The Modelling of Freezing Process in Saturated Soil Based on the Thermal-Hydro-Mechanical Multi-Physics Field Coupling Theory

Abstract

The freezing process of saturated soil is studied under the condition of water replenishment. The process of soil freezing was simulated based on the theory of the energy and mass conservation equations and the equation of mechanical equilibrium. The accuracy of the model was verified by comparison with the experimental results of soil freezing. One-side freezing of a saturated 10-cm-high soil column in an open system with different parameters was simulated, and the effects of the initial void ratio, hydraulic conductivity, and thermal conductivity of soil particles on soil frost heave, freezing depth, and ice lenses distribution during soil freezing were explored. During the freezing process, water migrates from the warm end to the frozen fringe under the actions of the temperature gradient and pore pressure. During the initial period of freezing, the frozen front quickly moves downward, the freezing depth is about 5 cm after freezing for 30 h, and the final freezing depth remains about 6 cm. The freezing depth of the soil column is affected by soil porosity and thermal conductivity, but the final freezing depth mainly depends on the temperatures of the top and lower surfaces. The frost heave is mainly related to the amount of water migration. The relationship between the amount of frost heave and the hydraulic conductivity is positively correlated, and the thickness of the stable ice lens is greatly affected by the hydraulic conductivity. With the increase of the hydraulic conductivity and initial void ratio, the formation of ice lenses in the soil become easier. With the increase of the initial void ratio and thermal conductivity of soil particles, the frost heave of the soil column also increases. With high-thermal-conductivity soil, the formation of ice lenses become difficult.

1. Introduction

Permafrost in China is mainly distributed in high latitudes. With the change of temperature, water and ice in the pores of frozen soil can transform into each other. Frozen soil is a dynamic system related to temperature, and the dynamic system involves heat transfer, material exchange, and soil deformation. Therefore, the structure, composition, and mechanical properties of frozen soils are more complicated than unfrozen soil [1,2]. The most important problem in permafrost engineering is frost heave, especially when the water content of soil is high. This has caused great damage to the infrastructure in cold regions. For example, frost heave can cause railroad undulation, the rupture of oil pipeline, and instability of the electric tower [3,4,5]. Therefore, exploring the mechanism of the freezing process in soils is of great significance to permafrost engineering. Many scholars have carried out a lot of studies on the soil freezing process by experimental and numerical simulation methods.

Taber [6] studied the soil freezing process by an experimental method to explain the frost heave of soil and the formation of ice lenses. Later, Everett [7] established the hydrothermal model of frozen soils through capillary theory, which is the first frost theory. This theory was used to explain the frost heave phenomenon and estimates frost heave. However, Everett did not verify the model through experimental data, and the theory cannot explain the formation of discontinuous ice lenses. When the soil freezes, the heat of soil migrates from the high-temperature region to thenlow-temperature region, and the content of unfrozen water in frozen soils drops sharply [8]. The soil freezing is a dynamic coupling process under the influence of the temperature gradient and hydraulic gradient [9]. Harlan [8] established a hydrothermal coupling model for partially frozen soil based on the analogy of the mechanism of water transport in partially frozen soils and unsaturated soils. However, this model cannot explain the formation of ice lenses. Kay and Groenevelt [10] proposed the appropriate energy equation and combined it with the Clapeyron equation to find a better transport coefficient. Accordingly, Kay and Groenevelt [10] proposed a theory of water transfer under a temperature gradient and heat transfer under a water pressure gradient. These early theories about soil freezing focused on the temperature field and the seepage field, and ignored the role of the stress field during the freezing process. Additionally, the formation of the ice lens has not been explained well.

Assuming that there is a region of low moisture content and low moisture conductivity between the bottom of the ice lens and the frozen front, Miller [11] established the second frost theory. With a large amount of studies about coarse-textured soil, Gilpin [12] developed a model of frost heave and ice lenses. The model assumed that the latency of pore water release at the isotherms where most of the pore water freezes or an ice lens is forming, and the water in the unfrozen water film is driven entirely by normal driven pressure. According to the experimental results, Konrad and Morgenstern [13] presented that seepage velocity is continuous across the region of the fringe under steady-state conditions, and used a model that distinguishes the passive and active regions to explain the frost heave characteristics of fine-textured soils. The model confirms the correctness of the Clausius–Clapeyron equation at the bottom of ice lenses. Nishimura et al. [14] used a thermo-hydro-mechanical finite element formulation to simulate the pipeline heave. Tan et al. [15] established a thermo-hydro coupling model and focused on the evolution of the temperature field and simulated the freezing process of the soil column without water replenishment. Wu et al. [16] believes that water activity is the inducement of phase change between water and ice. According to the coupling relationship between the formation and development of ice and the flow of water and heat, Wu et al. [16] established a kinetic model of ice development. Song et al. [17] and Xu et al. [18] used the lattice Boltzmann method (LBM) model to predict the distribution of water content and the unfrozen water content and hydraulic conductivity during freezing, respectively. The theory of soil freezing has been relatively well developed, and some models have also considered the ice lenses. However, most models do not consider the effects of soil deformation on heat transfer and water migration, and the effect of the horizontal displacement limit, which caused the model to describe the actual physical process of soil freezing inaccurately.

Soil freezing is a thermo-hydro-mechanical coupling process caused by a temperature gradient, and the coupling relationship of soil freezing can be shown in Figure 1. Due to the adsorption of the soil particles, the water film between the particles is kept in a liquid state when the temperature is lower than the freezing temperature. When the temperature continues to decrease, the water film begins to freeze, and with the ice wedging between the particles, the soil particles are separated. Under the appropriate condition of water replenishment and temperature, the ice layer between the soil particles becomes thicker to form the ice lens [19]. It has been widely accepted by scholars that phase change of the pore water and the formation of ice lens occur in the so-called ‘frozen fringe’ and the stress criterion is considered as the criteria for the formation of a new ice lens [12,20,21,22]. Thomas et al. [21] thought that a new ice lens forms when the pore pressure exceeds the sum of the separation strength and the overburden stress of the freezing soil. Cao and Liu [23] found that if the pore pressure criterion is used, the distribution of ice pressure will be discontinuous, and there are errors between the simulation results and experimental data. The simulation results using the ice pressure criterion are closer to the experimental data.

In this study, a new thermo-hydro-mechanical coupling model is established to describe the freezing process in soil. The coupling relationship between the physical fields is shown in Figure 1. The effects of the initial void ratio, hydraulic conductivity, and thermal conductivity of soil particles on soil freezing process are explored. A new condition for the formation of ice lenses is proposed. In this model, the following items are assumed:

• (1). The pores of soil are filled with water and ice during freezing, and the soil is regarded as isotropic and elastic mediums.

• (2). The pore ice is immobile relative to the soil skeleton.

• (3). The soil particles, pore ice, and water are incompressible.

• (4). A weighting algorithm is used to convert the soil of the three-phase system into a single-phase system for the calculation of heat transfer.

2. Theoretical Model of Soil Freezing

According to the theory of adsorption, the pore water is not frozen completely when the temperature is below the freezing temperature, and the unfrozen water adsorbs on the surface of the soil particles, and the thickness of the water film mainly depends on the temperature [24]. This process can be illustrated in Figure 2.

2.1. Heat Transfer

In the process of soil freezing, heat transfer and heat convection occur according to the thermodynamics theory. The evolution of the temperature field conforms to the law of conservation of energy. The heat conduction equation for soil freezing can be expressed as:

(1)$\frac{\partial \left(\mathsf{\Phi }dV\right)}{\partial t}+\nabla \left(-\lambda \nabla T\right)dV+{\rho }_l{c}_l\nabla T·\overrightarrow{v}dV=0$

where $t$ is time; T is temperature; $c_l$ and ${\rho }_l$ are the specific heat capacities and the density of pore water, respectively; $\nabla$ is the Hamiltonian operator and means $\nabla =\frac{\partial }{\partial x}+\frac{\partial }{\partial y}$ (in two dimensions); $\mathsf{\Phi }$ represents the heat content of soil per unit volume; $d$ stands for the differential operator; $dV$ refers to the volume element of the soils; $\lambda$ refers the thermal conductivity; and $\overrightarrow{v}$ is the seepage velocity of pore water:

(2)$\mathsf{\Phi }=H\left(T-T_r\right)-Ln{S}_i{\rho }_i$

where $L$ is the latent heat of water phase change; $n$ represents the porosity of the soil; $S_i$ is the ice volume percentage in pores; $T_r$ is the reference temperature; ${\rho }_i$ is the density of ice; $H$ is the volumetric heat capacity of soil, which can be expressed as:

(3)$H=\left(1-n\right){\rho }_s{c}_s+\left(1-S_i\right)n{\rho }_l{c}_l+S_{i}n{\rho }_i{c}_i$

where ${\rho }_s$ is the density of the soil particles, and $c_s$ and $c_i$ are the specific heat capacities of the soil particles and ice, respectively. Tice et al. [25] established a relationship between the ice percentage $S_i$ and temperature based on the results of experiments as follows:

(4)$S_i=\{\begin{array}{cc}1-{\left[1-\left(T-T_0\right)\right]}^{\alpha }& T\le T_0\\ 0& T>T_0\end{array}$

where T₀ is the freezing point temperature (273.15 K) of water in pores ignoring the effect of the water vapor sorption isotherm [21,22,25,26,27], and here $\alpha$ is a parameter depending on the pore size of soil:

(5)$dV=\left(1+e\right)d{V}_s$

where $e$ and $d{V}_s$ are the porosity ratio and the soil particle volume of soil element, respectively:

(6)$\lambda ={\lambda }_s{}^{\left(1-n\right)}{\lambda }_l{}^{\left(1-S_i\right)n}{\lambda }_i{}^{n{S}_i}$

Substituting Equations (2), (3), and (5) into Equation (1) and assuming that the solid particles are incompressible, the equations of heat transfer can be written as:

$C_{11}\frac{\partial T}{\partial t}+C_{12}\frac{1}{1+e}\frac{\partial e}{\partial t}+\nabla \left(-\lambda \nabla T\right)+{\rho }_l{c}_l\nabla T·\overrightarrow{v}=0$

(7)$C_{11}=\left[H+\frac{e}{1+e}\left(T-T_r\right)\left({\rho }_i{c}_i-{\rho }_l{c}_l\right)\frac{\partial S_i}{\partial T}-\frac{e}{1+e}L{\rho }_i\frac{\partial S_i}{\partial T}\right]$

$C_{12}=\left[\left(1-S_i\right){\rho }_l{c}_l+S_i{\rho }_i{c}_i\right]\left(T-T_r\right)-L{S}_i{\rho }_i$

2.2. Mechanical Equilibrium

During the freezing process of the soil, the skeleton of the soil is filled with pore ice and water. The total stress of frozen soil is equal the sum of the skeleton stress and pore pressure. According to Biot’s effective stress law [28], the equilibrium equation can be expressed as:

(8)${\sigma }_{ij,j}^e+{\alpha }^{*}{P}_{,i}+f_i=0$

where ${\sigma }_{ij}^e$ is the effective stress tensor; ${\alpha }^{*}$ is the Biot coefficient; $P$ is the pore pressure; and $f_i$ is the body force in the ith direction:

(9)${\alpha }^{*}=1-\frac{K}{K_s}$

where $K$ is the bulk modulus of the soil skeleton, and $K_s$ is the bulk modulus of soil particles.

The total stress ${\sigma }_{11}$ in the vertical direction can be expressed as:

(10)${\sigma }_{11}=P_{ob}+{{\displaystyle \int }}_x^h\gamma d{x}^{*}$

where $P_{ob}$ is the overburden pressure of the soil column; h is the height of the soil column; and $\gamma$ is the unit weight.

According to the composition of frozen soil, the unit weight $\gamma$ can be calculated by the following equation:

(11)$\gamma =g\left[\left(1-n\right){\rho }_s+n{S}_i{\rho }_i+n\left(1-S_i\right){\rho }_l\right],$

where $g$ is the gravitational acceleration.

Therefore, the initial total stress ${\sigma }_{11}^0$ can be written as:

(12)${\sigma }_{11}^0=P_{ob}+{{\displaystyle \int }}_x^h{\gamma }_{0}d{x}^{*}$

and the initial unit weight ${\gamma }_0$ can be given as:

(13)${\gamma }_0=g\left[\left(1-n_0\right){\rho }_s+n_0{\rho }_l\right]$

where $n_0$ is the initial porosity, $n_0=\frac{e_0}{1+e_0}$, and $e_0$ is the initial porosity ratio.

Based on the principle of effective stress, the relationship of the total stress ${\sigma }_{11}$, the effective stress ${\sigma }_{11}^e$, and pore pressure $P$ can be written as:

(14)${\sigma }_{11}={\sigma }_{11}^e+P$

For isotropic porous media materials, the relationship between the effective stress tensor and the strain tensor can be expressed as:

(15)${\sigma }_{ij}^e={\lambda }^{*}{\epsilon }_{kk}{\delta }_{ij}+2{\mu }^{*}{\epsilon }_{ij}+\left({\sigma }_{ij}^0-P_0\right)$

where ${\lambda }^{*}$ and ${\mu }^{*}$ are the Lame constant; ${\epsilon }_{kk}$ is the volume strain and ${\epsilon }_{kk}={\epsilon }_{11}+{\epsilon }_{22}+{\epsilon }_{33}$; ${\epsilon }_{11}, {\epsilon }_{22}, {\epsilon }_{33}$ are the normal strain; ${\delta }_{ij}$ is the Kronecker symbol; ${\sigma }_{ij}^0$ is the initial total stress tensor; and $P_0$ is the initial pore pressure.

The freezing process of the soil layer is simulated, and the schematic diagram of the soil column model is shown in Figure 3. When shallow soil freezes in the natural environment, frost heave does not occur in the horizontal direction. From Figure 3, it can be seen that the horizontal strains of the model meet the condition ${\epsilon }_{22}={\epsilon }_{33}=0$ and ${\epsilon }_{kk}={\epsilon }_{11}$. Equation (15) can be written as:

(16)${\sigma }_{11}^e={\lambda }^{*}{\epsilon }_{11}+2{\mu }^{*}{\epsilon }_{11}+{\sigma }_{11}^0$

Soil particles are assumed to be undeformed, and the frost heave is derived from the expansion of the pores. Before the deformation, the void ratio was $e_0$; the volume of soil particles is A; and the volume of pores is (1 − A). After the deformation, the void ratio is $e$; the volume of soil particles is A; and the volume of pores is Q. The following relationships can be obtained:

(17)$e_0=\frac{1-A}{A}$

(18)$e=\frac{Q}{A}$

The volumetric strain can be expressed as:

(19)${\epsilon }_{kk}={\epsilon }_{11}=\frac{Q-\left(1-A\right)}{1}$

Substituting Equations (17) and (18) into Equation (19), the ${\epsilon }_{11}$ can be written as:

(20)${\epsilon }_{kk}={\epsilon }_{11}=\frac{e-e_0}{1+e_0}$

(21)$P={\sigma }_{11}-{\sigma }_{11}^e$

Substituting Equations (10), (12), (16), and (20) into Equation (21), the pore pressure equation can be rewritten as:

(22)$P={{\displaystyle \int }}_x^h\left(\gamma -{\gamma }_0\right)d{x}^{*}+\frac{E\left(1-\nu \right)\left(e-e_0\right)}{\left(1+e_0\right)\left(1+\nu \right)\left(1-2\nu \right)}$

where $E$ is the compressive modulus and $\nu$ is Poisson’s ratio.

From the soil damage criterion, we can get the following formula:

(23)$P\le {\sigma }_{11}-{\sigma }_{sep}$

where ${\sigma }_{sep}$ is the separation strength of the freezing soil (negative value for tensile strength).

Substituting Equation (22) into Equation (23) and assuming that when $e$ is equal to $e_n$, Equation (23) takes the equal sign, the following relationships can be obtained:

(24)${{\displaystyle \int }}_x^h\left(\gamma -{\gamma }_0\right)d{x}^{*}+\frac{E\left(1-\nu \right)\left(e_n-e_0\right)}{\left(1+e_0\right)\left(1+\nu \right)\left(1-2\nu \right)}={\sigma }_{11}-{\sigma }_{sep}$

(25)$e_n=\frac{(P_{ob}+{{\displaystyle \int }}_x^h{\gamma }_{0}d{x}^{*}-{\sigma }_{sep})\left(1+\nu \right)\left(1-2\nu \right)\left(1+e_0\right)}{E\left(1-\nu \right)}+e_0$

Therefore, the Equation (22) can be expressed as:

(26)$P=\{\begin{array}{cc}{{\displaystyle \int }}_x^h\left(\gamma -{\gamma }_0\right)d{x}^{*}+\frac{E\left(1-\nu \right)\left(e-e_0\right)}{\left(1+e_0\right)\left(1+\nu \right)\left(1-2\nu \right)}& e<e_n\\ P_{ob}+{{\displaystyle \int }}_x^h\gamma d{x}^{*}-{\sigma }_{sep}& e\ge e_n\end{array}$

Frost heave can be calculated by integrating the strain in height and can be written as:

(27)$H_i={{\displaystyle \int }}_0^h{\epsilon }_{11}dx={{\displaystyle \int }}_0^h\frac{e-e_0}{1+e_0}dx$

2.3. Water Migration

The pore ice is assumed to be fixed during the freezing process, and the pore water is movable relative to the skeleton of soil. According to the law of mass conservation, the governing equation of seepage can be expressed as:

(28)$\frac{\partial \left(M_{a}dV\right)}{\partial t}+{\rho }_l\nabla \overrightarrow{v}dV=0$

where $M_a$ is the mass of pore water and ice in the unit volume of soils; $\overrightarrow{v}$ is the velocity relative to the soil particles:

(29)$M_a=n\left[S_i{\rho }_i+\left(1-S_i\right){\rho }_l\right]$

(30)$\overrightarrow{v}=-\frac{k}{{\gamma }_l}\nabla \left(\phi +{\gamma }_{l}z\right)$

where $k$ is the hydraulic conductivity; ${\gamma }_l$ is the weight of unit water; $\phi$ is the seepage pressure that consists of the pore pressure $P$ and the cryogenic suction $P_T$ (due to the ice/water interface tension) [29], as $\phi =\left(P+P_T\right)$:

(31)$\frac{k}{k_0}=\{\begin{array}{cc}{\left[1-\left(T-T_0\right)\right]}^{\beta }& T\le T_0\\ 1& T>T_0\end{array}$

where k₀ is the hydraulic conductivity when the soil is unfrozen, T₀ is the freezing point temperature (273.15 K) of water in pore ignoring the effect of the water vapor sorption isotherm [21,22,26,27], and $\beta$ is a parameter depending on the size and structure of pore.

During the soil freezing process, the experimental results showed that the temperature gradient will cause the pore water to migrate from the high-temperature region to the low-temperature region [9]. The Clapeyron equation is presented to describe the relationship between ice pressure, water pressure, and temperature [10]:

(32)$\frac{P_l}{{\rho }_l}+\frac{P_i}{{\rho }_i}=L\ln \frac{T}{T_0}$

where $P_l$ and $P_i$ are the water pressure and ice pressure, respectively.

The cryogenic suction $P_T$ can be expressed as [12]:

(33)$P_T=P_l-P_i\approx {\rho }_{l}L\frac{T-T_0}{T_0}$

Substituting Equations (29), (30), and (33) into Equation (28), the governing equation of water transfer can be written as:

$C_{21}\frac{\partial T}{\partial t}+C_{22}\frac{1}{1+e}\frac{\partial e}{\partial t}+{\rho }_l\nabla \left(-\frac{k}{{\gamma }_l}\nabla P\right)+{\rho }_l\nabla \left(-\frac{k{\rho }_{l}L}{{\gamma }_l{T}_0}\nabla T\right)+{\rho }_l\nabla \left(-k\nabla z\right)=0$

(34)$C_{21}=n\left({\rho }_i-{\rho }_l\right)\frac{\partial S_i}{\partial T}$

$C_{22}=S_i{\rho }_i+\left(1-S_i\right){\rho }_l.$

2.4. Ice Segregation

According to Miller’s second theory of frost heave, there is a frozen fringe in the process of soil freezing, and a new ice lens is formed in the frozen fringe. The development of the mechanical method is relatively complete because the mechanical method is relatively visualized and has a clear physical meaning. In 1961, Bishop gave an expression for pore pressure in unsaturated soils. Miller [30] applied the expression to saturated freezing soils as:

(35)$P=\chi P_l+\left(1-\chi \right)P_i$

where $\chi$ is a weighting coefficient, and can be expressed as the following Equation [31]:

(36)$\chi ={\left(1-S_i\right)}^{1.5}$

According to Equation (32), the pore water pressure can be obtained as:

(37)$P_l=\frac{{\rho }_l}{{\rho }_i}{P}_i+{\rho }_{l}L\ln \left(\frac{T}{T_0}\right)$

Substituting Equation (37) into Equation (35), the relationship between $P$ and $P_i$ can be obtained as:

(38)$P=\left[\frac{{\rho }_l}{{\rho }_i}\chi +\left(1-\chi \right)\right]P_i+\chi {\rho }_{l}L\ln \left(\frac{T}{T_0}\right)$

(39)$P\approx P_i+\chi {\rho }_{l}L\ln \frac{T}{T_0}$

The formation of ice lenses occurs when the ice pressure exceeds the sum of the external load and the soil separation strength:

(40)$P_i\ge {\sigma }_{11}-{\sigma }_{sep}$

Substituting Equation (39) into Equation (40), we have:

(41)$P-\chi {\rho }_{l}L\ln \frac{T}{T_0}\ge {\sigma }_{11}-{\sigma }_{sep}$

According to Equations (10), (22), and (25), the equivalent form of Equation (41) can be obtained

(42)$e\ge \frac{\left[{\left(1-S_i\right)}^{1.5}{\rho }_{l}L\ln \frac{T}{T_0}+P_{ob}+{{\displaystyle \int }}_x^h{\gamma }_{0}d{x}^{*}-{\sigma }_{sep}\right]\left(1+\nu \right)\left(1-2\nu \right)\left(1+e_0\right)}{E\left(1-\nu \right)}+e_0$

Equation (42) can be used as the criterion for the formation of ice lenses.

It can be seen from Equation (42) that the formation of ice lenses is related to the temperature, the overburden pressure, and the separation strength of the freezing soil.

3. Model Validation and Numerical Analysis

During the freezing process of soil, moisture will be replenished to the frozen fringe. In this section, the thermo-hydro-mechanical coupling model is used to describe the soil freezing process. After appropriate initial values and boundary conditions are given, the finite element method within the platform of COMSOL Multiphysics can be used to solve the equations to obtain the soil temperature field and displacement field during the freezing process.

3.1. Model Validation

Lai et al. [26] conducted one-side freezing experiments for silty clay column in an open system. The soil column was frozen from top to bottom and water was freely available at the base. The initial water content of the 20 cm soil column was 20.5%; the temperatures of the top and lower surfaces were controlled to 268.15 and 274.15 K; the initial temperature of the soil column was 274.15 K; and the temperature gradient of the soil column was 0.30 K/cm. The lower surface of the soil column can be used for water replenishment. Other parameters related to soil column from Lai et al. [26] are listed in Table 1. The simulation results and experimental data are compared in Figure 4. The frost heave and freezing depth of the soil column are shown in Figure 4a; the relationship between temperatures at the heights of 5, 10, 14, and 18 cm and time are shown in Figure 4b.

From Figure 4a, the frost heave of the soil column continues to increase with the freezing time. The frozen front descends very quickly at the beginning and then descends slowly. From Figure 4b, the temperature of the soil near the cold end descends quickly and then stabilizes. The simulation results of the frost heave, freezing depth, and temperature are almost consistent with the experimental data. As illustrated in Figure 4, the simulation results are consistent with the experimental results.

3.2. The Effect of Soil Parameters on the Soil Freezing Process

The types of soil are different in different regions or different depths. The porosity, thermal conductivity, and water conductivity of different soils are not identical. Previous scholars mainly studied the influence of external conditions on the freezing process of soils, such as the external temperature and overburden pressure [21,31]. There are few reports on the influence of soil properties on the freezing process. The effects of the initial void ratio ($e_0$), thermal conductivity of soil particles (${\lambda }_s$), and hydraulic conductivity ($k_0$) on the freezing process of soil were analyzed in this study. The parameters of the model required for calculation are listed in Table 2 and these parameters come from Yin et al. [27] and Zhou and Li [22]. In order to analyze the effect of soil properties on the freezing process of the soil column, $e_0$, $k_0$, and ${\lambda }_s$ assumed the different values listed in Table 3.

The relationships between temperatures at heights of 9.5, 9, 8, 7, 6, 4.5, and 2 cm and time in case 2 are demonstrated in Figure 5. The temperature drops rapidly near the cold end, which is mainly caused by a large temperature gradient. Near the temperature for 273.15 K, due to the release of latent heat by the phase transition of pore water, the decrease rate of temperature becomes smaller as illustrated in Figure 5. Zhou and Li [22] also obtained similar results.

The temperature distributions along the height of the soil column in case 2 are presented in Figure 6. The dotted line corresponds to 273.15 K. The temperature tends to be linear outside of the dotted line region. Near the temperature for 273.15 K, the change of the temperature gradient is mainly caused by the release of latent heat from the phase transition of pore water. The temperature fields at different times in case 2 are shown in Figure 7.

The effects of the initial void ratio ($e_0$), hydraulic conductivity ($k_0$), and thermal conductivity of soil particles (${\lambda }_s$) on frost heave and frozen front during the freezing process of soil are illustrated in Figure 8, Figure 9 and Figure 10, respectively. It can be seen from Figure 8 that the soil column with a high porosity ratio or a high water content has a more obvious frost heave and a smaller freezing depth at the same freezing time. A high initial void ratio means a high water content; therefore, more latent heat in the soil is released and the freezing depth of the soil column is also reduced. The final positions of the frozen front are basically unchanged.

The frost heave increases by 4 mm for every $0.5·{10}^{-11} \mathrm{m}/\mathrm{s}$ increase in hydraulic conductivity from $1·{10}^{-11} \mathrm{m}/\mathrm{s}$ to $2·{10}^{-11} \mathrm{m}/\mathrm{s}$ after freezing for 60 h, as presented in Figure 9. The seepage of soil is directly affected by the hydraulic conductivity, which causes more water to migrate from the unfrozen region to the frozen front. This leads to a greater release of latent heat. The freezing depth of the soil with high permeability conductivity is slightly smaller, but the change is not obvious. The final positions of the frozen front are basically the same.

The soil with large thermal conductivity causes a great amount of frost heave at the same freezing time. The heat transfer of soil is affected by the thermal conductivity of the soil particles. With the increase of the thermal conductivity, the depth of freezing increases as demonstrated in Figure 10. The freezing depth of the soil column is obviously affected by the change of thermal conductivity. The greater the thermal conductivity of the soil column, the greater the freezing depth. The final positions of the frozen front are basically unchanged. The freezing depth of the soil column with better thermal conductivity reaches a stable position faster.

As it is displayed in Figure 8, Figure 9 and Figure 10 that the frozen front quickly moved downward during the initial period of freezing, the freezing depth is about 5 cm after freezing for 30 h, and the final freezing depth remains about 6 cm. The final temperature is linearly distributed along the height of the soil column in Figure 6, and the finial position of the frozen front depends on the temperatures of the cold and warm ends. Through these seven cases, it can be found that frost heave increases almost linearly with freezing time.

The formation and distribution of ice lenses is important to study the freezing process of soils, the distribution of ice lenses of different cases after freezing 60 h is given in Figure 11, Figure 12 and Figure 13. The white portions represent the ice lenses, and the black portions represent the soil regions. Due to the final stabilization of the frozen front, a thick ice lens, which is called the ‘stable ice lens’ (region B in Figure 12), will be generated near the frozen front, and the ice lenses in region A is called as the ‘previous ice lenses’. Zhou and Li [22], Sheng et al. [32], and Lai et al. [26] obtained similar results about two types of ice lens.

After freezing for 60 h, the results of the cases under different initial void ratios are illustrated in Figure 11. The initial void ratio has little effect on the stable ice lens. This is because the final position of the frozen front will not be affected by the initial void ratio. With the increase of the initial void ratio, the number and thickness of previous ice lenses increase.

The formation of the ice lens is greatly affected by the hydraulic conductivity of soil, as presented in Figure 12. The number and thickness of previous ice lenses are increased in the case with large hydraulic conductivity, and the stable ice lens of soil becomes thicker. It is because the increase of the hydraulic conductivity causes the increase of the moisture of migration. However, there is no obvious change at the bottom of the stable ice lens. The reason is that the stable position of the frozen front is not affected by the hydraulic conductivity.

From Figure 13, the number of previous ice lenses is reduced, and stable ice lens has hardly changed with the increase in thermal conductivity. The large thermal conductivity of soil decreases the freezing time of frozen soil, which prevents further formation and increase of the ice lens in soil.

4. Conclusions and Discussion

In this study, a new thermo-hydro-mechanical coupling model describing the soil freezing process was established considering the formation of ice lenses. The model considered the heat transfer, heat convection, phase change of pore water, mechanical equilibrium, and formation of ice lenses. A new criterion for the formation of ice lenses was derived. Seven cases of soil freezing were simulated to study the effects of the initial void ratio, hydraulic conductivity, and thermal conductivity of soil particles on the freezing process of soils. The following conclusion is made:

The final stable position of the frozen front depends mainly on the temperature of the cold and warm ends, but the freezing front of the soil with a lower initial void ratio and hydraulic conductivity and a large thermal conductivity reaches the stable position faster, and there are fewer previous ice lenses in this case. There was a positive correlation between the frost heave and initial void ratio, hydraulic conductivity, and thermal conductivity. The stable ice lens is mainly affected by the hydraulic conductivity, and as it increases, the thickness of the stable ice lens becomes larger.

Because of the compressibility of soil particles, pore ice and pore water are ignored, and the applicability of the model will be reduced when used to study the freezing of deep soil. The research object of the theory proposed in this paper was saturated soil. For unsaturated soil, the degree of saturation needs to be introduced into the theory to adjust the model.

Author Contributions

D.L. Methodology, Investigation, Software, Writing—original draft. Y.Y. Conceptualization, Supervision, Writing—review and editing. C.C. Conceptualization, Writing—review and editing. Y.C. Data curation, Formal analysis. S.W. Resources, Formal analysis. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Fundamental Research Funds for the Central Universities, grant number 2019ZDPY18.

Acknowledgments

The authors are grateful to the editor and reviewers for providing constructive and positive comments for this paper.

Conflicts of Interest

The authors of this manuscript declare not to have any conflict of interest regarding this manuscript. None of the authors have any financial interests with any commercial entity that has interest in the subject or outcome of this manuscript including consultancy, paid expert consultancy, or honoraria, patent application, as well as other forms of conflict of interest, including personal and academic issues. The authors to the best of their knowledge conducted the study and reported the conclusions independently without any interference from partial or full funding sources or other entities.

Abbreviations

$e$	pore ratio
$n$	porosity
$P$	pore pressure, Pa
$t$	time, s
$T$	temperature, K
${\sigma }_{ij}, {\sigma }_{ij}^e$	total stress tensor and effective stress tensor, Pa
${\sigma }_{11}$	normal stress in vertical direction, Pa
${\epsilon }_{ij}$	strain tensor
${\epsilon }_{11}, {\epsilon }_{22}, {\epsilon }_{33}$	normal strain Constants and symbols
$c_s, c_l, c_i$	specific heat capacities of solid particles, pore water and pore ice, respectively, J·kg−1·K−1
$dV$	unit volume element of soils, m3
$e_0$	initial porosity ratio
$E$	compressive modulus, Pa
$f_i$	body force in the ith direction, Pa·m−1
$g$	gravitational acceleration, m·s−2
$h$	height of soil model, cm
$H$	volumetric heat capacity of the soils, J·K−1·m−3
$H_i$	frost heave, m
$\mathrm{K},{\mathrm{K}}_{\mathrm{S}}$	the bulk modulus of dry soil and soil particles, respectively, Pa
$k, k_0$	hydraulic conductivity of freezing soils and unfrozen soils, respectively, m·s−1
$L$	latent heat of fusion for water, J·kg−1
$M_a$	mass of pore water and ice in the unit volume of soils, kg·m−3
$P_i$, $P_l$	ice and water pressures, respectively, Pa
$P_T$	cryogenic suction due to the ice/water interface tension, Pa
$P_{ob}$	overburden pressure, Pa
$S_i$	volume percentage of ice in pores
$T_0$	freezing point temperature of water in pore, K
$T_r$	reference temperature, K
$\overrightarrow{v}$	seepage velocity of pore water, m·s−1
$\alpha$	parameter depending on the size of pore
${\alpha }^{*}$	Biot coefficient
$\beta$	parameter dependent on the size and structure of pore
$\mathsf{\Phi }$	heat content of soil per unit volume, J·m−3
$\lambda$	thermal conductivity of freezing soils, J·m−1·s−1·K−1
${\lambda }_s, {\lambda }_l, {\lambda }_i$	thermal conductivity of soil particles, water and ice, respectively, J·m−1·s−1·K−1
${\rho }_s, {\rho }_l, {\rho }_i$	density of soil particles, water and ice, respectively, kg·m−3
$\gamma$	unit weight of soil, kg·m−2·s−2
${\gamma }_l$	unit weight of water, kg·m−2·s−2
${\sigma }_{sep}$	separation strength of the freezing soil, Pa
${\sigma }_{11}^0$	initial normal stress in vertical direction, Pa
${\lambda }^{*}, {\mu }^{*}$	Lame constant
$\mathsf{\nu }$	Poisson’s ratio of soil
$\nabla$	Hamiltonian operator
$d$	differential operator
$\phi$	driving force that causes seepage
$\chi$	weighting coefficient
${\delta }_{ij}$	Kronecker symbol

Figures and Tables

Figure 1. Coupling relationship of physical fields in soil freezing.

Figure 2. Microscopic schematic of soil freezing.

Figure 3. Schematic diagram of the soil column mechanical model.

Figure 4. Comparison between simulation results and experimental data. (a: Frost heave profiles and frozen front profiles during freezing process; b: Temperature profiles at the heights of 5, 10, 14, and 18 cm during freezing process).

Figure 5. The temperature evolution curve of the soil column with different heights.

Figure 6. Temperature distribution at different times of case 2.

Figure 7. Temperature distribution at 0, 5, 10, and 35 h of case 2.

Figure 8. Freezing process curves of cases 1, 2, and 3.

Figure 9. Freezing process curves of cases 1, 4, and 5.

Figure 10. Freezing process curves of cases 5, 6, and 7.

Figure 11. Distribution of ice lenses at different initial void ratios. (Region A is called the ‘previous ice lenses’; Region B is called the ‘stable ice lens’.).

Figure 12. Distribution of ice lenses at different hydraulic conductivity.

Figure 13. Distribution of ice lenses at different thermal conductivity of solid particles.

Table 1

List of parameters used in the model validation.

Parameter	Value
Density of solid particles, ${\rho }_s$ (kg/m3)	2360
Density of water, ${\rho }_l$ (kg/m3)	1000
Density of ice, ${\rho }_i$ (kg/m3)	917
thermal conductivity of water, ${\lambda }_l$ (J·m−1·s−1·K−1)	0.58
thermal conductivity of ice, ${\lambda }_i$ (J·m−1·s−1·K−1)	2.22
thermal conductivity of solid particles, ${\lambda }_s$ (J·m−1·s−1·K−1)	1.5
Specific heat capacities of solid particles, $c_s$ (J·kg−1·K−1)	2360
Specific heat capacities of water, $c_l$ (J·kg−1·K−1)	4180
Specific heat capacities of ice, $c_i$ (J·kg−1·K−1)	1874
Gravitational acceleration, $g$ (m/s2)	9.81
Height of soil model, $h$ (cm)	20
Latent heat of fusion for water, $L$ (J/kg)	334,560

Table 2

List of parameters used in the calculation examples.

Parameter	Value
Density of solid particles, ${\rho }_s$ (kg/m3)	2700
Density of water, ${\rho }_l$ (kg/m3)	1000
Density of ice, ${\rho }_i$ (kg/m3)	917
Compressive modulus, E (MPa)	0.8
thermal conductivity of water, ${\lambda }_l$ (J·m−1·s−1·K−1)	0.58
thermal conductivity of ice, ${\lambda }_i$ (J·m−1·s−1·K−1)	2.22
Specific heat capacities of solid particles, $c_s$ (J·kg−1·K−1)	2360
Specific heat capacities of water, $c_l$ (J·kg−1·K−1)	4180
Specific heat capacities of ice, $c_i$ (J·kg−1·K−1)	1874
Gravitational acceleration, $g$ (m/s2)	9.8
Height of soil model, $h$ (cm)	10
Latent heat of fusion for water, $L$ (J/kg)	334,000
Overburden pressure, $P_{ob}$ (kPa)	200

Table 3

Values of the variables in the simulated cases.

Case No.	Initial Void Ratio	Initial Hydraulic Conductivity(10−11m s−1)	Temperatures of Top and Lower Surfaces (K)	Thermal Conductivity of Soil Particles(W·m−1 K−1)
1	0.4	1.0	271.65/274.15	0.70
2	0.6	1.0	271.65/274.15	0.70
3	0.8	1.0	271.65/274.15	0.70
4	0.4	1.5	271.65/274.15	0.70
5	0.4	2.0	271.65/274.15	0.70
6	0.4	2.0	271.65/274.15	1.20
7	0.4	2.0	271.65/274.15	1.70