Stress–Strain Model for Freezing Silty Clay under Frost Heave Based on Modiﬁed Takashi’s Equation

Featured Application: Potential applications of this work include predicting the stable stress and strain ﬁelds and analyzing the inﬂuence factors of frost heave force and deformation of soil surfaces after freezing processes. Abstract: In analyzing frost heave, researchers often simplify the compressive modulus of freezing soil by considering it as a constant or only as a function of temperature. However, it is a critical parameter characterizing the stress–strain behavior of soil and a variable that is inﬂuenced by many other parameters. Hence, herein several one-dimensional freezing experiments are conducted on silty clay in an open system subjected to multistage freezing by considering the compressive modulus as a variable. First, freezing soil under multistage freezing is divided into several layers according to the frozen fringe theory. Then, the correlation between the freezing rate and temperature gradient within each freezing soil layer is investigated. Takashi’s equation for frost heave analysis is modiﬁed to extend its application conditions by replacing its freezing rate term with a temperature gradient term. A mechanical model for the stress–strain behavior of freezing soil under the action of frost heave is derived within the theoretical framework of nonlinear elasticity, in which a method for determining the compressive modulus of freezing soil with temperature gradient, overburden pressure, and cooling temperature variables is proposed. This study further enhances our understanding of the typical mechanical behavior of saturated freezing silty clay under frost heave action.


Introduction
Frost heave is a primary cause of the failure of critical infrastructure in cold regions, such as underground gas pipelines and roadways, including bridges-especially those with their footings in frost-susceptible soils [1,2]. Therefore, frost heave research on freezing soils is relevant in today's scenario [3,4]. Additionally, in cold regions where this is a major problem, engineers are attempting to find practical solutions.
The efforts to explain the mechanism of frost heave began in the 1920s [5]. Since then, numerous experimental and theoretical efforts have been undertaken to systematically investigate the frost heave characteristics of a variety of soils, including the mechanism, amount, and pressure [4,[6][7][8]. Various methods have been proposed to describe the frost heave process, such as the capillary model [9], rigid ice model [10][11][12][13][14], segregation potential model [15], thermomechanical models [16], and semi-empirical approaches [17,18]. In addition, recent research interests in frost heave and frost action in frost-susceptible soils have been renewed in China [19,20].
Indeed, the frost heaving of frost-susceptible soils is a coupled hydro-thermomechanical (THM) and thermomechanical (TM) interaction process. The development of mathematical modeling and computing technology has promoted the advancement of the TM and THM methods. Recently, based on the discrete ice lens theory, Fremond and Mikkola (1991) [16] established a TM model to explain the water heat transfer and mechanical characteristics of freezing soils. Konrad and Shen (1996) [21] used the newly developed two dimensional (2D)-segregation potential (SP) frost action model and provided distributions of ice contents, stresses, and displacements. Michalowski and Zhu (2006) [22] proposed a TM model using a porosity rate function within the theoretical framework of heaving soil depending on the stress state and temperature gradient (∇T). Liu and Yu (2011) [23] adopted the THM theory to simulate an unsaturated porous medium during frost heave. Kanie et al. (2012) [24] proposed a practical method for three-dimensional frost heave simulation based on Takashi's equation. Lai et al. (2014) [25] took the displacement, temperature, and porosity of the soil as variables and developed a coupled THM model to characterize the one-dimensional (1D) freezing process of saturated soil. Zheng and Kanie (2014) [26] developed a two-dimensional (2D) numerical THM model by expanding Takashi's equation from 1D to 2D to evaluate the frost heave process.
For simplification, several analysis methods of the TM and THM models of freezing soil under frost heave are based on the assumption that the compressive modulus (E s ) of the soil is a constant, or an empirical formula is formulated in which it is a function of temperature only [27]. As such, from a macroscopic mechanical perspective, the expansion associated with phase change during freezing causes compressive stresses to develop in the freezing zone of the soil. Accordingly, freezing-induced stress and strain are important features because they cause the upward swelling of soil and damage the underlying infrastructure [4]. Meanwhile, numerous factors, including the temperature gradient, overburden pressure, and cooling temperature, significantly influence the E s of freezing soil.
In engineering, a more effective analysis method is one that uses an appropriate stress-strain model with an accurate E s for freezing soil to better predict the maximum potential deformations and stresses resulting from frost heave at a satisfactory level of confidence. In this regard, the key research questions of a stress-strain analytical model should be investigated and an accurate E s value for freezing soil should be determined.
This study describes a series of one-dimensional multistage freezing experiments with silty clay columns. The results highlight the correlation between the freezing rate (U) and ∇T for each freezing soil layer, divided according to multistage freezing. In the regard, Takashi's equation is modified by substituting U with ∇T for operation over wider application conditions. Then, the determination of the compression modulus of freezing soil under the action of frost heave is proposed and a nonlinear elastic mechanical model is developed to map the stress-strain relationship of the freezing soil. Finally, the conclusions of the study are drawn.

Freezing Experiments
A set of one-dimensional freezing test experiments for silty clay columns were carried out in an open frost heave test system. The one-dimensional frost heave test system ( Figure 1) consisted of an insulated environmental chamber, a freezing cell, two cryostats with a temperature accuracy of ±0.05 • C, and various sensors necessary for measuring the soil temperature, frost heave, and water intake. The sample was placed in an organic cylindrical glass container with inner dimensions of 10 × 20 cm (diameter × height). To enhance the cooling efficiency, the circumference of the container was insulated by a 50-mm-thick layer of Styrofoam, and the container with the prepared specimens was placed in a temperature-controlled chamber at a constant temperature of −2 • C to avoid heat loss. Alcohol was used as the working fluid circulating through the plate's bottom and top ( Figure 2), with temperatures ranging from −40 • C to +40 • C, with an accuracy of ±0.1 • C. In addition, the water supply system was composed of a Mariotte bottle connected to the bottom plate through a plastic tube. In all experiments, an unpressurised water supply process was employed.   Frost-susceptible silt clay obtained from the city of Harbin in Northeast China was used to prepare the reconstituted samples. The grain size distribution of this clay is shown in Figure 3. The grain size characteristics of silt clay are as follows: curvature coefficient = 1.44; nonuniformity index    Frost-susceptible silt clay obtained from the city of Harbin in Northeast China was used to prepare the reconstituted samples. The grain size distribution of this clay is shown in Figure 3. The grain size characteristics of silt clay are as follows: curvature coefficient = 1.44; nonuniformity index = 11.85; plasticity index = 10.67%. All circular solid samples were 15 cm in height and 10 cm in diameter. These soil samples were compacted under slight vibration at an optimal moisture content of approximately 14% to achieve a density of approximately 1670 kg/m 3 . The soil samples were then saturated through the vacuum method, bringing the water content of the sample close to 23.59%.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 4 of 22 = 11.85; plasticity index = 10.67%. All circular solid samples were 15 cm in height and 10 cm in diameter. These soil samples were compacted under slight vibration at an optimal moisture content of approximately 14% to achieve a density of approximately 1670 kg/m 3 . The soil samples were then saturated through the vacuum method, bringing the water content of the sample close to 23.59%.  Table 1 lists the five freezing tests and their experimental conditions (i.e., test FT1 to FT5) in this study. The saturated soil column was subjected to a range of overburden stress. The tops of the samples were kept cold by subjecting them to a negative temperature, while the bottoms were kept warm by subjecting the samples to a positive temperature. After 24 h of consolidation, a cooling temperature (Tc) was applied to the top of the soil column, whereas the bottom of the sample was subjected to a warm temperature (Tw) of 2 °C. To characterize the natural freezing environment and the sequence of freezing temperatures, the temperature of the top portion (i.e., Tc) varied with time and underwent three freezing stages (i.e., stages 1, 2, and 3). The Tc values in stages 1, 2, and 3 were −5, −10, and −15 °C, respectively, and each freezing stage lasted 72 h. These identical thermal conditions are illustrated in Figure 4. For each sample, the temperature gradient was calculated from the temperature difference between the warm end and the cold end divided by the height of the sample. The initial temperatures were the temperatures of the samples before testing. A Data Taker 80 instrument (a robust data logger   Table 1 lists the five freezing tests and their experimental conditions (i.e., test FT1 to FT5) in this study. The saturated soil column was subjected to a range of overburden stress. The tops of the samples were kept cold by subjecting them to a negative temperature, while the bottoms were kept warm by subjecting the samples to a positive temperature. After 24 h of consolidation, a cooling temperature (T c ) was applied to the top of the soil column, whereas the bottom of the sample was subjected to a warm temperature (T w ) of 2 • C. To characterize the natural freezing environment and the sequence of freezing temperatures, the temperature of the top portion (i.e., T c ) varied with time and underwent three freezing stages (i.e., stages 1, 2, and 3). The T c values in stages 1, 2, and 3 were −5, −10, and −15 • C, respectively, and each freezing stage lasted 72 h. These identical thermal conditions are illustrated in Figure 4.

Experimental Results and Analysis
The evolution of frost heave (S) in the five samples under various overburden pressures recorded over time is illustrated in Figure 5. In the beginning of each freezing stage, a similar evolutionary trend of rapid frost heaving occurs against various overburden pressures, before becoming more gradual. As Tc reduces from one stage to the next, the effect of frost heave exhibits a slight yet abrupt increase at the beginning of stages 2 and 3. By the end of the freezing (at t = 216 h), the sizes of the soil columns due to frost heaving are 43.39, 36.78, 30.43, 26.25, and 23.46 mm under the overburden stresses of 0, 0.015, 0.028, 0.041, and 0.054 MPa, respectively. From the graph, as the overburden stress (σv) increases, the effect of frost heave decreases. The reason for this is that an increasing σv inhibits the amount of water entering the soil. The temperature Tsoil appears to be insensitive to the overburden stress σv at a low level. As an example, the Tsoil at various depths in test FT1 is shown in Figure 6. The temperature field within the soil column under an identical thermal condition is similar for various overburden stresses. For each freezing stage, the temperature Tsoil at various depths follows a similar trend and requires different spans of freezing time t to attain a stable state.  For each sample, the temperature gradient was calculated from the temperature difference between the warm end and the cold end divided by the height of the sample. The initial temperatures were the temperatures of the samples before testing. A Data Taker 80 instrument (a robust data logger featuring USB memory stick support, built-in display, and extensive communications capability, which can be used for environmental monitoring) was employed to collect temperature, displacement, and stress data. During the testing, a displacement transducer was installed at the plate to monitor the effects of frost heave on the specimen's surface. The side of the container included four holes with four thermistors for each sample. Each thermistor was placed 3 cm apart along the depth to measure the internal soil temperature (T soil ), and two additional thermistors were placed on the top and bottom surfaces of the sample.

Experimental Results and Analysis
The evolution of frost heave (S) in the five samples under various overburden pressures recorded over time is illustrated in Figure 5. In the beginning of each freezing stage, a similar evolutionary trend of rapid frost heaving occurs against various overburden pressures, before becoming more gradual. As T c reduces from one stage to the next, the effect of frost heave exhibits a slight yet abrupt increase at the beginning of stages 2 and 3. By the end of the freezing (at t = 216 h), the sizes of the soil columns due to frost heaving are 43.39, 36.78, 30.43, 26.25, and 23.46 mm under the overburden stresses of 0, 0.015, 0.028, 0.041, and 0.054 MPa, respectively. From the graph, as the overburden stress (σ v ) increases, the effect of frost heave decreases. The reason for this is that an increasing σ v inhibits the amount of water entering the soil.

Experimental Results and Analysis
The evolution of frost heave (S) in the five samples under various overburden pressures recorded over time is illustrated in Figure 5. In the beginning of each freezing stage, a similar evolutionary trend of rapid frost heaving occurs against various overburden pressures, before becoming more gradual. As Tc reduces from one stage to the next, the effect of frost heave exhibits a slight yet abrupt increase at the beginning of stages 2 and 3. By the end of the freezing (at t = 216 h), the sizes of the soil columns due to frost heaving are 43.39, 36.78, 30.43, 26.25, and 23.46 mm under the overburden stresses of 0, 0.015, 0.028, 0.041, and 0.054 MPa, respectively. From the graph, as the overburden stress (σv) increases, the effect of frost heave decreases. The reason for this is that an increasing σv inhibits the amount of water entering the soil. The temperature Tsoil appears to be insensitive to the overburden stress σv at a low level. As an example, the Tsoil at various depths in test FT1 is shown in Figure 6. The temperature field within the soil column under an identical thermal condition is similar for various overburden stresses. For each freezing stage, the temperature Tsoil at various depths follows a similar trend and requires different spans of freezing time t to attain a stable state.   The temperature T soil appears to be insensitive to the overburden stress σ v at a low level. As an example, the T soil at various depths in test FT1 is shown in Figure 6. The temperature field Appl. Sci. 2020, 10, 7753 6 of 22 within the soil column under an identical thermal condition is similar for various overburden stresses. For each freezing stage, the temperature T soil at various depths follows a similar trend and requires different spans of freezing time t to attain a stable state.

Modified Takashi's Equation
In these experiments, the freezing front in stage 1 gradually moves downward with time until the freezing rate U decreases to zero and the temperature T soil becomes stable. In this stage, the soil samples can be divided into the freezing and unfrozen zones, as shown in Figure 8a. The temperature (T fre,c1 ) at the top of the first freezing soil layer equals −5 • C (i.e., T fre,c1 = T c = −5 • C at that moment).

Modified Takashi's Equation
In these experiments, the freezing front in stage 1 gradually moves downward with time until the freezing rate U decreases to zero and the temperature Tsoil becomes stable. In this stage, the soil samples can be divided into the freezing and unfrozen zones, as shown in Figure 8a. The temperature (Tfre,c1) at the top of the first freezing soil layer equals −5 °C (i.e., Tfre,c1 = Tc = −5 °C at that moment). Similarly, the freezing front continues to move downward with time. The soil samples in stage 2 can be divided into frozen, freezing, and unfrozen zones (Figure 8b). The frozen zone remains unchanged with time, the second freezing soil layer gradually moves downward, and the unfrozen zone gradually shrinks. Meanwhile, the temperature (Tfre,c2) at the top of the second freezing soil layer under Tc = −10 °C is −1.6 °C.
The soil samples in stage 3 can also be divided into frozen, freezing, and unfrozen zones, as shown in Figure 8c. Similarly, the temperature (Tfre,c3) at the top of the third freezing soil layer under Tc = −15 °C is −0.6 °C. Overall, three successive freezing soil layers (i.e., the first, second, and third freezing soil layers) appeared from stage 1 to stage 3 during the testing.  Similarly, the freezing front continues to move downward with time. The soil samples in stage 2 can be divided into frozen, freezing, and unfrozen zones (Figure 8b). The frozen zone remains unchanged with time, the second freezing soil layer gradually moves downward, and the unfrozen zone gradually shrinks. Meanwhile, the temperature (T fre,c2 ) at the top of the second freezing soil layer under T c = −10 • C is −1.6 • C.
The soil samples in stage 3 can also be divided into frozen, freezing, and unfrozen zones, as shown in Figure 8c. Similarly, the temperature (T fre,c3 ) at the top of the third freezing soil layer under T c = −15 • C is −0.6 • C. Overall, three successive freezing soil layers (i.e., the first, second, and third freezing soil layers) appeared from stage 1 to stage 3 during the testing.
When the temperature field changes from t to t + dt, the increment (dS) of frost heave at the top of the soil sample is equal to the sum of the frost heave increments of the frozen soil layer (dS fro ) and the freezing soil layer (dS fre ). Therefore: dS = dS fro + dS fre (1) As a result of the existence of the frozen fringe, which has a scant amount of unfrozen water and low permeability, water cannot migrate from the unfrozen area to the frozen area. Therefore, we assumed that frozen soils do not permeate water migration and that the water migration occurs only in the freezing soil layers. Thus, dS fro = 0 in Equation (1). Then: This means the increment of frost heave is mainly derived from the degree of frost heave in the freezing soil layer.
The increments of frost heave in the freezing soil layer (i.e., frost heave increments at the surfaces of soil samples) during each freezing stage are represented as dS fre,1 , dS fre,2 , and dS fre,3 , respectively. These values are calculated with Equation (2) and the experimental results in Figure 8, and are plotted in Figure 9. From this calculation, the increments of frost heave generated in the freezing zone eventually stabilize with time.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 8 of 22 When the temperature field changes from t to t + dt, the increment (dS) of frost heave at the top of the soil sample is equal to the sum of the frost heave increments of the frozen soil layer (dSfro) and the freezing soil layer (dSfre). Therefore: As a result of the existence of the frozen fringe, which has a scant amount of unfrozen water and low permeability, water cannot migrate from the unfrozen area to the frozen area. Therefore, we assumed that frozen soils do not permeate water migration and that the water migration occurs only in the freezing soil layers. Thus, dSfro = 0 in Equation (1). Then: This means the increment of frost heave is mainly derived from the degree of frost heave in the freezing soil layer.
The increments of frost heave in the freezing soil layer (i.e., frost heave increments at the surfaces of soil samples) during each freezing stage are represented as dSfre,1, dSfre,2, and dSfre,3, respectively. These values are calculated with Equation (2) and the experimental results in Figure 8, and are plotted in Figure 9. From this calculation, the increments of frost heave generated in the freezing zone eventually stabilize with time.   The increments of hf during stages 1, 2, and 3 (dh1, dh2, and dh3, respectively), which can be calculated using the results in Figure 6, are presented in Figure 10. The frost heave strain (ε) can be defined as the ratio of frost heave to the freezing soil layer thickness (China Construction Ministry, China, 1998). As such: where εi, dSfre,i, dhi are the frost heave strain, frost heave increment, and thickness of the ith freezing soil layer during stage i (where i = 1, 2 and 3). Thus, εi can be calculated with Equation (3), the calculation results for which are shown in Figure 11.   The increments of h f during stages 1, 2, and 3 (dh 1 , dh 2 , and dh 3 , respectively), which can be calculated using the results in Figure 6, are presented in Figure 10. The frost heave strain (ε) can be defined as the ratio of frost heave to the freezing soil layer thickness (China Construction Ministry, China, 1998). As such: where ε i , dS fre,i , dh i are the frost heave strain, frost heave increment, and thickness of the ith freezing soil layer during stage i (where i = 1, 2 and 3). Thus, ε i can be calculated with Equation (3), the calculation results for which are shown in Figure 11. The increments of hf during stages 1, 2, and 3 (dh1, dh2, and dh3, respectively), which can be calculated using the results in Figure 6, are presented in Figure 10. The frost heave strain (ε) can be defined as the ratio of frost heave to the freezing soil layer thickness (China Construction Ministry, China, 1998). As such: where εi, dSfre,i, dhi are the frost heave strain, frost heave increment, and thickness of the ith freezing soil layer during stage i (where i = 1, 2 and 3). Thus, εi can be calculated with Equation (3), the calculation results for which are shown in Figure 11.       The ∇Ti within the ith freezing soil layer during stage i can be calculated by ∇Ti = (T0,i -Tfre,ci)/dhi (4) where T0,i is the temperature of the freezing front in the ith freezing soil layer in stage i. Therefore, the relationship between εi and ∇Ti can be further determined, which is presented in Figure 12. From Figures 11 and 12, the ε of the freezing soil layer under various overburden stresses and the ∇Ti tend to become stable with time. We found that a lower Tfre,ci induces a higher ε.   The ∇T i within the ith freezing soil layer during stage i can be calculated by where T 0,i is the temperature of the freezing front in the ith freezing soil layer in stage i. Therefore, the relationship between ε i and ∇T i can be further determined, which is presented in Figure 12. From Figures 11 and 12, the ε of the freezing soil layer under various overburden stresses and the ∇T i tend to become stable with time. We found that a lower T fre,ci induces a higher ε. The ∇Ti within the ith freezing soil layer during stage i can be calculated by ∇Ti = (T0,i -Tfre,ci)/dhi (4) where T0,i is the temperature of the freezing front in the ith freezing soil layer in stage i. Therefore, the relationship between εi and ∇Ti can be further determined, which is presented in Figure 12. From Figures 11 and 12, the ε of the freezing soil layer under various overburden stresses and the ∇Ti tend to become stable with time. We found that a lower Tfre,ci induces a higher ε.   Because of its satisfactory performance and simplicity, Takashi's equation [28] is widely adopted as the theoretical framework to estimate ε, which is related to σv and U, as shown below: where U is in m/s, and ε0, σ0, and U0 are constants representing the soil characteristics obtained by a regulated experiment of the Japanese Geotechnical Society (JGS; JGS 0171-2003); U and σv in the freezing direction are two critical variables for determining ε. According to Equation (5), if σv is infinitesimally small, ε will have an infinitely large value. Therefore, the value of σv is limited between 1 and 15 kg/cm 2 (i.e., 98 and 1470 kPa) for Takashi's equation. Beyond this limit, the equation ceases to work. Based on their definition, the values for Ui and ∇Ti with the time obtained from Figure 6 are shown in Figure 13. It is noted that the tendency between Ui and ∇Ti is highly consistent. In view of this consistency between Ui and ∇Ti under the different Tfre,ci values, we introduced Tfre,ci and ∇Ti into the Takashi model and accordingly modified σv. Takashi's equation was re-constructed to allow for σv at a low level (even zero, i.e., without overburden pressure), as shown below.
where parameters A, B, C, and D are coefficients that are back calculated using data from Figure 12.   Because of its satisfactory performance and simplicity, Takashi's equation [28] is widely adopted as the theoretical framework to estimate ε, which is related to σ v and U, as shown below: where U is in m/s, and ε 0 , σ 0 , and U 0 are constants representing the soil characteristics obtained by a regulated experiment of the Japanese Geotechnical Society (JGS; JGS 0171-2003); U and σ v in the freezing direction are two critical variables for determining ε. According to Equation (5), if σ v is infinitesimally small, ε will have an infinitely large value. Therefore, the value of σ v is limited between 1 and 15 kg/cm 2 (i.e., 98 and 1470 kPa) for Takashi's equation. Beyond this limit, the equation ceases to work. Based on their definition, the values for U i and ∇T i with the time obtained from Figure 6 are shown in Figure 13. It is noted that the tendency between U i and ∇T i is highly consistent. In view of this consistency between U i and ∇T i under the different T fre,ci values, we introduced T fre,ci and ∇T i into the Takashi model and accordingly modified σ v . Takashi's equation was re-constructed to allow for σ v at a low level (even zero, i.e., without overburden pressure), as shown below.
where parameters A, B, C, and D are coefficients that are back calculated using data from Figure 12. The outline of the procedure is as follows: (1) the three corresponding values of ε can be obtained once the three discrete moments (i.e., ∇T1 = 0.052, ∇T1 = 0.064, and ∇T1 = 0.080 for the first freezing soil layer; ∇T2 = 0.078, ∇T2 =0.085, and ∇T2 =0.087 for the second freezing soil layer; and ∇T3 = 0.110, ∇T3 = 0.112, and ∇T3 = 0.113 for the third freezing soil layer) have been assigned; (2) the values of εi at the specified intervals are respectively substituted into Equations (7)-(9); (3) for each freezing soil layer, the parameters A, B, C, and D can be determined using the least squares method.
Therefore, for the first freezing soil layer: The outline of the procedure is as follows: (1) the three corresponding values of ε can be obtained once the three discrete moments (i.e., ∇T 1 = 0.052, ∇T 1 = 0.064, and ∇T 1 = 0.080 for the first freezing soil layer; ∇T 2 = 0.078, ∇T 2 =0.085, and ∇T 2 =0.087 for the second freezing soil layer; and ∇T 3 = 0.110, ∇T 3 = 0.112, and ∇T 3 = 0.113 for the third freezing soil layer) have been assigned; (2) the values of ε i at the specified intervals are respectively substituted into Equations (7)-(9); (3) for each freezing soil layer, the parameters A, B, C, and D can be determined using the least squares method. Therefore, for the first freezing soil layer: Similarly, for the second freezing soil layer: In addition, for the third freezing soil layer: In conclusion, Equations (7)-(9) can be further expressed as Equation (6) for the ith freezing soil layer: where A i , B i , C i and D i are parameters related to T fre,ci : Here, ε i can be calculated with Equation (10). As such, Figure 14 shows the comparison between the experimental and estimated results, indicating that Equation (10) In conclusion, Equations (7)-(9) can be further expressed as Equation (6) for the i th freezing soil layer: where Ai, Bi, Ci and Di are parameters related to Tfre,ci: Here, εi can be calculated with Equation (10). As such, Figure 14 shows the comparison between the experimental and estimated results, indicating that Equation (10)

Mechanical Stress-Strain Relationships in the Freezing Soil Layers
In this study, the following items were assumed in order to formulate the mathematical expressions: 1. The self-weight consolidation is finished before the testing and the compressibility of the soil particles is negligible; 2. The density change of the soil is negligible in the processes of soil freezing; 3. The soil is isotropic and behaves as a nonlinear elastic material; 4. Water migration occurs only when the temperature field in the soil changes. When the temperature field is stable, the water migration stops.

Definition of Free Frost Heaving Strain
Herein, the concept of free frost heaving strain, εe,i, in the freezing soil layer in stage i is

Mechanical Stress-Strain Relationships in the Freezing Soil Layers
In this study, the following items were assumed in order to formulate the mathematical expressions: 1.
The self-weight consolidation is finished before the testing and the compressibility of the soil particles is negligible; 2.
The density change of the soil is negligible in the processes of soil freezing; 3.
The soil is isotropic and behaves as a nonlinear elastic material; 4.
Water migration occurs only when the temperature field in the soil changes. When the temperature field is stable, the water migration stops.

Definition of Free Frost Heaving Strain
Herein, the concept of free frost heaving strain, ε e,i , in the freezing soil layer in stage i is introduced, which becomes ε when σ v equals 0. Hence, ε e,i from Equation (10) may be expressed as follows: As a result, the reduction of ε (i.e., ∆ε i ) due to the increase in σ v can obtained from: By combining Equations (10) and (11) with Equation (12), we derived ∆ε i as: A similar hyperbolic curve is utilized to fit the relationship between ∆ε i and σ v within the freezing soil layer during stage i.

Potential Frost Heaving Stress
The concept of the potential frost heaving stress (σ e,i ) in stage i is introduced, which is equal to the critical σ v required to fully prevent the generation of frost heave in the freezing soil layer, i.e., ε i = 0. Therefore, the value of σ e,i can be determined using Equation (10) as follows: Equation (14) gives the relationship between σ e,i , T fre,ci , and ∇T i within the ith freezing soil layer, which is indicated in Figure 15. Here, σ e,i gradually reduces as T fre,ci decreases, but it appears to be independent of σ v .
Appl. Sci. 2020, 10, x FOR PEER REVIEW 16 of 22 As a result, the reduction of ε (i.e., Δεi) due to the increase in σv can obtained from: By combining Equations (10) and (11) with Equation (12), we derived Δεi as: A similar hyperbolic curve is utilized to fit the relationship between Δεi and σv within the freezing soil layer during stage i.

Potential Frost Heaving Stress
The concept of the potential frost heaving stress (σe,i) in stage i is introduced, which is equal to the critical σv required to fully prevent the generation of frost heave in the freezing soil layer, i.e., εi = 0. Therefore, the value of σe,i can be determined using Equation (10) as follows: Equation (14) gives the relationship between σe,i, Tfre,ci, and ∇Ti within the i th freezing soil layer, which is indicated in Figure 15. Here, σe,i gradually reduces as Tfre,ci decreases, but it appears to be independent of σv.

Decomposition of Potential Frost Heaving Stress
Next, σe,i can be decomposed into two components and is expressed as: where σf,i is the effective frost heaving stress that reflects the actual amount of frost heave.
Accordingly, εi can be divided into two parts, as follows: where εr,i is defined as the resilient strain, which is generated due to σv in equilibrium with a portion of σe,i. Here, εf,i is termed as the real ε induced only by σf,i.

Decomposition of Potential Frost Heaving Stress
Next, σ e,i can be decomposed into two components and is expressed as: where σ f,i is the effective frost heaving stress that reflects the actual amount of frost heave. Accordingly, ε i can be divided into two parts, as follows: where ε r,i is defined as the resilient strain, which is generated due to σ v in equilibrium with a portion of σ e,i . Here, ε f,i is termed as the real ε induced only by σ f,i . Based on this concept, the resilience modulus (E r,i ) is expressed as: Similarly, the compression modulus (E c,i ), which only relates to σ f,i , can be defined as follows: Thus, the following formula can easily be obtained by combining Equations (16)- (18): Equation (15) is then substituted into Equation (19): By defining: Equation (20) can be rewritten as:

Determination of Compression Modulus
According to the definition of ε e,i and Equation (22), ε i = ε e,i can be obtained when σ v = 0, and Equation (22) can be rewritten as: Combining Equations (11), (14), and (23) yields: Equation (24) represents the final expression of E c,i as a function of T fre,ci and ∇T i .

Determination of the Resilience Modulus
By substituting Equation (23) into Equation (22), ε i can be further expressed as: Furthermore, Equation (25) can be rewritten as: By substituting Equations (10) and (11) into Equation (26), we get: The value of E r,i can be easily evaluated by combined Equations (21), (23), and (27). Accordingly, E r,i depends on T fre,ci , ∇T i , and σ v .

Mathematical Equations for Stress-Strain Model of the Freezing Soil
By combining the model with Equations (15) and (22), these equations can be rewritten as: Assuming:

Discussion
The frost heaving process of silty clay is a complex problem. The amount of moisture transfer in the soil is directly related to the frost heave increment on the specimen's surface. Firstly, a set of one-

Discussion
The frost heaving process of silty clay is a complex problem. The amount of moisture transfer in the soil is directly related to the frost heave increment on the specimen's surface. Firstly, a set of one-dimensional freezing experiments were carried out to research the fundamental law of frost heaving of saturated clay. Temperature and load are the main factors that determine the frost heaving characteristics of the soil. In order to reveal the development law and mechanism of soil frost heaving, it is necessary to deeply analyze the internal relationship between soil frost heaving characteristics and temperature distribution during soil freezing.
Takashi's equation was generalized so as to play a significant role in this research. It has been widely used as based on Japanese Geotechnical Society (JGS) regulations and is valuable for engineering purposes. According to Equation (5), if σ v is infinitesimally small, ε will have an infinitely large value. Therefore, Takashi's model is limited to cases with overburden stress larger than 98 kPa and smaller than 1470 kPa. It can be seen that the model does not take the influencing factors of the frost heave rate into account adequately. Necessary modifications should be carried out in order to operate in wider application conditions. In Takashi's equation, the freezing rate and effective constraining stress in the freezing direction are the factors dominating frost heave prediction. As shown in Figures 6 and 13, the U i is directly related to the ∇T i . Thus, the T fre,ci and ∇T i were introduced into the Takashi model, as shown in Equation (6). Then, Takashi's equation can be used in cases with low σ v , even without overburden pressure.
On the other hand, the analysis of soil deformation induced by frost heaving is another important target in the study. For this, the method of determining the compression modulus of freezing soil is given and a nonlinear elastic mechanical model is developed to map the stress-strain relationship of the freezing soil.
In future studies, the stress-strain of freezing soil must be researched in action conditions, such as in freezing-thawing cycles. Then, a more realistic stress-strain model that includes elaborate factors could be developed.

Conclusions
A stress-strain model for freezing silty clay under frost heave based on a modified Takashi's equation was proposed. The conclusions of the study and recommendations for future studies are as follows: (1) Surface heaving significantly increases with freezing time in each freezing stage. Furthermore, frost heaving decreases with the increase of σ v , which inhibits frost heave. Therefore, the influence of σ v should be considered when estimating the stress-strain behavior and mechanical properties of freezing soil. In addition, frost heave increases slightly with the decrease in the short-duration cooling temperature from one stage to the other; (2) By using the frost heave data, the correlation between the freezing rate (U i ) and temperature gradient (∇T i ) within each freeing soil layer was investigated. Then, Takashi's equation was modified by substituting U i with ∇T i and evaluating the frost heave strain within each freezing subsoil layer; (3) A nonlinear mechanical elasticity model with the compressive modulus as a function of the variables of temperature gradient, overburden pressure, and cooling temperature was proposed to characterize the stress-strain relationship of freezing soil resulting from the effect of frost heave; (4) Future studies should address other important factors and introduce them into this stress-strain model to achieve a more general understanding of the realistic stress-strain behavior of freezing soil resulting from the effects of frost heave. Meanwhile, researchers can use the proposed model with the finite element method to simulate more complex practical tests.