Establishment and Experimental Verification of Stress- Temperature Coupled Damage Model of Warm Frozen Soil

Under the background of rising environmental temperature, the state of warm frozen soil near the phase transition is extremely unstable. In order to explore the relationship between temperature and mechanical properties of warm frozen soil, a damage model of warm frozen soil structure under the coupling of stress and temperature is established based on the strain equivalent theory of damage mechanics. Based on the Mohr-Coulomb criterion, the nominal stress is used to represent the stress damage of frozen soil elements, the initial elastic modulus is used to represent the temperature damage, and a composite damage factor is introduced to describe their coupled relationship. Through a triaxial compression experiment of frozen soil, the experimental data and stress-strain curve are obtained. The full-fitting method based on the experimental data (method 1) and the semitheoretical semifitting method based on the characteristic points of the stress-strain curve (method 2) are used to obtain the shape parameters and scale parameters of the stress-temperature coupled damage model corresponding to different fitting methods. Based on the triaxial compression tests of frozen sand and frozen silty clay, the reliability of the stress-temperature coupled damage model results obtained by the two parameter determination methods under the conditions of strain softening and strain hardening is verified. The results show that both methods are applicable under the condition of strain softening and strain hardening and method 2 is better than method 1 under the condition of strain softening. Compared with the prediction results of the single stress damage model, the stress-temperature coupled damage model can effectively reduce the influence of the parameter estimation error on the results and improve the overall stability of the model.


Introduction
The alpine and polar cryosphere are retreating and are expected to continue melting in the future [1]. With global warming and increasing ground temperatures, permafrost gradually degenerates, and the extent of warm frozen soil further increases [2,3]. There is a close dynamic equilibrium relationship between unfrozen water content and temperature. As a cemented phase, ice exerts a crucial influence on the mechanical properties of frozen soil. Unfrozen water is highly sensitive to temperature changes in the near-phase transition region of warm frozen soil. Small temperature differences directly lead to changes in the unfrozen water membrane thicknesses of frozen soil particles. The difference is reflected in the mechanical index of frozen soil elements and is expected to increase to a certain extent [4][5][6]. Compared with normal frozen soil, warm frozen soil is likened to a mixed state of thawed soil and frozen soil [7]. This type of material strongly intensifies the significance of internal permafrost defects compared with normal permafrost in terms of both the relative quantity and influence on the soil mechanical properties [8][9][10]. The elastic modulus of frozen soil is assumed to be constant in numerical calculations owing to its high instantaneous strength [11][12][13][14]. However, for warm frozen soil, the deterioration of the frozen soil stiffness is more apparent because of the increased unfrozen water content at higher temperatures [15]. Large calculation uncertainties are introduced if mechanical parameters under a certain fixed temperature are used to numerically analyze the stress-strain relationship of frozen soil under variable temperature conditions, which can affect the safety of regional facilities [16].
Mechanical studies of warm frozen soil have mainly focused on strength and deformation characteristics based on a constitutive model; however, the overarching study of constitutive models is less constrained [7]. Lai et al. [8] showed that the strength and elastic modulus distribution of warm frozen soil is discrete. The distribution of the strength and elastic modulus of frozen soil can be described by a Weibull distribution. Strain is taken as the damage degree, and statistical damage mechanics have been used to develop the earliest damage model of warm frozen soil. According to Li et al. [9], the concept of mechanical damage represented by strain is not sufficiently clear and lacks theoretical support. Therefore, the Mohr-Coulomb criterion was introduced to assess the damage state of warm frozen soil elements, and a 1-D parameter expression was derived and achieved good results. The constitutive study of warm frozen soil has consistently focused on deformation development under the action of a single stress factor, and temperature has been considered to be a fixed condition rather than an action factor, which limits the application scope of the constitutive model ( [17]; Zhang, D et al., 2018). Previous studies on the coupled effects of stress and temperature have mainly focused on the level of physical field control relations, whereas investigations of a deeper constitutive model have occurred less frequently [14,18]. According to the strain equivalent theory of damage mechanics, Zhang et al. (2010) obtained a composite damage model of freeze-thaw rock under temperature and stress action. Huang et al. [19] deduced the statistical damage model of rock freezingthawing and stress coupling in cold regions and verified it by a tunnel project in cold regions. Yu et al. [20] introduced the H-B (Hoek-Brown) criterion to extend the application scope of the model to describe the stress-strain relationship of rock under warm and confining pressure conditions.
In view of the above problems and shortcomings, this study introduces the composite damage theory into a stress-strain simulation prediction of warm frozen soil, establishes a stress-temperature coupled damage model characterized by the Mohr-Coulomb criterion for 3-D (three-dimension) stress damage of warm frozen soil elements [21], verifies the model prediction effects with relevant experiments, and discusses its effectiveness and stability.

Stress-Temperature Coupled Damage Model
The initial internal state of warm frozen soil is damaged under the combined action of stress and temperature. Some existing damage defects will further develop, and new damage defects will gradually form [8,9].
According to the damage theory of a continuous medium, the development of this type of damage can reduce the actual stress element in the body under external action, which can be understood in terms of material structure damage [22]. This damage is reflected in the physical parameters of the material, which are generally reduced compared with the original parameters and can be described by a dam-age factor. The nominal strain due to material degradation under nominal stress is described as follows: where σ is the nominal stress, ε is the nominal strain,Ẽ is the effective elastic modulus, E is the initial state of the elastic modulus, and D is the damage factor. A coupled damage model based on Equation (1) where m and a are the shape and scale parameters of the distribution, respectively. The Mohr-Coulomb theoretical analysis is used to determine the extent of damage characterization F [9]. In triaxial compression tests under low confining pressure, the Mohr-Coulomb shear strength formulation can be expressed as follows by axial stress σ 1 and confining pressure σ 3 [24].
where c is the cohesive force and φ is the friction angle. Upon transferring terms σ 1 and σ 3 , the following equation is obtained: whereσ 1 andσ 3 represent the real stress of the material with the development of warm frozen soil damage. From the generalization of Equation (1), the following equation is obtained: Thus, the damage characterization F expressed by nominal stress can be established as follows: According to the stress-strain relationship of materials 2 Geofluids under 3-D stress, it can be concluded that In triaxial compression experiments, σ 2 = σ 3 , and the expression of 1 − D m can be obtained by the following expression: Equation (8) can be substituted into Equation (6) to obtain the following equation: When F corresponding to different axial strains ε 1 is calculated from Equation (9), σ 1 is an unknown quantity to be solved.
The extent of the damage evolution of frozen soil elements is assumed to be small over a small range of strains. The increment of F over a small range can be ignored and therefore expressed as follows: where n represents the accumulated step number of strain ε when accumulated from ε = 0 over a certain step. When n = 0, this represents the stress value corresponding to the initial state before the axial load begins after the end of the confining pressure load. At this time, σ 1 = σ 3 , and Equation (10) can be converted into the following expression: 2.2. Temperature Damage. The cementation force of ice on soil particles and its own strength in frozen soil provide rigidity values higher than those in the melting state [10]. The ice strength and cementation ability decrease with increasing temperature. This type of composite material degradation due to temperature effects can also be reflected by the concept of damage.
Equation (1) shows that the effects of material damage on mechanical parameters are most directly reflected in a reduction in the elastic modulus. Therefore, the temperature damage factor D T is defined, which represents the degree of deterioration from the initial state caused by an increase in temperature and a decrease in ice strength and cementation force in frozen soil.
where E T is the elastic modulus at a certain temperature after heating from the initial state. The elastic modulus of frozen soil can be expressed as a power function with temperature. Previous experimental studies have shown that an empirical relationship can be expressed by the following equation [13]: where β and γ are experimentally determined parameters.

Coupled Effect of Stress and Temperature. Zhang and
Yang's research shows that the composite damage caused by temperature and stress can be deduced and characterized by the composite damage factor D c [25].
The strain caused by the effective stress corresponding to the stress action acting on the stress element after the temperature damage is the same as the strain caused by the effective stress corresponding to the temperature action acting on the stress element after the stress damage. Further derivation shows that the real stress of the material is equal to the product of the strain after stress damage reduction and the elastic modulus after temperature damage reduction. Its essence is to treat the material after one damage as a nondestructive state and then calculate the influence of another damage on the state.
Equations (2) and (12) can be substituted into Equation (14) to obtain the following expression: Substituting Equation (13) into Equation (15) can be simplified to obtain the following expression.
where T 0 is the initial temperature and λ is defined as the temperature deterioration coefficient.
Then, Equation (16) can be introduced into Equation

Geofluids
(1), which gives the stress-temperature coupled damage model of warm frozen soil as follows:

Model Parameters
The stress damage of warm frozen soil is in accordance with a Weibull distribution [26]. The shape parameter m and scale parameter a in this distribution are two important parameters in the stress-temperature coupling model and are directly affected by the stress-strain curve. The curves of deviator stress σ 1 − σ 3 and axial strain ε 1 obtained from triaxial tests of warm frozen soil take different forms owing to different confining pressures and soil properties [16,27].
In general, the stress-strain curves of frozen sand, frozen silty sand, and frozen silty clay with high moisture content show typical strain softening, and the deviator stress shows a peak value. After peak strain, the deviator stress decreases with increasing strain ( [28][29][30]. No peak value is observed in the stress-strain curves of fine-grained soil with low water content, and the deviator stress gradually increases with increasing axial strain, but the growth rate gradually decreases and tends toward a stable value [31]. For two different stress-strain curves, there are two typical two-parameter solutions of a Weibull distribution. One solution is to obtain m and a by linear fitting (method 1), and the other solution is to further deduce Equation (9) through characteristic data points in the stress-strain curve and obtain an expression for m and a expressed by the values of these characteristic points (method 2) [32,33].
Using two logarithm operations at both ends of the equation yields, The left-hand terms and right-hand F can be obtained from triaxial experiments, which gives Then, the linear relationship between the following three can be obtained.
A set of X and Y values corresponding to the experimental data under a certain confining pressure and temperature can be obtained by introducing the triaxial experimental values. The expression of the final m and a values can be obtained by linear fitting and further calculation.
The detailed calculation method of relevant parameters is further explained in the experimental verification section. (18) is first modified to obtain

Method 2. Equation
The relationship between ε 3 , axial stress σ 1 , and confining pressure σ 3 can be established according to the generalized Hooke's law.
The following expression is obtained if the above equation is modified and substituted into Equation (2).
By taking the peak stress of the strain-softening curve as the characteristic point, the following relationship can be obtained: Determining the total differential of Equations (27) and (29) yields The selected damage characterization quantity is determined based on the Mohr-Coulomb theory. When Equation (27) is substituted into Equation (4), F can be expressed by ε 1 and σ 3 . For the convenience of distinction, F in this case is marked as F 1 .
In the application of the frozen soil damage model, m and a are often regarded to only be related to confining pressure σ 3 [9,24], and the following relationship can be obtained: If Equations (36), (38), and (39) are substituted into Equation (32), the total differential relationships of σ 1 , σ 3 , and ε 1 can be obtained as follows: Among them, If Equations (37), (38), and (39) are introduced into Equation (33), the total differential relations of σ 1 , σ 3 , and ε 3 can be obtained as follows: Among them, The total differential relationships among σ 1 , ε 1 , and ε 3 can be obtained by combining Equations (40) and (42).
Among them, The numerical conditions of the characteristic points of the stress-strain curve are introduced to solve the problem. Equation (30) is first introduced into Equation (27) Among them, Then, Equation (31) is introduced into Equation (45) and simplified to obtain After each value is introduced, it can be simplified to obtain The expression of m and a can be obtained by combining Equations (47) and (50) as follows: The values in the above two equations can be obtained from the deviator stress vs. axial strain curves obtained from the triaxial compression experiments. If the confining pressure σ 3 = 0, Equations (51) and (52) are reduced to a 1-D 5 Geofluids uniaxial compression.
A comparison shows a similar result of uniaxial compression to that derived in previous studies [9]. Specific values of the shape and scale parameters can be determined by the above derivation.

Experimental Method and Parameter Determination
According to the above information, the frozen triaxial compression test is used to determine the parameters of the constitutive model and verify the model. At the same time, in order to verify the applicability of the model, the strain softening and strain hardening curves are verified, respectively. Frozen sand (strain softening type) and frozen silty clay (strain hardening type) are selected for the experimental soil samples.  Table 1.
In order to ensure that the internal and external temperature of the specimens is constant and the internal moisture distribution is uniform, the low-temperature accelerated freezing is adopted to reduce the moisture migration, and the confining pressure chamber and the specimens are cooled together for a long time to reduce the test temperature fluctuation. The specimens had a diameter of 39.1 mm, height of 80.0 mm, dry weight of 15.3 kN/m 3 , and relative compactness of 77.84% after compaction preparation and were saturated and frozen for 48 h in a lowtemperature tank at −25°C. Three sets of triaxial compression experiments were carried out at a loading rate of 1 mm/min and at constant temperatures of −0.3, −1, or −2°C for 24 h in a target-temperature chamber before testing. The confining pressures of each group were 0.05, 0.15, or 0.25 MPa. The parameter experiment and verification experiment are described as Figure 1.
Nine groups of triaxial compression experiments can be divided into two parts: parameter experiment and verification experiment. The parameter experiment consists of 5 groups. The stress damage factor is determined by 3 groups of triaxial compression tests with different confining pressures under the same temperature conditions, and the temperature damage factor is determined by 3 groups of triaxial compression tests with different temperatures under the same confining pressure conditions. The other 4 groups can be used as validation experiments.
The experimental results are shown in Figure 2.
It can be seen from Figure 2 that the peak stress increases with the decrease of temperature at the same confining pressure. When the temperature is the same, the peak stress increases with the increase of confining pressure. The ice cementation of warm frozen sand is the main source of strength, and with the increase of temperature, the ice strength decreases, resulting in the decrease of peak stress. At the same time, when the confining pressure increases, the connection between ice and sand particles is closer, and the peak stress increases.
The stress-strain curve shows the characteristics of strain softening, which is mainly due to the failure of the connection between ice and sand particles providing main strength after the partial stress reaches the peak stress.
The stress-strain curve is approximately linear during the initial strain stage (Figure 2). The elastic modulus of Equation (12) characterizes the initial undamaged state of the frozen soil, and values closer to the initial elastic modulus are associated with more precise model solutions. Therefore, the modulus of the linear elastic section in the range of 0%-0.5% strain is approximated as the elastic modulus E of the initial damage state.
The basic mechanical parameters calculated from the experimental results are listed in Table 2.
The results in Table 2 show that the expression of the temperature damage modulus E T in the stress-temperature coupled damage model can be obtained from Equation (13). The above analysis shows that in theory, only the confining pressure must be controlled, and β and γ in Equation The goodness of fit R 2 = 0:9998 reflects high accuracy.

Parameter
Determination of the Damage Model. As shown in Figure 1, the nondestructive temperature is -2.0°C [7]. Figure 2 shows the solutions of m and a of the stress damage model using the above two methods. According to the data in Figure 2, the shape parameters and scale parameters of the stress damage model are solved by the above two methods. The final calculated results are shown in Tables 3 and 4. Table 3 includes the full fitting of Equations (9) and (21) Table 5. The particle analysis curve is shown in Figure 3. A cylindrical specimen with a dry density of 1.75 g/cm 3 , diameter of 10 cm, and height of 20 cm was prepared. The moisture content of the specimen was 21% (saturated soil). Then, the specimen was frozen for 48 h in a -25°C low temperature chamber, and the experiment was conducted in a target temperature chamber for 24 h. Then, three groups of triaxial compression experiments were carried out at a loading rate of 1 mm/min at -0.5°C, -1°C, and -2°C. The confining pressures of each group were 0.05 MPa, 0.15 MPa, and 0.25 MPa. The parameter experiment and verification experiment are shown as Figure 4.
The experimental results are shown in Figure 5. It can be found from Figure 5 that, unlike frozen sand, frozen silty clay does not appear strain softening. With the increase of strain, the growth rate of deviatoric stress decreases, but maintains an increasing trend.
There is unfrozen water in warm frozen silty clay. The strength of the three-phase system formed by ice, water, and silty clay particles not only comes from the cementation of ice phase but also provides some strength by the interaction between particles. Therefore, it does not have brittle failure, but in a certain degree of deformation, the stress increases slowly with strain.
The elastic modulus of the linear elastic section in the range of 0-0.5% strain is taken as the elastic modulus of the initial damage state. The basic mechanical parameters obtained from Figure 5 are shown in Table 6.
Based on the above analysis, the stress damage of warm frozen soil can be determined by three groups of triaxial compression tests with different experimental confining pressures when the experimental temperature is -2°C. The temperature damage can be determined by three groups of triaxial compression tests with different temperatures under the same confining pressure.
The initial elastic modulus of three groups of different confining pressures corresponding to a confining pressure of 0.25 MPa is selected to determine the temperature damage.
The goodness of fit R 2 = 0:9712 can more accurately reflect the changes in data.

Parameter Determination of the Damage Model.
According to the stress-strain curve obtained from the triaxial test of frozen silty clay, the parameters of the damage model are determined. Taking -2°C as the nondestructive temperature and using the above two methods to determine the shape parameters and scale parameters of the stress damage model, the parameter values are obtained when the experimental confining pressures are 0.05 MPa, 0.15 MPa, and 0.25 MPa. In method 2, the peak stress and peak strain data are needed to determine the model parameters, while the stress-strain curve of frozen silty clay is the strain hardening type, and there is no significant stress peak. Therefore, when the strain is 15%, the corresponding deviatoric stress is taken as the peak stress, and the parameters are calculated accordingly. The final results are shown in Tables 7 and 8.
The above parameters can be brought into the constitutive model to predict and verify the experimental stressstrain curve to verify the effectiveness of the constitutive model.

Experimental Verification of the Strain Softening Curve.
According to the above parameters, the model was verified. As shown in Figure 1, four groups of confining pressures of 0.05 MPa and 0.15 MPa and experimental temperatures of -0.3°C and -1.0°C were used to verify the effectiveness of the model, and the results are shown in Figure 6. Figure 6 shows that the overall simulation results are in good agreement with the experiments. The results are most accurate when the confining pressure is 0.05 MPa, and a certain deviation is observed when the bias stress simulation is 0.15 MPa. By comparing the simulation accuracy of the four  7 Geofluids groups of confining pressure curves to the original experimental curves, an increase in confining pressure is associated with a gradually increasing degree of deviation between the calculated and experimental curves. Under the same confining pressure conditions, the model results more accurately reproduce the experimental simulations at higher temperatures.
The influence of the two methods on the simulation results is shown in Figure 6. The simulation results of the parameters obtained by method 2 are notably closer to the experimental curve than those obtained by method 1. The difference between the two methods is not substantial under low confining pressures; however, the calculated curve using method 1 deviates more from the experimental curve with increasing confining pressure than that of method 1.
The above verification process shows that the stresstemperature coupled damage model is more accurate in predicting the strain softening stress-strain curve when the stress damage under a nondestructive temperature state and the temperature damage corresponding to the target temperature are known.

Experimental Verification of the Strain Hardening Curve.
Through the above parameters, the stress-strain curves of four groups of nonparametric determination experiments       Figure 7. As shown in Figure 7, when the confining pressure is 0.05 MPa, the prediction of the stress-strain curve by method 1 and method 2 is more accurate. Method 1 is more accurate in predicting the phase of rapid stress growth than method 2. However, when the growth rate of stress gradually slows, the prediction results of the two methods have some deviation compared with the experimental curves, and the prediction of the change trend of the stress-strain curve of method 1 is slightly worse than that of method 2. When the confining pressure is 0.15 MPa, there is also a certain deviation between the predicted results and the experimental results in the slow stress growth section (the strain is greater than 9%), and the deviation is greater than that when the confining pressure is 0.05 MPa. Method 2 is usually used for strain softening curves. Since there is no peak stress or peak strain in the strain hardening curve, the parameters required for the solution of method 2 must be calculated by changing the strain to peak strain when the stress is approximately unchanged. The above information shows that the strain hardening curve of warm frozen soil conforms to the Weibull distribution, while the two methods only describe the stress-strain relationship from different perspectives. When the artificial peak strain and peak stress are used to solve the model parameters through method 2, the small increase in stress after the peak strain will be ignored, but this part will not affect the main part of the model application. Simultaneously, due to the extremely slow growth of stress, the error can be regarded as a negligible high-order small amount of predicted value.
Based on the above analysis, method 1 and method 2 are still more accurate. Method 1 is better than method 2 in the prediction of the stress increasing section, and method 2 is better than method 1 in the prediction of the stress increasing section. Method 2, which is generally applicable to the strain softening curve, can also effectively predict the strain hardening curve.
In conclusion, the constitutive model described in this paper can effectively predict strain softening and strain hardening materials.

Discussion
Confining pressure and temperature are two important factors that affect the mechanical properties of frozen soil, which must be addressed in practical applications of the damage model ( [34,35]. In most cases, the experimental conditions under which parameters are determined do not fully correspond to the actual conditions, so fitting predictions of model parameters are often necessary.
The stress-temperature coupled damage model only considers the effect of a single factor in the prediction of model parameters. The parameters of the stress damage factor only consider the change in confining pressure, while the parameters of the temperature damage factor only consider the change in temperature. The prediction formula of the two is the corresponding monomial formula of a single factor. However, when the stress damage model is used to predict the damage factors under a certain temperature and confining pressure, it is necessary to consider the effects of temperature and confining pressure at the same time because the damage factors are not distinguished. The prediction formula of the related parameters is the polynomial of the corresponding two factors.
The strain softening curve and strain hardening curve are analyzed, and the influence of the parameter error of the damage factor on the predicted stress-strain curve is discussed.
First, the stress-strain curve of frozen sand is used to analyze the strain hardening curve. In order to explore the influence of model parameter deviation on the prediction results, the influence of 10% increase or decrease of each parameter on the calculation curve of damage model   Figure 8.
When each parameter is increased or decreased by 10% from the original parameter, Poisson's ratio error has a negligible effect on the final calculation curve. The influence of the error of φ and a on the calculation results is mainly reflected in the peak value. With increasing φ and a, the peak strain and peak deviator stress increase. An increase in φ by 10% leads to a 2.41% increase in peak strain and a 2.55% increase in peak deviator stress. A reduction in φ by 10% leads to a 2.41% reduction in peak strain and a 2.71% reduction in peak deviator stress. When a is increased/decreased by 10%, the peak strain increases/decreases by 10.25% compared with the reference value, and the peak deviator stress increases/decreases by 10.01%. Therefore, the model results are more sensitive to changes in a.
The elastic modulus E and shape parameter m affect the stress growth, peak stress value, and stress decline stage of the calculation curve and show similar overall change trends. During the stress growth stage, larger parameter values are associated with faster stress growth, and the stress growth caused by a change in E is more significant than that caused by a change in m. With an increasing parameter value, the peak strain decreases and the peak deviator stress increases. A change in E more strongly affects the peak strain, whereas m exerts a stronger influence on the peak deviator stress. During the stress descent stage, smaller parameter values are associated with slower descent rates, and the calculated curves are more sensitive to a changing E.
Based on the above analysis, the four other parameters (except Poisson's ratio) affect the reconstruction curve. Under conditions of an identical amplitude change, the influence of a is the most apparent. When multiple parameters simultaneously change, the difference between the reconstruction curve and reference curve is not discussed      10 Geofluids because ν has little influence on the reconstruction curve.
The final results are shown in Figure 9. Figure 9 shows that when multiple parameters simultaneously change, the influence of each parameter exerts a composite superposition effect, and the three reconstruction stages of the stress-strain curve largely deviate from the original data curve. The variation in peak deviator stress is mainly due to the error of a, and the deviation of the stress growth stage is mainly due to the changes in E and m. A larger number of fitting parameters are associated with a greater influence of the model parameter error on the reconstructed curve accuracy.
Frozen silty clay is also compared with the above method. Based on the calculation parameters of the model  11 Geofluids when the confining pressure is 0.05 MPa and a temperature of -0.5°C, the prediction curve obtained after each parameter is increased or decreased by 10% as Figure 10.
The above figure shows that the overall trend of the influence of each parameter of frozen silty clay on the final prediction results is consistent with that of frozen sand, and the main difference is reflected in the shape parameters. With the increase or decrease in the parameters, the influence of the final curve is mainly reflected in the peak stress section and the stress drop section. There is an intersection point in this part. Before this point, with the increase in parameters, the stress predicted by the final model increases at the same time; the increase in parameters will reduce the stress predicted by the model. After this point, the change direction of the parameters and stress estimation are opposite. With the increase in parameters, the estimated stress decreases. In the strain softening curve, this point is located in the stress drop section, so its influence on the final result is reflected in the first half of the curve, and the stress prediction changes in the same direction as the parameters. For the strain hardening curve, the point moves forward, which is at the intersection of the stress increasing segment and the stress peak segment. This position creates the difference in the first half being hidden and the difference in the second half being enlarged.
Similar to the frozen sand, further comparison of the influence on the final results when the parameters change at the same time is shown in Figure 11. Figure 11(a) shows the final stress-strain curve when the other four parameters increase or decrease by 10% at the same time except for Poisson's ratio. From the results, it seems that the effect of multiparameter error on the final result is not significant. However, based on the above analysis, the change in the strain hardening curve and shape parameter is reversed, while the changes in other parameters are in the same direction, which forms a "counteracting" effect.
When the other three parameters are the same but the shape parameters are opposite, the final result is shown in Figure 11(b). The a curve is the result obtained when the parameters increase by 10% and the shape parameters decrease by 10%, while the b curve is the opposite. The figure shows that the error of each parameter is significantly enlarged after superposition. In the case of a single parameter error, the change in the shape parameter has a great influence on the final result. When the parameter increases  by 10%, the peak stress increases by 16.88%. However, the peak stress of the a curve increases by 33.23% compared with the original curve, which is approximately twice the difference in the final result caused by a single parameter error. The data in Tables 4 and 8 show that a small change in confining pressure and temperature leads to a relatively small range of m and a values; however, the final result is highly sensitive to a change in these two parameters [24]. When the binary polynomial fitting is used, the error of the parameter simulation values further increases compared with those of the single variable polynomial fitting. When the two parameters simultaneously have large errors, the model uncertainty sharply increases, which affects the reliability of the final result.
In summary, it is necessary to predict the elastic modulus E, shape parameter m, scale parameter a, Poisson's ratio ν, and friction angle φ under certain target temperature and confining pressure conditions (i.e., nonparametric experi-mental state) when simulating the stress-strain curve of warm frozen soil using the stress damage model. Among them, E, m, and a are simultaneously affected by changes in temperature and confining pressure, and a binary function is required to fit them. When using the coupled stresstemperature damage model, only m, a, ν, and E are estimated at the target temperature. In contrast to the single stress damage model, the parameters required in the stresstemperature coupled damage model can be fitted by a single element polynomial, which substantially reduces the parameter prediction uncertainty.
Through the reconstructed stress-strain curve after the parameters of strain hardening curve and strain softening curve change by 10%, it can be found that the shape parameters and scale parameters have the most significant influence on the final result. The scale parameter mainly acts on the peak point of the stress-strain curve, and the scale parameter changes in the same direction as the peak stress.    Figure 11: Comparison of the influence of a 10% increase or decrease in multiple parameters of the model on the calculation curve (frozen silty clay).

Geofluids
The shape parameter mainly acts on the peak stress, and the shape parameter changes inversely with the peak stress.
One can further compare the minimum amount of experiment required for parameter experiments of the two models. As shown in Figures 1 and 4, to predict any point under certain temperature and confining pressure conditions, a single stress damage model will require at least nine groups of experimental data to fit the parameters. However, if the stress-temperature coupled damage model is used, the temperature damage factor can be determined by only the experimental data at (i) three different confining pressures and the lowest temperature in the condition area and (ii) three different temperatures and one of the three confining pressures. Therefore, the temperature damage factor can be determined using only five sets of experimental data. At the same time, the error of parameter prediction can also be avoided by checking the experimental curve, and the stability of the model is higher.
The stress-temperature coupled damage model has obvious advantages over the single stress-strain model in terms of the prediction accuracy of the stress-strain curve and the experimental quantity of the parameter experiment.

Conclusions
When the frozen soil is in the near phase transition temperature range, the unfrozen water content increases sharply and the ice phase strength decreases, which belongs to the unique temperature damage mechanism of warm frozen soil, which is different from normal frozen soil. The temperature damage is incorporated into the constitutive model of warm frozen soil as a variable, and the stress temperature coupling damage model is established. It is an extension of the traditional frozen soil stress-damage model and introduces temperature variables to describe the degradation of warm frozen soil under the combined action of stress and temperature. Based on Hooke's law and the strain equivalent theory of damage mechanics, a warm frozen soil damage model under the coupling of 3-D axisymmetric stress (σ 2 = σ 3 ) and temperature (stress-temperature coupling damage model) is established, and the damage of warm frozen soil element is described by Weibull distribution. Through theoretical analysis, mathematical derivation, and experimental verification, the following conclusions are obtained.
(1) Composite damage theory is introduced into a stress-temperature coupled damage model of warm frozen soil. Based on the Mohr-Coulomb criterion, the nominal stress is used to represent the stress damage in the permafrost element, the initial elastic modulus is used to represent the temperature damage, and a composite damage factor is introduced to better describe the coupled relationship between them (2) The triaxial compression tests of frozen sand and frozen silty clay in the temperature range of 0~-2°C were carried out under 3 confining pressure condi-tions. The results show that the increase of temperature has a significant deterioration effect on the strength of warm frozen soil. With the increase of confining pressure, the peak stress of the material increases. When the confining pressure is the same, the peak stress decreases rapidly when the temperature changes from -2°C to 0°C (3) There are two general methods used to solve the shape parameter m and scale parameter a of the stress-temperature coupled damage model. The experimental results show that for strain-softening materials, the difference between the two methods is small when the confining pressure level is low. However, under high confining pressures, the reconstruction curve obtained by the semitheoretical semifitting method based on characteristic points of the stress-strain curve is closer to the experimental curve, and the method is also suitable for strain hardening materials (4) The results of triaxial compression experiments at constant temperature and confining pressure show that the stress-temperature coupled damage model offers good simulation and prediction abilities. The model's overall prediction accuracy is high in the stress development stage and decreases slightly with increasing confining pressure. At the same confining pressure, the model prediction accuracy increases with increasing temperature in the frozen soil (5) The stress-temperature coupled damage model requires fewer model parameters and only a 1-D function fitting to complete the prediction. This model effectively reduces the influence of system errors on the final results caused by parameter estimation, effectively reduces the workload of parameter experiments, and improves the model stability

Data Availability
All data are included in the manuscript.

Disclosure
The original version of this paper, Study on Stress-Temperature Coupled Damage Model of Warm Frozen Soil, has been published in ResearchGate as a preprint (10.1002/ essoar.10502573.1).

Conflicts of Interest
The authors declare that they have no conflicts of interest.