Enhancing Constitutive Models for Soils: Adding the Capability to Model Nonlinear Small Strain in Shear

. The deviatoric stress-deviatoric strain relationship in soils is highly nonlinear, especially in the small strain range. However, the constitutive models which aim to replicate the small strain nonlinearity are often complex and rarely used in geotechnical engineering practice. The goal of this study is to oﬀer a simple way for updating the existing constitutive models, widely used in geotechnical practice, to take into account the small strain shear modulus changes. The study uses an existing small strain relationship to derive a yield surface. When the yield surface is introduced to an existing soil model, it enhances the model with the nonlinear deviatoric stress-deviatoric strain relationship in the small strain range. The paper also gives an example of such a model enhancement by combining the new yield surface with the Modiﬁed Cam Clay constitutive model. The validation simulations of the undrained triaxial tests on London Clay and Ham River sand with the upgraded constitutive models replicate the experiments clearly better than the base models, without any changes to existing model parameters and the core source code associated with the base model.

Despite the extensive investigations on the nonlinear shear stress-strain behaviour of soil, the most common soil constitutive models in engineering practice ignore this aspect of soil behaviour.ey employ linear or simple nonlinear elasticity instead, which fails to describe the nonlinear shear stress-shear strain behaviour in the small strain range [13].Furthermore, advanced constitutive models capable of describing nonlinear behaviour of soils in the small strain range are often complex and require special soil investigation and extra parameters beyond those related to small strain behaviour [10], which prevents engineers from applying them in routine practice.Many advanced constitutive models, which simulate special aspects of soil behaviour, such as Barcelona basic model (BBM) [14] and constitutive models used for simulation of large strain in soils [15,16], are based on commonly used soil constitutive models.us, they do not consider nonlinear shear stress-strain behaviour accurately.erefore, there is a need for simple constitutive models capable of capturing realistic nonlinear shear stress-shear strain behaviour in the small strain range.
Benz [17] introduced a small strain overlay model which can replicate small strain behaviour of soil with only 2 extra parameters.Combining this overlay model with elastoplastic constitutive models enhances them with nonlinear shear stress-shear strain behaviour.However, the combination with the overlay small strain model requires alternation of some other aspects of elastoplastic constitutive models and may lead to thermodynamic inconsistency.
is paper updates the procedure for enhancing elastoplastic constitutive models with nonlinear shear stressshear strain behaviour in the small strain range without changing other characteristics of the base models, first introduced in [18].e presented algorithm preserves other features of the base constitutive model, including thermodynamic consistency.e proposed procedure can upgrade constitutive models with a constant and nonlinear shear modulus.Furthermore, implementation of the proposed algorithm typically does not require any changes of the code associated with the base model and therefore can be added to an existing code as a patch.
e paper first derives the constitutive equations for small strain shearing of a base model (i), shows the thermodynamic compatibility of that small strain shearing model (ii), gives incremental stress-strain equations for the small strain shearing and a procedure to combine the base constitutive model and small strain shearing (iii), shows the improvement in replication of laboratory tests on clayey and sandy soil by Modified Cam Clay and Mohr-Coulomb constitutive models enhanced with the proposed small strain shear nonlinear behaviour (iv), and finally discusses the required steps for enhancing other base constitutive models with the capability to model the small strain nonlinearity (v).

Enhancing Constitutive Models with Small
Strain Shear Nonlinearity is section describes derivation of a yield surface for describing small strain shearing behaviour which can be introduced into simple constitutive models, with no other change in the constitutive behaviour of the base model.e section starts with (1) brief review of the small strain models followed by (2) the derivation of the constitutive equations for small strain shearing yield surface, (3) the introduction of hyperplasticity formulations for small strain shearing, (4) the development of the incremental stress-strain equations during small strain shearing, and (5) the coupling between the base constitutive model and the small strain one.is section will also show an example of enriching the constitutive model with small strain nonlinearity.e base model used in this example, the Modified Cam Clay model with constant shear modulus, is a very common model in geotechnical engineering practice.It is also a base for some advanced constitutive model (e.g.BBM) and can be formulated in thermodynamically consistent fashion [19].e same procedure can upgrade constitutive models with isotropic hardening as well as those with constant yield surface, as long as their formulation uses a constant elastic shear modulus.Upgrading constitutive models with a nonconstant shear modulus requires modification to the given procedure and is beyond the scope of this paper.Nevertheless, the paper will discuss steps required for enhancing such constitutive models in the next section.

Brief Review of the Hyperbolic Models for the Small Strain Behaviour in Shear.
e limited selection of models discussed here aim to alter the shear modulus in small strain range in order to capture the significant nonlinearity of the shear strain-shear stress curve.First, Hardin and Drnevich [1] proposed one of the first hyperbolic formulations: where G is the secant shear modulus, G 0 represents the elastic shear modulus of soil before shearing, c corresponds to the shear strain, and c r is the reference shear strain which is equal to the maximum shear stress divided to the elastic shear modules (τ max /G 0 ). is formulation requires calibration of c r which is often difficult as the maximum shear stress changes with mean effective stress level [10].In order to allow a stress independent calibration, Darendeli [8] suggested a slightly altered formulation for sand: where a is a curvature parameter and c r is the reference shear strain parameter at which the shear modulus reduces to 50% of G 0 value.To address both clayey and sandy soils, Correia et al. [9] and dos Santos et al. [11,12] introduced a formulation to predict secant shear modulus in the small strain range for both clayey and sandy soils: where c r is the reference shear strain and equal to shear strain in which shear modulus reduces to 70% of its maximum value.is equation is only applicable in the small strain region and can replicate the nonlinear behaviour of both clayey and sandy soils.Due to that comprehensiveness, equation ( 3) is selected to be introduced into constitutive models.

Derivation of Constitutive
Equations.Shearing of soils may cause both recoverable (elastic) and irrecoverable (plastic) deformations.erefore, its modelling requires, for example, an elastoplastic formulation.Such an elastoplastic formulation needs definitions of a hardening law, a yield surface, and a plastic potential function.
Equations ( 1)-(3), and similar expressions, provide nonlinear relation between shear stress and total shear strain: where q is the deviatoric shear stress, ε q is the deviatoric shear strain, and G is the secant shear modulus which is a function of strain.Elastic and plastic parts are the two composing parts of ε q , although only the elastic strain leads to stress increment: where ε e q is the elastic part of deviatoric strain and G e ss is the secant elastic shear modulus during small strain shearing.
is secant elastic modulus controls soil behaviour when subjected to very small shear strain at moderate rate.Experiments, based on the measurements of the shear wave velocity with bender elements, show that G e ss in the beginning of shearing has its maximum (initial) value which 2 Advances in Civil Engineering changes with the density and the stress state of soil [20,21].
During shearing, the value of the secant elastic shear modulus decreases.Combining equation ( 4) and ( 5) leads to the calculation of plastic strain, based on deviatoric shear stress level: where ε p q is the plastic part of deviatoric strain.Equation ( 6) also provides a way to find shear stress based on the amount of plastic strain.is equation can also lead to the hardening law for small strain if G and G e ss of small strain shearing are specified.
e secant shear modulus G depends on the type of soil and the base constitutive model.Several formulations are available for defining G. ese equations describe different soils, vary in accuracy, and depend on different variables and parameters.Although it is theoretically possible to use almost any formulation for the secant shear modulus G, the choice affects the complexity of the hardening law and may help to keep the base model intact.is paper uses equation (3) for describing nonlinear small strain behaviour of soil due to its simplicity and applicability for both clayey and sandy soils.Equation (3) may be recast into the deviatoric stress space: where ε qr is the reference deviatoric shear strain and is equal to deviatoric shear strain in which shear modulus reduces to 70% of its maximum value.e base constitutive model in this study uses a constant elastic shear modulus (G e MCC ) and ignores variation of elastic shear modulus during shearing.Taking a different elastic shear modulus for small strain shearing formulation (i.e., assuming G e ss ≠ G e MCC ) is a more realistic approach but it would lead to deviation from behaviour of the base model and is against the main goal of this section.Hence, the study assumes the constant G e MCC of base model for small strain shearing (G e ss � G e MCC ). is constant value cannot represent variation of shear modulus during small strain shearing.Inevitably, it leads to erroneous calculation of elastic strain and the equation (5) changes to where ε Error q stands for error caused by constant G e .Combining this equation with equation (4) results in recalculation of equation (6) as is paper introduces an internal variable as where α q is a scalar variable.is variable is negative when G has higher value than G e MCC and is similar to kinematic internal variable in Houlsby and Puzrin hyperplasticity framework [22], which will be investigated in next subsection.e internal variable α q replaces ε p q in calculation of hardening law to include errors caused by constant G e MCC .After defining G and G e ss for small strain shearing, finding an equation relating ε q to α q leads to calculation of the hardening law.Combining equations ( 4), (7), and ( 9) provides a function for α q based on ε q : where A 1 is equal to G 0 /G e MCC and A 2 is 0.385/ε qr .Inverse of (11) defines ε q based on α q : where the value of α q is negative and ε q is between 0 and . Small strain hardening law is calculated after combining equations ( 4), (7), and ( 12): where q 0 is the shear hardening parameter.e added −1 in this equation is an extra change of the variable and helps to simplify the hyperplastic derivation of the model.Equation ( 13) is independent from the volumetric plastic strain.erefore, the Modified Cam Clay yield surface remains unchanged during small strain shearing leading to preservation of the elastoplastic behaviour of the base model.
e small strain yield function separates the elastic and the elastoplastic regions.
is function should remain 0 during small strain shearing and generate negative value beyond it.Defining small strain yield surface function as satisfies these criteria.Plastic potential function defines direction of plastic strain during elastoplastic behaviour.As this study assumes associated flow rule, the plastic potential function is the same as the yield surface.

Hyperplastic Formulation of the Small Strain Model.
is subsection employs Houlsby and Puzrin [22] hyperplastic framework to derive the small strain model.e hyperplastic derivation provides an alternative way to look at the model enhancement.It also ensures thermodynamic consistency, as the first and second laws of thermodynamics are automatically enforced in the hyperplastic framework.
e following derivation utilizes a procedure similar to derivation of Modified Cam Clay model with constant G e [22].It first introduces the Gibbs free energy function for small strain model.en, it defines the dissipation function and uses it for finding the yield equation.Finally, Ziegler's orthogonality condition leads to recalculation of the shear small strain yield surface given in the previous subsection (F ss ). is procedure assumes that soil thermodynamic state Advances in Civil Engineering is dependent on the internal kinematic variable (α q ) as well as the state variables (p, q). e definition of α q is the same as in the last subsection (equation ( 10)).
Gibbs free energy is one of the possible options for defining the soil energy equation.is function depends on state and internal kinematic variables g � g(p, q, α q ) and allows for the calculation of the strains and the generalized stress: where ε p is the volumetric strain and χ q is the generalized stress conjugate of deviatoric internal variable.e study assumes the Gibbs free energy function as follows: where p 0 is the Modified Cam Clay hardening parameter and m(α q ) is a function of α q defined as e first two terms in equation ( 18) relate to the elastic behaviour of small strain model, whereas the third term adds the deviatoric internal variable, ensuring correct calculation of the deviatoric strain.e remaining part of the function introduces the hardening behaviour during the small strain shearing.
e Gibbs free energy function g(p, q, α q ) and equations ( 15)-( 17) provide ε p , ε q , and χ q as follows: e dissipation function is another function required in hyperplasticity.
is function must be nonnegative, homogenous of first order, and depend on the changes of internal variable d � d(δα q ).e dissipation function allows for calculation of where χ q is a dissipative generalized stress, conjugate to the deviatoric internal variable change.e study chooses dissipation function as which satisfies requirements of the dissipation function and leads to the calculation of χ q as χ q � zd zδα q � sgn δα q  .
In the hyperplastic framework, the Lagrange transformation of the dissipation function leads to calculation of the yield function. is yield function has to be a function of χ q and is found as Ziegler's orthogonality condition states that χ i � χ i .erefore, the hyperplastic yield surface of equation ( 26) becomes Combining equations ( 22) and ( 27) leads to the shear small strain yield surface definition which is the same as small strain yield surface of equation (14). is shows that the small strain shearing model is well defined within the hyperplastic framework.It also confirms that the proposed model satisfies the first and second laws of thermodynamics.is conservative behaviour is inherited from the base model since the small strain shearing uses the elastic law of the base model.

Incremental Stress-Strain Equations.
is subsection focuses on derivation of incremental stress-strain equations in the small strain range.It first recasts the derived yield surface into the general stress space.Next, it introduces a new internal variable to compensate using the elastic matrix of the Modified Cam Clay model instead of the small strain shearing tangent matrix.Finally, it uses the consistency condition to calculate plastic multiplier and the increment of stress.Combining these derived equations with the base model allows for replication of the small strain nonlinearity during shearing.
e enhancement of the Modified Cam Clay model with the small strain shear nonlinearity requires the small strain model to be formulated in the general stress space.Recasting the yield surface into the general stress space leads to 4 Advances in Civil Engineering where σ 11 , σ 22 , σ 33 , σ 12 , σ 13 , and σ 23 are components of the stress tensor (σ).e differential of this yield surface is and it should remain equal to 0 during small strain shearing in order to satisfy the Prager consistency condition.In this equation, dq 0 is the change of hardening parameter and where dα q is change of internal soil variable.Equation (11) allows to calculate dα q : where dε q is the nonnegative deviatoric strain increment.is equation becomes undefined when ε q is equal to and is thus only valid for . Equation ( 32) also leads to a negative value for dα q increment when the deviatoric strain increases and results in positive value for dα q if the deviatoric strain decreases.
e stress increment dσ in equation ( 30) is associated with change of elastic strain: where D el ss is the tangent elastic matrix for the small strain shearing and dε e is the increment of the elastic strain.e elastic matrix of the base model D where dε Error represents error caused by using D el MCC as elastic matrix for small strain shearing.is error is similar to the error calculated in equation (8).Similarly, the new internal variable for soil is where β is a tensorial variable similar to plastic strain.is variable replaces plastic strain in the incremental stressstrain equations and leads to the correct calculations of strain change: e formulation uses an associated flow rule.erefore, the increments of β are computed as where λ is the scalar plastic multiplier.Combining (34), (36), and (37) gives the stress increment as Introducing ( 31), (32), and ( 38) into (30) results in calculation of the plastic multiplier: Calculation of the plastic multiplier allows for calculation of dσ when the stress state satisfies the yield equation of small strain shearing.
Figure 1 shows how the yield surface evolve during small strain shearing.Before shearing starts, the small strain yield surface is a horizontal line located on the deviatoric strain axis.Increases of deviatoric strain result in negative values for dα q .ese negative dα q result in calculation of positive values for dq 0 .e change of stress dσ, which is calculated using equations ( 38) and (39), has a deviatoric part (dq) equal to dq 0 .erefore, the small strain yield surface moves upward during the small strain shearing.e movement of the small strain yield surface does not produce any volumetric plastic strain.erefore, the volumetric part of strain change (dε p ) results in the mean stress change according to the elastic law of the base model.Furthermore, application of a strain increment consisting only from the volumetric strain leads to no changes in the small strain yield surface.
Figure 2 compares the variation of the secant shear modulus calculated with the developed elastoplastic formulation with that of equation ( 3), assuming G 0 , G e , and c r to be 30 MPa, 1 MPa, and 0.0525%, respectively.e figure shows that the proposed formulation replicates the nonlinear shear stress-shear strain behaviour defined by equation (3).

Coupling Procedure.
is subsection describes coupling between the small strain shear yield surface and the base constitutive model.is coupling is shown on an example of a constitutive model which employs the introduced yield surface to enhance the Modified Cam Clay model with the small strain behaviour.
e obtained model ensures a smooth transition between the small strain and the higher strain range shearing and presents a consistent unloadingreloading behaviour at any shear strain.e main goal of this coupling is to leave the base model unaffected, except for adding the small strain shearing yield surface.
e elastoplastic behaviour on the small strain shearing yield surface starts as soon as shearing of the soil begins.is elastoplastic behaviour produces deviatoric elastic, deviatoric Advances in Civil Engineering plastic, and volumetric elastic strains.As the associated flow rule is assumed, the elastoplastic small strain shearing behaviour leads to no volumetric plastic strain.erefore, it does not affect the Modified Cam Clay yield surface.Formulations derived in the previous subsections (equations 29-39) allow for calculation of incremental stress-strain behaviour on this yield surface. Figure 1 shows a schematic of yield surface changes during small strain loading, assuming no yielding on the Modified Cam Clay yield locus.
Shearing behaviour of soil should change between small strain and higher strain range in a smooth manner.is smooth change happens automatically as shear strain reaches to ). e small strain formulation is not valid for ε q higher than (1/A 2 )( �� � A 1  − 1), and the slope of q − ε q curve has reduced to 3G e MCC at this point, which is the slope of Modified Cam Clay elastic behaviour.erefore, q 0 does not change when the deviatoric strain increases further, leading to the stress state no longer being on the small strain yield surface and a smooth change of behaviour to the purely elastic (Figure 3).
If the considered stress state is not on any of the yield surfaces (which requires some initial yielding on the introduced small strain yield surface), the elastic rule controls the soil behaviour until the stress state reaches one or both yield surfaces again.e enhanced model uses the elastic law of the base model.Figure 3 shows a schematic of transition between the two shearing mechanisms, elastic shearing, and position of yield surface in different situations.
After reaching point 3, the behaviour is elastic until reaching the Modified Cam Clay model surface (the base model).At that time, in general, during yielding, the deviatoric elastic, volumetric elastic, deviatoric plastic, and volumetric plastic strains are generated.e introduced internal variable increment (δα q ) becomes equal to the increment of the deviatoric plastic strain.erefore, the deviatoric plastic strain of the Modified Cam Clay model affects the small strain yield surface (moving it downwards until it reaches 0 and becoming nonactive).Afterwards, the changes of the deviatoric plastic strain lead to no changes in the deviatoric hardening parameter of the small strain shearing.Figure 4 shows a schematic of the small strain yield surfaces variation during loading on the Modified Cam Clay yield surface.Notably, changes in the volumetric plastic strain affect only the base model yield surface.However, both the Modified Cam Clay and the small strain shearing yield surface can change when loading happens on the Modified Cam Clay yield surface.
If the stress state satisfies both yield equations of the Modified Cam Clay and the small strain shearing, the plastic deviatoric strain obtained from the Modified Cam Clay model will lead to the stress state moving away from the small strain yield surface, hence yielding will occur only on the base model yield surface. erefore, the plastic behaviour of the Modified Cam Clay model remains intact on all the possible stress paths.
e unloading can be elastoplastic or purely elastic.at depends on the stress state at the beginning of the unloading.If the unloading starts or reaches the still active small strain yield surface, at that point, the elastoplastic deformations occur.
e unloading would move the small strain yield surface downward leading to the reversal of the small strain loading.However, the unloading is purely elastic if the previous loading on the Modified Cam Clay yield surface deactivated the small strain yield surface, i.e., it has moved the small strain yield surface to q � 0. Figure 5 shows a schematic of unloading from different points in the deviatoric strain-deviatoric stress space.In this figure, the continuous line represents the loading of the soil while the dashed lines show different cases of unloading.

Validation, Results, and Discussion
is section examines the performance of the proposed model based on the triaxial tests on clayey and sandy soils by Small strain shearing Figure 1: Schematic of small strain yield surface and its movements.F ssi denotes small strain yield surface position at point i, where i � 0, 1, and 2.   [23,24].Furthermore, the section shows simulations of Jardine et al.'s experiments on a Ham River sand (HRS) using a linear elastic Mohr-Coulomb model with an associated flow rule.Simulations on HRS use parameters from Jovcic [25].e model driver uses explicit stress integration with error control and NICE technique [26].e common triaxial tests conventionally use external strain measurement.
e external measurement includes sitting, alignment, bedding, and compliance errors in strain ranging from 0.001% to 0.1% [27,28].erefore, conventional triaxial test cannot realistically investigate small strain shearing behaviour in soil.However, local strain measurement of samples allows for avoiding errors in the small strain range [2,6].As the accuracy of local strain measurement is limited, strain measurement of below 0.0001% requires a different measuring technique.Dynamic soil testing is a suitable tool for very small strain measurement [6].
erefore, the value of G 0 in the simulations estimated using bender elements, a dynamic soil testing method.Viggiani and Atkinson [29] performed bender element tests on London Clay.Jovcic [25] provided similar data for estimation of Ham River sand G 0 .
Figure 6 illustrates the variation of secant shear modulus in two triaxial tests on intact London Clay.e figure shows an initial increase of secant shear modulus in both tests.We concluded that factors such as tilting of the overconsolidated samples and defects inside intact samples are likely reason for the increase.Figure 6 also shows hyperbolic reduction of secant shear modulus after the initial increase in both tests.Equation and data from Viggiani and Atkinson [29] on London Clay led to estimation of 36.82MPa and 32.21 MPa for G 0 of London Clay 1 (LC1) and London Clay 2 (LC2), respectively.ese G 0 values and triaxial data in Figure 6 allow for estimation of ε qr and G e MCC of enhanced model to be 0.035% and 2.5 MPa.Table 1 summarizes all parameters of the Modified Cam Clay and the enhanced Modified Cam Clay model.e enhanced model predicts the variation of the secant shear modulus as shown in Figure 6. is figure shows that the enhanced model cannot replicate the initial increase of the secant shear modulus (likely the artifact of some experimental inaccuracy) but can capture the remaining hyperbolic reduction of the modulus.
Figure 7 presents the variation of the deviatoric stress in the laboratory tests and compares them with prediction of the Modified Cam Clay and the enhanced Modified Cam Clay models.Results of the triaxial tests in Figure 7 show the highly nonlinear behaviour of soil during the initial stages of shearing.is behaviour is followed by a semilinear shear stress-shear strain relationship which continues upon reaching the peak deviatoric shear stress.e shear stress then reduces in both tests.
e simulation with the Modified Cam Clay model, with the constant shear modulus (dashed line in Figure 7) is elastic before reaching the peak shear stress in this heavily overconsolidated soil.e soil then becomes elastoplastic, leading to the reduction of shear stress after the peak.e elastic behaviour of the Modified Cam Clay model cannot replicate the highly nonlinear behaviour in the initial stages of shearing.However, the simulation based on the Modified Cam Clay model provides a good description of postpeak behaviour of soil.e simulation with the enhanced Modified Cam Clay model results in the correct behaviour during small strain shearing and enforces a smooth transition to the elastic behaviour of Modified Cam Clay in higher strain range.Figure 7 shows that the initial elastoplastic behaviour of the enhanced model replicates the nonlinear behaviour in the small strain range well.e enhanced model shows the same behaviour of Modified Cam Clay in the higher strain range.
e improvements in the small strain range lead to overall improved prediction of clayey soil behaviour.
Figure 8 shows the secant shear modulus in the two triaxial tests on Ham River sand and compares them with the prediction of the enhanced Mohr-Coulomb model. is figure shows an initial increase of the secant shear modulus in both tests, which is followed by a hyperbolic reduction.Data of bender element tests in [25] allow to estimate 95 MPa and 147 MPa for G 0 for Ham River sand 1 (HRS1) and Ham River sand 2 (HRS2) tests, respectively.ese G 0 values and the triaxial data of Figure 8 lead to an estimate ε qr and G e MC for the enhanced Mohr-Coulomb to be 0.04% and 11.5 MPa. Figure 8 illustrates that the enhanced model cannot predict the initial variations of secant shear modulus (again, being likely an artifact of the experimental technique) but can capture the remaining variation well.
Wanatowski and Chu's [30] investigation on the undrained behaviour of Changi sand shows a continuous  Advances in Civil Engineering strain-hardening behaviour for the medium-dense sand, which results in effective stress path approaching a constant stress ratio (asymptotic behaviour).Similar behaviour is also reported by other researchers and in different conditions [31][32][33][34].In an undrained triaxial test on medium-dense Changi sand, this constant stress ratio mobilizes a friction angle of 35 °which is slightly higher than 33 °critical state friction angle [30].Undrained triaxial shearing of mediumdense Ham River sand results in similar behaviour.e effective stress path of the HRS1 and HRS2 approaches a constant stress path with mobilized friction angle of 35 °that is slightly higher than 33 °critical state friction angle of Ham River sand reported in [35][36][37].Jardine et al. [2] reported that both HRS1 and HRS2 tests failed by cavitation of pore water and without reaching the critical state.erefore, the modelling uses mobilized friction angle of constant stress path (instead of critical state) to simulate these tests.Table 2 provides a summary of parameters used for the Mohr-Coulomb and enhanced Mohr-Coulomb constitutive models.Figure 9 illustrates the variation of the deviatoric shear stress in undrained triaxial tests on the medium-dense Ham River sand and compares them with the predictions of the Mohr-Coulomb and enhanced Mohr-Coulomb models.e sand behaves in a highly nonlinear fashion in the initial stages of shearing.is nonlinear behaviour changes to a linear one as shearing progresses.Both tests exhibit linear behaviour after reaching approximately 0.3% deviatoric shear strain.e figure shows that the simulation with the linear elastic Mohr-Coulomb model is incapable of capturing the initial stage of shearing, which leads to poor prediction in higher strain range.Figure 9 demonstrates how adding the small strain yield surface improves the Mohr-Coulomb model and allows for a better prediction of the soil behaviour.Furthermore, this figure shows that behaviour of the enhanced model and the base model is only different in the small strain range.
Shown results demonstrate the capabilities of the proposed framework.Figures 6 and 8 show that the enhanced models can replicate hyperbolic degradation of the secant shear modulus in clayey and sandy soils.e results also show that the proposed formulation cannot predict any increase of the secant shear modulus.e increases of experimental data in the very first stages of shearing are likely to be caused by experimental inaccuracy and cannot be predicted since equation (3) disregards them.
Shown results in Figures 7 and 9 confirm the capability of the proposed framework for the prediction of the highly nonlinear shear stress-shear strain behaviour of clayey and sandy soils. is improves the base constitutive models and allows them to model the soil behaviour more accurately.Furthermore, the results demonstrate that the proposed model enhancement does not affect the predicted soil behaviour except adding the small stain shearing yield surface. erefore, the procedure satisfies the main goal of this study.

Nonconstant Shear Modulus.
e shown procedure, which introduces an extra yield surface, can enhance the base constitutive models with nonconstant elastic shear modulus after some modifications.First, equation (12) and its range of applicability should be redefined since A 1 is no longer a parameter but rather would change with variations of base model elastic shear modulus.en, the hardening law and the small strain shearing yield surface should be adjusted accordingly.Furthermore, the differential of the yield surface may need to be adjusted depending on how the elastic shear modulus of the base constitutive model is defined.For example, if the elastic shear modulus of base is defined as a function of the mean pressure, equation ( 30 where dp denotes the change of the mean stress.is variation in the differential of the yield surface would lead to alteration in all the following incremental stress-strain equations.Furthermore, taking the base constitutive model as a model with a nonconstant elastic shear modulus may lead to a thermodynamic incompatibility during the small strain shearing as the introduced procedure uses the elastic law of base model and inherits the thermodynamic compatibility (or lack of it) from it.erefore, in case the base constitutive model with a nonconstant elastic shear modulus being nonconservative [38], the enhanced model with small strain shearing will not be derivable within the hyperplasticity framework.

Conclusions
is study addresses the need for simple constitutive models capable of capturing realistic nonlinear shear stress-shear strain behaviour in the small strain range.e paper first introduces a procedure for enhancing the base constitutive models, which use a constant shear modulus with the ability to replicate the nonlinear small strain behaviour.e paper later shows the thermodynamic compatibility of a model enhanced with the extra yield surface.Finally, the study illustrates the capabilities of the proposed approach by comparing the results of simulations with the enhanced models to the laboratory results on different type of soils and briefly discusses the required steps for enhancing constitutive models that use nonconstant shear modulus.
e introduced procedure can upgrade the capability of many elastoplastic constitutive models.
e upgraded models have better ability to describe the shearing of soils in the small strain region without affecting any other aspect of the base model behaviour.Furthermore, in the numerical implementation, the procedure does not require any changes in the original constitutive model as the yield surfaces of the base model and the small strain enhancement are independent.erefore, the procedure may be used to enhance constitutive models without the need of altering the part of the source code associated with the base model.
A possible further research may develop similar enhancement related to the volumetric behaviour of soils in the small strain range.
is study focuses on small strain shearing behaviour of soils and thus volumetric stress- q (kPa)  Advances in Civil Engineering volumetric strain behaviour of developed formulation is inherited from the base model and remains unchanged.Current treatments of the small strain volumetric behaviour in literature lead to change in other aspects of soil behaviour.erefore, development of a formulation for enhancement of volumetric stress-volumetric small strain behaviour without changing other aspect of soil behaviour is a potential expansion of the current study.

F ss3 and F ss4 4 Figure 3 :
Figure3: Schematic of small strain loading.F ssi denotes small strain yield surface position at point i, where i � 1, 2, 3, and 4.

Figure 5 :
Figure 5: Schematic of unloading from different points.

FF ss7 Figure 4 :
Figure4: Changes of small strain yield surface during loading on the Modified Cam Clay.F ssi denotes small strain yield surface position at point i, where i � 3-7.

Figure 6 :
Figure 6: Variation of secant shear modulus in London Clay based on data from Jardine et al. (1984)

Figure 8 :
Figure 8: Variation of secant shear modulus in Ham River sand based on data from Jardine et al. [2].

Figure 7 :
Figure 7: Variations of deviatoric stress in London Clay based on data from Jardine et al. [2].

Figure 9 :
Figure 9: Variations of deviatoric stress in Ham River sand based data from Jardine et al. [2].

Table 1 :
Parameters of Modified Cam Clay and enhanced Modified Cam Clay model.

Table 2 :
Parameters of Mohr-Coulomb and enhanced Mohr-Coulomb model.