A Nonlinear Constitutive Model for Remoulded Fine-Grained Materials Used under the Qinghai–Tibet Railway Line

Using undrained triaxial shear tests, this study investigates the mechanical properties of fine-grained materials (silty clay and sand) which are extensively used for China’s Qinghai–Tibet Railway (QTR) under different confining pressures (σ3) and freezing temperatures (T). The results show that a reduction in T causes an increase in the shear strength and elastic modulus of all the materials tested in the present study. In addition, the freezing of the silty clay has no significant effect on the type of soil behaviour (strain-hardening), whereas the freezing of the sand changes its strain-hardening behaviour to strain-softening. Supposing that the deviatoric stress–strain curves of the silty clay and sand can be divided into two segments due to a reverse bending point, it was assumed that the first segment follows a hyperbolic function. Meanwhile, the second segment is also a hyperbola, with the reverse bending point as the origin and the residual strength as the asymptote. Accordingly, a nonlinear relation constitutive model that considers σ3 and T is derived. All model parameters are identified. The reasonability of the new model was verified using the test results of the materials. A comparison of the predicted and test results shows that this model can well simulate the deviatoric stress–strain response in the failure process of the tested materials. In particular, it can reflect the residual deviatoric stress after the materials’ failure.


Introduction
The territory covered by frozen soil in China is the third largest in the world and includes 2.15 × 10 6 km 2 of permafrost areas [1]. Many key projects in China, such as railways, highways, and distance pipelines were built in these regions [2]. Among them, the Qinghai-Tibet Railway (QTR) is the highest and longest plateau railway in the world at present [3][4][5]. However, it has been reported that a large amount of frost damage has occurred since the QTR was built [6]. Therefore, it was essential for engineering design and maintenance to systematically investigate the mechanical properties of the frozen remoulded fine-grained materials.
In previous studies, the strength and deformation properties of frozen soil under static conditions considering different influencing factors have been researched. The influence of freezing temperatures (T) on the physical and mechanical characteristics of silty clay has been studied by many researchers [7][8][9][10][11][12][13]. It has been shown that T could significantly improve the physical and mechanical properties of soils due to the formation of a rigid ice-soil matrix [14]. Li et al. [15], Lat et al. [16], and Xu et al. [17] investigated the effects of the confining pressure (σ 3 ), T, and moisture content on the static strength, stiffness,

Study Area and Soil Properties
The silty clay and sand of the remoulded samples was collected within the range of 0.5 m under the shoulder at section K1013 along the Qinghai-Tibet Railway, the highest altitude railway in China's permafrost regions (see Figure 1). The gradation curves of the soil sample are shown in Figure 2. Through laboratory tests, the optimal moisture content of the silty clay was found to be 17.4%, and the corresponding maximum dry density was 1.74 g/cm 3 (see Figure 3). The optimal moisture content of the sand was 11.8%, and the corresponding maximum dry density was 1.88 g/cm 3 (see Figure 3). The two fillings basically accorded with the typical characteristics of the density and optimal water content of silty clay and sand, respectively.

Sample Preparation
According to the Code for Soil Tests of Railway Engineering in China (T [28], the method of sample compaction at different layers was adopted i study. The silty clay and sand obtained from the Qinghai-Tibet railway dried, and sieved, and only the particles under a 2 mm diameter sieve we make the remoulded samples. The samples were divided into five layer where the mass of each layer could be obtained considering 95% of maximu (ρd, max). Subsequently, purified water was added to reach the optimum mo (w). Later on, the prepared soil mixtures were kept in enclosed bags for 24 evaporation. Lastly, the cylindrical material samples with a height of 200 mm eter of 100 mm were compacted layer-by-layer using a standard proctor h compacting one layer, the layer interface was made sufficiently coarse to e

Sample Preparation
According to the Code for Soil Tests of Railway Engineering in China (TB10102-2010) [28], the method of sample compaction at different layers was adopted in the current study. The silty clay and sand obtained from the Qinghai-Tibet Railway was cleaned, dried, and sieved, and only the particles under a 2 mm diameter sieve were collected to make the remoulded samples. The samples were divided into five layers to compact, where the mass of each layer could be obtained considering 95% of maximum dry density (ρ d, max ). Subsequently, purified water was added to reach the optimum moisture content (w). Later on, the prepared soil mixtures were kept in enclosed bags for 24 h to prevent evaporation. Lastly, the cylindrical material samples with a height of 200 mm and a diameter of 100 mm were compacted layer-by-layer using a standard proctor hammer. After compacting one layer, the layer interface was made sufficiently coarse to ensure the two layers were integrated. Then, the specimens were wrapped with rubber sleeves, and the top and bottom were covered with epoxy resin platen to prevent water evaporation.

Test Procedures
The dynamic triaxial tests were conducted on a fully automated Global Digital System (GDS), a cryogenic triaxial apparatus, illustrated in Figure 4. The system is a digitally controlled servo pneumatic system that controls two parameters: axial stress and confining pressure (σ 3 ). The system incorporated a control and data acquisition system, which can maintain an auxiliary air receiver with a servo valve for cell pressure control, local deformation measurement in the vertical direction, and a submersible load cell measuring the applied axial load. The stable confining pressure ranged from 20 kPa to 4 MPa, and the maximum axial load and displacement were 40 kN and 85 mm, respectively.  According to Wang et al. [29], a closed system was adopted for the test program. A one-dimensional freezing-thawing model was employed to simulate the direction of freezing and thawing from the top to the bottom. To ensure one-dimensional freezing and thawing, only the sample top surface was exposed to the external environment. A 50 mm thick layer of insulating polystyrene was used to protect the perimeter and bottom of the cylinder.
The variation range in ground temperature along the Qinghai-Tibet Railway in China is roughly between −20 °C and 20 °C , with a large temperature difference. Therefore, the temperatures under the static loading condition were set to −10 °C , −5 °C , and −1 °C . The time needed to freeze the sample was determined by a pre-experiment test. A short resistance temperature detector probe was inserted at the centre of a pre-experiment test sample during the freezing process to show the variation in temperature. It revealed that 12 h was sufficient to freeze the sample to reach −10 °C with slight fluctuations, plus or minus 0.2 °C .
The test sample was isotropically consolidated at σ3 = 100 kPa, 150 kPa, or 200 kPa after the sample temperature was under the preset temperature for 12 h. After consolidation, the test specimen was then subjected to shearing at an axial strain rate (1.25 mm per min) until the strain of 15% was reached. The particular test scheme is shown in Table 1.  According to Wang et al. [29], a closed system was adopted for the test program. A onedimensional freezing-thawing model was employed to simulate the direction of freezing and thawing from the top to the bottom. To ensure one-dimensional freezing and thawing, only the sample top surface was exposed to the external environment. A 50 mm thick layer of insulating polystyrene was used to protect the perimeter and bottom of the cylinder.
The variation range in ground temperature along the Qinghai-Tibet Railway in China is roughly between −20 • C and 20 • C, with a large temperature difference. Therefore, the temperatures under the static loading condition were set to −10 • C, −5 • C, and −1 • C. The time needed to freeze the sample was determined by a pre-experiment test. A short resistance temperature detector probe was inserted at the centre of a pre-experiment test sample during the freezing process to show the variation in temperature. It revealed that 12 h was sufficient to freeze the sample to reach −10 • C with slight fluctuations, plus or minus 0.2 • C.
The test sample was isotropically consolidated at σ 3 = 100 kPa, 150 kPa, or 200 kPa after the sample temperature was under the preset temperature for 12 h. After consolidation, the test specimen was then subjected to shearing at an axial strain rate (1.25 mm per min) until the strain of 15% was reached. The particular test scheme is shown in Table 1.

Stress-Strain Behaviour
To illustrate the mechanical behaviour of the remoulded fine-grained materials that were influenced by the freezing temperature and σ 3 , Figure 5 shows the typical stress-strain curves for the samples. The results show that the shear strength of the frozen silty clay was significantly greater than that of the frozen sand at the same conditions, especially at low temperatures. It is seen that all the stress-strain curves for silty clay were of the strain-hardening type, and the principal stress deviation had nonlinear growth with the increase in the axial strain, in which the shape of curves was hyperbolic. The stress-strain curves of the silty clay show that the shape of the stress-strain curve did not change with the variation in the temperature, but the degree of hardening was weakened. Decreasing the freezing temperature to −5 • C and −10 • C increased the shear strength of the silty clay by 1374% and 2258% with respect to the silty clay of −1 • C when the σ 3 = 150 kPa.
As for the sand, when the freezing temperature was the same, the stress-strain curves' variation tendency was similar with the different σ 3 . However, the stress-strain curves gradually changed to strain-softening as the freezing temperature decreased. Moreover, the reduction in the temperature from −1 • C to −10 • C caused an increase in the shear strength of 661% for the frozen sand under a σ 3 of 150 kPa. It can be concluded that the influence of freezing on the increase in shear strength of the silty clay was much greater than that for sand.

Cohesion and Internal Friction Angle
The results in Figure 5 were processed using the Mohr-Coulomb strength criterion expressed by the "p-q" method to obtain the internal friction angle ϕ and cohesion c in Figure 6. According to Figure 6, the samples of silty clay with more water content than that of sand strengthened the bite friction during the freezing process; so, the cohesion and angle of the internal friction in silty clay were slightly greater than those in sand. Here, the cohesion increased by approximately 5900% and 3100% for silty clay and sand, respectively, comparing the samples at the temperatures of −1 • C and −10 • C. Similar to the cohesion, the reduction in temperature from −1 • C to −10 • C led to an 11.1% and 8.4% increase in the internal friction angle of silty clay and sand, respectively. It is noted that the rate of increase in the cohesion and angle of internal friction due to the reduction in temperature was almost identical for the silty clay and sand.

Elastic Modulus
The effect of σ3 and freezing temperature reduction on the modulus of elasticity (Ee) of the silty clay is presented in Figure 7a. The values of Ee for the silty clay under the temperature of −1 °C for σ3 values of 100, 150, and 200 kPa were approximately 7.5, 9, and 10 MPa, respectively. After the freezing temperature decreased to −5 °C, the growth in Ee was approximately 900%, 1567%, and 1880%, respectively. As the freezing temperature continued to decrease to −10 °C, the values of Ee increased by 1793%, 1844%, and 2230% with respect to the silty clay at −1 °C for the three σ3 values, respectively.

Elastic Modulus
The effect of σ3 and freezing temperature reduction on the modulus o of the silty clay is presented in Figure 7a. The values of Ee for the silty c temperature of −1 °C for σ3 values of 100, 150, and 200 kPa were approxima

Elastic Modulus
The effect of σ 3 and freezing temperature reduction on the modulus of elasticity (E e ) of the silty clay is presented in Figure 7a. The values of E e for the silty clay under the temperature of −1 • C for σ 3 values of 100, 150, and 200 kPa were approximately 7.5, 9, and 10 MPa, respectively. After the freezing temperature decreased to −5 • C, the growth in E e was approximately 900%, 1567%, and 1880%, respectively. As the freezing temperature continued to decrease to −10 • C, the values of E e increased by 1793%, 1844%, and 2230% with respect to the silty clay at −1 • C for the three σ 3 values, respectively. Moreover, a decrease in freezing temperature led to a significant increase in the Ee of sand. Figure 7b presents the effect of the freezing temperature reduction on the Ee of the sand tested under σ3 values of 100, 150, and 200 kPa at the freezing temperature of −10 °C , −5 °C , and −1 °C . As expected, the Ee of sand generally increased with the decreasing freezing temperature. The growth in the Ee of sand under the freezing temperature of −10°C relative to the sand at −1 °C was approximately 775%, 743%, and 715%, respectively, for the three σ3 values.

Establishment of the Model
As shown in Figure 8, a softening-type stress-strain model with a high degree of accuracy can be used to approximate the nonlinear stress-strain curves for the remoulded fine-grained materials under different freezing temperatures [30]. After careful inspection, the piecewise continuous functions for the proposed constitutive model are expressed as follows: where Dβ, b1, b2, b3, b4, and b5 are undetermined coefficients that can be obtained from the experimental results. The parameters Rσ, Rε, Rσt, and Rεt are defined by Equation Moreover, a decrease in freezing temperature led to a significant increase in the E e of sand. Figure 7b presents the effect of the freezing temperature reduction on the E e of the sand tested under σ 3 values of 100, 150, and 200 kPa at the freezing temperature of −10 • C, −5 • C, and −1 • C. As expected, the E e of sand generally increased with the decreasing freezing temperature. The growth in the E e of sand under the freezing temperature of −10 • C relative to the sand at −1 • C was approximately 775%, 743%, and 715%, respectively, for the three σ 3 values.

Establishment of the Model
As shown in Figure 8, a softening-type stress-strain model with a high degree of accuracy can be used to approximate the nonlinear stress-strain curves for the remoulded fine-grained materials under different freezing temperatures [30]. After careful inspection, the piecewise continuous functions for the proposed constitutive model are expressed as follows: where D β , b 1 , b 2 , b 3 , b 4 , and b 5 are undetermined coefficients that can be obtained from the experimental results. The parameters R σ , R ε , R σt , and R εt are defined by Equation (2a) through (2d) as follows: where (σ 1 − σ 3 ) is the deviator stress; σ 1 and σ 3 are the major principal stress and minor principal stress, respectively; ε is the axial strain; (σ 1 − σ 3 ) p and ε p are the deviator stress and deviator strain, respectively, at the peak point (i.e., point P in Figure 8); and (σ 1 − σ 3 ) t and ε p are the deviator stress and the deviator strain, respectively, at the inflection point (i.e., point T in Figure 8).
where (σ1 − σ3) is the deviator stress; σ1 and σ3 are the major principal stress and minor principal stress, respectively; ε is the axial strain; (σ1 − σ3)p and εp are the deviator stress and deviator strain, respectively, at the peak point (i.e., point P in Figure 8); and (σ1 − σ3)t and εp are the deviator stress and the deviator strain, respectively, at the inflection point (i.e., point T in Figure 8). (1) Coefficient determination for Equation (1a).
To accurately represent the relationship between (σ1 − σ3) and Rε, the terms on both sides of Equation (1a) were multiplied by (σ1 − σ3)p and Equation (1a) was rewritten as follows: The compressive strengths of the fillers are assumed to satisfy the Mohr-Coulomb failure criterion. The relationship between the compressive strength (σ1 − σ3)p and the confining pressure σ3 is described by the following equation: in which cp and φp are the Mohr-Coulomb strength parameters. Next, Equation (4) is substituted for (σ1 − σ3)p in Equation (3) to yield ( ) (1) Coefficient determination for Equation (1a).
To accurately represent the relationship between (σ 1 − σ 3 ) and R ε , the terms on both sides of Equation (1a) were multiplied by (σ 1 − σ 3 ) p and Equation (1a) was rewritten as follows: The compressive strengths of the fillers are assumed to satisfy the Mohr-Coulomb failure criterion. The relationship between the compressive strength (σ 1 − σ 3 ) p and the confining pressure σ 3 is described by the following equation: in which c p and ϕ p are the Mohr-Coulomb strength parameters.
Next, Equation (4) is substituted for (σ 1 − σ 3 ) p in Equation (3) to yield where D β = D × β, with β as a modified coefficient and D = E max /E p . E max and E p (i.e., E p = (σ 1 − σ 3 ) p /ε p ) are the initial tangent modulus at the original point (i.e., point O in Figure 8) and the secant modulus at the peak point, respectively. The determination of E max is equivalent to the determination of E max in the Duncan-Chang model [23]. Performing the differentiation on Equation (1a) produces the following expression In addition, R σ can be considered as a function of the parameter R ε in Equation (1a) and attains its extreme value at the peak. Thus, Combined with Equations (6) and (7), In this case, R σ = 1 and R ε = 1 at the peak point based on their definitions in Equation (2). Thus, these values are substituted into Equation (1a) to obtain According to Equation (1a), R εt and R σt may be expressed at the inflection point as follows: By combining Equations (9) and (10), b 2 and b 3 can be expressed as ., b 1 , b 2 , and b 3 ) are included in Equation (11). When the original b 1 is given, b 2 and b 3 are evaluated using Equations (11a) and (11b). Subsequently, b 2 and b 3 are substituted into Equation (8) to obtain a new value for b 1 . This procedure can be repeated by using the original b 1 until the given and recalculated b 1 are identical.
Similarly, the terms on both sides of Equation (1b) are simultaneously multiplied by (σ 1 − σ 3 ) p . Next, Equation (1b) is combined with Equation (2a) to (2d) as follows: then, Equation (12) is rewritten as The tangent modulus (E) at any point in the first segment of the curve (i.e., R ε ≤ R εt ) can be directly represented as The expressions of R σ and R ε can be rewritten as Equation (16) is substituted into Equation (15) to obtain Equation (6) is substituted into Equation (17), and the tangent modulus in the first segment of the curve can be expressed as In addition, D β = D × β, D = E max /E p , Equation (18) can be rewritten as Therefore, the tangent modulus (E t,1 ) at the inflection point in the first segment of the curve can be obtained from Equation (19).
In addition, the tangent modulus (E t,2 ) at the inflection point in the second segment of the curve can be expressed from Equation (12) by using Equation (21) in terms of the definition of the derivative.
To ensure a continuous condition at the inflection point, the tangent modulus at the inflection point in the two segments of the curve should be equivalent, which indicates that E t,1 = E t,2 at the inflection point. Therefore, Equation (21) is substituted into Equation (20) to obtain Consequently, b 4 can be evaluated when the value of A is substituted into Equation (13a). When the deviatoric stress (σ 1 − σ 3 ) r , which is the residual strength and its corresponding strain ε r at any point in the second segment of the curve, is substituted into Equation (2), R σr and R εr can be evaluated, and b 5 can be obtained using Equation (23). In this study, the residual deviatoric stress (σ 1 − σ 3 ) r is selected when the strain ε r is equivalent to 20%.

Determination of Model Parameters
Six parameters need to be determined from Equation (1) in the proposed model, including the parameters related to the basic elastic properties (D β ), the peak stress state (b 1 , b 2 , and b 3 ), and the residual stress state (b 4 and b 5 ). Most of these model parameters can be conveniently obtained through the monotonic triaxial tests incorporating σ 3 and T. The initial small stress-strain data can be adopted to calculate the D β . The residual stress state parameters (b 4 and b 5 ) can be evaluated through the monotonic triaxial tests results at the residual stress state using Equations (13) and (23). The peak state parameters (b 1 , b 2 , and b 3 ) can be determined by a trial-and-error process using Equations (8) and (11), i.e., by comparing the model predictions and laboratory data. Furthermore, the values of the abovementioned undetermined coefficients with the corresponding values are shown in Table 1.

Model Verification
The test results of the sand and silty clay subjected to different σ 3 and T are used to validate the proposed model. Figure 9 shows the predicted stress-strain of the sand and silty clay compared to the laboratory observations. Note that the model predictions match well with the test data, and the strain-hardening and strain-softening of the samples subjected to different σ 3 and T have been captured very well.

Sensitivity Analysis of the Parameters
In Equation (1a), let b 1 = 0.78, b 2 = −0.44, and b 3 = 0.21; thus, Figure 10a demonstrates the deviatoric stress-strain curves of the proposed constitutive model for different values of the parameter D β . It can be seen from Figure 10a that the change in the value D β has no obvious influence on the shape of the deviatoric stress-strain curve. However, with the increase in the value D β , the initial elastic modulus of the model increases gradually, indicating that the parameter D β mainly reflects the initial elastic modulus of the materials; indeed, the larger the parameter D β is, the higher the initial elastic modulus will be.  b1, b2, and b3, respectively. Figure 10b,d show that the change in the values b1, b2, and b3 has a great influence on the shape of the deviatoric stress-strain curve.  Figure 10b,d show that the change in the values b 1 , b 2 , and b 3 has a great influence on the shape of the deviatoric stress-strain curve. With a decrease in the values b 1 and b 3 and an increase in the value b 2 , the stressstrain curves gradually change from strain-softening to strain-hardening. This indicates that the parameters b 1 , b 2 , and b 3 mainly reflect the peak strength and brittleness degree of the materials' failure. The larger the parameters b 1 and b 3 are, and the smaller the parameter b 2 is, the more obvious the brittleness failure characteristic of the materials will be.
In addition, it can be seen from Figure 10 that the variation in b 4 and b 5 only changes the shape of the residual stress state. With an increase in b 4 and b 5 , the residual stress increases, indicating that the failure characteristic changes from brittleness to ductility.

Comparison with the Existing Model
In order to further illustrate the rationality of the proposed model, the triaxial compression test results of silty clay and sand, with σ 3 = 100 kPa and T = −5 • C, were employed to validate it. Figure 11 presents the comparison of the theoretical curves from the proposed model and the Duncan-Chang model. It can be seen from Figure 11 that the theoretical curves from the Duncan-Chang model had good agreement with the test data at the prepeak region but poor agreement at the postpeak region for the mechanical behaviours of the materials. In particular, it could not reflect the residual deviatoric stress after the materials' failure, while the theoretical curves from the model proposed in the current study conform well with the test data in both the prepeak region and postpeak region, no matter what the shape of the stress-strain curves. Therefore, the proposed model has better adaptability. curves gradually change from strain-softening to strain-hardening. This indicates that the parameters b1, b2, and b3 mainly reflect the peak strength and brittleness degree of the materials' failure. The larger the parameters b1 and b3 are, and the smaller the parameter b2 is, the more obvious the brittleness failure characteristic of the materials will be. In addition, it can be seen from Figure 10 that the variation in b4 and b5 only changes the shape of the residual stress state. With an increase in b4 and b5, the residual stress increases, indicating that the failure characteristic changes from brittleness to ductility.

Comparison with the Existing Model
In order to further illustrate the rationality of the proposed model, the triaxial compression test results of silty clay and sand, with σ3 = 100 kPa and T = −5 °C , were employed to validate it. Figure 11 presents the comparison of the theoretical curves from the proposed model and the Duncan-Chang model. It can be seen from Figure 11 that the theoretical curves from the Duncan-Chang model had good agreement with the test data at the prepeak region but poor agreement at the postpeak region for the mechanical behaviours of the materials. In particular, it could not reflect the residual deviatoric stress after the materials' failure, while the theoretical curves from the model proposed in the current study conform well with the test data in both the prepeak region and postpeak region, no matter what the shape of the stress-strain curves. Therefore, the proposed model has better adaptability.

Conclusions
This paper presented the results of a comprehensive experimental investigation to study the effect of confining pressure (σ3) and freezing temperature (T) on the mechanical behaviour of remoulded fine-grained materials (silty clay and sand) used under the Qinghai-Tibet Railway line. Accordingly, a nonlinear constitutive model, which could reflect the postpeak softening behaviour was established. Conventional monotonic triaxial tests incorporating σ3 and T carried out in the current study were used to evaluate all the model parameters. The salient outcomes of the model are summarized below.
(1) All the silty clay exhibited a strain-hardening type of stress-strain curve, but the sand under the temperatures of −5 °C and −10 °C showed strain-softening. Under the same test conditions, the shear strength of the silty clay was greater than that of the sand. In all cases, the cohesion (c) and angle of internal friction (φ) of the silty clay were Figure 11. Comparison of the predicted and tested data of the tested remoulded fine-grained materials under T = −5 • C and σ 3 = 100 kPa.

Conclusions
This paper presented the results of a comprehensive experimental investigation to study the effect of confining pressure (σ 3 ) and freezing temperature (T) on the mechanical behaviour of remoulded fine-grained materials (silty clay and sand) used under the Qinghai-Tibet Railway line. Accordingly, a nonlinear constitutive model, which could reflect the postpeak softening behaviour was established. Conventional monotonic triaxial tests incorporating σ 3 and T carried out in the current study were used to evaluate all the model parameters. The salient outcomes of the model are summarized below.
(1) All the silty clay exhibited a strain-hardening type of stress-strain curve, but the sand under the temperatures of −5 • C and −10 • C showed strain-softening. Under the same test conditions, the shear strength of the silty clay was greater than that of the sand. In all cases, the cohesion (c) and angle of internal friction (ϕ) of the silty clay were greater than that of sand. Furthermore, the modulus of elasticity of the materials tested increased due to freezing and the temperature reduction. (2) A practical constitutive model was developed to represent the nonlinear, stressdependent, and inelastic stress-strain behaviours of the fillers subjected to freezing and thawing. This model incorporated three important aspects of the stress-strain behaviour, including nonlinearity, strain-dependency softening, and inelasticity. A simple technique was used to interpret the test results and conveniently determine the six parameters in the model. (3) The triaxial test results of the remoulded fine-grained materials were employed to evaluate the reasonability of the proposed model established in this paper. A comparison of the predicted and test results showed that this model could well simulate the deviatoric stress-strain response in the failure process of the tested materials.
In particular, it could reflect the residual deviatoric stress after materials' failure. (4) This study analysed the behaviour of the fillers with optimum water content that were exposed to the freeze-thaw cycles to develop a constitutive model. If appropriate experimental results are available, the parameter values in the proposed model can be derived from the triaxial test results. Therefore, additional experiments should be conducted to investigate other parameters, such as the temperature, duration of freezing and thawing, and the moisture content and compactness of the fillers, which are important characteristics of fillers in regions that are subjected to seasonal freezing and thawing.