# A Review of the Research on Thermo-Hydro-Mechanical Coupling for the Frozen Soil

Authors: Yun Chao Teng; Zhen Chao Teng; Yu Liu; Xiao Yan Liu; Ya Dong Zhou; Jia Lin Liu; Bo Li
Journal: Geofluids
 (2022)

---

## Abstract

This paper reviews the history of the research development on the coupling mechanism of the multiphysical field, e.g., thermo-hydro-mechanical (THM), for frozen soil. The objective is to deepen the current understanding of the theories and mechanism of multiphysical field coupling in the frozen soil and the dynamic changes in the temperature, moisture, and stress fields during soil freezing. A new differential equation of the coupling of temperature field and moisture field is proposed. Based on the DiscreteFrechetDist algorithm, a fitting method of evaluating a curve is proposed. The paper is expected to help understand the soil freezing process in cold regions and enhance the innovativeness of the research methodologies dealing with multifield coupling for the frozen soil.


---

## Body

## 1. Introduction

Frost heave and thaw settlement are the cause of damage to infrastructure in the seasonal frozen region. The occurrence of frost heave and thaw settlement is primarily due to the changes in temperature, moisture, stress, and concentration fields of the frozen ground and interactions between multiphysical fields [1–4]. The soil stability in the seasonally frozen ground depends on the interactions and mutual influence between heat transfer, moisture migration, and phase change during the freeze-thaw cycles. Therefore, it is necessary to understand the mechanisms behind the multiphysical field coupling of the frozen soil.The coupling process of the multiphysical field mainly comprises thermo-hydro-mechanical (THM) coupling, thermo-hydro-vapor-mechanical (THVM) coupling [5–9], and thermo-hydro-salt-mechanical (THSM) coupling [10, 11]. Bai et al. [12] built the THVM coupling model and verified the model experimentally. It was found that when the initial water content was below the critical value, the water vapor migration was significant, and the frost heave took place at a slow rate. When the initial water content was above the critical value, liquid water migration dominated. The frost heave took place at an accelerating rate. The model lays the foundation for studying the three-dimensional THVM coupling model for unsaturated soil. Zhang et al. [10] proposed the mathematical model describing the THSM multiple fields that satisfied the solute conservation equation. However, this model did not consider the effects of salt on frozen soil deformation. The crystallization dynamic under phase change was revealed. The THSM coupling model was built in COMSOL for the saturated frozen sulfate soil. The phase change pore pressure drove water-salt migration and soil deformation.The soil freezing process involves the interaction of several physical fields. The moisture, temperature, stress, and salt concentration fields interact with each other and mutually influence each other, and these effects are coupled and combined. Some experts attempt to understand the multiphysical areas on three different levels through frozen soil simulation: multifields, multiregional, and multiscales. Specifically, the study of the response and excitation of multiphysical fields covers the frozen soil system, the direct interactions between frozen soil continua with different features via the boundaries, and the continuous variation from microscopic to macroscopic behaviors of the frozen soil system under different scales.A thorough study on the multiphysical field coupling for the frozen soil is highly significant for the safety and stability of lifeline engineering in cold regions, including building construction, railway operation, and oil pipelines. Researchers have conducted extensive 3D simulations and experiments on the multiphysical field coupling for the frozen soil but made only limited achievements given the broad range of engineering sectors involved in the multiphysical field coupling for the frozen soil. Therefore, it was crucial to probe deeper into the simulation and experiments of the multiphysical field coupled systems for the frozen soil.Based on the existing applied cases of the THM coupling theories, mathematical models, and finite element software for the frozen soil to solve the problems of THM coupling, we summarize the current research status and conceptualize the future direction of research on THM coupling for the frozen soil. This paper reviews the existing studies on the differential equations of multiphysical field coupling and experimental simulation of multiphysical field coupling for the frozen soil. This summary not only provides some reference for the future research work but also clarifies the research direction and focus. Frozen soil coupling mode and characteristics are given in Table1.Table 1
Intensity of multiphysical field coupling of the frozen soil and the features of coupling of different intensities.
The intensity of multiphysical field couplingCoupling methodCoupling featuresStrong couplingElement matrix methodIt is difficult to control the calculation accuracy. Nevertheless, this method has advantages when used for theoretical analysis of the coupled field of frozen soilWeak couplingSuperposition methodThe calculation is highly complex and involves many time steps. Therefore, this approach is more suitable for the numerical analysis of multiphysical fields for frozen soilNote: element matrix method: use an element matrix (or load vector) to fit the coupling to governing equations entirely. Then, the system of governing equations is solved directly. Superposition method: the calculation results of the first physical field are imposed onto the second and third physical fields as an external load to realize coupling between the three fields.


## 2. Research on the Temperature, Moisture, and Stress Field Equations for the Frozen Soil

As already known, the physical fields of the frozen soil can be solved using a system of partial differential equations. Considering the complexity of differential equations and partial differential equations, the coupling process is divided into the following types based on the features of the physical model, as shown in Figure1.Figure 1
Coupling classification.The driving force due to temperature changes during soil freezing leads to moisture migration in the soil. Meanwhile, the thermal expansion during soil freezing also generates temperature stress. The latent heat of ice-water phase change poses an inhibitory effect on cold energy transfer. Besides, moisture migration and the formation of ice lens in the soil can alter the original pore size, and the consequent volumetric strain further causes the frost heave of the soil. Given the above, the frost heave itself is the result of multiphysical field coupling behavior. The coupling theory of the multiphysical field remains one of the most extensively studied theories related to frozen soil at home and abroad. Considering the complexity of multiphysical field coupling, it is recommended that the temperature, moisture, and stress fields should be analyzed separately at the preliminary stage.


### 2.1. Research on the Mathematical Model of the Frozen Soil Temperature Field

In the differential equation of the temperature field, the heat transfer parameters relate to atmospheric factors, including evaporation, radiation, wind velocity, and humidity. Therefore, Wang et al. [13] proposed equation (1) to describe the phase-change nonsteady temperature field considering the above factors.
(1)K+NDtTt=Pt+NDtTt−Dt,(2)k1⋅n=∫−11∫−11λ16JD2∂Hn∂ξ+D1∂Hn∂ηdξdη+Aαsl,n=i,j,k,m,(3)n1⋅n=∫−11∫−11ρCpJHlHndξdηl,n=i,j,k,m,whereP is the synthetic array and the subscript t indicated the time [14]. According to Wang et al. [13], λ=4.17w2+1504100.25ρd−3.9; ρCp is the volumetric heat capacity and specific heat capacity [15] C=ρd1.27+0.021w103. A=0 for the non-Robin boundary, and
(4)A=1/3,l=n,1/6,else,for the Robin boundary.Some researchers considered THM coupling for the frozen soil but did not couple the heat release due to crystallization-related phase change of salts into the model. Zhang et al. [16] proposed a two-dimensional differential equation of THS coupling for the first time for salinized frozen soil. Based on the conservation of energy and Fourier’s law, the equation took the exothermic effect of the crystallization-related phase change of Glauber’s salt and the latent heat of ice-water phase change as the heat sources. Then, the differential equation (5) of heat transfer in the salinized soil was derived.
(5)ρCθ∂T∂t=λθ∇2T+Lcρc∂θc∂t+Lmρi∂θi∂t.In seasonally frozen ground, the temperature field around the buried oil pipeline is influenced not only by the thermophysical property but also by the soil temperature within the pipeline, freeze-thaw cycle, and isothermal layer. Teng et al. [17] established the transient state heat conduction:
(6)Cρs∂T∂t=divλgradT⟶−Lρw∂θw∂y+βtT∂εv∂t.Based on the above research, a new coupling equation is proposed. In the seasonal freezing zone, the influence of steam on the latent heat of phase change, the latent heat of salt crystallization, and the volume strain effect under temperature changes should be considered at the same time, and the dimensions of the differential terms are the same.(7)ρCθ∂T∂t=λθ∇2T+Liρi∂θi∂t+Lvρv∂θv∂t+Lcρc∂θc∂t+βtT∂εν∂t,whereρ, ρc, ρi, and ρv are the density of soil, salt, ice, and vapor, respectively (kg/m3); θi, θv, and θc are the volumetric content of ice, vapor, and salt, respectively (cm3/cm3); C is the specific heat capacity of frozen soil (J/(kg·°C)); Lc, Li, and Lv are the latent heat of salt crystallization, water freezing, and water evaporation, respectively (J/g); T is the temperature (°C); βt is the thermal stress coefficient of the frozen soil (MPa/k); εv is the volumetric strain (1); and λ is the thermal conductivity (W/(cm·°C)).The temperature field theory for frozen soil develops from conservation of energy, Fourier’s law, and the law of thermodynamics. The temperature field theory was first proposed for molten soil and later extended to frozen soil. While the TH coupling theory has been intensively studied at home and abroad, the THM coupling receives less attention. For seasonally frozen ground, the temperature gradients in the temperature field are large during one freeze-thaw cycle of soil. Therefore, the differential equation for the temperature field needs to incorporate the differential term of the volumetric strain effect to consider the temperature changes in frozen saline soil, the joint action of exothermic effect caused by salt crystallization phase change (e.g., Glauber’s salt), and the latent heat from ice-water phase change. In addition, the differential equation of the temperature field in the unsaturated frozen soil should consider the influence of ice-water phase change and the impact of water vapor on the migration and freezing of liquid water. Most past studies focus on the two-dimensional models of the temperature field of the frozen soil, but rarely on the three-dimensional models.


### 2.2. Research on the Mathematical Model of the Frozen Soil Moisture Field

The unfrozen water in the salinized soil exists below 0°C, and the pore water migration satisfies Darcy’s law [18]. Based on the Richards equation, Zhang et al. [19] considered the influence of frozen water in the pores on the migration of unfrozen water and the effect of latent heat of ice-water phase change and phase change of Glauber’s salt in reducing the content of unfrozen water in the soil. Hence, the following equation was proposed [16]:The differential equation of the liquid water migration in the roadbed made of unsaturated saline soil takes the following form:(8)∂θu∂t+ρiρw∂θi∂t+ρcM180ρwMmx∂θc∂t=∇Dθu∇θu+kgθu,wherekgθu=∂k/∂z is the permeability coefficient of the unsaturated soil in the gravitational direction and Dθu=kθu/cθu⋅I is the diffusion rate of water. I is the impedance factor, representing the impeding effect of ice on liquid water migration I=1010θi; M180 is the molar mass of 10 water molecules (180 g/mol); Mmx is the molar molecular mass of Glauber’s salt (322 g/mol). Thus, the moisture content of the saline water is composed of three parts, ice, liquid water, and crystal water of Glauber’s salt. Therefore, the volume content of water in the soil is θ=θu+ρI/ρwθi+θcρcM180/ρwMmx.The differential equation of the water field for unsaturated frozen water should consider the influence of ice-water phase change and the impact of water vapor on the migration and freezing of liquid water. Zhang et al., Teng et al., and Zhang et al. [20–22] demonstrated that the model proposed by Philip and de Vries [23] was able to simulate the water-heat-water vapor behavior of unsaturated soil. Bai et al. [24] established the following equation for the flow of unfrozen water and water vapor based on the modified Darcy’s law and the equation of mass conservation:
(9)∂θl∂t+∂θv∂t+ρiρw∂θi∂t=∇⋅Klh∇h+z+KlT∇T+∇⋅Kνh∇h+KνT∇T.On the basis of the above research on soil moisture field [16, 24], a new coupling equation of moisture field is proposed. In seasonally frozen areas with high-salt soils, pore steam flow, the effect of ice on the transport of unfrozen water, the influence of the latent heat of ice-water phase change and concentration changes on the reduction of soil unfrozen water content should be considered at the same time. The dimension of each differential term is the same.
(10)∂θl∂t+∂θv∂t+ρiρw∂θi∂t+ρcM180ρwMmx∂θc∂t=∇⋅Klh∇h+z+KlT∇T+∇⋅Kνh∇h+KνT∇T,whereθl, θi, θv, and θc are the volumetric content of soil, ice, vapor, and salt, respectively (cm3/cm3); ρw, ρi, and ρc are the densities of liquid water, ice, and salt, respectively (kg/m3); Klh/Kvh and KlT/KvT are the isothermal and thermal hydraulic/vapor conductivities for liquid phase fluxes due to gradients in h and T, respectively; z is the spatial coordinate; M180 is the molar mass of 10 water molecules; and Mmx is the molar molecular mass of Glauber’s salt.The differential equation of the water field for unsaturated frozen soil should consider TH coupling and the influence of phase change due to the crystallization of Glauber’s salt on the unfrozen water content for the salinized soil. Apart from this, it is also necessary to consider the influence of water vapor on the migration and freezing of liquid water.


### 2.3. Research on the Mathematical Model of the Stress Field

#### 2.3.1. Discriminant Method for Formation of Ice Lens

The formation and development of ice lens are shown in Table2.Table 2
The criterion of ice lens formation.
CriterionFormulaSeparation pressure:icepressure>separationpressure [25]Psep=Pob+2σiw/RfPR (∗)Pore ice pressure:icepressure=totalexternalload [26]PiT=ρiPwT/ρw−LT/T0 (∗∗)Pore pressure:icepressure>sumofexternalload+soilseparationstrength [12]e⩾1−Si1.5ρlLlnT/T0+Pob+∫xhγ0dx∗−σsep1+ν1−2v1+e0/E1−ν+e0 (∗∗∗)Premelting film theory: the molecular interactions between ice generates net thermal molecular force and soil particles [27, 28]FT=−∫TPTdΓ=∫Tγwiκ−ρiLTm−T/TmdΓ (∗∗∗∗)The microscopic criterion of ice lens formation:averagestress>tensilestressoffrozensoil [1]Pri⩾σt (∗∗∗∗∗)Note: equations (∗,∗∗,∗∗∗,∗∗∗∗,∗∗∗∗∗) provided the judgment criteria for the formation of ice lenses. Ice lens formation is not only related to temperature and overburden pressure but also to the separation strength of the frozen soil.


#### 2.3.2. Research on Constitutive Relation

The frozen soil undergoes a time-varying deformation under the action of load. According to the seepage test of the frozen soil [29–32], the Bingham model could well describe the rheological properties of the frozen soil. The constitutive relation for a time step per unit time is from the equation below [1, 30, 33]:
(11)Δσ=DTΔε−Δεvp+Δεfh,whereDT is the elastic modulus array:
(12)DT=ET1−vT1+vT1−2vTETvT1+vT1−vT0ETvT1+vT1−v2TETvT1+vT1−v2TET1−vT1+vT1−v2T0ETvT1+vT1−v2T00ET21+vT0ETvT1+vT1−v2TETvT1+vT1−v2T0ET1−vT1+vT1−2vT,Er is the elastic modulus relative to soil temperature; νT is the Poisson’s ratio relative to soil temperature; Δεvp is the viscoplastic strain:ε̇vp=γTΦF∂Q/∂σ; Δεfh is the strain caused by frost heave; γT is the viscosity parameter relative to soil temperature; Φ is any function: ΦF=F−F0/F0,
(13)ΦF=ΦF,ifF>F0,ΦF=0,ifF⩽F0,whereQ is the plastic potential in the law of flow and Q is equivalent to the yield function of F [31, 34]. Since adopting the D-P yield criterion in this literature [31, 34], there is ∂Q/∂σ=∂F/∂σ=∂F/∂σr∂F/∂σz∂F/∂τrz∂F/∂σ0T; as the soil freezes, the resulting frost heave is caused by in situ transformation with moisture migration and freezing. Thus, the increase in soil volume is given by Liu et al. [1], Li et al. [35], Lai et al. [31], and Zienkiewicz and Taylor [34]. ΔεfhV=θit+Δt+θwt+Δt−nst, where nst is the volume porosity. If ΔεfhV<0, then ΔεfhV=0.THM coupling behavior is a complex and highly nonlinear problem. Due to the complexity of analytical solution, Li et al. [36] developed MHMIP-2D based on the FORTRAN program, which found practical applications [29, 30, 37].In 2020, Zhang et al. [38] considered the influence of soil gravity and expansion caused by ice-water phase change on soil deformation. The governing equation of the stress-strain relation of the frozen soil becomes
(14)∇Dε−εth+Fν=0,whereD is the elasticity matrix:
(15)D=E1+ν1−2v1−ννν000ν1−νν000νν1−ν0000001−2vv0000001−2vv0000001−2vv,E is the elastic modulus (Pa), ν is Poisson’s ratio, ε is the strain matrix, and εth is the thermal strain matrix. For the thermal expansion of ordinary material, a linear thermal expansion coefficient represents the thermal strain of frozen soil caused by frost heave related to ice-water phase change.εth=εxνεyνεzνγxνγyνγzνT, where εth is the total strain matrix; εvx, εvy, and εvz are the typical strain components of thermal expansion in x, y, and z directions, respectively. For frozen isotopic soil, εxν=εyν=εzν=−αT−Tref, where α is the coefficient of linear thermal expansion (l/°C), which is related to the mass content and temperature; Tref is the reference temperature (°C); and γvx, γvy, and γvz are the components of shear stress induced by thermal expansion in the x, y, and z directions, respectively. For isotropic frozen soil, γvx=γvy=γvz=0, Fv is a volume force matrix.


### 2.4. Numerical Simulation-Theoretical Solution-Test Solution Method and Process

#### 2.4.1. Simulation-Test Solution Process

(1)
Mathematical equations for the temperature field, moisture field, and stress field of frozen soil(2)
Input simulation frozen soil software, such as COMSOL, ANSYS, and ABAQUS (using MATLAB language, ABAQUS as solver development THM multifield coupling program, and analysis of frozen soil multifield coupling effect)(3)
Obtain a regular model consistent with the experiment


#### 2.4.2. Theoretical Solution-Experimental Solution Process

(1)
Mathematical equations for temperature field, moisture field, and stress field of frozen soil(2)
When unknown variables are greater than the number of equations, supplementary equations need to be added(3)
Software MATLAB calculates the theoretical value(4)
Obtain the law theory consistent with the experiment


## 2.1. Research on the Mathematical Model of the Frozen Soil Temperature Field

In the differential equation of the temperature field, the heat transfer parameters relate to atmospheric factors, including evaporation, radiation, wind velocity, and humidity. Therefore, Wang et al. [13] proposed equation (1) to describe the phase-change nonsteady temperature field considering the above factors.
(1)K+NDtTt=Pt+NDtTt−Dt,(2)k1⋅n=∫−11∫−11λ16JD2∂Hn∂ξ+D1∂Hn∂ηdξdη+Aαsl,n=i,j,k,m,(3)n1⋅n=∫−11∫−11ρCpJHlHndξdηl,n=i,j,k,m,whereP is the synthetic array and the subscript t indicated the time [14]. According to Wang et al. [13], λ=4.17w2+1504100.25ρd−3.9; ρCp is the volumetric heat capacity and specific heat capacity [15] C=ρd1.27+0.021w103. A=0 for the non-Robin boundary, and
(4)A=1/3,l=n,1/6,else,for the Robin boundary.Some researchers considered THM coupling for the frozen soil but did not couple the heat release due to crystallization-related phase change of salts into the model. Zhang et al. [16] proposed a two-dimensional differential equation of THS coupling for the first time for salinized frozen soil. Based on the conservation of energy and Fourier’s law, the equation took the exothermic effect of the crystallization-related phase change of Glauber’s salt and the latent heat of ice-water phase change as the heat sources. Then, the differential equation (5) of heat transfer in the salinized soil was derived.
(5)ρCθ∂T∂t=λθ∇2T+Lcρc∂θc∂t+Lmρi∂θi∂t.In seasonally frozen ground, the temperature field around the buried oil pipeline is influenced not only by the thermophysical property but also by the soil temperature within the pipeline, freeze-thaw cycle, and isothermal layer. Teng et al. [17] established the transient state heat conduction:
(6)Cρs∂T∂t=divλgradT⟶−Lρw∂θw∂y+βtT∂εv∂t.Based on the above research, a new coupling equation is proposed. In the seasonal freezing zone, the influence of steam on the latent heat of phase change, the latent heat of salt crystallization, and the volume strain effect under temperature changes should be considered at the same time, and the dimensions of the differential terms are the same.(7)ρCθ∂T∂t=λθ∇2T+Liρi∂θi∂t+Lvρv∂θv∂t+Lcρc∂θc∂t+βtT∂εν∂t,whereρ, ρc, ρi, and ρv are the density of soil, salt, ice, and vapor, respectively (kg/m3); θi, θv, and θc are the volumetric content of ice, vapor, and salt, respectively (cm3/cm3); C is the specific heat capacity of frozen soil (J/(kg·°C)); Lc, Li, and Lv are the latent heat of salt crystallization, water freezing, and water evaporation, respectively (J/g); T is the temperature (°C); βt is the thermal stress coefficient of the frozen soil (MPa/k); εv is the volumetric strain (1); and λ is the thermal conductivity (W/(cm·°C)).The temperature field theory for frozen soil develops from conservation of energy, Fourier’s law, and the law of thermodynamics. The temperature field theory was first proposed for molten soil and later extended to frozen soil. While the TH coupling theory has been intensively studied at home and abroad, the THM coupling receives less attention. For seasonally frozen ground, the temperature gradients in the temperature field are large during one freeze-thaw cycle of soil. Therefore, the differential equation for the temperature field needs to incorporate the differential term of the volumetric strain effect to consider the temperature changes in frozen saline soil, the joint action of exothermic effect caused by salt crystallization phase change (e.g., Glauber’s salt), and the latent heat from ice-water phase change. In addition, the differential equation of the temperature field in the unsaturated frozen soil should consider the influence of ice-water phase change and the impact of water vapor on the migration and freezing of liquid water. Most past studies focus on the two-dimensional models of the temperature field of the frozen soil, but rarely on the three-dimensional models.


## 2.2. Research on the Mathematical Model of the Frozen Soil Moisture Field

The unfrozen water in the salinized soil exists below 0°C, and the pore water migration satisfies Darcy’s law [18]. Based on the Richards equation, Zhang et al. [19] considered the influence of frozen water in the pores on the migration of unfrozen water and the effect of latent heat of ice-water phase change and phase change of Glauber’s salt in reducing the content of unfrozen water in the soil. Hence, the following equation was proposed [16]:The differential equation of the liquid water migration in the roadbed made of unsaturated saline soil takes the following form:(8)∂θu∂t+ρiρw∂θi∂t+ρcM180ρwMmx∂θc∂t=∇Dθu∇θu+kgθu,wherekgθu=∂k/∂z is the permeability coefficient of the unsaturated soil in the gravitational direction and Dθu=kθu/cθu⋅I is the diffusion rate of water. I is the impedance factor, representing the impeding effect of ice on liquid water migration I=1010θi; M180 is the molar mass of 10 water molecules (180 g/mol); Mmx is the molar molecular mass of Glauber’s salt (322 g/mol). Thus, the moisture content of the saline water is composed of three parts, ice, liquid water, and crystal water of Glauber’s salt. Therefore, the volume content of water in the soil is θ=θu+ρI/ρwθi+θcρcM180/ρwMmx.The differential equation of the water field for unsaturated frozen water should consider the influence of ice-water phase change and the impact of water vapor on the migration and freezing of liquid water. Zhang et al., Teng et al., and Zhang et al. [20–22] demonstrated that the model proposed by Philip and de Vries [23] was able to simulate the water-heat-water vapor behavior of unsaturated soil. Bai et al. [24] established the following equation for the flow of unfrozen water and water vapor based on the modified Darcy’s law and the equation of mass conservation:
(9)∂θl∂t+∂θv∂t+ρiρw∂θi∂t=∇⋅Klh∇h+z+KlT∇T+∇⋅Kνh∇h+KνT∇T.On the basis of the above research on soil moisture field [16, 24], a new coupling equation of moisture field is proposed. In seasonally frozen areas with high-salt soils, pore steam flow, the effect of ice on the transport of unfrozen water, the influence of the latent heat of ice-water phase change and concentration changes on the reduction of soil unfrozen water content should be considered at the same time. The dimension of each differential term is the same.
(10)∂θl∂t+∂θv∂t+ρiρw∂θi∂t+ρcM180ρwMmx∂θc∂t=∇⋅Klh∇h+z+KlT∇T+∇⋅Kνh∇h+KνT∇T,whereθl, θi, θv, and θc are the volumetric content of soil, ice, vapor, and salt, respectively (cm3/cm3); ρw, ρi, and ρc are the densities of liquid water, ice, and salt, respectively (kg/m3); Klh/Kvh and KlT/KvT are the isothermal and thermal hydraulic/vapor conductivities for liquid phase fluxes due to gradients in h and T, respectively; z is the spatial coordinate; M180 is the molar mass of 10 water molecules; and Mmx is the molar molecular mass of Glauber’s salt.The differential equation of the water field for unsaturated frozen soil should consider TH coupling and the influence of phase change due to the crystallization of Glauber’s salt on the unfrozen water content for the salinized soil. Apart from this, it is also necessary to consider the influence of water vapor on the migration and freezing of liquid water.


## 2.3. Research on the Mathematical Model of the Stress Field

### 2.3.1. Discriminant Method for Formation of Ice Lens

The formation and development of ice lens are shown in Table2.Table 2
The criterion of ice lens formation.
CriterionFormulaSeparation pressure:icepressure>separationpressure [25]Psep=Pob+2σiw/RfPR (∗)Pore ice pressure:icepressure=totalexternalload [26]PiT=ρiPwT/ρw−LT/T0 (∗∗)Pore pressure:icepressure>sumofexternalload+soilseparationstrength [12]e⩾1−Si1.5ρlLlnT/T0+Pob+∫xhγ0dx∗−σsep1+ν1−2v1+e0/E1−ν+e0 (∗∗∗)Premelting film theory: the molecular interactions between ice generates net thermal molecular force and soil particles [27, 28]FT=−∫TPTdΓ=∫Tγwiκ−ρiLTm−T/TmdΓ (∗∗∗∗)The microscopic criterion of ice lens formation:averagestress>tensilestressoffrozensoil [1]Pri⩾σt (∗∗∗∗∗)Note: equations (∗,∗∗,∗∗∗,∗∗∗∗,∗∗∗∗∗) provided the judgment criteria for the formation of ice lenses. Ice lens formation is not only related to temperature and overburden pressure but also to the separation strength of the frozen soil.


### 2.3.2. Research on Constitutive Relation

The frozen soil undergoes a time-varying deformation under the action of load. According to the seepage test of the frozen soil [29–32], the Bingham model could well describe the rheological properties of the frozen soil. The constitutive relation for a time step per unit time is from the equation below [1, 30, 33]:
(11)Δσ=DTΔε−Δεvp+Δεfh,whereDT is the elastic modulus array:
(12)DT=ET1−vT1+vT1−2vTETvT1+vT1−vT0ETvT1+vT1−v2TETvT1+vT1−v2TET1−vT1+vT1−v2T0ETvT1+vT1−v2T00ET21+vT0ETvT1+vT1−v2TETvT1+vT1−v2T0ET1−vT1+vT1−2vT,Er is the elastic modulus relative to soil temperature; νT is the Poisson’s ratio relative to soil temperature; Δεvp is the viscoplastic strain:ε̇vp=γTΦF∂Q/∂σ; Δεfh is the strain caused by frost heave; γT is the viscosity parameter relative to soil temperature; Φ is any function: ΦF=F−F0/F0,
(13)ΦF=ΦF,ifF>F0,ΦF=0,ifF⩽F0,whereQ is the plastic potential in the law of flow and Q is equivalent to the yield function of F [31, 34]. Since adopting the D-P yield criterion in this literature [31, 34], there is ∂Q/∂σ=∂F/∂σ=∂F/∂σr∂F/∂σz∂F/∂τrz∂F/∂σ0T; as the soil freezes, the resulting frost heave is caused by in situ transformation with moisture migration and freezing. Thus, the increase in soil volume is given by Liu et al. [1], Li et al. [35], Lai et al. [31], and Zienkiewicz and Taylor [34]. ΔεfhV=θit+Δt+θwt+Δt−nst, where nst is the volume porosity. If ΔεfhV<0, then ΔεfhV=0.THM coupling behavior is a complex and highly nonlinear problem. Due to the complexity of analytical solution, Li et al. [36] developed MHMIP-2D based on the FORTRAN program, which found practical applications [29, 30, 37].In 2020, Zhang et al. [38] considered the influence of soil gravity and expansion caused by ice-water phase change on soil deformation. The governing equation of the stress-strain relation of the frozen soil becomes
(14)∇Dε−εth+Fν=0,whereD is the elasticity matrix:
(15)D=E1+ν1−2v1−ννν000ν1−νν000νν1−ν0000001−2vv0000001−2vv0000001−2vv,E is the elastic modulus (Pa), ν is Poisson’s ratio, ε is the strain matrix, and εth is the thermal strain matrix. For the thermal expansion of ordinary material, a linear thermal expansion coefficient represents the thermal strain of frozen soil caused by frost heave related to ice-water phase change.εth=εxνεyνεzνγxνγyνγzνT, where εth is the total strain matrix; εvx, εvy, and εvz are the typical strain components of thermal expansion in x, y, and z directions, respectively. For frozen isotopic soil, εxν=εyν=εzν=−αT−Tref, where α is the coefficient of linear thermal expansion (l/°C), which is related to the mass content and temperature; Tref is the reference temperature (°C); and γvx, γvy, and γvz are the components of shear stress induced by thermal expansion in the x, y, and z directions, respectively. For isotropic frozen soil, γvx=γvy=γvz=0, Fv is a volume force matrix.


## 2.3.1. Discriminant Method for Formation of Ice Lens

The formation and development of ice lens are shown in Table2.Table 2
The criterion of ice lens formation.
CriterionFormulaSeparation pressure:icepressure>separationpressure [25]Psep=Pob+2σiw/RfPR (∗)Pore ice pressure:icepressure=totalexternalload [26]PiT=ρiPwT/ρw−LT/T0 (∗∗)Pore pressure:icepressure>sumofexternalload+soilseparationstrength [12]e⩾1−Si1.5ρlLlnT/T0+Pob+∫xhγ0dx∗−σsep1+ν1−2v1+e0/E1−ν+e0 (∗∗∗)Premelting film theory: the molecular interactions between ice generates net thermal molecular force and soil particles [27, 28]FT=−∫TPTdΓ=∫Tγwiκ−ρiLTm−T/TmdΓ (∗∗∗∗)The microscopic criterion of ice lens formation:averagestress>tensilestressoffrozensoil [1]Pri⩾σt (∗∗∗∗∗)Note: equations (∗,∗∗,∗∗∗,∗∗∗∗,∗∗∗∗∗) provided the judgment criteria for the formation of ice lenses. Ice lens formation is not only related to temperature and overburden pressure but also to the separation strength of the frozen soil.


## 2.3.2. Research on Constitutive Relation

The frozen soil undergoes a time-varying deformation under the action of load. According to the seepage test of the frozen soil [29–32], the Bingham model could well describe the rheological properties of the frozen soil. The constitutive relation for a time step per unit time is from the equation below [1, 30, 33]:
(11)Δσ=DTΔε−Δεvp+Δεfh,whereDT is the elastic modulus array:
(12)DT=ET1−vT1+vT1−2vTETvT1+vT1−vT0ETvT1+vT1−v2TETvT1+vT1−v2TET1−vT1+vT1−v2T0ETvT1+vT1−v2T00ET21+vT0ETvT1+vT1−v2TETvT1+vT1−v2T0ET1−vT1+vT1−2vT,Er is the elastic modulus relative to soil temperature; νT is the Poisson’s ratio relative to soil temperature; Δεvp is the viscoplastic strain:ε̇vp=γTΦF∂Q/∂σ; Δεfh is the strain caused by frost heave; γT is the viscosity parameter relative to soil temperature; Φ is any function: ΦF=F−F0/F0,
(13)ΦF=ΦF,ifF>F0,ΦF=0,ifF⩽F0,whereQ is the plastic potential in the law of flow and Q is equivalent to the yield function of F [31, 34]. Since adopting the D-P yield criterion in this literature [31, 34], there is ∂Q/∂σ=∂F/∂σ=∂F/∂σr∂F/∂σz∂F/∂τrz∂F/∂σ0T; as the soil freezes, the resulting frost heave is caused by in situ transformation with moisture migration and freezing. Thus, the increase in soil volume is given by Liu et al. [1], Li et al. [35], Lai et al. [31], and Zienkiewicz and Taylor [34]. ΔεfhV=θit+Δt+θwt+Δt−nst, where nst is the volume porosity. If ΔεfhV<0, then ΔεfhV=0.THM coupling behavior is a complex and highly nonlinear problem. Due to the complexity of analytical solution, Li et al. [36] developed MHMIP-2D based on the FORTRAN program, which found practical applications [29, 30, 37].In 2020, Zhang et al. [38] considered the influence of soil gravity and expansion caused by ice-water phase change on soil deformation. The governing equation of the stress-strain relation of the frozen soil becomes
(14)∇Dε−εth+Fν=0,whereD is the elasticity matrix:
(15)D=E1+ν1−2v1−ννν000ν1−νν000νν1−ν0000001−2vv0000001−2vv0000001−2vv,E is the elastic modulus (Pa), ν is Poisson’s ratio, ε is the strain matrix, and εth is the thermal strain matrix. For the thermal expansion of ordinary material, a linear thermal expansion coefficient represents the thermal strain of frozen soil caused by frost heave related to ice-water phase change.εth=εxνεyνεzνγxνγyνγzνT, where εth is the total strain matrix; εvx, εvy, and εvz are the typical strain components of thermal expansion in x, y, and z directions, respectively. For frozen isotopic soil, εxν=εyν=εzν=−αT−Tref, where α is the coefficient of linear thermal expansion (l/°C), which is related to the mass content and temperature; Tref is the reference temperature (°C); and γvx, γvy, and γvz are the components of shear stress induced by thermal expansion in the x, y, and z directions, respectively. For isotropic frozen soil, γvx=γvy=γvz=0, Fv is a volume force matrix.


## 2.4. Numerical Simulation-Theoretical Solution-Test Solution Method and Process

### 2.4.1. Simulation-Test Solution Process

(1)
Mathematical equations for the temperature field, moisture field, and stress field of frozen soil(2)
Input simulation frozen soil software, such as COMSOL, ANSYS, and ABAQUS (using MATLAB language, ABAQUS as solver development THM multifield coupling program, and analysis of frozen soil multifield coupling effect)(3)
Obtain a regular model consistent with the experiment


### 2.4.2. Theoretical Solution-Experimental Solution Process

(1)
Mathematical equations for temperature field, moisture field, and stress field of frozen soil(2)
When unknown variables are greater than the number of equations, supplementary equations need to be added(3)
Software MATLAB calculates the theoretical value(4)
Obtain the law theory consistent with the experiment


## 2.4.1. Simulation-Test Solution Process

(1)
Mathematical equations for the temperature field, moisture field, and stress field of frozen soil(2)
Input simulation frozen soil software, such as COMSOL, ANSYS, and ABAQUS (using MATLAB language, ABAQUS as solver development THM multifield coupling program, and analysis of frozen soil multifield coupling effect)(3)
Obtain a regular model consistent with the experiment


## 2.4.2. Theoretical Solution-Experimental Solution Process

(1)
Mathematical equations for temperature field, moisture field, and stress field of frozen soil(2)
When unknown variables are greater than the number of equations, supplementary equations need to be added(3)
Software MATLAB calculates the theoretical value(4)
Obtain the law theory consistent with the experiment


## 3. Current Research Status of THM Coupling for the Frozen Soil

### 3.1. Mechanism of Multiphysical Field Coupling for the Frozen Soil

At present, the theory of multiphysical field coupling is one of the most discussed topics related to frozen soil. It is also a cutting-edge field of international research: moisture migration, stress, and heat exchange jointly influences the frozen soil environment. If more factors of the frozen soil environment are considered, the calculation results will be more accurate [23, 39, 40].Some scholars [41–46] developed the THM coupling model by combining the TH coupling model with the porosity-strain relationship and conducted a further study on the THM coupling mechanism for the frozen soil by integrating Galerkin discretization, fluid dynamics theory, and an actual case. The phase change equation for the temperature field takes porosity as the variable, and the moisture migration equation is in Figure 2. Phenomenological thermodynamics, Fourier’s law, and the law of energy conservation form the basis for the above temperature, whereas the Richards equation, the law of conservation of energy, and Darcy’s law form the basis for the moisture field. The stress field depends on the Navier-Stokes equation and elastoplasticity theory.Figure 2
THM coupling diagram of frozen soil.For multiphysical field coupling, most scholars prefer the weak form theory of partial differential equation and attempt to connect the weight function in the weak state and the shape function in the finite element method. The existing governing equations of temperature, moisture, and stress fields weaken their derived forms. Besides, the tasks and governing equations are simplified and translated into the finite element model. The models have evolved from a simple form to a much more complex version. The algorithms are constantly modified so that it is more likely for the finite element simulation to achieve a convergent solution.


### 3.2. Development of THM Coupling for the Frozen Soil

Based on research on hydro-thermal (TH) coupling [47], the theory of hydro-thermal-mechanical (THM) coupling has been developed. A rigid ice model presumes the existence of the frozen fringe. Y. Zhou and G. Q. Zhou. [48] simplified the heat-fluid coupling model based on equivalent water pressure and explained the mechanism of the growth of a single ice lens (fittingdegree=0.3202). Despite the abundance of TH coupling models proposed so far, few of them deal with the problems caused by stress during the freezing process. Considering the stress change and the corresponding deformation, Li et al. [36] studied the phase change, convective heat transfer, and development of ice lens (fittingdegree=0.8930). Thomas et al. [49] built the heat-fluid coupling model (fittingdegree=0.2270). However, the two-field coupling is insufficient to characterize the stress pattern of pipelines in the frozen soil or realistically study the multiphysical coupling problem in the frozen soil. Based on the principle of effective stress, Nishimura et al. [50] proposed a THM coupling model to characterize the interactions in the saturated soil during the freeze-thaw cycles. Since unsaturated soils are prevalent in cold regions, the above model needs revision. Lai et al. [51] proposed a THM coupling model for frost heave based on the Clapeyron equation. They believed that the negative pore-water pressure was the primary driving force for moisture migration (fittingdegree=0.4058).The value (fitting degree=) is the fitting degree of the curve (test-simulation; test-calculated value) for time-frost heave amount/frost heave force. The smaller the value, the better the fit. Specific processes and code are described in the appendix. And the overall process is shown in Figure3.Figure 3
The flow chart for calculating the matching degree of the mathematical physics algorithm.The heat-fluid coupling model proposed by Y. Zhou and G. Q. Zhou [48] did not consider the mutual influence between stress development and heat transfer during the frost heave. Zhou et al. [52] described the THM multiphysical coupling model (fittingdegree=1.0643). Koniorczyk [53] conducted a theoretical study of the hydrothermal transfer process related to material deformation based on water phase-change dynamics. The volume averaging theory opens up an even broader horizon for THM coupling [54] . Based on this theory, Tong et al. and Na and Sun (2017) [55, 56] proposed the multiphase flow model for THM coupling and the porous medium model for stable THM coupling. Based on the volume average theory, the above models can reasonably simulate the freezing process where moisture migration, heat conduction, stress, and strain interact with each other in the soil in engineering construction.The interactions between pressures caused by frost heave and frosting have been rarely studied. Ji et al. [57] discussed the role of frosting-caused pressure in the frost heave equation (fittingdegree=0.7634). They were devoted to developing a one-dimensional device. The frosting-induced pressure affects the equivalent water pressure on the warm side of the ice lenses, which further influences the moisture flow velocity. The frosting-induced pressure between the soil particles can lead to the consolidation of liquid water. This pressure can also inhibit ice crystal formation, thus preventing frost heave. Therefore, new ice crystals are formed when the segregation pressure at the ice-water interface is larger than the sum of the frost heave pressure and the critical pressure.Considering the complexity of the nonlinear THM coupling model during the freezing process, Li et al. [36] employed the FORTRAN program to develop MHMIP-2D code, which provides new perspectives for the coupling study of the frozen soil.Research on the mechanism of complex frost heave path has been rarely reported before 2019. Bai et al. [58] studied the role of frost heave path in the frost heave mechanism based on the existing observations of ice lens distribution and the corresponding numerical simulation (fittingdegree=0.2655). They further proposed the three-dimensional THM coupling model with the presence of segregated ice.Lei et al. [12] simplified the three-phase heat transfer into single-phase heat transfer using the weighting algorithm. The relationship between soil freezing, porosity, water conductivity, and thermal conductivity of soil was discussed, and the new conditions for the formation of ice lenses were clarified (fittingdegree=1.6288). However, the defects of the above study lie in the following aspects. (1) The COMSOL simulation of soil freezing does not consider the compressibility of pore water, pore ice, and soil particles. Neither can this model be used to simulate the soil freezing process at a deeper layer. (2) The proposed model is specific for saturated soil, while unsaturated soil is prevalent in most cold regions. Therefore, the engineering applicability of the model is improved by changing the saturation when necessary. (3) No details are provided for the COMSOL simulation to solve the temperature field and displacement field. Therefore the coupling between the temperature field and the stress field may not be accurately solved.The THM coupling model will find more applications for boundary conditions, stress-strain relationships, and hydrothermal coupling processes.


## 3.1. Mechanism of Multiphysical Field Coupling for the Frozen Soil

At present, the theory of multiphysical field coupling is one of the most discussed topics related to frozen soil. It is also a cutting-edge field of international research: moisture migration, stress, and heat exchange jointly influences the frozen soil environment. If more factors of the frozen soil environment are considered, the calculation results will be more accurate [23, 39, 40].Some scholars [41–46] developed the THM coupling model by combining the TH coupling model with the porosity-strain relationship and conducted a further study on the THM coupling mechanism for the frozen soil by integrating Galerkin discretization, fluid dynamics theory, and an actual case. The phase change equation for the temperature field takes porosity as the variable, and the moisture migration equation is in Figure 2. Phenomenological thermodynamics, Fourier’s law, and the law of energy conservation form the basis for the above temperature, whereas the Richards equation, the law of conservation of energy, and Darcy’s law form the basis for the moisture field. The stress field depends on the Navier-Stokes equation and elastoplasticity theory.Figure 2
THM coupling diagram of frozen soil.For multiphysical field coupling, most scholars prefer the weak form theory of partial differential equation and attempt to connect the weight function in the weak state and the shape function in the finite element method. The existing governing equations of temperature, moisture, and stress fields weaken their derived forms. Besides, the tasks and governing equations are simplified and translated into the finite element model. The models have evolved from a simple form to a much more complex version. The algorithms are constantly modified so that it is more likely for the finite element simulation to achieve a convergent solution.


## 3.2. Development of THM Coupling for the Frozen Soil

Based on research on hydro-thermal (TH) coupling [47], the theory of hydro-thermal-mechanical (THM) coupling has been developed. A rigid ice model presumes the existence of the frozen fringe. Y. Zhou and G. Q. Zhou. [48] simplified the heat-fluid coupling model based on equivalent water pressure and explained the mechanism of the growth of a single ice lens (fittingdegree=0.3202). Despite the abundance of TH coupling models proposed so far, few of them deal with the problems caused by stress during the freezing process. Considering the stress change and the corresponding deformation, Li et al. [36] studied the phase change, convective heat transfer, and development of ice lens (fittingdegree=0.8930). Thomas et al. [49] built the heat-fluid coupling model (fittingdegree=0.2270). However, the two-field coupling is insufficient to characterize the stress pattern of pipelines in the frozen soil or realistically study the multiphysical coupling problem in the frozen soil. Based on the principle of effective stress, Nishimura et al. [50] proposed a THM coupling model to characterize the interactions in the saturated soil during the freeze-thaw cycles. Since unsaturated soils are prevalent in cold regions, the above model needs revision. Lai et al. [51] proposed a THM coupling model for frost heave based on the Clapeyron equation. They believed that the negative pore-water pressure was the primary driving force for moisture migration (fittingdegree=0.4058).The value (fitting degree=) is the fitting degree of the curve (test-simulation; test-calculated value) for time-frost heave amount/frost heave force. The smaller the value, the better the fit. Specific processes and code are described in the appendix. And the overall process is shown in Figure3.Figure 3
The flow chart for calculating the matching degree of the mathematical physics algorithm.The heat-fluid coupling model proposed by Y. Zhou and G. Q. Zhou [48] did not consider the mutual influence between stress development and heat transfer during the frost heave. Zhou et al. [52] described the THM multiphysical coupling model (fittingdegree=1.0643). Koniorczyk [53] conducted a theoretical study of the hydrothermal transfer process related to material deformation based on water phase-change dynamics. The volume averaging theory opens up an even broader horizon for THM coupling [54] . Based on this theory, Tong et al. and Na and Sun (2017) [55, 56] proposed the multiphase flow model for THM coupling and the porous medium model for stable THM coupling. Based on the volume average theory, the above models can reasonably simulate the freezing process where moisture migration, heat conduction, stress, and strain interact with each other in the soil in engineering construction.The interactions between pressures caused by frost heave and frosting have been rarely studied. Ji et al. [57] discussed the role of frosting-caused pressure in the frost heave equation (fittingdegree=0.7634). They were devoted to developing a one-dimensional device. The frosting-induced pressure affects the equivalent water pressure on the warm side of the ice lenses, which further influences the moisture flow velocity. The frosting-induced pressure between the soil particles can lead to the consolidation of liquid water. This pressure can also inhibit ice crystal formation, thus preventing frost heave. Therefore, new ice crystals are formed when the segregation pressure at the ice-water interface is larger than the sum of the frost heave pressure and the critical pressure.Considering the complexity of the nonlinear THM coupling model during the freezing process, Li et al. [36] employed the FORTRAN program to develop MHMIP-2D code, which provides new perspectives for the coupling study of the frozen soil.Research on the mechanism of complex frost heave path has been rarely reported before 2019. Bai et al. [58] studied the role of frost heave path in the frost heave mechanism based on the existing observations of ice lens distribution and the corresponding numerical simulation (fittingdegree=0.2655). They further proposed the three-dimensional THM coupling model with the presence of segregated ice.Lei et al. [12] simplified the three-phase heat transfer into single-phase heat transfer using the weighting algorithm. The relationship between soil freezing, porosity, water conductivity, and thermal conductivity of soil was discussed, and the new conditions for the formation of ice lenses were clarified (fittingdegree=1.6288). However, the defects of the above study lie in the following aspects. (1) The COMSOL simulation of soil freezing does not consider the compressibility of pore water, pore ice, and soil particles. Neither can this model be used to simulate the soil freezing process at a deeper layer. (2) The proposed model is specific for saturated soil, while unsaturated soil is prevalent in most cold regions. Therefore, the engineering applicability of the model is improved by changing the saturation when necessary. (3) No details are provided for the COMSOL simulation to solve the temperature field and displacement field. Therefore the coupling between the temperature field and the stress field may not be accurately solved.The THM coupling model will find more applications for boundary conditions, stress-strain relationships, and hydrothermal coupling processes.


## 4. Conclusions and Future Prospects

### 4.1. Conclusions

The latent heat of ice-water phase change, the influence of steam on the latent heat process of phase change, the latent heat of salt crystallization phase change, and the volume strain effect during temperature change are considered at the same time. TM coupling (equation (7)) has certain guiding significance for the study of the temperature field of saline soil in a seasonal frozen region. The TH coupling (equation (10)) considering the concentration change is also proposed. Quantitative evaluation of a matching degree about the test-model curve from a mathematical and physical perspective is proposed, which has certain guiding significance for the generalization of civil engineering.When there are more unknown variables of differential equations in multiphysical field coupling than differential equations, empirical equations and ice-water ratios are usually introduced to solve the differential equations. The empirical coefficients calculated from experiments may vary from one study to another, pointing to the necessity of improving the empirical coefficients for different soils in different regions.Freeze-thaw cycles are common occurrences in soils in seasonally frozen ground. It is necessary to introduce the differential term of the volumetric strain effect under temperature changes into the differential equation of the temperature field to realize thermo-mechanical coupling in the real sense.The COMSOL is a powerful tool for weighting mechanical, physical, heat transfer, and fluid parameters. This software fits actual engineering data to calculate the weights of different parameters for different engineering projects and thus achieve accurate THM coupling for the frozen soil. The refinement of test equipment, modern nondestructive detection methods, for example, real-time CT scanning technology, SEM, TEN, X-ray stereo ultrasonography, and nuclear magnetic resonance imaging (NMRI), and the improvement of the coupling between the stress boundary conditions and TH boundary conditions have laid a solid foundation for the multiphysical field coupling for the frozen soil.


### 4.2. Future Prospects

The coupling of the multiphysical field for the frozen soil is an integral part of structure and oil pipeline design and railway operation in cold regions. Multiphysical field coupling for frozen soil is a complex, multidisciplinary topic, and more work needs to be done to improve the prediction of the soil freezing process, coupling mechanism, coupling method, and solving approach [59–62]. However, many issues remain unsolved concerning the multiphysical field coupling of the frozen soil.
(1)
THM coupling for the frozen soil has been extensively studied. However, most of these studies only achieve weak coupling on the hydrothermal and stress boundaries. Although pairwise couplings have existed between the mechanical equation, heat transfer equation, and fluid equation, they were not coupling in the real sense(2)
There seems to be a lack of academic interest in the THVM and THCM coupling models for unsaturated soil under the freeze-thaw cycles. Besides, in most studies, the upper boundary is generally an open boundary. Therefore, one cannot depict the liquid water migration, vaporous water migration, ice content, heat distribution, and stress and strain status in the overlying roadbed in cold regions. There is a scarcity of theories and quantitative analysis tools for the coupling of water, ice, water vapor, temperature, stress, and strain when studying the damage mechanism of roadbed related to frost heave in these regions


## 4.1. Conclusions

The latent heat of ice-water phase change, the influence of steam on the latent heat process of phase change, the latent heat of salt crystallization phase change, and the volume strain effect during temperature change are considered at the same time. TM coupling (equation (7)) has certain guiding significance for the study of the temperature field of saline soil in a seasonal frozen region. The TH coupling (equation (10)) considering the concentration change is also proposed. Quantitative evaluation of a matching degree about the test-model curve from a mathematical and physical perspective is proposed, which has certain guiding significance for the generalization of civil engineering.When there are more unknown variables of differential equations in multiphysical field coupling than differential equations, empirical equations and ice-water ratios are usually introduced to solve the differential equations. The empirical coefficients calculated from experiments may vary from one study to another, pointing to the necessity of improving the empirical coefficients for different soils in different regions.Freeze-thaw cycles are common occurrences in soils in seasonally frozen ground. It is necessary to introduce the differential term of the volumetric strain effect under temperature changes into the differential equation of the temperature field to realize thermo-mechanical coupling in the real sense.The COMSOL is a powerful tool for weighting mechanical, physical, heat transfer, and fluid parameters. This software fits actual engineering data to calculate the weights of different parameters for different engineering projects and thus achieve accurate THM coupling for the frozen soil. The refinement of test equipment, modern nondestructive detection methods, for example, real-time CT scanning technology, SEM, TEN, X-ray stereo ultrasonography, and nuclear magnetic resonance imaging (NMRI), and the improvement of the coupling between the stress boundary conditions and TH boundary conditions have laid a solid foundation for the multiphysical field coupling for the frozen soil.


## 4.2. Future Prospects

The coupling of the multiphysical field for the frozen soil is an integral part of structure and oil pipeline design and railway operation in cold regions. Multiphysical field coupling for frozen soil is a complex, multidisciplinary topic, and more work needs to be done to improve the prediction of the soil freezing process, coupling mechanism, coupling method, and solving approach [59–62]. However, many issues remain unsolved concerning the multiphysical field coupling of the frozen soil.
(1)
THM coupling for the frozen soil has been extensively studied. However, most of these studies only achieve weak coupling on the hydrothermal and stress boundaries. Although pairwise couplings have existed between the mechanical equation, heat transfer equation, and fluid equation, they were not coupling in the real sense(2)
There seems to be a lack of academic interest in the THVM and THCM coupling models for unsaturated soil under the freeze-thaw cycles. Besides, in most studies, the upper boundary is generally an open boundary. Therefore, one cannot depict the liquid water migration, vaporous water migration, ice content, heat distribution, and stress and strain status in the overlying roadbed in cold regions. There is a scarcity of theories and quantitative analysis tools for the coupling of water, ice, water vapor, temperature, stress, and strain when studying the damage mechanism of roadbed related to frost heave in these regions



---

*Source: 8274137-2022-03-21.xml*