A finite volume-based model for the hydrothermal behavior of soil under freeze–thaw cycles

Freeze–thaw cycles in soil are driven by water migration, phase transitions, and heat transfer, which themselves are closely coupled variables in the natural environment. To simulate this complex periglacial process at different time and length scales, a multi-physics model was established by solving sets of equations describing fluid flow and heat transfer, and a dynamic equilibrium equation for phase changes in moisture. This model considers the effects of water–ice phase changes on the hydraulic and thermal properties of soil and the effect of latent heat during phase transition. These equations were then discretized by using the finite volume method and solved using iteration. The open-source software OpenFOAM was used to generate computational code for simulation of coupled heat and fluid transport during freezing and thawing of soil. A set of laboratory freezing tests considering two thermal boundary conditions were carried out, of which the results were obtained to verify the proposed model. In general, the numerical solutions agree well with the measured data. A railway embankment problem, incorporating soil hydrothermal behavior in response to seasonal variations in surface temperature, was finally solved with the finite volume-based approach, indicating the algorithm’s robustness and flexibility.


Introduction
Unsaturated soil at or near to the ground surface may experience frequent freeze-thaw cycles driven by seasonal climate change. The associated change of soil temperature, fluid flow, and formation of ice all affect the long-term performance of infrastructure, posing problems for engineers [1][2][3]. As such, being able to predict fluid flow patterns and temperature gradients in soil is of great significance for geotechnical projects in cold regions.
The successive freezing and thawing of soil mass are essentially a heat and mass transfer process, which is governed by thermal conduction, the migration, and phase transformation of water in the subsurface environment. The migration and freezing of water are closely coupled processes, and jointly influence the temperature field and the rate of heat transfer within soil. Two types of mathematical models have previously been developed to describe this coupling strategy: mechanistic models based on the viscous flow of liquid water in porous media and the thermal balance of soil [4], and thermodynamic models that apply principles of irreversible thermodynamics for water and heat fluxes within soil [5,6]. A variety of numerical codes use mechanistic models when applied to geotechnical engineering problems, which is a technique that we adopt in this study. Many previous studies have laid a solid foundation for the parameterization of complex numerical models. Harlan [7] presented the first finite difference solution to a one-dimensional hydrothermal problem regarding the freezing and thawing of a homogeneous, porous medium. Subsequently, Newman et al. [8] presented an approach for predicting heat and mass transfer in freezing, unsaturated soil. In their work, characteristics for soil-ice and soil-water mixtures were used to combine heat and mass transfer relationships into a single equation for freezing (or frozen) soils. Hansson et al. [9,10] presented a robust, energy-and mass-conservative solution to a coupled heat transport and variably saturated water flow, which was designed to apply to a temperature range across 0˚C. In addition, several unusual water-migration phenomena that occur during soil freezing have also been reported and studied comprehensively, such as the frozen fringe and its suction effect [11,12], the canopy effect [13,14], and the excess pore-water pressure effect [15].
Based on these investigations, several codes have been developed for modeling the freezethaw cycles, including the HYDRUS model [16], the FROST model [17], the SHAW model [18], the SOIL model [19], and the VAPOR model [20]. These codes are sometimes unsuitable for addressing more complex problems, as they are limited to 1D or 2D problems. In addition, because of their geometric adaptability, they can only be applied to structured grids, and these codes do not support parallel computing, which limits the size of problems to be investigated. In general, simulation capabilities accounting for complex geometries and high-resolution problems are essentially limited.
In light of this, we have adopted an unstructured mesh-based code to addresses these gaps in modeling capabilities, namely the finite element method (FVM). In the FVM, volume integrals in a partial differential equation that contain a divergence term are converted to surface integrals, using the divergence theorem. These terms are then evaluated as fluxes at the surfaces of each finite volume. In contrast to the finite element method and finite difference method, the FVM evaluates exact expressions for the average value of the solution over volume, and uses this data to construct approximations of the solution within cells. The FVM-based approach shows notable advantages of automatic matrix construction, solution capabilities for partial differential equations, and ease of implementation in some open-source platforms, such as OpenFOAM [21]. The FVM code developed in this study was validated by laboratory freezing test data, and further tested with seasonal simulation of a railway embankment case involving heat transport and water flow in response to seasonal variations in surface temperature.
The remainder of paper is organized as follows: Section 2 presents the fundamentals for the coupled heat transfer and water migration process in the soil domain; Section 3 describes the finite volume method and developing strategy for the coupled model; Section 4 presents the validation and application of the FVM code; Finally, conclusions and some proposals for future investigations are given in Section 5.

Heat transfer
Transitions between unfrozen, frozen, and thawing soil are necessarily accompanied by conductive heat transfer and changes in the physical state of H 2 O. According to heat transfer theory, the control equation for unsteady heat transfer can be formulated in terms of the latent heat of phase changes being either a heat source or sink [22], given as where C is the volumetric heat capacity of the soil, J�m −3 �K −1 ; T is the soil temperature, K; λ is the soil thermal conductivity, W�m −1 �K −1 ; L is the latent heat released as water transforms into ice or vice versa, J�kg −1 ; ρ i is the ice density, kg�m −3 ; and θ i is the volumetric ice content. During the freezing or thawing of soil, pore water is redistributed and θ i changes simultaneously. The physical and mechanical properties of the soil also change due to variations in soil-water-ice proportions. The volumetric heat capacity of soil can be expressed as the sum of the volumetric heat capacity of the solid matrix, liquid water, and ice phases multiplied by their volumetric fractions, expressed as [23] where θ s is the volumetric soil matrix content, C s is the volumetric heat capacity of soil matrix, J�m −3 �K −1 ; C i is the volumetric heat capacity of ice, J�m −3 �K −1 ; θ u is the volumetric liquid water content, and C u is the volumetric heat capacity of liquid water, J�m −3 �K −1 .
The effective thermal conductivity of soil depends on the thermal conductivity of all its components (i.e., soil matrix, liquid water, and ice), their volumetric fractions, and the spatial distribution of each phase. Thermal conductivity is generally obtained by a logarithmic law [23] from the physical properties of the soil matrix, ice, and water contents as follows where λ s is the thermal conductivity of the soil matrix, W�m −1 �K −1 ; λ i is the thermal conductivity of ice, W�m −1 �K −1 ; and λ u is the thermal conductivity of water, W�m −1 �K −1 . Soil moisture consists of liquid water located in pores or fractures within the soil, and ice located around the grain surface. Because of the difference in densities of ice and water, the volumetric moisture content of frozen soil is calculated by where ρ w is the density of liquid water, kg�m −3 .

Moisture migration
The water migration model is derived from the application of Darcy's law to moisture movement in a permeable media. The migration of liquid water in soil at temperatures across 0˚C is described by the modified Richards equation [23] @y u @t where D(θ u ) is the soil water diffusivity, m 2 �s −1 ; K(θ u ) is the hydraulic conductivity of unsaturated frozen soil, m�s −1 . The right first term in Eq (5) accounts for the fluid flow due to water concentration gradients, while the right second term represents the fluid flow by gravity. The water diffusivity in soil is then calculated by [9,10]: where C(θ u ) is the specific water capacity, m −1 , I(θ i ) represents the blocking effect of ice on water migration in frozen soil [9,10], calculated by In this study, the hydraulic properties are described by [24,25] where S is the relative saturation of the frozen soil; a 0 , l, and m are empirical parameters; K s is the saturated hydraulic conductivity, m�s −1 . The normalized saturation is given as [25] S ¼ where the θ r is the residual volumetric water content, and θ s is the volumetric water content of saturated soil.

Coupling strategy
In the heat transfer and water migration models, there are three unknown variables, namely temperature, ice content, and liquid water content. An additional equation is required for a closed-form solution. In freezing soil, the unfrozen water content is seen as a function of temperature, generally determined from soil freezing characteristic curve [8]. If liquid water and the sub-zero temperature of frozen soil are in a dynamic equilibrium state, the relationship between liquid water contents and below-freezing temperatures is given by [26] w 0 where w 0 is the initial water content of unfrozen soil, w u is the liquid water content of frozen soil, T f is the freezing point, K, and B is a parameter related to the salt content and soil type, which can be determined by either using the one-point method [26] or by adopting typical values for silty sand (0.58) and sandy gravels (0.65). This relationship between liquid water contents and sub-zero temperatures can be represented as Eqs (11), (12) and (4) are combined to describe a relationship between the ice content, liquid water content, and temperature, expressed as Combining Eqs (1), (5) and (13) leads to the framework of a hydrothermal model. With proper boundary conditions, the temperature (T), volumetric liquid water content (θ u ), and volumetric ice content (θ i ) of a soil can be calculated. However, since analytical solutions to these equations are only available in special cases, numerical techniques must be employed to solve them.

Numerical solution
Equations for heat transfer and water flow are both nonlinear. In this study, the Picard iteration method [27] is adopted to linearize both equations by allowing the nonlinearities to lag by one step behind. More specifically, both heat transfer and water flow equations are mathematically coupled through their phase change component, meaning that a mutual solution must be generated by iteration. In this section, we introduce a calculation method to address the coupling issue and the finite volume-based procedure for the subsequent simulation.

Solution strategy
Eq (14) is used to eliminate the ice content variable in the definitions of heat and mass transfer, leaving both with only two undetermined variables.
where the parameter r is expressed as Then, the time derivative of θ i can be rewritten as Substituting Eq (16) into Eq (5) produces In numerical calculations, the iterative process for each time step is given as: 1. The temperature T i k+1 at the end of the predicted time step T i k is approximated by the temperature at the beginning of this step. The i k+1 component of the volumetric ice content (θ i ) and liquid water content (θ u ) i k+1 term at the end of this period is obtained by the dynamic equilibrium defined in Eq (13).
2. The (θ i ) i k+1 expression obtained by the water diffusion coefficient D i k and the result produced in Step (1) allows Eq (5) to be solved, and so the (θ u ) i k+1 of each node at the end of this step is obtained. Furthermore, Eq (14) is used with this expression for (θ i ) i k+1 to solve , and the results of step (2) allows the temperature control Eq (1) to be solved, and the real T i k+1 of each node at the end of the period can be determined. This is both the correction value of this iteration step and the forecast value of the next iteration step.
The calculation is continued until the difference between unfrozen water content and temperature from two consecutive cycles satisfies the predefined tolerance.
The general process of this approach is shown graphically in Fig 1.

Finite volume method
All equations in this model have been spatially discretized using the finite volume method, which is widely used for numerical simulation of problems involving fluid flow and heat and mass transfer. As the finite volume method can be implemented on unstructured polygonal meshes, it is well-suited for modeling the hydrothermal process in frozen soil. In this study, the open-source C++ library OpenFOAM was selected for code execution, and the following parts describe the implementation of this solver. All governing equations presented in Section 2 can be written in the general form: The parameters of interest are given in Table 1. The first step of the discretization process is to partition the whole domain into non-overlapping elements, also known as finite volumes. The partial differential equations are then discretized into algebraic equations by integration over each discrete element. The parameters used in this discretization process are shown in Fig 2. For example, consider two adjacent control volumes P and N, which share the plane f. Here, the vector S f is the normal vector to f and the module is the area of f. By integrating Eq (18) over the control volume, Eq (18) Replacing the volume integrals of the Laplace terms by surface integrals based on the Gauss theorem, Eq (19) Replacing the surface integral by a summation over the control volume faces, the discretization equation on each cell can be described by

PLOS ONE
Finite volume-based model for hydrothermal behavior of soil Then, the Laplacian terms in Eq (21) are discretized by the central difference scheme as described in Eq (22), which is always of second-order accuracy: Calculation of the diffusion coefficient at the cell surface requires a conservative interpolation based on a harmonic mean, allowing the diffusion coefficients to guarantee continuity of both flux and current. For the two control volumes P and N shown in Fig 2, each of which has homogeneous diffusion coefficients (Γ P and Γ N , respectively), the harmonic mean of the diffusion coefficient can be calculated by Eq (21) must be solved by a time step procedure. In addition, due to convergence and stability requirements in the calculation, a fully implicit backward difference scheme is used to

PLOS ONE
Finite volume-based model for hydrothermal behavior of soil discretize the time derivative term in Eq (21). The time variables are discretized, shown in Fig  3. This allows the discrete-time set and the corresponding series of time intervals [t n − 1 , t n ], n = 1, . . ., N to be derived, and the time step is: The time derivative at time t n can then be approximated by Eq (25), and all equations are solved step-by-step for each time interval. @r� @t Finally, the discretization of the above equations produces a set of algebraic equations that are solved to obtain discrete values.

Validation and applications
Validation Experiment description. The theory and numerical solution presented in Sections 2 and 3 were validated against the results of one-dimensional soil freezing tests. The soil sample used in this experiment was collected from an embankment along the Baotou-Lanzhou Railway (BLR), North China, which is recognized as the first railway crossing a desert area of China. The BLR was constructed in the 1950s and has since suffered from severe frost damage due to the underlying natural surface soil, which was not well treated with respect to the frost issues. The soil utilized in our experiment was excavated at the BLR design milepost DK72+612.
Particle analysis of this soil showed the following composition: 0.74% fine gravel, 39.87% sand, 45.74% silt, and 13.65% clay. This soil also has a maximum dry density of 1.89 g/cm 3 , optimal moisture content of 13.0%, specific density of 2.74, liquid limit of 26.6%, and plastic limit of 13.4%. The sample was molded into a cylinder with 10 cm in height and 15 cm in diameter. Standard procedures for sample preparation were adhered to, such that the relative compaction was 90%, and the moisture content was 17.5%.
Freeze-thaw testing was conducted in the Frozen Soil Engineering Laboratory of Beijing Jiaotong University [28], as shown in Fig 4. The soil samples were conditioned by maintaining a constant temperature of +1˚C before each test, which produced a uniform water content of 17.5%. During the experiment, the sides of the samples were thermally insulated and the temperature at the bottom of each soil column was held at +1˚C. The top surfaces of the soil samples were exposed to a circulating fluid with a temperature of either −1.5˚C or −3˚C. This circulation kept the top surfaces of the soil samples at a constant temperature. After 48 h, the

PLOS ONE
Finite volume-based model for hydrothermal behavior of soil temperatures and moisture contents of the soils were measured as functions of the distances from the bottom of each column.
Numerical simulations of these experiments utilized the soil properties and water migration model parameters given in Tables 2 and 3. All simulations considered a mesh consisting of 50 control volumes with 0.2 cm spacing in the axial direction of the soil column. A fully implicit backward difference time-stepping routine was adopted with a time step of 30 s. The code implementation was deemed convergent if the liquid water contents and temperatures did not change by more than 1% between successive iterations at each time step for all test cases.
In each simulation, the initial and boundary thermal conditions for the heat equation were specified as follows: T (t = 0 s) = 1˚C and T (t > 0 s, y = 10 cm) = −0.75˚C/−3˚C, T (t > 0 s, y = 0) = 1˚C. In addition, the initial and boundary conditions for the water equation were defined by: θ| t = 0 = 0.175, @y @n j t>0 s;y ¼ 0:0 ¼ 0, and @y @n j t>0 s;y ¼ 10 cm ¼ 0. Comparison of experimental and numerical results. The calculated total moisture content and vertical temperature distribution within the soil column were compared against the experimental results to verify the code program proposed.   Table 2. Physical properties of the tested soil sample.

PLOS ONE
Finite volume-based model for hydrothermal behavior of soil

PLOS ONE
Finite volume-based model for hydrothermal behavior of soil moisture contents show reasonable agreement with the measured values. Experiments conducted with a top temperature of −0.75˚C formed a freezing front at the cold end (Fig 5a), and measured moisture contents in the adjacent unfrozen soil always exceeded those in other locations concerning the freezing front. Water diffusivity in the unfrozen zone thus drives fluid flow towards the cold face. Experiments using an upper temperature of −3.0˚C (Fig 4b) produced a different moisture redistribution profile, whereby the measured moisture content was highest near to the cold end (top) of the column. Moving downwards through the column, the moisture content first decreases, then increases, and finally, the freezing point sharply decreases. Consequently, two peak moisture contents occur, one on the top surface of the soil sample and the second located at the position of the freezing edge. Formation of the latter can be explained by the freezing front moving relatively fast compared to diffusive fluid flow, such that the maximum velocity of water migration occurs at the freezing front, although further migration through the column is significantly reduced by the ice formation. These data also demonstrate a strong correlation between the measured and computed soil temperatures (Fig 6), both decrease linearly with depth in each experiment. Fig 7 presents calculated accumulation of ice content in the experiment using a soil-top temperature of −3˚C. In this scenario, when the soil sample was subjected to a sub-freezing temperature, part of the pore water within the cold portion of the sample was the first to solidify into ice. Simultaneously, a gradient in water potential developed in the same direction as the temperature gradient. Water then diffused up from the warm portion to feed the accumulation of pore ice, resulting in water migration and an increase in ice content in the frozen area. At the same time, the liquid water content and the permeability of soil close to the freezing front both decreased. Therefore, it became increasingly difficult for water to flow to the frozen area behind the freezing front. With continued freezing, the freezing front advanced down the soil column to form a frozen area characterized by a relatively homogenous distribution of ice.
The distribution of water and ice within this soil column has also been predicted using our numerical model for boundary conditions ranging between −0.75˚C and −5˚C (Fig 8). These results all show an increase in total moisture content behind the advancing of freezing front, namely, a second peak in moisture content and the depletion in the unfrozen domain. The modeled speed at which the freezing front advances downwards through the column controls the moisture content increase in the frozen layer. For example, a slow-moving freezing front causes a larger increase in moisture content in the frozen soil (Fig 8a). In addition, the extent of the frozen area as a function of depth was modeled to increase as temperatures dropped, and the peak ice content also showed a slight reduction. However, compared with the degree of redistribution of moisture content at different temperatures, predicted ice contents had similar ranges. This indicates that the speed of the advancing freezing front primarily affects liquid water migration. These modeled results are reasonable and have been validated in the literature [29,30]. The close correlation between observed (experimental) and calculated (simulated) results implies that our theoretical formulation can accurately describe the physical processes involved in freeze-thaw weathering.

PLOS ONE
Finite volume-based model for hydrothermal behavior of soil Applications Railway embankment problem. The application of this model was illustrated by simulating the variation of temperature and moisture profiles in a railway embankment over a period. The site belongs to the Baotou-Lanzhou railway in the North China (Fig 9), where the maximum frost depth can reach 1-1.2 m below the natural ground. Field observations show that frost heave exceeds 20 mm in certain area, posing a threat to the safe operation of trains. The hydrothermal method was applied to develop a two-dimensional model of an embankment to evaluate the extent of frost heave defects and related long-term stability of the earth work.
Model and parameters. The spatial finite-volume meshes adopted for modeling of the railway embankment is shown in Fig 10. Here, the upper portion of the mesh is unstructured and the lower part is structured, simulating artificial and natural ground, respectively. The computational domain for both the embankment and natural ground used the soil parameters listed in Tables 2 and 3. This is a defensible assumption since the natural ground was once artificially compacted, and the embankment was constructed using the same soil without

PLOS ONE
Finite volume-based model for hydrothermal behavior of soil geotechnical pre-treatment or modification. The water content of the soil was set to 9%, based on the arid continental climate and deep-seated groundwater table in the region.
Boundary and initial conditions. In Fig 10, the upper boundary EF represents the top surface of embankment fills, and boundaries BE and FI are the embankment's sloped surfaces. Horizontal boundaries AB and IJ represent the natural ground surfaces. Long-term observation of the test site allowed estimation of simplified boundary temperature conditions, as defined by Eq (26) and reported in Table 4.
where T 0 is the mean annual temperature, A is the annual amplitude of temperature variation, p is a vibration cycle equal to 8,760 h, t is the time in hours, and R 0 is the heating rate, taken here to be 0.04˚C/a. Table 4. Thermal boundary conditions used in railway embankment modeling.

Model boundary T 0 (˚C) A (˚C)
Natural ground surfaces (AB and IJ in Fig

PLOS ONE
Finite volume-based model for hydrothermal behavior of soil A constant temperature of T = 7˚C was applied to the bottom surface (LK) of the computational model, which was situated 20 m below the natural ground surface. The left-and righthand boundaries of the model were assumed to be adiabatic.
This numerical simulation assumed that the embankment was constructed on August 15, which typically records the highest annual temperature in this region. The initial temperature profile of the natural ground was obtained by calculating a 50-year transient solution for North China without considering long-term climate changes. The initial temperature distribution within the embankment fills was then determined by the top-surface temperature of natural ground at the point in time when the embankment was constructed. The spatiotemporal characteristics of the temperature and moisture fields within the embankment were then forward-modeled based on its construction 10 years later.

PLOS ONE
Finite volume-based model for hydrothermal behavior of soil Results and analysis. Problems of frost heave in embankments are associated with the frost depth. For simplicity, 0˚C was taken as the soil-freezing temperature in this simulation. The effects of the freeze-thaw process calculated by this model are shown in Fig 11. The upper layer of the embankment is predicted to freeze at the exposed surface on December 11 in the first year after construction, and the maximum freezing depth occurs at around March 27 of the following year. Interestingly, the subsequent thawing process simultaneously propagates from both the top and bottom of the embankment column central line, creating a two-direction mode (Fig 11). The predicted duration of this thawing period is relatively short, lasting approximately 20 days. By contrast, the predicted freezing period lasts from December of the first year to April of the next year, spanning almost 130 days. Fig 12 shows the temperature distribution of the embankment in the 10 th year after construction. During annual cooling associated with decreasing air temperature, the cooling zone of the embankment section is enlarged continuously. The freezing depth on April 1 is situated −1.16 m below the natural ground surface, while the freezing depth of the embankment under the centerline is 1.29 m and 0.13 m thicker than that under the natural ground surface. Moreover, the minimum temperature occurs on the shoulders of the embankment, which suffers from a two-way invasion of cold fronts, causing them to experience the greatest freezing depth. This result implies that preventative measures should be adopted to avoid asymmetrical frost heave and crack formation on embankments, which can otherwise severely affect the stability of railway foundations.
Spatiotemporal moisture redistribution is an important parameter for evaluating the longterm stability of railway embankments. Fig 13 shows the distribution of volumetric moisture content predicted during the freezing period defined above. In the early freezing period, unfrozen water in the shallow layers of the embankment decreases in volume, while free water within the foundation migrates upward towards the frozen fringe (Fig 13a). As a result, moisture accumulated in the frozen fringe, and a layer of ice gradually formed (Fig 13b and 13c). Thereafter, a frozen portion of the modeled soil core appeared beneath the ground surface due to initiation of a two-direction thawing mode (Fig 13d), which shrunk over time. The modeled ice content near to the top and shoulder of the railway embankment increased significantly during this period, although was not predicted to be abundant at the foot of the embankment due to the protective effect of the revetment. Hence, the construction of a revetment is suggested to be a cost-effective measure for preventing frost accumulation, and although additional precautionary measures should be adopted at other locations on the surface of the embankment to prevent ice formation. Critically, the spatiotemporal effect of frost heave within the subgrade can be predicted via our numerical code, purely based on the calculated temperature profile and ice content distribution. Fig 14 shows how volumetric ice content is predicted to vary directly beneath the centerline of the railway embankment section for 10 years after its construction. These data reveal that ice should accumulate progressively closer to the upper part of the embankment over time, suggesting that moisture migration inside the embankment should also continue on an annual basis due to annual temperature changes. As such, moisture-blocking layers should be installed within railway embankments to ensure their long-term stability.

Conclusions
In this paper, a numerical model was presented for the hydrothermal process in unsaturated soils upon freezing and thawing. The model was developed using the open-access OpenFOAM platform. Numerical simulations of laboratory freezing experiments show close correlations between predicted and measured values for spatiotemporal changes in moisture content and temperature, indicating the robustness of our code. Following this validation, this model was applied to predict how temperature and moisture redistribution varies in a railway embankment over time, which are important physical parameters that govern its long-term stability. Our preliminary results demonstrate that the numerical model and associated code presented herein may be useful for design, maintenance, and research on other embankments in seasonally frozen regions.
The data and results presented in this study confirm that our model effectively considers multiple aspects of temperature and fluid flow, and successfully resolves the interdependence between heat and moisture transfer in soil. Future research projects will quantify the magnitude of deformation of soil columns due to freeze-thaw cycles, and field-based data will be collected and integrated into our code to better describe the hydraulic and thermal properties of unsaturated soil.