# Numerical Simulation of the Thermal-Hydro-Mechanical Characteristics of High-Speed Railway Roadbeds in Seasonally Frozen Regions

Authors: Yuzhi Zhang; Jianghui Bei; Pei Li; Xiaojie Liang
Journal: Advances in Civil Engineering
 (2020)

---

## Abstract

A multiphysics mathematical model of high-speed railway (HSR) roadbeds is necessary to facilitate a good level of understanding of the frost heaving mechanism. Based on the classical hydrodynamic model and fundamental thermoelasticity theories, we propose a thermo-hydro coupled model, based on the soil-water characteristic curve and solid-liquid ratio as the relation equations, with the effects of the ice-water phase change and water migration due to temperature change considered. With the linear expansion coefficient related to the temperature and the mass of ice content in roadbeds as the relation equation, we establish a macroscopic thermal-hydro-mechanical model for unsaturated soil to calculate the roadbed deformations. Based upon the field data of a typical cross section of the Harbin-Dalian HSR roadbed, the variation of the thermal-hydro-mechanical characteristics is simulated and studied. The results demonstrate that the increase of water content in the roadbed’s central line mainly appears in soil layers at depths less than 1.2 m and most ice-containing soil layers are at depths less than 0.6 m. Under the driving force of thermal and hydraulic migration, the vertical displacement of the west shoulder is increased to 18 mm. Then the settled maximum surface unevenness reaches 16 mm between the shoulder and centre line.


---

## Body

## 1. Introduction

The coarse-grained fillings and well-graded broken stone widely used in HSR roadbeds are generally regarded as materials with little or no frost heaving. However, field monitoring of the Harbin-Dalian HSR roadbed in seasonally regions of China reveals that frost heaving is commonly present [1, 2]. Hence, a complete understanding of the mechanism of freeze-thaw deformations of HSR roadbeds in cold regions is essential. Currently, the mechanism of freeze-thaw deformations in roadbeds is most often studied inside a laboratory by considering the frost heaving behaviour of unsaturated coarse-grained soil with respect to macroscopic factors such as particle size and grade [3, 4], but the heat and water migration in HSR roadbeds and the variation of their thermal-hydro-mechanical characteristics have seldom been studied.Due to limitations of the testing conditions and long duration, neither laboratory tests nor field monitoring can completely clarify the freeze-thaw deformations of coarse-grained fillings. Therefore, a thermal-hydro-mechanical numerical model is useful and facilitates an understanding of the mechanism of freeze-thaw deformations of coarse-grained fillings. Typical mechanical models of frozen soil already presented in the literature include a hydrodynamic model [5, 6], a segregation potential model [7, 8], a rigid ice model [9–12], and a thermodynamic model [13]. Based on the classical theory of hydrodynamics, various thermal-hydro-mechanical coupled models that consider the effects of volumetric strain caused by phase transformations between ice and water that occur in situ and the water migration in soil deformation have been proposed [14–16]. For example, Lai et al. [17] proposed a thermal-hydro-mechanical coupled model of saturated soil in which the ice content function is introduced and the Clapeyron equation is used as the equilibrium condition of the ice/water phase transformation. Similarly, thermal-hydro-mechanical coupled models have been established based on the partial expression of frozen soil components [18, 19]. Also, based on the similarity to frost heaving due to the ice/water phase transformation and the thermal linear expansion of materials, unsaturated thermal-hydro-mechanical frost heaving models have been proposed [20, 21].In previous studies, thermal-hydro-mechanical models were constructed by considering the Clapeyron equation generally used to describe phase equilibrium conditions. Also, the applicability of the Clapeyron equation under unsteady-state conditions is controversial, and moreover, it cannot distinguish between different soil characteristics such as particle composition, water content, and saturation degree, which are very important for practical engineering. The segregation criterion in the thermal-hydro-mechanical models is often empirical and cannot reach a unified judgment. A practical and realistic thermal-hydro-mechanical coupled model has not been put forward yet, since the existing thermal-hydro-mechanical models are often very complicated and some physical parameters are difficult to be measured. At the same time, previous studies explored moisture migration, frost heaving mechanism, and frost heaving variation of frozen soils in different ways, but most of the existing calculation models have only been verified through laboratory tests. Moreover, few of these consider the variation of the thermal-hydro-mechanical characteristics of HSR roadbeds in seasonally frozen regions, especially those with sunny-shady slopes.Therefore, based on the hydrodynamic model and the basic theory of elasticity, a unified soil-water characteristic curve (SWCC) is adopted to represent the soil water potential in unsaturated soil. The permeability coefficient related to the volumetric water content is introduced, and the parameter changes due to the water content in unsaturated soil are considered. The specific solid-liquid ratio is used as the coupling equation, and thus, the soil freezing law is looked at from the macroscopic perspective, instead of the microscopic study on the segregation criterion and the ice lens initiation mechanism. Meanwhile, the influence of the transformation of ice into water and water migration on the temperature field is considered, so as to realize the full hydrothermal coupling. The linear expansion coefficient, which is related to the soil temperature and mass of water content, is introduced into the mechanics of porous media equilibrium equation as the coupling equation between the hydrothermal state and the stress-strain relationship. The step function of hydrothermal characteristic parameters for frozen and unfrozen soil is used to accelerate the convergence. Subsequently, the thermal-hydro-mechanical mathematical model for unsaturated soil is constructed. Also, based on the field data of a typical section of the Harbin-Dalian HSR roadbed, a model considering the sunny-shady slope effect is developed and verified. Additionally, the variation of thermal-hydro-mechanical characteristics of HSR roadbeds in seasonally frozen regions is investigated, as well as the variations of the temperature field, moisture field, and deformation field.The organization of the remaining sections of this article is as follows. In Section2, the development of the thermal-hydro-mechanical model is presented, followed by an application of the model to the HSR roadbed in Section 3. In this latter section, the achieved results are presented and discussed. In Section 4, a conclusion is made.


## 2. Development of the Thermal-Hydro-Mechanical Model

In this section, the development of the thermal-hydro-mechanical model is presented.


### 2.1. Basic Assumptions

In deriving the model, the following basic assumptions were considered:(1)
The moisture migration mainly occurs at the frost front and unfrozen area, while that in the frozen area is negligible(2)
Moisture migration is only in the form of liquid water. In other words, migration of vapors (i.e., sublimation and condensation of moisture) and moisture evaporation at roadbed boundaries are excluded(3)
Moisture migration can be described by Darcy’s law(4)
Soil particles, ice, and water are incompressible substances with constant densities(5)
The change of solute is neglected, namely, the osmotic pressure caused by concentration difference is neglected(6)
Soil is an isotropic linear elastic continuum material that is completely consolidated and no account of creep is taken into consideration(7)
Once the mass of ice content exceeds a critical level, the linear expansion coefficient of soil follows a linear relationship with respect to the mass content of ice


### 2.2. Governing Equations

#### 2.2.1. Moisture Migration

According to the principle of mass conservation and Darcy’s law, the governing equation of moisture for unsaturated soil is given as follows [5]:(1)∂θu∂t+ρiρw∂θi∂t+∇−Kθu∇ψθu=0,where θu is the volumetric content of unfrozen water; t is the time (s); ρi is the ice density (kg/m3);ρw is the water density (kg/m3); θi is the volumetric content of ice; Kθu is the permeability coefficient (m/s); and ψ is the soil water potential.The soil water potential can be calculated by considering the generalized Clapeyron equation and the hydrodynamic model. Since the Clapeyron equation does not include any soil properties, it is only applicable to a state of phase equilibrium and can hardly measure the ice-water phase pressure. Hence, the soil water potential is calculated by the hydrodynamic model in this study. The soil water potential of unsaturated soil mainly consists of the matrix potential and the gravity potential. By introducing the diffusion coefficientDθu, equation (1) can be rewritten as follows:(2)∂θu∂t+ρiρw∂θi∂t+∇−Dθu∇θu−Kθu=0,(3)Dθu=Kθucm=Kθu∂θu/∂φpm=Kθu∂φpm∂θu.where Dθu is the water diffusion coefficient in soil (m2/s); Kθu is the soil permeability coefficient related with θu; Cm is the specific water capacity (m−1); and φpm is the matric potential (m). In most cases, the relation of matrix potential φpm and θu is described by the unsaturated soil-water characteristic curve (SWCC).


#### 2.2.2. Temperature Field

Since the temperature field distribution is affected by both the moisture field and ice/water phase transformation, the heat transfer induced by migrating water in the hydrodynamic model cannot be ignored [6]. Considering the effects of migrating water and the ice/water phase transformation on heat transfer in soil, the control equation for frozen soil temperature is given by(4)C∂T∂t+Cw−Kθu∇ψθu∇T+∇−λ∇T−Lρi∂θi∂t=0,where C is the soil volumetric heat capacity (kJ·m−3·°C −1); T is the temperature (°C); Cw is the volumetric heat capacity of water (kJ·m−3·°C−1); λ is the thermal conductivity (W·m−1·°C−1); and L is the latent heat of the ice (kJ/kg).


#### 2.2.3. Stress-Strain Relationship

In order to simplify the model and according to the basic assumptions, the effects of soil gravity and expansion induced by ice/water phase transformation on soil deformation are considered, but the effects of train loads on soil deformation are neglected. The stress-strain control equation of frozen soils is given by(5)∇Dε−εth+Fv=0,where D is the elasticity matrix, ε is the strain matrix, εth is the thermal strain matrix, and Fv is the volumetric force matrix:(6)D=E1+v1−2v1−vvv000v1−vv000vv1−v0000001−2v20000001−2v20000001−2v2,where E is the elastic modulus (Pa) and v is Poisson’s ratio.As for the thermal expansion of common materials, linear expansion coefficients are introduced to describe thermal strains in frozen soils induced by the ice/water phase transformation frost heaving:(7)εth=εxvεyvεzvγxvγyvγzvT,where εth is the total strain matrix, εxv, εyv, and εzv are components of the normal strain of thermal expansion in x, y, and z directions, respectively. For isotropic frozen soils, εxv=εxv=εxv=−αT−Tref, where α is the linear thermal expansion coefficient (1/°C) and it is related to the mass content of ice and temperature; Tref is the reference temperature (°C). γxv,γyv,andγzv are components of the shear strain of thermal expansion in x, y, and z directions, respectively. For isotropic frozen soils, γxv=γyv=γzv=0.


#### 2.2.4. Coupling Equations

From equations (2) and (4), two differential equations and three unknown basic variables (θu, θi, and T) are identified for coupling the moisture and temperature fields, and so one more coupling equation is required to close and solve the system of equations. Indeed, liquid water cannot completely transform into ice during freezing, and consequently, unfrozen water is present around the soil particles. The dynamic equilibrium of the volumetric content of unfrozen water and temperature can be characterized by the soil freezing characteristic curve (SFC) or other empirical equations. Based on the model for the prediction of unfrozen water [22], the solid-liquid ratio is introduced to describe the relation of the volumetric content of ice and volumetric content of unfrozen water and temperature, thus establishing a coupling equation. The ice-containing state of soil can be described by considering the temperature and the content of unfrozen water as shown in equations (8) and (9) to effectively improve the efficiency of the calculation:(8)w0wu=TTfb,where w0 is the initial gravimetric moisture content (w0=wi+wu); wi is the mass content of ice (wi=ρiθi/ρd); wu is the mass content of unfrozen water (wu=ρwθu/ρd); ρd is the soil dry density; T is the soil temperature (°C); b is an empirical constant related to soil texture [22]; and Tf is the freezing temperature of water (°C). Also,(9)Bi=θiθu=ρwρiTTfb−ρwρi,T<Tf,0,T≥Tf.where Bi is the solid-liquid ratio and θi can be described by θi and T. A complete hydrothermal coupling analysis can be achieved by considering equations (2), (4), and (9). Based on the solution of the moisture and temperature fields, the thermal-hydro-mechanical coupling analysis can be realized by combining equation (5) with the linear thermal expansion coefficient (α) as the coupling equation.


#### 2.2.5. Model Parameters

The ice content calculated by the hydrothermal model is combined with the hydrodynamic frost heaving model. The determined frost heaving coefficient based on the ice content was substituted by the thermal expansion coefficient in COMSOL (a Multiphysics software). The linear thermal expansion coefficient of clay loess in the foundation was set to be −0.0015, while the linear thermal expansion coefficients of well-graded broken stone and coarse-grained fillings in roadbeds [20] were determined by the following equation:(10)αjp=0.0156wi−0.0042,wi>0.027,0,wi≤0.027,αAB=0.076wi−0.0028,wi>0.037,0,wi≤0.037,where αjp is the linear thermal expansion coefficient of the well-graded broken stone layer (1/K); αAB is the linear expansion coefficient (1/K) of coarse-grained fillings; and wi is the mass content of ice.To facilitate the model’s numerical convergence, the hydrothermal parameters of soil [22] were selected by considering step functions so that these are continuous and differentiable:(11)C=ρdCs+ρwρdθuCw+ρiρdθiCi,Cs=Csf+Csu−CsfHTf,ΔT,where C is the volumetric thermal capacity of the soil (kJ·kg−1·°C−1); Cs is the specific heat capacity of soil skeleton (kJ·kg−1·°C−1); Cw is the specific heat capacity of water (kJ·kg−1·°C−1); Ci is the specific heat capacity of ice (kJ·kg−1·°C−1); Csf is the specific heat capacity of frozen soil skeleton (kJ·kg−1·°C−1); Csu is the specific heat capacity of thaw soil skeleton (kJ·kg−1·°C−1); H (Tf, ΔT) is the Heaviside step function, as shown in Figure 1; and ΔT is the 1/2 phase transition interval. Thus, the soil thermal conductivity λ (W·m−1·°C−1) is considered to vary according to the following equation:(12)λ=λf+λu−λfHTf,ΔT,where λf is the frozen soil thermal conductivity (W·m−1·°C−1) and λu is the thaw soil thermal conductivity (W·m−1·°C−1). Also, the soil permeability K (m/s) is considered to vary according to the following equation:(13)K=Kf+Ku−KfHTf,ΔT,where Kf is the frozen soil permeability (m/s) and Ku is the thaw soil permeability (m/s).Figure 1
Heaviside step function.Based on the model proposed by Van Genuchten [23], which is suitable for the calculation of the permeability coefficient of unsaturated soil, Ku and the specific water capacity Cm are calculated by the following equation:(14)Ku=KsKr=Ks1−ohk−11+ohk−m21+ohkm/2,where Ks is the saturated permeability coefficient (m/s); Kr is the relative permeability; o, m, and k are empirical constants related to soil properties [23]; and h is the pressure head (m), namely, φpm (m). And(15)Cm=−omθs−θr1−mθu−θrθs−θr1/m1−θu−θrθs−θr1/mm,where θs is the saturated volumetric content of water and θr is the volumetric content of residual water. Also,(16)Kf=KuIwhere I is the impedance coefficient. According to Taylor’s model [6, 14], I=10−10θi.In summary, the complete hydrothermal coupling analysis can be achieved by considering equations (2), (4), (9), and (11)–(16), and then, the thermal-hydro-mechanical analysis can be obtained by solving simultaneously equations (5) and (10).


## 2.1. Basic Assumptions

In deriving the model, the following basic assumptions were considered:(1)
The moisture migration mainly occurs at the frost front and unfrozen area, while that in the frozen area is negligible(2)
Moisture migration is only in the form of liquid water. In other words, migration of vapors (i.e., sublimation and condensation of moisture) and moisture evaporation at roadbed boundaries are excluded(3)
Moisture migration can be described by Darcy’s law(4)
Soil particles, ice, and water are incompressible substances with constant densities(5)
The change of solute is neglected, namely, the osmotic pressure caused by concentration difference is neglected(6)
Soil is an isotropic linear elastic continuum material that is completely consolidated and no account of creep is taken into consideration(7)
Once the mass of ice content exceeds a critical level, the linear expansion coefficient of soil follows a linear relationship with respect to the mass content of ice


## 2.2. Governing Equations

### 2.2.1. Moisture Migration

According to the principle of mass conservation and Darcy’s law, the governing equation of moisture for unsaturated soil is given as follows [5]:(1)∂θu∂t+ρiρw∂θi∂t+∇−Kθu∇ψθu=0,where θu is the volumetric content of unfrozen water; t is the time (s); ρi is the ice density (kg/m3);ρw is the water density (kg/m3); θi is the volumetric content of ice; Kθu is the permeability coefficient (m/s); and ψ is the soil water potential.The soil water potential can be calculated by considering the generalized Clapeyron equation and the hydrodynamic model. Since the Clapeyron equation does not include any soil properties, it is only applicable to a state of phase equilibrium and can hardly measure the ice-water phase pressure. Hence, the soil water potential is calculated by the hydrodynamic model in this study. The soil water potential of unsaturated soil mainly consists of the matrix potential and the gravity potential. By introducing the diffusion coefficientDθu, equation (1) can be rewritten as follows:(2)∂θu∂t+ρiρw∂θi∂t+∇−Dθu∇θu−Kθu=0,(3)Dθu=Kθucm=Kθu∂θu/∂φpm=Kθu∂φpm∂θu.where Dθu is the water diffusion coefficient in soil (m2/s); Kθu is the soil permeability coefficient related with θu; Cm is the specific water capacity (m−1); and φpm is the matric potential (m). In most cases, the relation of matrix potential φpm and θu is described by the unsaturated soil-water characteristic curve (SWCC).


### 2.2.2. Temperature Field

Since the temperature field distribution is affected by both the moisture field and ice/water phase transformation, the heat transfer induced by migrating water in the hydrodynamic model cannot be ignored [6]. Considering the effects of migrating water and the ice/water phase transformation on heat transfer in soil, the control equation for frozen soil temperature is given by(4)C∂T∂t+Cw−Kθu∇ψθu∇T+∇−λ∇T−Lρi∂θi∂t=0,where C is the soil volumetric heat capacity (kJ·m−3·°C −1); T is the temperature (°C); Cw is the volumetric heat capacity of water (kJ·m−3·°C−1); λ is the thermal conductivity (W·m−1·°C−1); and L is the latent heat of the ice (kJ/kg).


### 2.2.3. Stress-Strain Relationship

In order to simplify the model and according to the basic assumptions, the effects of soil gravity and expansion induced by ice/water phase transformation on soil deformation are considered, but the effects of train loads on soil deformation are neglected. The stress-strain control equation of frozen soils is given by(5)∇Dε−εth+Fv=0,where D is the elasticity matrix, ε is the strain matrix, εth is the thermal strain matrix, and Fv is the volumetric force matrix:(6)D=E1+v1−2v1−vvv000v1−vv000vv1−v0000001−2v20000001−2v20000001−2v2,where E is the elastic modulus (Pa) and v is Poisson’s ratio.As for the thermal expansion of common materials, linear expansion coefficients are introduced to describe thermal strains in frozen soils induced by the ice/water phase transformation frost heaving:(7)εth=εxvεyvεzvγxvγyvγzvT,where εth is the total strain matrix, εxv, εyv, and εzv are components of the normal strain of thermal expansion in x, y, and z directions, respectively. For isotropic frozen soils, εxv=εxv=εxv=−αT−Tref, where α is the linear thermal expansion coefficient (1/°C) and it is related to the mass content of ice and temperature; Tref is the reference temperature (°C). γxv,γyv,andγzv are components of the shear strain of thermal expansion in x, y, and z directions, respectively. For isotropic frozen soils, γxv=γyv=γzv=0.


### 2.2.4. Coupling Equations

From equations (2) and (4), two differential equations and three unknown basic variables (θu, θi, and T) are identified for coupling the moisture and temperature fields, and so one more coupling equation is required to close and solve the system of equations. Indeed, liquid water cannot completely transform into ice during freezing, and consequently, unfrozen water is present around the soil particles. The dynamic equilibrium of the volumetric content of unfrozen water and temperature can be characterized by the soil freezing characteristic curve (SFC) or other empirical equations. Based on the model for the prediction of unfrozen water [22], the solid-liquid ratio is introduced to describe the relation of the volumetric content of ice and volumetric content of unfrozen water and temperature, thus establishing a coupling equation. The ice-containing state of soil can be described by considering the temperature and the content of unfrozen water as shown in equations (8) and (9) to effectively improve the efficiency of the calculation:(8)w0wu=TTfb,where w0 is the initial gravimetric moisture content (w0=wi+wu); wi is the mass content of ice (wi=ρiθi/ρd); wu is the mass content of unfrozen water (wu=ρwθu/ρd); ρd is the soil dry density; T is the soil temperature (°C); b is an empirical constant related to soil texture [22]; and Tf is the freezing temperature of water (°C). Also,(9)Bi=θiθu=ρwρiTTfb−ρwρi,T<Tf,0,T≥Tf.where Bi is the solid-liquid ratio and θi can be described by θi and T. A complete hydrothermal coupling analysis can be achieved by considering equations (2), (4), and (9). Based on the solution of the moisture and temperature fields, the thermal-hydro-mechanical coupling analysis can be realized by combining equation (5) with the linear thermal expansion coefficient (α) as the coupling equation.


### 2.2.5. Model Parameters

The ice content calculated by the hydrothermal model is combined with the hydrodynamic frost heaving model. The determined frost heaving coefficient based on the ice content was substituted by the thermal expansion coefficient in COMSOL (a Multiphysics software). The linear thermal expansion coefficient of clay loess in the foundation was set to be −0.0015, while the linear thermal expansion coefficients of well-graded broken stone and coarse-grained fillings in roadbeds [20] were determined by the following equation:(10)αjp=0.0156wi−0.0042,wi>0.027,0,wi≤0.027,αAB=0.076wi−0.0028,wi>0.037,0,wi≤0.037,where αjp is the linear thermal expansion coefficient of the well-graded broken stone layer (1/K); αAB is the linear expansion coefficient (1/K) of coarse-grained fillings; and wi is the mass content of ice.To facilitate the model’s numerical convergence, the hydrothermal parameters of soil [22] were selected by considering step functions so that these are continuous and differentiable:(11)C=ρdCs+ρwρdθuCw+ρiρdθiCi,Cs=Csf+Csu−CsfHTf,ΔT,where C is the volumetric thermal capacity of the soil (kJ·kg−1·°C−1); Cs is the specific heat capacity of soil skeleton (kJ·kg−1·°C−1); Cw is the specific heat capacity of water (kJ·kg−1·°C−1); Ci is the specific heat capacity of ice (kJ·kg−1·°C−1); Csf is the specific heat capacity of frozen soil skeleton (kJ·kg−1·°C−1); Csu is the specific heat capacity of thaw soil skeleton (kJ·kg−1·°C−1); H (Tf, ΔT) is the Heaviside step function, as shown in Figure 1; and ΔT is the 1/2 phase transition interval. Thus, the soil thermal conductivity λ (W·m−1·°C−1) is considered to vary according to the following equation:(12)λ=λf+λu−λfHTf,ΔT,where λf is the frozen soil thermal conductivity (W·m−1·°C−1) and λu is the thaw soil thermal conductivity (W·m−1·°C−1). Also, the soil permeability K (m/s) is considered to vary according to the following equation:(13)K=Kf+Ku−KfHTf,ΔT,where Kf is the frozen soil permeability (m/s) and Ku is the thaw soil permeability (m/s).Figure 1
Heaviside step function.Based on the model proposed by Van Genuchten [23], which is suitable for the calculation of the permeability coefficient of unsaturated soil, Ku and the specific water capacity Cm are calculated by the following equation:(14)Ku=KsKr=Ks1−ohk−11+ohk−m21+ohkm/2,where Ks is the saturated permeability coefficient (m/s); Kr is the relative permeability; o, m, and k are empirical constants related to soil properties [23]; and h is the pressure head (m), namely, φpm (m). And(15)Cm=−omθs−θr1−mθu−θrθs−θr1/m1−θu−θrθs−θr1/mm,where θs is the saturated volumetric content of water and θr is the volumetric content of residual water. Also,(16)Kf=KuIwhere I is the impedance coefficient. According to Taylor’s model [6, 14], I=10−10θi.In summary, the complete hydrothermal coupling analysis can be achieved by considering equations (2), (4), (9), and (11)–(16), and then, the thermal-hydro-mechanical analysis can be obtained by solving simultaneously equations (5) and (10).


## 2.2.1. Moisture Migration

According to the principle of mass conservation and Darcy’s law, the governing equation of moisture for unsaturated soil is given as follows [5]:(1)∂θu∂t+ρiρw∂θi∂t+∇−Kθu∇ψθu=0,where θu is the volumetric content of unfrozen water; t is the time (s); ρi is the ice density (kg/m3);ρw is the water density (kg/m3); θi is the volumetric content of ice; Kθu is the permeability coefficient (m/s); and ψ is the soil water potential.The soil water potential can be calculated by considering the generalized Clapeyron equation and the hydrodynamic model. Since the Clapeyron equation does not include any soil properties, it is only applicable to a state of phase equilibrium and can hardly measure the ice-water phase pressure. Hence, the soil water potential is calculated by the hydrodynamic model in this study. The soil water potential of unsaturated soil mainly consists of the matrix potential and the gravity potential. By introducing the diffusion coefficientDθu, equation (1) can be rewritten as follows:(2)∂θu∂t+ρiρw∂θi∂t+∇−Dθu∇θu−Kθu=0,(3)Dθu=Kθucm=Kθu∂θu/∂φpm=Kθu∂φpm∂θu.where Dθu is the water diffusion coefficient in soil (m2/s); Kθu is the soil permeability coefficient related with θu; Cm is the specific water capacity (m−1); and φpm is the matric potential (m). In most cases, the relation of matrix potential φpm and θu is described by the unsaturated soil-water characteristic curve (SWCC).


## 2.2.2. Temperature Field

Since the temperature field distribution is affected by both the moisture field and ice/water phase transformation, the heat transfer induced by migrating water in the hydrodynamic model cannot be ignored [6]. Considering the effects of migrating water and the ice/water phase transformation on heat transfer in soil, the control equation for frozen soil temperature is given by(4)C∂T∂t+Cw−Kθu∇ψθu∇T+∇−λ∇T−Lρi∂θi∂t=0,where C is the soil volumetric heat capacity (kJ·m−3·°C −1); T is the temperature (°C); Cw is the volumetric heat capacity of water (kJ·m−3·°C−1); λ is the thermal conductivity (W·m−1·°C−1); and L is the latent heat of the ice (kJ/kg).


## 2.2.3. Stress-Strain Relationship

In order to simplify the model and according to the basic assumptions, the effects of soil gravity and expansion induced by ice/water phase transformation on soil deformation are considered, but the effects of train loads on soil deformation are neglected. The stress-strain control equation of frozen soils is given by(5)∇Dε−εth+Fv=0,where D is the elasticity matrix, ε is the strain matrix, εth is the thermal strain matrix, and Fv is the volumetric force matrix:(6)D=E1+v1−2v1−vvv000v1−vv000vv1−v0000001−2v20000001−2v20000001−2v2,where E is the elastic modulus (Pa) and v is Poisson’s ratio.As for the thermal expansion of common materials, linear expansion coefficients are introduced to describe thermal strains in frozen soils induced by the ice/water phase transformation frost heaving:(7)εth=εxvεyvεzvγxvγyvγzvT,where εth is the total strain matrix, εxv, εyv, and εzv are components of the normal strain of thermal expansion in x, y, and z directions, respectively. For isotropic frozen soils, εxv=εxv=εxv=−αT−Tref, where α is the linear thermal expansion coefficient (1/°C) and it is related to the mass content of ice and temperature; Tref is the reference temperature (°C). γxv,γyv,andγzv are components of the shear strain of thermal expansion in x, y, and z directions, respectively. For isotropic frozen soils, γxv=γyv=γzv=0.


## 2.2.4. Coupling Equations

From equations (2) and (4), two differential equations and three unknown basic variables (θu, θi, and T) are identified for coupling the moisture and temperature fields, and so one more coupling equation is required to close and solve the system of equations. Indeed, liquid water cannot completely transform into ice during freezing, and consequently, unfrozen water is present around the soil particles. The dynamic equilibrium of the volumetric content of unfrozen water and temperature can be characterized by the soil freezing characteristic curve (SFC) or other empirical equations. Based on the model for the prediction of unfrozen water [22], the solid-liquid ratio is introduced to describe the relation of the volumetric content of ice and volumetric content of unfrozen water and temperature, thus establishing a coupling equation. The ice-containing state of soil can be described by considering the temperature and the content of unfrozen water as shown in equations (8) and (9) to effectively improve the efficiency of the calculation:(8)w0wu=TTfb,where w0 is the initial gravimetric moisture content (w0=wi+wu); wi is the mass content of ice (wi=ρiθi/ρd); wu is the mass content of unfrozen water (wu=ρwθu/ρd); ρd is the soil dry density; T is the soil temperature (°C); b is an empirical constant related to soil texture [22]; and Tf is the freezing temperature of water (°C). Also,(9)Bi=θiθu=ρwρiTTfb−ρwρi,T<Tf,0,T≥Tf.where Bi is the solid-liquid ratio and θi can be described by θi and T. A complete hydrothermal coupling analysis can be achieved by considering equations (2), (4), and (9). Based on the solution of the moisture and temperature fields, the thermal-hydro-mechanical coupling analysis can be realized by combining equation (5) with the linear thermal expansion coefficient (α) as the coupling equation.


## 2.2.5. Model Parameters

The ice content calculated by the hydrothermal model is combined with the hydrodynamic frost heaving model. The determined frost heaving coefficient based on the ice content was substituted by the thermal expansion coefficient in COMSOL (a Multiphysics software). The linear thermal expansion coefficient of clay loess in the foundation was set to be −0.0015, while the linear thermal expansion coefficients of well-graded broken stone and coarse-grained fillings in roadbeds [20] were determined by the following equation:(10)αjp=0.0156wi−0.0042,wi>0.027,0,wi≤0.027,αAB=0.076wi−0.0028,wi>0.037,0,wi≤0.037,where αjp is the linear thermal expansion coefficient of the well-graded broken stone layer (1/K); αAB is the linear expansion coefficient (1/K) of coarse-grained fillings; and wi is the mass content of ice.To facilitate the model’s numerical convergence, the hydrothermal parameters of soil [22] were selected by considering step functions so that these are continuous and differentiable:(11)C=ρdCs+ρwρdθuCw+ρiρdθiCi,Cs=Csf+Csu−CsfHTf,ΔT,where C is the volumetric thermal capacity of the soil (kJ·kg−1·°C−1); Cs is the specific heat capacity of soil skeleton (kJ·kg−1·°C−1); Cw is the specific heat capacity of water (kJ·kg−1·°C−1); Ci is the specific heat capacity of ice (kJ·kg−1·°C−1); Csf is the specific heat capacity of frozen soil skeleton (kJ·kg−1·°C−1); Csu is the specific heat capacity of thaw soil skeleton (kJ·kg−1·°C−1); H (Tf, ΔT) is the Heaviside step function, as shown in Figure 1; and ΔT is the 1/2 phase transition interval. Thus, the soil thermal conductivity λ (W·m−1·°C−1) is considered to vary according to the following equation:(12)λ=λf+λu−λfHTf,ΔT,where λf is the frozen soil thermal conductivity (W·m−1·°C−1) and λu is the thaw soil thermal conductivity (W·m−1·°C−1). Also, the soil permeability K (m/s) is considered to vary according to the following equation:(13)K=Kf+Ku−KfHTf,ΔT,where Kf is the frozen soil permeability (m/s) and Ku is the thaw soil permeability (m/s).Figure 1
Heaviside step function.Based on the model proposed by Van Genuchten [23], which is suitable for the calculation of the permeability coefficient of unsaturated soil, Ku and the specific water capacity Cm are calculated by the following equation:(14)Ku=KsKr=Ks1−ohk−11+ohk−m21+ohkm/2,where Ks is the saturated permeability coefficient (m/s); Kr is the relative permeability; o, m, and k are empirical constants related to soil properties [23]; and h is the pressure head (m), namely, φpm (m). And(15)Cm=−omθs−θr1−mθu−θrθs−θr1/m1−θu−θrθs−θr1/mm,where θs is the saturated volumetric content of water and θr is the volumetric content of residual water. Also,(16)Kf=KuIwhere I is the impedance coefficient. According to Taylor’s model [6, 14], I=10−10θi.In summary, the complete hydrothermal coupling analysis can be achieved by considering equations (2), (4), (9), and (11)–(16), and then, the thermal-hydro-mechanical analysis can be obtained by solving simultaneously equations (5) and (10).


## 3. Application of the Developed Model

### 3.1. Application to the Roadbed

Currently, in frost heaving models, the boundary temperatures are usually simplified in practical applications by a symmetric boundary condition. Considering a section of the Harbin-Dalian HSR as an object, the thermal-hydro-mechanical relation of the HSR roadbeds in seasonally frozen regions is investigated using field data and with the sunny-shady slope effect taken into consideration, the effect of which was previously reported by Zhang et al. [24] in the north-south oriented part of the Harbin-Dalian HSR.Steady-state monitored field data were obtained, and the results were used for a transient analysis of the initial temperature field. Once the initial temperature field is established, the thermal-hydro-mechanical coupling is divided into two steps. First, the hydrothermal coupling analysis of the roadbeds is performed by secondary development of the mathematical module in COMSOL, to obtain the distributions and variations of moisture and temperature fields. Second, stress and deformation fields are calculated by the geomechanics module in COMSOL, which is based on the proposed model and the results of the hydrothermal coupling analysis.As shown in Figure2(a), the target section is located near the Shuangcheng Station in Harbin and its altitude is 166.4 m above mean sea level [25]. The left part of Figure 2(b) shows the composition of each roadbed component [25]. The width and the height of the roadbed are 13.6 m and 5.433 m, respectively, and the roadbed upper layer is a 0.4 m layer of well-graded broken stone. A composite geomembrane insulating layer is laid just below the upper layer with a 0.05 m medium and coarse sand cushion below and above, respectively. The part under the insulating layer consists of a 1.0 m antifreezing layer with coarse-grained fillings (High quality and good filling materials in China railway roadbed design code [26]), while the lowest layer is composed of normal coarse-grained fillings. An antifrost heaving berm is set at the foot of the embankment slope. The foundation soil is clayey loess, silty clay, and sand, reinforced by CFG (cement fly-ash gravel) piles arranged in equilateral triangles. No surface water spillover was observed at the test site. The groundwater type is Quaternary porous phreatic water, which is mainly recharged by atmospheric precipitation and exposed to severe seasonal variation. The water depth is relatively large (between 13 m and 18.5 m), and the groundwater level has an amplitude approximately between 1 m and 2 m.Figure 2
Location, composition, and measurement arrangement of the observed section. (a) Location of the observed section. (b) Composition and arrangement of measured points of roadbed [25].
(a)(b)The right part of Figure2(b) shows the arrangement of the roadbed measuring points considered [25]. Slope foots and temperature measuring holes are only arranged on the west part of the roadbed and other measuring points are arranged symmetrically with respect to the central line of the roadbed. Temperature measuring points are arranged at intervals of 0.5 m between 0.8 m and 3.8 m at intervals of 1 m between 3.8 m and 10.8 m below the roadbed surface. Five measuring lines are arranged at the east shoulder, west shoulder, central line, west slope foot, and 20 m away from the west slope foot, respectively, and each measuring line consists of 14 measuring points.With the thermistor as the main component, the temperature sensor has a measuring range of −40°C to +60°C. Its accuracy is ±0.03°C in the range of −20°C to +20°C.At the target section, well-graded broken stone and tracks were laid in August 2010 and September-October 2010, respectively. The temperature was measured every five days from August 01, 2010, to August 01, 2014. The ice densityρi and water density ρw were 917 and 1000 kg m−3. The ice thermal conductivity λi and the water thermal conductivity λw were 2.22 and 0.58 W m−1·°C−1, respectively. The latent heat of the ice L was 334.56 kJ kg−3. And the specific heat capacity of ice Ci and the specific heat capacity of water Cw were 2.09 and 4.18 kJ kg−1·°C−1. Tables 1 and 2 summarize the basic characteristic parameters of the roadbed and foundation layers [22, 27].Table 1
Thermal characteristic parameters of soil.
Csf (kJ·kg−1·°C−1)Csu (kJ·kg−1·°C−1)Tf (°C)λf (W·m−1·°C−1)λu (W·m−1·°C−1)Well-graded broken stone0.710.79−0.562.531.86Coarse-grained fillings0.620.66−0.561.401.15Clay loess0.760.84−0.452.401.54Silty clay0.760.84−1.002.121.42Silty sand0.730.84−0.601.241.86Table 2
Characteristic parameters of water for different soils.
w0 (%)α (m−1)bKs (m·s−1)lmnρd (kg·m−3)θsθrWell-graded broken stone100.660.612E−40.50.141.161 7800.400.05Coarse-grained fillings90.660.611E−40.50.141.162 1220.400.05Clayey loess252.650.568E−80.50.261.351 3600.530.08Silty clay252.650.562E−80.50.261.351 6000.510.08Silty sand251.000.476E−70.50.261.351 5000.410.06Considering the effects of temperature on the elastic modulus and Poisson’s ratio, the mechanical parameters of roadbeds are determined by the method reported elsewhere (see Table3) [24].Table 3
Mechanical parameters of roadbeds.
Temperature (°C)−20−10−5−1−0.5030E (MPa)380.81272.00200.21114.0095.9761.0061.00v0.210.280.3150.3450.3450.350.35By fitting a function of the ground temperature obtained from field monitoring data, as shown in Figure3, the Dirichlet temperature boundary conditions considering the sunny-shady slopes effect for different locations are designed [27]:(17)Tx,y,t=T0+A0sin2πtp+ω+R0t,x,y∈AB¯,BC¯,CD¯,DE¯,EF¯,FG¯,GH¯,HI,¯andIJ¯,∂Tx,y∂n=0,x,y∈AL¯andJK¯,Tx,y=T1,x,y∈KL¯.Figure 3
Sketch map of roadbed model [27].The initial phase angleω is determined according to the initial data obtained from the monitored area. Based on the boundary layer principle (boundary layer thickness was considered as 0.5 m) [24], the average ground temperatures T0 at the east shoulder, central line, and west shoulder of the roadbed were determined to be 9.55°C, 7.71°C, and 7.80°C, respectively. A0 is the temperature amplitude of the boundary layer bottom, t is the time from the initial moment, measured in hours, and p is the vibration period, adopted as 8760 hours. R0 is the thermal gradient of the boundary layer bottom temperature increased by global warming, adopted as 0.04°C/8760 h [28], n is the outer normal direction of the boundary, and T1 is the temperature at the lower boundary of the calculated region.The dimensions of the model are shown in Figure3 [27]. The height of clayey loess, silty clay, and sand was determined to be 11 m, 11 m, and 3 m, respectively, by drilling.East shoulder to central line (the depth from the roadbed surfacez range = 0–5.433 m):(18)Tx,y,t=9.55−1.84x6.8+18.560−1.920x6.8exp−0.2957zsin2πtp+π2−πx10×6.8−0.2957z.Central line to west shoulder (the depth from the roadbed surfacez range = 0–5.433 m):(19)Tx,y,t=7.71+0.09x−6.86.8+16.640−2.70x−6.86.8exp−0.2957zsin2πtp+2π5−πx−6.810×6.8−0.2957z,where x is the distance from the east shoulder (0–13.8 m).To discretize the variations of the moisture and temperature fields of roadbeds in different seasons, the total time step and substep of the hydrothermal coupling transient analysis considered are 365 days and 0.5 days, respectively. The starting date of the transient analysis is August 1 (temperature peak moment).


### 3.2. Results and Discussion

#### 3.2.1. Model Verification

Figure4 shows a comparison between the simulated frost depth (SFR) and the monitored frost depth (MFR) of the west shoulder and a comparison between the simulated vertical displacement (SVD) of the roadbed and the monitored vertical displacement (MVD) within 2.7 m from the adjacent roadbed section [29]. As it is evident, simulated and monitored frost depths follow consistent trends. Simulation of the deformation demonstrates a reasonable match to the monitored deformation basically. The region with frost depth less than 1.2 m is exposed to severe soil frost heaving. When the frost depth exceeded 1.2 m, frost heaving developed slowly even though the frost depth increased further. The duration of thawing subsidence was short. As shown in Figure 4, field monitoring of the adjacent roadbed section indicated 2 mm residual deformation after thawing subsidence of the roadbed [29], while previous studies claimed that roadbed deformations almost returned to the initial condition after thawing subsidence. In some cases, a small amount of settlements were observed after thawing subsidence [27].Figure 4
Simulated and monitored frost depth and vertical displacement.Therefore, from the results, it can be concluded that the proposed numerical model is effective and applicable for a comprehensive analysis of the temperature field, moisture field, and deformation field of the roadbed.


#### 3.2.2. Analysis of the Simulation Results

Figure5 illustrates the temperature field distribution of the section at certain moments in time. As observed, the freezing process is unidirectional and lasted for a long period, while the thawing process is bidirectional and has a short period. Due to the sunny-shady slope effect, the temperature fields of the roadbed are asymmetrically distributed in the transverse direction during the freezing-thawing process. The maximum frost depth occurred at the west shoulder (right shoulder, sunny shoulder), followed by the east shoulder (left shoulder, shady shoulder) and that at the central line. The freezing sequence is west side slope, east side slope, and top surface. During thawing, frozen soil nuclei are observed in the roadbed and the thawing sequence is west side slope, east side slope, and roadbed top surface. Additionally, thawing of shoulders occurred earlier than that of the slope foot.Figure 5
Temperature field distribution of the section at a typical moment (°C). Temperature field distribution on (a) August 01, 2013, (b) December 10, 2013, (c) March 10, 2014, (d) April 10, 2014, and (e) August 01, 2014.
(a)(b)(c)(d)(e)Figure6 illustrates the distribution of the mass content of ice at the roadbed section. Due to the presence of the temperature field, ice-containing soil layers are observed at depths less than 1.2 m. At the same moment, the distribution of ice content of the roadbed section is uneven; the ice content of the west part is relatively large during freezing and melted rapidly during thawing.Figure 6
Distribution of the mass content of ice at the roadbed section on March 10, 2014 (%).Figure7 shows time history curves of the temperature field and moisture fields of the central line of the roadbed. As observed, the ice-containing soil layer at the central line has a thickness of about 1.2 m, mainly within the range 0 to 0.6 m. Moisture migration is observed in both positive and negative temperature periods, but most especially in the negative temperature period. The moisture migration is concentrated in the region at depths less than 2.5 m. Here, the moisture distribution in positive and negative temperature periods is mainly affected by self-equilibrium of the moisture and temperature fields, respectively. During freezing, the mass content of moisture in the roadbed increased with depths less than 1.2 m decreased with depths between 1.2 m and 2.5 m and remained almost constant with depths greater than 2.5 m. Here, the mass content of moisture increased near the top surface of the roadbed and within the layer at a depth of 1 m. The maximum variation of the mass content of moisture is 2%, which was about 20% higher than that before freezing. In other words, moisture migration has a significant effect. The critical depth of moisture variation is 1.2 m. At this depth, the variation range of temperature is −0.7∼0°C (with the phase transition temperature zone) during freezing. Also, a frost front is presented, resulting in the migration of moisture to the lower part of the roadbed from an unfrozen region to a frozen one.Figure 7
Trends of temperature field and moisture field of the central line. (a) Temperature (°C). (b) Mass content of ice (%). (c) Mass content of unfrozen water (%). (d) Mass content of moisture (%).
(a)(b)(c)(d)The distribution of the deformation of roadbeds is investigated based on the above analysis. In this work, the vertical displacement of the roadbed surface was the major concern. Figure8 shows the time history curve, temperature curve, frost depth curve of the roadbed deformations. Figure 9 shows vertical displacements of roadbed section at different moments. Figure 10 shows the horizontal and vertical displacement fields of the roadbed section on March 10, 2014. As it is evident, vertical displacement of the roadbed surface followed the sequence of west shoulder, east shoulder, and central line, while the vertical displacement followed the sequence of west shoulder which is higher that of the east shoulder which is higher than that at the central line. The maximum vertical displacement of the roadbed surface is approximately 11∼18 mm (15 mm is allowed according to Chinese railway roadbed design code). At the early stage of freezing, the vertical displacement increased with the frost depth. Once the frost depth exceeded a critical level, the vertical displacement did not increase further with the frost depth. The roadbed vertical displacements are asymmetrically distributed in the transverse direction and also an uneven deformation of the roadbed surface appeared. In this study, the left and right shoulders are exposed to ambient conditions, while the central line is covered by the track slab. Hence, roadbed deformations at the shoulders and the central line are uneven in the early negative temperature period (up to 14 mm) and the difference subsequently decreases.Figure 8
Time history curves of vertical displacement and frost depth.Figure 9
Vertical displacements of roadbed section at different moments.Figure 10
Distribution of roadbed section on March 10, 2014. (a) Vertical displacement field. (b) Horizontal displacement field.
(a)(b)


#### 3.2.3. Discussion

For frost heaving of roadbeds in cold regions, changes of the temperature, moisture, and deformation are the cause, the medium, and the result, respectively. Roadbed deformations are macroscopic effects of hydrothermal migrations in roadbed soils.Field monitoring revealed that frost heaving of the Harbin-Dalian HSR roadbed is present from the surface to a frost depth of 1.2 m. Here, frost heaving of the well-graded broken stone dominated the frost heaving of the roadbed to a frost depth of 0.5 m [27]. Additionally, the moisture content varied drastically within the depth of 0.6 m.The results of this study demonstrated that moisture migration at the freezing front is observed in the roadbed of the Harbin-Dalian HSR under the driving force of the temperature field. Here, the moisture content increased significantly in the range between 0 m and 1.2 m, and the ice-containing soil layers are concentrated in this region especially in the range between 0 and 0.6 m. Therefore, phase transition and moisture migration in Harbin-Dalian HSR roadbed driven by an uneven distribution of the temperature field are the main causes of roadbed deformations.The transverse variation of the ground temperature of roadbeds can be attributed to the differences in the solar radiation intensity and also the surface covering layer [24, 30]. These variations of the ground temperature lead to an unbalanced development of the frost depth time and value, thus resulting in an unbalanced moisture migration. In this case, uneven transverse deformations are inevitable [30, 31]. Field monitoring data revealed differences in deformations of the east and west shoulders of the Harbin-Dalian HSR roadbed [32]. Simulation results demonstrated that the distribution of deformation field of roadbed in the transverse direction is asymmetric and uneven due to the effects of temperature and mass content of ice. Indeed, the deformation is maximized at the shoulders. The deformation of a specific part of roadbed is affected by a variation of the ground temperature and moisture at this part and of the surrounding soil. Hence, roadbed deformations are highly complicated. Frost heaving was first observed at the west shoulder, since freezing was first observed over there. However, deformation was initiated before freezing, since the freezing of the side slope occurred earlier than roadbed freezing. The uneven transverse deformation of the roadbed results in severe frost heaving at the shoulders and side slope. In fact, slope foots should be protected as weak spots due to severe ice accumulation results in horizontal displacement.The results demonstrated that a reduction of the roadbed freeze-thaw deformations should be achieved based on the analysis of thermo-hydro migration characteristics of roadbeds, instead of a thermal insulation layer on the roadbed surface. Elimination of the uneven transverse deformations should be obtained by minimizing the difference of the ground temperature in the transverse direction, which can be achieved by asymmetric transitional layout of roadbed soils [30] or by an asymmetric antifrost heaving berm [33]. Based on a theoretical study, this study provides references to the design and maintenance of HSR roadbeds in seasonally frozen regions.


## 3.1. Application to the Roadbed

Currently, in frost heaving models, the boundary temperatures are usually simplified in practical applications by a symmetric boundary condition. Considering a section of the Harbin-Dalian HSR as an object, the thermal-hydro-mechanical relation of the HSR roadbeds in seasonally frozen regions is investigated using field data and with the sunny-shady slope effect taken into consideration, the effect of which was previously reported by Zhang et al. [24] in the north-south oriented part of the Harbin-Dalian HSR.Steady-state monitored field data were obtained, and the results were used for a transient analysis of the initial temperature field. Once the initial temperature field is established, the thermal-hydro-mechanical coupling is divided into two steps. First, the hydrothermal coupling analysis of the roadbeds is performed by secondary development of the mathematical module in COMSOL, to obtain the distributions and variations of moisture and temperature fields. Second, stress and deformation fields are calculated by the geomechanics module in COMSOL, which is based on the proposed model and the results of the hydrothermal coupling analysis.As shown in Figure2(a), the target section is located near the Shuangcheng Station in Harbin and its altitude is 166.4 m above mean sea level [25]. The left part of Figure 2(b) shows the composition of each roadbed component [25]. The width and the height of the roadbed are 13.6 m and 5.433 m, respectively, and the roadbed upper layer is a 0.4 m layer of well-graded broken stone. A composite geomembrane insulating layer is laid just below the upper layer with a 0.05 m medium and coarse sand cushion below and above, respectively. The part under the insulating layer consists of a 1.0 m antifreezing layer with coarse-grained fillings (High quality and good filling materials in China railway roadbed design code [26]), while the lowest layer is composed of normal coarse-grained fillings. An antifrost heaving berm is set at the foot of the embankment slope. The foundation soil is clayey loess, silty clay, and sand, reinforced by CFG (cement fly-ash gravel) piles arranged in equilateral triangles. No surface water spillover was observed at the test site. The groundwater type is Quaternary porous phreatic water, which is mainly recharged by atmospheric precipitation and exposed to severe seasonal variation. The water depth is relatively large (between 13 m and 18.5 m), and the groundwater level has an amplitude approximately between 1 m and 2 m.Figure 2
Location, composition, and measurement arrangement of the observed section. (a) Location of the observed section. (b) Composition and arrangement of measured points of roadbed [25].
(a)(b)The right part of Figure2(b) shows the arrangement of the roadbed measuring points considered [25]. Slope foots and temperature measuring holes are only arranged on the west part of the roadbed and other measuring points are arranged symmetrically with respect to the central line of the roadbed. Temperature measuring points are arranged at intervals of 0.5 m between 0.8 m and 3.8 m at intervals of 1 m between 3.8 m and 10.8 m below the roadbed surface. Five measuring lines are arranged at the east shoulder, west shoulder, central line, west slope foot, and 20 m away from the west slope foot, respectively, and each measuring line consists of 14 measuring points.With the thermistor as the main component, the temperature sensor has a measuring range of −40°C to +60°C. Its accuracy is ±0.03°C in the range of −20°C to +20°C.At the target section, well-graded broken stone and tracks were laid in August 2010 and September-October 2010, respectively. The temperature was measured every five days from August 01, 2010, to August 01, 2014. The ice densityρi and water density ρw were 917 and 1000 kg m−3. The ice thermal conductivity λi and the water thermal conductivity λw were 2.22 and 0.58 W m−1·°C−1, respectively. The latent heat of the ice L was 334.56 kJ kg−3. And the specific heat capacity of ice Ci and the specific heat capacity of water Cw were 2.09 and 4.18 kJ kg−1·°C−1. Tables 1 and 2 summarize the basic characteristic parameters of the roadbed and foundation layers [22, 27].Table 1
Thermal characteristic parameters of soil.
Csf (kJ·kg−1·°C−1)Csu (kJ·kg−1·°C−1)Tf (°C)λf (W·m−1·°C−1)λu (W·m−1·°C−1)Well-graded broken stone0.710.79−0.562.531.86Coarse-grained fillings0.620.66−0.561.401.15Clay loess0.760.84−0.452.401.54Silty clay0.760.84−1.002.121.42Silty sand0.730.84−0.601.241.86Table 2
Characteristic parameters of water for different soils.
w0 (%)α (m−1)bKs (m·s−1)lmnρd (kg·m−3)θsθrWell-graded broken stone100.660.612E−40.50.141.161 7800.400.05Coarse-grained fillings90.660.611E−40.50.141.162 1220.400.05Clayey loess252.650.568E−80.50.261.351 3600.530.08Silty clay252.650.562E−80.50.261.351 6000.510.08Silty sand251.000.476E−70.50.261.351 5000.410.06Considering the effects of temperature on the elastic modulus and Poisson’s ratio, the mechanical parameters of roadbeds are determined by the method reported elsewhere (see Table3) [24].Table 3
Mechanical parameters of roadbeds.
Temperature (°C)−20−10−5−1−0.5030E (MPa)380.81272.00200.21114.0095.9761.0061.00v0.210.280.3150.3450.3450.350.35By fitting a function of the ground temperature obtained from field monitoring data, as shown in Figure3, the Dirichlet temperature boundary conditions considering the sunny-shady slopes effect for different locations are designed [27]:(17)Tx,y,t=T0+A0sin2πtp+ω+R0t,x,y∈AB¯,BC¯,CD¯,DE¯,EF¯,FG¯,GH¯,HI,¯andIJ¯,∂Tx,y∂n=0,x,y∈AL¯andJK¯,Tx,y=T1,x,y∈KL¯.Figure 3
Sketch map of roadbed model [27].The initial phase angleω is determined according to the initial data obtained from the monitored area. Based on the boundary layer principle (boundary layer thickness was considered as 0.5 m) [24], the average ground temperatures T0 at the east shoulder, central line, and west shoulder of the roadbed were determined to be 9.55°C, 7.71°C, and 7.80°C, respectively. A0 is the temperature amplitude of the boundary layer bottom, t is the time from the initial moment, measured in hours, and p is the vibration period, adopted as 8760 hours. R0 is the thermal gradient of the boundary layer bottom temperature increased by global warming, adopted as 0.04°C/8760 h [28], n is the outer normal direction of the boundary, and T1 is the temperature at the lower boundary of the calculated region.The dimensions of the model are shown in Figure3 [27]. The height of clayey loess, silty clay, and sand was determined to be 11 m, 11 m, and 3 m, respectively, by drilling.East shoulder to central line (the depth from the roadbed surfacez range = 0–5.433 m):(18)Tx,y,t=9.55−1.84x6.8+18.560−1.920x6.8exp−0.2957zsin2πtp+π2−πx10×6.8−0.2957z.Central line to west shoulder (the depth from the roadbed surfacez range = 0–5.433 m):(19)Tx,y,t=7.71+0.09x−6.86.8+16.640−2.70x−6.86.8exp−0.2957zsin2πtp+2π5−πx−6.810×6.8−0.2957z,where x is the distance from the east shoulder (0–13.8 m).To discretize the variations of the moisture and temperature fields of roadbeds in different seasons, the total time step and substep of the hydrothermal coupling transient analysis considered are 365 days and 0.5 days, respectively. The starting date of the transient analysis is August 1 (temperature peak moment).


## 3.2. Results and Discussion

### 3.2.1. Model Verification

Figure4 shows a comparison between the simulated frost depth (SFR) and the monitored frost depth (MFR) of the west shoulder and a comparison between the simulated vertical displacement (SVD) of the roadbed and the monitored vertical displacement (MVD) within 2.7 m from the adjacent roadbed section [29]. As it is evident, simulated and monitored frost depths follow consistent trends. Simulation of the deformation demonstrates a reasonable match to the monitored deformation basically. The region with frost depth less than 1.2 m is exposed to severe soil frost heaving. When the frost depth exceeded 1.2 m, frost heaving developed slowly even though the frost depth increased further. The duration of thawing subsidence was short. As shown in Figure 4, field monitoring of the adjacent roadbed section indicated 2 mm residual deformation after thawing subsidence of the roadbed [29], while previous studies claimed that roadbed deformations almost returned to the initial condition after thawing subsidence. In some cases, a small amount of settlements were observed after thawing subsidence [27].Figure 4
Simulated and monitored frost depth and vertical displacement.Therefore, from the results, it can be concluded that the proposed numerical model is effective and applicable for a comprehensive analysis of the temperature field, moisture field, and deformation field of the roadbed.


### 3.2.2. Analysis of the Simulation Results

Figure5 illustrates the temperature field distribution of the section at certain moments in time. As observed, the freezing process is unidirectional and lasted for a long period, while the thawing process is bidirectional and has a short period. Due to the sunny-shady slope effect, the temperature fields of the roadbed are asymmetrically distributed in the transverse direction during the freezing-thawing process. The maximum frost depth occurred at the west shoulder (right shoulder, sunny shoulder), followed by the east shoulder (left shoulder, shady shoulder) and that at the central line. The freezing sequence is west side slope, east side slope, and top surface. During thawing, frozen soil nuclei are observed in the roadbed and the thawing sequence is west side slope, east side slope, and roadbed top surface. Additionally, thawing of shoulders occurred earlier than that of the slope foot.Figure 5
Temperature field distribution of the section at a typical moment (°C). Temperature field distribution on (a) August 01, 2013, (b) December 10, 2013, (c) March 10, 2014, (d) April 10, 2014, and (e) August 01, 2014.
(a)(b)(c)(d)(e)Figure6 illustrates the distribution of the mass content of ice at the roadbed section. Due to the presence of the temperature field, ice-containing soil layers are observed at depths less than 1.2 m. At the same moment, the distribution of ice content of the roadbed section is uneven; the ice content of the west part is relatively large during freezing and melted rapidly during thawing.Figure 6
Distribution of the mass content of ice at the roadbed section on March 10, 2014 (%).Figure7 shows time history curves of the temperature field and moisture fields of the central line of the roadbed. As observed, the ice-containing soil layer at the central line has a thickness of about 1.2 m, mainly within the range 0 to 0.6 m. Moisture migration is observed in both positive and negative temperature periods, but most especially in the negative temperature period. The moisture migration is concentrated in the region at depths less than 2.5 m. Here, the moisture distribution in positive and negative temperature periods is mainly affected by self-equilibrium of the moisture and temperature fields, respectively. During freezing, the mass content of moisture in the roadbed increased with depths less than 1.2 m decreased with depths between 1.2 m and 2.5 m and remained almost constant with depths greater than 2.5 m. Here, the mass content of moisture increased near the top surface of the roadbed and within the layer at a depth of 1 m. The maximum variation of the mass content of moisture is 2%, which was about 20% higher than that before freezing. In other words, moisture migration has a significant effect. The critical depth of moisture variation is 1.2 m. At this depth, the variation range of temperature is −0.7∼0°C (with the phase transition temperature zone) during freezing. Also, a frost front is presented, resulting in the migration of moisture to the lower part of the roadbed from an unfrozen region to a frozen one.Figure 7
Trends of temperature field and moisture field of the central line. (a) Temperature (°C). (b) Mass content of ice (%). (c) Mass content of unfrozen water (%). (d) Mass content of moisture (%).
(a)(b)(c)(d)The distribution of the deformation of roadbeds is investigated based on the above analysis. In this work, the vertical displacement of the roadbed surface was the major concern. Figure8 shows the time history curve, temperature curve, frost depth curve of the roadbed deformations. Figure 9 shows vertical displacements of roadbed section at different moments. Figure 10 shows the horizontal and vertical displacement fields of the roadbed section on March 10, 2014. As it is evident, vertical displacement of the roadbed surface followed the sequence of west shoulder, east shoulder, and central line, while the vertical displacement followed the sequence of west shoulder which is higher that of the east shoulder which is higher than that at the central line. The maximum vertical displacement of the roadbed surface is approximately 11∼18 mm (15 mm is allowed according to Chinese railway roadbed design code). At the early stage of freezing, the vertical displacement increased with the frost depth. Once the frost depth exceeded a critical level, the vertical displacement did not increase further with the frost depth. The roadbed vertical displacements are asymmetrically distributed in the transverse direction and also an uneven deformation of the roadbed surface appeared. In this study, the left and right shoulders are exposed to ambient conditions, while the central line is covered by the track slab. Hence, roadbed deformations at the shoulders and the central line are uneven in the early negative temperature period (up to 14 mm) and the difference subsequently decreases.Figure 8
Time history curves of vertical displacement and frost depth.Figure 9
Vertical displacements of roadbed section at different moments.Figure 10
Distribution of roadbed section on March 10, 2014. (a) Vertical displacement field. (b) Horizontal displacement field.
(a)(b)


### 3.2.3. Discussion

For frost heaving of roadbeds in cold regions, changes of the temperature, moisture, and deformation are the cause, the medium, and the result, respectively. Roadbed deformations are macroscopic effects of hydrothermal migrations in roadbed soils.Field monitoring revealed that frost heaving of the Harbin-Dalian HSR roadbed is present from the surface to a frost depth of 1.2 m. Here, frost heaving of the well-graded broken stone dominated the frost heaving of the roadbed to a frost depth of 0.5 m [27]. Additionally, the moisture content varied drastically within the depth of 0.6 m.The results of this study demonstrated that moisture migration at the freezing front is observed in the roadbed of the Harbin-Dalian HSR under the driving force of the temperature field. Here, the moisture content increased significantly in the range between 0 m and 1.2 m, and the ice-containing soil layers are concentrated in this region especially in the range between 0 and 0.6 m. Therefore, phase transition and moisture migration in Harbin-Dalian HSR roadbed driven by an uneven distribution of the temperature field are the main causes of roadbed deformations.The transverse variation of the ground temperature of roadbeds can be attributed to the differences in the solar radiation intensity and also the surface covering layer [24, 30]. These variations of the ground temperature lead to an unbalanced development of the frost depth time and value, thus resulting in an unbalanced moisture migration. In this case, uneven transverse deformations are inevitable [30, 31]. Field monitoring data revealed differences in deformations of the east and west shoulders of the Harbin-Dalian HSR roadbed [32]. Simulation results demonstrated that the distribution of deformation field of roadbed in the transverse direction is asymmetric and uneven due to the effects of temperature and mass content of ice. Indeed, the deformation is maximized at the shoulders. The deformation of a specific part of roadbed is affected by a variation of the ground temperature and moisture at this part and of the surrounding soil. Hence, roadbed deformations are highly complicated. Frost heaving was first observed at the west shoulder, since freezing was first observed over there. However, deformation was initiated before freezing, since the freezing of the side slope occurred earlier than roadbed freezing. The uneven transverse deformation of the roadbed results in severe frost heaving at the shoulders and side slope. In fact, slope foots should be protected as weak spots due to severe ice accumulation results in horizontal displacement.The results demonstrated that a reduction of the roadbed freeze-thaw deformations should be achieved based on the analysis of thermo-hydro migration characteristics of roadbeds, instead of a thermal insulation layer on the roadbed surface. Elimination of the uneven transverse deformations should be obtained by minimizing the difference of the ground temperature in the transverse direction, which can be achieved by asymmetric transitional layout of roadbed soils [30] or by an asymmetric antifrost heaving berm [33]. Based on a theoretical study, this study provides references to the design and maintenance of HSR roadbeds in seasonally frozen regions.


## 3.2.1. Model Verification

Figure4 shows a comparison between the simulated frost depth (SFR) and the monitored frost depth (MFR) of the west shoulder and a comparison between the simulated vertical displacement (SVD) of the roadbed and the monitored vertical displacement (MVD) within 2.7 m from the adjacent roadbed section [29]. As it is evident, simulated and monitored frost depths follow consistent trends. Simulation of the deformation demonstrates a reasonable match to the monitored deformation basically. The region with frost depth less than 1.2 m is exposed to severe soil frost heaving. When the frost depth exceeded 1.2 m, frost heaving developed slowly even though the frost depth increased further. The duration of thawing subsidence was short. As shown in Figure 4, field monitoring of the adjacent roadbed section indicated 2 mm residual deformation after thawing subsidence of the roadbed [29], while previous studies claimed that roadbed deformations almost returned to the initial condition after thawing subsidence. In some cases, a small amount of settlements were observed after thawing subsidence [27].Figure 4
Simulated and monitored frost depth and vertical displacement.Therefore, from the results, it can be concluded that the proposed numerical model is effective and applicable for a comprehensive analysis of the temperature field, moisture field, and deformation field of the roadbed.


## 3.2.2. Analysis of the Simulation Results

Figure5 illustrates the temperature field distribution of the section at certain moments in time. As observed, the freezing process is unidirectional and lasted for a long period, while the thawing process is bidirectional and has a short period. Due to the sunny-shady slope effect, the temperature fields of the roadbed are asymmetrically distributed in the transverse direction during the freezing-thawing process. The maximum frost depth occurred at the west shoulder (right shoulder, sunny shoulder), followed by the east shoulder (left shoulder, shady shoulder) and that at the central line. The freezing sequence is west side slope, east side slope, and top surface. During thawing, frozen soil nuclei are observed in the roadbed and the thawing sequence is west side slope, east side slope, and roadbed top surface. Additionally, thawing of shoulders occurred earlier than that of the slope foot.Figure 5
Temperature field distribution of the section at a typical moment (°C). Temperature field distribution on (a) August 01, 2013, (b) December 10, 2013, (c) March 10, 2014, (d) April 10, 2014, and (e) August 01, 2014.
(a)(b)(c)(d)(e)Figure6 illustrates the distribution of the mass content of ice at the roadbed section. Due to the presence of the temperature field, ice-containing soil layers are observed at depths less than 1.2 m. At the same moment, the distribution of ice content of the roadbed section is uneven; the ice content of the west part is relatively large during freezing and melted rapidly during thawing.Figure 6
Distribution of the mass content of ice at the roadbed section on March 10, 2014 (%).Figure7 shows time history curves of the temperature field and moisture fields of the central line of the roadbed. As observed, the ice-containing soil layer at the central line has a thickness of about 1.2 m, mainly within the range 0 to 0.6 m. Moisture migration is observed in both positive and negative temperature periods, but most especially in the negative temperature period. The moisture migration is concentrated in the region at depths less than 2.5 m. Here, the moisture distribution in positive and negative temperature periods is mainly affected by self-equilibrium of the moisture and temperature fields, respectively. During freezing, the mass content of moisture in the roadbed increased with depths less than 1.2 m decreased with depths between 1.2 m and 2.5 m and remained almost constant with depths greater than 2.5 m. Here, the mass content of moisture increased near the top surface of the roadbed and within the layer at a depth of 1 m. The maximum variation of the mass content of moisture is 2%, which was about 20% higher than that before freezing. In other words, moisture migration has a significant effect. The critical depth of moisture variation is 1.2 m. At this depth, the variation range of temperature is −0.7∼0°C (with the phase transition temperature zone) during freezing. Also, a frost front is presented, resulting in the migration of moisture to the lower part of the roadbed from an unfrozen region to a frozen one.Figure 7
Trends of temperature field and moisture field of the central line. (a) Temperature (°C). (b) Mass content of ice (%). (c) Mass content of unfrozen water (%). (d) Mass content of moisture (%).
(a)(b)(c)(d)The distribution of the deformation of roadbeds is investigated based on the above analysis. In this work, the vertical displacement of the roadbed surface was the major concern. Figure8 shows the time history curve, temperature curve, frost depth curve of the roadbed deformations. Figure 9 shows vertical displacements of roadbed section at different moments. Figure 10 shows the horizontal and vertical displacement fields of the roadbed section on March 10, 2014. As it is evident, vertical displacement of the roadbed surface followed the sequence of west shoulder, east shoulder, and central line, while the vertical displacement followed the sequence of west shoulder which is higher that of the east shoulder which is higher than that at the central line. The maximum vertical displacement of the roadbed surface is approximately 11∼18 mm (15 mm is allowed according to Chinese railway roadbed design code). At the early stage of freezing, the vertical displacement increased with the frost depth. Once the frost depth exceeded a critical level, the vertical displacement did not increase further with the frost depth. The roadbed vertical displacements are asymmetrically distributed in the transverse direction and also an uneven deformation of the roadbed surface appeared. In this study, the left and right shoulders are exposed to ambient conditions, while the central line is covered by the track slab. Hence, roadbed deformations at the shoulders and the central line are uneven in the early negative temperature period (up to 14 mm) and the difference subsequently decreases.Figure 8
Time history curves of vertical displacement and frost depth.Figure 9
Vertical displacements of roadbed section at different moments.Figure 10
Distribution of roadbed section on March 10, 2014. (a) Vertical displacement field. (b) Horizontal displacement field.
(a)(b)


## 3.2.3. Discussion

For frost heaving of roadbeds in cold regions, changes of the temperature, moisture, and deformation are the cause, the medium, and the result, respectively. Roadbed deformations are macroscopic effects of hydrothermal migrations in roadbed soils.Field monitoring revealed that frost heaving of the Harbin-Dalian HSR roadbed is present from the surface to a frost depth of 1.2 m. Here, frost heaving of the well-graded broken stone dominated the frost heaving of the roadbed to a frost depth of 0.5 m [27]. Additionally, the moisture content varied drastically within the depth of 0.6 m.The results of this study demonstrated that moisture migration at the freezing front is observed in the roadbed of the Harbin-Dalian HSR under the driving force of the temperature field. Here, the moisture content increased significantly in the range between 0 m and 1.2 m, and the ice-containing soil layers are concentrated in this region especially in the range between 0 and 0.6 m. Therefore, phase transition and moisture migration in Harbin-Dalian HSR roadbed driven by an uneven distribution of the temperature field are the main causes of roadbed deformations.The transverse variation of the ground temperature of roadbeds can be attributed to the differences in the solar radiation intensity and also the surface covering layer [24, 30]. These variations of the ground temperature lead to an unbalanced development of the frost depth time and value, thus resulting in an unbalanced moisture migration. In this case, uneven transverse deformations are inevitable [30, 31]. Field monitoring data revealed differences in deformations of the east and west shoulders of the Harbin-Dalian HSR roadbed [32]. Simulation results demonstrated that the distribution of deformation field of roadbed in the transverse direction is asymmetric and uneven due to the effects of temperature and mass content of ice. Indeed, the deformation is maximized at the shoulders. The deformation of a specific part of roadbed is affected by a variation of the ground temperature and moisture at this part and of the surrounding soil. Hence, roadbed deformations are highly complicated. Frost heaving was first observed at the west shoulder, since freezing was first observed over there. However, deformation was initiated before freezing, since the freezing of the side slope occurred earlier than roadbed freezing. The uneven transverse deformation of the roadbed results in severe frost heaving at the shoulders and side slope. In fact, slope foots should be protected as weak spots due to severe ice accumulation results in horizontal displacement.The results demonstrated that a reduction of the roadbed freeze-thaw deformations should be achieved based on the analysis of thermo-hydro migration characteristics of roadbeds, instead of a thermal insulation layer on the roadbed surface. Elimination of the uneven transverse deformations should be obtained by minimizing the difference of the ground temperature in the transverse direction, which can be achieved by asymmetric transitional layout of roadbed soils [30] or by an asymmetric antifrost heaving berm [33]. Based on a theoretical study, this study provides references to the design and maintenance of HSR roadbeds in seasonally frozen regions.


## 4. Conclusions

With SWCC, solid-liquid ratio, and the linear thermal expansion coefficient as coupling equations, considering the influence of the transformation of ice into water and water migration on the temperature field, a thermal-hydro-mechanical model applicable for unsaturated soils was proposed. The hydrothermal characteristic parameters of frozen and unfrozen soils were expressed by step functions. Then the thermal-hydro-mechanical model is introduced in this study to explore the thermal field, water migration, and deformation variations of the HSR roadbed.From the work described in this article, it is concluded that(1)
The physical property fields of the roadbed varied significantly in soil layers with depths less than 2.5 m. Redistribution of the moisture field was triggered by a variation of the temperature field. Due to the significant effect of moisture migration, the mass content of moisture increased within the roadbed at depths greater than 1.2 m and the maximum variation of the mass content of moisture was 2 %. Most ice-containing soil layers were at depths less than 0.6 m.(2)
Affected by the boundary of the sunny-shady slopes, distributions of the temperature field, moisture field, and deformation of the roadbed in the transverse direction were asymmetric and uneven. The maximum frost depth was greatest at the west shoulder (right shoulder, shady shoulder) followed by the east shoulder (left shoulder, sunny shoulder) and least at the central line. The maximum uneven settlement between the shoulder and the central line was 16 mm, while that between the centre of the two tracks was 4 mm. Here, the maximum vertical displacements of the west shoulder, east shoulder, and central line were 18 mm, 16 mm, and 12 mm, respectively.(3)
The roadbed deformations are affected by the distribution and variation of the temperature and moisture fields. Therefore, reduction of freeze-thaw deformations and uneven deformations of roadbeds should be achieved by controlling the roadbed thermo-hydro migration characteristics.



---

*Source: 8849754-2020-09-03.xml*