A frost heave model of unsaturated coarse-grained soil considering vapour transfer

Substantial frost heave has been observed in coarse fills in high-speed railway embankments. These coarse fills have low fine contents and very low water content. The groundwater table is located below the coarse fills. The coarse fills were considered not susceptible to frost heave. Recent experimental results in the literature showed that vapour transfer has a considerable influence on the frost heaving of unsaturated coarse-grained soil. But vapour transfer has been rarely considered in the modelling of frost heave. This study presents a new frost heave model with considering vapour transfer and its contribution to ice formation. The rigid ice theory is applied to initiate an ice lens formation in the frozen fringe. An updated computer programme PCHeave is developed by considering the vapour transfer. The results of the proposed model are compared with laboratory test results, which show reasonable agreement. The prediction of the model agrees well with the measured frost heave and frost depth, which indicates that the proposed model can reasonably reflects the process of frost heave in unsaturated coarse soil.


INTRODUCTION
In seasonally frozen soil areas, the effects of frost heave on high-speed railways are significant because of stringent deformation requirement of truck foundation. The Harbin-Dalian Passenger Dedicated Railway is the first dedicated passenger line in the world to operate in cold regions. These coarse fills for high-speed railways, classified as Group A and Group B in the Chinese standard, have low fines contents (typically, <10%) and very low initial water content (typically, 5% in mass). The groundwater table is located far below the coarse fills. The coarse fills were considered not susceptible to frost heave, but substantial frost heave has been observed in the coarse fills in high-speed railway embankments [1][2][3][4]. Therefore, it is important to explain the unusual phenomenon and predict the frost heave of coarse-grained soils.
A number of theories in the literature have been developed to predict frost heave. Five main streams of the theories can be identified: (1) the capillary theory based on the Young-Laplace equation for primary frost heave [5]; (2) the rigid ice model referred to as the secondary frost heave that assumes the existence of a frozen fringe [6][7][8]; (3) the hydrodynamic models or the thermo-hydro-mechanical models [9][10][11][12]; (4) the segregation potential models, which are based on the proportional relation between the velocity of water entering the unfrozen soil and the temperature gradient [13][14]; (5) the theory of premelting dynamics that considers the intermolecular forces causing premelted fluid to migrate and supply segregated ice growth and frost heave [15][16][17].
In all the models and theories discussed above, including the more recent premelted ice model, frost heave is considered to be a consequence of freezing and migration of liquid water in fine grained soil. Indeed, the frozen fringe is assumed to be fully saturated by water and ice, both in the rigid ice model, in the segregational potential theory and in the premelted ice model. According to the above theories, liquid water in coarsegrained soils above the groundwater table is usually discontinuous. It is then difficult for liquid water to flow to the freezing front or the frozen fringe. The frozen fringe in coarse-grained soils such as sand or gravel is very thin or zero [6], which prevents the formation of ice lens. According to existing models and theories, freezing of a coarse-grained soil does not cause the formation of ice lenses or much frost heave.
The experimental results in the literature show that vapour transfer has an important effect on the frost heave of coarse-grained soil [4]. Teng et al. performed a series of laboratory tests to investigate the factors that will affect vapour transfer in a clean sand, including initial water content, dry density, boundary temperature, and water supply pattern [18]. They found that vapour transfer can be the primary mechanism of moisture migration in coarse-grained soil, especially when the initial water contents of the soils are low. Under subfreezing temperatures, ice crystals in coarse-grained soils will further promote the migration of vapour, and the ice lens will grow continuously, leading to frost heave in the soils [12,19]. The ice crystal is equivalent to a water suction pump, which increases the matric suction of the soil and reduces the concentration of vapour, thus continuously driving the migration of vapour to the freezing zone [18]. While the physical process of frost heave in coarse-grained soils is relatively clear and validated by experiments, the corresponding mathematical model is still lacking. It is relatively unknown: how to model frost heave in coarse-grained soils and how much frost heave can form due to vapour transfer.
The objective of this study is to develop a new frost heave model that takes into account the vapour transfer and phase change in coarse-grained soils. The numerical model is described in detail and further implemented into the program PCHeave. The results of the numerical modelling are then compared with laboratory test results. The results show that the model simulation matches well with the test results. The significance of vapour transfer in coarse-grained soils is analysed based on the proposed model, and finally, some concluding remarks are summarized.

Basic Assumptions
PCHeave is a computer program that developed by Sheng et al. [7] that can simulate the formation of discrete ice lenses in soil and predict the onedimensional frost heave. In PCHeave, a frozen fringe is assumed to exist that is always saturated by water or ice. This is in accordance with the rigid ice concept proposed in Miller [6]. Although PCHeave is one of the few models that can be used to predict frost heave in engineering practice, it cannot explain the frost heave in coarse-grained soil caused by vapour transfer. This study further develops the PCHeave model into a new version to account for the effect of vapour transfer on the accumulation of pore ice. A basic assumption that the vapour inflow at the freezing front will condense to water and fill the pores of the frozen fringe is added in the updated model.

Ice Lens Initiation
A new ice lens is assumed to be initiated when and where the effective stress in the soil approaches zero [7]. Through the Clapeyron equation in thermodynamics, the water pressure at the warm side of the lens, which appears in suction, is affected by the segregational temperature and the overburden pressure. Assuming the ice-water interface is at equilibrium, the ice pressure, water pressure and temperature are related through the integrated Clapeyron equation [1][2]7]. It is shown that the effective stress approaches zero when and where the new ice lens is initiated.

Heat and Mass Transfer
In order to establish the heat and mass balance equations, we first study the mass transfer within the frozen fringe. At a time level t n , we assume that the frozen fringe is located between x f and x b , where x f is the frost front and x b is the warm boundary of the latest ice lens, Fig. 1. Within the fringe, we denote the soil solid fraction by the area filled with gray block, liquid water by the area filled with yellow block and air by the area filled with sparse grid lines.
The volumetric ice content is 1 above x b , I b (not necessarily equal 1) at x b and zero at x f . If the rigid ice body moves upwards at a velocity V i , a space of the size ሺͳܸ ݀‫ݐ‬ሻ will be left during time period dt and must be filled by new formed ice. This volume of ice is represented by the area filled with white block in Fig. 1, which includes the ice formed at x b , i.e. ‫ܫ‬ ܸ ‫,ݐ݀‬ and that within the fringe, i.e. ሺͳ െ ‫ܫ‬ ሻܸ ‫ݐ݀‬ . In addition, the penetration of the frost front ሺെ݀‫ݔ‬ ሻ will also cause an ice formation of the volume of ‫ܫ‬ሺെ݀‫ݔ‬ ሻ, where ‫ܫ‬ is the average ice content in the frozen fringe. This amount of ice is represented by the area filled with close grid lines in Fig. 1. Assume the location of the frozen fringe is known at the time t. The aim is to determine the rate of heave during a small time step ǻt and the new location of the frozen fringe at time t+ǻt.
The heat balance is considered for each soil layer. At the warm surface of the latest ice lens (x b in Figure 1), the heat flow rate to x b plus the release rate of latent heat by phase change must equal the rate of heat removal upward through the frozen zone. For a steady thermal state on each side of x b , the heat balance can be expressed as˖ ( ) where ߣ is the effective thermal conductivity of the layer above x b , ߣ is the effective thermal conductivity of the current frozen fringe, x c is either the cold boundary of the soil profile for coarse calculation or the location of the cold boundary of the growing ice lens, T s is the temperature at the ice-water interface or the socalled segregational temperature, T c is the temperature at location x c , T f is the temperature at location x f , ȡ i is the density of ice, L is the specific latent heat of water.
Within the frozen fringe, the rate of ice formation is V i due to rigid ice motion, in addition to ‫ܫ‬ሺെ݀‫ݔ‬ Ȁ݀‫ݐ‬ሻ due to frost front penetration. The vapour inflow at the freezing front will condense to water and fill the pores of the frozen fringe. Therefore, the overall heat balance for the frozen fringe can be expressed as follows: where dx f /dt is the advancing rate of the frost front x f or the frost penetration rate during ǻt and is negative when frost advances downwards, ߣ ୳ is the effective thermal conductivity of the unfrozen zone between x f and x w (Figure 1), x w is the boundary of the soil layer below the frost front, T w is the temperature at warm side x w , s r is the degree of saturation, ‫ܫ‬ is the mean ice content in the current frozen fringe, n is the soil porosity, ȡ w is the density of liquid water, L vap is the specific latent heat of vapour. The first term on the right-hand side represents the latent heat for the new ice mass formed in the frozen fringe. The second term represents the latent heat for the new liquid water formed in the frozen fringe. The unfrozen soil is unsaturated, but the frozen fringe is saturated. So, the vapour inflow at the frost front will condense to water and fill the pores of frozen fringe. The mass conservation at x b requires that the water flow to x b equals the ice mass formed there: where v ff denotes the rate of water flow in the frozen fringe.
The overall mass balance within the frozen fringe requires that the outflow of ice mass at x b must balance the inflow of water mass at x f plus the mass loss in the frozen fringe as follows: where v u is the rate of water flow at the frost front (in the unfrozen soil). The second term on the right-hand side is the mass loss due to phase change, and the third term is the water mass needed to fill the initially unsaturated soil.
In the case where the frost front is above the groundwater table, the rate of water flow in the unfrozen soil (v u ) can be determined from Darcy's law as follows: where k u is the effective permeability of the unfrozen soil between the frost front and the groundwater table, g is the acceleration of gravity, x gw is the location of the groundwater table, u w (x gw ) is the pore water pressure at the ground water table, u w (x f ) is the pore water pressure at the frost front x f , and v vap is the vapour flux in the unfrozen zone. The pore water pressure at the groundwater table is usually zero, unless an excess pore pressure is considered.
The vapour flux is computed as follows: where K vh and k vT are the isothermal and thermal vapour hydraulic conductivities, respectively. They can be estimated by the method of Zhang et al. [20].
If the water flow rate in the frozen fringe (v ff ) is assumed to be constant, the pore water pressure in the frozen fringe can be found by integrating from Darcy's law, with the assumed permeability function as follows: Substituting Eq. (14) into Eq. (13) to eliminate k ff , and integrating for u w leads to the following: The integration constant C can be determined from the boundary values at x b, where the pore water pressure is related to the pore ice pressure (overburden) and the temperature via the generalised Clapeyron equation as follows: where ı is the overburden pressure at location x b , k u is the saturated permeability of the soil with in the frozen fringe (before freezing), and b is a parameter that will be determined below. The pore water pressure at the position of x b and x f can be given as follows: Equations (1)-(6) and (10) form the basis of the frost heave model. Providing the location of the frozen fringe (x f , x b ), the ice content (I) within the frozen fringe and the material parameters involved are known, and they can be solved by iteration for the six unknowns T s , V i , dx f /dt, v ff , v u and u w (x).
If the frost front penetrates below the groundwater table, the pore water pressure at the frost front is then determined by the following: In this case, Eqs. (1), (2), (3) and (10) are still valid. Together with Eq. (18), the five equations are solved for the five unknowns T s , V i , dx f /dt, v ff and u w (x).

Steady State of Growth
In the case where the frozen fringe is too thin for computational significance, the base of the latest ice lens is set to the frost front. The rate of water flow in the unfrozen zone is then in equilibrium with the rate of ice lensing; thus, the frost front will remain stationary. The following equation is used instead of Eqs. (1) and (2) : The following mass balance equation is used instead of Eqs. (3) and (4) : Equations (5) and (10) are still valid. Together with Eqs. (13) and (14), the four equations can be solved for the four unknowns: T s , V i , v u , and u w (x).
So far, the governing equations have been established based on knowledge of the location of the frozen fringe. The next step is to locate the frozen fringe at each time step. If the frozen fringe (x f , x b ) is known at time t, the frost penetration rate dx f /dt can be determined and the new position of the frost front at time t+ǻt is x f +(dx f /dt)ǻt. The position of the warm surface of the latest ice lens x b remains fixed unless a new ice lens appears within the frozen fringe. A new ice lens appears when and where the effective stress vanishes.

Computation Strategy
The governing equations contain a number of material parameters that must be determined. Such parameters include the effective thermal conductivity and the permeability of a multi-layer structure, the thermal conductivity and the permeability of a multiphase material, and the ice content of the frozen fringe. The calculation of these parameters use the method given by Sheng et al. [7]. The final parameters used in PCHeave include the bottom position of the soil layer, the thermal conductivity of the soil solids, the water content by dry weight, the dry density, the degree of saturation, the saturated permeability and the unfrozen water content at -1 o C (in percentage of the initial water content). This parameter is kept as an input so that nonsoil materials such as insulation, snow and concrete can also be simulated in PCHeave.
In addition to the material parameters, the boundary conditions, groundwater table, external load and geometric make-up of the soil profile must also be specified. The general input parameters in PCHeave include the number of material layers, the total depth of the soil profile, the number of temperature intervals at the cold and warm boundaries, the overburden pressure at the ground surface, the distance to the groundwater table from the ground surface, the total time length and the time step. The temperatures at the cold and warm boundaries can be specified step-wise in a number of intervals. The overburden pressure is the extra load applied at the cold boundary. The soil self-weight is calculated automatically in PCHeave. The size of the time step can affect the results, but the solution converges as the step size is sufficiently small. As the computation in PCHeave is extremely fast (usually in seconds), an appropriate step size can be found by trial and error.
The calculation procedure in PCHeave is illustrated in Fig. 2. The position of the first ice lens is determined according to the theory discussed in the section location of the frozen fringe. The heat and mass balance Eqs. (1) -(6) and (10) are then solved by iteration, for T s , V i , dx f /dt, and u w (x). The length of the soil column is elongated by V i ǻt. The frost front moves (-dx f /dt)ǻt downwards. The position of the new frozen fringe is determined according to the ice lens formation criterion. The system is then ready for the next time step.

Model validation
A series of one-side freezing experiments of a coarse-grained soil were conducted by Gao et al. [19] to test the effect of vapour transfer on the frost heave of a coarse-grained soil. Graded crushed stones containing 5% fines by weight are adopted as the test material. The non-uniformity coefficient of the cumulative curve of the test material is 15.78, and the coefficient of curvature is 2.63, which indicates that the test material belongs to a well graded soil. The soil specimen is 198 mm in diameter and 200 mm in height. The dry density is controlled at 2.087 g/cm 3 , and the initial water content by weight is 8%. During the whole tests progress, nopressure water supplement were performed. The top temperatures and the bottom temperatures are kept constant at -12ႏ and 2ႏ, respectively. The sidewall of the sample is wrapped with thermal insulation material to achieve a unidirectional freezing state from top to bottom. The test period is 12 days. The experimental results are used to validate the proposed model in PCHeave. The material parameters, the general input parameters and the boundary temperature specification in PCHeave are given in Table 1. The predicted results without considering vapour transfer are generated by the original model in PCHeave, while the results of considering vapour transfer is computed by the proposed model (Fig. 3). The predicted frost depth and frost heave after 288 h are 183 mm and 5.2 mm, respectively, while the measured values for them are 181 mm and 5.5 mm, respectively. The predicted results of the proposed model agree well with the measured data. When using the original model of PCHeave, the predicted frost heave is 1.9 mm after the test, which is considerably less than the measured result of 5.5 mm. This finding is consistent with the results of Zhang et al. [3] that the contribution of vapour flow to total water flux can be significant in the unsaturated freezing coarse soils [20]. It also highlights the importance of vapour transfer in modelling frost heave in coarse-grained soils. The predicted frost depth by the original PCHeave model shows a faster descending rate compared to the experimental values, which may be caused by neglecting the vapour transfer. Phase change from vapour into ice is an exothermic process, and can delay the descent of the freezing front. The updated PCHeave code can simulate ice lens formation in soils. Figure 4(a) presents the forms of ice lenses at different times when considering vapour transfer, and Figure 4(b) presents the simulated result without considering vapour transfer. Figure 4(a) shows that a few ice lenses are distributed in the frozen zone of the soil specimen, and the thickest ice lenses occur near the final freezing front. However, if the vapour transfer is E3S Web of Conferences 195, 02017 (2020) E-UNSAT 2020 https://doi.org/10.1051/e3sconf/202019502017 neglected, the ice lens only occurs at th freezing front, and there are no ice le positions. The reason is that the liq bottom of the sample can hardly mov vapour will flow upwards driven by the humidity gradients. As the permeability to the ice lens decreases and the free downwards, the ice lens will generate at (a) considering vapour trans (b) without considering vapour to describe frost ain feature of this r flow, which is es. The proposed rogram PCHeave. ment is used to posed model. The frost depth show lues of laboratory reasonably reflect ined soils.
The updated PCHeave prog many input parameters for properties, it is very simple to extremely fast. The proposed investigate frost heave in coars