Study on triaxial test and constitutive prediction model of frozen silty clay

With the increasing demand for engineering construction in the seasonal frozen area and the background of the Belt and Road Initiative, the frozen soil constitutive model should be studied in depth. At present, the constitutive prediction model of frozen silty clay has many problems, such as complex formula, single model application and poor prediction ability. Random forest optimal model hyperparameter input was very difficult. Particle Swarm Optimization (PSO) was used to optimize the parameters of the number of neurons, dropout and batch_size in the Long-term and Short-Term Memory network (LSTM) structure. The optimization results were 61, 0.09 and 95 respectively. The results showed that the strength tended to be stable after 6,9,6,9 and 9 freeze-thaw cycles under initial moisture content = 25, 22.5, 20, 17.5, and 15%, respectively. After 18 freeze-thaw cycles, the strength decreased by 2.66%, 11.85%, 18.83%, 16.79, and 29.02%, respectively. The predicted values of frozen soil binary medium model (BM), random forest model (RF) and PSO-LSTM model were compared with the measured values under different working conditions, and good accuracy was obtained. The R2 of the PSO-LSTM model test set was trained to more than 98%, and RMSE, MAE and MAPE were also trained to the lowest under the same working conditions. The influencing factors of deviator stress of frozen silty clay were given in order from strong to weak: initial moisture content>strain>confining pressure>number of freeze-thaw cycles. The LSTM optimal combination input parameters were searched by PSO, and the parameter adjustment speed of the model for the data learning process of frozen silty clay was greatly increased, which was conducive to the promotion of other soil constitutive prediction models. A new constitutive prediction model of frozen silty clay was developed using PSO-LSTM algorithm. 15 working conditions had been verified, and the optimal model had high accuracy in the constitutive prediction of frozen silty clay, which provided a good reference for the application of frozen soil engineering in cold regions.


Introduction
The world's cold regions can be divided by air temperature, snow thickness and frozen soil depth. The distribution area of global frozen soil accounts for about 24% of the total land. Seasonal and intermittent frozen soil regions account for 50.6% and 6.6% of the land area in the Northern Hemisphere, respectively (Zhang et al., 2003), indicating that more than half of the land in the Northern Hemisphere experiences a freeze-thaw cycle every year (Li et al., 2021). The freezing damage of soil around building foundation, road foundation and pipeline is controlled by soil mechanical properties under freezing state. Frozen silty clay is widely distributed in cold regions of the Northern Hemisphere. In order to study the engineering properties of soil in cold regions, it is of great significance to study the constitutive relationship of frozen silty clay. Constitutive model is the key to analyze the stress-strain relationship of frozen soil, and also the basis for studying the stability of landslide Li et al., 2022).
Soil elastic-plastic definite solution needs to be solved simultaneously with equilibrium differential equation, geometric equation, deformation coordination equation and known boundary. The constitutive model, as the deformation coordination equation, is an important theoretical basis for analyzing the soil stress-strain relationship. So far the most studied soil constitutive models can be divided into: elastic-plastic constitutive model (Zhang and Liu, 2019;Lai et al., 2010;Wang et al., 2020;Yu et al., 2022), nonlinear constitutive model (Yang et al., 2013;Loria et al., 2017), critical state constitutive model (Wu et al., 1996;Yao et al., 2004;Yang et al., 2010;Yin et al., 2010;Lai et al., 2016;Jin et al., 2017;Sun et al., 2020), hypoplastic constitutive model (Lai et al., 2014;Mašín, 2015;Chang et al., 2019) and the micromechanical constitutive model (Yin et al., 2010;Yin et al., 2014;Xiong et al., 2017;Qu et al., 2021;Qu et al., 2021;. However, there are some problems in the above soil constitutive model: 1) the stress-strain relationship between frozen soil and thawed soil is proposed based on certain assumptions (Zhang et al., 2022;Yin et al., 2010). 2) Frozen soil and thawed soil models are generally only applicable to one soil type . 3) The mathematical formula of stress and strain of frozen soil and thawed soil comes from data fitting, but the model prediction effect is not good (Yin et al., 2011). 4) With the increase of fitting parameters (Yin and Jin, 2021), the constitutive formula of frozen soil becomes increasingly complex.
The stiffness of frozen soil is between soft soil and rock, so a binary medium model can be used to describe the freeze-thawcompression failure mechanism and stress evolution of frozen silty clay. In terms of the research on the frozen soil binary medium model and freeze-thaw damage mechanism, few articles have done relevant analysis.  verified the data obtained from the triaxial test of frozen soil through homogenization theory and binary medium theory.
In recent years, frozen soil constitutive theory has been being studied by some methods, such as experimental research, theoretical derivation, formula fitting, numerical simulation of particle flow, and data-driven deep learning. Numerous ML algorithms have been used to study thawing and frozen soil constitutive, such as Evolutionary Polynomial Regression (EPR) (Nassr et al., 2018), Support Vector Machine (SVM) Kohestani and Hassanlourad, 2016), Back Propagation Neural Network (BPNN) (Shahin and Indraratna, 2006;Johari et al., 2011;Rashidian and Hassanlourad, 2014;Stefanos and Gyan, 2015;Lin et al., 2019), radial basis function (RBF) neural network (Peng et al., 2008), recurrent neural network (RNN) (Zhu et al., 1998;Romo et al., 2001), long short-term memory (LSTM) neural network . Problems such as gradient explosion or gradient disappearance can be better avoided by LSTM neural network. However, the selection of the optimal model input parameter combination is still an urgent problem to be studied. However, for different frozen silty clay datasets, the optimal neuron numbers, dropout probability and batch_size combination corresponding to the minimum stress loss function should be different. Aiming at the problems of difficult parameter adjustment and poor fitting accuracy, a method based on long short-term memory neural network (LSTM) and particle swarm optimization (PSO) is proposed to predict the stress-strain evolution law of frozen silty clay.
The structure of this paper is as follows. In the first part of the experimental study, the stress-strain evolution process of frozen silty clay is qualitatively analyzed by mechanism. In the second part, the binary medium frozen silty clay constitutive model based on the simplification of freeze-thaw damage parameters is studied, and the fitted mathematical formula provides a reference for engineering design. The third part, frozen silty clay constitutive prediction model principle is introduced. In the fourth part, the relationship between the RF model, the frozen soil binary medium model and the PSO-LSTM model and the generalization value of the optimal model are discussed, and the mechanism of the influence of the three factors on the strength of frozen silty clay is deeply analyzed. Finally, this paper summarizes and puts forward the future work plan.

Test overview
The test instrument adopts the frozen soil dynamic and static triaxial tester FST-250 designed by Xi'an Lichuang Material Testing Technology Co., Ltd. (see Figure 1). The instrument adopts displacement control mode and is mainly composed of triaxial pressure system, XT5718ULT-E2000 cooling water circulating device, data acquisition system, cooling system of testing machine, etc. The temperature control range is  -70°C~+90°C, the temperature control accuracy is 0.5°C, the confining pressure range is 0～40MPa, the accuracy is 1%, and the maximum axial pressure is 250KN.
The microcomputer-controlled electro-hydraulic servo triaxial testing machine for frozen soil is mainly used for dynamic and static triaxial tests of conventional soil and frozen soil. The basic physical properties of the test soil are shown in Table 1. In order to reduce the heterogeneity and dispersion of natural silty clay, 1.25 mm sieve is used to screen out impurities such as plant fiber and gravel of natural silty clay. The distribution of soil with particle size less than 1.25 mm is shown in Table 2. Direct shear test is used to measure cohesion and internal friction angle of soil samples. The liquid limit and plastic limit in Table 1 are measured by LP-100D digital display soil liquid and plastic limit tester. The standard consolidation method is adopted in the consolidation test, and the compressibility factor and modulus are obtained. The maximum and minimum values are obtained from experiments, and the average values are obtained from calculation. According to the liquid plastic limit and particle size distribution, the test soil can be classified as silty clay.

Soil preparation
The dry soil that has passed through 1 mm sieve is sprayed and mixed with soil according to the initial moisture content (15%, 17.5%, 20%, 22.5%, and 25%). Then the soil samples are sealed with plastic wrap for 24 h. The dry density of soil sample is controlled to 1.62 g/ cm3. Soil samples with different water contents are divided into five layers and compacted into a three-part mold with an inner diameter of 61.8 mm and a height of 125 mm. After standing for a certain time, soil samples with different initial moisture contents are demolded. Then, the soil samples are packaged and labelled with plastic wrap. Finally, the soil samples are put into the high and low temperature alternating test chamber TEMP880.

Freeze-thaw cycles test
The high and low temperature alternating test chamber TEMP880 designed by Shanghai-Huayi Equipment Co., Ltd. is used for the freeze-thaw cycle test. The test chamber is set at -35°C-35°C. The specific temperature time history curve in the high and low temperature alternating test chamber during the experiment is shown in Figure 2 below. The soil samples undergoe 4 h of freezing and 4 h of thawing, that is, a freeze-thaw cycle is performed every 8 h . Then, the soil samples are packaged with a thin latex sleeve, and place in a constant temperature test chamber with the same design temperature for 24 h after sticking the label. YH-10 aviation hydraulic oil in the triaxial pressure chamber is used to apply confining pressure and maintain temperature stability. Next, the soil samples with sufficient temperature control in the constant temperature test chamber are placed in the triaxial pressure chamber. When the temperature of the soil sample and hydraulic oil is coupled to the design temperature, the design confining pressure is applied for the test.

Time scaling principle of freeze-thaw cycle
Because the time scale of soil freezing and thawing process in seasonally frozen soil area is considerably large, the reduced similarity ratio is adopted in the test. The similarity criteria based on the similarity theoretical model are as follows (Guo et al., 2010): Where L represents the characteristic length, m. Γ represents the period of surface temperature change, h a represents the proportional coefficient.

FIGURE 2
Temperature variation of soil sample.
Frontiers in Earth Science frontiersin.org 04 The physical simulation characteristics of frozen silty clay are consistent with the test process and the actual process. The experiment lasted 8 h (4 h of freezing and 4 h of thawing), and the 1-year change process under the action of freezing and thawing cycles was simulated. The similarity calculation process is shown in Table 3:

Shear test
Constant deformation loading is adopted during the test, and the shear rate is 1 mm·min-1. The test is terminated when the specimen strain reached 20% or when it failed. At -15°C, Triaxial tests are carried out on frozen silty clay nuder five freeze-thaw cycles, five initial moisture contents and three confining pressures. The change of soil sample stress is monitored in real time and recorded. After the triaxial test, the frozen soil sample is convex in the upper part and has expansion deformation in the radial direction. The frozen silty clay generally presents an oblique crack failure form (Xu et al., 2020;Niu et al., 2022;Zhang et al., 2022).

Long short-term memory neural network structure and implementation
Long short-term memory neural network (LSTM) was created primarily to solve the problem of gradient disappearance and explosion. The biggest difference between LSTM and traditional neural network is the introduction of Forget gate, Input gate and Output gate. Another key is the change in the internal state, which is the core of each activated neuron and can be seen as a carrier for adding information or deleting information.
Step 1: Some information is selectively discarded by the forget gate, and LSTM receives the input of C t and h t−1 . And output a number between 0 and 1 to the internal state S t .
Step 2: Some information in the internal state is determined to be stored. The input gate determines which values are updated, and then a candidate state S (t) is built by the input node g (t) .
The original internal state is updated to S (t) by Eqs 2-4.
Step 3: The output information is determined by the output gate.
Step 4: Based on the output gate and the new internal state, the new hidden layer state h (t) is output.
where W and b represent the weight and bias, respectively, W fC is the weight of C (t) in the forget gate, and b f is the forget gate bias.

Introducing dropout to prevent overfitting
L1 and L2 regularization methods are often used in neural networks to prevent model overfitting. Some parameters in the loss function are limited by this method (Jiang et al., 2018;Huang et al., 2020). Neurons and fully connected layers are randomly deleted through Dropout technology to form new samples (Srivastava et al., 2014). The overfitting problem is effectively solved by this sample replacement effect. At the same time, the entire network is no longer very sensitive to the specific weights of neurons. Therefore, the generalization ability of the model has been greatly improved Chang et al., 2022).

Particle swarm Optimization(PSO)-LSTM process
The steps of PSO-LSTM are as follows, and the flow chart is shown in Figure 3.
1) The dataset is divided into LSTM training set, validation set and test set. 2) Initialize the PSO algorithm: parameters such as the number of iterations, inertia weight, learning factor and iteration step range are initialized. Then, the position and velocity of each particle are randomly initialized.   Frontiers in Earth Science frontiersin.org 07 Teng et al. 10.3389/feart.2022.1069182

FIGURE 5
Stress-strain curves of frozen silty clay with different freeze-thaw cycles.
Frontiers in Earth Science frontiersin.org 08 3) The three-dimensional particle swarm composed of Neurons, dropout and batch_size is randomly generated. 4) In each iteration, the LSTM training set is used as the training set, and the PSO optimization set is used as the test set. Therefore, the LSTM prediction process is simulated and the solution (particle) with the smallest prediction error in iterations is selected. 5) When the end condition is reached, the historical optimal solution of the swarm is output. The optimal number of neurons, dropout probability and batch_size are input to the next LSTM prediction. Then MSE is used to evaluate the fitting effect of stress-strain curve of frozen silty clay.,Finally, the optimized prediction results are obtained.

Optimization process analysis
According to the optimization parameters and constraints set above, the particle swarm optimization algorithm operation is performed. The fitness change curve is shown in Figure 4. It can be seen from Figure 4A that when the iteration exceeds 4 times, the fitness value quickly reaches the maximum. Therefore, the optimal solution is selected when the number of iterations (search steps) is 0-4, and the local optimal solution and the global optimal solution are output for each cycle. After 20 iterations, the global optimal solution can be obtained: [Number of neurons in the first layer, dropout, batch_size]= [61, 0.09, 95]. Therefore, PSO-LSTM structure optimization combination parameters are obtained and input.
It can be seen from the MSE attenuation curve of PSO-LSTM model in Figure 4B above that when the sample reaches 7.5 epochs, the MSE value tends to be stable. It shows that the validation set of the model has reached a good accuracy at this time. After determining the structural optimal parameters of PSO-LSTM, the weights and thresholds of each model are iteratively trained using the training set to obtain the final model. The test set is used to predict and compare each model (Chang et al., 2020;Huang et al., 2022).

Test results analysis
75 groups of triaxial tests are conducted on frozen silty clay with initial moisture content(IMC) of 25%, 22.5%, 20%, 17.5%, 15%, freeze-thaw cycles (FTC) of 0, 3, 6, 9, 18 times and confining pressure(CP) of 0.1Mpa, 0.2Mpa, 0.3Mpa. It can be seen from Figure 5 that freeze-thaw cycles and confining pressure do not change the strain characteristics of frozen silty clay. With the increase of initial moisture content, there is hysteresis phenomenon in the failure strain of frozen silty clay.
It can be seen from Figure 5 that when the initial moisture content is 22.5% and 25%, the deviatoric stress shows three development stages with the increase of the axial strain: linear elastic growth stage, elastic-plastic development stage, and strain hardening stage. In the linear elastic growth stage, the deviatoric  stress is dominated by cementation elements composed of pore ice and soil particles. Due to the strain is small, the pore ice does not crush and melt, and the structure does not have cracks. At this time, the overall deformation modulus is large, and the deviatoric stress increases linearly.
In the elastic-plastic stage, due to the large strain at this time, micro-cracks are generated in the pore ice. Accompanied by the pressure-melting phenomenon at the ice crack, part of the cementation element structure is destroyed and transformed into friction elements. Since the

Frontiers in Earth Science
frontiersin.org 11 deformation modulus of friction element is smaller than that of cementation element, the overall deformation modulus decreases. The deviatoric stress increases slowly with the increase of axial strain. During the strain hardening stage, the number of cementation elements remains stable, and the friction element, as the main bearing structure, show linear growth trend with the increase of strain.
It can be seen from Figure 5 that when the initial moisture content is 15%, 17.5%, and 20%, the frozen silty clay shows strain softening properties. The pore ice is broken and local shear bands are formed. The number of cementing elements decreases rapidly, the number of friction elements increases, and the integrity of the structural block is destroyed. The soil samples shows strain softening state. And with the decrease of initial water content, the strain softening phenomenon is more obvious. And under the same working condition, different confining pressures do not change the evolution law of deviatoric stress.

Binary medium constitutive model of frozen silty clay
The simplified frozen soil binary medium model based on damage parameters is shown in formula (8). The symbolic interpretation and derivation process of formula (8) are given (Zhang et al., 2022).
Since the research object in this paper is frozen silty clay, the local strain coefficient c is obtained from the stress-strain curve. The comparison between the fitted value of the formula and the actual test value is shown in Figure 6.
It can be seen from Figure 6 that the theoretical values of the constitutive model of frozen soil binary medium are very close to the measured values under different parameters. In addition, there is a high matching degree in the linear elastic stage, elasticplastic stage and strain hardening stage of the stress strain curve. Therefore, it is further verified that the model can better simulate the curve relationship between stress and axial strain of silty clay in the low-temperature triaxial test. The process parameters of the binary medium model fitted by the frozen silty clay test are shown in Table 4.

Optimal model selection
In order to compare the calculation accuracy of the model and select the optimal model, the frozen silty clay binary medium model, random forest model, PSO-LSTM model proposed in this paper and the frozen soil triaxial test are compared and analyzed. The stress sharing curve of the cementation element and the friction element is made by Formula (1), and the specific comparison is shown in Figure7. The deviatoric stress-strain curve under IMC25-FTC18-CP0.2 is shown in Figure 7. It can be seen from Figure 7 and formula (1) that the deviatoric stress curve of frozen silty clay is formed by superposition of the stress sharing curve of cementation element and the stress sharing curve of friction element.
In the linear elastic stage, the evolution trend of deviatoric stress of frozen silty clay is more consistent with the sharing curve of cementation element. At this time, the deviatoric stress shared by the friction elements is less than 0.039 MPa. Therefore, the pore ice inside the cementation element is not crushed and melted in the linear elastic stage. The integrity of cementation element is good, and the cementation element mainly plays the role of bearing force under low temperature triaxial action.
In the elastic-plastic stage, the stress sharing curve of friction element shows a linear growth trend, and the stress sharing curve of cementation element begins to grow slowly. When the strain exceeds 2.7%, the difference between the stress sharing curve of cementation element and the total stress curve gradually increases. In the strain hardening stage, the stress sharing curve of cementation element tends to be stable. It shows that cementation element has yielded and stabilized at this time, and the stress sharing curve of the friction element still maintains linear growth trend. It shows that pore ice in the cemented element is crushed and melted under the action of low

FIGURE 11
Three dimensional changes of freeze-thaw cycles, initial moisture content and strength under different confining pressures.

FIGURE 10
Three dimensional changes of confining pressure, initial moisture content and strength under different freeze-thaw conditions.  Figure 7 that in the range of axial strain of 20%, the stress sharing curve of cementation element is always larger than stress sharing curve of the friction element. The peak value of the total stress is also the maximum value of the superposition of the stress sharing curve of the cementation element and the friction element. It can be seen from Figure 7 that the calculated values of the frozen silty clay binary medium model and PSO-LSTM model are completely close to the test values. The constitutive relation of frozen silty clay can be better simulated by PSO-LSTM model. Because excellent matching degree can be achieved in linear elastic stage, elastic-plastic stage and strain hardening stage. The calculated values of the RF model and the test values show a high consistency in the strain hardening stage, but there is a large difference between the linear elastic stage and the elastic-plastic stage, which indicates that the model has a large error and still needs to be adjusted.

Model evaluation
Based on the random forest model test set, the stress-strain relationship of frozen silty clay with different moisture content after the 18th freeze-thaw cycles is obtained In the constitutive relation curve obtained by selecting the random forest model, the working condition with lower R 2 is only shown in Figure 8. Therefore, the worst fitting condition of random forest model is selected as the test set of PSO-LSTM model. Finally, the universality and applicability of PSO-LSTM model can be well verified.
It can be seen from the results in Figure 8 that the strains corresponding to the maximum deviatoric stress errors of the RF model are 18.18%, 7.92%, 4.79%, and 15.04% under the initial moisture contents of 22.5, 20, 17.5, and 15%. The deviatoric stresses corresponding to the RF model are 6.16, 6.50, 5.76, and 2.07 MPa respectively, and the deviatoric stresses corresponding to the test are 5.80, 6.10, 5.45, and 2.41MPa respectively. The errors are 5.8%, 6.1%, 5.38%, and 16.4% respectively. Therefore, the trend of the random forest model curve is consistent with that of the test curve, but there is always a large deviation in the whole evolution process of the partial stress. The evaluation indexes of the initial training set and test set of the RF model are significantly different, and the training set is more excellent. It shows that RF model has over-fitting phenomenon in the training process, and the unsatisfactory test results are generated. Therefore, cross-validation is added to avoid overfitting. In the end, the maximum error of the model is still 16.4%.
The maximum error of PSO-LSTM model is only 1.56% under IMC20FTC18CP0.3 working condition. Compared with RF, the predicted value of the PSO-LSTM model constructed in

FIGURE 13
Influencing factors of soil particle and pore ice force.

FIGURE 12
Three dimensional changes of freeze-thaw cycles, confining pressure and strength under different initial moisture content.
Frontiers in Earth Science frontiersin.org this paper fits better with the actual measured value, and the relative error is smaller and more stable. As can be seen from Figure 8, the evolution of deviatoric stress under each condition can be well simulated. The constitutive relationship predicted by the random forest model will have fluctuation segments, and it is prone to the problem of low strain sensitivity. The PSO-LSTM model can well simulate the evolution trend of deviatoric stress in each working condition.
Analysis of the data in Table 5 shows that the R 2 of the deviatoric stress prediction data of the frozen silty clay obtained by the RF model, the binary medium model and the PSO-LSTM model are all above 94%. Therefore, in the process of iterative training, the internal test laws of initial moisture content, freezethaw degradation effect and confining pressure are well excavated, and more reasonable constitutive prediction values are output. However, compared with the RF model and frozen soil binary medium model, the R 2 obtained by the PSO-LSTM model reaches 98.44%. At the same time, MSE and MAE are the lowest, which are 0.0078 and 0.01778 respectively, indicating that the model prediction effect is the best. The evaluation index of the initial moisture content of 25 is only shown in Table 5. Therefore, the evaluation index of the model test set under other conditions is given in the following Figure 9.
It can be seen from Table 6 that the accuracy of the optimal model under each working condition is above 90%. It can be seen from the comparison model that the accuracy of the PSO-LSTM model in each working condition is higher than that of the RF model. At the same time, RMSE and MAE errors are lower than RF model in all conditions, which is enough to show that the optimal model can well simulate the constitutive relation of frozen silty clay in all conditions.
The evaluation indicators of the RF and PSO-LSTM models on the validation set and the test set are relatively consistent, so the evaluation indicators of the validation set are omitted here. At the same time, the RF model and the PSO-LSTM model have not been fitted. In addition, the R 2 of PSO-LSTM for frozen silty clay constitutive prediction in the test set are all greater than 0.99. At the same time, RMSE, MSE and MAE are the smallest under each water content working condition, and they are the best in all networks. The better fitting accuracy and prediction performance of this model have been fully proved, and it is very accurate for the constitutive prediction of frozen silty clay. It also provides a new method for the prediction of frozen silty clay constitutive model.

FIGURE 15
Distribution of feature importance.

FIGURE 14
Structural model of frozen silty clay.
Frontiers in Earth Science frontiersin.org 5 Discussion 5.1 Effect of freeze-thaw cycles, confining pressure and initial moisture content on strength of frozen silty clay As shown in Figure 10, the strength tends to be stable after 6, 9, 6, 9, and 9 freeze-thaw cycles at initial moisture content of 25, 22.5, 20, 17.5, and 15%, respectively. After 18 freeze-thaw cycles, the shear strength of frozen silty clay decreased by 2.66%, 11.85%, 18.83%, 16.79, and 29.02%, respectively. Under the same confining pressure, the shear strength of frozen silty clay decreases first and then tends to be stable with the increase of freeze-thaw cycles. With the increase of freeze-thaw cycles, the bearing structure of frozen silty clay changes from initial coarse particles and pore ice to more deteriorated fine particles and pore ice. That is, the positions of the pore ice in frozen silty clay are occupied by more fine particles. On the one hand, the thickness of pore ice is reduced, and thin pore ice is more likely to expand cracks under low temperature triaxial stress. On the other hand, the unfrozen water film on the particle surface shows lubrication effect. And with the increase of freeze-thaw cycles, the rolling friction and sliding friction between particles increase more, which makes the shear strength decrease more obviously. Under a certain confining pressure, the rolling friction and sliding friction between particles tend to be stable after 6 and 9 freeze-thaw cycles. So the connection, arrangement and stress history of the soil particles are changed by the freeze-thaw cycles. This deterioration effect decreases with the increase of freezethaw cycles and tends to be stable after a certain number of freeze-thaw cycles (6 or 9 times). Figure 11, the deviatoric stress of frozen silty clay increases with the increase of confining pressure. This is very consistent with the actual project. And with the increase of the initial moisture content, the effect of confining pressure on the strength increase is no longer obvious. This is because in the third stage, the stress sharing curve of the friction element shows a linear growth trend. The stress sharing of the cemented element keeps the yield stable, and with the increase of the confining pressure, the strength of the cemented element for the yield stability is also greater. When the optimum moisture content is reached, the pore ice of cementation element tends to be saturated. At this time, with the increase of the initial moisture content, the expansion of the cracks generated by the thin pore ice is not obvious. Therefore, the influence of the confining pressure on the yield stability of the cemented element is gradually reduced.

As shown in
As shown in Figure 12, the strength of frozen silty clay is the highest near the initial moisture content of 20%. With the increase of freeze-thaw times and the decrease of confining pressure, the strength of frozen silty clay decreases. Therefore, the most unfavorable working condition is IMC15-FTC18-CP0.1, and the highest strength working condition is IMC20-FTC18-CP0.3.
Compared with Figures 10, 11, it is obvious that the strength difference of frozen silty clay is the largest under different initial moisture content. That is, the initial moisture content has greater influence on the strength of frozen silty clay (Cao et al., 2022;. This is because the effects of freezing and thawing and confining pressure on the connection and arrangement of soil particles and the stress history are not as significant as the initial moisture content. The bearing structures of soils with different initial moisture contents are quite different. As a result, the difference of pore ice thickness is large, and the difference of crack germination and expansion mode is increased. On the other hand, the difference of rolling friction and sliding friction between frozen silty clay particles under different initial moisture content is large, resulting in a large difference in shear strength.

Frozen soil structure model
The structural composition of unsaturated frozen soil  which influences factors of soil particle and pore ice force and soil strength change process is shown in Figure 13 below.
The structural model of frozen silty clay can be seen in Figure 14. According to the concept of the binary medium model (Zhang et al., 2022), cementation element is a stable ice-soil skeleton composed of clay particles and pore ice. The deformation modulus of cementation element is large, and it usually shows better integrity in case of small deformation. Under triaxial stress or low temperature, the pore ice inside There are a large number of cemented elements and friction elements in the frozen silty clay. And under the action of triaxial stress and low temperature, cementation element is gradually transformed into friction element. The strength characteristic difference of frozen soil is closely related to the quantity difference between cementation elements and friction elements.
The cementation element plays a major role in linear elastic stage, and the stress sharing curve of friction element can be ignored. In the elastic-plastic stage, the bearing capacity of the cementation element reaches the yield stability and continuously transforms into the friction element. During the strain hardening stage, the stress sharing curve of the cementation element remains stable, and the stress sharing curve of the friction element continues to increase linearly.

Contribution weight analysis of stressstrain curve of frozen silty clay
Compared with confining pressure and freeze-thaw cycles, the initial moisture content has the greatest influence on the evolution trend of deviatoric stress. Therefore, machine learning method is used to study the importance of features. It provides a new method for the study of the contribution weight of different factors to frozen silty clay constitutive curve. The contribution weight of deviatoric stress is shown in Figure 15.
The contribution weight of the four factors to the deviatoric stress can be seen from the importance map of random forest characteristics. According to the above experimental study, the trend of the stress-strain relationship curve will be changed by the different initial moisture content. The state of frozen silty clay before the failure will be determined by the initial moisture content, and the position of peak strain will also be changed. With the increase of initial moisture content, the failure stage changes from strain softening to strain stabilization, and then to strain hardening. At the same time, the failure stress will also lag. This also verifies the strong connection between the moisture water content and the stress evolution law.
In a certain range (IMC = 25,22.5,20,17.5, 15%, FTC = 0, 3, 6, 9, 18 times, CP = 0.1, 0.2, 0.3 MPa), the influence factors of frozen silty clay stress from strong to weak are: initial moisture content > strain > confining pressure > freeze-thaw cycles. The contribution of confining pressure to the stress evolution law is small, which is because the freezing depth of the seasonal frozen area is less than 2 m. Therefore, the gradient and range of the confining pressure setting in this paper are small, so the contribution weight is small.

Comparison and promotion of frozen silty clay binary constitutive model, RF model and PSO-LSTM model
Although the binary constitutive model of frozen silty clay has obtained a good prediction effect. And the reference formula is provided, which is beneficial to the application in practical engineering. However, there are some shortcomings in the model. For example, the binary constitutive model of frozen silty clay in this paper is only suitable for the fitting of the constitutive relationship under the condition of IMC=25%, but it can not achieve a good fitting effect for the IMC=22.5%, 20%, 17.5%, and 15%. Therefore, the formula fitting model still needs to find the constitutive formula of similar curve law. The formula fitting software 1STOPT can also be used to perform a random search of constitutive relations. But some parameters in the formula cannot be explained mathematically and physically. Therefore, in-depth research and promotion are limited.
Based on the shortcoming that the above model is applicable to a single working condition, the RF model is established in this paper. The RF model can be applied to all initial moisture content conditions, but the model prediction accuracy is not satisfactory.
The PSO-LSTM model is established in this paper to solve the problems of the above models. There are many advantages, such as high prediction accuracy, wide application range, and no need to propose based on certain assumptions. The LSTM optimal combination input parameters are searched by PSO, and the parameter adjustment speed of the model for the data learning process of other soil types is greatly increased, which is conducive to popularization. Therefore, it is only necessary to carry out frozen soil triaxial tests of different soil types in various regions. PSO-LSTM model has extremely high accuracy in constitutive prediction, which can provide reference for frozen soil engineering application in cold regions.
Moreover, the deep learning algorithm has a strong fitting ability for multiple operating conditions and multiple factors. The data driven deep learning method will have an excellent fitting effect for more dimensional sample data. For other soil types, the triaxial test data are imported, and the high-precision constitutive relation curve can be obtained by modifying the training set.
In the Code for Design of Soil and Foundation of Building in Frozen Soil Region, the characteristic value fa of foundation bearing capacity is selected based on temperature and soil type. Only four types of frozen soil with ground bearing capacity characteristic values at -3°C~-0.5°C are given in the specification. However, it is not suitable for saline frozen soil and frozen peat soil, and a large number of experiments need to be carried out to revise the specification. Therefore, the optimal model proposed in this paper has important generalization significance.
That is, a more comprehensive reference  Figure 16, the predicted strength law is consistent with the actual situation. It provides a new method for the constitutive study of frozen soil engineering in cold regions, and provides a favorable reference for the strength of remolded soil in seasonal frozen regions. Therefore, it is of great significance to carry out the model study in this paper. Random forest and PSO-LSTM can provide a new method for frozen soil constitutive research. The constitutive model of frozen soil is developed based on mathematical skills, which can deeply learn the data of frozen soil triaxial test. Machine learning model has strong nonlinear mapping ability, and is not limited by frozen soil constitutive formula In the process of deep mining of multi-dimensional large sample data, there is usually a better fitting effect. Therefore, machine learning shows extremely high accuracy and wide applicability in frozen soil constitutive model prediction. It can be used to modify the characteristic value fa of foundation bearing capacity of the Code for Design of Soil and Foundation of Building in Frozen Soil Region.

Conclusion
The study of constitutive model is of great significance to the construction of frozen soil engineering. In this paper, 75 sets of experiments are carried out, and a comparative study of the constitutive models of frozen silty clay is carried out. The following conclusions are drawn: 1) Under the same confining pressure, the shear strength of frozen silty clay decreases first and then tends to be stable with the increase of freeze-thaw cycles. The strength tends to be stable after 6, 9, 6, 9, and 9 times freeze-thaw cycles at initial moisture content of 25, 22.5, 20, 17.5, and 15%, respectively. After 18 freeze-thaw cycles, the shear strength of frozen silty clay decreased by 2.66%, 11.85%, 18.83%, 16.79, and 29.02%, respectively. 2) Compared with RF model and frozen soil binary medium model, the R 2 obtained by the PSO-LSTM model reaches 98.44%. At the same time, MSE and MAE are the lowest, which are 0.0078 and 0.01778 respectively, indicating that the model prediction effect is the best. In the process of iterative training, the internal test laws of initial moisture content, freeze-thaw degradation effect and confining pressure are well excavated, and more reasonable constitutive prediction values are output. 3) In a certain range (IMC = 25,22.5,20,17.5,15%,FTC = 0,3,6,9,18 times, CP = 0.1, 0.2, 0.3 MPa), the influence factors of frozen silty clay stress from strong to weak are: initial moisture content > strain > confining pressure > freeze-thaw cycles. Tongming et al., 2021, Yuanming et al., 2010.

Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions
Y-CT: Design, Experiment, Algorithm, Data Processing, Writing; Z-CT: writing, translation; J-LL: Experiment, Writing; Y-DZ: Experiment, Writing, X-YL: Writing, Z-WL: Experiment, C-YT: Experiment. All authors contributed to the article and approved the submitted version.