A Unified Bounding Surface Plasticity Model for Both High- and Low-Plasticity Clays Under Low Stress Ratios

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 models originally developed to describe the behavior of normally consolidated (NC) clays based on the critical state soil mechanic (CSSM). It is evident that in subsequent research, numerous modifications of the MCC have been proposed to enhance its applicability to geomaterials [e.g., 11–22]. For OC clays, the yield function (f) of the Modified Cam Clay model often overestimates the peak strength of heavily overconsolidated materials on the dry side of the critical state [23,24,25,26,27]. Therefore, it is commonly revised to better capture the OC behavior. A key feature in reformulations is the incorporation of the Hvorslev envelope [28] into the yield function, resulting the familiar bullet- or teardrop-shaped yield surfaces. Embedding the Hvorslev envelope within the yield function of MCC not only improves the accuracy of peak-strength predictions for moderated to heavily overconsolidated clay but also ensures smooth variation of the surface gradient across the full range of 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 Figures 1a–b) and performs well for clays that exhibit relatively low plastic deformation under low stress ratios, η = q/p (hereafter referred to as LPLS). 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 Figures 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., Ref#) with the adaptation of a non-associated flow rule, in which the yield function differs from the plastic potential function (f ≠ 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 D= (∂g/∂p) / (∂g/∂q). It should be noted that, when applying D in constitutive formulations, a paradox arises because D is often assumed to equal (∂f/∂p)/(∂f/∂q) instead of (∂g/∂p)/(∂g/∂q) [e.g., 29-31]. 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 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 paths. To this end, two key aspects are introduced: (1) the formulation of an appropriate yield function that facilitates such modeling, and (2) the adoption of a suitable plastic potential function that eliminates the paradox from the model formulation.

2. Formulations of the Constitutive Model

2.1. Yield and Plastic Potential Functions

As mentioned earlier, the yield surface in the MCC model is known to overestimate the shear strength on the dry side, which reduces the predictive accuracy of the model. To correct this limitation, the bullet/teardrop-shaped yield surfaces are often introduced. In this study, the proposed model to be applicable to clays exhibiting both low and high plasticity under low stress ratios η, the yield function must be sufficiently versatile. Specifically, under low stress ratios, it should allow a right-skewed teardrop form which is suitable for low-plasticity overconsolidated (OC) clays, as well as a left-skewed teardrop form which is more appropriate for high-plasticity OC clays. The proposed yield function is presented in Equation 1 which contain additional model parameter ($\mathsf{\Psi }$), and effect of the corresponding variation of parameter $\mathsf{\Psi }$ to yield function is depicted in Figure 2.

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

(1)

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)}^{\mathsf{\Omega }}\right]+{\left({\displaystyle \frac{\mathsf{\eta }}{\mathrm{M}}}\right)}^{\mathsf{\Psi }}- \mathrm{H} = 0$$

(2)

and

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

(3)

As shown in Figure 2, a left-skewed teardrop form can be reproduced with $\mathsf{\Psi }$ =1.0, whereas a right-skewed teardrop form is obtained with $\mathsf{\Psi }$ =2.0. In practice, the parameter $\mathsf{\Psi }$ may be varied within the range of 1.0–2.0 to best capture the stress path in the p–q space for both normally consolidated (NC) and overconsolidated (OC) clays. It should also be noted that values of $\mathsf{\Psi }$ greater than 2.0 are permissible; such values accentuate the height of the teardrop shape on the right-hand side, which is consistent with the behavior of lower-plasticity clays under low stress ratios.

In the study of [32], the shape of the plastic potential surface g was found to resemble the left-skewed teardrop form of without a specific equation. Motivated by this observation, this study adopted a closed-form expression of the yield function of the original Cam clay (OCC) [9] as the plastic potential function g, represented as a left-skewed teardrop similar to [32], as expressed in Equation 4.

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

(4)

Figure 3. Plastic potential function g of the proposed model compared with that of the MCC [10] model.

Figure 3. Plastic potential function g of the proposed model compared with that of the MCC [10] model.

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

2.2. Dilatancy Equation

An important advantage of employing stress dilatancy is that it enables the simulation of both volumetric and deviatoric plastic strains, even when the plastic potential function is not explicitly defined. In the classical bounding-surface plasticity framework, the plastic potential function g is not readily obtained in closed form. Accordingly, the literature commonly adopts the stress-dilatancy relation—defined as the ratio of the volumetric to the deviatoric plastic strain increments—which can be identified more straightforwardly from experiments. This ratio, referred to as stress dilatancy, is typically expressed as shown in Equation 5.

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

(5)

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,33-35]. Although higher-order stress dilatancy relations have also been proposed to capture the behavior of both normally consolidated (NC) and OC clays [30-31], 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 5, which presents the relationship between D and η for Fujinomori Clay and Lower Cromer till. The results confirm that, within the practical range of η, a linear approximation D-η 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 4, thereby eliminating the paradox entirely. The expressions for ∂g/∂P and ∂g/∂Q can be derived as presented in Equations 6–7.

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

(6)

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

(7)

Substituting Equation 6 and Equation 7 into Equation 5 yields:

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

(8)

It is evident that the expression for D in Equation 8 is linear, which is consistent with the experimental data for Fujinomori clay [36] and Lower Chrome till [37] presented in Figure 5. Accordingly, the plastic potential function is shown to be in agreement with the experimental results without introducing any paradox within the model.

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 1. The bounding surface employed in this model is expressed as follow:

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

(9)

where $\mathsf{\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 ($\mathsf{\Psi }$ =1) or a right-skewed form ($\mathsf{\Psi }$ =2). In general, $\mathsf{\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 $\mathsf{\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 F}{\partial P}}$, ${\displaystyle \frac{\partial F}{\partial Q}}$ and ${\displaystyle \frac{\partial F}{\partial P_0}}$ are presented in Equations 10-12, respectively.

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

(10)

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

(11)

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

(12)

Another fundamental aspect in modeling overconsolidated (OC) clays is the application of the mapping rule. In general, the mapping rule establishes a correspondence between the stress state located 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 transition from the OC state to the normally consolidated (NC) state through the spacing ratio $\mathrm{R}=\mathrm{l}/𝓁=\mathrm{p}/\mathrm{P}=\mathrm{q}/\mathrm{Q}$, as illustrated in Figure 6. 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.

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-log(p) plane with a slope equal to κ, where conventional calculations give only purely elastic strains though a constant κ/(1+e0). 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 λ. 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 Cp=(λ–κ)/(1+e0) 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}\mathsf{\epsilon }}_{\mathrm{v}}^{\mathrm{p}}}{\mathrm{d}{\mathrm{P}}_0}}={\displaystyle \frac{\mathrm{Cp}}{{\mathrm{P}}_0}}={\displaystyle \frac{1}{{\mathrm{P}}_0}}{\displaystyle \frac{\mathsf{\lambda }-\mathsf{\kappa }}{1+\mathrm{e}0}}$$

(13)

where e0 and${d\epsilon }_v^p$ are the initial void ratio, volumetric plastic strain increment, respectively. From Equation 13, for a single sub-step computation, $d{P}_0$ can be determined as ${d{P}_0=\left({\displaystyle \frac{P_0}{\mathrm{C}\mathrm{p}}}\right)d\epsilon }_v^p.$For the plastic modulus Kp which is employed together with the proportional constant $\mathrm{d}\Lambda$ to determine the change of the hardening parameter in consistency condition which expressed as $\mathrm{dH} = \mathrm{d}\Lambda \cdot \mathrm{Kp}$. By substituting this relation into Equation 3 and rearranging, the expression for Kp can be obtained as presented in Equation 14.

$$\mathrm{Kp}={\displaystyle \frac{1}{\mathrm{Cp}}}{\displaystyle \frac{\partial \mathrm{g}}{\partial \mathrm{P}}}={\displaystyle \frac{1}{\mathrm{Cp}}}\mathrm{D}{\displaystyle \frac{\partial \mathrm{g}}{\partial \mathrm{Q}}}={\displaystyle \frac{1+\mathrm{e}0}{\mathsf{\lambda } - \mathsf{\kappa }}}\left(\mathrm{M} - \mathsf{\eta }\right){\displaystyle \frac{\partial \mathrm{g}}{\partial \mathrm{Q}}}$$

(14)

Upon closer examination Equation 14, using stress point (P, Q) is highly restricted to OCR=1, because P and Q are located on the bounding surface. To allow Kp 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 Kp 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 Kp is that Kp must be a positive number at the critical state. However, as given by Equation 14, Kp 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) [38] proposed the concept of a virtual peak stress ratio, Mk, to replace M in Equation 14, thereby providing a smoother and more accurate description of OC behavior. The expression for Mk is shown in Equation 15.

$$\mathrm{Mk}=\mathrm{M}{\mathrm{R}}^{-\mathrm{n}}$$

(15)

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, Mk​ 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 15 is introduced to adjust the virtual stress ratio Mk​, and consequently the plastic modulus Kp​. 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 Mk​ and Kp​. 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 16.

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

(16)

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 (Δ) decreases with increasing OCR. Figure 7 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 7 also reveals that the critical state stress ratio M exerts an additional influence on Δ; specifically, larger values of M correspond to greater values of Δ. Hence, the relationship can be expressed as Δ = f(OCR, M).

In this study, the relationship among Δ, OCR, and M is presented in Figure 4 and expressed in Equation 17. From Equation 17, it can be observed that Δ is more precisely defined as Δ = f(OCR, M²). Accordingly, it may be inferred that the parameter M² 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 18.

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

(17)

$$\mathrm{n}=0.5+M^2+{\displaystyle \frac{2}{\mathrm{OCR}}} (for Drained Loading)$$

(18)

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

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

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 π-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 9. 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 ${\tilde{\mathsf{\sigma }}}_{ij}$

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

$${\tilde{\mathsf{\sigma }}}_{ij}= \tilde{p}{\delta }_{ij}+{\tilde{s}}_{ij}=p{\delta }_{ij}+{\displaystyle \frac{q_{SMP}}{q_c}}{s}_{ij}$$

(19)

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

$$\tilde{p}={\displaystyle \frac{1}{3}}{\tilde{\mathsf{\sigma }}}_{ii}$$

(20)

$${\tilde{s}}_{ij}={\tilde{\mathsf{\sigma }}}_{ij}-\tilde{p}{\delta }_{ij}$$

(21)

$q_{SMP}$ and $q_c$ are the radius of the circular cone and the Extended Mises cone which can be calculated using

$$q_{SMP}=\sqrt{{\displaystyle \frac{2}{3}}}{\displaystyle \frac{2I_1}{3\sqrt{\left(I_1{I}_2-I_3\right)\left(I_1{I}_2-9I_3\right)-1}}}$$

(22)

$$q_c=\sqrt{s_{kl}{s}_{kl}}$$

(23)

where $I_1$, $I_2$ and $I_3$ are the first, second, and third stress invariants, respectively.

3.2. Stress–strain Modelling for Elastoplastic Constitutive Model

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

$$d\epsilon = d{\epsilon }^e+d{\epsilon }^p=\left(d{\epsilon }_v^e+d{\epsilon }_d^e\right)+\left(d{\epsilon }_v^p+d{\epsilon }_d^p\right)$$

(24)

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

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

(25)

Thus,

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

(26)

It is worth noting that $\mathrm{d}\Lambda$ in Equation 24 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 $d{\epsilon }_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}{\mathsf{\epsilon }}_v^e\\ \mathrm{d}{\mathsf{\epsilon }}_d^e\end{array}\right\}$$

(27)

Where

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

(28)

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

$$K={\displaystyle \frac{1+e_0}{\kappa }}p$$

(29)

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

(30)

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

$$\left\{\begin{array}{c}d{\epsilon }_v^e\\ d{\epsilon }_d^e\end{array}\right\}=\left\{\begin{array}{c}d{\epsilon }_v\\ d{\epsilon }_d\end{array}\right\}-\left\{\begin{array}{c}d{\epsilon }_v^p\\ d{\epsilon }_d^p\end{array}\right\}$$

(31)

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

(32)

thus,

$$\mathrm{d}\Lambda ={\displaystyle \frac{\mathrm{K}\frac{\partial \mathrm{F}}{\partial \mathrm{P}}\mathrm{d}{\mathsf{\epsilon }}_{\mathrm{v}}+3\mathrm{G}\frac{\partial \mathrm{F}}{\partial \mathrm{Q}}\mathrm{d}{\mathsf{\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}}}}$$

(33)

$${\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}(\mathsf{\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]$$

(34)

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

$$\left\{\begin{array}{c}dp\\ dq\end{array}\right\}=D_{2\times 2}^{ep}\left\{\begin{array}{c}d{\epsilon }_v\\ d{\epsilon }_d\end{array}\right\}$$

(35)

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 Ψ. This parameter constitutes a key factor enabling the model to reproduce clay behavior as either HPLS or LPLS across different OCR values. Figure 10 depicts the normalized stress paths (p′/p′c – q/p′c) for OCR values of 1, 2, 6, and 20 conjunctions with Ψ 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 Ψ values. This outcome is consistent with the aforementioned aim of this study.

on normalized stress path p’/p’c - q/p’c for various OCR.

4. Model Verifications

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.

4.1. Undrained Analysis

Figure 11, Figure 12, Figure 13 and Figure 14 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 15 and Figure 16.

It is clearly evident from the triaxial compression test results of London clay and Lower Cromer Till, shown in Figure 11, 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 13 and Figure 14, 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.

For the undrained triaxial extension tests on Lower Cromer Till shown in Figure 15, 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 16), 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 17, Figure 18 and Figure 19 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 17, Figure 18 and Figure 19 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.

5. Conclusions

This study advances a unified bounding-surface plasticity formulation that reproduces both high-plasticity and low-plasticity clay responses under low stress ratios across a wide range of OCRs. The model’s core ingredients—a teardrop yield surface governed by two shape parameters Ψ and $\mathsf{\Omega }$, an explicit plastic potential function g that removes the long-standing dilatancy paradox, a linear stress-dilatancy relation D = M−η for practical calibration, a virtual peak ratio $Mk=M{R}^{-n}$linked to OC behavior, and SMP-based strength correction via transformed stresses—work in tandem to deliver a coherent, implementation-friendly framework. Together, these elements provide a smooth NC↔OC transition, preserve numerical robustness, and avoid ad-hoc switching between distinct models when layered deposits mix HPLS and LPLS characteristics.

Verification against undrained triaxial compression data shows that the proposed model captures the atypically steep, near-linear NC stress paths of London Clay and Lower Cromer till (HPLS) where two recent alternatives misclassify the response as LPLS—while maintaining high fidelity for Boston Blue and Kaolin clays (LPLS). Drained comparisons further demonstrate correct reproduction of contractive responses at low OCR and the contractive-to-dilative transition at higher OCR (with appropriate residual strength trends) in Lower Cromer Till, Shanghai clay, and Fujinomori clay. These results collectively highlight that a single Ψ-controlled family of surfaces, coupled with an explicit plastic potential function and OCR-aware hardening, is sufficient to recover the salient features of both LPLS and HPLS clays without inflating parameter sets or resorting to higher-order dilatancy forms.

From an engineering standpoint, the proposed model yields substantial, concrete benefits: (i) a single, unified framework that reproduces both HPLS and LPLS behaviors without ad-hoc model switching; (ii) a smooth NC–OC transition governed by the shape parameters Ψ and $\mathsf{\Omega }$, preserving behavioral continuity across OCRs; (iii) parameter parsimony and straightforward calibration from standard triaxial and oedometer data, reducing practitioner burden; (iv) faithful recovery of steep, near-linear NC stress paths (HPLS) as well as curved LPLS trajectories over OCR ≈ 1–20; (v) consistent capture of contractive responses and contractive-to-dilative transitions together with realistic residual-strength trends; and (vi) interpretability—key parameters tied to measurable indices and critical-state descriptors—facilitating quality control and communication.