Experimental study and numerical modeling of the thermo-hydro-mechanical processes in soil freezing with different frost penetration directions

This research work presents an experimental and numerical study of the coupled thermo-hydro-mechanical (THM) processes that occur during soil freezing. With focusing on the artificial ground freezing (AGF) technology, a new testing device is built, which considers a variety of AGF-related boundary conditions and different freezing directions. In the conducted experiments, a distinction is made between two thermal states: (1) The thermal transient state, which is associated with ice penetration, small deformations, and insignificant water suction. (2) The thermal (quasi-) steady state, which has a much longer duration and is associated with significant ice lens formation due to water suction. In the numerical modeling, a special focus is laid on the processes that occur during the thermal transient state. Besides, a demonstration of the micro-cryo-suction mechanism and its realization in the continuum model through a phenomenological retention-curve-like formulation is presented. This allows modeling the ice lens formation and the stiffness degradation observed in the experiments. Assuming a fully saturated soil as a biphasic porous material, a phase-change THM approach is applied in the numerical modeling. The governing equations are based on the continuum mechanical theory of porous media (TPM) extended by the phase-field modeling (PFM) approach. The model proceeds from a small-strain assumption, whereas the pore fluid can be found in liquid water or solid ice state with a unified kinematics treatment of both states. Comparisons with the experimental data demonstrate the ability and usefulness of the considered model in describing the freezing of saturated soils.


Introduction
The need for efficient, environmentally friendly and modern infrastructures, especially in big cities with a high population density, has led to an increase in the dependency on underground constructions such as underground traffic routes. Moreover, in areas with difficult geological conditions, such as hills, mountains, and river crossings, relying on underground constructions, such as tunnels, provides feasible economical and ecological solutions. In the construction of underground structures especially in inner-city areas, artificial ground freezing (AGF) is a widely used technology. It utilizes the advantages of frozen soils, which can bear static loading and have a waterproofing functionality. Nevertheless, the accompanying volumetric expansions might lead to soil deformations and frost heave at the ground surface and, thus, damaging the nearby constructions. This occurs frequently in finegrained soils due to ice lens formation. In particular, ice lenses are pure ice layers in the soil, which evolve at negative temperatures near the frost line by the accumulation of sucked groundwater. This suction occurs during the phase-change process in the partially frozen area (frozen fringe) and leads to water transport from the unfrozen area toward the ice lenses (see, e.g., [7,21,57,96,102,110,117]). Different case studies in connection with projects related to the AGF technique could be found in [3,12,74].
Early studies of ice lens formation in open systems were carried out by [109,110] who investigated the frost penetration into the soil due to low ambient temperatures. He considered different types of soils, investigated the effect of overburden pressure and studied several factors that contribute to the water suction. Additionally, Taber [109] pointed out that the grain sizes, the water availability, the size and percentage of voids, as well as the cooling rate play a major role in the frost heave and the ice lens formation. To quantify the frost heave onset due to soil freezing, tests on a macroscopic level were carried out by, e.g., [4,55,57,60,62,[64][65][66]69] and [93]. Moreover, microscopic tests and studies using X-ray images [5,61,92], particle image velocimetry [9,29], and time lapse photography using fluorescent tracer [118] helped in better understanding the formation process of discrete ice lenses and the frozen fringes characteristics.
While most of the studies on soil freezing that can be found in the literature are motivated by weather-induced freezing and considered the corresponding boundary conditions, they overlooked the differences in the soil behavior between naturally induced and artificially induced soil freezing. Of particular importance in this regard are the effects of freezing direction and the different thermal boundary conditions. In this work, frost heave experiments are carried out, which include an investigation of the variation of the thermal state during AGF, the different freezing directions, the temperature conditions, the water pressure change and the overburden loading. Therefore, the outlined results in this paper focus on investigations of thermal states and freezing directions at different temperature boundary conditions. For an appropriate evaluation of the experimental results, the different thermal states during freezing have to be considered. The ice lenses can initiate during the thermal transient state. However, with the decrease in the frost penetration rate, the ice lens thickness increases until the frost penetration rate becomes almost zero and the thermal (quasi-) steady state is reached. At a given water availability, the heave due to ice lens formation during the thermal (quasi-) steady state can significantly exceed the heave due to ice formation during the thermal transient state. The consideration of the thermal (quasi-) steady state is important since the artificial ground freezing is usually used for long construction periods. For this reason, both thermal states are considered separately in the analysis of frost heave tests. Previous studies on frost heave with ice lens formation mostly refer to the natural frost penetration in cold regions and, therefore, focus on the transient state [63].
Due to gravity effects, soil freezing is direction-dependent, which is addressed in the underlying study. For instance, in tunnel construction with horizontal or slightly inclined freezing pipes, ice accumulates around the pipes, however, more rapidly upward as the gravity enhances the water flow toward the freezing area. In this connection, it is worth mentioning that [8,25] and Niggemann [81] reported increasing permeability in the vertical direction as a result of longitudinal shrinkage cracks.
In the underlying work, a numerical modeling framework that takes the multi-physical thermo-hydro-mechanical (THM) processes of ground freezing into account is presented. In this, a saturated soil is treated as a nonisothermal, deformable, fully saturated biphasic porous material with a single fluid that can change depending on the thermal conditions between a solid ice and a liquid water state. The model is based on a coupled phase-fieldporous media approach [108], where the main focus is laid on the temperature-driven processes that lead to the phase transition between water and ice and the freezing-related deformations that can be observed in the experiments. The governing equations of the macroscopic model are based on the well-founded theory of porous media (TPM) extended by the phase-field modeling (PFM), whereas the combined PFM-TPM approach is implemented within a finite element method (FEM) framework to reliably simulate the coupled process of soil freezing.
The macroscopic TPM [19,32,36,40,70,71], which extends the Theory of Mixtures (TM) [18,113] by the concept of volume fractions is used as a continuum-mechanical framework for the description of the deformable, heterogeneous porous solid. For a more detailed discussion of the general frame of the TPM, its extensions, solution schemes, applications and historical development, the interested reader is referred to [20,24,31,35,36,49,50,52,71,85,101] among others. Besides, the phase-field model (PFM) is employed for deriving a unified kinematics formulation for the phase changing fluid constituent. Details and references related to the energetic PFM for phase-change materials description can be found in, e.g., [16,59,105,115].
In the literature, different THM models within porous media frameworks have been developed to capture the multi-physical phenomena related to soil freezing, see, e.g., [28,56,68,79,80,82]. Additionally, several THM numerical models were developed within the AGF context, to capture the volume expansion of water transforming into ice and the contribution of the micro-cryo-suction mechanism responsible for driving liquid water toward the frozen zone, see, e.g., [23,112,124,125]. Related to the physical aspects and modeling of the ice lens formation during the freezing process, several theories and models that utilize different ice lens formation criteria were proposed in [10,67,76,77,88,123]. Recently, Gawin et al. [45] presented a new poro-mechanical model of the coupled THM processes in partially saturated porous materials that are exposed to temperatures below the freezing point of the pore water. Their model considered freezing-thawing cycles, and can capture important events, such as the effects of hysteresis, the water supercooling, and the cyclic freezing strains.
In related research works, thermodynamically consistent modeling of the phase-change process in porous materials were developed [30,33]. Considering liquid/gas phase transitions within saturated elasto-plastic porous media, TPM models were later presented in [37,38]. Related to the freezing and thawing in porous media including the accompanied phase transition process and volumetric deformation, several models were proposed [14,15,98,99]. In the aforementioned publications, the phase-change process was modeled by imposing a mass production term in the mass balance relation to describe the mass transfer from one phase to another. This method, however, added some complexity and numerical stability challenges to their models. To overcome such problems, an artificial restriction for the volume fraction of the liquid during freezing has to be provided; otherwise, the volume fraction tends to become negative. In the present work, this challenge is avoided by employing the energetic phasefield method to describe the phase-change. Consequently, we consider that ice and water phases belong to the same pore-fluid and are described using PFM-based unified kinematics, thus a more stable and easier computational treatment is provided.
Utilizing the energy-based PFM [11,16,95,105], explicit tracking of the freezing interface is completely avoided, which makes it more attractive among the other interface capturing methods, e.g., the level-set [89] and the volume-of-fluid [54], and easy to implement in finite-element codes. In the context of PFM, a phase-field variable / is employed to indicate the phases of the material, i.e., / ¼ 0 for the solid and / ¼ 1 for the liquid states. The PFM allows for a smooth variation between these distinct values across the interface, resulting in a finite interface thickness. A certain evolution equation for / can rigorously be derived from thermodynamically consistent theories of continuum phase transitions (see, e.g., [115]). The development and applications of phase-field models in the context of continuum models of phase transitions have been presented by many researchers [16,27,43,58,59,78,105,115,119,122]. Away from the phase-change applications, the PFM can also be found in recent contributions within fracture mechanics and crack propagation as, e.g., [6, 39, 51-53, 84, 90, 91, 94]. The PFM is being recently applied to tackle phase-change problems on the macroscale [107,120,121], especially for studying the behavior of phase-change materials (PCMs). However, none of the above mentioned phase-change-PFM-based articles considered the challenging problem of porous structure deformation due to the volume changes of the fluid during solidification or melting. In the current contribution, this challenge will be tackled using an efficient thermodynamically consistent combined TPM-PFM approach.
In a nutshell, the major objective of this work is to gain new insight into the challenging coupled thermo-hydromechanical processes that occur during soil freezing via applying an advanced experimental and numerical study. This is achieved by (1) building a new multi-directional soil freezing testing device that allows the application of different mechanical and thermal boundary conditions, (2) measuring the time-dependent thermo-hydro-mechanical changes in the domain and at the boundaries of the specimens, (3) implementation and validation of a simplified unified water/ice kinematics porous media approach to model the coupled phase-change thermo-hydro-mechanical processes, and (4) extending the modeling framework via a retention-curve-like formulation to consider the microcryo-suction mechanism and, thus, allowing to capture the deformations in the transient and the (quasi-) steady state due to ice lens formation.
To give an overview, Sect. 2 describes the experimental method with specifications of the test setup, the existing boundary conditions and the properties of the test material. The evaluation of the experimental data is explained in Sect. 3, including the measured experimental data as well as the distinction of thermal transient and thermal (quasi-) steady state. In Sect. 4, the results of the two test series with different thermal boundary conditions as well as considering the different freezing directions, are explained. Section 5 provides a description of the fundamental kinematics and basic principles of the TPM and PFM and introduces the governing balance relations. Moreover, an extension of the model for the description of ice lens formation is presented in this section. Section 6 presents the numerical treatment of the governing equations, including the weak formulations, as well as the spatial and the temporal discretization. Verification of the numerical model for the thermal transient state based on experimental data is carried out in Sect. 7.1. Section 7.2 presents the numerical results related to the ice lens formation and a comparison with the experimental results for the thermal steady state.
Finally, a brief summary, conclusions and an outlook of future works are presented in Sect. 8.

Setup of the test device
This project focuses on frost heave analysis due to AGF, which includes different freezing directions. This mimics real AGF scenarios as, e.g., via using horizontal freezing pipes, in which the vertical freezing direction toward the ground surface is decisive in the resulting maximum heave. To investigate the maximum frost heave, an experimental testing device is constructed, which allowed for the examination of the one-dimensional upward freezing (bottom freezing -BF). The results are then compared to one-dimensional downward-freezing tests (top freezing-TF) with a more common setup used for frost heave tests with regard to pavement freezing in cold regions (Fig. 1).
The dimensions of the cylindrical specimen are d ¼ 25 cm in diameter, and h ¼ 20 cm in height. Each specimen is placed in a special test cylinder to achieve a uniaxial deformation during freezing. It consists of eight PVC ring elements with a 2:5 cm thickness. During the freezing test, the ring elements can be pulled apart, so that the soil samples do not undergo any frictional resistance due to the vertical expansion. The top and the bottom of the samples, respectively, are cooled by a freezing element À, which is connected to a cooling circulation thermostat (CC-508, Huber Kaeltemaschinenbau AG, Germany). A water reservoir Þ is connected to a filter plate`at the warm side of the sample, providing a water supply under atmospheric pressure for ice lens formation. The warm end of the soil sample is additionally heated by a heating plate´, that is connected to a warming circulation thermostat (Kiss 104A, Huber Kaeltemaschinenbau AG, Germany) to avoid freezing of the whole sample. All tests are carried out in a climatic chamber with a constant temperature of 15 C and humidity of 40% to exclude influences such as seasonal fluctuations.
For data recording, temperature sensors (MWT1/10, TMH GmbH, Germany)ˆare distributed over the sample height. In order to calculate the temperature profile as a function of time, the vertical displacements of the temperature sensors during frost heave are measured using displacement transducers (LRW2-C-50, WayCon Positionsmesstechnik GmbH, Germany). The uniaxial deformations at the top of the specimen are also recorded by displacement transducers (LRW2-C-50, WayCon Positionsmesstechnik GmbH, Germany)˜. To determine the water inflow into the sample due to suction, the mass of the water supply is recorded by a load cell (S2, HBM, Austria) Þ.

Boundary conditions
In the following, the thermal, the hydraulic, and the mechanical boundary conditions in connection with the test device setup are defined, which are also going to be considered in the related numerical models. Figure 2 shows the boundary conditions for the bottom as well as for the top freezing setups. The freezing element cools down the surface of the soil sample and leads to axial soil deformations. The PVC cylinder avoids a radial expansion of the freezing sample. After reaching the set temperature, the freezing element maintains a constant temperature. The radial heat transfer can be neglected, which is ensured by an insulating PVC cylinder and an additional extruded polystyrene resistance foam covering the test cylinder. The environmental conditions are controlled by the climatic chamber. The ambient temperature for the test series presented here is h a ¼ 15 C, which is also the initial temperature h i of the soil sample. During the entire experiment, water can be sucked from an external supply, which is located at the warm side of the sample. For the presented study, the overburden load results only from construction elements (i.e., the porous heating plate for BF and the freezing plate for TF) and an additional static 10 kg load to ensure a frictional connection between the plates and the sample. This results in a total load of p ¼ 3:6 kPa.

Test material and procedure
The experiments are carried out by considering a lean clay (CL, classified by the Unified Soil Classification System), see, Fig. 3 for the particle size distribution, with a liquid limit of 32:9% and a plasticity index of 14:5%.
The material is produced artificially to ensure an equal composition for all tests and, thus, to improve the comparability of the results. The material is prepared with a water content of 20% to ensure a high saturation up to 95% as well as an appropriate consistency. For saturation, the dry soil is mixed with the required amount of water and stored for 24 hours to ensure that water can penetrate the pores. Afterward, the soil is installed in the test cylinder. The water content present is checked for each soil by a drying procedure.
The essential material properties are summarized in Table 1. These include the mechanical properties, which are determined via tests carried out at the Chair of Geotechnical Engineering, RWTH Aachen University, and the thermal properties, which are determined via tests at the Institute of Applied Geophysics and Geothermal Energy, RWTH Aachen University.
The prepared soil is stored for approx. 24 hours in the climatic chamber at 15 C to ensure a uniform material temperature at the beginning of each test. After acclimatization, the freezing process is turned on. At the same time, the tempering process also starts on the warm side of the sample. During the freezing process, liquid water through a reservoir is provided at the permeable warm boundary, so that the suction from external source is not barriered. The duration of the freezing process is about 10 days per experiment.

Test series
Two test series are carried out focusing on the effects of freezing directions (BF and TF) and the temperature boundary conditions.
In test series 1, the frost penetration is not regulated by the warm side temperature but only by the cold side Young's modulus E S Ã (kPa) 670 Shrinkage limit SL (%) 16.8 Plastic limit PL (%) 18.4 Liquid limit LL (%) 32.9 Plasticity index PI (%) 14.5 Thermal conductivity H SR (W=ðm Á KÞ) 3 Specific heat capacity C S p ( J=ðkg Á KÞ) 776 temperature. With temperature-induced freezing at the bottom or top boundary, the applied temperature varies between approx. À5 C to approx. À15 C. On the warm side, there is only a minimal heat penetration to prevent freezing of the whole sample. The temperature at the warm boundary changes between approx. 5 C to approx. 10 C depending on the set cold side temperature. In this series, the final frost depth at thermal (quasi-) steady state changes with the applied temperature at the freezing boundary. Thus, the temperature gradients over the sample height at thermal (quasi-) steady state are constants (grad h ¼ const:).
To investigate the influence of different temperature gradients, test series 2 includes experiments carried out with adjusted warm side temperatures to reach the same frost penetration depth at different freezing temperatures (grad h 6 ¼ const:). The set temperature gradient in the frozen area grad h f varies between 0:6 C= cm and 1:6 C= cm.

Extracted parameters
The results of each frost heave test (BF and TF) include the measured • total frost heave, • temperatures at six measuring points (incl. motion of sensors), • sucked water mass during ice lens formation.
The total frost heave is measured at the top of the sample to record uniaxial deformations due to pore water expansion and ice lens formation. To determine the position of the frost line during the freezing test, the temperature along the sample height is measured. Figure 4 shows the time-dependent measured temperature along the soil sample and the resulting calculated temperature surface. This spacetime temperature contour plot is used to determine the frost line at 0 C (see Sect. 3.2 for more details). The recorded sucked water from the reservoir provides information about the reservoir water mass used for ice lens formation and, therefore, the influence of gravitational forces (only at BF), which favor the inflow of water.
After each freezing test, the water content is determined in the longitudinal direction of the soil sample to examine the water distribution after freezing. For this purpose, the sample is divided into approx. 2:5 cm thick layers.

Definition of thermal states
During artificial ground freezing, two thermal states result from different operating phases. Initially, the frost body freezes to the required static thickness with a continuously changing frost line position (''moving ice/water boundary''). The penetration velocity depends on the temperature gradient. This state is referred to as the thermal transient state with pore water expansion and less ice lens formation due to freezing pore channels and therefore a restricted water supply. In the subsequent thermal steady state, the thickness of the frozen body is maintained. Since a last ice lens is still formed in this state, it is also called thermal (quasi-) steady state. With unrestricted water supply, the last ice lens causes a very large heave over time. Therefore, this state is especially decisive for long-term projects. In order to investigate the ice lens formation at the thermal (quasi-) steady state, the beginning of this state must be defined in the experiments.
The definition of the thermal (quasi-) steady state, when no further frost penetration occurs, is based on the course of the 0 C-isotherm (Fig. 5, left), which is determined by the temperature space-time distribution in Fig. 4 defining the 0 C over time. The time derivative of the 0 C-isotherm ( Fig. 5, right) gives the frost penetration velocity v fp , which approaches zero at (quasi-) steady state conditions. The beginning of the thermal (quasi-) steady state is defined by a critical frost penetration velocity v fp;crit less than 0:01 cm= h. From this point onward, the water flow to the ice lens is not prevented by freezing pore spaces due to a moving frost line and the frost heave significantly exceeds the heave at the transient state. Figure 6 shows the results of frost heave H at thermal transient and (quasi-) steady state. In the transient state, the pore water heave H PW is determined by the one-dimensional pore water expansion of 9% due to the frost penetration x(t), under consideration of the of the volumetric porosity n F and the saturation degree S r . In the following equation, it is assumed that the water freezes at 0 C and does not expand into air filled voids. The difference between total heave and pore water expansion is the calculated ice lens formation at the transient state. Compared to the ice lens formed at the (quasi-) steady state, the heave due to ice lens formation in the transient state is minor. This is particularly the case for long periods under (quasi-) steady state conditions.

Water content after freezing
The measured water mass sucked in from the reservoir does not reflect the total water mass used for ice lens formation. A comparison of the sucked in reservoir water with the total water mass, calculated back from the total heave due to ice lens formation, shows a difference, which reflects the sucked in pore water. The amount of sucked in reservoir and pore water depends on the flow path length between the water reservoir and the ice lens. A longer flow path leads to a lower pressure gradient over the unfrozen part and thus a lower water accumulation of reservoir water. With a decreasing water supply from the reservoir, the amount of pore water of the unfrozen soil, which contributes to the ice lens formation, increases. Figure 7 (left) shows an example of the ratio of sucked water from the reservoir and pore water to the total water mass used to form ice lenses. At the beginning of ice lens formation, most of the water that forms the ice lens is pore water. This is due to the low permeability of the soil and the long flow path from the reservoir to the ice lens, which restricts the flow rate of the reservoir water. The unfrozen soil shrinks as water drains off, resulting in longitudinal cracks along the sample, whereby the permeability of the unfrozen soil increases and additional reservoir water flows through the unfrozen part to the arising ice lens. This phenomenon of crack formation was also observed earlier by other scientists (e.g., [8,25,109]) and is attributed to particularly impermeable soils. Figure 7 (right) shows the water content along the soil sample (about 2:5 cm layers) after freezing. Compared to the initial water content (20%) the unfrozen part shows a lower water content, especially near the last ice lens. In the frozen area it is higher compared to the initial state caused by slight ice lens formation during the thermal transient state.

Test series 1
In test series 1, only the cold side temperature is varied between À5 C and À15 C to reach different frost penetrations. The ice lens heave at the thermal (quasi-) steady state H IL depends on the length of the unfrozen part or flow path d uf between the ice lens and the sample surface with water supply (see, Fig. 8). With a lower temperature at the freezing boundary h f , greater frost penetration and, therefore, a decrease in the flow path occurs. The shorter flow path enhances the water supply for ice lens formation due to an increased water flow velocity from the water reservoir to the ice lens. Thus, more frost heave due to ice lens formation occurs.
The frost heave H results in Figs. 9 and 10 show differences between BF and TF at thermal transient and (quasi-) steady states. Considering the thermal transient state (Fig. 9), the differences between the TF and the BF result from the different frost depth due to the freezing temperature. A comparison between the two different freezing directions shows the influence of overburden soil pressure. The higher the overweight during BF, the lower is the heave gradient at the beginning, which results in a lower total heave at the end of the thermal transient state.
A comparison of the ice lens heave H IL with different freezing directions shows that the BF exceeds the TF with smaller flow paths d uf (Fig. 10, left). A similar result is shown for the sucked water mass m w (Fig. 10, right). The increasing heave difference between BF and TF is also attributed to better water accumulation in the unfrozen part according to the formation of longitudinal cracks during freezing (Fig. 11). At BF, water accumulation increases due to gravity, especially when the frost line is near the sample surface, causing the cracks to extend to the warm sample surface where the water from the reservoir can flow into the sample. The better water supply results in larger heave due to ice lens formation. If the frost line is not close to the warm side, the cracks do not reach the water source (warm side surface). In this case, the heave difference between BF and TF is primarily due to the overburden load. Depending on the position of the frost line, the overburden load also changes slightly. This effect cannot be excluded, thus further research needs to be carried out to investigate the influence of overburden pressure.

Test series 2
In test series 2, the freezing and heating temperatures were set to achieve the frost line at thermal (quasi-) steady state approximately in the center of the sample. The set temperature gradient in the frozen area grad h f varies between 0:6 C= cm and 1:6 C= cm.
During thermal transient state, the ice lens formation depends on the temperature gradient, because a slow frost penetration rate improves water accumulation. For the researched temperature gradients the heave differences were negligible in comparison to the final ice lens that is formed at (quasi-) steady state conditions. Figure 12 (left) shows the frost heave 150 hours after the beginning of the thermal (quasi-) steady state for different temperature gradients in the frozen area grad h f . The  Fig. 12 (right) and are controlled to be similar by the temperature boundary conditions. Despite different temperature gradients in the frozen area grad h f , the same frost heaves H IL occur. It can be concluded that the temperature gradient in the frozen area grad h f have no decisive impact on the formation of ice lenses at thermal (quasi-) steady state in the investigated range.

Multiphase continuum porous media frost action model
We start in this section with a mathematical model that describes the main physical mechanisms occurring in the thermal transient state during the saturated-soil freezing. In particular, there are the temperature-field distribution, the phase-field evolution, the volume expansion, and the accompanying frost heave. In an additional section, we will briefly demonstrate and include the micro-cryo-suction effect, which is the main driving force for the water suction toward the frozen area and the ice lens formation, which occur mainly during the thermal (quasi-) steady state.  In the current study, the saturated soil is modeled as a twophase material consisting of immiscible solid (a ¼ S) and fluid (a ¼ F) constituents. In the macroscopic TPM-based description of the microscopically heterogeneous twophase material, representative volume elements (RVEs) are chosen and homogenized so that each macroscopic portion of the material u consists of superimposed and interacting constituents u a . While the solid phase u S undergoes linear elastic deformations, the fluid phase u F undergoes a phase transition, where it could be found either in a liquid state (u L ), an ice state (u I ), or a liquid-ice mixed state representing the interface zone separating (u L ) and (u I ). Therefore, the following relations are assigned for any macroscopic spatial point xðtÞ at any time t 2 ½t; T: The volume fraction n a , which describes the local composition of the single constituent, is defined as the ratio of the volume element dv a of the constituent u a with respect to the volume element dv of the overall aggregate u with 0\n a :¼ dv a =dv\1. Analogously, volume fractions representing the liquid (b ¼ L) and the frozen (b ¼ I) states of the pore-fluid, can also be defined as n b :¼ dv b =dv with 0 n b n F . Having both the solidity (n S ) and the porosity (n F ), the saturation constraint for the case of a fully saturated medium can be given as X a n a ¼ n S þ n F ¼ 1 with n F ¼ n L : liquid state; n I : ice state; n L þ n I : interfacial state: Moreover, for each constituent, one defines a material (intrinsic) density q aR :¼ dm a =dv a and a partial density q a : ¼ dm a =dv with dm a being the local mass element of u a . Both densities are related to each other via q a ¼ n a q aR .

Phase-field modeling of pore-fluid freezing
The phase-change process during the pore-fluid freezing is modeled based on the PFM. The PFM is an interfacecapturing approach that treats the occurring mushy zone between the liquid and ice states as a diffusive interface, where an auxiliary phase-field variable / F ðx; tÞ is used to indicate the phases. In particular, this can be expressed as Using / F , a constitutive interpolation polynomial Pð/ F Þ representing an averaging function that interpolates between the bulk properties and free energies of the liquid and ice phases is employed to allow for a phase-field-dependent smooth variation of the different thermo-mechanical properties between the liquid and frozen phases, e.g., [95]. It can be expressed as For instance, the fluid material density is expressed as The phase-field-dependent variation of other kinematic or material properties across the diffusive interface will briefly be discussed in the subsequent sections. More details can also be found in [108].

Kinematics of multi-phase materials
In a fluid-solid aggregate, the deformation of the solid matrix is given within a Lagrangian description using the solid displacement vector u S and velocity v S , whereas the fluid constituents are described in a modified Eulerian setting by use of the seepage velocities w F , i.e., where v F is the fluid velocity. Additionally, the material time derivative of an arbitrary vector quantity ðÞ following the motion of u a is expressed with respect to the motion of the solid constituent as In the present contribution, the model development is restricted to linear isotropic thermo-elasticity. Thus, the linearized infinitesimal solid strain tensor e can be expressed as

Governing balance relations
The present work proceeds from the following assumptions and simplifications: • There is no mass production or mass exchange between the solid and the fluid constituents, i.e., in Eq. (9) q a ¼ 0. • Non-dissipative phase transition, i.e., the capillary hysteresis during freeze-thaw cycles is ignored. • The freezing-point depression of the pore water is neglected.
The discussed mathematical description for the freezing process in a fluid-saturated porous medium results in the following coupled partial differential equations (PDE) with primary variables {p F , w F , u F , u S , h, / F }, where p F is the pore-fluid (ice or water) pressure. More details about the derivation of the balance and constitutive relations can be found in [108].

Mass balance of u
This mass balance could be reformulated based on Eq. (7), the derivative of the saturation condition (3), and the volume balance of the solid, ðn S Þ 0 Assuming a barotropic pore fluid as discussed in [17,53,106], a relation among the material density q FR , the pore-fluid pressure p F and the fluid compressibility j F : Inserting this relationship into Eq. (10) yields the governing equation of fluid mass balance as It is worth mentioning here that the volume balance of the materially incompressible solid phase under the smallstrains assumption yields the deformation-dependent solid volume fraction as n S ¼ n S 0S ð1 À div u S Þ with n S 0S representing the initial value of n S .

Momentum balance of the overall aggregate u
The momentum balance of the overall aggregate results from the sum of the momentum balances of the individual constituents u a yielding div Herein, T a E with a 2 fS; Fg represent the effective stress tensors and b denotes the mass-specific body force. The effective stress tensors of the solid and fluid constituents can be expressed according to [108] as Therein, e D a represents the deviatoric strain tensor, j a and l a are the compression and shear moduli, respectively, and m a h ¼ À3 j a a a h is the stress-temperature modulus with a a h being the coefficient of thermal expansion. a I v denotes the volumetric expansion coefficient due to phase change of u F and Dh a 0 represents the temperature variation compared to the reference temperature h a 0 . The subscript {0 a} refers to the quantity in the reference configuration. The quantities defined in Eq. (14) are obtained based on the restrictions following the evaluation of the entropy inequality and the definition of the free energies w a where T a E :¼ q a o w a =o e a . For the detailed explanation and derivations the interested reader could refer to [46,48,108]. Noting here that the fluid displacement, needed in the computation of e F in the case of ice, is computed based on the fluid velocity through integrating the equation of the material time derivative together with the governing equations at each time step, i.e.,

Momentum balance of u F
Neglecting the dynamic effects and assuming isotropic permeability throughout the freezing and deformation processes, the momentum balance for u F can expressed as Herein, l FR [ 0 is the fluid dynamic viscosity and K S represents the isotropic intrinsic permeability. The deformation dependence of K S is being considered following [41] via where K S 0S is the intrinsic permeability of the undeformed porous solid and j is a material parameter that governs the dependency of the permeability on the solid deformation. For more information regarding the permeability with anisotropic properties, the reader is referred to [24]. The term K S =l FR in Eq. (16) represents the specific permeability K F , which decreases to a very low value upon ice formation as the fluid dynamic viscosity is given as l FR ¼ Pð/ F Þl LR þ 1 À Pð/ F Þ À Á l IR with l IR ) l LR . This resembles the decrease in the hydraulic conductivity or Darcy flow coefficient as ice forms and also ensures that the fluid velocity v F reduces to values almost equal to the solid velocities v S when the fluid is in the frozen state.

Energy balance of the overall aggregate u
The local energy balance for the homogenized medium is reformulated using the Legendre transformation to determine the internal energy via e a ¼ w a þ h g a with g a being the specific entropy. Considering the thermal-equilibrium assumption in this work, the temperatures of the constituents at each material point are approximately the same. Thus, it is suitable to consider the energy balance for the overall aggregate, which is concluded via the sum of energy balances of the individual constituents. Moreover, since the current numerical study is restricted to the small deformation regime, it is convenient to neglect the influences of the deformations on the internal energy. The nonlinear terms were proved in [108] to have a negligible effect based on the numerical simulations. Therefore, the simplified version of the energy balance for the overall aggregate u could be expressed as In this equation, h a ¼ h, cf. Sect. 5.3, C a p represents the specific heat capacity, L F is the latent heat of fusion of ice at the freezing point, H aR represents the scalar-valued heat conductivity coefficient, andp F is the fluid direct momentum production. Related to the PFM, is the gradient energy coefficient that controls the width of the transition zone between the solid and liquid zones, W is an energy hump parameter and Gð/ F Þ ¼ ð/ F Þ 2 ð1 À / F Þ 2 is a double-well energy function. G 0 ð/ F Þ and P 0 ð/ F Þ are the derivatives of Gð/ F Þ and Pð/ F Þ with respect to / F , respectively.

Phase-field evolution equation
The evolution equation for / F is obtained as constitutive equation based on the evaluation of the entropy inequality as discussed in [108]: In this, M [ 0 is the interface mobility parameter. According to [16,116], the parameters , M and W could be obtained by relating them to the material parameters via with c IL , m F k and d being the ice-water surface tension, kinetic coefficient and interface width of the fluid, respectively. The derivation of such parameters is based on solving the one-dimensional phase-field equation for isothermal coexistence of liquid and solid phases at a planar interface of a pure substance, see, [115,116].

Ice-water equilibrium in pores
Following well-established premelting dynamics theories, which describe liquid freezing in the pores of a porous solid material, a thin unfrozen liquid film is assumed to exist at the interface between the frozen liquid and the curved solid grains (pore wall), see, e.g., [97]. The existence of this premelted liquid film can be explained based on the minimization of the surface energies at equilibrium, where the surface energy of ice-grains interface (c IS ) is greater than the sum of the surface energies of the iceliquid (c IL ) and liquid-grains (c LS ) interfaces. This liquid film is characterized by a low pressure, which induces suction of liquid from the surrounding unfrozen ambient leading to additional ice lens formation and frost heave. In this work, to take this suction effect in the soil freezing continuum model into account, a phenomenological formulation similar to the soil/water retention curve (SWRC) in unsaturated soils is implemented. More details about frost heave physics and modeling within a continuum mechanical framework can also be found in [79,87,88,104]. Starting with the mechanical equilibrium at the microscopic pore scale, the relationship between the pressures of the liquid water and the ice crystals can be described by the Young-Laplace's law as where c IL stands for the interface energy between ice and liquid water ( J=m 2 ) and r denotes the mean curvature radius of the interface separating the crystal from the adjacent liquid water. Similarly to Eq. (21), the mechanical equilibrium between vapor (gas) and liquid water of unsaturated porous materials obeys also the Young-Laplace's law as where c GL is the liquid-air interface energy. The similarity between Eqs. (21) and (22) paves the path to implement similar macroscopic formulations that consider the suction effects. In particular, we proceed from the van Genuchten model [114], which is used to express the relationship between the pore capillary pressure and the liquid saturation v L during the drying process of unsaturated porous materials and is expressed as In this, p G is the gas pressure, whereas N Ã and m Ã are adjustable model parameters related to the material pore structure and pore-fluid properties. In analogy, the cryogenic suction due to the ice/water interface tension, i.e., S cryo :¼ p I À p L , can be expressed by the following phenomenological relationship where N and m are retention model parameters, adjusted depending on the solid and fluid phases properties. The function in (24) is introduced as a non-hysteretic phenomenological characteristic function to relate the liquid water saturation v L with the pressure difference between the crystal ice and the liquid water ( S cryo ¼ p I À p L [ 0) [82,124]. In our treatment of unified water/ice kinematics, the phase-field variable / F is utilized as the degree of water saturation and, thus, corresponds to v L . In order to include the suction effect in the macroscopic model presented in Sect. 5.3, one has to distinguish between p L and p I , whereas p I À p L can explicitly be determined from Eq. (24). For this purpose, the definition of the net pore pressure p :¼ p F is exploited, which is equal to the weighted sum of the pore ice pressure p I and the pore water pressure p L [88]. Based on the PFM formulation with v L % / F , this can be described as This reformulation is applied to eliminate the ice phase pressure p I from the primary variables by expressing it as a function of the liquid phase pressure p L and the cryo-suction S cryo . Note that in the current reformulation, p L becomes the primary pressure variable instead of p F . Therefore, to overcome the instability challenges, a quasicompressibility method is used to enforce the stability via adding a stabilization term of the form ''Àa div ðgrad p L Þ'' as mentioned in Sect. 6.

Ice segregation criterion
A stress criterion for the initiation of a new ice lens and, thus, degradation of the solid stiffness, was adopted by many researchers [47,83,88]. The stress-based method considers that in a one-dimensional configuration, the soil stiffness degrades and a new ice lens forms when the pore pressure exceeds the sum of the overburden stress p ob and the separation strength of the soil p sep . The frost heave then results from the moisture accumulation and the volumetric expansion due to ice lens formation. In this paper, an alternative criterion approach presented by [123] and [77] is adopted, where a separating porosity or void ratio, as a counterpart to the critical separation stress, determines the formation of ice lenses. In particular, when the void ratio e :¼ n F =n S is greater than or equal to the separation void ratio e sep , i.e., e ! e sep , ice crystals interconnect and ice lenses begin to form. The separation void ratio can be expressed in one-dimensional space according to [77,123] as Herein, c :¼ g n S q SR þ n W q WR þ n I q IR ð Þ is the unit weight of the soil, and the subscript 0 refers to the initial state. To represent the damage to the solid phase once ice segrega-tion takes place, the stiffness of soil is allowed to degrade to a very low value by multiplying the effective stress by a degradation function ! d . This method is chosen in analogy to degradation of the stiffness of the solid matrix during the crack growth in porous media [52,53]. Moreover, this method agrees with described methods in the literature as in [103,111], which utilize other approaches for ice lens formation. In particular, a phenomenological parameter is presented in the sense of an indicator function for ice lens formation. It can be expressed for a smooth transition as where j D is the residual stiffness after separation, which is set 10 À2 in the simulations so that no numerical stability problems occur. Moreover, e max is a maximum value of the void ratio, after which the stiffness of the solid phase is considered to be completely lost. In the ice lenses simulation, the degradation function in Eq. (27) is multiplied by the solid effective stress T S E defined in Eq. (14) 1 . Therefore, the momentum balance of the overall aggregate u in Eq. (13) can be written as which is applied in the ice-lens modeling.

Numerical Implementation
In the numerical solution of the coupled partial-differential equations (PDEs), the finite element method (FEM) within a Bubnov-Galerkin approach is applied. As a representative case, the original model presented in Sect. 5.3 will be considered, whereas the modified model in Sect. 5.4 can be treated in a similar fashion. In this, the PDEs with primary variables nðx; tÞ :¼ ½p F w F u F u S h / F T are multiplied by the corresponding independent test functions dnðxÞ :¼ ½dp F dw F du F du S dh d/ F and integrated over the spatial domain X T to derive the weak formulation. Thereafter, the product rule and the Gaussian integral theorem are utilized to extract the boundary terms C n ¼ o X n , where C n is split into Dirichlet (essential) and Neumann (natural) boundaries, see, e.g., [49,50,73,126] for analogous treatment. For the case of a nearly incompressible pore fluid, i.e., a very high bulk modulus j F , the term corresponding to ðp F Þ 0 S in the fluid mass balance Eq. (12) tends to vanish, and Eq. (12) becomes independent of the pressure variable. Thus, numerical instabilities appear in the monolithic solution due to the violation of the inf-sup condition. To overcome this instability challenge, a quasi-compressibility method is used to enforce the stability via adding a stabilization term of the form ''Àa div ðgrad p F Þ'' to Eq. (12) with a as a small stabilization parameter. Following this, the weak form of the governing equations read: where w ¼ q F w F Á n is the mass efflux of the fluid draining through the Neumann boundary C w with outward-oriented unit surface normal n.
Herein, t F ¼ T F E À n F p F I À Á n represents the external fluid load vectors acting on the Neumann boundary C t F . The extra equation induced by the use of the relation where t ¼ r n represents the external solid load vector acting on Neumann boundary C t .
• Energy balance of the mixture u Z where H SF represents H S þ H F and q SF ¼ H SF grad h ð ÞÁn is the heat flux acting on Neumann boundary C q SF .
In this, the term 2 grad / F À Á Á n vanishes on the Neumann boundary C / F . The TPM-PFM coupled problem is solved using the FEM package FlexPDE (professional, ver. 7.15.), where the time discretization of the equations is carried using the secondorder implicit backward difference formula. For the spatial discretization, the Galerkian method with quadratic shape functions is applied. Having the time-dependent nodal values vector of the primary variables nðtÞ :¼ ½p F w F u F u S h / F T , the coupled PDE system after spatial discretization can be written in an abstract form as where G p F , G w F , G u F , G u S , G h and G / F represent the spacediscrete fluid mass balance, fluid momentum balance, fluid displacement, mixture momentum balance, mixture energy balance, and the phase-field evolution equation, respectively.
Related to the time integration, the volumetrically coupled thermo-mechanical problem is solved using a monolithic implicit scheme with the quasi-compressibility stabilization scheme for the case of a low-compressible fluid. It is worth mentioning that the computational efficiency gained by the staggered strategy in thermo-elasticity is not guaranteed for the case where the differential equations are strongly coupled (many nonlinear dependencies) [42]. In the current implementation, time-stepping adaptivity is applied to keep the solution within the specified limits of accuracy. Regarding the system of nonlinear differential equations, a modified Newton-Raphson iteration procedure is utilized by FlexPDE to attain the numerical solution, where the Jacobian matrix of the current system equations is computed and Newton iterations are performed at each time step until the convergence is reached. Useful and detailed information to the solution of coupled problems can be found in [1,2,72,73,85,107] and the included references.
In dealing with the phase-change process, the mesh size should be small enough so the solid-liquid interface is efficiently covered. For an accurate representation of the physical processes at the moving, diffusive front, it is very important to have small values of the parameter since it plays a central role in the sense that the diffusive interface converges to a sharp one for ! 0. However, from the computational cost viewpoint, it is desirable for to be as large as possible so that a solution of the phase-field equation can be obtained for a practical computational effort. In this sense, interfaces are taken artificially wide in most of the phase-field simulations to decrease the computational costs. In several studies it was shown that when the PFM is applied to model the phase-change process on the macroscale, without the interest in capturing the microstructural evolution and dendrite growth, the interface width could be increased by orders of magnitude without significantly affecting the result [22,108,120]. This serves significantly in decreasing the computational costs by eliminating the need for a very fine mesh size.  Table 1, while those of ice and water are adopted from [13]. Other material parameters related to the PFM are chosen as follows: c IL ¼ 0:065 J= m 2 , m F k ¼ 0:001 m=s Á K, d ¼ 0:5 mm, and ¼ 0:8 ð J= mÞ 1 2 . Additionally, triangular grids with mesh size of ca. 4 mm was chosen for all the simulations. The initial conditions for both BF and TF are the same, where water is found in a liquid state at an ambient temperature of h i ¼ 15 C. The drained boundary conditions are prescribed by imposing zero pore pressure p F ¼ 0 Pa at the heating boundary of each case to allow for water supply, while impermeable boundary conditions are set at the other boundaries. Moreover, a total mechanical load of p ob ¼ 3:6 kPa corresponding to the overburden pressure is applied to the top boundary, see, Fig. 2 for illustration. In both BF and TF, the sample is being cooled along the cold boundary, while the warm boundary is slightly heated to prevent freezing of the complete sample as discussed in Sect. 2. The cooling and heating temperature loads for the BF and TF tests, which are also applied in the numerical models, are illustrated in Fig. 13.
In discussing the numerical results, the relevant phasefield and temperature contour plots for the BF and TF tests after 20 and 25 hours are shown in Figs. 14 and 15, respectively.
The different positions of the temperature sensors (six sensors) in the experiments are illustrated in Fig. 16. A comparison between the experimental and numerical results of the temperature-time history at different sensor positions is given in Fig. 17. In both BF and TF cases, the temperature starts decreasing due to cooling until reaching 0 C, where the water starts freezing. The temperature continues decreasing with time while the thickness of the ice layer grows continuously as shown in Fig. 18. Figure 19 shows the frost heave evolution extracted from the FE simulations and compared to the measured data. These simulations are performed proceeding from elastic solid skeleton deformations, whereas the heave evolution is measured at the top surface. The resulting frost heave in the transient phase is mainly due to the volume expansion of the sample due to ice formation. These results, which cover the phase change process and the thermo-mechanical behavior of soil during freezing, show a good agreement between the numerical and the experimental date. However, simplification in the applied model, such as ignoring the plastic deformations of soil, have presumably contributed to the deviation in the displacement field in Fig. 19 between the numerical and experimental results.
To provide a percentage-based deviation measure of the predicted numerical results from the reference experimental data, several methods could be found in the literature. In this, the mean absolute percentage error (MAPE) [34] provides a simple dimensionless measure that could easily be computed and interpreted Here, Y t denotes the measured data (experimental data) at time t, F t is the forecast data (numerical predicted data) and n is the number of different times for which the considered variable is forecast. However, this method has inherent difficulty in handling low and zero values in the denominator. Moreover, having a single case with large error can dominate the overall result of the MAPE. To overcome these challenges, especially related to the 'infinite error', the weighted mean absolute percentage error (WMAPE) [75,86] can alternatively be used. This is equivalent to dividing the sum of absolute differences of the measured and forecast values by the sum of the measured values yielding Another alternative for the MAPE is the symmetric mean absolute percentage Error (SMAPE) [44,100], where it expresses the absolute error as a percentage of the average of the forecast and the measured values. However, this measure is found not to be suitable when interpreting percentage errors of data including mixed positive and negative values (as, e.g., temperatures on a C-scale in our case). Moreover, undefined values could appear when the sum of the measured and forecast values is zero in the denominator. An overview of the different types of measures, including scale-dependent measure and relativebased measures, can be found in [26,44]. According to this, we employ the WMAPE (Eq. (37)) for quantitative comparisons between the experimental and the numerical results related to the different outputs, i.e., the frost heave, the temperature field, and the ice/water interface position.

(a) (b)
Fig. 14 Contour plots of a spatial temperature distribution and b phase-field evolution for the BF case after 25 hours (a) (b) Fig. 15 Contour plots of a temperature distribution and b phase-field evolution for TF case after 20 hours As shown in Table 2, the lowest WMAPE values are recorded for the interface position in each test. In the TF tests, the temperatures show the largest deviations followed by the frost heave, while, for the BF tests the frost heave shows the largest deviations among the other outputs. This could be attributed, as discussed previously, to neglecting the plastic deformations of the solid phase in the freezing process. In general, these deviation errors remain in acceptable range, especially that the previously discussed figures showed almost similar trends between the simulation and experimental curves for the different studied cases.   Based on these results and comparisons between the numerical and the experimental data for the different studied cases, we conclude that the proposed TPM-PFM approach is powerful in the simulation of phase-change processes and in the description of the non-isothermal behavior during freezing of water-saturated soils.

Numerical results related to the ice lens formation
For the current numerical study, the suction-related parameters N and m are chosen as 0:1 MPa and 0.5, respectively, which gives a good fit with the experimental data. Here, we consider BF test 3 as a representative case. The boundary and initial conditions are the same as described in the previous sections. We assume a separation strength of the soil p sep of 200 kPa, which corresponds according to Eq. (26) to e sep % 1:1. The temperature loads for this representative case are shown in Fig. 20.

Frost heave
A comparison of the frost heave resulting from the simulations with that of the experimental data is presented in Fig. 21. During the transient freezing process and frozen fringe formation, the heave occurs due to the volumetric expansion of the frozen water (ice-water phase change). As the temperature gradient becomes steady, the rate of frozen front penetration decreases while the frost heave continues its gradual increase. In this phase, the pore water migrates from the unfrozen zone and accumulates at the freezing front, which results in the ice lens formation. This indicates  that the heave due to ice lens formation in the thermal (quasi-) steady state is of great importance. The results also show the capability of the proposed numerical model in predicting the total heave due to pore-water freezing.
Regarding the percentage error related to numerical and experimental results of the frost heave, the WMAPE results in a value of 15:57%. Referring to Fig. 11 and the obtained results in the related section, at the early stages of the experiments (during the thermal transient state), the vertical cracks are insignificant, and the frost heave is mainly due to the ice expansion upon phase-change. In the later steady-state phase, the frost heave is mainly due to the sucked water that comes from the unfrozen surrounding and the external supply through the surrounding, which accumulates and freezes at the freezing front. The vertical cracks grow significantly at the late stages of the steady-state phase, where they play a bigger role in altering the sample's permeability and increasing the frost heave. Figure 21 shows that the frost heave in the numerical result at the late stage of the steady state is less compared to the experimental result. This could be due to the vertical cracks that formed at that stage and the relevant vertical permeability enhancement, which has not been taken into consideration in the model.
Ice lens contour plot at two different times of the freezing and scaled deformations by factor 7. The red color refers to a complete loss of solid stiffness (separation) and a continuous ice lens formation

Ice lens distribution
The predicted ice lens, which fulfills the condition e ! e sep , is illustrated in Fig. 22, where a scaling factor of 7 is used. As long as the time 40 h, the height of the soil column is increasing, but no visible ice lenses are formed, i.e., the void ratio in the domain is less than e sep . As the rate of the frozen front penetration decreases in the thermal (quasi-) steady state, water in the unfrozen zone has enough time to flow to the frozen zone (suction effect) and to form a localized ice lens. In other words, as time goes by, the ice lens becomes larger and forms at greater spacing progressively (e ) e sep ) causing an increase in the frost heave. The numerical analysis agrees with the experimental results regarding the position of the final ice lens. This shows that the proposed model could serve well in the description of the ice lens formation.

Conclusions and future work
Ground frost and the resulting heave-induced damages are mainly known from road construction in cold regions. However, they can also occur during artificial ground freezing. The boundary conditions resulting from artificial ground freezing differ in some aspects from the boundary conditions in road construction, for example, concerning the decisive freezing direction, which runs toward the ground surface during ground freezing. This and other influences on frost heave during artificial ground freezing are investigated in this research work using a newly designed experimental device.
In artificial ground freezing, the frost body first freezes to the static required thickness and is maintained on this thickness thereafter. Due to the large temperature gradients at the beginning of freezing (thermal transient state), less ice lens formation occurs. In this, the experimental results showed that the frost heave due to the formation of the last ice lens in the thermal (quasi-) steady state significantly exceeds the frost heave at the thermal transient state. The  (quasi-) steady state is therefore decisive in the process of artificial ground freezing.
The outlined results show that the unfrozen area, which provides a flow path to the ice lens, has a major influence on ice lens growth. A shorter flow path leads to an increase in the pressure gradient and, therefore, an increase in the pore-water velocity toward the ice lens. In contrast, no significant influence of different temperature gradients on ice lens formation in the (quasi-) steady state could be observed. A comparison of the freezing directions BF and TF showed that especially with short flow paths, the frost heave of the BF experiments exceeds those of the TF. This could be attributed to the effect of gravity and the longitudinal crack formation during freezing. In particular, the hydraulic conductivity increases due to the initiation of the cracks, so that an improved water supply is provided.
The results presented here are part of a comprehensive study to determine frost heave for the boundary conditions of artificial ground freezing. Further experimental investigations concerning the overburden load and the influence of water pressure are pending.
In the numerical simulation of soil freezing and the related instances, a thermo-hydro-mechanical model that augments the macroscopic theory of porous media (TPM) with the phase-field modeling (PFM) for simulating the freezing process in a fluid-saturated porous medium is presented. The water-saturated soil is described within a continuum mechanical framework as a biphase porous material, consisting of an elastic solid (soil skeleton), and pore-fluid (pore-water/ice). In this context, thermodynamically consistent constitutive relationships and phase-field evolution obtained based on satisfying the entropy inequality are utilized. A special feature of the current work lies in presenting a unified phase-field kinematic framework for treating the ice and water phases, considering that both phrases refer to the same pore-fluid. Moreover, the PFM ensures a smooth transition for the different thermo-mechanical properties during the phasechange process.
The mathematical description of the coupled porousmedia phase-field problem results in six partial-differential equations. The numerical solution of these equations is carried out using a monolithic scheme within the finite element framework. The numerical solutions of the proposed TPM-PFM approach together with the qualitative comparisons with experimental results showed the capability of the proposed model to efficiently capture the coupled processes of temperature field development in space and time, the phase-change process, and the volumetric expansion due to ice formation during freezing, and the formation of ice lenses together with solid stiffness degradation.
To this end, the presented numerical model in this research work can serve as a base for future studies and real applications in the field of soil freezing. However, improvement areas are to be addressed in the future work of this research project. Future works will focus on providing a more realistic description of the soil freezing phenomenon by extending the model toward unsaturated porous media. Additionally, the initiation and formation of ice lens will be modeled in future works based on the phase-field diffusive fracture approach. This would require extending the energy formulation through providing a thermodynamically consistent derivation of the fracturerelated PF-equation while considering the associated THM changes in the domain.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this Fig. 22 Ice lens contour plot at two different times of the freezing and scaled deformations byfactor 7. The red color refers to a complete loss of solid stiffness (separation) and a continuousice lens formation article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.