Thermo-mechanical simulation of frost heave in saturated soils

Abstract

Roads are exposed to various degradation mechanisms during their lifetime. The pavement deterioration caused by the surrounding environment is particularly severe in winter when the humidity and subfreezing temperatures prevail. Frost heave-induced damage is one of the winter-related pavement deterioration. It occurs when the porewater in the soil is exposed to freezing temperatures. The study of frost heave requires conducting a multiphysics analysis, considering the thermal, mechanical, and hydraulic fields. This paper presents the use of a coupled thermo-mechanical approach to simulate frost heave in saturated soils. A function predicting porosity evolution is implemented to couple the thermal and mechanical field analyses. This function indirectly considers the effect of the water seepage inside the soil. Different frost heave scenarios with uniform and non-uniform boundary conditions are considered to demonstrate the capabilities of the method. The results of the simulations indicate that the thermo-mechanical model captures various processes involved in the frost heave phenomenon, such as water fusion, porosity variation, cryogenic suction force generation, and soil expansion. The characteristics and consequences of each process are determined and discussed separately. Furthermore, the results show that non-uniform thermal boundaries and presence of a culvert inside the soil result in uneven ground surface deformations.

Graphical abstract

Keywords

frost heave / multiphysics analysis / thermo-mechanical approach / saturated soils

1 Introduction

Meteorological factors, such as atmospheric temperature and incoming and outgoing radiation, determine the temperature at the ground surface. These parameters fluctuate daily, annually, or even over a longer period. The temperature cycling results in heating and cooling of the ground until a certain depth. The water in the soil skeleton undergoes a phase change because of the freezing and thawing processes. This phenomenon, known as frost action, alters the physical properties of the soil and leads to its volumetric expansion, which manifests itself as upward deformations at the ground surface. The swelling of the soil due to frost action is called frost heave. This phenomenon can cause severe damage to structures, either resting on or buried in the ground, such as roads, bridges, and pipelines [1–3].

Taber [4] reported the earliest attempt in investigating the frost heave phenomenon. He provided valuable facts about the frost heave mechanisms asserting that the water seepage to the freezing front is a decisive factor in frost heave occurrence. The heave that occurs because of water flow to the freezing front is known as secondary heave. It is much greater than primary heave caused by the volumetric expansion of pore water. Mitchell [5] attempted to elucidate the force that draws water to the frozen fringe. He argued that it stems from the freezing-point depression (a drop in the freezing point of porewater), thermal gradient, and thermodynamic condition of the system. Takagi [6] studied the generation of this force by introducing the concept of segregation freezing. Wilen and Dash [7] examined the water flow due to the temperature gradient by considering an isolated crystal of ice. They reported the presence of a thin surface-melted layer at temperatures close to the freezing point and indicated that the water transfer from the surrounding environment to the freezing region is because of the existence of this layer. Rempel et al. [8] also confirmed that the existence of a thin layer of water segregating solid particles from the surrounding ice causes water transfer during frost heave. Notably, different soil types demonstrate different behaviors when subjected to freezing temperatures. In general, fine-grained soils are more sensitive to frost action as they provide a favorable condition for the capillary movement of water in the soil. However, the frost susceptibility of clay soils varies over time because of variations in their hydraulic conductivity. This type of soil has less frost susceptibility in the short term but undergoes a large heave in the long term if the prerequisites are satisfied [9,10].

Various frost heave models have been proposed by implementing either macro- or micro-level analysis. The interaction between various mechanisms, such as heat and mass transfer occurring simultaneously during frost action, adds to the complexity of the problem. This issue leads to a lack of consensus regarding frost heave simulation. The macro-level models neglect or oversimplify some mechanisms, such as ice lens formation involved in frost heave. In contrast, micro-level investigations come up with complicated models requiring a broad range of parameters that must be determined experimentally. Takashi et al. [11] proposed a formula to determine the upward movement of soils by conducting simple frost heave experiments. In this equation, the rate of freezing penetration and the effective stress in the freezing direction are two chief parameters. Zheng and Kanie [12] developed a thermo-hydro-mechanical (THM) approach adopting Takashi’s equation to simulate multidimensional frost heave in saturated soils. They proposed an anisotropic parameter determining the frost heave ratio in different directions. Konrad and Morgenstern [13] proposed a semi-empirical approach that determines frost heave based on the segregation potential (SP) theory. They also developed a finite difference procedure to calculate the amount of heave based on heat and mass transfer [14]. O’Neill and Miller [15] introduced the “rigid ice model” in saturated granular soils, which can predict the ice lens distribution and amount of soil deformation due to frost action. The equations implemented in this model are derived from thermo-mechanical considerations and laboratory investigations. Michalowski [16] proposed a constitutive model to simulate frost heave in saturated soils. He presented a model at the macroscopic level by studying the porosity increase because of ice lens growth. In this approach, the frost susceptibility of soils is embedded in a constitutive property defined by the porosity rate function. Subsequently, Michalowski and Zhu [17] improved this model by modifying the porosity rate function. Li et al. [18,19] established a coupled heat−moisture−deformation framework for studying the frost heave in saturated soils by employing energy, continuity, and equilibrium principles. By including the air phase in the soil, Li et al. [20] further improved this model to simulate frost heave in unsaturated soils. Nishimura et al. [21] proposed a fully coupled THM model to simulate freezing and thawing in saturated soils based on the study of the individual thermal, hydraulic, and mechanical processes, and their various interactions. Liu and Yu [22] proposed a THM model, which can simulate frost heave in unsaturated soils. They implemented Fourier’s law for thermal field analysis, Richard’s equation to consider the water migration, and mechanical constitutive relations as the mechanical field representation. Zhou and Li [23] introduced the separating void ratio as a criterion in a THM framework to evaluate the formation of ice lenses. The segregated pore ice particles start to connect and form ice lenses when the porosity of the soil exceeds the separating void ratio. Yin et al. [24] extended this model and considered the vapor phase in unsaturated soils if vapor exists only in the unfrozen region and condenses to water in the freezing front. Subsequently, Weng et al. [25] improved the THM framework proposed by Zhou and Li [23] acknowledging the inelastic behavior of the freezing soils. In this model for saturated freezing soils, the viscoplasticity of the soil is reflected using Norton–Hoff’s law. By considering the coexistence of water and ice in the frozen fringe, Lai et al. [26] devised a THM model to simulate one-side freezing in saturated soils. This model implements a comprehensive criterion to predict the location and thickness of ice lenses. Zhang and Michalowski [27] developed another constitutive model to simulate frost heave and thaw settlement in saturated soils. The model also determines the changes in the strength of the soil associated with the freezing and thawing processes.

Determining the differential heave is a crucial factor as regards studying damages caused by frost heave in structures either built on or buried in the soil. The constitutive thermo-mechanical frost heave model developed by Michalowski and Zhu [17] can capture non-uniform frost heave and simulate physical processes that takes place inside the soil during frost heave. The other advantage of this method is its compatibility with continuum mechanics, making ‎it possible to be coupled with continuum damage models.‎ In this approach, rather than performing a detailed hydraulic field analysis, A porosity growth function is introduced to account for the effect of water seepage during frost heave. Nevertheless, this feature of the method does not compromise the accuracy.

1.1 Aim and structure of the paper

Although many frost heave models were proposed in the literature, only few works have been conducted on simulating uneven frost heave. This study aims to simulate uneven frost heave and frost action in saturated soils by implementing a thermo-mechanical approach developed by Michalowski and Zhu [17]. Moreover, this approach, is extended to obtain a three-dimensional frost heave model in the present work. It confers the possibility to simulate frost heave under more general conditions (e.g., when the plane strain condition is violated).

The paper first describes the theoretical background of the study. Afterward, the geometries, boundary conditions, and material properties of the generated models are described. Then, the model is calibrated and verified based on experimental frost heave test results reported in the literature. Following this, the results are presented and discussed. Finally, the main conclusions drawn from the simulations and recommendations for future studies are provided.

2 Thermo-mechanical approach

As aforementioned, simulating the frost heave phenomenon requires a multiphasic analysis. In this section, the thermo-mechanical approach for frost heave simulation is explained in detail. First, the energy balance during heat conduction inside the soil is studied. Subsequently, the porosity growth because of water migration and ice formation is discussed. Then, the mechanical equilibrium within the soil is formulated. Finally, we present the soil freezing characteristic curve (SFCC). Notably, in the following sections, bold symbols \mathit{x} represent tensor parameters and the nabla symbol \mathrm{\nabla } denotes the vector differential operator.

2.1 Energy balance

Assuming that energy transfer occurs only by heat conduction while penetration of the frost inside the soil with transient water flow, the energy balance is defined as [28]

(1)C_{\text{p}}\frac{\mathrm{\partial }T}{\mathrm{\partial }t}-L_{\text{f}}{\rho }_{\text{i}}\frac{\mathrm{\partial }{\theta }_{\text{i}}}{\mathrm{\partial }t}=\mathrm{\nabla }\cdot \left(\lambda \mathrm{\nabla }T\right),

where T and C_{\text{p}} are the temperature and the heat capacity per unit volume of soil, respectively, in the order given, L_{\text{f}} is the latent heat of water fusion per unit mass, {\rho }_{\text{i}} and {\theta }_{\text{i}} are the density and volumetric content of ice, respectively, and \lambda is the thermal conductivity of the soil composite.

Considering the soil in a saturated condition, the heat capacity per unit volume of the soil composite can be determined as

(2)C_{\text{p}}=C_{\text{s}}{\theta }_{\text{s}}+C_{\text{w}}{\theta }_{\text{w}}+C_{\text{i}}{\theta }_{\text{i}},

where the volumetric content and volumetric heat capacity of the soil constituents are denoted by \theta and C, respectively. The indices ‘s’, ‘w’, and ‘i’ denote the soil grains, unfrozen water, and ice, in the order given. The schematic phase diagram of a freezing soil is illustrated in Fig.1.

Using the chain rule, the first two terms in Eq. (1) can be combined as

(3)C_{\text{p}}\frac{\mathrm{\partial }T}{\mathrm{\partial }t}-L_{\text{f}}{\rho }_{\text{i}}\frac{\mathrm{\partial }{\theta }_{\text{i}}}{\mathrm{\partial }t}=C_{\text{p}}\frac{\mathrm{\partial }T}{\mathrm{\partial }t}-L_{\text{f}}{\rho }_{\text{i}}\frac{\mathrm{\partial }{\theta }_{\text{i}}}{\mathrm{\partial }T}\frac{\mathrm{\partial }T}{\mathrm{\partial }t}=\frac{\mathrm{\partial }\left(C_{\text{p}}-L_{\text{f}}{\rho }_{\text{i}}{\displaystyle \frac{\mathrm{\partial }{\theta }_{\text{i}}}{\mathrm{\partial }T}}\right)T}{\mathrm{\partial }t}.

Thus, Eq. (1) takes the following form:

(4)C_{\text{a}}\frac{\mathrm{\partial }T}{\mathrm{\partial }t}=\mathrm{\nabla }\cdot \left(\lambda \mathrm{\nabla }T\right),

(5)C_{\text{a}}=C_{\text{p}}-L_{\text{f}}{\rho }_{\text{i}}\frac{\mathrm{\partial }{\theta }_{\text{i}}}{\mathrm{\partial }T},

where C_{\text{a}} is the apparent volumetric heat capacity of the freezing soil.

The thermal conductivity of freezing soil is defined as

(6)\log \lambda =\left(1-n\right)\log {\lambda }_{\text{s}}+{\theta }_{\text{w}}\log {\lambda }_{\text{w}}+{\theta }_{\text{i}}\log {\lambda }_{\text{i}},

or

(7)\lambda ={\lambda }_{\text{s}}^{(1-n)}{\lambda }_{\text{w}}^{{\theta }_{\text{w}}}{\lambda }_{\text{i}}^{{\theta }_{\text{i}}},

where {\lambda }_{\text{s}}, {\lambda }_{\text{w}}, and {\lambda }_{\text{i}} are the thermal conductivity of soil constituents, and n denotes the soil porosity.

The soil temperature at the ground surface is not equal to the air temperature when simulating ground surface deformations due to the penetration of freezing temperature. Newton’s law of cooling must be used to determine the ground surface temperature as

(8)\left(\lambda \mathrm{\nabla }T\right)\cdot \mathit{e}=h_{\text{c}}\left(T_{\text{air}}-T\right),

where h_{\text{c}} is the coefficient of convection heat transfer, T_{\text{air}} is the ambient temperature, and e is the unit vector perpendicular to the boundary plane.

2.2 Porosity evolution

The porosity of soil varies during the frost action inside the soil. This is because of the water flow into the freezing region occurring in response to the generation of cryogenic suction force [29]. Michalowski and Zhu [17] introduced the porosity growth function, in which the water intake to the freezing front and the ice lens growth are considered, to determine the porosity of the soil during frost heave. In other words, instead of conducting an analytical study using the mass conservation law to simulate the water flow and adopting a method to capture the ice lens formation and growth, the product of these processes, namely the porosity evolution, is considered. Accordingly, the application of an appropriate technique to determine the porosity evolution within the soil implicitly guarantees the mass balance in the thermo-mechanical approach.

Michalowski [16] proposed a phenomenological function ascertaining porosity growth in freezing soils. This function was subsequently modified by Michalowski and Zhu [17] so as to achieve better agreement with experimental results. The adjusted porosity growth equation is expressed as

(9)

where n_{\text{m}} and T_{\text{m}} are the maximum rate of porosity evolution and the corresponding temperature, respectively, g_{\text{T}} is the temperature gradient to which the soil sample is subjected to (the ratio n_{\text{m}}/g_{\text{T}} is constant for a given soil), T_0 is the freezing point of water, \zeta is a material parameter accounting for stress dependency, vector l is the direction of heat flow (see Fig.2), and {\sigma }_{\text{J}} is the first invariant of the total stress tensor in the soil, which is defined as

(10){\sigma }_{\text{J}}={\sigma }_{ii}+{\sigma }_{jj}+{\sigma }_{kk},

where i, j, and k are principle directions.

The heat flow direction is obtained by

(11)\mathit{l}=\frac{\mathrm{\nabla }T}{\Vert \mathrm{\nabla }T\Vert }.

Thus, the temperature gradient along the heat flow direction for the given point is determined according to

(12)\frac{\mathrm{\partial }T}{\mathrm{\partial }\mathit{l}}=\mathrm{\nabla }T\cdot \mathit{l}.

2.3 Mechanical equilibrium

Using Newton’s second law, the equation of stress equilibrium in soil is expressed as [30]

(13)\rho \frac{{\mathrm{\partial }}^2\mathit{u}}{\mathrm{\partial }{t}^2}=\mathrm{\nabla }\cdot \mathit{\sigma }+{\mathit{f}}_{\text{v}},

where u and {\mathit{f}}_{\text{v}} are the displacement and body force vectors, respectively, \rho is the density of soil composition, and \mathit{\sigma } is the 2nd-order Cauchy stress tensor, which is determined as

(14)\mathit{\sigma }=\mathit{C}(E,v):{\mathit{\epsilon }}_{\text{el}},

where {\mathit{\epsilon }}_{\text{el}} is the 2nd-order elastic strain tensor and C is the 4th-order stiffness tensor whose components are functions of the Poisson’s ratio v and elastic modulus E.

The decomposition of the mechanical response of materials is an effective method for simulating their mechanistic behavior [31–33]. In this study, the total strain is additively decomposed into the strain induced by the frost heave {\mathit{\epsilon }}_{\text{fh}} and the elastic strain, such that

(15)\mathit{\epsilon }={\mathit{\epsilon }}_{\text{el}}+{\mathit{\epsilon }}_{\text{fh}}.

The strain increment is determined by differentiating Eq. (15) as

(16)\text{d}\mathit{\epsilon }=\text{d}{\mathit{\epsilon }}_{\text{el}}+\text{d}{\mathit{\epsilon }}_{\text{fh}}.

The strain increment caused by frost heave \text{d}{\mathit{\epsilon }}_{\text{fh}} is calculated based on the porosity evolution inside the soil as [17]

(17)\left[\begin{array}{c}\text{d}{\epsilon }_{\text{fh,11}}\\ \text{d}{\epsilon }_{\text{fh,22}}\\ \text{d}{\epsilon }_{\text{fh,33}}\end{array}\right]=\left[\begin{array}{c}\xi \text{d}n\\ \frac{1}{2}(1-\xi )\text{d}n\\ \frac{1}{2}(1-\xi )\text{d}n\end{array}\right],

where \xi is the expansion coefficient (\xi =1/3 for isotropic deformation and \xi =1 when one-dimensional growth occurs). Furthermore, while deriving the equation above, the coordinate system is assumed such that direction 1 locally coincides with the heat flow direction.

By implementing Hooke’s law, the elastic strain increment is expressed as

(18)\left[\begin{array}{c}\text{d}{\epsilon }_{\text{el,}xx}\\ \text{d}{\epsilon }_{\text{el,}yy}\\ \text{d}{\epsilon }_{\text{el,}\mathit{\text{z}}\mathit{\text{z}}}\\ \text{d}{\epsilon }_{\text{el,}xy}\\ \text{d}{\epsilon }_{\text{el,}y\mathit{\text{z}}}\\ \text{d}{\epsilon }_{\text{el,}\mathit{\text{z}}x}\end{array}\right]=\left[\begin{array}{c}\frac{1}{E}\left(\text{d}{\sigma }_{xx}-v(\text{d}{\sigma }_{yy}+\text{d}{\sigma }_{\mathit{\text{z}}\mathit{\text{z}}})\right)\\ \frac{1}{E}\left(\text{d}{\sigma }_{yy}-v(\text{d}{\sigma }_{xx}+\text{d}{\sigma }_{\mathit{\text{z}}\mathit{\text{z}}})\right)\\ \frac{1}{E}\left(\text{d}{\sigma }_{\mathit{\text{z}}\mathit{\text{z}}}-v(\text{d}{\sigma }_{xx}+\text{d}{\sigma }_{yy})\right)\\ \frac{(1+\upsilon )\text{d}{\tau }_{xy}}{E}\\ \frac{(1+\upsilon )\text{d}{\tau }_{y\mathit{\text{z}}}}{E}\\ \frac{(1+\upsilon )\text{d}{\tau }_{\mathit{\text{z}}x}}{E}\end{array}\right],

where x, y, and z are axes of an arbitrary Cartesian coordinate system.

The plane strain condition is common while determining ground surface deformations due to frost heave. Considering this, the frost heave strain increment in any arbitrary direction can be determined using the transformation rule as

(19)\left[\begin{array}{c}\text{d}{\epsilon }_{\text{fh,}xx}\\ \text{d}{\epsilon }_{\text{fh,}yy}\\ \text{d}{\epsilon }_{\text{fh,}xy}\end{array}\right]=\left[\begin{array}{ccc}\text{co}{\text{s}}^2\varphi & \text{si}{\text{n}}^2\varphi & -\text{cos}\varphi \text{sin}\varphi \\ \text{si}{\text{n}}^2\varphi & \text{co}{\text{s}}^2\varphi & \text{cos}\varphi \text{sin}\varphi \\ 2\text{sin}\varphi \text{cos}\varphi & -2\text{sin}\varphi \text{cos}\varphi & \text{co}{\text{s}}^2\varphi -\text{si}{\text{n}}^{\text{2}}\varphi \end{array}\right]\left[\begin{array}{c}\text{d}{\epsilon }_{\text{fh,}11}\\ \text{d}{\epsilon }_{\text{fh,}22}\\ 0\end{array}\right],

where \varphi is the angle between the heat flow direction and axis x.

The strain increment induced by frost heave is calculated in a local coordinate system associated with the direction of heat flow (see Fig.3). Therefore, it must be transferred to the global coordinate system using the following transformation rule below to be used in a three-dimensional simulation.

(20)\text{d}{\mathit{\epsilon }}_{\text{FH}}=R_{3,(\frac{\text{π}}{2}-\alpha )}\cdot \left(R_{2,\omega }\cdot \text{d}{\mathit{\epsilon }}_{\text{fh}}\cdot R_{2,\omega }^{\text{T}}\right)\cdot R_{\text{3,(}\frac{\text{π}}{\text{2}}-\alpha )}^{\text{T}},

where \text{d}{\mathit{\epsilon }}_{\text{FH}} and \text{d}{\mathit{\epsilon }}_{\text{fh}} are the strain increments induced by frost heave in the global and local coordinate systems, respectively, \alpha is the angle between the project of the heat flow direction in the x–y plane and axisy, and \omega is the angle between the heat flow direction and x–y plane.

The rotation matrices (R_{\text{3,(}\frac{\text{π}}{\text{2}}-\alpha \text{)}}, R_{\text{2,}\omega }) are defined as

(21)R_{3,(\frac{\text{π}}{2}-\alpha )}=\left[\begin{array}{ccc}\text{sin}\alpha & \text{cos}\alpha & 0\\ -\text{cos}\alpha & \text{sin}\alpha & 0\\ 0& 0& 1\end{array}\right],R_{2,\omega }=\left[\begin{array}{ccc}\text{cos}\omega & 0& \text{sin}\omega \\ 0& 1& 0\\ -\text{sin}\omega & 0& \text{cos}\omega \end{array}\right].

2.4 Soil freezing characteristic curve

As the freezing temperature penetrates the soil, only a portion of the water in the soil voids freezes and the other part remains in the liquid state at a temperature below the freezing point. The analogy of drying and freezing processes in the soil (in both cases water is removed from soil voids) can help define SFCC based on the concept of soil−water characteristic curve (SWCC) [34,35]. SFCC describes the relationship between the amount of unfrozen water and the cryogenic force during frost action inside the soil. In this approach, the negative suction in unsaturated soils is treated as the cryogenic force in freezing soils. The cryogenic suction (\psi) can be obtained as [36]

(22)\psi =-L_{\text{f}}{\rho }_{\text{w}}\text{ln}\left(\frac{T+273.15}{T_0+273.15}\right),

where {\rho }_{\text{w}} is the water density. It is worth mentioning that the pore ice pressure is assumed to be equal to the atmospheric pressure and the solute effect is neglected while determining the cryogenic suction.

Considering the van Genuchten SWCC model [37], the amount of unfrozen water while frost penetration into the soil is expressed as

(23){\theta }_{\text{w}}={\theta }_{\text{r}}+\left({\theta }_{\text{sat}}-{\theta }_{\text{r}}\right){\left(1+{\left(\eta \psi \right)}^{\beta }\right)}^{\frac{1-\beta }{\beta }},

where {\theta }_{\text{r}} and {\theta }_{\text{sat}} represent the residual and the saturated water contents, respectively, and \eta and \beta are fitting parameters.

3 Numerical simulations

In this section, the geometry, thermal and mechanical boundary conditions as well as initial states of the models are explained first. Afterward, the properties of the materials implemented in the models are presented.

3.1 Geometry, boundary, and initial conditions

The finite element software COMSOL Multiphysics^® 5.6 was used to solve the coupled system of differential equations representing the thermal field, mechanical field, and porosity growth during frost action inside the soil. The implicit backward differentiation formula (BDF) with variable orders was used as a time-stepping method in this study.

The governing equations, detailed in the previous section, are defined as

(24)\{\begin{array}{rl}& C_{\text{p}}\frac{\mathrm{\partial }T}{\mathrm{\partial }t}-L_{\text{f}}{\rho }_{\text{i}}\frac{\mathrm{\partial }{\theta }_{\text{i}}}{\mathrm{\partial }t}=\mathrm{\nabla }\cdot \left(\lambda \mathrm{\nabla }T\right),\\ & \frac{\text{d}n}{\text{d}t}=n_{\text{m}}{\left(\frac{T-T_0}{T_{\text{m}}}\right)}^2\exp \left(1-{\left(\frac{T-T_0}{T_{\text{m}}}\right)}^2\right)\frac{\left|{\displaystyle \frac{\mathrm{\partial }T}{\mathrm{\partial }\mathit{l}}}\right|}{g_{\text{T}}}\\ & \cdot \exp \left(-\frac{\left|{\sigma }_{\text{kk}}\right|}{\zeta }\right)\exp \left(-\frac{\mathrm{\partial }{\theta }_{\text{i}}}{\mathrm{\partial }{\theta }_{\text{w}}}\right),\\ & \mathrm{\nabla }\cdot \mathit{\sigma }+{\mathit{f}}_{\text{v}}=0.\end{array}

In simulation (A), the analysis was carried out within a soil domain with an initial height of 300 mm and a square section of side 100 mm. Fig.4 shows the thermal and mechanical boundary conditions. The displacement is free only along the z axis (the heat flow direction), while in two other perpendicular directions, the translational displacements are restricted. The cooling process begins on the upper side, where the temperature drops to −5 °C. On the bottom side, it is maintained at the initial temperature, which is set to 1 °C. The side faces are assumed to be fully insulated, such that

(25)\mathrm{\nabla }T\cdot \mathit{e}=0.

Simulation (B) intends to capture the uneven deformation of the ground surface with non-uniform thermal boundary conditions by conducting a plane strain analysis. In frost heave scenario (C), a three-dimensional simulation is carried out to determine the frost heave around a culvert buried in the soil. Fig.5 shows the thermal and mechanical boundary conditions for study scenarios (B) and (C). In these simulations, the created models are subjected to Stockholm temperature from mid-December 2009 until mid-January 2010 (see Fig.6) [38]. During this period, the thawing effect is negligible.

3.2 Material properties, calibration of the soil freezing characteristic curve, and verification of the model

Tab.1 presents the mechanical and thermal characteristics of the soil. They were extracted by Williams and Smith [39] and Selvadurai et al. [40]. As shown in Fig.7, the elastic modulus of the soil composite is temperature-dependent [41].

A numerical simulation was carried out based on the experimental test conducted by Fukuda et al. [42] to verify the model. In this test, a cylindrical soil specimen with a radius of 50 mm and a height of 70 mm was exposed to a controlled temperature gradient. While the initial temperature of the sample was 5 °C, the test started by abruptly reducing the temperature on the bottom side to −5 °C.

As illustrated in Fig.8, the simulated curve and experimental data indicate a good level of matching. In addition, the SFCC is calibrated using the data in the work by Fukuda et al. [42], wherein the amount of unfrozen water at subfreezing temperatures was reported (see Fig.9). The parameters to be used in SFCC and porosity rate function are provided in Tab.2 and Tab.3.

4 Results and discussion

In this section, the outcomes of the analyses are presented. Simulation (A) elaborates on the whole process of frost action within the soil, while simulations (B) and (C) determine the differential ground surface deformations due to frost heave.

4.1 Frost heave scenario (A), frost action study, and uniform frost heave

The freezing process starts in regions of the soil within which the temperature falls below 0 °C. Therefore, as a prerequisite for frost action study, determining the temperature profile inside the soil in advance is necessary. Fig.10 shows the evolution of the temperature distribution inside the soil. The heat conduction inside the soil tends to attain a steady-state since the temperature at the upper side of the model is kept constant after reaching −5 °C. As the temperature falls below 0 °C, the phase transition of the water from liquid to solid (ice) begins. This process is a type of transition to a less energetic state since the atoms in liquid are moving much faster than in solids. The energy is released while turning the water into ice. This can be observed in Fig.11, wherein the apparent heat capacity increases while ice formation. However, Fig.11 shows the energy increase caused by the phase change of the water does not occur only at a certain level where the temperature is 0 °C, rather it happens within a larger region. The interpretation of this phenomenon is that the porewater gradually turns into ice at subfreezing temperatures and water fusion does not only take place at a temperature of 0 °C. This can also be inferred by considering Fig.12 where the ice content is illustrated at different time steps. The growth of ice content even on the upper side of the model at different time steps confirms the gradual fusion of water.

Fig.13 shows the amount of cryogenic suction. This negative force is generated because of the transformation of liquid porewater to ice. It can be observed that with the penetration of the freezing point into the soil, the amount of suction force and the area in which this force is generated simultaneously increase. Water seepage from the surrounding environment to the freezing area is the main ramification of cryogenic suction force generation. The availability of porewater in liquid form at temperatures far below 0 °C also confirms the water migration, as shown in Fig.14. Notably, the water flows through capillary channels inside the soil under cryogenic suction. In soils with larger grains where these channels are seldom, the water flow to frozen fringe is not possible. The grain size is a decisive factor in assessing frost susceptibility of soils.

The water migration to the freezing region as well as water fusion yields in porosity growth inside the soil. Fig.15 illustrates the porosity evolution during frost heave. The expected consequence of frost action will be the increase of the soil volume since the volume of the soil is a function of its porosity. The upward deformation of the soil due to porosity growth and freezing front penetration is presented in Fig.16. Furthermore, the volumetric expansion of the soil occurs in all directions even in one-dimensional heat flow conditions. However, the assumed mechanical boundary conditions restrict lateral strains. This matter leads to lateral stress development inside the soil, as shown in Fig.17. The amount of stress exceeds the yield stress of the unfrozen soil at some points. This is because the modulus of elasticity of the frozen soil is significantly higher than that of the soil in unfrozen conditions.

4.2 Frost heave scenarios (B) and (C) and non-uniform frost heave simulation

As aforementioned, in simulation (B), a region with non-uniform boundaries is subjected to freezing temperatures. Fig.18 illustrates the heat flow directions and the temperature field at a certain time. In this figure, it is observed that the heat flow direction varies inside the region because of non-uniform boundaries. As previously explained, the strain increment induced by frost heave is larger in the direction of frost penetration than in other directions. The consequence of this characteristic of frost heave is illustrated in Fig.19, where within region (2) the amount of upward deformation is larger at places with the vertical heat flow direction. Fig.19 also indicates that the amount of upward deformations occurred in region (2) are greater compared to region (1). It is because the area experiencing subfreezing temperatures in region (2) is larger than that region (1), as shown in Fig.18(a).

The heat flow directions and temperature distribution at the last time step in simulation (C) are illustrated in Fig.20. It shows that around the culvert, the direction of freezing temperature penetration is not vertical, and the rate of freezing temperature penetration above the culverts is greater compared to other parts. This is because less thermal energy is required for heat conduction in the vicinity of the culvert because of insulation. Fig.21 illustrates ground surface deformations at different time steps. Because of faster heat conduction, the ground above the culvert experiences a “bump”. However, the effect of the heat conduction rate abates, and the frost heave grows faster in other places than above the culvert. This is the time when the effect of the cavity, inside which frost heave does not occur, begins to become dominant. This can be observed in Fig.21, where the height of the bump above the culvert decreases over time.

5 Conclusions and future work

A thermo-mechanical framework has been implemented to simulate the frost heave phenomenon. In this approach, the porosity growth function was used to couple the thermal and mechanical fields. This function implicitly considers the water flow and ensures mass conservation during the heaving process. In other words, rather than simulating the water seepage to the frozen region inside the soil, the consequence of this action, namely porosity evolution, is considered.

The capabilities of the thermo-mechanical model in simulating frost heave were presented using three different frost heave scenarios. In simulation (A), where one-dimensional heat flow occurs, the whole process of the frost action within the soils was studied. It was shown that subjecting a saturated soil into freezing temperature results in water fusion and generation of a negative pressure, which is known as cryogenic suction. This force, in the presence of capillary tubes inside the soil, drives the water from unfrozen parts of the soil to the frozen region where it turns into ice. The seepage of water feeds the formation and growth of ice lenses and leads to swelling of the soil by increasing its porosity.

In simulations (B) and (C), the non-uniform frost heave was investigated. It was shown that the non-uniform thermal boundaries and the presence of the culvert inside the soil yields in uneven ground surface deformations. This is either because of the varying heat flow directions or different rates of the freezing temperature penetration since frost heave is a function of the temperature gradient and it is greater along the heat flow direction.

As regards future research, conducting experimental investigations to calibrate porosity rate function for other soil types with different initial porosities is highly recommended. This would afford the opportunity to study the effect of the initial porosity variation inside the soil on frost heave.

References

[1] Chen X, Shi D. Influence of freeze–thaw cycles on apparent dynamic tensile strength, apparent dynamic fracture toughness and microstructure of concrete under impact loading. European Journal of Environmental and Civil Engineering, 2019, 25(11): 1977–2001

[2] Sylvestre O, Bilodeau J P, Doré G. Effect of frost heave on long-term roughness deterioration of flexible pavement structures. International Journal of Pavement Engineering, 2019, 20(6): 704–713

[3] Sun B, Yang Y, Li X, Yan M, Xie M, Bao Y. Full-scale investigation of post-tensioned prestressed concrete bridge girders subjected to frost heaving in cold regions. Engineering Structures, 2022, 250: 113413

[4] Taber S. Frost heaving. Journal of Geology, 1929, 37(5): 428–461

[5] MitchellJ K. Fundamentals of the Soil Behavior. New York: John Wiley & Sons, 1976

[6] TakagiS. Segregation freezing as the cause of suction force for ice lens formation. Engineering Geology, 1979, 13(1−4): 93−100

[7] Wilen L A, Dash J G. Frost heave dynamics at a single crystal interface. Physical Review Letters, 1995, 74(25): 5076–5079

[8] Rempel A W, Wettlaufer J S, Worster M G. Interfacial premelting and the thermomolecular force: Thermodynamic buoyancy. Physical Review Letters, 2001, 87(8): 088501

[9] DagliD. Laboratory Investigations of Frost Action Mechanisms in Soils. Luleå: Luleå University of Technology, 2017

[10] ChamberlainE J. Frost Susceptibility of Soil, Review of Index Tests. Hanover: Cold Regions Research and Engineering Lab, 1981

[11] TakashiTYamamotoHOhraiTMasudaM. Effect of penetration rate of freezing and confining stress on the frost heave ratio of soil. In: Proceedings of the 3rd International Conference on Permafrost. Proceedings of the 3rd International Conference on Permafrost, 1978

[12] Zheng H, Kanie S. Combined thermal-hydraulic-mechanical frost heave model based on Takashi’s equation. Journal of Cold Regions Engineering, 2015, 29(4): 04014019

[13] Konrad J M, Morgenstern N R. The segregation potential of a freezing soil. Canadian Geotechnical Journal, 1981, 18(4): 482–491

[14] Konrad J M, Morgenstern N R. Frost heave prediction of chilled pipelines buried in unfrozen soils. Canadian Geotechnical Journal, 1984, 21(1): 100–115

[15] O’Neill K, Miller R D. Exploration of a rigid ice model of frost heave. Water Resources Research, 1985, 21(3): 281–296

[16] Michalowski R L. A constitutive model of saturated soils for frost heave simulations. Cold Regions Science and Technology, 1993, 22(1): 47–63

[17] Michalowski R L, Zhu M. Frost heave modeling using porosity rate function. International Journal for Numerical and Analytical Methods in Geomechanics, 2006, 30(8): 703–722

[18] Li N, Chen B, Chen F, Xu X. The coupled heat−moisture−mechanic model of the frozen soil. Cold Regions Science and Technology, 2000, 31(3): 199–205

[19] Li N, Chen F, Su B, Cheng G. Theoretical frame of the saturated freezing soil. Cold Regions Science and Technology, 2002, 35(2): 73–80

[20] Li N, Chen F, Xu B, Swoboda G. Theoretical modeling framework for an unsaturated freezing soil. Cold Regions Science and Technology, 2008, 54(1): 19–35

[21] Nishimura S, Gens A, Olivella S, Jardine R J. THM-coupled finite element analysis of frozen soil: Formulation and application. Geotechnique, 2009, 59(3): 159–171

[22] Liu Z, Yu X. Coupled thermo-hydro-mechanical model for porous materials under frost action. Acta Geotechnica, 2011, 6(2): 51–65

[23] Zhou J, Li D. Numerical analysis of coupled water, heat and stress in saturated freezing soil. Cold Regions Science and Technology, 2012, 72: 43–49

[24] Yin X, Liu E, Song B, Zhang D. Numerical analysis of coupled liquid water, vapor, stress and heat transport in unsaturated freezing soil. Cold Regions Science and Technology, 2018, 155: 20–28

[25] Weng X, Sun Y, Yang Z, Wang D, Yu H. A thermal-hydromechanical coupled model for poro-viscoplastic saturated freezing soil. European Journal of Environmental and Civil Engineering, 2020, 26(9): 4220–4236

[26] Lai Y, Pei W, Zhang M, Zhou J. Study on theory model of hydro-thermal-mechanical interaction process in saturated freezing silty soil. International Journal of Heat and Mass Transfer, 2014, 78: 805–819

[27] Zhang Y, Michalowski R L. Thermal-hydro-mechanical analysis of frost heave and thaw settlement. Journal of Geotechnical and Geoenvironmental Engineering, 2015, 141(7): 04015027

[28] ÇengelY A. Heat and Mass Transfer: A Practical Approach. New York: McGraw-Hill, 2007

[29] Hohmann M. Soil freezing—The concept of soil water potential. State of the art. Cold Regions Science and Technology, 1997, 25(2): 101–110

[30] BowerA. Applied Mechanics of Solids. Boca Raton: CRC Press, 2009

[31] Vu-Bac N, Duong T X, Lahmer T, Areias P, Sauer R A, Park H S, Rabczuk T. A NURBS-based inverse analysis of thermal expansion induced morphing of thin shells. Computer Methods in Applied Mechanics and Engineering, 2019, 350: 480–510

[32] Vu-Bac N, Rabczuk T, Park H S, Fu X, Zhuang X. A NURBS-based inverse analysis of swelling induced morphing of thin stimuli-responsive polymer gels. Computer Methods in Applied Mechanics and Engineering, 2022, 397: 115049

[33] Areias P, Vu-Bac N, Rabczuk T. A multisurface constitutive model for highly cross-linked polymers with yield data obtained from molecular dynamics simulations. International Journal of Mechanics and Materials in Design, 2018, 14(1): 21–36

[34] Ren J, Vanapalli S K, Han Z. Soil freezing process and different expressions for the soil-freezing characteristic curve. Sciences in Cold and Arid Regions, 2018, 9(3): 221–228

[35] Koopmans R W R, Miller R D. Soil freezing and soil water characteristic curves. Soil Science Society of America Journal, 1966, 30(6): 680–685

[36] Kurylyk B L, Watanabe K. The mathematical representation of freezing and thawing processes in variably-saturated, non-deformable soils. Advances in Water Resources, 2013, 60: 160–177

[37] van Genuchten M T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal, 1980, 44(5): 892–898

[38] MobergA. Stockholm Historical Weather Observations—Daily Mean Air Temperatures Since 1756. Stockholm: Bolin Centre Database, 2020

[39] WilliamsP JSmithM V. The Frozen Earth, Fundamentals of Geocryology. Cambridge: Cambridge University Press, 1989

[40] Selvadurai A P S, Hu J, Konuk I. Computational modeling of frost heave induced soil–pipeline interaction: I. Modeling of frost heave. Cold Regions Science and Technology, 1999, 29: 215–228

[41] LadanyiBShenM. Freezing pressure development on a buried chilled pipeline. In: Proceedings of International Symposium on Frost Geotechnical Engineering. Rotterdam: A.A. Balkema, 1993

[42] FukudaMKimHKimY. Preliminary results of frost heave experiments using standard test sample provided by TC8. In: Proceedings of Ground Freezing 97: Frost Action in Soils. Boca Raton: CRC Press, 1997

RIGHTS & PERMISSIONS

Higher Education Press