A Robust Constitutive Model for Clays over a Wide Range of Plasticity

1. Introduction

In practical geotechnical engineering, clay deposits rarely consist of a homogeneous material with a uniform overconsolidation ratio (OCR). Instead, they often contain strata exhibiting variations in OCR [1,2,3,4,5,6,7,8]. Consequently, geotechnical analyses should adopt a robust constitutive model capable to rigorously predicts the mechanical response of both normally consolidated (NC) and overconsolidated (OC) clays over a wide range of OCRs formulated within a unified theoretical framework.

The original Cam Clay (OCC) [9] and Modified Cam Clay (MCC) [10] models were the first two elastoplastic constitutive models originally developed to describe the behavior of normally consolidated (NC) clays within the framework of Critical State Soil Mechanics (CSSM). It is evident that, in subsequent studies, numerous modifications of the MCC model have been proposed to improve its applicability to geomaterials (e.g., [11,12,13,14,15,16,17,18,19,20,21,22]). For overconsolidated (OC) clays, the yield function (f) of the MCC model—characterized by its elliptical shape—tends to overestimate the peak strength of heavily overconsolidated materials on the dry side of the critical state [23,24,25,26,27]. Therefore, the dry side portion of the MCC yield surface is commonly revised to better capture the OC behavior. A key feature in such reformulations is the incorporation of the Hvorslev envelope [28] into the yield function, resulting in the familiar bullet- or teardrop-shaped yield surfaces. Embedding the Hvorslev envelope in the MCC yield function enhances peak-strength predictions for moderately to heavily overconsolidated clays while maintaining a smooth surface-gradient transition across all stress ratios.

On the wet side of the bullet/teardrop yield surface, the formulation generally remains similar to the Modified Cam Clay (MCC) model. This portion of the surface provides a convex stress path for normally consolidated (NC) clay (as illustrated in Figure 1a–b) and performs well for clays that exhibit relatively low plastic deformation under low stress ratios, $\mathrm{\eta }$ = q/p (hereafter referred to as LPLS). In Figure 1 and the related text, the mean effective stress p and deviatoric stress q are defined following the fundamental principle of soil mechanics as p = ($\mathrm{\sigma }1 + \mathrm{\sigma }2 + \mathrm{\sigma }3$)/3 and q = $\mathrm{\sigma }1 - \mathrm{\sigma }3$ where $\mathrm{\sigma }1$, $\mathrm{\sigma }2$ and $\mathrm{\sigma }3$ are the principal stresses satisfying $\mathrm{\sigma }1 \ge \mathrm{\sigma }2 \ge \mathrm{\sigma }3$. The parameter ${\mathrm{p}}_0$ denotes the preconsolidation pressure, representing the maximum historical mean effective stress that controls the current yield and bounding surfaces. The Critical State Line (CSL), defined by $\mathrm{q} = \mathrm{Mp}$, marks the ultimate stress condition where plastic volumetric strain ceases. The normalized form $\mathrm{p}/{\mathrm{p}}_0- \mathrm{q}/{\mathrm{p}}_0$ was used to remove dimensional influence and enable direct comparison among clays with different consolidation histories. This normalization also highlights the relative position of the current stress state with respect to its yield surface, thus facilitating a unified interpretation for both normally and overconsolidated clays.

However, certain clays deviate from this behavior: while the models may predict responses reasonably well for overconsolidated (OC) clays with OCR > 1, they often fail to reproduce NC behavior with comparable fidelity. This limitation diminishes their practical applicability in geotechnical analysis, as inadequate representation of even a single clay stratum can compromise the reliability of the overall boundary-value solution. Such clays frequently follow nearly linear stress paths in the p–q plane, rather than the curved trajectories shown in Figure 1c–d, reflecting the onset of substantial plastic straining from the very early stages of loading—even under low stress ratios (denoted as HPLS in this study).

In addition, in modeling normally consolidated (NC) clay, it is generally assumed that no plastic strain develops inside the yield surface. However, for overconsolidated (OC) clays, plastic strains may occur even within the yield surface. To capture this behavior, a bounding surface plasticity framework was introduced for OC clays (e.g., 29-30,32-33) with the adaptation of a non-associated flow rule, in which the yield function differs from the plastic potential function $(\mathrm{f} \ne \mathrm{g})$. Although the exact form of g is not explicitly known, it can be evaluated based on the principle of stress–dilatancy D, defined as the ratio of plastic volumetric strain to plastic deviatoric strain, given by $\mathrm{D} = (\partial \mathrm{g}/\partial \mathrm{p})/(\partial \mathrm{g}/\partial \mathrm{q})$. It should be noted that, when applying D in constitutive formulations, a paradox arises because $\mathrm{D}$ is often assumed to equal $(\partial \mathrm{f}/\partial \mathrm{p})/(\partial \mathrm{f}/\partial \mathrm{q})$ instead of $(\partial \mathrm{g}/\partial \mathrm{p})/(\partial \mathrm{g}/\partial \mathrm{q})$ (e.g., [29,30,32]). Therefore, to eliminate this paradox in models employing the bounding surface plasticity framework, an appropriate specification of the plastic potential function g is required.

The primary objective of this study Is to develop a unified constitutive model capable of capturing the behavior of both normally consolidated (NC) and overconsolidated (OC) clays, with particular emphasis on NC clays that exhibit both low- and high-plastic deformation under low stress ratios. To achieve this objective, several new theoretical features are introduced. First, a dual-parameter control ($\mathrm{\Psi }$, $\mathrm{\Omega }$) is embedded into the yield function to enable flexible and explicit regulation of the yield-surface skewness and curvature. Second, an explicitly defined plastic potential function is formulated to eliminate the long-standing stress–dilatancy paradox in bounding-surface-type models. Third, a linear stress–dilatancy relation is adopted to ensure physical transparency and facilitate practical parameter calibration. Through these features, the proposed formulation provides a unified framework that can consistently capture both LPLS and HPLS responses within a single constitutive structure, thereby advancing the state of the art in constitutive modeling of clays under low stress ratios.

2. Formulations of the Constitutive Model

2.1. Yield and Plastic Potential Functions

In this section, for the sake of clarity, we begin with the plastic potential function g, followed by the yield function f, noting that certain elements of the latter have evolved from the former. We then proceed to examine appropriate pairs of plastic-potential and yield functions that enable the model to reproduce the behavior of both low- and high-plasticity clays under a consistent and compatible f-g formulation.

2.1.1. Plastic Potential Function

As noted earlier, the explicit form of the plastic potential function g is generally unknown. However, the use of the stress–dilatancy relation D within the bounding-surface plasticity framework provides a practical $\mathrm{D}–\mathrm{\eta }$ relationship that is easier to determine experimentally, albeit unavoidably accompanied by a well-known paradox. In this study, we establish an explicit form of g. Inspired by the work in [33], which proposed a model applicable to both NC and OC clays and presented the plastic potential surface as a left-skewed teardrop without a closed-form equation, we adopt the closed-form expression of the original Cam-Clay (OCC) yield function [9] as the plastic potential function g. This formulation yields a left-skewed teardrop surface similar to that in [33], as given in Equation 1. The resemblance in shape provides an appropriate gradient for computations within the bounding-surface plasticity framework, and having an explicit form of g completely eliminates the paradox inherent in the theory. Consequently, the plastic potential surface g used in this study can be directly compared in size and shape with the elliptical yield surface of MCC, as illustrated in Figure 2.

$$\mathrm{G}=\ln \left({\displaystyle \frac{\mathrm{p}}{{\mathrm{p}}_0}}\right)+\left({\displaystyle \frac{\mathrm{\eta }}{\mathrm{M}}}\right)=0$$

(1)

2.1.2. Yield Function

Once, the Modified Cam Clay (MCC) model [10] was developed by improving the yield function of the original Cam Clay (OCC) model [9]. This modification was motivated by the observation that the OCC yield function did not fully reproduce the experimental behavior of normally consolidated (NC) clays, whose yield surface on the wet side should be nearly vertical to the p-axis at the point ${\mathrm{p}}_0$ (the low-stress-ratio region). Such behavior corresponds to the LPLS condition discussed earlier in this study.

However, later literature clearly reported that some NC clays do not exhibit this LPLS-type response; instead, their stress paths are inclined linearly, consistent with the yield function of the OCC model—referred to here as HPLS-type. It is therefore evident that neither the OCC nor the MCC model can fully capture the mechanical behavior of both clay types within a single constitutive framework. Accordingly, an appropriate yield function should be capable of representing both responses within one unified model. This consideration provides the fundamental motivation for developing a new unified yield function in the present study.

In constructing this unified function, the OCC yield function was adopted as the baseline because it already reproduces the HPLS-type response satisfactorily. Further modifications were then introduced so that the model could also describe the LPLS-type behavior.

The OCC yield function inherently exhibits a left-skewed shape. To enable the model to also represent the right-skewed form associated with the LPLS condition, the yield-surface shape must be controlled through a suitable parameter, here denoted as $\mathrm{\Psi }$, introduced in the $\mathrm{\eta }/\mathrm{M}$ term of the OCC yield function. Figure 3a illustrates the influence of $\mathrm{\Psi }$ on the yield surface shape, showing its transition from left- to right-skewed forms for values of $\mathrm{\Psi }$ = 1 to 3. Specifically, $\mathrm{\Psi }$ = 1 represents a left-skewed shape, $\mathrm{\Psi }$ = 2 corresponds to a right-skewed one, and intermediate values (1 < $\mathrm{\Psi }$ < 2) capture transitional clay behavior between these two extremes.

Additionally, cases with $\mathrm{\Psi }$ > 2, also included in Figure 3a, allow the model to reproduce clay behavior exhibiting high-strain softening or a distinct stress-path peak. Furthermore, it is essential to accurately control the shear strength at the critical state ($\mathrm{\eta }$ = M), corresponding to the intersection of the yield function with the Critical State Line (CSL) (e.g., [27]). Therefore, another parameter, denoted as $\mathrm{\Omega }$, was introduced in the $\mathrm{p}/{\mathrm{p}}_0$ term of the OCC yield function to precisely govern this critical-state shear strength. The corresponding equation is presented in Equation 2.

$$\mathrm{F}= \ln \left[{\left({\displaystyle \frac{\mathrm{p}}{{\mathrm{p}}_0}}\right)}^{\mathrm{\Omega }}\right]+{\left({\displaystyle \frac{\mathrm{\eta }}{\mathrm{M}}}\right)}^{\mathrm{\Psi }}= 0$$

(2)

Accordingly, Equation 2 is formulated as a revised Cam-clay-type yield function, expanded into a general logarithmic form with two shape parameters ($\mathrm{\Psi }$ and $\mathrm{\Omega }$). These parameters provide flexible control of both the wet- and dry-side curvature, enabling the surface to reproduce right- or left-skewed teardrop shapes. The conceptual basis was inspired by Hvorslev-envelope-based yield functions (e.g., Yao et al., 2009; Chen & Yang, 2017) but re-parameterized here to achieve symmetric and continuous curvature control.

It is worth noting before proceeding that the physical meaning of $\mathrm{\Psi }$ lies in controlling the skewness of the teardrop-shaped yield function, whereas $\mathrm{\Omega }$ serves as a shear-strength adjuster across different OCRs. These interpretations provide a convenient basis for subsequent referencing throughout this study and for establishing practical guidelines for selecting these two parameters.

A general form of the yield function, once plastic strain develops beyond zero after loading has commenced, incorporates a hardening parameter H, which is defined as a function of plastic strain as follows:

$$\mathrm{f} = \ln \left[{\left({\displaystyle \frac{\mathrm{p}}{{\mathrm{p}}_0}}\right)}^{\mathrm{\Omega }}\right]+{\left({\displaystyle \frac{\mathrm{\eta }}{\mathrm{M}}}\right)}^{\mathrm{\Psi }}- \mathrm{H} = 0$$

(3)

and

$$\mathrm{H} =\int \mathrm{dH}=\int {\displaystyle \frac{1}{{\mathrm{c}}_{\mathrm{p}}}}{\mathrm{d}\mathrm{\epsilon }}_{\mathrm{v}}^{\mathrm{p}}= \int {\displaystyle \frac{1+{\mathrm{e}}_0}{\mathrm{\lambda } - \mathrm{\kappa }}}{\mathrm{d}\mathrm{\epsilon }}_{\mathrm{v}}^{\mathrm{p}}$$

(4)

where ${\mathrm{c}}_{\mathrm{p}} = (\mathrm{\lambda } - \mathrm{\kappa })/(1 + {\mathrm{e}}_0)$, ${\mathrm{e}}_0$is the initial void ratio, $\mathrm{\lambda } \mathrm{and} \mathrm{\kappa }$ denote the slope of virgin compression line (VLC) and unloading-reloading line in the one-dimensional Oedometer test, respectively.

Figure 3b presents the model predictions obtained in this study using the plastic potential function g in Equation 1 together with the yield function f in Equation 2. The figure illustrates how variations in the state parameter $\mathrm{\Psi }$ lead to distinct degrees of curvature or straightness in the stress paths for NC clays. These variations in $\mathrm{\Psi }$ provide a preliminary indication that an appropriate and feasible pair of f and g can accommodate both HPLS and LPLS behaviors within a single yield function.

In the proposed model, the yield function f and the plastic potential function g are formulated to operate in tandem, thereby enhancing the predictive capability for both normally consolidated (NC) and overconsolidated (OC) clays, encompassing low- and high-plasticity soils under low stress ratios. The framework can be categorized into two representative forms, as illustrated in Figure 4. In Figure 4a, f exhibits a right-skewed curve while g adopts a left-skewed curve, a combination that is particularly suited to modeling low-plasticity clays under low stress ratios (LPLS). Conversely, Figure 4b presents a pair of left-skewed surfaces for both f and g, which proves more appropriate for capturing the behavior of high-plasticity clays at low stress levels (HPLS).

A general guideline for selecting the shape parameters $\mathrm{\Psi }$ and $\mathrm{\Omega }$ can be established by linking their roles to conventional clay indices obtained from one-dimensional oedometer tests. For NC clays exhibiting markedly curved p′–q stress paths, $\mathrm{\Psi }$ ≈ 1.2–1.4 with $\mathrm{\Omega }$ = 1.0 is recommended, whereas nearly linear responses can be reproduced using $\mathrm{\Psi }$ = $\mathrm{\Omega }$ = 1.0. Additional refinement can be achieved by adjusting $\mathrm{\Omega }$ slightly above or below unity to increase or decrease the predicted shear strength. Consistently, the calibrated $\mathrm{\Omega }$ values in Table 1 remain close to one (0.93–1.17), providing a practical starting range.

Further analyses show the practical guideline that $\mathrm{\Psi }$ and $\mathrm{\Omega }$ correlate strongly with ($\mathrm{\lambda } - \mathrm{\kappa }$), where $\mathrm{\lambda }$ and $\mathrm{\kappa }$ represent the slopes of the virgin compression line and the unloading–reloading line, respectively, in the e–ln(p) plane obtained from oedometer tests as shown in Figure 5. These empirical relationships (Equations 5–6) provide an objective basis for parameter estimation and enhance the reproducibility of the model.

$$\mathrm{\Psi }=3.28(\mathrm{\lambda }-\mathrm{\kappa })+0.83$$

(5)

$$\mathrm{\Omega }=26.82{(\mathrm{\lambda }-\mathrm{\kappa })}^2 - 6.18(\mathrm{\lambda }-\mathrm{\kappa })+1.28$$

(6)

2.2. Dilatancy Equation

An important advantage of employing stress–dilatancy D is that it enables the simulation of both volumetric (${\mathrm{d}\mathrm{\epsilon }}_{\mathrm{v}}^{\mathrm{p}}$) and deviatoric ${(\mathrm{d}\mathrm{\epsilon }}_{\mathrm{d}}^{\mathrm{p}}$) plastic strain increments, even when the plastic potential function g is not explicitly defined. In the classical bounding-surface plasticity framework, the plastic potential function is often not readily available in closed form. Accordingly, the stress–dilatancy relation—defined as the ratio of the volumetric to the deviatoric plastic strain increments—has been widely adopted in the literature, as it can be more straightforwardly identified from experimental observations. This ratio, referred to as stress–dilatancy, is expressed as shown in Equation 7.

From the theory of plasticity, the plastic strain increments are computed as the product of a proportionality constant $\mathrm{d}\mathrm{\Lambda }$ and the gradient of the plastic potential function g, such that ${\mathrm{d}\mathrm{\epsilon }}_{\mathrm{v}}^{\mathrm{p}} = \mathrm{d}\mathrm{\Lambda }\partial \mathrm{g}/\partial \mathrm{P}$ and ${\mathrm{d}\mathrm{\epsilon }}_{\mathrm{d}}^{\mathrm{p}}=\mathrm{d}\mathrm{\Lambda }\partial \mathrm{g}/\partial \mathrm{Q}$. As a result, the stress–dilatancy ratio becomes independent of the plastic multiplier $\mathrm{d}\mathrm{\Lambda }$, and can be written directly in terms of the gradients of g, leading to the following expression:

$$\mathrm{D} = {\displaystyle \frac{{\mathrm{d}\mathrm{\epsilon }}_{\mathrm{v}}^{\mathrm{p}}}{{\mathrm{d}\mathrm{\epsilon }}_{\mathrm{d}}^{\mathrm{p}}}} = {\displaystyle \frac{(\mathrm{d}\mathrm{\Lambda })\frac{\partial \mathrm{g}}{\partial \mathrm{P}}}{(\mathrm{d}\mathrm{\Lambda })\frac{\partial \mathrm{g}}{\partial \mathrm{Q}}}} = {\displaystyle \frac{\frac{\partial \mathrm{g}}{\partial \mathrm{P}}}{\frac{\partial \mathrm{g}}{\partial \mathrm{Q}}}}$$

(7)

For OC clays, the most widely adopted expression of D is the second-order dilatancy formulation of the Modified Cam-Clay (MCC) type (e.g., [29,34,35,36]). Although higher-order stress–dilatancy relations have also been proposed to capture the behavior of both normally consolidated (NC) and OC clays [30,32], closer examination reveals that, for stress ratios η not exceeding the critical state stress ratio M (approximately 1.3–1.4), a first-order or linear dilatancy relation provides satisfactory predictions and yields results comparable to the second-order formulation. Moreover, as highlighted in [30], although a high-order dilatancy expression was proposed, the parameter values adopted in practice often reduce the equation to, or make it nearly equivalent to, a linear approximation. This observation is supported by Figure 6, which presents the relationship between D and η for Fujinomori Clay and Lower Cromer till. The results confirm that, within the practical range of $\mathrm{\eta }$, a linear approximation D-$\mathrm{\eta }$ remains appropriate. In addition, opting for such a simplified linear form facilitates practical implementation while maintaining predictive adequacy.

In the proposed model, an explicit expression for g is formulated, as shown in Equation 1, thereby eliminating the paradox entirely. The expressions for $\partial \mathrm{g}/\partial \mathrm{P}$ and $\partial \mathrm{g}/\partial \mathrm{Q}$ can be derived as presented in Equations 8–9.

$${\displaystyle \frac{\partial \mathrm{g}}{\partial \mathrm{P}}}={\displaystyle \frac{1}{\mathrm{P}}}\left(1 - \left({\displaystyle \frac{\mathrm{\eta }}{\mathrm{M}}}\right)\right)$$

(8)

$${\displaystyle \frac{\partial \mathrm{g}}{\partial \mathrm{Q}}}={\displaystyle \frac{1}{\mathrm{MP}}}$$

(9)

Substituting Equation 8 and Equation 9 into Equation 7 yields:

$$\mathrm{D} = {\displaystyle \frac{\frac{\partial \mathrm{g}}{\partial \mathrm{P}}}{\frac{\partial \mathrm{g}}{\partial \mathrm{Q}}}} = \mathrm{M} - \mathrm{\eta }$$

(10)

It is evident that the expression for D in Equation 10 is linear, which is consistent with the experimental data for Fujinomori clay [37] and Lower Cromer till [38] presented in Figure 6. This confirms that the explicit derivation of the plastic potential function in this study leads to a physically consistent stress–dilatancy response observed in experiments. Accordingly, the plastic potential function is shown to be in agreement with the experimental results without introducing any paradox within the model. A limitation of the linear stress–dilatancy formulation adopted in Equation 10 is that, as shown in Figure 6, each experimental D–$\mathrm{\eta }$ curve exhibits a practical upper limit. For the clays examined in this study, the stress ratio $\mathrm{\eta }$ remains within this range, and the linear approximation is therefore appropriate. If stress ratios approach such limits in practice, higher-order dilatancy relations may be required; however, these generally preclude an explicit plastic potential function and may reintroduce the stress–dilatancy paradox.

2.3. Bounding Surface

Once the yield surface has been properly defined, the bounding surface can be established by employing the principle of a smooth transition from the OC state to the NC state. Hence, the shape of the bounding surface is identical to that of the yield surface defined in Equation 2. The bounding surface employed in this model is expressed as follow:

$$\mathrm{F} = \ln \left[{\left({\displaystyle \frac{\mathrm{P}}{{\mathrm{P}}_0}}\right)}^{\mathrm{\Omega }}\right]+{\left({\displaystyle \frac{\mathrm{\eta }}{\mathrm{M}}}\right)}^{\mathrm{\Psi }}= 0$$

(11)

where $\mathrm{\Psi }$ is the model parameter controlling the shape of the bounding surface. The distinctive feature of this formulation is that the bounding surface can take a left-skewed form ($\mathrm{\Psi }$ = 1) or a right-skewed form ($\mathrm{\Psi }$ = 2). In general, $\mathrm{\Psi }$ = 1.0 is appropriate for soils whose stress paths for NC clay on the p–q plane is nearly linear, corresponding to high plastic deformation under low stress ratios (HPLS). By contrast, values of $\mathrm{\Psi }$ in the range of approximately 1.2–1.4 are suitable for soils exhibiting low plastic deformation under low stress ratios (LPLS). The equations required for elastoplastic computations ${\displaystyle \frac{\partial \mathrm{F}}{\partial \mathrm{P}}}$, ${\displaystyle \frac{\partial \mathrm{F}}{\partial \mathrm{Q}}}$ and ${\displaystyle \frac{\partial \mathrm{F}}{\partial {\mathrm{P}}_0}}$ are presented in Equations 12-14, respectively.

$${\displaystyle \frac{\partial \mathrm{F}}{\partial \mathrm{P}}}={\displaystyle \frac{1}{\mathrm{P}}}\left(\mathrm{\Omega }-{\mathrm{\Psi }\left({\displaystyle \frac{\mathrm{\eta }}{\mathrm{M}}}\right)}^{\mathrm{\Psi }}\right)$$

(12)

$${\displaystyle \frac{\partial \mathrm{F}}{\partial \mathrm{Q}}}={\displaystyle \frac{\mathrm{\Psi }}{\mathrm{MP}}}{\left({\displaystyle \frac{\mathrm{\eta }}{\mathrm{M}}}\right)}^{\mathrm{\Psi }-1}$$

(13)

$${\displaystyle \frac{\partial \mathrm{F}}{\partial {\mathrm{P}}_0}}=-{\displaystyle \frac{\mathrm{\Omega }}{{\mathrm{P}}_0}}$$

(14)

Another fundamental aspect in modeling overconsolidated (OC) clays is the application of the mapping rule. In general, when clay initially resides in an OC state, the increasing applied load causes its stress state to gradually transition toward an NC-like condition, depending on the magnitude of loading and the soil’s stress history [31]. The mapping rule establishes a correspondence between the stress state at the current stress point—commonly referred to as the loading surface—and the stress state defined on the bounding surface. This linkage facilitates a smooth mechanical transition from the OC state to the normally consolidated (NC) state through the spacing ratio $\mathrm{R} = \mathrm{l}/\mathrm{L} = \mathrm{p}/\mathrm{P} = \mathrm{q}/\mathrm{Q}$, as illustrated in Figure 7. When $\mathrm{R}<1$, the stress point lies inside the bounding surface and plastic strains are gradually mobilized through the mapping mechanism, enabling progressive yielding in OC clays. As $\mathrm{R}\to 1$, the stress state approaches the bounding surface and the response naturally converges to the NC behavior governed by the hardening law. Typically, both surfaces are assumed to share a similar shape, with the current loading surface being smaller in size and nested within the bounding surface. In practical terms, the spacing ratio $\mathrm{R}$ controls the degree of plasticity, stiffness degradation, and peak strength development, which is particularly important for soils exhibiting complex stress–strain behavior.

2.4. Hardening Law and Plastic Modulus

Although, for the OC state, the stress of clay lies on a portion of the unloading–reloading line on e-ln(p) plane with a slope equal to κ, where conventional calculations give only purely elastic strains though a constant $\mathrm{\kappa }/(1 + {\mathrm{e}}_0).$ However, in the bounding surface plasticity framework, it is common to introduce an alternative approach to enable OC clays to develop plastic strains in addition to elastic strains. Specifically, it assumes that the stress point lies on the isotropic compression line which has a slope of $\mathrm{\lambda }$. This conceptual framework not only allows the model to capture plastic strains in OC clays but also provides a rational basis for employing the parameter ${\mathrm{c}}_{\mathrm{p}}$in the hardening law. Accordingly, the isotropic hardening law can be consistently applied to both the NC and OC states within a unified framework. The hardening law employed in this model is expressed as follows:

$${\displaystyle \frac{{\mathrm{d}\mathrm{\epsilon }}_{\mathrm{v}}^{\mathrm{p}}}{\mathrm{d}{\mathrm{P}}_0}} = {\displaystyle \frac{{\mathrm{c}}_{\mathrm{p}}}{{\mathrm{P}}_0}} = {\displaystyle \frac{1}{{\mathrm{P}}_0}}{\displaystyle \frac{\mathrm{\lambda } - \mathrm{\kappa }}{1 + {\mathrm{e}}_0}}$$

(15)

where${\mathrm{d}\mathrm{\epsilon }}_{\mathrm{v}}^{\mathrm{p}}$ is the volumetric plastic strain increment. From Equation 15, for a single sub-step computation, $\mathrm{d}{\mathrm{P}}_0$ can be determined as ${\mathrm{P}}_0 = \left({\displaystyle \frac{{\mathrm{P}}_0}{{\mathrm{c}}_{\mathrm{p}}}}\right){\mathrm{d}\mathrm{\epsilon }}_{\mathrm{v}}^{\mathrm{p}}$ For the plastic modulus ${\mathrm{K}}_{\mathrm{p}}$ which is employed together with the proportional constant $\mathrm{d}\mathrm{\Lambda }$ to determine the change of the hardening parameter in consistency condition which expressed as $\mathrm{dH} = \mathrm{d}\mathrm{\Lambda }·{\mathrm{K}}_{\mathrm{p}}$. By substituting this relation into Equation 3 and rearranging, the expression for ${\mathrm{K}}_{\mathrm{p}}$ can be obtained as presented in Equation 16.

$${\mathrm{K}}_{\mathrm{p}}={\displaystyle \frac{1}{{\mathrm{c}}_{\mathrm{p}}}}{\displaystyle \frac{\partial \mathrm{g}}{\partial \mathrm{P}}}={\displaystyle \frac{1}{{\mathrm{c}}_{\mathrm{p}}}}\mathrm{D}{\displaystyle \frac{\partial \mathrm{g}}{\partial \mathrm{Q}}}={\displaystyle \frac{1+{\mathrm{e}}_0}{\mathrm{\lambda } - \mathrm{\kappa }}}\left(\mathrm{M} - \mathrm{\eta }\right){\displaystyle \frac{\partial \mathrm{g}}{\partial \mathrm{Q}}}$$

(16)

Upon closer examination Equation 16, using stress point (P, Q) is highly restricted to OCR = 1, because P and Q are located on the bounding surface. To allow ${\mathrm{K}}_{\mathrm{p}}$ valid for OCR > 1—i.e., for stress states on the loading surface that is nested within the bounding surface—it is therefore necessary to revise ${\mathrm{K}}_{\mathrm{p}}$ so it remains applicable to both normally and overconsolidated clays.

A key mechanism enabling plastic strain for stress points located inside the bounding surface via ${\mathrm{K}}_{\mathrm{p}}$ is that ${\mathrm{K}}_{\mathrm{p}}$must be a positive number at the critical state. However, as given by Equation 16, ${\mathrm{K}}_{\mathrm{p}}$ vanishes at the critical state which is a condition appropriate for normally consolidated clay but unsuitable for overconsolidated clay. To address this, Dafalias and Herrmann (1986) [39] proposed the concept of a virtual peak stress ratio, ${\mathrm{M}}_{\mathrm{k}}$, to replace M in Equation 16, thereby providing a smoother and more accurate description of OC behavior. The expression for ${\mathrm{M}}_{\mathrm{k}}$ is shown in Equation 17.

$${\mathrm{M}}_{\mathrm{k}}=\mathrm{M}{\mathrm{R}}^{-\mathrm{n}}$$

(17)

where M is the stress ratio at the critical state, R is the spacing ratio, and n is a positive parameter. Since the spacing ratio R equals 1 when the stress state lies on the bounding surface, ${\mathrm{M}}_{\mathrm{k}}$becomes identical to M, but takes values greater than M when the stress point lies inside the bounding surface. Consequently, for cases with OCR > 1, the stress–strain response may exhibit a peak at a stress ratio exceeding M.

The parameter n in Equation 17 is introduced to adjust the virtual stress ratio ${\mathrm{M}}_{\mathrm{k}}$​​, and consequently the plastic modulus ${\mathrm{K}}_{\mathrm{p}}$​. When OCR = 1, the parameter n has no effect because the spacing ratio R​ equals unity. Although many models employ a constant value for n. However, for OCR > 1, a value of n should lead to a different value of both ${\mathrm{M}}_{\mathrm{k}}$​​ and ${\mathrm{K}}_{\mathrm{p}}$​. Intuitively, this suggests that n may not be a constant but rather a quantity related to the soil’s overconsolidation history, since its influence on model responses varies with OCR. In [29] a power-law regression for n-OCR has been proposed. Nevertheless, this equation may be compatible only with specific constitutive models, depending on the mapping rule adopted. In this study, a new n–OCR relationship is proposed as presented in Equation 18.

$$\mathrm{n}=0.5+{\displaystyle \frac{2}{\mathrm{OCR}}}$$

(18)

Another crucial issue in considering the parameter n is whether it should take the same value for both undrained and drained loading conditions. It is well known that the drained shear strength is generally higher than the undrained shear strength due to the pore pressure that develops during undrained loading. However, the difference between drained and undrained shear strength ($\mathrm{\Delta }$) decreases with increasing OCR. Figure 8 depicts the typical paired stress paths for drained and undrained conditions across various OCR values, clearly demonstrating that the shear strength difference is OCR-dependent. Moreover, Figure 8 also reveals that the critical state stress ratio M exerts an additional influence on $\mathrm{\Delta }$; specifically, larger values of M correspond to greater values of $\mathrm{\Delta }$. Hence, the relationship can be expressed as $\mathrm{\Delta }$ = f(OCR, M).

In this study, the relationship among $\mathrm{\Delta }$, OCR, and M is presented in Figure 4 and expressed in Equation 19. From Equation 19, it can be observed that $\mathrm{\Delta }$ is more precisely defined as $\mathrm{\Delta }$ = f(OCR, ${\mathrm{M}}^2$). Accordingly, it may be inferred that the parameter ${\mathrm{M}}^2$ plays a key role in differentiating the value of n between drained and undrained cases. The specific expression for n under drained loading is provided in Equation 20.

$$\mathrm{\Delta }={ \mathrm{M}}^2\left[0.013\left(\mathrm{OCR}\right)+{\displaystyle \frac{1.2}{\mathrm{OCR}}} - 0.26\right]$$

(19)

Figure 9. Comparison of the disparity between drained and undrained shear strengths at the critical state for various OCR values.

Figure 9. Comparison of the disparity between drained and undrained shear strengths at the critical state for various OCR values.

$$\mathrm{n}=0.5+{\mathrm{M}}^2+{\displaystyle \frac{2}{\mathrm{OCR}}}\left(\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{D}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{d}\mathrm{L}\mathrm{o}\mathrm{a}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{g}\right)$$

(20)

3. Elastoplastic Constitutive Relations

The constitutive models employing the stress p and q implicitly adopt the Extended Mises failure criterion, which forms a circular failure envelope on the $\mathrm{\pi }$-plane and typically leads to an overestimation of shear strength. This limitation can be eliminated by employing the SMP failure criterion [13] as shown in Figure 10. The present model also adopts this approach by employing a transformed stress formulation to integrate the SMP criterion to the proposed model. A brief description of the transformed stress and the SMP criterion is provided prior to presenting the governing equations of the proposed model.

3.1. Transformed Stress Tensor ${\stackrel{\mathrm{~}}{\mathrm{\sigma }}}_{\mathrm{ij}}$

A transformed stress tensor ${\stackrel{\mathrm{~}}{\mathrm{\sigma }}}_{\mathrm{ij}}$ is employed instead of the ordinary stress tensor ${\mathrm{\sigma }}_{\mathrm{ij}}$ in the analyses to improve shear strength prediction by introducing the SMP criterion [13]. Once the ${\mathrm{\sigma }}_{\mathrm{ij}}$ is given, the ${\stackrel{\mathrm{~}}{\mathrm{\sigma }}}_{\mathrm{ij}}$ can be calculated as follows:

$${\stackrel{\mathrm{~}}{\mathrm{\sigma }}}_{\mathrm{ij}} = \stackrel{\mathrm{~}}{\mathrm{p}}{\mathrm{\delta }}_{\mathrm{ij}} + {\stackrel{\mathrm{~}}{\mathrm{s}}}_{\mathrm{ij}} = \mathrm{p}{\mathrm{\delta }}_{\mathrm{ij}} + {\displaystyle \frac{{\mathrm{q}}_{\mathrm{SMP}}}{{\mathrm{q}}_{\mathrm{c}}}}{\mathrm{s}}_{\mathrm{ij}}$$

(21)

where ${\mathrm{\delta }}_{\mathrm{ij}}$ is the Kronecker’s delta, $\mathrm{p}$ and ${\mathrm{s}}_{\mathrm{ij}}$ are the mean stress and deviatoric stress tensor calculated from the ordinary stress tensor, respectively, $\stackrel{\mathrm{~}}{\mathrm{p}}$ and ${\stackrel{\mathrm{~}}{\mathrm{s}}}_{\mathrm{ij}}$ are the mean stress and deviatoric stress tensor calculated from the transformed stress tensor, respectively, which can be written as

$$\stackrel{\mathrm{~}}{\mathrm{p}}={\displaystyle \frac{1}{3}}{\stackrel{\mathrm{~}}{\mathrm{\sigma }}}_{\mathrm{ii}}$$

(22)

$${\stackrel{\mathrm{~}}{\mathrm{s}}}_{\mathrm{ij}}={\stackrel{\mathrm{~}}{\mathrm{\sigma }}}_{\mathrm{ij}}-\stackrel{\mathrm{~}}{\mathrm{p}}{\mathrm{\delta }}_{\mathrm{ij}}$$

(23)

${\mathrm{q}}_{\mathrm{SMP}}$ and ${\mathrm{q}}_{\mathrm{c}}$ are the radius of the circular cone and the Extended Mises cone which can be calculated using

$${\mathrm{q}}_{\mathrm{SMP}}=\sqrt{{\displaystyle \frac{2}{3}}}{\displaystyle \frac{2{\mathrm{I}}_1}{3\sqrt{\left({\mathrm{I}}_1{\mathrm{I}}_2-{\mathrm{I}}_3\right)\left({\mathrm{I}}_1{\mathrm{I}}_2-9{\mathrm{I}}_3\right)-1}}}$$

(24)

$${\mathrm{q}}_{\mathrm{c}}= \sqrt{{\mathrm{s}}_{\mathrm{kl}}{\mathrm{s}}_{\mathrm{kl}}}$$

(25)

where ${\mathrm{I}}_1$, ${\mathrm{I}}_2$ and ${\mathrm{I}}_3$ are the first, second, and third stress invariants, respectively.

3.2. Stress–Strain Modelling for Elastoplastic Constitutive Model

The total strain increment $\mathrm{d}\mathrm{\epsilon }$is decomposed into elastic $\mathrm{d}{\mathrm{\epsilon }}^{\mathrm{e}}$ and plastic $\mathrm{d}{\mathrm{\epsilon }}^{\mathrm{p}}$components as follows:

$$\mathrm{d}\mathrm{\epsilon } = \mathrm{d}{\mathrm{\epsilon }}^{\mathrm{e}} + \mathrm{d}{\mathrm{\epsilon }}^{\mathrm{p}} = \left(\mathrm{d}{\mathrm{\epsilon }}_{\mathrm{v}}^{\mathrm{e}} + \mathrm{d}{\mathrm{\epsilon }}_{\mathrm{d}}^{\mathrm{e}}\right) + \left(\mathrm{d}{\mathrm{\epsilon }}_{\mathrm{v}}^{\mathrm{p}} + \mathrm{d}{\mathrm{\epsilon }}_{\mathrm{d}}^{\mathrm{p}}\right)$$

(26)

In the calculation of $\mathrm{d}\mathrm{\Lambda }$, it can be computed by applying the consistency condition (dF = 0) conjunction with application of chain rule, as shown in Equation 27.

$${\displaystyle \frac{\partial \mathrm{F}}{\partial \mathrm{P}}}\mathrm{dP}+{\displaystyle \frac{\partial \mathrm{F}}{\partial \mathrm{Q}}}\mathrm{dQ}+{\displaystyle \frac{\partial \mathrm{F}}{\partial \mathrm{H}}}\mathrm{dH}={\displaystyle \frac{\partial \mathrm{F}}{\partial \mathrm{P}}}\mathrm{dP}+{\displaystyle \frac{\partial \mathrm{F}}{\partial \mathrm{Q}}}\mathrm{dQ} - \mathrm{d}\mathrm{\Lambda }{\mathrm{K}}_{\mathrm{p}}=0 \left(\because \partial \mathrm{F}/\partial \mathrm{H} =-1\right)$$

(27)

Thus,

$$\mathrm{d}\mathrm{\Lambda }={\displaystyle \frac{\frac{\partial \mathrm{F}}{\partial \mathrm{P}}\mathrm{dP}+\frac{\partial \mathrm{F}}{\partial \mathrm{Q}}\mathrm{dQ}}{{\mathrm{K}}_{\mathrm{p}}}}$$

(28)

It is worth noting that $\mathrm{d}\mathrm{\Lambda }$in Equation 28 is highly beneficial for conducting drained analyses, since the pore pressure is zero and the input data can therefore be directly expressed in terms of the stress increments dP and dQ. In contrast, for undrained analyses, the input for the computation (especially in the numerical analyses) must be given in terms of strain increments with the condition of $\mathrm{d}{\mathrm{\epsilon }}_{\mathrm{v}} = 0$. The corresponding calculations can then be carried out as follows.

$$\left\{\begin{array}{c}\mathrm{dp}\\ \mathrm{dq}\end{array}\right\}={\mathrm{D}}_{2\times 2}^{\mathrm{e}}\left\{\begin{array}{c}\mathrm{d}{\mathrm{\epsilon }}_{\mathrm{v}}^{\mathrm{e}}\\ \mathrm{d}{\mathrm{\epsilon }}_{\mathrm{d}}^{\mathrm{e}}\end{array}\right\}$$

(29)

Where

$${\mathrm{D}}_{2\times 2}^{\mathrm{e}}=\left[\begin{array}{cc}\mathrm{K}& 0\\ 0& 3\mathrm{G}\end{array}\right]$$

(30)

and the elastic bulk and shear modulus can be expressed in term of the initial void ratio ${\mathrm{e}}_0$, the Poisson’s ratio $\mathrm{\upsilon }$ and the slope of unloading-reloading line in e-ln(p) plane ($\mathrm{\kappa })$ as

$$\mathrm{K}= {\displaystyle \frac{1+{\mathrm{e}}_0}{\mathrm{\kappa }}}\mathrm{p}$$

(31)

$$\mathrm{G}=\mathrm{K}{\displaystyle \frac{3\left(1 - 2\mathrm{\upsilon }\right)}{2\left(1+\mathrm{\upsilon }\right)}}$$

(32)

The elastic and plastic strain increment can be written as follows:

$$\left\{\begin{array}{c}\mathrm{d}{\mathrm{\epsilon }}_{\mathrm{v}}^{\mathrm{e}}\\ \mathrm{d}{\mathrm{\epsilon }}_{\mathrm{d}}^{\mathrm{e}}\end{array}\right\} = \left\{\begin{array}{c}\mathrm{d}{\mathrm{\epsilon }}_{\mathrm{v}}\\ \mathrm{d}{\mathrm{\epsilon }}_{\mathrm{d}}\end{array}\right\} - \left\{\begin{array}{c}\mathrm{d}{\mathrm{\epsilon }}_{\mathrm{v}}^{\mathrm{p}}\\ \mathrm{d}{\mathrm{\epsilon }}_{\mathrm{d}}^{\mathrm{p}}\end{array}\right\}$$

(33)

$$\left\{\begin{array}{c}\mathrm{D}{\mathrm{\epsilon }}_{\mathrm{v}}^{\mathrm{p}}\\ \mathrm{d}{\mathrm{\epsilon }}_{\mathrm{d}}^{\mathrm{p}}\end{array}\right\}= 〈\mathrm{d}\mathrm{\Lambda }〉{\displaystyle \frac{\partial \mathrm{G}}{\partial \mathrm{Q}}}\left\{\begin{array}{c}\mathrm{D}\\ 1\end{array}\right\}$$

(34)

thus,

$$\mathrm{d}\mathrm{\Lambda } ={\displaystyle \frac{\mathrm{K}\frac{\partial \mathrm{F}}{\partial \mathrm{P}}\mathrm{d}{\mathrm{\epsilon }}_{\mathrm{v}}+3\mathrm{G}\frac{\partial \mathrm{F}}{\partial \mathrm{Q}}\mathrm{d}{\mathrm{\epsilon }}_{\mathrm{q}}}{{\mathrm{K}}_{\mathrm{p}}+\mathrm{KD}\frac{\partial \mathrm{F}}{\partial \mathrm{P}}\frac{\partial \mathrm{g}}{\partial \mathrm{Q}}+3\mathrm{G}\frac{\partial \mathrm{F}}{\partial \mathrm{Q}}\frac{\partial \mathrm{g}}{\partial \mathrm{Q}}}}$$

(35)

$${\mathrm{D}}_{2\times 2}^{\mathrm{ep}}=\left[\begin{array}{cc}\mathrm{K}& 0\\ 0& 3\mathrm{G}\end{array}\right] - {\displaystyle \frac{\mathrm{h}(\mathrm{\Lambda })}{{\mathrm{K}}_{\mathrm{p}}+\mathrm{KD}\frac{\partial \mathrm{F}}{\partial \mathrm{P}}\frac{\partial \mathrm{g}}{\partial \mathrm{Q}}+3\mathrm{G}\frac{\partial \mathrm{F}}{\partial \mathrm{Q}}\frac{\partial \mathrm{g}}{\partial \mathrm{Q}}}}\left[\begin{array}{cc}{\mathrm{K}}^2\mathrm{D}{\displaystyle \frac{\partial \mathrm{F}}{\partial \mathrm{P}}}{\displaystyle \frac{\partial \mathrm{g}}{\partial \mathrm{Q}}}& 3\mathrm{KGD}{\displaystyle \frac{\partial \mathrm{F}}{\partial \mathrm{Q}}}{\displaystyle \frac{\partial \mathrm{g}}{\partial \mathrm{Q}}}\\ 3\mathrm{KG}{\displaystyle \frac{\partial \mathrm{F}}{\partial \mathrm{P}}}{\displaystyle \frac{\partial \mathrm{g}}{\partial \mathrm{Q}}}& 9{\mathrm{G}}^2{\displaystyle \frac{\partial \mathrm{F}}{\partial \mathrm{Q}}}{\displaystyle \frac{\partial \mathrm{g}}{\partial \mathrm{Q}}}\end{array}\right]$$

(36)

where h($\mathrm{\Lambda }$) = 1 for $\mathrm{d}\mathrm{\Lambda }$ > 0, and h($\mathrm{\Lambda }$) = 0 when $\mathrm{d}\mathrm{\Lambda }$ < 0. The stress-strain relationship can be calculated as

$$\left\{\begin{array}{c}\mathrm{dp}\\ \mathrm{dq}\end{array}\right\}={ \mathrm{D}}_{2\times 2}^{\mathrm{ep}}\left\{\begin{array}{c}\mathrm{d}{\mathrm{\epsilon }}_{\mathrm{v}}\\ \mathrm{d}{\mathrm{\epsilon }}_{\mathrm{d}}\end{array}\right\}$$

(37)

Before presenting a comparison between the analytical results of the proposed model and the experimental data in the following section, the final part of this section illustrates the variation of the model’s predictions, with particular emphasis on the influence of the parameter $\mathrm{\Psi }$. This parameter constitutes a key factor enabling the model to reproduce clay behavior as either HPLS or LPLS across different OCR values. Figure 11 depicts the normalized stress paths (p′/${\mathrm{p}\text{'}}_0$ – q/${\mathrm{p}\text{'}}_0$) for OCR values of 1, 2, 6, and 20 conjunctions with $\mathrm{\Psi }$ in the range of 1.0–3.0. It can be observed that Ψ effectively captures the characteristic responses of HPLS and LPLS for normally consolidated clays (OCR = 1) as well as for slightly to moderately overconsolidated clays (OCR ≈ 2). However, in the case of heavily overconsolidated clays, the results converge, showing little distinction across different $\mathrm{\Psi }$ values. This outcome is consistent with the aforementioned aim of this study.

4. Model Verifications and Discussions

In this section, the results obtained from the proposed model are presented for both drained and undrained analyses. The results are compared against experimental data for four clays i.e., Boston Blue clay (OCR = 1, 2, 4 and 8), Lower Cromer till (OCR = 1, 2, 4 and 10), Kaolin clay (OCR = 1, 1.2, 1.45, 1.8, 2.5, 3, 4.5 and 6.5) and London clay (OCR = 1, 2, 6 and 20). Furthermore, for each of these clays, the experimental results are systematically compared with predictions from three constitutive models: (I) the proposed model, (II) the model introduced by Xu et al., 2024 [30], and (III) the model developed by Chen and Yang, 2017 [29]. The soil parameters used in the analyses for 3 models mentioned above are summarized in Table 1.

To ensure the reliability of the proposed formulation, model validation was performed against well-documented triaxial test data under both drained and undrained compression and extension conditions. Because any constitutive model represents the stress–strain behavior of soil, it is rational to assess its accuracy through characteristic stress paths (e.g., $\mathrm{p}–\mathrm{q}$ or normalized $\mathrm{p}/{\mathrm{p}}_0– \mathrm{q}/{\mathrm{p}}_0$) and strain-based relationships such as $\mathrm{q}/{\mathrm{p}}_0–{ \mathrm{\epsilon }}_1$and ${\mathrm{\epsilon }}_{\mathrm{v}}–{\mathrm{\epsilon }}_{\mathrm{d}}$, where ${\mathrm{\epsilon }}_1$, ${\mathrm{\epsilon }}_{\mathrm{v}}$, and ${\mathrm{\epsilon }}_{\mathrm{d}}$denote the axial, volumetric, and deviatoric strains, respectively. All these quantities are typically obtained from triaxial compression and extension tests, and the corresponding comparisons between calculated and experimental results are presented in Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18, Figure 19 and Figure 20.

4.1. Undrained Analysis

Figure 12, Figure 13, Figure 14 and Figure 15 present comparisons between the triaxial compression test results and the predictions of the three aforementioned models, shown separately for each model to highlight their differences. In particular, the proposed model demonstrates superior capability in modeling the high plasticity under low stress ratio (HPLS) behavior compared with the other models. The corresponding comparison for the undrained triaxial extension test is provided in Figure 16 and Figure 17.

It is clearly evident from the triaxial compression test results of London clay and Lower Cromer Till, shown in Figure 12 and Figure 13, respectively. Both London clay and Lower Cromer till can be classified as exhibiting high plasticity under low stress levels (HPLS), as the stress paths for OCR = 1 display a pronounced inclination from the very onset of loading—a characteristic not commonly observed. Nevertheless, the proposed model provides predictions that align closely with the experimental results. By contrast, in both figures, the other two models of Xu et al. (2024) and Chen and Yang (2017) fail to capture this behavior adequately. Their predicted stress paths instead correspond to low plasticity under low stress levels (LPLS). In addition, it is observed that for both soils, all the models yield comparable predictions when OCR > 2.

For the experimental results of Boston Blue Clay and Kaolin Clay, shown in Figure 14 and Figure 15, respectively, the stress paths at OCR = 1 exhibit a pronounced curvature consistent with the bullet/teardrop-shaped yield surface in bounding surface models. This behavior is characteristic of clays classified as having low plasticity under low stress ratios (LPLS). In both figures, all three models successfully reproduce stress paths consistent with LPLS behavior. Among them, the proposed model and the model of Xu et al. 2024) demonstrate relatively high accuracy across all OCR values, whereas the model developed by Chen and Yang (2017) provides reasonable accuracy only for OCR > 1, while at OCR = 1 it tends to overestimate the response.

The ability of the proposed model to reproduce both the nearly linear HPLS-type and the curved LPLS-type stress paths directly reflects the combined roles of the yield surface parameter Ψ and the plastic potential function g introduced in Section 2.1, as illustrated earlier in Figure 2, Figure 3 and Figure 4.

For the undrained triaxial extension tests on Lower Cromer Till shown in Figure 16, the simulated stress paths are in good agreement with the experimental results, and the model also reproduces shear strengths comparable to the tests at OCR = 1, 2, 4, and 10. This is attributed to the use of the SMP failure criterion, which accurately captures the shear strength envelope, thereby reflecting the well-known fact that the undrained triaxial extension shear strength is lower than that under undrained triaxial compression. It should be noted that Lower Cromer till under undrained loading exhibits high-plasticity behavior at low stress ratios (HPLS) in both compression and extension, underscoring the existence of clays that truly fall into the HPLS category. By contrast, for the undrained triaxial extension tests on Kaolin clay (Figure 17), the stress path exhibits low-plasticity behavior at low stress ratios (LPLS). In this case, the model underestimates the shear strength at OCR = 1 and 1.2 but provides close agreement with the experimental data at OCR = 6 and 10.

4.2. Drained Analysis

To further evaluate the predictive capability of the proposed model, three drained triaxial compression test results were selected for comparison. In addition to the Lower Cromer till, the parameters and test results of two other clays—Shanghai clay and Fujinomori clay—were also considered in the investigation, as summarized in Table 2. Figure 18, Figure 19 and Figure 20 present the test results and corresponding model predictions for Lower Cromer till (OCR = 1, 1.25, 1.5, 2.0, 4.0, and 10), Fujinomori clay (OCR = 1, 2, 4, and 8), and Shanghai clay (OCR = 1, 2, 4, 8, and 12).

The comparisons in Figure 18, Figure 19 and Figure 20 show a close agreement between the experimental data and model predictions. For Lower Cromer till, the normalized stress–strain curves reveal contractive behavior at OCR = 1–4, characterized by curved stress paths without a distinct peak, along with residual shear strengths that vary systematically with OCR. In contrast, Shanghai clay exhibits contractive behavior without a peak at OCR = 1–4, but transitions to a contractive response followed by dilation at OCR = 8 and 10. Notably, for all OCR values, Shanghai clay develops essentially the same residual shear strength. Fujinomori clay demonstrates responses similar to Shanghai clay: contractive behavior without a peak at OCR = 1–4, followed by contractive behavior and subsequent dilation at OCR = 8, with residual shear strength remaining nearly constant across all OCR values. The proposed model successfully captures this behavior in close agreement with the experimental results.

From a practical geotechnical engineering perspective, the proposed model can be readily implemented in numerical analyses for problems involving soft to stiff clays under low stress ratios, such as shallow foundations, staged embankment construction, and consolidation-induced shear. The unified description of both normally consolidated and overconsolidated behaviors within a single framework allows engineers to simulate loading–unloading–reloading processes without switching constitutive models. In particular, the ability to reproduce both high- and low-plasticity clay responses provides a practical advantage for site conditions where clay plasticity varies spatially. These features make the model suitable for predicting stress–strain and volumetric responses in complex field conditions while maintaining a relatively simple parameter calibration procedure.

5. Conclusions

This study presents a unified bounding-surface plasticity formulation that reproduces both high-plasticity (HPLS) and low-plasticity (LPLS) clay responses under low stress ratios across a broad OCR range. The model combines a teardrop-shaped yield surface controlled by the shape parameters $\mathrm{\Psi }$ and $\mathrm{\Omega }$, an explicit plastic-potential $g$that resolves the dilatancy paradox, a practical linear stress–dilatancy relation $\mathrm{D} = \mathrm{M}-\mathrm{\eta }$, and an OCR-sensitive virtual peak ratio. Together, these ingredients produce a coherent, implementation-friendly framework with smooth NC↔OC transitions and no ad-hoc model switching.

Although two additional shape parameters are introduced, they can be directly determined from standard clay indices ($\mathrm{\lambda }$ and $\mathrm{\kappa }$), preserving the parameter parsimony typical of MCC while enhancing predictive capability. Validation against undrained and drained tests shows that the model captures the steep, near-linear NC paths characteristic of HPLS clays and the curved trajectories of LPLS clays, as well as correct contractive behavior at low OCR and contractive-to-dilative transitions at higher OCR—without inflating the parameter set or requiring higher-order dilatancy forms.

From an engineering perspective, the formulation offers a single, interpretable constitutive model that is straightforward to calibrate from standard triaxial and oedometer data, preserves numerical robustness, and is readily implementable in finite-element analyses for practical geotechnical design, enabling consistent prediction of clay behavior under coupled loading, unloading, and consolidation processes within a unified constitutive framework.