Constitutive Modeling of Normally and Over-Consolidated Clay with a High-Order Yield Function

: In this paper, a simple sub-loading yield surface model for both normally consolidated and over-consolidated clay is proposed with emphasis on the effect of the yield surface shape. Compared with the modiﬁed Cam-clay model, only one additional material parameter is introduced to reﬂect geometry features of the yield surface. A higher-order stress–dilatancy relation is given in the current study, leading to a new yield function capable of offering an adequate description of the yield surface of soil samples in the p – q plane. By introducing the concept of the sub-loading yield surface and the uniﬁed hardening parameter, the proposed model can capture the main features of the overconsolidated clay with dilatancy and strain-softening behavior and the main features of the normally consolidated clay with contraction and strain-hardening behavior. The results show that adjusting the yield surface leads to more accurate predictions than the modiﬁed Cam-clay model. The proposed model can also reasonably describe its mechanical behavior for clay samples.


Introduction
The family of Cam-clay models developed on the basis of critical state soil mechanics [1] can describe a broad range of constitutive behaviors of saturated normally consolidated clay with only five parameters. In general, more fitting parameters are needed to capture the mechanical behavior of the over-consolidated clay. Figure 1 shows the typical stress-strain, deformation, and strength properties of clay samples with different over-consolidation ratios (OCR) under the drained compression condition, wherein the solid lines denote the behavior of over-consolidated clay and the dotted lines denote the behavior of normally consolidated clay. For an over-consolidated clay sample, the current stress point reaches the peak strength state before entering the critical state, resulting in the strain-softening behavior in the stress-strain curve. Meanwhile, the clay sample is compacted at first and then expanded to show the dilatancy in the volumetric deformation.
The stress-dilatancy relation plays a crucial role in developing critical state models since dilatancy is an essential factor in characterizing both strength and deformation of soils. One of the distinct features of over-consolidated clay is that the dilatancy closely depends on the degree of over-consolidation [2][3][4]. The modified Cam-clay (MCC) model will underestimate the volume expansion of heavily over-consolidated clays under the drained compression condition because its stress-dilatancy relation only depends on the stress ratio [5]. To overcome this limitation, Yin and Chang [6] developed an empirical stress-dilatancy relation in terms of e/e cs to account for the effect of over-consolidation ratio (OCR) on the dilatancy of clay, where e and e cs are the current and the critical state void ratios at the same stress level, respectively. Similarly, the state parameter proposed by The stress-dilatancy relation plays a crucial role in developing critical state models since dilatancy is an essential factor in characterizing both strength and deformation of soils. One of the distinct features of over-consolidated clay is that the dilatancy closely depends on the degree of over-consolidation [2][3][4]. The modified Cam-clay (MCC) model will underestimate the volume expansion of heavily over-consolidated clays under the drained compression condition because its stress-dilatancy relation only depends on the stress ratio [5]. To overcome this limitation, Yin and Chang [6] developed an empirical stress-dilatancy relation in terms of / to account for the effect of over-consolidation ratio (OCR) on the dilatancy of clay, where and are the current and the critical state void ratios at the same stress level, respectively. Similarly, the state parameter proposed by Been and Jefferies [7], which can be denoted as = − , can also be used to capture the dilative response of an over-consolidated clay [8][9][10]. The formulas of those stressdilatancy relations are quite different from that of the MCC model, which may have some difficulties in developing a practical model for the over-consolidated clay. To simplify the procedure, Gao et al. [11] extended the stress-dilatancy relation of the MCC model by introducing a new parameter to quantify the degree of over-consolidation, which was defined as the ratio between the distances of the current and image stress points from the mapping center. The new relation can be used in either the bounding surface or the subloading surface models for over-consolidated clay.
Another limitation of the MCC model is its elliptical yield surface. The "dry" side of the yield surface above the critical state line has a large elastic zone that will overestimate the undrained shear strength of the over-consolidated clay [12]. In addition, for some geomaterials, experimental observations indicate that the shapes of the yield surface may range from the bullet to the drop [13][14][15]. The shape of the yield surface can also affect the undrained shear strength of the soil. To achieve better performance, the elliptical yield surface of the MCC model has been modified in the literature. For example, Collins [16] developed the -yield function using the thermomechanical approach, where and are two parameters that control the shape of the surface. The -yield function was adopted by Gao et al. [11] to develop a bounding surface model for over-consolidated soil. Another limitation of the MCC model is its elliptical yield surface. The "dry" side of the yield surface above the critical state line has a large elastic zone that will overestimate the undrained shear strength of the over-consolidated clay [12]. In addition, for some geomaterials, experimental observations indicate that the shapes of the yield surface may range from the bullet to the drop [13][14][15]. The shape of the yield surface can also affect the undrained shear strength of the soil. To achieve better performance, the elliptical yield surface of the MCC model has been modified in the literature. For example, Collins [16] developed the α-γ yield function using the thermomechanical approach, where α and γ are two parameters that control the shape of the surface. The α-γ yield function was adopted by Gao et al. [11] to develop a bounding surface model for over-consolidated soil. Panteghini and Lagioia [15] generalized the yield function of the MCC model for geomaterials to describe the arbitrary meridional and deviatoric shapes of yield surface while retaining the features of entire convexity double homothety.
To further improve the performance of the MCC model in modeling the behaviors of over-consolidated clay, the sub-loading surface concept [17,18] was developed, which can capture the irreversible plastic deformation when the current stress point is inside the yield surface [19]. Thus, a smooth transition between the elastic and plastic responses can be obtained. In sub-loading surface models, the change in the size of the sub-loading surface is usually controlled by a similarity ratio R, which is defined as the ratio between the current effective mean stress and the reference effective mean stress [20][21][22]. The concepts of the unified hardening parameter [5,23,24] and the void ratio difference [25][26][27][28] are the alternative hardening parameters that can be used to describe the movement of the sub-loading surface in the principal stress space.
In this study, a high-order stress-dilatancy relation is proposed to describe the dilatancy of normally and over-consolidated clays, which performs better than the MCC model. On this basis, a demonstrative sub-loading surface model is then developed. The yield function representing the sub-loading surface can produce different shapes, ranging from the bullet to the drop. The proposed model introduces the void ratio difference concept to quantify the degree of over-consolidation, which serves as a hardening parameter. The proposed model provides a simple scenario and can predict many essential features of over-consolidated clay, including the strain-softening and dilatancy features.

High-Order Yield Function
Constitutive equations were developed on the basis of critical state soil mechanics and defined in terms of the mean pressure p = tr(σ)/3 and the equivalent von Mises stress q = 3J s,2 where J s,2 = s:s/2 is the second invariant of the deviatoric stress tensor s = σ − pI and σ and I are the stress and the second-order identity tensors, respectively.
In those models, the stress-dilatancy relation defines the relationship between the stress ratio η and the dilatancy ratio d, which can determine the plastic flow direction. In the MCC model, the stress-dilatancy relation is given by where η is the ratio between the mean effective stress p and the deviatoric stress q (η = q/p), d is the ratio between the plastic volumetric strain increment d where p 0 is the intersection of the yield surface with the p-axis in the initial condition. The predicted yield surface is an ellipse in the p-q plane. However, the real shapes of the yield surface may vary from a bullet to a drop for different types of soils. In this study, a simple stress-dilatancy relation with only one additional parameter compared with that of the MCC model is proposed, which is defined as where β is a positive material constant that can be calibrated using the measured dilatancy of soils. The proposed high-order stress-dilatancy relation leads to a generalized MCC-type yield function that is expressed as Figure 2a shows the measured stress-dilatancy relation and yield surface of Bourke silt by Uchaipichat and Khalili [29]. It can be found that the MCC-type stress-dilatancy relation predicts a larger dilatancy ratio d than the measured value at a low stress level. However, the high-order stress-dilatancy relation in Equation (3) can fit the test data well when β = 2.6. In addition, the corresponding yield function in Equation (4) also has a better performance in fitting the measured yield points of Bourke silt than the MCC-type yield function. The comparison between those two yield functions in the three-dimensional stress space is shown in Figure 2b. relation predicts a larger dilatancy ratio than the measured value at a low stress level. However, the high-order stress-dilatancy relation in Equation (3) can fit the test data well when = 2.6. In addition, the corresponding yield function in Equation (4) also has a better performance in fitting the measured yield points of Bourke silt than the MCC-type yield function. The comparison between those two yield functions in the three-dimensional stress space is shown in Figure 2b.

Yield Surfaces
Unlike normally consolidated clay, the initial stress state of over-consolidated clay before shearing is located inside the normal yield surface due to the loading history. It may pass through the critical state line and move backward before the current stress state reaches the yield surface, following a specific stress path (i.e., = 3 for a drained triaxial test). Thus, the stress-strain relation of over-consolidated clay exhibits softening behavior during the post-peak strength loading process, accompanied by dilative volumetric deformation. Note that plastic deformation always occurs in the soil even an infinitesimal loading is applied. This suggests that the current stress point remains on a subloading surface, and the soil sample is yielded during the loading process. As illustrated in Figure 3, in the -plane, the current stress state A ( , ) (i.e., the red circles) and the reference stress state B ( , ) (i.e., the grey circles) are located on the sub-loading surface and the normal yield surface, respectively. We can use the similarity ratio to determine the current size of the sub-loading surface to that of the yield surface. is the inverse of the over-consolidation ratio and is defined as = / ̅ = / . For normally consolidated clay, the sub-loading surface coincides with the normal yield surface as = 1. On the other hand, is less than 1.0 when the soil is over-consolidated, and the sub-loading surface is located inside the normal yield surface.

Yield Surfaces
Unlike normally consolidated clay, the initial stress state of over-consolidated clay before shearing is located inside the normal yield surface due to the loading history. It may pass through the critical state line and move backward before the current stress state reaches the yield surface, following a specific stress path (i.e., dq = 3 dp for a drained triaxial test). Thus, the stress-strain relation of over-consolidated clay exhibits softening behavior during the post-peak strength loading process, accompanied by dilative volumetric deformation. Note that plastic deformation always occurs in the soil even an infinitesimal loading is applied. This suggests that the current stress point remains on a sub-loading surface, and the soil sample is yielded during the loading process. As illustrated in Figure 3, in the p-q plane, the current stress state A (p,q) (i.e., the red circles) and the reference stress state B (p,q) (i.e., the grey circles) are located on the sub-loading surface and the normal yield surface, respectively. We can use the similarity ratio R to determine the current size of the subloading surface to that of the yield surface. R is the inverse of the over-consolidation ratio OCR and is defined as R = p/p = q/q. For normally consolidated clay, the sub-loading surface coincides with the normal yield surface as R = 1. On the other hand, R is less than 1.0 when the soil is over-consolidated, and the sub-loading surface is located inside the normal yield surface.
According to the high order yield function, the reference surface that passes through the reference stress point B(p, q) is expressed as where p x0 is a constant that corresponds to the length of the major axis of the reference yield surface at the initial state of shearing, e 0 is the initial void ratio, and λ and κ are the slopes of the isotropic compression line and the compression/swelling line in the e-p plane. Substitution of the definition of the similarity ratio R leads to  According to the high order yield function, the reference surface that passes through the reference stress point ( ̅ , ) is expressed as where ̅ is a constant that corresponds to the length of the major axis of the reference yield surface at the initial state of shearing, is the initial void ratio, and and are the slopes of the isotropic compression line and the compression/swelling line in the e-p plane.
Substitution of the definition of the similarity ratio leads to Equation (6) shows that the similarity ratio is a function of the current mean effective stress , stress ratio , and plastic volumetric strain . In the sub-loading surface models, plays an essential role because it can reflect the relationship between the subloading surface and the normal yield surface. Usually, the evolution law of is assumed to be a function of the norm of the incremental plastic strain tensor p with d = u( )‖d ‖, where the function ( ) is a monotonously decreasing function of satisfying the conditions of (0) = ∞, ( ) > 0 for 0 < < 1, and (1) = 0. However, this study calculates the value of explicitly, which can be used to further determine the hardening parameter [5,23,24]. Figure 4 shows the evolution of the similarity ratio for over-consolidated clay in a drained triaxial test. It can be found that the stress-strain relation of the soil only exhibits the hardening behavior before the current stress point reaches the peak strength state and is then followed by softening behavior. The phase transformation state occurs earlier than the peak strength state, which is a turning point for the sign of the incremental volumetric strain. The soil is continuously compacted before the phase transformation state and then Equation (6) shows that the similarity ratio R is a function of the current mean effective stress p, stress ratio η, and plastic volumetric strain p v . In the sub-loading surface models, R plays an essential role because it can reflect the relationship between the sub-loading surface and the normal yield surface. Usually, the evolution law of R is assumed to be a function of the norm of the incremental plastic strain tensor p with dR = u(R) d p , where the function u(R) is a monotonously decreasing function of R satisfying the conditions of u(0) = ∞, u(R) > 0 for 0 < R < 1, and u(1) = 0. However, this study calculates the value of R explicitly, which can be used to further determine the hardening parameter , where the function u(R) is a monotonously decreasing function of u satisfying the conditions of u(0) = ∞, u(R) > 0 for 0 < R < 1, and u(1) = 0. However, this study calculates the value of R explicitly, which can be used to further determine the hardening parameter H [5,23,24]. Figure 4 shows the evolution of the similarity ratio R for over-consolidated clay in a drained triaxial test. It can be found that the stress-strain relation of the soil only exhibits the hardening behavior before the current stress point reaches the peak strength state and is then followed by softening behavior. The phase transformation state occurs earlier than the peak strength state, which is a turning point for the sign of the incremental volumetric strain. The soil is continuously compacted before the phase transformation state and then gradually expanded. The soil behaves like a fluid without volumetric deformation at the critical state. Figure 4 also indicates that the initial value of R, in this case, is smaller than 1.0 due to over-consolidation. Continuous loading leads to a monotone increase in R until it reaches 1.0.
Similar to the reference surface, the sub-loading surface that passes through the current stress point A(p, q) is expressed as where p x0 is the length of the major axis of the subloading surface at the initial state before shearing. The hardening parameter H for overconsolidated clays is defined as a function of the current stress ratio η, plastic volumetric strain p v , and potential failure stress ratio M f as follows: Mathematics 2022, 10, x FOR PEER REVIEW 6 of 13 gradually expanded. The soil behaves like a fluid without volumetric deformation at the critical state. Figure 4 also indicates that the initial value of , in this case, is smaller than 1.0 due to over-consolidation. Continuous loading leads to a monotone increase in until it reaches 1.0. Similar to the reference surface, the sub-loading surface that passes through the current stress point ( , ) is expressed as where is the length of the major axis of the subloading surface at the initial state before shearing. The hardening parameter for overconsolidated clays is defined as a function of the current stress ratio , plastic volumetric strain , and potential failure stress ratio as follows: According to Yao et al. [23], the determination of is elaborated in conjunction with the Hvorslev envelope in the subsequent subsection, and is a function of that can be expressed as follows: where parameter is defined as = ( ) .
Because the associated flow rule is adopted, the high-order yield function can also serve as the plastic potential. Thus, the plastic strain increment can be easily calculated as According to Yao et al. [23], the determination of M f is elaborated in conjunction with the Hvorslev envelope in the subsequent subsection, and M f is a function of R that can be expressed as follows: where parameter k is defined as k = M 2 12(3−M) . Because the associated flow rule is adopted, the high-order yield function can also serve as the plastic potential. Thus, the plastic strain increment can be easily calculated as To ensure the current stress point remains on the sub-loading surface during shearing, the consistency condition should be satisfied as follows: with ∂ f ∂p The conventional solution of the incremental stress-strain relationship is based on explicit integration of the following rate constitutive equation: where . σ and .
ε are the incremental stress and strain tensors, respectively. The symbol · represents the Macauley brackets which imply that x = x when plasticity is in effect; otherwise, x = 0.
The elastic stress-strain tensor C e ijkl is given by where δ ij is the Kronecker delta, and K and G are the elastic bulk and shear modulus, respectively. The plastic stress-strain tensor C p ijkl is given by where K p is the plastic modulus with the following definition:

Model Performance
Compared with the original and modified Cam-clay models, the proposed model in this study only introduces one extra model constant as β to adjust the shapes of the yield surface. The other critical state parameters Γ, λ, κ, M, and ν can be calibrated through conventional laboratory tests following a standard procedure. Figure 5 shows the influence of parameter β on the stress-strain relation and volumetric deformation for over-consolidated clay in a drained triaxial test. The relationship between the mean effective stress increment dp and deviatoric stress increment dq is given as dq = 3 dp. It can be found that the tangent modulus is increased with increasing β value, and the soil becomes stiffer, resulting in a higher peak strength. Meanwhile, the soil experiences expansive volumetric deformation after compaction, and the final dilation is also increased with a larger β value.
We further evaluated the influence of the OCR on the constitutive behaviors of clays under the undrained compression condition. In this case, the soil volume remains unchanged, i.e., d υ = 0, and negative pore pressure is developed. Figure 6 gives the model predictions on the stress-strain relation and volumetric deformation for clays with different initial over-consolidated ratios under the drained triaxial compression condition. This model can adequately predict strain hardening and volumetric contraction for normally consolidated or lightly over-consolidated clay; meanwhile, for over-consolidated clay, the strain-softening and volumetric expansion behaviors are also captured by the proposed model. We further evaluated the influence of the on the constitutive behaviors of clays under the undrained compression condition. In this case, the soil volume remains unchanged, i.e., d = 0, and negative pore pressure is developed. Figure 6 gives the model predictions on the stress-strain relation and volumetric deformation for clays with different initial over-consolidated ratios under the drained triaxial compression condition. This model can adequately predict strain hardening and volumetric contraction for normally consolidated or lightly over-consolidated clay; meanwhile, for over-consolidated clay, the strain-softening and volumetric expansion behaviors are also captured by the proposed model.

Model Verification
An essential feature of the proposed model is that it can describe the various shapes of the yield surface with a high-order yield function. By incorporating the sub-loading surface concept, a smooth transition between elastic and plastic deformation is obtained directly since the irreversible plastic strain developed inside the yield surface is captured. To verify the model, experimental data of both normally consolidated and over-consolidated clays in the literature were adopted, i.e., Weald clay (Bishop and Henkel [30]), black  We further evaluated the influence of the on the constitutive behaviors of clays under the undrained compression condition. In this case, the soil volume remains unchanged, i.e., d = 0, and negative pore pressure is developed. Figure 6 gives the model predictions on the stress-strain relation and volumetric deformation for clays with different initial over-consolidated ratios under the drained triaxial compression condition. This model can adequately predict strain hardening and volumetric contraction for normally consolidated or lightly over-consolidated clay; meanwhile, for over-consolidated clay, the strain-softening and volumetric expansion behaviors are also captured by the proposed model.
(a) (b) Figure 6. Effect of OCR on the behaviors of the over-consolidated clay during undrained shearing: (a) the stress ratio-axial strain relation; (b) the pore water pressure-axial strain relation.

Model Verification
An essential feature of the proposed model is that it can describe the various shapes of the yield surface with a high-order yield function. By incorporating the sub-loading surface concept, a smooth transition between elastic and plastic deformation is obtained directly since the irreversible plastic strain developed inside the yield surface is captured. To verify the model, experimental data of both normally consolidated and over-consolidated clays in the literature were adopted, i.e., Weald clay (Bishop and Henkel [30]), black Figure 6. Effect of OCR on the behaviors of the over-consolidated clay during undrained shearing: (a) the stress ratio-axial strain relation; (b) the pore water pressure-axial strain relation.

Model Verification
An essential feature of the proposed model is that it can describe the various shapes of the yield surface with a high-order yield function. By incorporating the sub-loading surface concept, a smooth transition between elastic and plastic deformation is obtained directly since the irreversible plastic strain developed inside the yield surface is captured. To verify the model, experimental data of both normally consolidated and over-consolidated clays in the literature were adopted, i.e., Weald clay (Bishop and Henkel [30]), black Kaolin clay (Zervoyanis [31]), white clay (Biarez and Hicher [32]), and Fujinomori clay (Nakai and Hinokio [25]). All parameters were calibrated and are listed in Table 1.

Weald Clay
Bishop and Henkel [30] conducted a series of undrained triaxial tests on Weald clay, which is a typical normally consolidated clay. During the test, the soil sample was first isotropically consolidated up to 207 kPa of pre-consolidated pressure and then sheared under the undrained condition. The model parameters are given in Table 1. Figure 7 shows the comparison between the measured and predicted results calculated by the proposed model, Cam-clay (CC) model, and modified Cam-clay (MCC) model. Figure 7 indicates that both the original and the modified Cam-clay models predicted a smaller undrained peak strength (i.e., approximately 80 kPa and 110 kPa) than the measured value of 122 kPa. A discrepancy between the measured and calculated pore pressure can also be observed, as shown in Figure 7. However, the proposed model predictions were in good agreement with the measured undrained peak strength and pore pressure. In addition, the calculated effective stress path was quite close to the measured one. (Nakai and Hinokio [25]). All parameters were calibrated and are listed in Table 1.

Weald Clay
Bishop and Henkel [30] conducted a series of undrained triaxial tests on Weald clay, which is a typical normally consolidated clay. During the test, the soil sample was first isotropically consolidated up to 207 kPa of pre-consolidated pressure and then sheared under the undrained condition. The model parameters are given in Table 1. Figure 7 shows the comparison between the measured and predicted results calculated by the proposed model, Cam-clay (CC) model, and modified Cam-clay (MCC) model. Figure 7 indicates that both the original and the modified Cam-clay models predicted a smaller undrained peak strength (i.e., approximately 80 kPa and 110 kPa) than the measured value of 122 kPa. A discrepancy between the measured and calculated pore pressure can also be observed, as shown in Figure 7. However, the proposed model predictions were in good agreement with the measured undrained peak strength and pore pressure. In addition, the calculated effective stress path was quite close to the measured one.  Figure 8 shows the comparisons between the measured and predicted results of the over-consolidated black Kaolin clay subjected to drained triaxial compression. In the tests, the remolded black Kaolin clay samples were first isotropically consolidated up to 800 kPa and then unloaded to 400, 200, and 100 kPa to obtain specific over-consolidation ratios OCR = 1, 2, 4, and 8,. After that, axial loading was applied to the soil samples until they failed, so as to investigate the influence of the over-consolidation ratio. All the model parameters are listed in Table 1. As indicated in Figure 8, the model predictions fit the measured data reasonably well for different OCRs. Furthermore, the stress-strain relations for OCR = 1 and 2 showed a continuous hardening characteristic due to the contractive void ratio change in Figure 8. By contrast, the stress-strain relations for OCR = 4 and 8 showed a slightly softening feature accompanied by a dilative volumetric deformation.  Figure 8 shows the comparisons between the measured and predicted results of the over-consolidated black Kaolin clay subjected to drained triaxial compression. In the tests, the remolded black Kaolin clay samples were first isotropically consolidated up to 800 kPa and then unloaded to 400, 200, and 100 kPa to obtain specific over-consolidation ratios = 1, 2, 4, and 8,. After that, axial loading was applied to the soil samples until they failed, so as to investigate the influence of the over-consolidation ratio. All the model parameters are listed in Table 1. As indicated in Figure 8, the model predictions fit the measured data reasonably well for different . Furthermore, the stress-strain relations for = 1 and 2 showed a continuous hardening characteristic due to the contractive void ratio change in Figure 8. By contrast, the stress-strain relations for = 4 and 8 showed a slightly softening feature accompanied by a dilative volumetric deformation.  Figure 9 shows the comparison between the measured and predicted results of the over-consolidated white clay subjected to undrained triaxial compression. During the tests, three samples were first isotropically consolidated up to 800 kPa, and two of them were then unloaded to 400 and 67 kPa to obtain the values of = 4 and 12, respectively. All the model parameters can be found in Table 1. As calibrated by Yin et al. [6], the reference critical state void 0 is 1.21, and the compression index and the swelling index are 0.089 and 0.034, respectively. Therefore, the initial void ratio 0 at different overconsolidated ratios ( = 1, 2, and 12) was 0.653, 0.677, and 0.737, respectively.  Figure 9 shows the comparison between the measured and predicted results of the over-consolidated white clay subjected to undrained triaxial compression. During the tests, three samples were first isotropically consolidated up to 800 kPa, and two of them were then unloaded to 400 and 67 kPa to obtain the values of OCR = 4 and 12, respectively. All the model parameters can be found in Table 1. As calibrated by Yin et al. [6], the reference critical state void e cr0 is 1.21, and the compression index λ and the swelling index κ are 0.089 and 0.034, respectively. Therefore, the initial void ratio e 0 at different over-consolidated ratios (OCR = 1, 2, and 12) was 0.653, 0.677, and 0.737, respectively.

White Clay
As indicated in Figure 9, the model predictions were in good agreement with the measured values. The measured effective stress paths in Figure 9b show that the stress paths did not go beyond the critical state line for the normally consolidated and lightly over-consolidated samples. In those cases, the developed pore pressure would be positive. On the other hand, for a heavily over-consolidated sample (i.e., OCR = 12), the current stress point would pass through the critical state line and then move backward, resulting in negative pore pressure. Overall, the proposed model could predict the constitutive behaviors of undrained tests for specimens with different over-consolidated ratios. As indicated in Figure 9, the model predictions were in good agreement with the measured values. The measured effective stress paths in Figure 9b show that the stress paths did not go beyond the critical state line for the normally consolidated and lightly over-consolidated samples. In those cases, the developed pore pressure would be positive. On the other hand, for a heavily over-consolidated sample (i.e., = 12), the current stress point would pass through the critical state line and then move backward, resulting in negative pore pressure. Overall, the proposed model could predict the constitutive behaviors of undrained tests for specimens with different over-consolidated ratios. Figure 10 shows the comparison between the measured and predicted results of the over-consolidated Fujinomori clay subjected to drained triaxial compression. In the tests, the remolded samples were prepared with = 1, 2, 4, and 8, respectively. All the model parameters are listed in Table 1. It can be seen from Figure 10 that the predicted peak strength was close to the measured data, while the corresponding axial strain was overestimated by the proposed model. At the same time, a discrepancy in the volumetric deformation can also be observed. To some degree, the proposed model could capture the mechanical behavior of Fujinomori clay, especially when the sample was normally consolidated with OCR = 1.  Figure 10 shows the comparison between the measured and predicted results of the over-consolidated Fujinomori clay subjected to drained triaxial compression. In the tests, the remolded samples were prepared with OCR = 1, 2, 4, and 8, respectively. All the model parameters are listed in Table 1. It can be seen from Figure 10 that the predicted peak strength was close to the measured data, while the corresponding axial strain was overestimated by the proposed model. At the same time, a discrepancy in the volumetric deformation can also be observed. To some degree, the proposed model could capture the mechanical behavior of Fujinomori clay, especially when the sample was normally consolidated with OCR = 1. As indicated in Figure 9, the model predictions were in good agreement with the measured values. The measured effective stress paths in Figure 9b show that the stress paths did not go beyond the critical state line for the normally consolidated and lightly over-consolidated samples. In those cases, the developed pore pressure would be positive. On the other hand, for a heavily over-consolidated sample (i.e., = 12), the current stress point would pass through the critical state line and then move backward, resulting in negative pore pressure. Overall, the proposed model could predict the constitutive behaviors of undrained tests for specimens with different over-consolidated ratios. Figure 10 shows the comparison between the measured and predicted results of the over-consolidated Fujinomori clay subjected to drained triaxial compression. In the tests, the remolded samples were prepared with = 1, 2, 4, and 8, respectively. All the model parameters are listed in Table 1. It can be seen from Figure 10 that the predicted peak strength was close to the measured data, while the corresponding axial strain was overestimated by the proposed model. At the same time, a discrepancy in the volumetric deformation can also be observed. To some degree, the proposed model could capture the mechanical behavior of Fujinomori clay, especially when the sample was normally consolidated with OCR = 1. As shown in Figures 7-10, the proposed model could provide a unified approach for describing the constitutive behaviors of both normally consolidated and over-consolidated clays.

Conclusions
This paper developed a new unified sub-loading surface model for both normally and over-consolidated clay, with particular attention paid to the effect of yield surface shape.
Firstly, a new and straightforward dilatancy equation with high order was proposed by introducing only one parameter compared with the MCC model. The yield surface function was, therefore, obtained by integrating the dilatancy equation. Such a yield function could describe the yield surface with varying shapes.
Secondly, a sub-loading constitutive model was developed on the basis of the newly proposed yield surface function and the concept of a unified hardening parameter. It was demonstrated that the proposed model can perform well in describing the mechanical behaviors of both normally and over-consolidated clay.