Experimental Testing and Numerical Modelling of Mechanical Behaviors of Silty Clay under Freezing-Thawing Cycles

Seasonal freezing-thawing cycle is one of the most common physical weathering processes in cold regions, which can significantly affect the mechanical behaviors of soil. In this paper, a series of freezing-thawing (F-T) cycle and consolidated drained triaxial tests have been carried out on silty clay samples collected fromTibetan Plateau. To do so, amodified numericalmodel was developed taking into F-T effect. Test results showed that the stress-strain curves of original soil specimens presented strain hardening characteristics, accompanied with shear shrinkage. In F-Texperienced specimens, volumetric strain in triaxial loading stage was gradually increased, while failure strength was decreased. Elliptic and parabolic functions were selected in numerical modelling to describe volume and shear yield surfaces on a p-q plane, respectively. Moreover, a double-yield surface constitutive model was developed to describe relationships among deviatoric stress, axial strain, and volumetric strain. Furthermore, equations for model parameters with the number of F-Tcycles as variable were derived based on the triaxial test results which were then substituted into the established model to take into account the effects of F-T cycles. Finally, numerical results were validated with experimental findings.


Introduction
In cold regions, ground bases, and soil foundations such as subgrades, generally suffer seasonal and diurnal freezingthawing (F-T) cycles due to atmospheric temperature changes. F-T process causes soil structure degradation, strength reduction, settlement, etc. Numerical constitutive models can quantitatively predict the mechanical behaviors of soil [1]. erefore, numerical modelling of soil mechanical behaviors is one of the main geotechnical engineering research issues in cold regions.
Currently, soil constitutive models mainly include nonlinear elastic and elastoplastic types. Nonlinear elastic models include Duncan E-v, E-B, and K-G models, none of which can describe soil dilatancy and the influences of both compress and shear loadings [2]. Elastoplastic models mainly refer to the Cam-clay model which was first proposed by Roscoe and Schofield [3] and a series of modified single-yield surface models developed on the basis of the original one. Among them, Asaoka et al. [4] established a modified model capable of reflecting historical stress path and soil structure characteristics, by embedding super-consolidation and structural parameters into the Cam-clay model yield function to form upper and lower loading yield surfaces in the p-q plain. Based on state-related dilatancy theory, Li and Dafalias [5] suggested several parameters for the modification of the Cam-clay model to make it capable of describing shear contraction and dilatancy characteristics of soil. Moreover, Yin and Graham [6] developed a viscoelastic-plastic constitutive model considering time effect based on equivalent time concept. Furthermore, Huang et al. [7] established the Tsinghua model, and Li [8] improved it. ey also experimentally determined the direction of plastic strain increment under each stress state and found yield surface according to the associated flow law which is only appropriate for soil strain hardening process. Yao [9] proposed UH models suitable for both super-and normal consolidated soils, based on the Cam-clay model and lower loading surface concept. To sum up, the single-yield surface model is a "hat" model. erefore, it can only reflect soil shear shrinkage. However, there may still be yielding phenomena within yield surfaces in stress space.
To address the shortcomings of single-yield surface models, two-yield and multiple-yield surface models were developed. Lade [10] first elaborated two-yield surface concept and further defined volumetric and shear yield surfaces as failure and plastic expansion yield surfaces, respectively; the latter one being a bullet-shaped cone. Shen [11] combined the advantages of Duncan-Chang and Cam-clay models and proposed a two-yield surface model called Nanshui model, which could simulate the stress path characteristics of decreasing confining or average pressures. Zheng and Chen [12] developed multiple-yield surfaces based on the vector principle and introduced the concepts of volume yield surface, shear yield surface, and shear yield surface of θσ direction. Yin et al. [13,14] proposed yield and hardening functions relating to compress and shear loadings, respectively, and then established another two-yield surface model using correlated flow law. Huang et al. [15] introduced a shape parameter into the Cam-clay model and developed a two-yield surface model which was able to reflect soil strain softening. Several other yield surface models and boundary theory-based models have also been developed [16].
At present, one of the popular research fields is to improve models to make them applicable for different soil properties, environments, and boundary conditions based on the above models [17]. Particularly, for soil constructions in cold regions, F-T process of soil causes drastic changes in its physical and mechanical characteristics. Many forms of viscoelastoplastic and elastoplastic models [18,19] have been developed based on the mechanical behaviors of frozen soil. Meanwhile, the thawing of frozen soil is a major source of engineering problems. Regarding soil thawing, Liu et al. [20] proposed two new factors based on the Duncan model, namely, deviatoric stress reduction ratio and the number of F-T cycles, to predict the stress-strain relationship of soil under F-T cycles. Cui et al. [21] fitted volume and shear hardening functions with logarithmic and exponential functions, respectively, and evaluated F-T effect by describing the variation law of hardening function coefficient. Chang et al. [22,23] established a two-yield surface model considering F-T effect by fitting the equations of hardening function coefficients as a cubic polynomial with the number of F-T cycles. e above studies showed that soil F-T effect could be simply determined by selecting reasonable yield surfaces and equations of functional coefficients based on typical two-yield surface constitutive models. However, the applicability of this model considering F-T effect remains to be further analyzed.
In this paper, a series of F-T cycle and consolidated drained triaxial tests were carried out on the silty clay samples collected from Tibetan Plateau, where soil ground suffers continuous F-T cycles. A two-yield surface constitutive model was first established based on the conventional Nashui Model. en, model parameter equations by considering the number of F-T cycles as a variable were derived according to triaxial test results. Finally, a modified model was proposed and validated to predict the deviatoric stress, axial strain, and volumetric strain of soil under F-T cycles.

Soil Materials.
e soil samples employed in the experiments were collected from the newly constructed Gonghe-Yushu Expressway, which locates in Tibetan Plateau, China, and is the first expressway in the world to be constructed in a permafrost region, as shown in Figure 1.
e collected soil samples were used as a subgrade filler for expressway and showed obvious F-T effects on both strength and deformation behaviors. e grain composition of the soil sample is shown in Figure 2. Soil was turned over, dried, and passed through a 2 mm sieve before testing, according to the preparation requirements of soil samples as advised in reference [24]. e soil sample had maximum dry density of 1.828 g/cm 3 , optimal moisture content of 14.8%, relative density of 2.64, liquid limit of 28.0%, and plasticity index of 10.3. Based on the above values, the soil sample was specified as silty clay.

Preparation of Specimens and Test
Scheme. According to subgrade compaction requirements, the specimen preparation standards for moisture content and compaction degree were considered to be 14.8% and 95%, respectively. Figure 3 shows the details of the simulation of subgrade fillers in real operational conditions. Cylinder specimens with diameter 39.1 mm and height 80 mm were prepared through the layered compaction method. First, F-T cycle tests were performed in an incubator with alternate constant temperatures of − 5°C and 20°C with a lasting time of 12 hours at each temperature to ensure the complete freezing and thawing of specimens, as shown in Figure 3(a). F-Texperienced specimens are shown in Figure 3(b), which were covered with plastic layers to prevent water evaporation or absorption. en, in order to test the pore water pressure accurately to further study the mechanical and deformation characteristics of tested soils and to obtain the potential most unfavorable results, the thawed specimens were taken out of the incubator for vacuum saturation in a sealed vessel, as shown in Figure 3(c). Finally, F-T-experienced saturated specimens underwent consolidated drained triaxial tests on a triaxial loading apparatus made by GDS Instrument Company, UK, as shown in Figure 3(d). Specifically, the F-Texperienced sample was further saturated through a backpressure saturation process, which was conducted by a stepwise (20 kPa) method and continued until the ratio of pore water pressure increment to confining pressure increment was larger than 98%. e drain valve of the GDS apparatus was then opened, and the consolidation process was carried out under a predetermined confining pressure, until the pore water pressure dissipated by more than 95%. At last, the loading test started with a constant axial shearing rate of 0.01%/min, and it was continued until the ultimate axial strain of each sample reached 20%. All the above testing methods are referred to the provisions of reference [24].
In this study, the number of F-T cycles and confining pressure were selected as test variables in F-T cycle and consolidated drained triaxial tests, respectively. F-T cycle numbers were selected as 0, 1, 3, 6, 9, and 12 times, and confining pressures for both consolidating and triaxial loading process were set at 100, 200, 300, and 400 kPa. e axial strain rate was set at 0.01%/min, and the test was terminated when the axial strain reached 20%. A total of 24 specimens were tested.

Volume Change of Specimens Subjected to F-T Cycles.
e volume change of specimens with the number of F-Tcycles is illustrated in Figure 4. It can be seen that the specimen volume was increased by increasing the number of F-T cycles. e ultimate volume increasing rate can be 3.63% at the F-T cycle numbers of 12. It is because the volume of liquid water in soil pores increases by 9% after freezing [25]. Pore water freezing process destroys the original soil structure and the connection of soil particles which cannot be restored during subsequent thawing. As a result, continuous degradation on soil mechanical behaviors was observed in specimens.

Deviatoric Stress-Axial Strain Characteristics.
Deviatoric stress-axial strain curves under different test conditions are shown in Figure 5. It can be seen that all stress-strain curves showed strain hardening and the deviatoric stress was gradually increased with the increase of axial strain. In addition, at low confining pressures, stressstrain curves showed certain weak strain softening characteristics. Meanwhile, the spatial position of the stress-strain curve was gradually decreased with the increase of F-T cycle numbers.
From the perspective of practice, it is feasible to use the peak value of the deviatoric stress as the failure criterion. For simplicity, the failure strength of each specimen was obtained based on its deviatoric stress value corresponding to axial strain of 15% as advised in reference [24]. e change rate of failure strength is displayed in Figure 6. It can be seen in the figure that failure strength was decreased with the increase of the number of F-T cycles. Meanwhile, the change rate of failure strength was decreased by increasing confining pressure, which indicated that confining pressure could weaken F-T effect.

Volumetric Strain-Axial Strain Characteristics.
Consolidated drained triaxial test results for pore water drainage volume after different F-Tcycle numbers are shown   Mathematical Problems in Engineering 3 in Table 1. It can be seen that drainage volume was increased by increasing both F-T cycle number and confining pressure. It can be seen from Figure 6 and Table 1 that confining pressure at both triaxial loading stage and consolidation stage is beneficial to reduce the structural deterioration of the specimen in the F-T cycling test. Meanwhile, drainage volume increase due to F-T cycles revealed that the unfavorable degradation of soil structure was increased with the number of F-T cycles. Nonetheless, soil structure was recovered to some extent during consolidating process, and recovery degree was increased by confining pressure. e volumetric strain-axial strain relationships of the loading stage under different test conditions are shown in Figure 7. It can be seen that specimens basically showed shear shrinkage property, and volumetric strain was gradually increased by increasing axial strain under similar conditions. Except for non-F-T-experienced specimens, they showed a dilatation behavior in the early loading stage under the confining pressure of 100 kPa. Meanwhile, with the increase of F-T cycle numbers and confining pressures, volumetric strain corresponding to similar axial strain was also increased. However, at a confining pressure of 400 kPa, the volumetric strain corresponding to the same axial strain showed an opposite decreasing trend compared with that at 300 kPa, which was because of more significant consolidation drainage under the former test condition.

Two-Yield Surface Constitutive Model
In elliptic-parabolic two-yield surface constitutive model for soil [13,14], yield surface is regarded as the boundary of elastic deformation regions which are composed of elliptic and power functions. Soil deformation consists of elastic dε e   Figure 8, f 1 and f 2 represent the yield locus and N is the current point. Hence, the p-q plane is divided into four regions of A 0 , A 1 , A 2 , and A 3 , where A 0 is only characterized by elastic strain, A 1 is related to plastic compression yield, A 2 is related to plastic expansion yield, and A 3 is related to both compression and expansion yields simultaneously. At point N, plastic increment direction was the vector sum of different plastic strain increments. As a result, the composition of total strain dε was

Yield Associated with Compression and Expansion.
e yield surface corresponding to compression adopts elliptic shape in the modified Cam-clay model [2] and improved yield equation f 1 is still in the form of the elliptic function.
e main equations and parameters were expressed as follows [13,14]: where M 1 presents the elliptic shape change, p is the mean normal stress, q is the generalized shear stress, ε p-q plane, and H a1 is the hardening function of compression yield surface. e calculation method for p r was stated as Under isostatic triaxial testing conditions, the relationship between p 0 and ε p v could be expressed in a hyperbolic form. erefore, the hardening function corresponding to the first yield surface was written as [13,14] where p 0 is the horizontal coordinate of the intersection point of yield locus f 1 and axis p and h and t are the dimensionless parameters to describe the relationship between confining pressure and plastic volumetric strain. Yield equation f 2 corresponding to shear expansion was a parabolic function expressed as [13,14] where a represents the strain ratio of yield surface f 2 to total plastic strain; therefore, high values of a correspond to dilation. G is the elastic shear modulus, ε p s is the hardening parameter of parabolic yield surface, and H a2 is the hardening function of expansion yield surface.
G was calculated as where for the elastic modulus, we had E = 2.0E i . Generally, the proportion of elastic strain in the total strain is small; therefore, v is always set at 0.3. e initial tangent modulus E i is related to confining pressure, and its empirical equation was derived as where p a is the standard atmospheric pressure, which is 101.3 kPa, and k and n are the intercept and slope of the line lg(E i /p a )∼lg(σ 3 /p a ). e hardening function corresponding to the second yield surface was obtained to be [13,14]

Constitutive Equations in Increment
Form. e three strain components of equation (1) where K is the bulk modulus. By differentiating both sides of equation (2), we found Assuming both yield surfaces to be applicable to associated flow law, the following equation was obtained:      was obtained by combining equations (11) and (12). e plastic compression components of the second yield surface included volumetric elastic strain dε p2 v and shear strain dε p2 s increments. Using the same derivation method as the first yield surface, the following equation was obtained for dε p2 s : . (14) e calculation equation of dε p2 v was obtained by combining equations (13) and (14).
According to typical elastic-plastic theory, dε v and dε s were defined as where A, B, C, and D were expressed as follows:

Equations of Triaxial Tests.
In conventional consolidated drained triaxial tests, stress state is generally in A 3 region, as illustrated in Figure 8. e relationships between p and q and increments dp and dq are given by equations (20) and (21), respectively: In strain-controlled triaxial tests, the relationship of axial dε a and volumetric dε v strain increments were By substituting equation (21) into equation (15), the volumetric strain increment was found to be By substituting equations (22) and (23) into equation (15), axial strain increment was derived as

Mathematical Problems in Engineering
Finally, the incremental constitutive equation of the relationships among deviatoric stress, axial strain, and volumetric strain of soil were obtained by combining equations (16)∼ (19), (23), and (24). ere are 9 parameters in the above numerical model including c, φ, M 1 , h, t, k, n, M 2 , and a.

Functional Expressions for Model Parameters
According to the test results of consolidated drained triaxial tests, the model parameter variation rules for the number of F-T cycles N and corresponding function expressions were determined, namely, c(N), φ(N), M 1 (N), h(N), t(N), k(N), n(N), M 2 (N), and a(N).

Determination of c(N) and φ(N).
e internal friction angle φ and cohesion c were calculated according to the slope and intercept of the tangent line of Mohr's stress circle as where α is the slope angle and d is the intercept of the envelope. e tangent lines of Mohr's stress circles for different F-T cycle numbers are shown in Figure 9, and shear strength parameters are summarized in Table 2. It was seen that both c and φ were decreased with the number of F-T cycles. e logistic function, as given by equation (26), was used to fit the regression relationships among φ, c, and number of F-T cycles, as shown in Figure 10: where x 0 is the shear strength parameters of non-F-T-experienced specimens and A 1 ∼A 3 are the fitting parameters.

Determination of M 1 (N).
e empirical equation of M 1 (N) is derived as [13,14] where β is the ratio of volumetric strain to axial strain at stress level 75% and β values vary with confining pressure; therefore, its average value was adopted. M is the slope of failure line as shown in Figure 8, which can be obtained by equation (28) combined with triaxial compression test results: e values of β, M, and M 1 are defined in Table 3. It can be seen that β was gradually increased with the number of F-T cycles, while M and M 1 were decreased. e regression relationship between M 1 and the number of F-T cycles was described by the logistic function, as shown in Figure 11.

Determination of h(N) and t(N). h(N) and t(N)
were obtained through modifying equation (4) [13,14]: When B p � dp 0 /dε According to the previous experiences, plastic volumetric strain ε p v accounts for about 65% of volumetric strain by increasing the stress level from 0 to 50%. Based on equation (2), p 0 was obtained through p and q at the stress level of 50%. en, B p was calculated according to ε p v and p 0 , specifically. e intercept and slope of the curve of (B p /p a ) 0.5 ∼p 0 /p a were h 0.5 and t/h 0.5 , respectively; therefore, h and t were also obtained. Curves of (B p /p a ) 0.5 ∼p 0 /p a for different F-T cycle numbers are shown in Figure 12, and h and t values are given in Table 4. It was seen that h was gradually decreased with the number of F-Tcycles while t was increased. e relationships among h, t, and the number of F-T cycles were fitted by the logistic function, as shown in Figure 13.

Determination of k(N) and n(N).
e curve of lg (σ 3 /p a )∼ lg (E i /p a ) based on triaxial test results is shown in Figure 14. It was seen that the spatial positions of fitting lines moved downwards by increasing the number of F-T cycles. e values of intercept k and slope n are given in Table 5. As can be seen, k was decreased with the number of F-T cycles while n was increased. e relationships among k, n, and the number of F-T cycles were described by the logistic function, as illustrated in Figure 15. (N) and a(N). M 2 and a were determined by transforming equations (5) and (8) [13,14]:

Determination of M 2
In equation (31), ε p s was obtained using the following empirical equation: where d is the slope of the curve of ε v ∼ε a at a stress levels of 75%∼95%. e above valuing method of d is actually an   empirical one, and its reliability has been verified [13,14]. Negative values of d meant that there was dilation. e values of d in each test curve were different, and therefore, its mean values were adopted. e intercept of the curve of (p + p r )/q∼(q/ε p s G) 2 was 1/M 2 , and its slope was a 2 /M 2 ; therefore, M 2 and a were also obtained. e values of M 2 and a for different F-T cycle numbers are summarized in Table 6. It was seen that M 2 and a were decreased with the number of F-T cycles and the relationships among M 2 , a, and the number of F-T cycles were described by the logistic function, as shown in Figure 16.
Based on the two-yield surface model described in Section 4, model parameters were fitted as functional expressions with the number of F-T cycles N as variable. A modified model was then developed considering the influences of confining pressure and F-T cycles. Equations (23) and (24)

Validation of the Modified Two-Yield
Surface Model e above incremental constitutive equations were calculated using a spreadsheet. Firstly, according to strain-controlled triaxial test results, equation (24) was used to calculate deviator stress increment dq corresponding to a constant axial strain increment value of 0.125%. en, according to the obtained dq and equation (23), the volumetric strain corresponding to each level of axial strain was obtained.
e relationships among deviatoric stress, axial strain, and volumetric strain under different conditions were obtained through the calculation program shown in Figure 17. Among them, the results obtained after 4 F-T cycles were a supplementary verification test. It was seen that numerical results were basically consistent with experimental findings. For the same axial strain, deviatoric stress was decreased with the number of F-T cycles, while the volumetric strain was increased. It was demonstrated that the modified     Mathematical Problems in Engineering numerical model could predict the variations of deformation and strength characteristics of silty clays exposed to F-T cycles.

Conclusions
(1) e volume of silty clay was increased with the number of F-T cycles, which degraded soil mechanical behaviors to some extent. In the consolidated drained triaxial test, the discharge amount of pore water in the consolidating stage was increased with the number of F-T cycles. In the triaxial loading stage, specimens basically showed a shear contraction property. Meanwhile, volumetric strain was increased with the number of F-T cycles. Confining pressure in both consolidating and triaxial loading stages weakened the unfavorable effects of F-T cycles. (2) e incremental constitutive equations of deviatoric stress, axial strain, and volumetric strain increments have been derived in this paper by introducing the associated flow law to a typical two-yield surface model. Under the effect of F-T cycles, model parameters c, φ, h, k, M 1 , M 2 , and a were decreased while t and n were increased. Moreover, the logistic function was used to fit the regression relationship between the above model parameters and the number of F-T cycles. Furthermore, a modified numerical model considering F-T effects was developed by substituting model parameter expressions to incremental constitutive equations. (3) A program for the calculation of modified constitutive equations was developed. e comparison between numerical and test results showed that the proposed modified two-yield surface model well predicted the effect of F-T cycles on soil mechanical behaviors.

Data Availability
All data, models, and code generated or used during the study appear in the submitted article.

Conflicts of Interest
e authors declare no conflicts of interest.