# A Numerical Model of Vapour Transfer and Phase Change in Unsaturated Freezing Soils

Authors: Sihao Liang; Jidong Teng; Feng Shan; Sheng Zhang
Journal: Advances in Civil Engineering
 (2020)

---

## Abstract

In recent studies, vapour transfer is reported to lead to remarkable frost heave in unsaturated soils, but how to better model this process has not been answered. In order to avoid the great uncertainty caused by the phase change term of vapour-water-ice in the numerical iteration process, a new numerical model is developed based on the coupled thermal and hydrological processes. The new model avoids using the local equilibrium assumption and the hydraulic relations that accounts for liquid water flow, which provides a new way for the water-heat coupling movement problem. The model is established by using COMSOL Multiphysics, which is a multiphysics simulation software through finite element analysis. The model is evaluated by comparing simulated results with data from column freezing experiments for unsaturated coarse-grained soils. Simulated values of the total water content compare well with experimental values. The model is proved to be applicable and numerically stable for a high-speed railway subgrade involving simultaneous heat and moisture transport. An agreement can be found between the predicted and measured frost/thawed depth and soil moisture profiles, demonstrating that the model is able to simulate rapidly changing boundary conditions and nonlinear water content profiles in the soil.


---

## Body

## 1. Introduction

Many engineering problems, including cracking of pavements, damage to the foundation of structures, and fracture of pipelines, are caused by the freezing and thawing process in cold regions. [1–3]. Previous studies mostly focused on the problem of frost heave caused by the liquid water transfer from warm to cold side in freezing soils [4–7]. In recent studies, it is shown that vapour transfer can result in ice accumulation in unsaturated freezing soil and, thus, cause frost damage to infrastructures. For example, Eigenbrod and Kennepohl [8] contended that vapour flow led to water accumulation at the base of pavements in the northern parts of Canada and the United States. Severe frost heave problems were reported to occur at an airport in Northwestern China with 20 m-deep groundwater table and limited annual rainfall. The specific phenomenon was named as the canopy effect by Zhang et al. [9]. The researchers attributed the water accumulation beneath the impervious cover to the vapour transfer in unsaturated freezing soils. Niu et al. [10] performed field monitoring of soil temperature, water content, and deformation in the railway embankment and concluded that water vapour diffusion might be an essential influence factor to the formation of near-surface ice.The study on frozen soil models began in the 1970s, which aims to solve the coupled heat and mass conservative equation, for example, the works of Harlan [11] and Talyor and Luthin [12]. In the earlier studies, vapour flow is usually neglected comparing to the liquid water transfer in freezing soil. This assumption is made on the basis of the cognition that (1) neglecting the vapour flow could largely simplify the numerical computation and (2) it is truly insignificant for the vapour flow in the soil of having a continuous liquid phase. Considering that more evidence has revealed the complexity in soil freezing process, a large number of sophisticated models or parameterizations have been presented during the following several decades [13–17]. Some of these studies take the contribution of vapour flow and its phase change into account, while others choose to neglect yet. The proposed models have laid a good foundation for understanding the soil freezing process.In the most popular models for frozen soils, it is common to set up two governing equations with four undetermined variables, i.e., soil temperature, matric suction of soil, liquid water content, and ice content. In order to numerically solve the governing equations, another two relations are required. The one is a relationship between matric suction and liquid water content called the Soil Water Characteristic Curve (SWCC), such as the VG model [18] and BC model [19]. The other is a formula that relates liquid water content and ice content to temperature, which is referred to as the Soil Freezing Characteristic Curve (SFCC) [20, 21]. The differences in the existing models are the form of governing equations and the functions for the SWCC and SFCC. It is noted that the studies that take vapour flow into account usually relate soil vapour content (or vapour density) to temperate and matric suction on the basis of a local equilibrium theory, which was initially proposed by Philip and de Vries [22]. But, here, the characterising parameter for vapour flow is a nonindependent variable.When applying the pervious models to analyse the frost heave in unsaturated soils, it is found they are incapable to reveal the mechanism of vapour flow, due to the following reasons. Firstly, when soil water content is relatively low, specifically near to the residual water content, it is inaccurate to describe liquid water flow by using the SWCC and the derived unsaturated permeability function [23, 24]. Besides, liquid water flow may be smaller several magnitude orders than vapour flow at this stage. The vapour flow may dominate the mass transfer when the soil is relatively dry [25]. Secondly, the theory proposed by Philip and de Vries [22] assumes that vapour concentration in air is always in equilibrium with liquid water, i.e., the vapour concentration is determined by the curved liquid/vapour interface. But, the recent study suggests that the equilibrium assumption works at the low water content in soil because the equilibrium establishment is not instantaneous at that time [26–28]. Moreover, when the soil temperature reaches to the freezing point, the formation of ice phase in the soil will change the liquid/vapour interface. Thus, the local equilibrium theory would be challenged in case of being used to describe the moisture transfer in frozen soil. Thirdly, the mechanism of multiphase transfer in unsaturated freezing soils has been less studied. They remain unknown, for example, the phase changes among pore ice, vapour, and liquid water and the influence of ice phase on vapour and liquid water flow [29]. The existing model introduced a few empirical equations in their models to clarify the role of ice phase in unsaturated soil while these fitting parameters are hard to be determined by laboratory experiment. Therefore, although many different simulations are dealing with freezing soil issue, few works have been reported to model the vapour flow in unsaturated freezing soil.Aiming to better understand the mechanism of vapour flow and phase change in unsaturated freezing soils, this study presents a theoretical framework to formulate the coupled thermal and hydrological process. In the new model, vapour flow governs the mass transfer process, and liquid water flow is neglective. The coupled model is solved numerically by the finite element method. A series of laboratory test results are used to validate the numerical formulations/codes. Finally, some conclusions are drawn based on the results and discussion.


## 2. Mathematical Model

Figure1 depicts the physical system which is a vertical, one-dimensional soil profile covered by an impermeable plate. The plate represents the infrastructures that covered on the soil surface. The soil is considered as a homogeneous porous medium, where vapour diffuses from the warm and moisture side to the cold and dry side. Ice lens is formed beneath the impervious cover where the soil temperature is below the 0 °C due to vapour-ice desublimation. It is noted that some experimental research has reported the generation of the ice lens zone [8, 9, 30]; however, less attention has been paid from the theoretical approach to model this phenomenon. As shown in Figure 1, beyond the ice lens zone, the vapour flow is governed by the vapour concentration in soil and will condense into liquid water. The condensed water will continue freezing into pore ice in the soil where the temperature drops below 0°C. The transient profiles of soil temperature (solid line in Figure 1) and total water saturation (dashed line) are the outcomes of this problem.Figure 1
(a) Schematic diagram of the vapour-ice desublimation process in soil; (b) the transient profiles of water content and soil temperature.
(a)(b)Some assumptions are made to simplify the quantitative description of the coupled heat and moisture process, as follows:(a)
The deformation of the soil matrix due to variations of temperature and pore water pressure or ice formation can be neglected, which is consistent with the assumption in the work of Sheng et al. [3] and Teng et al. [31].(b)
The freezing front is the depth where the soil temperature is 0°C.(c)
The ice lens zone is the region where the soil is saturated. The soil cannot be oversaturated, i.e., the total amount of water in the ice lens zone can be represented by the thickness of the ice lens zone, which can be supported by the experimental observation in the work of Teng et al. [32].


### 2.1. Governing Equations for the Vapour Transfer Zone

The one-dimensional governing equation for vapour flow in the unsaturated freezing soil is given by the following mass conservation equation:(1)∂n1−Sρv∂t=−∂qv∂z−∂nSρw∂t,where z(m) and t(s) represent the position and time, respectively, n is the porosity of soil (unitless), and ρ is the density in soil (kg m−3). The subscripts w, v, and i denote the liquid water, vapour, and pore ice, respectively. qv is vapour flux (kg m−2 s−1), and S is the total saturation of liquid water and pore ice (unitless); here, S = Sw + Siρi/ρw.The vapour movement in unsaturated freezing soil is recognised to be driven by both the pressure gradient and vapour concertation gradient [33]. The vapour mass flux in this region according to Darcy’s law and Fick’s law is(2)qv=−ρvKμv∂P∂z−1−SnτD∂ρv∂z,where K is the air permeability of soil (m2), μv is the dynamic viscosity of vapour (kg m−1 s−1), P is the vapour pressure (kPa), D is the vapour diffusivity in soil (m2 s−1), and τ is the tortuosity factor (unitless). The parameter τ is the ratio of real length of transfer path to apparent length, which is assigned to be 1.2 in this study [34].The vapour phase in unsaturated freezing soil can be considered as ideal gas. It follows the ideal gas law,P = ρvRT. Substituting this equation into equation (2) leads to(3a)qv=−KT∂T∂z−Ks∂S∂z,(3b)KT=KμvRρv2,(3c)KS=KμvRρvT+1−SnτD∂ρv∂S,where KS and KT are the effective diffusivity by the saturation gradient (kg m−1 s−1) and effective mass conductivity by the temperature gradient (kg m−1 s−1 K−1), respectively. R is the specific gas constant of water vapour (461.89 J kg−1 K−1).The mass transfer equation can be obtained by substituting equation (3a) into equation (1):(4)∂n1−Sρv∂t+∂nSρw∂t=∂∂zKT∂T∂z+Ks∂S∂z.It has been revealed that the heat conduction is two orders of magnitude larger than the sensible heat of vapour flow [35]. Thus, the heat conservative equation neglects the term of the sensible heat of vapour flow in the soil and is written as(5a)ρC∂T∂t=∂∂zλ+KTLi,v∂T∂z+∂∂zKSLi,v∂S∂z−Li,w∂nSwρw∂t,(5b)ρC=1−nρsCs+nSiρiCi+nSwρwCw+n1−SρvCv+n1−SLi,v∂ρv∂T,(5c)λ=1−nλs+nSiλi+nSwλw+n1−Sλv,where ρ is the total density of soil and C is the volumetric heat capacity and is given in equation (5b). λ is the thermal conductivity (W m−1 K−1) that is defined in equation (5c). L donates the latent heat between different phases. The subscripts s in the abovementioned equations indicate the solid phase of soil particles.The saturation of the unfrozen water is an empirical function of temperature in the unsaturated freezing soil. Anderson and Tice [36] and Anderson and Morgenstern [37] found a reasonable approximate equation with a power law:(6)Suw=ρsρw1−n100nαT0−Tβ,T≤T0,where T0 is 273.15 K, the freezing point of free water, and α and β are empirical fitting parameters that relate to the specific surface area of soil. Anderson and Tice [36]; Anderson and Morgenstern [37]; and Blanchard and Fremond [38] tested the empirical formulas of parameters α and β for different types of soils and listed their recommended values. Although the power relationship as equation (6) has been proposed several decades ago, it is still commonly given in reference texts as a valid approach to estimating the unfrozen water content in freezing soils [39, 40].Suw in equation (6) presents the maximum liquid water content that will not freeze at a subzero temperature T. Based on the computed values of S and Suw, the criteria for determining the saturation of pore ice and liquid water can be obtained [7], as follows:(7a)Si=0,T>T0,0,T≤T0andS≤Suw,ρwS−SuwρiT≤T0andS≤Suw,(7b)Sw=S−Siρiρw.In the hygroscopic porous material, the three phases of water, i.e., unfrozen liquid water, ice, and vapour coexist in the freezing material pores. The presence of unfrozen liquid water complicates the process of phase change. In this case, the thermodynamic equilibrium relationship (the Clapeyron equation) fails to describe the truth because the Clapeyron equation is a way only for characterising the discontinuous phase transition between two phases of matter of a single constituent [41, 42]. It has been recognised from chemical engineering that an adsorption-desorption equilibrium relationship should be replaced to describe the desublimation or sublimation process of a hygroscopic porous media with bound moisture [43, 44]. Such a relation is used to describe the equilibrium state among vapour, liquid water, and pore ice in unsaturated freezing soil. Wang et al. [45] tested serval kinds of adsorption-desorption equilibrium relations in the literature and found that Kelvin’s style in the exponential form proposed by Rajniak and Yang [46] could generate a better performance. The equation is expressed as(8)PT,SPvsT=e−γ1−S/SRT,where Pvs(T) is the saturated vapour pressure at temperature T (Pa), given in Table 1. P(T, S) is the vapour pressure in unsaturated freezing soil (Pa), which is a function of temperature T and saturation S. γ is an empirical parameter, a suggested value is 5000 according to Wang et al. [45].Table 1
The parameters used in computation.
SymbolValue and equationsReferenceρw (kg m−3)1000[47]ρi (kg m−3)913[48]Ci (J kg−1K−1)1930[48]Cs (J kg−1K−1)800[49]Cv (J kg−1K−1)1886[48]Cw (J kg−1K−1)4180[25]λi (W m−1K−1)2.22[48]λs (W m−1K−1)2.68[25]λv (W m−1K−1)0.022[48]λw (W m−1K−1)0.54[25]μv (kg m−1 s−1)0.011 × (T/273.15)1.5/(T + 961)[50]Li,v (J kg −1)2.839 × 106[51]Li,w (J kg −1)0.334 × 106[15]Pvs(T) (Pa)lnPvs(T) = 9.55 − 5723.265/T + 3.53ln(T) − 0.007283 T[52]Da (m2 s−1)Da = 2.12 × 10−5(T/273.15)2[49]The simulated model consists of mass conservation equation (4) and energy balance equation (5), which are highly nonlinear and coupled. There are four unknown variables in the two equations, Sw, Si, T, and ρv. Equations (6)–(8) provide the other two additional relations by adding the saturation of unfrozen water Suw, such that the simultaneous equations can be solved mathematically.


### 2.2. Governing Equations for the Ice Lens Zone

It is assumed that there is no moisture movement in the ice lens zone, i.e., the density of the frost layer is uniform. The mass balance at this region can be expressed by(9)ρ=ρi,0<z<zfst,where zfs (t) represents the position of the frost surface of the ice lens, which is a function of time.The governing equation for heat transfer in the ice lens zone can be written as(10)ρiCi∂T∂t=λi∂2T∂z2,0<z<zfst,where Ci is the specific heat of ice lens (J kg−1 K−1) and λi is the effective thermal conductivity of the frost layer (W m−1 K−1).


### 2.3. Heat and Mass Balance at the Frost Surface

The heat and moisture mass balance at the moving front can be expressed by(11)λi∂T∂zz=zfs−−λ∂T∂zz=zfs+=−Li,vqvz=zfs+.The water vapour on the frost surface is assumed to be saturated.(12)ρv=ρvs=PvsTfsRTfs,atz=zfs.Two new parameters are added into the model,zfs at the frost surface and T at the ice lens zone, and ρvs is the saturated vapour density, which is a function of temperature [53].


### 2.4. Soil Hydraulic Properties

In complete dry soil, the air permeabilityK is equal to inherent permeability ks. As for the unsaturated freezing soil, an impedance factor is added to the air permeability function [9, 12, 19].(13)KS=10−cSiks1−S21−S2+b/b,where b and c are empirical parameters and the impedance factor 10−cSi denotes the obstruction of pore ice to vapour flows. When the soil temperature is greater than the freezing point, Si becomes 0, and equation (13) reverts to the formula as defined by Brooks and Corey [19].The vapour transfer in the porous medium described by Fick’s law can be divided into two types, molecular diffusion and Knudsen diffusion. Molecular diffusion is caused by the collisions with the gas molecules, while Knudsen diffusion is caused by the collisions between gas molecules and the wall of the transport channel. Molecular diffusion dominates the process when the mean free path of the molecules is smaller than the average pore radius. Otherwise, Knudsen diffusion does [54]. The harmonic averaging method is adopted here to compute the vapour diffusivity D based on the molecular diffusion coefficient Dm and the Knudsen diffusion coefficient Dk [55]:(14)D=11/Dm+1/Dk,where Dm can be written as follows [56]:(15)Dm=0.66n1−SDa,where Da is the diffusivity of water vapour in air (m2 s−1), as given in Table 1. The expression of Dk was given by Geankoplis [57] as follows:(16)Dk=97.0d2Tmw,where d is the averaged soil pore diameter (m) and mw is the molecular weight of water (kg mol−1). The values of related parameters are presented in Table 1.


### 2.5. Numerical Implementation

Equations (4), (5a), (5b), and (5c) are highly nonlinear since primary variables vary with time and primary variables. The finite element method for spatial discretization and the finite difference method for temporal discretization are used to numerically solve the governing equations, which is performed in the COMSOL Multiphysics package. A detailed description for the numerical approach can be found in the work of Teng et al. [32].In order to achieve a numerical stability of the solution, the following strategies are adopted: the domain is divided into 10000 elements, and the total time duration is divided into 12000 steps.


## 2.1. Governing Equations for the Vapour Transfer Zone

The one-dimensional governing equation for vapour flow in the unsaturated freezing soil is given by the following mass conservation equation:(1)∂n1−Sρv∂t=−∂qv∂z−∂nSρw∂t,where z(m) and t(s) represent the position and time, respectively, n is the porosity of soil (unitless), and ρ is the density in soil (kg m−3). The subscripts w, v, and i denote the liquid water, vapour, and pore ice, respectively. qv is vapour flux (kg m−2 s−1), and S is the total saturation of liquid water and pore ice (unitless); here, S = Sw + Siρi/ρw.The vapour movement in unsaturated freezing soil is recognised to be driven by both the pressure gradient and vapour concertation gradient [33]. The vapour mass flux in this region according to Darcy’s law and Fick’s law is(2)qv=−ρvKμv∂P∂z−1−SnτD∂ρv∂z,where K is the air permeability of soil (m2), μv is the dynamic viscosity of vapour (kg m−1 s−1), P is the vapour pressure (kPa), D is the vapour diffusivity in soil (m2 s−1), and τ is the tortuosity factor (unitless). The parameter τ is the ratio of real length of transfer path to apparent length, which is assigned to be 1.2 in this study [34].The vapour phase in unsaturated freezing soil can be considered as ideal gas. It follows the ideal gas law,P = ρvRT. Substituting this equation into equation (2) leads to(3a)qv=−KT∂T∂z−Ks∂S∂z,(3b)KT=KμvRρv2,(3c)KS=KμvRρvT+1−SnτD∂ρv∂S,where KS and KT are the effective diffusivity by the saturation gradient (kg m−1 s−1) and effective mass conductivity by the temperature gradient (kg m−1 s−1 K−1), respectively. R is the specific gas constant of water vapour (461.89 J kg−1 K−1).The mass transfer equation can be obtained by substituting equation (3a) into equation (1):(4)∂n1−Sρv∂t+∂nSρw∂t=∂∂zKT∂T∂z+Ks∂S∂z.It has been revealed that the heat conduction is two orders of magnitude larger than the sensible heat of vapour flow [35]. Thus, the heat conservative equation neglects the term of the sensible heat of vapour flow in the soil and is written as(5a)ρC∂T∂t=∂∂zλ+KTLi,v∂T∂z+∂∂zKSLi,v∂S∂z−Li,w∂nSwρw∂t,(5b)ρC=1−nρsCs+nSiρiCi+nSwρwCw+n1−SρvCv+n1−SLi,v∂ρv∂T,(5c)λ=1−nλs+nSiλi+nSwλw+n1−Sλv,where ρ is the total density of soil and C is the volumetric heat capacity and is given in equation (5b). λ is the thermal conductivity (W m−1 K−1) that is defined in equation (5c). L donates the latent heat between different phases. The subscripts s in the abovementioned equations indicate the solid phase of soil particles.The saturation of the unfrozen water is an empirical function of temperature in the unsaturated freezing soil. Anderson and Tice [36] and Anderson and Morgenstern [37] found a reasonable approximate equation with a power law:(6)Suw=ρsρw1−n100nαT0−Tβ,T≤T0,where T0 is 273.15 K, the freezing point of free water, and α and β are empirical fitting parameters that relate to the specific surface area of soil. Anderson and Tice [36]; Anderson and Morgenstern [37]; and Blanchard and Fremond [38] tested the empirical formulas of parameters α and β for different types of soils and listed their recommended values. Although the power relationship as equation (6) has been proposed several decades ago, it is still commonly given in reference texts as a valid approach to estimating the unfrozen water content in freezing soils [39, 40].Suw in equation (6) presents the maximum liquid water content that will not freeze at a subzero temperature T. Based on the computed values of S and Suw, the criteria for determining the saturation of pore ice and liquid water can be obtained [7], as follows:(7a)Si=0,T>T0,0,T≤T0andS≤Suw,ρwS−SuwρiT≤T0andS≤Suw,(7b)Sw=S−Siρiρw.In the hygroscopic porous material, the three phases of water, i.e., unfrozen liquid water, ice, and vapour coexist in the freezing material pores. The presence of unfrozen liquid water complicates the process of phase change. In this case, the thermodynamic equilibrium relationship (the Clapeyron equation) fails to describe the truth because the Clapeyron equation is a way only for characterising the discontinuous phase transition between two phases of matter of a single constituent [41, 42]. It has been recognised from chemical engineering that an adsorption-desorption equilibrium relationship should be replaced to describe the desublimation or sublimation process of a hygroscopic porous media with bound moisture [43, 44]. Such a relation is used to describe the equilibrium state among vapour, liquid water, and pore ice in unsaturated freezing soil. Wang et al. [45] tested serval kinds of adsorption-desorption equilibrium relations in the literature and found that Kelvin’s style in the exponential form proposed by Rajniak and Yang [46] could generate a better performance. The equation is expressed as(8)PT,SPvsT=e−γ1−S/SRT,where Pvs(T) is the saturated vapour pressure at temperature T (Pa), given in Table 1. P(T, S) is the vapour pressure in unsaturated freezing soil (Pa), which is a function of temperature T and saturation S. γ is an empirical parameter, a suggested value is 5000 according to Wang et al. [45].Table 1
The parameters used in computation.
SymbolValue and equationsReferenceρw (kg m−3)1000[47]ρi (kg m−3)913[48]Ci (J kg−1K−1)1930[48]Cs (J kg−1K−1)800[49]Cv (J kg−1K−1)1886[48]Cw (J kg−1K−1)4180[25]λi (W m−1K−1)2.22[48]λs (W m−1K−1)2.68[25]λv (W m−1K−1)0.022[48]λw (W m−1K−1)0.54[25]μv (kg m−1 s−1)0.011 × (T/273.15)1.5/(T + 961)[50]Li,v (J kg −1)2.839 × 106[51]Li,w (J kg −1)0.334 × 106[15]Pvs(T) (Pa)lnPvs(T) = 9.55 − 5723.265/T + 3.53ln(T) − 0.007283 T[52]Da (m2 s−1)Da = 2.12 × 10−5(T/273.15)2[49]The simulated model consists of mass conservation equation (4) and energy balance equation (5), which are highly nonlinear and coupled. There are four unknown variables in the two equations, Sw, Si, T, and ρv. Equations (6)–(8) provide the other two additional relations by adding the saturation of unfrozen water Suw, such that the simultaneous equations can be solved mathematically.


## 2.2. Governing Equations for the Ice Lens Zone

It is assumed that there is no moisture movement in the ice lens zone, i.e., the density of the frost layer is uniform. The mass balance at this region can be expressed by(9)ρ=ρi,0<z<zfst,where zfs (t) represents the position of the frost surface of the ice lens, which is a function of time.The governing equation for heat transfer in the ice lens zone can be written as(10)ρiCi∂T∂t=λi∂2T∂z2,0<z<zfst,where Ci is the specific heat of ice lens (J kg−1 K−1) and λi is the effective thermal conductivity of the frost layer (W m−1 K−1).


## 2.3. Heat and Mass Balance at the Frost Surface

The heat and moisture mass balance at the moving front can be expressed by(11)λi∂T∂zz=zfs−−λ∂T∂zz=zfs+=−Li,vqvz=zfs+.The water vapour on the frost surface is assumed to be saturated.(12)ρv=ρvs=PvsTfsRTfs,atz=zfs.Two new parameters are added into the model,zfs at the frost surface and T at the ice lens zone, and ρvs is the saturated vapour density, which is a function of temperature [53].


## 2.4. Soil Hydraulic Properties

In complete dry soil, the air permeabilityK is equal to inherent permeability ks. As for the unsaturated freezing soil, an impedance factor is added to the air permeability function [9, 12, 19].(13)KS=10−cSiks1−S21−S2+b/b,where b and c are empirical parameters and the impedance factor 10−cSi denotes the obstruction of pore ice to vapour flows. When the soil temperature is greater than the freezing point, Si becomes 0, and equation (13) reverts to the formula as defined by Brooks and Corey [19].The vapour transfer in the porous medium described by Fick’s law can be divided into two types, molecular diffusion and Knudsen diffusion. Molecular diffusion is caused by the collisions with the gas molecules, while Knudsen diffusion is caused by the collisions between gas molecules and the wall of the transport channel. Molecular diffusion dominates the process when the mean free path of the molecules is smaller than the average pore radius. Otherwise, Knudsen diffusion does [54]. The harmonic averaging method is adopted here to compute the vapour diffusivity D based on the molecular diffusion coefficient Dm and the Knudsen diffusion coefficient Dk [55]:(14)D=11/Dm+1/Dk,where Dm can be written as follows [56]:(15)Dm=0.66n1−SDa,where Da is the diffusivity of water vapour in air (m2 s−1), as given in Table 1. The expression of Dk was given by Geankoplis [57] as follows:(16)Dk=97.0d2Tmw,where d is the averaged soil pore diameter (m) and mw is the molecular weight of water (kg mol−1). The values of related parameters are presented in Table 1.


## 2.5. Numerical Implementation

Equations (4), (5a), (5b), and (5c) are highly nonlinear since primary variables vary with time and primary variables. The finite element method for spatial discretization and the finite difference method for temporal discretization are used to numerically solve the governing equations, which is performed in the COMSOL Multiphysics package. A detailed description for the numerical approach can be found in the work of Teng et al. [32].In order to achieve a numerical stability of the solution, the following strategies are adopted: the domain is divided into 10000 elements, and the total time duration is divided into 12000 steps.


## 3. Application to a Laboratory Freezing Experiment

### 3.1. Testing Program

A series of laboratory tests were carried out on the basis of a specifically designed new device. A detailed illustration can be found in the work of Teng et al. [32]. The temperatures of the top and bottom ends of the specimen are controlled accurately by a pair of plates. The top end of the specimen is sealed to prevent moisture source or sink. A Mariotte bottle is used to supply distilled water, leaving a void approximately 1 cm in height between the bottom end of the specimen and water surface, such that only vapour can enter the specimen.A kind of silica sand is applied in this study, which has a particle size range from 0.5 mm to 1 mm. The physical and hydraulic properties of the sample are presented in Table2. The measured Soil Freezing Characteristic Curve is shown in Figure 2, where the curve is fitted by equation (6). The soil specimen is placed into the cylinder with a controlled dry density of 1.40 g/cm3. The final height of the sample is 13.5 cm. Three Time-Domain-Reflectometer transducers and seven thermistors are inserted into the cylinder to measure the transit temperatures and water contents. In addition, the specimen is divided into 1 cm-high columns to check water contents at different depths after the completion of the test.Table 2
Physical and hydraulic properties of the sample.
SampleSpecific gravityMaximum dry densityInherent permeabilityAverage pore diameter (m)Fitting parameterbFitting parameter cSilica sand2.63 g/cm31.60 g/cm33.47 × 10−10 m20.000753.708.00Figure 2
The simulated and measured SFCC of the sample. The best-fit results forα and β in equation (4) are 30.880 and −1.855, respectively.Five experimental cases are designed as listed in Table3. In all the 5 cases, the top end of the sample is subjected to a subfreezing temperature, while the bottom end is subjected to a superfreezing different temperature. Case 1 and case 2 correspond to two different initial water contents. Case 3 to case 5 change the test period in order to study the transient total water content profile in time.Table 3
Test conditions.
CaseInitial saturation (%)Bottom temperature (°C)Top temperature (°C)Water supply pattern and test period131.8110−10Control test, vapour, 7d2010−10Vapour, 7d331.8110−10Vapour, 1d431.8110−10Vapour, 3d531.8110−10Vapour, 5d


### 3.2. Simulation Result

The computed and measured water content profiles for cases 1 and 2 are shown in Figure3. As for the result of case 1, it shows that the peak water content appears at the top surface and the freezing front. A large amount of ice is accumulated at the top surface, which is the ice lens zone. The simulated result has a good general agreement with the measured data and the simulated saturation of ice and water at the top surface, and the freezing front is quite close to the measured result. A minor disagreement can be observed in the unfrozen zone, which may attribute to the deviation of adsorption-desorption equilibrium relations. The simulated result has a quite good match with the tested data, where only one peak value at the top surface can be observed.Figure 3
Comparison between the measured and simulated saturation of water and ice for cases 1 and 2.Figure4 shows a comparison between the measured and computed result for cases 3, 4, and 5. In both the simulated and measured results, the peak values can be observed at both the top surface and the freezing front where the ice content increases with the test period. At the same time, the water content in the unfrozen zone seems to be independent of the test period. It is apparent that liquid water goes downward within the first day of freezing due to gravity, which implies that the formation of ice is mainly induced by vapour flow.Figure 4
Comparison between the measured and simulated saturation of water and ice for cases 3, 4, and 5. The symbols are the measured result, while the solid lines are the numerical simulation.


## 3.1. Testing Program

A series of laboratory tests were carried out on the basis of a specifically designed new device. A detailed illustration can be found in the work of Teng et al. [32]. The temperatures of the top and bottom ends of the specimen are controlled accurately by a pair of plates. The top end of the specimen is sealed to prevent moisture source or sink. A Mariotte bottle is used to supply distilled water, leaving a void approximately 1 cm in height between the bottom end of the specimen and water surface, such that only vapour can enter the specimen.A kind of silica sand is applied in this study, which has a particle size range from 0.5 mm to 1 mm. The physical and hydraulic properties of the sample are presented in Table2. The measured Soil Freezing Characteristic Curve is shown in Figure 2, where the curve is fitted by equation (6). The soil specimen is placed into the cylinder with a controlled dry density of 1.40 g/cm3. The final height of the sample is 13.5 cm. Three Time-Domain-Reflectometer transducers and seven thermistors are inserted into the cylinder to measure the transit temperatures and water contents. In addition, the specimen is divided into 1 cm-high columns to check water contents at different depths after the completion of the test.Table 2
Physical and hydraulic properties of the sample.
SampleSpecific gravityMaximum dry densityInherent permeabilityAverage pore diameter (m)Fitting parameterbFitting parameter cSilica sand2.63 g/cm31.60 g/cm33.47 × 10−10 m20.000753.708.00Figure 2
The simulated and measured SFCC of the sample. The best-fit results forα and β in equation (4) are 30.880 and −1.855, respectively.Five experimental cases are designed as listed in Table3. In all the 5 cases, the top end of the sample is subjected to a subfreezing temperature, while the bottom end is subjected to a superfreezing different temperature. Case 1 and case 2 correspond to two different initial water contents. Case 3 to case 5 change the test period in order to study the transient total water content profile in time.Table 3
Test conditions.
CaseInitial saturation (%)Bottom temperature (°C)Top temperature (°C)Water supply pattern and test period131.8110−10Control test, vapour, 7d2010−10Vapour, 7d331.8110−10Vapour, 1d431.8110−10Vapour, 3d531.8110−10Vapour, 5d


## 3.2. Simulation Result

The computed and measured water content profiles for cases 1 and 2 are shown in Figure3. As for the result of case 1, it shows that the peak water content appears at the top surface and the freezing front. A large amount of ice is accumulated at the top surface, which is the ice lens zone. The simulated result has a good general agreement with the measured data and the simulated saturation of ice and water at the top surface, and the freezing front is quite close to the measured result. A minor disagreement can be observed in the unfrozen zone, which may attribute to the deviation of adsorption-desorption equilibrium relations. The simulated result has a quite good match with the tested data, where only one peak value at the top surface can be observed.Figure 3
Comparison between the measured and simulated saturation of water and ice for cases 1 and 2.Figure4 shows a comparison between the measured and computed result for cases 3, 4, and 5. In both the simulated and measured results, the peak values can be observed at both the top surface and the freezing front where the ice content increases with the test period. At the same time, the water content in the unfrozen zone seems to be independent of the test period. It is apparent that liquid water goes downward within the first day of freezing due to gravity, which implies that the formation of ice is mainly induced by vapour flow.Figure 4
Comparison between the measured and simulated saturation of water and ice for cases 3, 4, and 5. The symbols are the measured result, while the solid lines are the numerical simulation.


## 4. Application to a High-Speed Railway Subgrade

### 4.1. Problem Description

In order to evaluate the numerical stability of the proposed model and to illustrate its application to high-speed railway subgrade, this section will perform a long-term test of the model by Niu et al. [10] who performed in situ tests to measure the ground temperatures, moisture content, and frost heave of the subgrade of the Harbin-Dalian Passenger Dedicated Line. The subgrade at one site (K977) is built on the undisturbed ground surface, while the other is in a cut section (K1004). The measured data start from Nov 1, 2013, to Oct 17, 2014, that last for 351 days. The upper layer of the subgrade structure is a well-graded gravel with cement (thickness is 70 cm). The materials in the second layer are group A/B fills (thickness is 230 cm). The lower layer is a well-graded crushed stone and sand (thickness is 50 cm). In order to simplify the numerical simulation, the subgrade is considered as a 300 cm-high embankment with the 70 cm cemented gravel at the top position, and the geomembrane is not considered.The main differences between the in situ measurement and the laboratory column experiment are that the soil is subject to freezing, as well as thawing, and the temperature variation in situ is much higher. The measured near-surface temperature is shown in Figure5, where a fitting function is obtained as the input upper boundary. The lower temperature of the soil specimen is kept constant 3°C. It is assumed that there is no flow at both the upper and lower boundaries. The thermal properties of soil can be found in Table 1. The other parameters are the same as those for silica sand, which can be found in Zhang et al. [9]. The hydraulic properties of the soil are the same as those in the work of Teng et al. [32].Figure 5
The measured temperature for the K977 site and K1004 site. The solid lines are the best-fit function, which is the input of the upper boundary.


### 4.2. Freezing-Thawing Simulations

Figure6 presents the predicted frost and thawed depth variation in the 160 days, which indicates one freezing and thawing cycle. The frost and thawed depth are usually defined as the depth where the soil temperature is zero. Thus, the predicted frost and thawed depth are obtained from the simulated transient temperature profiles. It can be observed that the frost depths increase roughly linearly and reach a peak at about 90 days and 80 days, respectively, for the site of K977 and K1004, and then, the depth gradually decreases with time. The thawed depth shows a linear increase tendency. A quite good agreement between the measured and predicted frost and thawed depth can be observed, which indicates that the proposed numerical code performs well in predicting the temperature field.Figure 6
The predicted and measured frozen/thawed depth versus time.Figure7 presents measured and predicted water content profiles in the subgrades of the two sites during the onset of freezing conditions and rapid warming period. It can be observed that the water content near the top position of the subgrade at both sites decreases rapidly as the freezing condition develops, which is caused by the ice accumulation as the freezing front progresses downward and the upward migration of vapour into the frozen zone. The predicted result can well capture the features of moisture profile. During the thawing period, the moisture in soil drains downward and the soil moisture becomes lower comparing to the freezing condition. The position of the peak moisture content during the thawing condition is relatively lower than that during the freezing condition. The proposed result can agree reasonably with the measured data. It is emphasized that a subgrade of the layered structure has much lower water contents in the upper part. The measured distributions may be different from the hypothetical condition. This case is hard to be modelled in a numerical code [15], while the proposed model can reasonably reproduce the tendency.Figure 7
The predicted and measured volumetric water content profiles at K977 and K1004 sites.
(a)(b)


## 4.1. Problem Description

In order to evaluate the numerical stability of the proposed model and to illustrate its application to high-speed railway subgrade, this section will perform a long-term test of the model by Niu et al. [10] who performed in situ tests to measure the ground temperatures, moisture content, and frost heave of the subgrade of the Harbin-Dalian Passenger Dedicated Line. The subgrade at one site (K977) is built on the undisturbed ground surface, while the other is in a cut section (K1004). The measured data start from Nov 1, 2013, to Oct 17, 2014, that last for 351 days. The upper layer of the subgrade structure is a well-graded gravel with cement (thickness is 70 cm). The materials in the second layer are group A/B fills (thickness is 230 cm). The lower layer is a well-graded crushed stone and sand (thickness is 50 cm). In order to simplify the numerical simulation, the subgrade is considered as a 300 cm-high embankment with the 70 cm cemented gravel at the top position, and the geomembrane is not considered.The main differences between the in situ measurement and the laboratory column experiment are that the soil is subject to freezing, as well as thawing, and the temperature variation in situ is much higher. The measured near-surface temperature is shown in Figure5, where a fitting function is obtained as the input upper boundary. The lower temperature of the soil specimen is kept constant 3°C. It is assumed that there is no flow at both the upper and lower boundaries. The thermal properties of soil can be found in Table 1. The other parameters are the same as those for silica sand, which can be found in Zhang et al. [9]. The hydraulic properties of the soil are the same as those in the work of Teng et al. [32].Figure 5
The measured temperature for the K977 site and K1004 site. The solid lines are the best-fit function, which is the input of the upper boundary.


## 4.2. Freezing-Thawing Simulations

Figure6 presents the predicted frost and thawed depth variation in the 160 days, which indicates one freezing and thawing cycle. The frost and thawed depth are usually defined as the depth where the soil temperature is zero. Thus, the predicted frost and thawed depth are obtained from the simulated transient temperature profiles. It can be observed that the frost depths increase roughly linearly and reach a peak at about 90 days and 80 days, respectively, for the site of K977 and K1004, and then, the depth gradually decreases with time. The thawed depth shows a linear increase tendency. A quite good agreement between the measured and predicted frost and thawed depth can be observed, which indicates that the proposed numerical code performs well in predicting the temperature field.Figure 6
The predicted and measured frozen/thawed depth versus time.Figure7 presents measured and predicted water content profiles in the subgrades of the two sites during the onset of freezing conditions and rapid warming period. It can be observed that the water content near the top position of the subgrade at both sites decreases rapidly as the freezing condition develops, which is caused by the ice accumulation as the freezing front progresses downward and the upward migration of vapour into the frozen zone. The predicted result can well capture the features of moisture profile. During the thawing period, the moisture in soil drains downward and the soil moisture becomes lower comparing to the freezing condition. The position of the peak moisture content during the thawing condition is relatively lower than that during the freezing condition. The proposed result can agree reasonably with the measured data. It is emphasized that a subgrade of the layered structure has much lower water contents in the upper part. The measured distributions may be different from the hypothetical condition. This case is hard to be modelled in a numerical code [15], while the proposed model can reasonably reproduce the tendency.Figure 7
The predicted and measured volumetric water content profiles at K977 and K1004 sites.
(a)(b)


## 5. Conclusions

The mechanism of vapour migration and phase transformation in unsaturated freezing soils has been less understood. In order to avoid the great uncertainty caused by the phase change term of vapour-water-ice in the numerical iteration process, this study presents numerical solution schemes for coupled vapour and heat transfer at above-zero and subzero temperature conditions, and hence freezing-thawing cycles. The new model is solved numerically by the finite element method, which is, then, validated by comparing to the laboratory measurement and in situ observation. The main findings are as follows:(1)
The new model avoids using the local equilibrium assumption that accounts for water-vapour-ice phase change, and an adsorption-desorption equilibrium relation is proposed to direct describe phase change between vapour and ice phase, which provides a novel way for modelling the water-heat coupling movement problem.(2)
The numerical model has a close match with the results of laboratory freezing experiment, which shows a good performance in numerical instability and mass or heat conservation.(3)
The numerical solutions performed quite well for the field application that considers rapidly varying surface temperatures. The proposed model can generate a reasonable result for the frost or thawed depth and soil moisture profiles.(4)
It should be noted that the new model should be tested eventually against a more realistic and complete field data set, which can provide detailed information on soil properties. The new expressions used for vapour transfer and the adsorption-desorption equilibrium relations are applicable to both frozen and unfrozen conditions. The new expressions can replenish the existing models, which should be tested for more kinds of soils.



---

*Source: 8874919-2020-10-31.xml*