#### Abstract

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 offer 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 Modified 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.

#### 1. Introduction

Soil shear stress-strain behaviour in small strain range is highly nonlinear and well documented in the literature of geotechnical engineering. Within the last 50 years, many researchers [1–7] have studied nonlinear shear stress-shear strain behaviour in different types of soils and proposed formulations describing the variation of shear modulus with strain in the small strain range [1, 8–12].

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. They 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. Thus, they do not consider nonlinear shear stress-strain behaviour accurately. Therefore, 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.

This paper updates the procedure for enhancing elastoplastic constitutive models with nonlinear shear stress-shear strain behaviour in the small strain range without changing other characteristics of the base models, first introduced in [18]. The presented algorithm preserves other features of the base constitutive model, including thermodynamic consistency. The 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.

The 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).

#### 2. Enhancing Constitutive Models with Small Strain Shear Nonlinearity

This 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. The 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.

This section will also show an example of enriching the constitutive model with small strain nonlinearity. The 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]. The 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.

##### 2.1. Brief Review of the Hyperbolic Models for the Small Strain Behaviour in Shear

The 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 is the secant shear modulus, represents the elastic shear modulus of soil before shearing, corresponds to the shear strain, and is the reference shear strain which is equal to the maximum shear stress divided to the elastic shear modules (). This formulation requires calibration of 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 is a curvature parameter and is the reference shear strain parameter at which the shear modulus reduces to 50% of 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 is the reference shear strain and equal to shear strain in which shear modulus reduces to 70% of its maximum value. This 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.

##### 2.2. Derivation of Constitutive Equations

Shearing of soils may cause both recoverable (elastic) and irrecoverable (plastic) deformations. Therefore, 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 is the deviatoric shear stress, is the deviatoric shear strain, and is the secant shear modulus which is a function of strain. Elastic and plastic parts are the two composing parts of , although only the elastic strain leads to stress increment:where is the elastic part of deviatoric strain and is the secant elastic shear modulus during small strain shearing. This 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 in the beginning of shearing has its maximum (initial) value which 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 is the plastic part of deviatoric strain. Equation (6) also provides a way to find shear stress based on the amount of plastic strain. This equation can also lead to the hardening law for small strain if and of small strain shearing are specified.

The secant shear modulus depends on the type of soil and the base constitutive model. Several formulations are available for defining . These 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 , the choice affects the complexity of the hardening law and may help to keep the base model intact. This 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 is the reference deviatoric shear strain and is equal to deviatoric shear strain in which shear modulus reduces to 70% of its maximum value.

The base constitutive model in this study uses a constant elastic shear modulus () and ignores variation of elastic shear modulus during shearing. Taking a different elastic shear modulus for small strain shearing formulation (i.e., assuming ) 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 of base model for small strain shearing (). This 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 towhere stands for error caused by constant . Combining this equation with equation (4) results in recalculation of equation (6) as

This paper introduces an internal variable aswhere is a scalar variable. This variable is negative when has higher value than and is similar to kinematic internal variable in Houlsby and Puzrin hyperplasticity framework [22], which will be investigated in next subsection. The internal variable replaces in calculation of hardening law to include errors caused by constant .

After defining and for small strain shearing, finding an equation relating to leads to calculation of the hardening law. Combining equations (4), (7), and (9) provides a function for based on :where is equal to and is . Inverse of (11) defines based on :where the value of is negative and is between 0 and .

Small strain hardening law is calculated after combining equations (4), (7), and (12):where is the shear hardening parameter. The 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. Therefore, the Modified Cam Clay yield surface remains unchanged during small strain shearing leading to preservation of the elastoplastic behaviour of the base model.

The small strain yield function separates the elastic and the elastoplastic regions. This function should remain 0 during small strain shearing and generate negative value beyond it. Defining small strain yield surface function assatisfies 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.

##### 2.3. Hyperplastic Formulation of the Small Strain Model

This subsection employs Houlsby and Puzrin [22] hyperplastic framework to derive the small strain model. The 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.

The following derivation utilizes a procedure similar to derivation of Modified Cam Clay model with constant [22]. It first introduces the Gibbs free energy function for small strain model. Then, 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 (). This procedure assumes that soil thermodynamic state is dependent on the internal kinematic variable () as well as the state variables (*p*, *q*). The definition of 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. This function depends on state and internal kinematic variables and allows for the calculation of the strains and the generalized stress:where is the volumetric strain and is the generalized stress conjugate of deviatoric internal variable. The study assumes the Gibbs free energy function as follows:where is the Modified Cam Clay hardening parameter and is a function of defined as

The 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. The remaining part of the function introduces the hardening behaviour during the small strain shearing. The Gibbs free energy function and equations (15)–(17) provide , , and as follows:

The dissipation function is another function required in hyperplasticity. This function must be nonnegative, homogenous of first order, and depend on the changes of internal variable . The dissipation function allows for calculation ofwhere is a dissipative generalized stress, conjugate to the deviatoric internal variable change. The study chooses dissipation function aswhich satisfies requirements of the dissipation function and leads to the calculation of as

In the hyperplastic framework, the Lagrange transformation of the dissipation function leads to calculation of the yield function. This yield function has to be a function of and is found as

Ziegler’s orthogonality condition states that . Therefore, the hyperplastic yield surface of equation (26) becomes

Combining equations (22) and (27) leads to the shear small strain yield surface definitionwhich is the same as small strain yield surface of equation (14). This 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. This conservative behaviour is inherited from the base model since the small strain shearing uses the elastic law of the base model.

##### 2.4. Incremental Stress-Strain Equations

This 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.

The 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 towhere , , , , , and are components of the stress tensor (). The differential of this yield surface isand it should remain equal to 0 during small strain shearing in order to satisfy the Prager consistency condition. In this equation, is the change of hardening parameter andwhere is change of internal soil variable. Equation (11) allows to calculate :where is the nonnegative deviatoric strain increment. This equation becomes undefined when is equal to and is thus only valid for below . Equation (32) also leads to a negative value for increment when the deviatoric strain increases and results in positive value for if the deviatoric strain decreases.

The stress increment in equation (30) is associated with change of elastic strain:where is the tangent elastic matrix for the small strain shearing and is the increment of the elastic strain. The elastic matrix of the base model is dependent on constant and should not be the same as which is constructed with . Nevertheless, to keep elastoplastic behaviour of the combined model same as the base model, the paper assumes that . Therefore, equation (33) changes towhere represents error caused by using as elastic matrix for small strain shearing. This error is similar to the error calculated in equation (8). Similarly, the new internal variable for soil iswhere is a tensorial variable similar to plastic strain. This variable replaces plastic strain in the incremental stress-strain equations and leads to the correct calculations of strain change:

The formulation uses an associated flow rule. Therefore, the increments of are computed aswhere 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 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 . These negative result in calculation of positive values for . The change of stress , which is calculated using equations (38) and (39), has a deviatoric part () equal to . Therefore, the small strain yield surface moves upward during the small strain shearing. The movement of the small strain yield surface does not produce any volumetric plastic strain. Therefore, the volumetric part of strain change () 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 , , and to be , , and , respectively. The figure shows that the proposed formulation replicates the nonlinear shear stress-shear strain behaviour defined by equation (3).

##### 2.5. Coupling Procedure

This subsection describes coupling between the small strain shear yield surface and the base constitutive model. This 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. The obtained model ensures a smooth transition between the small strain and the higher strain range shearing and presents a consistent unloading-reloading behaviour at any shear strain. The main goal of this coupling is to leave the base model unaffected, except for adding the small strain shearing yield surface.

The elastoplastic behaviour on the small strain shearing yield surface starts as soon as shearing of the soil begins. This elastoplastic behaviour produces deviatoric elastic, deviatoric plastic, and volumetric elastic strains. As the associated flow rule is assumed, the elastoplastic small strain shearing behaviour leads to no volumetric plastic strain. Therefore, 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. This smooth change happens automatically as shear strain reaches to . The small strain formulation is not valid for higher than , and the slope of curve has reduced to at this point, which is the slope of Modified Cam Clay elastic behaviour. Therefore, 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. The 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. The introduced internal variable increment () becomes equal to the increment of the deviatoric plastic strain. Therefore, 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. Therefore, the plastic behaviour of the Modified Cam Clay model remains intact on all the possible stress paths.

The unloading can be elastoplastic or purely elastic. That 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. The 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.

#### 3. Validation, Results, and Discussion

This section examines the performance of the proposed model based on the triaxial tests on clayey and sandy soils by Jardine et al. [2]. The shown simulations replicate Jardine et al.’s experiments on intact overconsolidated London Clay samples using (a) the Modified Cam Clay model with constant shear modulus and (b) the upgraded Modified Cam Clay model with small strain shearing capability, taking model parameters reported in [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]. The model driver uses explicit stress integration with error control and NICE technique [26].

The common triaxial tests conventionally use external strain measurement. The external measurement includes sitting, alignment, bedding, and compliance errors in strain ranging from 0.001% to 0.1% [27, 28]. Therefore, 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]. Therefore, the value of 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 .

Figure 6 illustrates the variation of secant shear modulus in two triaxial tests on intact London Clay. The 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.82 MPa and 32.21 MPa for of London Clay 1 (LC1) and London Clay 2 (LC2), respectively. These values and triaxial data in Figure 6 allow for estimation of and of enhanced model to be 0.035% and 2.5 MPa.

**(a)**

**(b)**

Table 1 summarizes all parameters of the Modified Cam Clay and the enhanced Modified Cam Clay model. The enhanced model predicts the variation of the secant shear modulus as shown in Figure 6. This 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. This behaviour is followed by a semilinear shear stress-shear strain relationship which continues upon reaching the peak deviatoric shear stress. The shear stress then reduces in both tests.

**(a)**

**(b)**

The 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. The soil then becomes elastoplastic, leading to the reduction of shear stress after the peak. The 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. The 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. The enhanced model shows the same behaviour of Modified Cam Clay in the higher strain range. The 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. This 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 for Ham River sand 1 (HRS1) and Ham River sand 2 (HRS2) tests, respectively. These values and the triaxial data of Figure 8 lead to an estimate and 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.

**(a)**

**(b)**

Wanatowski and Chu’s [30] investigation on the undrained behaviour of Changi sand shows a continuous 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–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 medium-dense Ham River sand results in similar behaviour. The 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–37]. Jardine et al. [2] reported that both HRS1 and HRS2 tests failed by cavitation of pore water and without reaching the critical state. Therefore, 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. The sand behaves in a highly nonlinear fashion in the initial stages of shearing. This nonlinear behaviour changes to a linear one as shearing progresses. Both tests exhibit linear behaviour after reaching approximately 0.3% deviatoric shear strain. The 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.

**(a)**

**(b)**

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. The results also show that the proposed formulation cannot predict any increase of the secant shear modulus. The 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. This 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. Therefore, the procedure satisfies the main goal of this study.

##### 3.1. Nonconstant Shear Modulus

The 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 is no longer a parameter but rather would change with variations of base model elastic shear modulus. Then, 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) should becomewhere denotes the change of the mean stress. This 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. Therefore, 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.

#### 4. Conclusions

This study addresses the need for simple constitutive models capable of capturing realistic nonlinear shear stress-shear strain behaviour in the small strain range. The 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. The 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.

The introduced procedure can upgrade the capability of many elastoplastic constitutive models. The 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. Therefore, 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. This study focuses on small strain shearing behaviour of soils and thus volumetric stress-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. Therefore, 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.

#### Data Availability

All the data used to support the findings of this study are included within the article. A matlab code for adding the introduced small strain shearing formulation to the constitutive models with constant shear modulus is provided in the Supplementary Material.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was partially funded by the Academy of Finland under the project “Progressive failure and postfailure modelling of slopes with Generalized Interpolation Material Point Method (GIMP)” (decision no. 286628).

#### Supplementary Materials

A matlab code for adding the introduced small strain shearing formulation to the constitutive models with constant shear modulus is provided.* ((Supplementary Materials))*