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

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. *e 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. *e model is established by using COMSOL Multiphysics, which is a multiphysics simulation software through finite element analysis. *e 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. *e 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.


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][2][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][5][6][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. e specific phenomenon was named as the canopy effect by Zhang et al. [9]. e 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. e 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. is 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][14][15][16][17]. Some of these studies take the contribution of vapour flow and its phase change into account, while others choose to neglect yet. e 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. e 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]. e 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]. e 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. e 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][27][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. us, the local equilibrium theory would be challenged in case of being used to describe the moisture transfer in frozen soil. irdly, the mechanism of multiphase transfer in unsaturated freezing soils has been less studied. ey 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]. e 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. erefore, 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. e 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. Figure 1 depicts the physical system which is a vertical, onedimensional soil profile covered by an impermeable plate. e plate represents the infrastructures that covered on the soil surface. e 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.

Mathematical Model
e condensed water will continue freezing into pore ice in the soil where the temperature drops below 0°C. e transient profiles of soil temperature (solid line in Figure 1) and total water saturation (dashed line) are the outcomes of this problem.
Some assumptions are made to simplify the quantitative description of the coupled heat and moisture process, as follows: (a) e 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) e freezing front is the depth where the soil temperature is 0°C. (c) e ice lens zone is the region where the soil is saturated. e 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].

Governing Equations for the Vapour Transfer
Zone. e one-dimensional governing equation for vapour flow in the unsaturated freezing soil is given by the following mass conservation equation: 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 ). e subscripts w, v, and i denote the liquid water, vapour, and pore ice, respectively. q v is vapour flux (kg m −2 s −1 ), and S is the total saturation of liquid water and pore ice (unitless); here, S � S w + S i ρ i /ρ w .
2 Advances in Civil Engineering e vapour movement in unsaturated freezing soil is recognised to be driven by both the pressure gradient and vapour concertation gradient [33]. e vapour mass flux in this region according to Darcy's law and Fick's law is where K is the air permeability of soil (m 2 ), μ 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 (m 2 s −1 ), and τ is the tortuosity factor (unitless). e parameter τ is the ratio of real length of transfer path to apparent length, which is assigned to be 1.2 in this study [34]. e vapour phase in unsaturated freezing soil can be considered as ideal gas. It follows the ideal gas law, P � ρ v RT.

Substituting this equation into equation (2) leads to
where K S and K T 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 ). e mass transfer equation can be obtained by substituting equation (3a) into equation (1): It has been revealed that the heat conduction is two orders of magnitude larger than the sensible heat of vapour flow [35].
us, the heat conservative equation neglects the term of the sensible heat of vapour flow in the soil and is written as 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  Advances in Civil Engineering 3 subscripts s in the abovementioned equations indicate the solid phase of soil particles. e 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: where T 0 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]. S uw 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 S uw , the criteria for determining the saturation of pore ice and liquid water can be obtained [7], as follows: In the hygroscopic porous material, the three phases of water, i.e., unfrozen liquid water, ice, and vapour coexist in the freezing material pores. e 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. e equation is expressed as where P vs (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. c is an empirical parameter, a suggested value is 5000 according to Wang et al. [45]. e simulated model consists of mass conservation equation (4) and energy balance equation (5), which are highly nonlinear and coupled.
ere are four unknown variables in the two equations, S w , S i , T, and ρ v . Equations (6)-(8) provide the other two additional relations by adding the saturation of unfrozen water S uw , such that the simultaneous equations can be solved mathematically.

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. e mass balance at this region can be expressed by where z fs (t) represents the position of the frost surface of the ice lens, which is a function of time. e governing equation for heat transfer in the ice lens zone can be written as where C i 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 ).

Heat and Mass
Balance at the Frost Surface. e heat and moisture mass balance at the moving front can be expressed by e water vapour on the frost surface is assumed to be saturated.
Two new parameters are added into the model, z fs 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].

Soil Hydraulic Properties.
In complete dry soil, the air permeability K is equal to inherent permeability k s . As for the unsaturated freezing soil, an impedance factor is added to the air permeability function [9,12,19].
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, S i becomes 0, and equation (13) reverts to the formula as defined by Brooks and Corey [19]. e 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]. e harmonic averaging method is adopted here to compute the vapour diffusivity D based on the molecular diffusion coefficient D m and the Knudsen diffusion coefficient D k [55]: where D m can be written as follows [56]: where D a is the diffusivity of water vapour in air (m 2 s −1 ), as given in Table 1. e expression of D k was given by Geankoplis [57] as follows: where d is the averaged soil pore diameter (m) and m w is the molecular weight of water (kg mol −1 ). e values of related parameters are presented in Table 1. (4), (5a), (5b), and (5c) are highly nonlinear since primary variables vary with time and primary variables. e 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].

Numerical Implementation. Equations
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.

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]. e temperatures of the top and bottom ends of the specimen are controlled accurately by a pair of plates. e 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. e physical and hydraulic properties of the sample are presented in Table 2.
e measured Soil Freezing Characteristic Curve is shown in Figure 2, where the curve is fitted by equation (6). e soil specimen is placed into the cylinder with a controlled dry density of 1.40 g/cm 3 . e final height of the sample is 13.5 cm.
ree 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.
Five experimental cases are designed as listed in Table 3. In all the 5 cases, the top end of the sample is subjected to a   [49] Advances in Civil Engineering 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.

Simulation Result.
e computed and measured water content profiles for cases 1 and 2 are shown in Figure 3. 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. e 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.
e simulated result has a quite good match with the tested data, where only one peak value at the top surface can be observed. Figure 4 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.

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. e subgrade at one site (K977) is built on the undisturbed ground surface, while the other is in a cut section (K1004). e measured data start from Nov 1, 2013, to Oct 17, 2014, that last for 351 days. e upper layer of the subgrade structure is a wellgraded gravel with cement (thickness is 70 cm). e materials in the second layer are group A/B fills (thickness is 230 cm). e 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. e 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. e measured near-surface temperature is shown in Figure 5, where a fitting function is obtained as the input upper boundary. e 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. e thermal properties of soil can be found in Table 1. e other parameters are the same as those for silica sand, which can be found in Zhang et al. [9]. e hydraulic properties of the soil are the same as those in the work of Teng et al. [32]. Figure 6 presents the predicted frost and thawed depth variation in the 160 days, which indicates one freezing and thawing cycle. e frost and thawed depth are usually defined as the depth where the soil temperature is zero.

Freezing-awing Simulations.
us, 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. e 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 7 presents measured and predicted water content profiles in the subgrades of the two sites during the       Advances in Civil Engineering 7 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. e 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. e position of the peak moisture content during the thawing condition is relatively lower than that during the freezing condition. e 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. e measured distributions may be different from the hypothetical condition. is case is hard to be modelled in a numerical code [15], while the proposed model can reasonably reproduce the tendency.

Conclusions
e 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 freezingthawing cycles. e new model is solved numerically by the finite element method, which is, then, validated by comparing to the laboratory measurement and in situ observation. e main findings are as follows: (1) e 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) e 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) e numerical solutions performed quite well for the field application that considers rapidly varying surface temperatures. e 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. e new expressions used for vapour transfer and the adsorption-desorption equilibrium relations are applicable to both frozen and unfrozen conditions. e new expressions can replenish the existing models, which should be tested for more kinds of soils.

Notations b:
Empirical parameter of the BC model, unitless c: Empirical parameter of air permeability in frozen soil, unitless C: Effective volumetric heat capacity, J m −3 K −1 C i : Volumetric heat capacity of ice, J m −3 K −1 C s : Volumetric heat capacity of soil particle, J m −3 K −1 C v : Volumetric heat capacity of vapour, J m −3 K −1 C w : Volumetric heat capacity of water, J m −3 K −1 d: Average

Data Availability
A theoretical model is established to formulate the coupled thermal and hydrological process by COMSOL, where the vapour flow governs the mass transfer process. A series of laboratory test results, which are published on Geotechnique [17], are used to validate the numerical formulations/codes.

Conflicts of Interest
e authors declare that they have no conflicts of interest.