# Numerical Simulation of Thawing Process in Frozen Soil

Authors: Zhaoyu Yan; Wei Pan; Junjie Fang; Zihui Liu
Journal: Geofluids
 (2020)

---

## Abstract

Permafrost has been thawing faster due to climate change which would release greenhouse gases, change the hydrological regimes, affect buildings above, and so on. It is necessary to study the thawing process of frozen soil. A water-heat coupling model for frozen soil thawing is established on Darcy’s law and Heat Transfer in Porous Media interfaces in Comsol Multiphysics 5.5. Three curves of total liquid water volume, minimum temperature, and total heat flux in the thawing process are obtained from a numerical simulation. The distributions of liquid water, temperature, and pressure based on time are simulated too. The liquid water distribution is consistent with the total liquid water volume curve. The temperature distribution is confirmed by the minimum temperature and total heat flux curve. The pressure distribution represents ice in the frozen soil that generates negative pressure during the melting process. The numerical simulation research in this article deepens the understanding of the internal evolution in the process of frozen soil thawing and has a certain reference value for subsequent experimental research and related applications.


---

## Body

## 1. Introduction

There are large areas of permafrost in boreal and high-elevation regions such as Qinghai-Tibet Plateau [1]. However, hydrological and ecological changes are exacerbating because of climate change [2–7]. Permafrost has been thawing faster affected by the global temperature rise. At the same time, greenhouse gases such as carbon dioxide and methane released by the thawing of permafrost will in turn exacerbate global warming [8–10]. This trend will continue in the foreseeable decades [11]. The distributions of surface water and subface water pathways can be changed by permafrost thaw or changing patterns of seasonal subsurface ice [12–14], so that permafrost thaw makes a big difference to hydrological regimes. In addition, the freeze-thaw cycles of soil can also affect the buildings and facilities above [15, 16]. Thus, soil under freezing and thawing conditions is the key point for studying the impact of climate change on boreal and high-elevation regions [17].It is a mainstream research direction to carry out experimental research for variations of soil under freezing and thawing conditions. Wang et al. [18] studied physicomechanical property changes of Qinghai-Tibet clay due to cyclic freezing and thawing. Zhang et al. [19] studied the variations of the temperatures and volumetric unfrozen water contents during a freezing-thawing process. Wang et al. [20–22] studied fracturing characteristics and mechanical behavior of rock. Darrow et al. [23] measured the zeta potential of cation-treated soils to find its implication on unfrozen water mobility. Cao et al. [24] studied the effects of cyclic freeze-thaw treatments on the fracture characteristics of sandstone. Zhou et al. [25] studied the mechanical behavior of frozen loess under freeze-thaw cycles. Han et al. [26] studied the effect of freeze-thaw cycles on the shear strength of saline soil. Zhao et al. [27–30] studied the dynamic strength of frozen sulfate saline silty clay under cyclic loading. De Guzman et al. [31] discovered that the shear strength of frozen soil samples was reduced by as much as 50% upon freezing and thawing conditions. In addition, some scholars studied the shear behavior of frozen soil-concrete interface under freezing and thawing conditions [32, 33]. Besides the experimental study, numerical simulation is also a feasible approach. Lin et al. [34–36] established the stress field of an open flaw tip under uniaxial compression and shear models. Li et al. [37] developed a THM coupling model to research freeze-thaw failure mechanisms and engineer safeguarding measures. Zhao et al. [38–41] established models about rock behavior which contain the THM model. Grenier et al. [42] studied groundwater flow and heat transport for systems undergoing freeze-thaw. He et al. [43] proposed a coupled heat-moisture-deformation model for saturated frozen soils. Considering the effects of water flow in porous resulted by Soret effect and segregation potential on seepage velocity and water pressure distribution, Tan et al. [44] established a TH coupling model for the construction of Galonglatunnel. Amiri et al. [45] extended the theory of soil freezing curve paradigm. Kurylyk and Watanabe [46] summarized and examined recent models simulating coupled heat and water transport in cold regions.Soil under freeze-thaw conditions has been studied from a variety of perspectives. However, the variation in the frozen soil cannot be displayed intuitively. Because it is difficult to observe the internal evolution of frozen soil during the thawing process, numerical simulation is the approach which is chosen to study the variations of water, temperature, heat, and pressure in the process. The research conclusion has a certain reference value for the subsequent experimental research and related applications.


## 2. Model Building

### 2.1. Model Definition

The numerical simulation of the frozen soil thawing process is carried out by the Comsol Multiphysics 5.5 software. The interfaces selected to simulate the thawing process of frozen soil are Darcy’s law and Heat Transfer in Porous Media.The numerical model is shown in Figure1. A flow channel 3 meters in length and 1 meter in width which is from left to right contains a squared ice inclusion of 0.333 meter in length. The initial temperature of the ice inclusion is Tini=−5°C. The temperature in the porous channel is Tinw=+5°C. The fluid flow is driven by a hydraulic head gradient of ΔH=3% over the length of the channel.Figure 1
Numerical model.In order to facilitate research work, the model is built on these assumptions: no thermal dispersion is included in the heat transfer equation, and water density and dynamic viscosity are considered constant with respect to the temperature (see Table1).Table 1
Material properties of ice, water, and solid matrix.
Material propertyDensity (kg/m3)Heat capacity at constant pressure (J/(kg·K))Thermal conductivity (W/(m·K))Ice92020602.14Water100041820.6Solid matrix26508359


### 2.2. Model Equations

Using the storage model node, Darcy’s law interface contains an implementation of Darcy’s law which explicitly includes an option to define the linearized storageS (SI unit: 1/Pa) using the compressibility of the fluid and the porous matrix:
(1)ρS∂p∂t+∇∙ρ−κμ∇p+ρg∇D=Qm.In this case, the storageS is defined as
(2)S=Swβεp,where Sw is the water saturation, εp is the porosity, and β is the effective compressibility, which is a combined value of water, ice, and solid matrix compressibility. For simplicity, the gravity term in Equation (1) is neglected. The variable Qm is a mass source which represents the additional liquid water due to the melting of the ice inclusion:
(3)Qm=Swεpρi−ρw∂Sw∂t.Here,ρi and ρw are the ice and liquid water densities. The liquid water saturation Sw depends on the phase change as
(4)Sw=Swres+1−Swresθ2.Here,Swres is the residual liquid water saturation, θ2 is a smooth step function defined in the phase change material node. The step function θ2T is zero for temperatures below the melting temperature Tpc, and it equals to 1 for temperatures above Tpc.The heat transfer in the porous medium is described by the following equation:(5)ρCeq∂T∂t+ρCwu∙∇T+∇∙−keq∇T=Q.Here,ρCeq is the equivalent value of density ρ (kg/m3) and heat capacity C at constant pressure (J/kg·K), keq is the effective thermal conductivity (W/m·K), T is the temperature (K), and Q is a heat source (W/m3). The variable Cw is the effective fluid’s heat capacity at constant pressure.


### 2.3. Model Data

The benchmark example describes the ice-to-water phase change of an ice inclusion within a porous channel. The values for thermal properties of water, ice, and solid matrix are presented in Table1.A latent heat of melting,L=334 kJ/kg is used. Additional physical parameters used in the example are summarized in Table 2.Table 2
Additional physical parameters.
Parameterμw (Pa·s)εβ (Pa-1)SwresΩkint (m2)DescriptionDynamic viscosity of waterPorosityEffective compressibilityResidual water saturationImpedance factorIntrinsic permeabilityValue1.793·10−30.3710-80.05501.3·10−10


## 2.1. Model Definition

The numerical simulation of the frozen soil thawing process is carried out by the Comsol Multiphysics 5.5 software. The interfaces selected to simulate the thawing process of frozen soil are Darcy’s law and Heat Transfer in Porous Media.The numerical model is shown in Figure1. A flow channel 3 meters in length and 1 meter in width which is from left to right contains a squared ice inclusion of 0.333 meter in length. The initial temperature of the ice inclusion is Tini=−5°C. The temperature in the porous channel is Tinw=+5°C. The fluid flow is driven by a hydraulic head gradient of ΔH=3% over the length of the channel.Figure 1
Numerical model.In order to facilitate research work, the model is built on these assumptions: no thermal dispersion is included in the heat transfer equation, and water density and dynamic viscosity are considered constant with respect to the temperature (see Table1).Table 1
Material properties of ice, water, and solid matrix.
Material propertyDensity (kg/m3)Heat capacity at constant pressure (J/(kg·K))Thermal conductivity (W/(m·K))Ice92020602.14Water100041820.6Solid matrix26508359


## 2.2. Model Equations

Using the storage model node, Darcy’s law interface contains an implementation of Darcy’s law which explicitly includes an option to define the linearized storageS (SI unit: 1/Pa) using the compressibility of the fluid and the porous matrix:
(1)ρS∂p∂t+∇∙ρ−κμ∇p+ρg∇D=Qm.In this case, the storageS is defined as
(2)S=Swβεp,where Sw is the water saturation, εp is the porosity, and β is the effective compressibility, which is a combined value of water, ice, and solid matrix compressibility. For simplicity, the gravity term in Equation (1) is neglected. The variable Qm is a mass source which represents the additional liquid water due to the melting of the ice inclusion:
(3)Qm=Swεpρi−ρw∂Sw∂t.Here,ρi and ρw are the ice and liquid water densities. The liquid water saturation Sw depends on the phase change as
(4)Sw=Swres+1−Swresθ2.Here,Swres is the residual liquid water saturation, θ2 is a smooth step function defined in the phase change material node. The step function θ2T is zero for temperatures below the melting temperature Tpc, and it equals to 1 for temperatures above Tpc.The heat transfer in the porous medium is described by the following equation:(5)ρCeq∂T∂t+ρCwu∙∇T+∇∙−keq∇T=Q.Here,ρCeq is the equivalent value of density ρ (kg/m3) and heat capacity C at constant pressure (J/kg·K), keq is the effective thermal conductivity (W/m·K), T is the temperature (K), and Q is a heat source (W/m3). The variable Cw is the effective fluid’s heat capacity at constant pressure.


## 2.3. Model Data

The benchmark example describes the ice-to-water phase change of an ice inclusion within a porous channel. The values for thermal properties of water, ice, and solid matrix are presented in Table1.A latent heat of melting,L=334 kJ/kg is used. Additional physical parameters used in the example are summarized in Table 2.Table 2
Additional physical parameters.
Parameterμw (Pa·s)εβ (Pa-1)SwresΩkint (m2)DescriptionDynamic viscosity of waterPorosityEffective compressibilityResidual water saturationImpedance factorIntrinsic permeabilityValue1.793·10−30.3710-80.05501.3·10−10


## 3. Results and Discussion

The simulation runs for 56 hours. After 19.2 hours, the ice completely melted, but the minimum temperature in the channel did not return to the initial 5°C until the 56th hour. With the melting of ice, the changes of various physical indicators in the channel can reflect the internal evolution of the thawing process in frozen soil, which is worthy of further discussion.


### 3.1. Numerical Simulation Curves

With the thawing of frozen soil, the total liquid water volume changed as shown in Figure2. After the simulation started running, the ice in the soil gradually melted, and the total liquid water volume in the channel began to increase. At about 19.2 h, the total liquid water volume in the channel rose to 1.11 and remained constant which demonstrates that the ice has completely melted.Figure 2
Total liquid water volume. Display the variation of total liquid water volume in the thawing process.At the beginning of the simulation, the temperature of the soil in the channel is +5°C, and the temperature of the ice is -5°C; at this time, the minimum temperature in the channel is the temperature of the ice, -5°C. After 56 hours, the minimum temperature in the channel returned to the initial 5°C, and the change curve of the minimum temperature in the channel is shown in Figure3. As the figure shows, the temperature rose rapidly to -1°C in the first two hours after the ice melt. From then on until about 19.2 h, the minimum temperature rose at a very slow rate. In the narrow interval from 19.2 h to 19.9 h, the minimum temperature increased rapidly and exceeded the freezing point. Then, the rise of minimum temperature has experienced three stages: slow, fast, and slow. The minimum temperature in the channel returned to the initial value at 56 h.Figure 3
The minimum temperature. The minimum temperature in the channel as a function of time.In the model, fluid flows from left to right in the channel, so the right side of the rectangle is the outflow boundary of the system. Total heat flux via the outflow boundary is shown as Figure4. Because the ice is in the left half of the channel, there is a certain distance from the right outflow boundary; the heat flux of the outflow boundary on the right side of the channel is not consistent with the heat change caused by ice melting; therefore, there is no obvious correspondence between Figures 3 and 4. The corresponding relationship between heat flux and temperature change will be described in the following content. The change of heat flux at the outflow boundary is driven by ice melting, so the heat flux is negative. The absolute value of total heat flux at the outflow boundary first increased and then decreased.Figure 4
Total heat flux. Total heat flux exiting the system via the outflow boundary as a function of time.


### 3.2. Liquid Water Distribution

The variation of liquid water distribution in the channel is shown in Figure5. In Figure 5, there is the liquid water distribution in numerical simulation for 0 h, 5 h, 10 h, 15 h, 20 h, and 25 h, respectively. It can be seen that the ice begins to melt from the four corners of the square. The left side of the square ice is facing the liquid flow direction, and the melting speed of the two adjacent corners is higher than that of the right side. In the process of melting, the shape of ice gradually changed from square to round, and the scope of ice gradually reduced until it completely melted. It can be seen from the figure at 20 h and 25 h that the ice in the channel has completely melted after 20 hours, the liquid water was full of the channel, and the total liquid water volume is constant with the curve in Figure 2.Figure 5
Liquid water saturation. Effectively, the areas with liquid water appear in blue whereas the areas of ice appear in white.


### 3.3. Temperature Distribution

The temperature distribution in the channel develops with the melting of ice which is shown in Figure6. After the ice melted and heated up, the low temperature diffused around and then moved to the right as the liquid flowed. Finally, the low temperature flowed out of the channel and the initial 5°C was restored in the channel. Obviously, the part below zero is ice; from the figures in the time range of 0 h to 20 h in Figure 6, we can find that the distribution of ice in the melting process inferred from Figure 6 is consistent with Figure 5. Compared with Figure 3, temperature figures at 0 h, 0.083 h, 0.25 h, and 1 h verify the rapid rise of the minimum temperature in the first hour of the simulation. The melting of ice caused the surrounding water temperature to cool down; low-temperature water flowed to the right with the hydraulic head which would change the heat flux of the outflow boundary. The flow of low-temperature water in Figure 6 also explains the variation curve of total heat flux in Figure 4. In Figure 6, there are blue low-temperature water passing through the outflow boundary between 20 h and 30 h, which is in line with the larger absolute value of heat flux in Figure 4.Figure 6
Temperature distribution. Display the variation of temperature distribution in the thawing process.


### 3.4. Pressure Distribution

The pressure distribution in the channel develops with the melting of ice which is shown in Figure7. The streamlines in the figure are Darcy’s velocity field, which represents the flow path of fluid in porous media. It can be seen that the pressure changes in the channel are mainly concentrated on ice. With the beginning of ice melting, the pressure in the square ice area began to decrease. First of all, it changed from positive pressure to negative pressure from four corners and then completely changed to negative pressure in the ice area. The ice melted from a square to a circle, and the negative pressure area also changed to a circle. Then, the area of ice became smaller and smaller, and the absolute value of negative pressure gradually increased until the ice melted completely, and the pressure in the channel remained constant. In addition, the change of streamlines indicates that the water flow path gradually occupied the ice melting area.Figure 7
Pressure distribution. Display the variation of pressure distribution in the thawing process.


## 3.1. Numerical Simulation Curves

With the thawing of frozen soil, the total liquid water volume changed as shown in Figure2. After the simulation started running, the ice in the soil gradually melted, and the total liquid water volume in the channel began to increase. At about 19.2 h, the total liquid water volume in the channel rose to 1.11 and remained constant which demonstrates that the ice has completely melted.Figure 2
Total liquid water volume. Display the variation of total liquid water volume in the thawing process.At the beginning of the simulation, the temperature of the soil in the channel is +5°C, and the temperature of the ice is -5°C; at this time, the minimum temperature in the channel is the temperature of the ice, -5°C. After 56 hours, the minimum temperature in the channel returned to the initial 5°C, and the change curve of the minimum temperature in the channel is shown in Figure3. As the figure shows, the temperature rose rapidly to -1°C in the first two hours after the ice melt. From then on until about 19.2 h, the minimum temperature rose at a very slow rate. In the narrow interval from 19.2 h to 19.9 h, the minimum temperature increased rapidly and exceeded the freezing point. Then, the rise of minimum temperature has experienced three stages: slow, fast, and slow. The minimum temperature in the channel returned to the initial value at 56 h.Figure 3
The minimum temperature. The minimum temperature in the channel as a function of time.In the model, fluid flows from left to right in the channel, so the right side of the rectangle is the outflow boundary of the system. Total heat flux via the outflow boundary is shown as Figure4. Because the ice is in the left half of the channel, there is a certain distance from the right outflow boundary; the heat flux of the outflow boundary on the right side of the channel is not consistent with the heat change caused by ice melting; therefore, there is no obvious correspondence between Figures 3 and 4. The corresponding relationship between heat flux and temperature change will be described in the following content. The change of heat flux at the outflow boundary is driven by ice melting, so the heat flux is negative. The absolute value of total heat flux at the outflow boundary first increased and then decreased.Figure 4
Total heat flux. Total heat flux exiting the system via the outflow boundary as a function of time.


## 3.2. Liquid Water Distribution

The variation of liquid water distribution in the channel is shown in Figure5. In Figure 5, there is the liquid water distribution in numerical simulation for 0 h, 5 h, 10 h, 15 h, 20 h, and 25 h, respectively. It can be seen that the ice begins to melt from the four corners of the square. The left side of the square ice is facing the liquid flow direction, and the melting speed of the two adjacent corners is higher than that of the right side. In the process of melting, the shape of ice gradually changed from square to round, and the scope of ice gradually reduced until it completely melted. It can be seen from the figure at 20 h and 25 h that the ice in the channel has completely melted after 20 hours, the liquid water was full of the channel, and the total liquid water volume is constant with the curve in Figure 2.Figure 5
Liquid water saturation. Effectively, the areas with liquid water appear in blue whereas the areas of ice appear in white.


## 3.3. Temperature Distribution

The temperature distribution in the channel develops with the melting of ice which is shown in Figure6. After the ice melted and heated up, the low temperature diffused around and then moved to the right as the liquid flowed. Finally, the low temperature flowed out of the channel and the initial 5°C was restored in the channel. Obviously, the part below zero is ice; from the figures in the time range of 0 h to 20 h in Figure 6, we can find that the distribution of ice in the melting process inferred from Figure 6 is consistent with Figure 5. Compared with Figure 3, temperature figures at 0 h, 0.083 h, 0.25 h, and 1 h verify the rapid rise of the minimum temperature in the first hour of the simulation. The melting of ice caused the surrounding water temperature to cool down; low-temperature water flowed to the right with the hydraulic head which would change the heat flux of the outflow boundary. The flow of low-temperature water in Figure 6 also explains the variation curve of total heat flux in Figure 4. In Figure 6, there are blue low-temperature water passing through the outflow boundary between 20 h and 30 h, which is in line with the larger absolute value of heat flux in Figure 4.Figure 6
Temperature distribution. Display the variation of temperature distribution in the thawing process.


## 3.4. Pressure Distribution

The pressure distribution in the channel develops with the melting of ice which is shown in Figure7. The streamlines in the figure are Darcy’s velocity field, which represents the flow path of fluid in porous media. It can be seen that the pressure changes in the channel are mainly concentrated on ice. With the beginning of ice melting, the pressure in the square ice area began to decrease. First of all, it changed from positive pressure to negative pressure from four corners and then completely changed to negative pressure in the ice area. The ice melted from a square to a circle, and the negative pressure area also changed to a circle. Then, the area of ice became smaller and smaller, and the absolute value of negative pressure gradually increased until the ice melted completely, and the pressure in the channel remained constant. In addition, the change of streamlines indicates that the water flow path gradually occupied the ice melting area.Figure 7
Pressure distribution. Display the variation of pressure distribution in the thawing process.


## 4. Conclusions

A 2D water-heat coupling model is established to analyse the thawing process of frozen soil. The numerical simulation in this paper reveals the internal evolution of frozen soil in the thawing process from four aspects of water, temperature, heat, and pressure. The main conclusions are as follows:(1)
The variation of liquid water distribution in the thawing process of frozen soil is simulated, which is in full compliance with the total liquid water volume curve, and can reflect the area of unmelted ice(2)
The variation of temperature distribution in the thawing process of frozen soil is simulated. On the one hand, the variation of temperature distribution explains the minimum temperature curve; on the other hand, the minimum temperature water passing through the outflow boundary is consistent with total heat flux(3)
The change of pressure and the flow path of water in frozen soil are simulated. It can be seen that the ice in the frozen soil generates negative pressure during the melting process(4)
Though this numerical simulation has a reference value, it is necessary to carry out more field observation and research on permafrost thawing so that the actual melting mechanism of permafrost can be clearly understood. One potential approach is deploying sensors on a large scale in permafrost to monitor the variation of permafrost thawing



---

*Source: 8822320-2020-10-23.xml*