An Elastic-Plastic Constitutive Model for Unsaturated Structural Loess

Abstract

The water sensitivity and structural characteristics of collapsible loess are two typical factors that significantly influence its mechanical behaviors. This paper presents a simple and practical elastic-plastic model based on the modified Cam-Clay model to well capture the essential behavior of unsaturated intact loess. The model employs deviator stress and spheric stress as the stress variables, with the water content serving as the moisture variable. The critical state surface of the model can be determined by utilizing the shear strength parameters of unsaturated soil under axisymmetric stress conditions. An initial yield surface equation is established by incorporating structural strength into the elliptical yield surface equation, which is used to determine the starting point for elastic-plastic deformation calculations under different humidity and stress combinations. The model comprises several parameters, each of which has a clear physical interpretation and can be conveniently obtained through conventional triaxial tests. The validity of the model for unsaturated intact loess is confirmed through a comparison with the stress–strain relationship of unsaturated intact loess in the axisymmetric stress state. This work has the potential to significantly enhance our ability to predict and mitigate potential geotechnical disasters, such as foundation deformation under axisymmetric conditions and slope stability problems under non-axisymmetric conditions. Ultimately, the application of this model could contribute to the safety and stability of infrastructure and construction projects in loess regions.

1. Introduction

Loess is extensively distributed worldwide, with a continuous coverage area of approximately 63 × 10⁴ km² in northwest China [1]. Due to the arid climate characteristics and regional geological processes, the intact loess has developed significant structural and water-sensitive characteristics [2,3]. The intact loess undergoes strength softening and collapse deformation once exposed to water, even with the acting pressure unchanged, which leads to a series of geological disasters [4,5,6]. The structure and unsaturation of intact loess are crucial factors for accurately understanding and evaluating its mechanical behavior, drawing significant attention from scientists [7,8,9,10,11]. Therefore, it is critical to illustrate the main feature of unsaturated structural loess using a simple and practical model.

In the case of unsaturated soil, it is essential to consider not only the stress variables but also the humidity variables when assessing strength and deformation behaviors. Many constitutive models utilize deviatoric stress and net stress as stress variables, such as the modified Cam-Clay (MCC) model. This model effectively links volume change and shear strength based on elastic-plastic theory, enabling the description of the critical state and stress–strain relationship of saturated soil, and it serves as the foundation for numerous elastic-plastic constitutive models [12]. Matric suction is commonly employed as a humidity variable in constitutive models of unsaturated soil. Alonso et al. [13] described the yield characteristics of unsaturated soil under the combined action of suction and net stress by establishing a loading-collapse (LC) curve, leading to the proposal of the well-known Barcelona basic model (BBM) to predict the mechanical behaviors of unsaturated soil. However, hydraulic hysteresis can lead to differences in the mechanical properties of soil with the same suction due to varying degrees of saturation [14]. The variation in saturation is dependent on the suction change and the soil deformation, serving as a state variable in some elastoplastic models to capture the hydro-mechanical behaviors of unsaturated soil [15,16,17]. Nonetheless, saturation changes with soil deformation caused by loading make it unsuitable for describing mechanical properties under the coupling of stress and moisture. To address the irreversible behavior of unsaturated soils under drying/wetting and various loading conditions, Sheng et al. [18] improved the Barcelona basic model by adopting the non-associated flow rules and incorporating the two suction-related yield surfaces to model unsaturated soil mechanical behaviors. While these models enhance the accuracy of simulating the hydro-mechanical behaviors of unsaturated soil, their complexity limits their application and promotion in engineering.

Intact soil typically possesses a natural structure that profoundly influences its mechanical behavior. The isotropic stress compression curve of intact soil exhibits obvious turning points, with the corresponding stress known as the structural compression yield stress or initial yield stress. This stress is widely recognized as the boundary parameter reflecting elastic compression deformation and elastic-plastic compression deformation [10]. Zhang et al. [19] incorporated the initial yield stress into the yield function of the modified Cam-Clay model to capture the structural properties of natural loess. Zhu and Yao [20] proposed a structured unified hardening model to account for soil structure, further extending the model by introducing an additional yielding mechanism to consider the effect of overconsolidation on the volume change in the microstructure during wetting [21]. In order to capture the effects of soil structure on the loading-collapse (LC) curve, Mu et al. [22] introduced an elastoplastic constitutive model that incorporates a normalized loading-collapse (LC) curve with structure degradation, enabling the simulation of its evolution and the initial soil structure under hydro-mechanical loads. Considering the water retention and evolution of loess structure, Weng et al. [23] proposed a coupling model based on the Barcelona basic model, and it can predict the hydro-mechanical behaviors of loess with different stress paths. Yao et al. [24] introduced evolution equations for loess structure during wetting and loading into the Barcelona basic model and proposed an elastoplastic model. To describe the behavior of unsaturated structured soils under hydro-mechanical loading, a boundary surface plastic model was constructed to achieve a smooth transition between structured and unstructured states by controlling the size of the boundary surface in reference [25,26]. However, the complexity of these models may limit their engineering applications. It is essential for a model to be simple, with few parameters, each with a clear meaning and being easily determinable, while effectively capturing the essential behaviors of soil.

The paper proposes a simple and practical model to describe the primary mechanical behaviors of intact loess. The model highlights the two key characteristics of intact loess: water sensitivity and structure, while minimizing the influence of other factors on its mechanical properties. By investigating the effects of humidity and structural strength on the mechanical response, a constitutive relationship tailored for unsaturated structured loess has been developed. This model utilizes the strength parameters of unsaturated soil to establish the critical state line, and it introduces the structural strength of unsaturated loess into the elliptical yield surface equation to determine the initial yield surface. Based on the modified Cam-Clay model, an elastic-plastic constitutive model for unsaturated structured loess is constructed. The effectiveness of the model is validated by comparing its predictions with triaxial test data.

2. Constitutive Model

2.1. The Critical State Line of Unsaturated Intact Loess under Axisymmetric Stress Conditions

The stress state variables, spheric stress p and deviator stress q, which are the same as those in the modified Cam-Clay model, are utilized in this study. The corresponding expressions are as follows:

$$p={\displaystyle \frac{{\sigma }_1+{\sigma }_2+{\sigma }_3}{3}}={\displaystyle \frac{{\sigma }_1+2{\sigma }_3}{3}}$$

(1)

$$q={\displaystyle \frac{1}{\sqrt{2}}}·\sqrt{{\left({\sigma }_1-{\sigma }_2\right)}^2+{\left({\sigma }_2-{\sigma }_3\right)}^2+{\left({\sigma }_3-{\sigma }_1\right)}^2}={\sigma }_1-{\sigma }_3$$

(2)

here, σ₁ is the major principal stress, and σ₃ is the minor principal stress; in an axisymmetric stress state, it is believed that σ₂ = σ₃.

The critical state indicates that the soil undergoes shear failure under the action of the stress state. In the case of normally consolidated saturated soil, the critical state line is a straight line that passes through the origin of the p-q stress coordinates [12]. The unsaturated intact loess consistently exhibits high cohesion and tensile strength owing to its structure [27]. Consequently, the critical state line of the unsaturated intact loess does not pass through the origin, but rather has an intercept ξ on the q axis, as illustrated in Figure 1. The equation for the critical state line is expressed as follows:

$$q=\xi +Mp$$

(3)

here, M is the slope of the critical state line.

For geotechnical materials, the Mohr–Coulomb failure criterion is utilized. The coordinates of the tangent points τ_f and σ can be expressed in terms of the major and minor principle stresses σ₁ and σ₃ of the critical stress circle, as depicted in Figure 2.

The Mohr–Coulomb shear strength formula is widely employed for the shear testing of all types of soil, and it is expressed as follows:

$${\tau }_f={\sigma }_f\tan \phi +c$$

(4)

Based on the geometric relationship in Figure 2, the following relationship can be derived:

$${\displaystyle \frac{{\sigma }_1-{\sigma }_3}{2}}={\displaystyle \frac{{\sigma }_1+{\sigma }_3}{2}}\sin \phi +c·\cos \phi$$

(5)

In three-dimensional space, the major principal stress σ₁ and minor principal stress σ₃ can be represented by the net stress p and deviator stress q according to Equations (1) and (2), and their expressions are as follows:

$${\sigma }_1=p+{\displaystyle \frac{2}{3}}q$$

(6)

$${\sigma }_3=p-{\displaystyle \frac{1}{3}}q$$

(7)

The relationship equation between the deviatoric stress q_f and net stress p_f in the failure state can be obtained by substituting Equations (6) and (7) into Equation (5) as follows:

$$q_f={\displaystyle \frac{6\sin \phi }{3-\sin \phi }}·p_f+{\displaystyle \frac{6c·\cos \phi }{3-\sin \phi }}$$

(8)

Hence, according to the Mohr–Coulomb failure criterion, under the axisymmetric stress state, the relationship equation between slope M and internal friction angle φ, and the relationship equation among the intercept ξ, internal friction angle φ, and cohesion c can be obtained from Equations (3) and (8) as follows:

$$M={\displaystyle \frac{6\sin \phi }{3-\sin \phi }}$$

(9)

$$\xi ={\displaystyle \frac{6c·\cos \phi }{3-\sin \phi }}$$

(10)

Based on the shear test data of unsaturated intact loess in different regions, it is evident that the shear strength parameters exhibit a consistent pattern of variation with the water content [28,29]. The test data demonstrate that the internal friction angle of the intact loess remains approximately constant across different water contents, as depicted in Figure 3. Meanwhile, the cohesion decreases with an increasing water content, as illustrated in Figure 4. The relationship between the water content and cohesion can be accurately modeled using a power function as follows:

$$c=A{w}^B$$

(11)

here, A and B are fitting parameters.

Figure 3 indicates that w has minimal impact on the internal friction angle (φ) of the intact loess. Consequently, the value of M for the intact loess with varying water contents can be considered constant. The intercept ξ can be determined by substituting Equation (11) into Equation (10), as follows:

$$\xi =A^{\prime }{w}^B$$

(12)

here, A′ = 6Acosφ/(3 − sinφ).

Assuming that the internal friction angle (φ) remains constant, parameter A′ also remains constant. Consequently, the reduction in the intercept ξ follows a non-linear relationship with the water content (w), which differs from the assumption in the Barcelona basic model [13].

Overall, the cohesion (c) and internal friction angle (φ) of the intact loess samples with varying water contents can be measured using conventional triaxial shear tests. By analyzing the variation in the internal friction angle (φ) and cohesion (c) with water content, it is determined that parameters A and B in Equation (12) can be obtained by fitting the trend of the variation in cohesion (c) with the water content using a power function. By substituting these values into Equations (9) and (12), the intercept and slope of the critical state line in Equation (3) can be obtained. Hence, the expression for the critical state surface can be determined.

2.2. The Initial Yield Surface of Unsaturated Intact Loess

The assumption that the yield surface of the unsaturated intact loess is elliptical is in line with the modified Cam-Clay model. The yield surface of the unsaturated intact loess in the p-q space is depicted in Figure 5, and it aligns with the yield surface of cemented clay described in the literature [30]. The equation for the elliptical yield surface in this study is as follows:

$${\left({\displaystyle \frac{q}{\beta }}\right)}^2+{\left({\displaystyle \frac{p-\gamma p_0}{\alpha }}\right)}^2={p_0}^2$$

(13)

here, p₀ is the intercept of the ellipse on the p axis; parameters α and β are removed during the derivation. The geometric relationship is shown in Figure 5.

The equations corresponding to the geometric relationships are as follows:

$$\alpha p_0={\displaystyle \frac{1}{2}}\left(p_0+\xi /M\right)$$

(14)

$$\beta p_0={\displaystyle \frac{1}{2}}\left(p_0+\xi /M\right)M$$

(15)

$$\gamma p_0={\displaystyle \frac{1}{2}}\left(p_0-\xi /M\right)$$

(16)

The yield function equation can be derived by substituting Equations (14)–(16) into Equation (13), as follows:

$$q^2+M^2{p}^2+Mp\xi =p_0\left(M^{2}p+M\xi \right)$$

(17)

here, p₀ is the stress state that produces the same hardening parameters under the isotropic compression stress path.

Intact loess is a typical underconsolidated soil due to its formation by wind accumulation. However, its compression characteristics during isotropic compression testing resemble those of overconsolidated soil, owing to its special cementation structure. The pressure corresponding to the turning point of its e-lnp curve is called the structural yield pressure p_y [9]. It is assumed that only elastic deformation occurs when the stress state of the intact loess is p ≤ p_y, and the elastic-plastic deformation begins when p > p_y.

The cementation structure of the intact loess partly consists of soluble salt, which can be weakened or damaged when exposed to water, leading to a decrease in the structural yield pressure p_y [31,32,33]. Consequently, even at the same sampling depth, samples with different water contents exhibit significant differences in the structural yield pressure p_y, unlike the preconsolidation pressure. Based on the compression test data of unsaturated intact loess in different regions [28,34], it is shown that the structural yield pressure p_y decreases with the water content w, as illustrated in Figure 6. The relationship between the structural yield pressure p_y and water content can be fitted using a power function as follows:

$$p_y=R{w}^T$$

(18)

Hence, the initial yield surface equation can be obtained via p₀ = p_y. The initial yield surface equation can be derived by substituting Equations (12) and (18) into Equation (17), as follows:

$$q^2+M^2{p}^2+Mp{A}^{\prime }{w}^B-R{w}^T\left(M^{2}p+M{A}^{\prime }{w}^B\right)=0$$

(19)

The initial yield surface of the intact loess in the p-q-w space is depicted in Figure 7. It is assumed that, within this range, only elastic deformation occurs, and no elastic-plastic deformation occurs under the combined action of humidity and stress.

2.3. Hardening Laws

The variation curve of the pore ratio (e) with loading pressure (p) for the unsaturated intact loess under isotropic compression or lateral compression conditions exhibits a clear turning point in the e-lnp coordinate system, as illustrated in Figure 8.

When p ≤ p_y, the loading pressure p does not exceed the structural yield pressure p_y, and it seems that the soil only undergoes elastic deformation without plastic deformation. When p = p_y, the porosity ratio decreases from e₀ to e_y. The equation for the porosity ratio e under the loading pressure p during this stage is as follows:

$$e=e_0-\kappa \left(\ln p-\ln p_a\right)$$

(20)

The equation for the corresponding volumetric strain ε_v is as follows:

$${\epsilon }_v={\displaystyle \frac{\kappa }{1+e_0}}\left(\ln p-\ln p_a\right)$$

(21)

here, κ is the elastic compression parameter, which is the slope of the elastic section of the compression curve before yielding, and e₀ is the initial porosity ratio of the intact loess with pressure p_a; to avoid being meaningless when taking logarithmic coordinates, the value of p_a can be taken as 1 kPa.

When p > p_y, the structure of the intact loess is damaged, plastic deformation begins to occur, and the deformation of the soil is elastic-plastic deformation, which includes two parts during this stage: elastic deformation and plastic deformation. The equation for the elastic-plastic deformation of the intact loess that occurs after yielding is as follows:

$$\Delta e=e-e_y=\lambda \left(\ln p_0-\ln p_y\right)$$

(22)

here, λ is the elastic-plastic compression parameter, which is the slope of the elastic-plastic section of the compression curve after yielding.

The equation for the corresponding elastic-plastic volumetric strain ε_v is as follows:

$${\epsilon }_v={\displaystyle \frac{\lambda }{e_0+1}}\left(\ln p_0-\ln p_y\right)$$

(23)

The equation for the elastic deformation in the elastic-plastic deformation is as follows:

$$\Delta e^{\mathrm{e}}=\kappa \left(\ln p_0-\ln p_y\right)$$

(24)

The equation for the corresponding elastic volumetric strain ε_v^e is as follows:

$${{\epsilon }_v}^{\mathrm{e}}={\displaystyle \frac{\kappa }{e_0+1}}\left(\ln p_0-\ln p_y\right)$$

(25)

Therefore, the equation for the corresponding plastic volumetric strain ε_v^p is as follows:

$${{\epsilon }_v}^p={\epsilon }_v-{{\epsilon }_v}^e={\displaystyle \frac{\lambda -\kappa }{e_0+1}}\left(\ln p_0-\ln p_y\right)$$

(26)

The plastic volumetric strain ε_v^p is selected as the hardening parameter in the model. Parameter p₀ is the stress state of the soil that produces the same hardening parameter ε_v^p under the isotropic compression stress path. The yield function equation can be derived by substituting Equation (26) into Equation (17), as follows:

$${\displaystyle \frac{\lambda -\kappa }{e_0+1}}\ln \left({\displaystyle \frac{q^2}{M^2{p}^2}}+1+{\displaystyle \frac{\xi }{Mp}}\right)-{\displaystyle \frac{\lambda -\kappa }{e_0+1}}\ln \left({\displaystyle \frac{p_y}{p}}+{\displaystyle \frac{p_y\xi }{M{p}^2}}\right)-{\epsilon }_v^p=0$$

(27)

A series of triaxial isotropic compression tests were conducted on unsaturated intact loess in the literature [35], and the variations in parameters λ and κ with water content w were studied, as shown in Figure 9.

Parameter κ, however, fluctuated in a small range with the change in the water content, and it could be assumed to be independent of the water content; thus, the average value could be taken. Parameter λ increased with the water content, which could be described by an empirical formula as follows:

$$\lambda ={\alpha }^{\prime }\ln w-{\beta }^{\prime }$$

(28)

The variation rule of parameters λ and κ with humidity variables is consistent with the conclusions in previous studies [17,19,36].

2.4. Flow Rules

The flow rule determines the direction of plastic strain, and the hardening law determines the relative magnitude of plastic strain. The model adopts the associative flow rule, and assuming that the plastic potential function and yield function are the same, the equation is as follows:

$$F=g={\displaystyle \frac{\lambda -\kappa }{e_0+1}}\ln \left({\displaystyle \frac{q^2}{M^2{p}^2}}+1+{\displaystyle \frac{\xi }{Mp}}\right)-{\displaystyle \frac{\lambda -\kappa }{e_0+1}}\ln \left({\displaystyle \frac{p_y}{p}}+{\displaystyle \frac{p_y\xi }{M{p}^2}}\right)-{\epsilon }_v^p=0$$

(29)

The equation for plastic deformation can be written as follows:

$$d{\epsilon }_v^P=\mathsf{\Lambda }{\displaystyle \frac{\partial F}{\partial p}},d{\epsilon }_s^q=\mathsf{\Lambda }{\displaystyle \frac{\partial F}{\partial q}}$$

(30)

here, Λ is the plastic multiplier.

$$\mathsf{\Lambda }=-{\displaystyle \frac{\left(\partial F/\partial p\right)\mathrm{d}p+\left(\partial F/\partial q\right)\mathrm{d}q}{\left(\partial F/\partial {{\epsilon }_v}^p\right)\left(\partial F/\partial p\right)}}$$

(31)

By taking the derivatives of the different variables in Equation (29), the following equations can be obtained:

$${\displaystyle \frac{\partial F}{\partial p}}={\displaystyle \frac{\lambda -\kappa }{1+e_0}}{\displaystyle \frac{1}{p+\frac{\xi }{M}}}{\displaystyle \frac{M^2{p}^2-q^2+2Mp\xi +{\xi }^2}{M^2{p}^2+q^2+Mp\xi }}$$

(32)

$${\displaystyle \frac{\partial F}{\partial q}}={\displaystyle \frac{\lambda -\kappa }{1+e_0}}{\displaystyle \frac{2q}{M^2{p}^2+q^2+Mp\xi }}$$

(33)

$${\displaystyle \frac{\partial F}{\partial {\epsilon }_v^p}}=-1$$

(34)

The plastic multiplier Λ equation can be obtained by substituting Equations (32)–(34) into Equation (31), as follows:

$$\mathsf{\Lambda }=dp+{\displaystyle \frac{2q\left(p+\frac{\xi }{M}\right)}{M^2{p}^2-q^2+2Mp\xi +{\xi }^2}}dq$$

(35)

2.5. Stress–Strain Relationships

The increment equations for plastic volumetric strain (dε_v^p) and plastic shear strain (dε_s^p) can be obtained by substituting Equation (35) into Equation (30), as follows:

$$d{{\epsilon }_v}^p={\displaystyle \frac{\lambda -\kappa }{1+e_0}}{\displaystyle \frac{1}{p+\frac{\xi }{M}}}{\displaystyle \frac{M^2{p}^2-q^2+2Mp\xi +{\xi }^2}{M^2{p}^2+q^2+Mp\xi }}dp+{\displaystyle \frac{\lambda -\kappa }{1+e_0}}{\displaystyle \frac{2q}{M^2{p}^2+q^2+Mp\xi }}dq$$

(36)

$$d{{\epsilon }_s}^p={\displaystyle \frac{\lambda -\kappa }{1+e_0}}{\displaystyle \frac{2q}{M^2{p}^2+q^2+Mp\xi }}dp+{\displaystyle \frac{\lambda -\kappa }{1+e_0}}{\displaystyle \frac{2q}{M^2{p}^2+q^2+Mp\xi }}{\displaystyle \frac{2q\left(p+\frac{\xi }{M}\right)}{M^2{p}^2-q^2+2Mp\xi +{\xi }^2}}dq$$

(37)

The matrix form can be written as follows:

$$\left[\begin{array}{c}d{\epsilon }_v^P\\ d{\epsilon }_s^p\end{array}\right]={\displaystyle \frac{\lambda -\kappa }{1+e_0}}{\displaystyle \frac{2q}{M^2{p}^2+q^2+Mp\xi }}\left[\begin{array}{cc}{\displaystyle \frac{\chi }{2pq+\frac{2q\xi }{M}}}& 1\\ 1& {\displaystyle \frac{2pq+\frac{2q\xi }{M}}{\chi }}\end{array}\right]\left[\begin{array}{c}dp\\ dq\end{array}\right]$$

(38)

here, χ = M²p² − q² + 2Mpξ + ξ²

According to the calculation methods in the literature [37], the elastic strain increment equations are as follows:

$$d{{\epsilon }_v}^e={\displaystyle \frac{\kappa }{1+e_0}}{\displaystyle \frac{dp}{p}}$$

(39)

$$d{{\epsilon }_s}^{\mathrm{e}}={\displaystyle \frac{2\left(1+\mu \right)}{9\left(1-2\mu \right)}}{\displaystyle \frac{\kappa }{1+e_0}}{\displaystyle \frac{dq}{p}}$$

(40)

The total deformation increment includes two parts: the elastic deformation increment and the plastic deformation increment. The increment equations are as follows:

$$\mathrm{d}{\epsilon }_v=\mathrm{d}{{\epsilon }_v}^e+\mathrm{d}{{\epsilon }_v}^p,\mathrm{d}{\epsilon }_s=\mathrm{d}{{\epsilon }_s}^e+\mathrm{d}{{\epsilon }_s}^p$$

(41)

Elastic deformation only occurs within the initial yield surface, and the total deformation of the soil can be obtained through Equations (39) and (40). Elastic-plastic deformation occurs when the combination of stress and humidity exceeds the initial yield surface, and the total deformation of the soil needs to be solved using Equation (41). Thus, the stress–strain relationship of unsaturated intact loess can be obtained.

3. Model Prediction and Validation

The model comprises eight parameters, with three parameters (A′, B, and M) derived from triaxial shear tests with various water contents to characterize the critical state surface, and the remaining five parameters (α′, β′, κ, R, and T) determined through isotropic compression tests with different moisture contents to define the hardening law and initial yield surface. The calculation results will be compared with triaxial shear data of intact loess across various regions to validate the effectiveness of the model.

3.1. Verification of Intact Loess in the Guanzhong Region

The test samples are representative unsaturated Q₃ intact loess, and they were sourced from a pit with a depth of 3.5 m in Xi′an City, Shaanxi Province [35]. The natural moisture content of the test samples is 13.6%, and the natural dry density is 1.35 g/cm³. The initial porosity ratio of the intact loess is 0.989, and the plasticity index is 13.5. The model parameter values in the literature are determined to be as follows: A’ = 82,185, B = −2.1, M = 0.88; α′ = 0.106, β′ = 0.25, κ = 0.005, R = 1931, and T = −0.73. The triaxial shear test results and model predictions are shown in Figure 10.

The comparison between the model predictions and the test data indicates that the model effectively simulates the stress-strain behavior of unsaturated intact loess, particularly the stress–strain relationship of the strain-hardening failure modes. However, the model’s prediction of the stress–strain relationship for the strain-softening failure type is deemed irrational. This failure type primarily occurs in samples with a low moisture content under low confining pressure.

3.2. Verification of Intact Loess in the Northern Shaanxi Region

The samples were collected from a slope in the Yan’an new area, Shaanxi Province. The soil sample is a typical Q₃ loess, excavated from a depth of 5 m below the top of the slope [38]. The natural moisture content of the test samples is 10.6%, and the natural dry density is 1.45 g/cm³. The initial porosity ratio of the intact loess is 1.06, and the plasticity index is 11.82. The model parameter values in the literature are determined to be as follows: A’ = 2004, B = −1.3, M = 0.98; α’ = 0.15, β’ = 0.29, κ = 0.01, R = 5509, and T = −1.46. The triaxial shear test results and model predictions are shown in Figure 11.

The predictions closely match the experimental data presented in Figure 11, indicating that the model can effectively capture the stress–strain relationship of unsaturated intact loess. However, there are slight deviations for weak strain hardening and strong strain hardening. The deviations for weak strain hardening mainly occur within the range of data points with strain less than 4%, while the deviations for strong strain hardening are primarily concentrated in the data points with strain greater than 12%.

3.3. Verification of Intact Loess in Gansu Province

The test samples are Q₃ intact loess, and they were sourced from a pit with a depth of 5 m in Lanzhou City, Gansu Province [39]. The natural moisture content of the test samples is 8%, and the natural dry density is 1.30 g/cm³. The initial porosity ratio of the intact loess is 1.075, and the plasticity index is 8.4. The model parameter values in the literature are determined to be as follows: A’ = 1418, B = −1.38, M = 1.35; α’ = 0.018, β’ = 0.053, κ = 0.04, R = 2295, and T = −1.26. The triaxial shear test results and model predictions are shown in Figure 12.

Figure 12 shows that the predictions for intact loess under varying moisture contents and stress conditions align well with the experimental data. According to the literature, the loess in Lanzhou is classified as silt-like loess, indicating that the model performs well in fitting this type of soil.

In summary, comparisons with experimental data from representative loess across various regions of China indicate that the model exhibits good applicability, especially for silt-like loess. However, there are notable deviations in fitting the stress–strain relationships for different failure types; the model performs well for strain hardening types but struggles with strain softening. These deviations are primarily concentrated in data points characterized by low moisture content and low stress conditions. This may be attributed to the relatively high structural strength of the intact loess at low water content, as the structural strength used in the calculations is derived from isotropic compression tests. The stress state under a triaxial shear stress path at low confining pressure shows significant deviation compared to that under an isotropic compression stress path. Hence, it will be essential to investigate the structural strength of intact loess under various stress paths in future research. Fortunately, intact loess with a low moisture content generally demonstrates favorable mechanical properties, and engineering disasters are rare under low stress conditions.

4. Conclusions

This study establishes an elastoplastic model with few parameters, and each parameter can be easily determined based on the main features of unsaturated intact loess. The model can predict the stress–strain relationship of intact loess with different humidity states. The conclusions drawn from this study are as follows:

(1) The critical state surface of unsaturated intact loess is constructed based on the variation in strength index with water content, which can be determined using three parameters obtained from conventional triaxial shear tests.

(2) A method for determining the initial yield surface is proposed, and it can determine the starting point of elastic-plastic deformation calculations under different humidity and stress combinations.

(3) The model contains eight parameters, and each parameter can be determined using conventional triaxial tests. The effectiveness of the model fitting results is verified by comparing them with existing triaxial test data.