A Multi‐Phase Heat Transfer Model for Water Infiltration Into Frozen Soil

In many regions, water infiltration into frozen soil is a critical process. Frozen soil is known to have a substantially reduced infiltration capacity, resulting in surface water runoff. This surface runoff is often associated with erosive behavior posing a potential natural hazard as it facilitates snow avalanches and debris flow, especially in springtime when rainfall and snowmelt coincide. As infiltration capacity depends on the ice content within the porous soil, thermal and hydraulic processes are strongly coupled. The involved phases during water infiltration into frozen soil are initially not in a local thermal equilibrium as the infiltrated water is warmer than the frozen soil. The duration of the temperature differences is yet to be determined. To adequately describe this thermal state, a local thermal non‐equilibrium (LTNE) model needs to be applied, which accounts for separate phase temperatures and describes the heat transfer between the phases explicitly. While LTNE models are common in more simple setups of saturated porous media flow, in this work, a multi‐phase LTNE model for water infiltration into a partly saturated, frozen soil is presented. All required parameters, such as heat transfer area and heat transfer coefficient are discussed and derived based on a capillary tube model. The theoretical model is implemented into a numerical model and compared to experimental data available from literature. Substantial differences of up to 0.5° C between phase temperatures during freezing and melting are observed after hours of infiltration and are especially pronounced around the freezing front.


Introduction
Water infiltration into frozen soil is a crucial process in many parts of the world as more than 50% of the exposed land in the Northern Hemisphere are experiencing seasonally frozen soils (T. Zhang et al., 2003). Frozen soil has a significantly reduced infiltration capacity depending on the ice content compared to unfrozen soil and therefore infiltration and surface runoff substantially depend on the ice content in the soil (Mohammed et al., 2018). If the liquid water cannot infiltrate into the soil, erosion and substantial surface water runoff might take place causing landslides and debris flow (Seyfried & Murdock, 1997). Also, the contribution of snowmelt to the groundwater resources can be an important factor in regions with otherwise little soil water storage (Hood & Hayashi, 2015). With the snow melt, contaminants originally stored in the snow are transported by the meltwater and infiltrate into groundwater resources (e.g., Feng et al., 2001;Jones et al., 1989). If water infiltrates into the soil and then freezes, frost heave might occur and damage roads, train tracks, pipelines and other critical infrastructure. In general, the vadose zone experiences complex interactions between soil, water, and air with chemical and biological activity within the vadose zone strongly effected by temperature and the availability of liquid water (Hayashi, 2013). Therefore, understanding the flow behavior under freezing conditions and the thermal interaction between soil and infiltrating water is crucial for hazard prevention and groundwater management. The thermal state of the involved phases is key for the hydraulic processes.
Water infiltration into frozen soil is a multi-phase scenario with an -at least initially -local thermal non-equilibrium (LTNE) between the involved phases because the infiltrating water has a temperature above the freezing point, while the soil and its pore filling are frozen. LTNE situations are well-known for a saturated porous matrix (Gossler et al., 2020;Hamidi et al., 2019) and the LTNE concept has been recently extended to partially saturated soils (Heinze & Blöcher, 2019). In the context of geothermal energy exploration, the LTNE concept has gained more and more attention due to the obvious limitations and the oversimplification of models assuming local thermal equilibrium (LTE) (e.g., Gelet et al., 2013;Shaik et al., 2011;C. Xu et al., 2015). LTE models are based on the volumetric weighting of thermal parameters of the involved Abstract In many regions, water infiltration into frozen soil is a critical process. Frozen soil is known to have a substantially reduced infiltration capacity, resulting in surface water runoff. This surface runoff is often associated with erosive behavior posing a potential natural hazard as it facilitates snow avalanches and debris flow, especially in springtime when rainfall and snowmelt coincide. As infiltration capacity depends on the ice content within the porous soil, thermal and hydraulic processes are strongly coupled. The involved phases during water infiltration into frozen soil are initially not in a local thermal equilibrium as the infiltrated water is warmer than the frozen soil. The duration of the temperature differences is yet to be determined. To adequately describe this thermal state, a local thermal non-equilibrium (LTNE) model needs to be applied, which accounts for separate phase temperatures and describes the heat transfer between the phases explicitly. While LTNE models are common in more simple setups of saturated porous media flow, in this work, a multi-phase LTNE model for water infiltration into a partly saturated, frozen soil is presented. All required parameters, such as heat transfer area and heat transfer coefficient are discussed and derived based on a capillary tube model. The theoretical model is implemented into a numerical model and compared to experimental data available from literature. Substantial differences of up to  0.5 E C between phase temperatures during freezing and melting are observed after hours of infiltration and are especially pronounced around the freezing front.

of 21
phases and assign all involved phases one common mixture temperature. LTE models are straight forward in the mathematical formulation and numerical implementation. On the other hand, LTNE models are based on separate heat equations for each phase but the principle can be extended to other scenarios of separate domains, e.g., separating pore and fracture fluid in a dual domain approach (Heinze & Hamidi, 2017). Heat exchange between the phases is described explicitly using a heat transfer term depending on the heat transfer coefficient, the specific contact area and the temperature difference between the phases (e.g., Nield & Bejan, 2013). The determination of the specific surface area for porous media is based on geometric assumptions of the porous matrix, e.g., with respect to grain size (Nield & Bejan, 2013). The heat transfer coefficient is often set constant based on experimental values, but experimental evidence shows that the heat transfer coefficient depends on flow velocity, among other parameters (e.g., Wakao & Funazkri, 1978;Zhao & Tso, 1993). Very few analytical models for the heat transfer coefficient exist. For fracture flow the velocity dependence of the heat transfer coefficient is derived to be first order (Heinze, 2021), while experimental data suggests dependence with the third root for the Darcy velocity in porous media (Gossler et al., 2020;Roshan et al., 2014;Wakao & Funazkri, 1978).
There is a number of theoretical and numerical models describing freeze and thaw processes in soil. Hansson et al. (2004) present a fully implicit numerical model for coupled water flow and heat transport in variably saturated porous media above and below the freezing/melting temperature. Based on the Richards equation and the assumption that soil freezing has similarities to soil drying (e.g., Spaans & Baker, 1996), a energy-and mass-conservative solution is derived and implemented in Hydrus-1d. The model is compared to experiments of frozen soil columns. Those laboratory experiments are used for comparison with numerical models in various studies (e.g., Dall'Amico et al., 2011;Kelleners, 2013). Dall'Amico et al. (2011) present a physical and numerical model for soil freezing in variably saturated soil using a splitting method, separating equations with constant temperature being only functions of the liquid water potential, and equations with constant water potential and variable temperature. They formulate the energy equation in terms of the internal energy to allow a formulation similar to a diffusion-advection equation like the mass balance equation. Models with various complexity for freezing and thawing in soil, as well as in snow, are presented by Kelleners et al. (2009);Kelleners (2013) and Kelleners et al. (2016). Those models further include root water uptake and water vapor flow in addition to freezing and melting. Similarity of mass and heat transport in snow and frozen soil is addressed in detail by Kelleners et al. (2016). Coupling between snow and soil, also with varying snow density and time varying thermal parameters has also been addressed in Rawlins et al. (2013). In the broader sense, permafrost hydrological models on a catchment scale incorporate vertical and lateral water transfer, e.g., during snowmelt (Y. Zhang et al., 2013). Multi-phase heat transport assuming thermal equilibrium between phases has been studied for a river-plain system over several months (Grenier et al., 2013). Popular major numerical frameworks for modeling freezing and melting in porous media include the USGS SUTRA model (McKenzie et al., 2007), the COUP model (Jansson & Karlberg, 2001) the SHAW model (Flerchinger et al., 1996) as well as FeFlow with the extension piFreeze (Langford et al., 2020;Magnin et al., 2017). A comparison of various numerical programs is given through the InterFrost initiative (Grenier et al., 2018).
There are also a couple of discontinuous models to explicitly describe the interface between frozen and unfrozen soil or water on different scales. A semi-analytical model based on dynamic interfaces between frozen and unfrozen as well as saturated and partly saturated layers has been developed for the application in permafrost degradation (Devoie & Craig, 2020). A discontinuous three-phase heat and mass transfer model is presented in F. Xu et al. (2018) and implemented using the Lattice-Boltzmann method (Y. Zhang et al., 2018). Their model describes the pre-melted film surrounding the soil grains and links film thickness to soil temperature. Interactive forces between soil, water and ice are used to describe the moisture movement (F. Xu et al., 2018). The model successfully reproduces experimental findings by Hansson et al. (2004) and others. In Peng et al. (2016) a thermal non-equilibrium between the ice and the liquid water is introduced. While the Richards equation is used to describe liquid water movement in frozen and unfrozen soil, it is modified with a sink term coupled to the freezing of liquid water. The numerical 1D model shows an even better agreement with the experimental data of Hansson et al. (2004) than many other continuous models (Peng et al., 2016).

of 21
All those continuous as well as discontinuous models are based on the LTE assumption. Therefore, a common temperature is assigned for all involved phases and a potential temperature difference between phases cannot be resolved. For infiltration scenarios, the LTE assumption results in limits of the model design as initial or boundary conditions have to be formulated for a common temperature, e.g., based on volumetric weighting. The duration of a temperature difference between phases during water infiltration into frozen soil has not been quantified so far and the applicability of the LTE assumption therefore remains open. Experimentally, water infiltration into a frozen soil is studied seldom compared to standard freezing or thawing experiments. An experimental investigation of water infiltration into frozen soil is presented by Watanabe and Kugisaki (2017) with special emphasis on the role of macropores in the soil.
In this work, a mathematical framework is presented for heat transfer between multiple phases with dynamic volume fractions. The possible phase change between liquid and frozen water in both ways, freezing of infiltrating water and melting of initially frozen pore water, pose special complexity to a heat transfer model due to the dynamic volume fractions of the involved phases. Additionally, air is replaced by infiltrating water in the described scenario. The different heat transfer terms between the phases and their parameterization are discussed in detail on the basis of a capillary bundle model. The mathematical model is then implemented into a numerical model and compared to experimental data from Watanabe and Kugisaki (2017).
Common models of soil freezing are mostly based on a modified Clausius-Clapeyron equation, that links hydraulic head and temperature. Subsequently the amounts of frozen and liquid water are determined (Kurylyk & Watanabe, 2013). For melting scenarios as considered in this work, thermal equilibrium between liquid water and ice required for the application of the Clausius-Clapeyron equation is not given. Therefore, the phase change needs to be described explicitly. Such an approach was developed by Kelleners et al. (2009) for the melting of snow. Once a phase temperature reaches the temperature of phase transition, any additional change in thermal energy, increase during melting and decrease during freezing, is transferred into the phase change of some amount of the material.
The melting of the existing frozen water can occur when its temperature f E T is equal to the temperature of phase transition ph E T . Any additionally available energy is used for phase transition. The energy required is determined by the latent heat of fusion f E L and given by (Kelleners et al., 2016) with the amount of melted ice , f m E  . Similarly, for freezing the energy of the liquid water is reduced by the latent heat given by (Kelleners et al., 2016) with the amount of melted ice , l f E  . As freezing and melting cannot occur at the same time at the same location, the change in volume fractions can be summarized as done above by the use of one common volume fraction , l fm E  .

Heat Transfer Area Between Phases
The heat transfer area and the heat transfer coefficient are required, besides the phase temperatures, to calculate the transferred heat. To describe the porous matrix, a bundle of capillary tubes is a common conceptional 6 of 21 model. The pore space is split in separate tubes with the shape of cylinders. Different pore radii are represented by a selection of tubes with varying radii of the respective cylinders. The porosity is then given as the ratio of the volume of the tubes t E V ( 3 m E ) and the total volume of the matrix m E V :   V V t m / . For a capillary bundle with E M different pore radii, and E N numbers of pores with the respective radius E R (m), the volume of the tubes is with length E L (m) and tortuosity  E (−) of the tubes. The surface area t E S ( 2 m E ) of all tubes is accordingly Subsequently, the specific surface area can be expressed as For a capillary bundle consisting of a homogeneous pore radius E R , this equation simplifies to It is noteworthy, that in such a formulation the specific surface area increases with increasing porosity because porosity depends on 2 E R . This is contradictory to a typically used formulation of the specific surface area for soils assuming spherical grains of homogeneous size (Dullien, 1992;Nield & Bejan, 2013). In a recent study, a similar relationship as presented here was derived for consolidated rock assuming spherical or elliptical inclusions embedded in a matrix, while the spherical grain approximation has been found valid for rather loose materials such as dune sand (Hussaini & Dvorkin, 2021). Considering fine grained soil, the derived relationship here is considered a suitable and realistic representation. The tube model has multiple benefits regarding its use for multi-phase models and has been used successfully in previous studies, e.g., for describing hydraulic conductivity in frozen soils (Watanabe & Flury, 2008). The tube radii distribution can be derived from the infiltration behavior and the respective parameters as described in Watanabe and Flury (2008).
After defining the heat transfer area of the porous matrix, the surface area of the frozen water can be described using the same approach. The volume fraction of the frozen water f E  is the ratio of the frozen water and the total volume. Similar to the tube model above, the specific surface area of this ice cylinder inside the tubes can be written as assuming that each tube with similar radius is also filled with a similar ice core with radius f E R (m). Again, for a medium with only one pore radius, this simplifies to ior is substantial for the consistency of the model. Further, as the radius of the ice core f E R is linked to the volumetric water content of the frozen phase by  , any change in frozen water content reflects on the heat transfer area. With the same relationships as above, the dependence of ice core radius f E R on the volume fraction of the frozen water can be derived The derivation of this relationship is purely geometric. Any effect of the van der Waals force or any other process inhibiting or restricting ice core grow is neglected (cf. Dash et al., 1995). However, such effects could be included by a representative description of the volume fractions, e.g., limiting the growth of f E  .
Calculating the heat transfer areas for soil * s E A through Equation 17 and for frozen water * f E A through Equation 19, the contact area between the liquid and the gaseous phase still needs to be defined. Considering a scenario in which air is replaced by infiltrating water, the contact area in the described capillary tube model is the cross section of these tubes. Compared to the length of the tubes, the cross section is negligible and it becomes even smaller in the presence of an ice core.

Heat Transfer Coefficients Between Phases
For a liquid or gaseous phase moving through a porous media, the most common expression for the heat transfer coefficient is based on the work of Wakao and Funazkri (1978) (Gossler et al., 2020;Hamidi et al., 2019;Roshan et al., 2014). This empirical expression links the heat transfer coefficient in dependence of the Prandtl and Reynolds numbers. The dimensionless Prandtl number Pr E is the ratio between viscosity and thermal conductivity for liquids and gases (Nield & Bejan, 2013)  The dimensionless Reynolds number Re E (−) is a quantity describing the flow behavior of liquids and gases over a system relevant length scale, such as the particle diameter The heat transfer coefficient E h between a stationary porous matrix and a flowing fluid is then expressed as (Wakao & Funazkri, 1978)   Besides the particle diameter, all quantities in this expression are the properties of the fluid/gas. Most experimental studies, resulting in comparable relationships, have been performed for gases flowing through porous media in the context of mechanical engineering applications (e.g., Singhal, 2017;Zanoni et al., 2017;Zhu et al., 2019). Geoscientific applications and representative materials are studied comparably seldom (e.g., Gossler et al., 2019;Roshan et al., 2014). For aquifers, a recent study proposed a relationship independent of the Prandlt number (Gossler et al., 2020). In constrast with engineering materials, such as between air and metal used in heat exchangers, the ratio of thermal parameters of water and sediment rocks is significantly smaller. Therefore, the influence of those parameters in negligible in geoscientific applications (Gossler et al., 2020).
Considering the multi-phase nature of this work, the more general approach by Wakao and Funazkri (1978) is applied because it enables the calculation of the heat transfer between soil, flowing water and gas with a common equation. Through the inclusion of the Reynolds number, the heat transfer coefficient is also velocity dependent. In the considered scenarios of this work, this is especially important as hydraulic conductivity is assumed to increase rapidly, when ice blocking the cores is melting. As all studies show a negligible dependence of the heat transfer coefficient on the thermal parameters of the solid phase, the same Equation 24 is also used for describing the heat transfer between the liquid and gaseous phase and the ice core. However, it needs to be noted that there have not been done any heat transfer experiments with ice up to this point. The heat transfer coefficient between air and liquid water inside a porous media has also not been studied experimentally so far. For open waters, several studies tried to obtain the heat transfer coefficient for the water-air interface and indicate a significantly lower value than for the air -rock interface (e.g., Kalinowska, 2019;Laguerre et al., 2017; h and fl E h .
Due to the dependence on flow velocity, all heat transfer coefficients can change over time according to the hydraulic dynamics of the system.

Heat Transfer Terms Between Phases
The heat transfer terms between the separate phases are all based on Newton's law of cooling given in Equation 3. In the most general sense, the heat transferred from one phase to another is a source term for one phase and a sink for the other phase. For the interaction with the soil, the heat transfer term can be defined comparably simple due to its constant volume fraction and transfer area. Due to the assumption that the ice core grows in the middle of the pore and that always a liquid layer will be between the ice and the solid matrix, there is no heat transfer between the rock and the frozen water (Bittelli et al., 2003;Lebeau & Konrad, 2012). For the liquid and gaseous phase, a homogeneous distribution is assumed, and both components interact with the soil matrix based on their volume fractions (Heinze & Blöcher, 2019). In such a case, the interaction between the soil matrix and the gaseous and liquid components can be defined as The contact area of the ice core with the liquid and the gaseous phase is calculated accordingly with respect to their volume fractions.
As the contact area between air and liquid is negligible, there is no heat transfer term  lg 0 E Q .

Freezing and Melting
Based on the previously introduced capillary tube model, the temperature of phase transition and the numbers and radii of pores affected by phase transition at given temperatures is calculated. The procedure described here is analogous to the work of Watanabe and Flury (2008). For a given water retention curve, the associated distribution of capillary radii is determined by dividing the water retention curve into equally spaced intervals of the water content.  (Watanabe & Flury, 2008). The effective porosity of the soil can then be expressed as As stated above, in a porous media an ice core will always be surrounded by a thin liquid water layer of thickness E d due to adsorptive forces along the capillary walls dependent on the Hamaker constant E A (J), which describes the average interactions between molecules at the solid-fluid interface due to van der Waals forces (Dash et al., 1995;Watanabe & Flury, 2008). Ice generation is constrained by the Gibbs-Thompson effect as only aggregates of sufficient size are stable to grow. For a cylindrical tube, the critical radius GT E r depends on the ice-water interfacial free energy  iw E (Watanabe & Flury, 2008). Subsequently, in a tube of radius J E R an ice core can be generated when the condition   J GT E R d r is true (Watanabe & Flury, 2008).
Larger tubes are prone to develop an ice core quickly and an annulus of unfrozen water will exist in these tubes. known temperature, the amount of pore space that could possibly freeze based on the above outlined conditions can be determined. Based on the thermal conditions of the liquid -and possibly existing frozen -water, the amount of actually freezing water can be calculated.

Hydraulic Flow Model
Water infiltration into a unsaturated soil, frozen or not, is described by the well-known Richards equation (Richards, 1931), e.g., in its head based form (Farthing & Ogden, 2017): In the presence of ice, the hydraulic conductivity is substantially reduced depending on the volumetric ice content. The hydraulic conductivity is then given as (Dall'Amico et al., 2011;Hansson et al., 2004)  with the hydraulic conductivity at saturated conditions sat E K (m/s) and an impedance factor Ω E (−), that is usually set to 7 (e.g., Dall'Amico et al., 2011;Hansson et al., 2004;Peng et al., 2016). The soil moisture capacity is given by (van Genuchten, 1980)  Besides the water infiltration, change in liquid water volume fraction can also be caused through melting or freezing as discussed above. This is included into the model through the source term E M , which can be associated with  , / l f m E  as the density contrast between liquid water and ice is neglected. Additionally, the amount of maximum achievable liquid water content , l sat E  is reduced or increased during melting or freezing. If the maximum achievable water saturation of an unfrozen soil is given by  , l sat E , in the presence of ice the maximum water fraction is given as:

Numerical Implementation
The possible change in volume fractions due to flow and phase change requires specific attention since hydraulic flow and temperature equations are strongly coupled. Further, unphysical temperatures for the separate phases, such as ice temperature above the melting temperature, may not occur. To address these problems, a splitting algorithm was introduced in Dall' Amico et al. (2011), to separate advective mass and heat fluxes from diffusive fluxes and phase change processes. In analogy to this procedure, a predictor-corrector scheme separating mass transport, diffusive heat transport and phase change.
In a first step, the volume fractions of gas and liquid will be calculated based on the hydraulic model outlined in Section 2.7. This includes the solution of Equations 30-34. Based on the volume fraction of the frozen water, the radius of the ice core in the pores f E R can be calculated for heat transfer (Equation 21).
Using all volume fractions and flow velocities, the heat transfer coefficient and the heat transfer area can be calculated as described in Sections 2.3 and 2.4, given in Equations 18, 20, and 24. Subsequently, the amount of heat transferred between the phases can be calculated (Equations 25-28), which then allows the calculation of the phase temperatures by solving the heat equations for each phase given in Section 2.2 (Equations 5, 7, 10, and 12) neglecting the phase transition terms in Equations 10 and 12. This calculation is the predictor step for the phase temperatures as phase change is considered next.
In a second step, the phase change as well as its influence on temperature and volume fractions is considered. If in the predictor step the calculated temperature  f ph E T T , then the amount of available energy per The amount of frozen water that can be melted with this available energy per discrete time step E dt (s) is calculated as with the latent heat of fusion f E L (J/kg). Based on the relationship given above    , f l m E   . In a corrector step, the phase temperature is increased by the latent heat described in Equation 13, which sets the frozen water temperature back to the temperature of phase transition ph E T . Freezing of liquid water is described in a similar way and outlined as pseudo-code in Algorithm 1. Once the temperature of the liquid water l E T is lower than the temperature of phase change ph E T , any additional decrease in thermal energy is transferred into a phase change similar to Equation 35 The change in liquid volume fraction due to freezing is then given as Also during the corrector step, the temperature of the liquid water is increased to ph E T in addition to the latent heat described in Equation 14 and the temperature of the frozen water is increased due to the addition of newly frozen water. Subsequently, the influence of phase change on the hydraulic properties such as saturation and hydraulic conductivity need to be considered. With an in-or decrease of maximum achievable saturation (Equation 34), there is a subsequent change of the effective saturation. The described approach is similar to a rock-mechanical approach with a dynamic porosity due to pore pressure changes in a partly saturated rock (Heinze et al., 2015;Rutqvist et al., 2002). The temperature of phase transition ph E T is determined based on the amount of frozen water, as depending on the derived pore size distribution only fractions of the available pore water will freeze at a specific temperature. The derived relationship between ph E T and the volume available to freezing or thawing has been assumed piece-wise linear for a suitable numerical implementation.
The governing equations are implemented in MATLAB ® using an one-dimensional explicit finite difference scheme. The Richards equation is implemented in the head based form as due to the presence of ice a strong heterogeneity in hydraulic conductivity might occur (Farthing & Ogden, 2017). Mass conserva-tion is secured through a small time step (  0.01 E dt s), which is required by the thermal calculations as well. The code of the Richards equation has been compared to Hydrus 1D solutions. Several infiltration scenarios on various column length, infiltration duration and soil type have been simulated and compared to the respective Hydrus 1D solutions. Both numerical simulations agree well with small deviations toward dry conditions. The thermal advection is implemented with a third-order upwind scheme. The thermal diffusion is calculated with a forward in time -centered in space finite difference (FTCS) scheme. Additional constraints on the maximum allowed time step occur through the heat transfer terms (Heinze & Hamidi, 2017). The results shown in this work have been tested to be independent from the numerical discretization within a sufficient range of spatial (0.05-1.5 cm) and temporal resolution (0.001-0.05 s). The results presented here are conducted with 0.1 cm grid distance. The used parameters to describe the four involved phases are given in Table 1. The value of the dry soil was estimated based on the work of Yun and Santamarina (2008). Thermal boundary conditions are formulated in terms of heat extraction following the approach of Hansson et al. (2004) and Dall'Amico et al. (2011). Hydraulic boundary conditions follow experimental setups of either constant hydraulic pressure or no flow. Δ 0.015 E , which is in a similar range, or slightly better, than the divergence between experimental results and simulations of freezing processes in column experiments previously presented (cf. Dall'Amico et al., 2011;Hansson et al., 2004;Peng et al., 2016). Additionally, intrinsic soil heterogeneity, which is not included in the numerical model, might cause some smaller divergences between simulation and experiment toward the bottom of the column. The upward movement of the liquid water toward the freezing front due to capillary suction is clearly visible, which has been experimentally observed several times (e.g., Hansson et al., 2004). This effect is more significant for the sandy loam than for the silt loam (Figure 1). The divergence between experimentally obtained and calculated volume fractions is smaller for the silt loam than for the sandy loam (Figures 1a and 1c). The difference between temperatures from simulation and experiment is more significant for silt loam than for sandy loam and for the temperatures in the unfrozen part of the column. While the simulation results follow an almost linear trend from the bottom toward the freezing front, the measured temperature distribution in the silt loam is quite non-linear. This is most likely due to the soil heterogeneity, that can also be seen in water content distribution and initial temperatures. For the sandy loam, there is a substantially better agreement between the calculated and measured temperatures. In the simulation, separate phase temperatures have been cal-  (Figures 1b and 1d). The amount of calculated ice for the last time step (48 hr) in the sandy loam is significantly higher in the simulated results than for the measured values (Figure 1c). However, the freezing fronts, which are determined as the 0 C limit, agree fairly well for both soil types (Figure 2). The freezing front is slightly underestimated in the simulation for the silt loam toward the end and at the start for the sandy loam.
Subsequently to the 48 hr of freezing, a liquid water column of 150 mm height with a temperature close to  1.0 E C was placed on top of the soil columns. The hydraulic head was kept constant and the experiments lasted until drainage was observed at the bottom. Again, liquid water content and temperature were monitored continuously. The simulation results reproduce the experimental findings quite well but the divergence between experimentally observed and calculated liquid water content is larger for the thawing than for the freezing experiments (Figure 3). Similar to the freezing experiments, the agreement between simulation and experiment is better for the silt loam than for the sandy loam. The small-scale heterogeneity has a much more significant influence on the results during the thawing than during the freezing experiments. The maximum liquid water content in the sandy loam columns is much higher in the upper part than in the lower part (Figure 3c). Through the assumption of a homogeneous sand column due to a lack of more detailed information about the experimental setup, such small scale heterogeneity is not reproduced in the numerical simulation. However, position and amount of the remaining frozen part of the soil is very well reproduced in all experiments, though in the simulation slightly closer toward the top than observed in the experiment (Figure 3). The infiltration resistance of the frozen soil is underestimated by the model, which is especially obvious for the sandy loam with a larger discrepancy between initial and saturated water content compared to the silt loam experiments (Figure 3). The initial continuation of freezing at the bottom end of the existing frozen part of the soil is reproduced in the simulations accordingly to the experimental results. The measured temperatures agree fairly well between laboratory measurements and simulations with a tendency of higher temperatures in the simulations (Figures 3b and 3d). Again, for the simulation mixture temperatures are calculated based on the volume fractions and the separate phase temperatures. The almost constant temperature value close to the melting point along the depth of the remaining frozen part of the soil is well reproduced for both soil types.

Differences in Phase Temperatures
With the newly derived model the separate phase temperatures can be studied in detail (Figures 1b, 1d,,  and 3d). In the freezing experiment, the separate phases have a common temperature initially (Figure 1). Due to the different thermal properties the thermal evolution between phases can differ if the heat transfer between the phases is not sufficiently large to maintain a local thermal equilibrium. Varying heat transfer area and heat transfer coefficient due to dynamic volume fractions inside the pores, can alter the heat transfer over time. Additionally, through freezing-point depression the liquid water can take on temperatures below  0 E C in smaller pores due to adhesive forces at the soil surface (e.g., Watanabe & Flury, 2008). As outlined in Section 2.6, frozen and unfrozen pores can be calculated based on a derived pore size distribution depending on the van Genuchten parameters. Depending on the calculated pore size distribution, the maximum possible volume fractions for frozen water increases for decreasing temperatures (Figure 4). Although both tested soil types have a comparably small mean pore size, the vast majority of the pore water would freeze at  0 E C. The freezing-point depression is even stronger for finer materials with larger clay content, and vanishes for coarser materials, such as sand.
Besides the initial thermal equilibrium between phases, a difference between phase temperatures evolves around the freezing front of around  0.3 E C after 6 hr of freezing ( Figure 5). With growing freezing depth, this difference between the separate phases becomes smaller and finally vanishes. Rock and air cool quicker than the liquid water, due to their lower heat capacity. Due to its higher thermal conductivity the rock phase is slightly colder than the air. Ice and liquid water have very similar temperatures due to their comparably large heat transfer between each other, which is in agreement with previous studies assuming thermal equilibrium between ice and water using the Clapeyron equation (e.g., Dall'Amico et al., 2011;Kurylyk & Watanabe, 2013). Through the energy required by the phase transition, and due to their high heat capacity, liquid and frozen water do not cool down as quickly as rock and air.
During the melting experiment with liquid water intrusion from the top, the liquid water is substantially warmer (more than  1 E C) at the top than all other phases (Figure 6c). Due to the numerical necessity to avoid very small volume fractions, the ice can be slightly warmer than the melting point and even slightly warmer than the liquid phase. However, this happens only by a small margin of less than  0.1 E C around the melting point (Figure 6b). Due to their lower heat capacity rock and air are the two phases warming quicker than ice and liquid water over many parts of the soil column based on the new boundary conditions and the phase transfer to the warmer inflowing water. Once melting started around the remaining frozen part of the soil, ice and liquid water have very similar temperatures (Figure 6b). However, a clear difference between Figure 2. Comparison between the freezing front as obtained from the experimental data presented in Watanabe and Kugisaki (2017) and using the presented model. the water temperatures (frozen and liquid) and the rock and the air phase of around  0.1 E C can be observed (Figures 6c and 6d). In regions without the presence of ice, a thermal equilibrium between the phases is reached (Figure 6).

Comparison to Soil Column Freezing Experiment
A common freezing experiment used by various researchers for code comparison is presented by Hansson et al. (2004). In this experiment four identical soil columns of a length of 20 cm and a diameter of 8 cm were filled with sandy loam. An initial homogeneous water distribution of ∼33% was achieved. From an initial temperature of  6.7 E C, the top was cooled to   6.0 E C while all other sides were thermally and hydraulically isolated. The cooling procedure was stopped after 12, 24, and 50 hr for one column each and the columns were sliced in 1 cm thick pieces for which the total water content (liquid + frozen) was determined. The previously presented model was applied to this scenario for comparison. No flow hydraulic boundary conditions were applied at top and bottom as well as a thermal isolation boundary condition at the bottom. At 3.5 10 E , 0.011 and 0.016. Around the freezing front a small divergence in phase temperatures can be observed in the simulations similar to the results presented above.

Coarse Grained Material
The experimental setup of Watanabe and Kugisaki (2017) has been applied to a sand with a porosity of 0.43, a saturated hydraulic conductivity of   5 8.25 10 E m/s and a residual water content of 0.045. The van Genuchten parameters were set to   14.5 E 1/m and  2.68 E n . This results in a calculated mean pore radius of   4 2.3 10 E m, which is approximately a factor of 10 larger than for the sandy loam and silt loam used above (cf. Table 2). The initial saturation was set to 0.3. All other physical and numerical parameters as well as boundary and initial conditions remained unchanged compared to the infiltration experiments of  Watanabe and Kugisaki (2017) described above. Due to the coarser material and the larger pore size, the heat transfer behavior in sand differs from the previous experiments as heat transfer area and heat transfer coefficient to the soil are reduced.
The freezing behavior in sand is different to the previous examples because of the differences in thermal and hydraulic parameters. The freezing front after 48 hr of freezing is much shallower and the temperature difference between liquid and frozen water is substantially larger in sand than for the loam studied above (Figure 8). Further, the frozen water is substantially cooler than the liquid water because once ice forms at the top the high thermal conductivity of the frozen water cools the ice from the top, while the other phases cool less quickly. The temperature difference between the phases is several degree Celsius (Figure 8). In the parts of the soil in which ice is present and a bit further down, the water is within the numerical range necessary around zero degree Celsius, indicating that further freezing would occur for longer freezing times (Figure 8). Small amounts of liquid water remain in the small pores in the frozen part of the column and there is no observable difference between air and rock temperature, similarly to the examples above due to freezing point depression and the heat transfer between rock and air. Once the infiltration starts from the top, the rock and air warm by heat transfer from the liquid and so does the ice. A drop in water temperature can be seen in the region where the water infiltrated in the upper centimeters and replaced the remaining air. Due to the small volume fractions of ice generated during 48 hr of freezing, here a snapshot after 2 hr of water infiltration is shown. For further times, the ice temperature is approaching the rock/gas temperature and melting starts once the temperature is above the melting point.

Reproduction of Experimental Results
The presented model has been compared to laboratory experiments reproducing freeze and thaw behavior presented in Watanabe and Kugisaki (2017) and Hansson et al. (2004). In the freezing experiments, the model agrees well with the experimental data, similar to other models (e.g., Dall'Amico et al., 2011;Hansson et al., 2004;Peng et al., 2016). Based on the different parameterization used in this model with respect to the hydraulic parameters strongly coupled to heat transfer, the obtained parameters to match the experiments of Hansson et al. (2004) slightly diverge from values used by other studies, though a spread in used parameters as well as a variation in thermal boundary conditions is common across studies (cf. Daanen et al., 2007;Dall'Amico et al., 2011;Hansson et al., 2004;Painter, 2011;Peng et al., 2016). The parameters derived in this work are consistent across soil types as both silt loam experiments result in comparable values and in agreement with literature values for those soil types. Reproduction of thawing experiments with water infiltration has not been found in the literature so far. The model shows a slightly larger divergence between simulation and experiment during thawing, probably due to the increased impact of small scale heterogeneity and of the hydraulic conductivity. The impedance factor Ω E describing the reduction of hydraulic conductivity in dependence on the present ice content has been set to 7 based on previous studies (Dall 'Amico et al., 2011;Hansson et al., 2004). Apparently this model underestimates the blocking of fluid pathways through the soil in the current model as outflow occurs earlier in the simulation than in the laboratory experiments. In this context, the extraordinarily high saturated hydraulic conductivity for the sandy loam should be noted as taken from Watanabe and Kugisaki (2017), which is significantly larger than expected for this soil type. Such a high hydraulic conductivity might be explained through the Figure 7. Comparison between total water content (frozen + liquid) for the experimental setup presented in Hansson et al. (2004) between experimental data (symbols) and simulation results (lines). filling process of the column and of course strongly influences the results of the thawing experiment. In general, the thawing processes seem much more sensitive to hydraulic changes than the freezing processes. With the involved processes of simultaneous freezing toward the bottom and thawing toward the top of the soil column and the strong coupling to the flow behavior, the thawing simulations are much more complex than the freezing simulations.
The experiments of water infiltration into the frozen soil show the methodological strength of the newly developed model. The separate phase temperatures can be included into the model without any further assumption, e.g., with respect to thermal parameters or the mixture temperature. For example, the boundary conditions can be described consistently without any limiting assumption regarding an instant thermal equilibrium between the phases at the inflow. Further, any change in volume fractions results in an automatic change of the thermal state through the heat transfer functions linking the separate phases in the presented model. A strong difference between melting and freezing point by several degree Celsius has been observed for Walla Walla silt loam (Bittelli et al., 2003). In the same experiments the freezing/ melting process occurred rapidly at one temperature (Bittelli et al., 2003). This is not in agreement with the experimental data of Watanabe and Kugisaki (2017), which shows that some part of the soil sample froze at colder temperatures than others. This is a known phenomenon depending on the pore size distribution and has been acknowledged as such in the presented model. Small differences between model and experiment remain because a quantification of this process is difficult and strongly dependent on soil type and its small scale heterogeneity. Due to a lack of data for the tested soils, a shift between freezing and melting point as proposed by Bittelli et al. (2003) has been neglected in this study.
During the freezing experiments, a substantial difference between phase temperatures can be observed around the freezing front due to the different thermal properties of the phases and the phase transition of the water. During the thawing experiments, the temperature difference between the phases is initiated through the infiltration of warmer liquid water as well as through the melting process of the ice. In both kinds of experiments, freezing and melting, rock and air temperatures are quite similar but distinct from water and ice temperatures, which again are quite similar to each other for small grained materials but can diverge by several degree Celsius in coarse grained materials. It is known, that soils with a larger mean grain size experience a more substantial divergence between phase temperatures between solid rock/soil and water (Gossler et al., 2020;Heinze & Blöcher, 2019). The calculation of mixture temperatures for comparison of experimental and numerical results therefore might cause some divergence, especially for the coarser sandy loam with a larger mean grain size.
The simulations of freezing and thawing in sand indicate that the temperature of the infiltrating water in coarse grained materials is less relevant for the melting process than ambient temperatures, such as surface temperature, because the heat transfer between the liquid water and the other phases is small compared to the effect of the heat conduction within the ice body. For a transfer toward field applications the specific laboratory setup needs to be considered as all phases are cooled at the top. In a natural setting thermal gradients in the soil and thermal boundary conditions might be different possibly emphasizing the effect of heat transfer for the melting of the ice. Strong differences in phase temperatures question the assumption of thermal equilibrium between phases during freezing and melting processes in coarse materials. However, the results might be considered with care and require further experimental investigation as the heat transfer coefficient becomes of great importance in coarse grained materials. The formula of Wakao and Funazkri (1978) has been used in this work due to its broad applicability but has been critically discussed for its applicability in aquifer systems (Gossler et al., 2020).

Future Research Perspective
The presented model is the first true multi-phase heat transfer model with more than two phases and applications to freezing and thawing processes. The derived mathematical scheme can easily be extended or altered to other multi-phase scenarios. The application to water infiltration into frozen soil is just one very obvious application of heat transfer in a truly non-equilibrium situation due to the initially diverging phase temperatures. The experiments reproduced are rather simple. The next step are inhomogeneous initial water content distributions, e.g., from true infiltration experiments. However, to this point such experiments have not been conducted and field observations have a comparably coarse sensor network. To achieve a more realistic representation of the soil, also preferential pathways and macropores should be considered. There is a heat transfer concept for dual domain models, which could be adopted (Heinze & Hamidi, 2017). It might be reasonable to parameterize the heat transfer coefficient E h differently along preferential pathways as E h is known to be different for the porous matrix and fractures in consolidated rock. The dependence of the heat transfer on flow velocity is more significant in fractures than in a porous media (Frank et al., 2021) and during the melting of ice cores in fractures melting rate strongly depends on the width of the available flow channel and flow velocity (Heinze, 2020). Similar effects could occur within micro-cracks and preferential pathways within the soil. Also for the porous matrix, the used empirical heat transfer coefficient has never been validated for ice and might need further testing. From the derived thermal model it would be straight forward to consider a possible density contrast between frozen and liquid water. However, a suitable flow model has not been presented so far.
A study of relevance of LTNE models in infiltration scenarios into frozen soil is required to constrain the limits of the LTE assumption. First results in this work indicate that for small grained soils the LTE assumption might be a valid assumption. From previous work, such as Hamidi et al. (2019); Heinze and Blöcher (2019); Gossler et al. (2020), several conclusions can be transferred for infiltration scenarios into frozen soil. It can be expected, that the temperature differences between phases persists longer for (a) larger grain sizes, (b) higher flow velocities and (c) larger initial temperature differences between phases. The effect of larger grain sizes has been studied in this work and agrees with findings from saturated media. The scenarios to be studied examining the limit of the LTE assumption should include complex temperature distributions in the soil as might occur during freeze and thaw cycles when upper soil layers get heated but deeper layers remain frozen.

Conclusions
A multi-phase heat transfer model has been derived and applied to freezing and thawing processes in partly saturated soils. The model has been implemented in a numerical framework and compared to freezing and thawing experiments with fluid intrusion in silt loam and sandy loam. The newly derived model matches the experimental results while demonstrating the benefits of an explicit description of the heat transfer between phases and partly questions the assumption of LTE in coarse grained soils. The boundary conditions of a comparably warm liquid water infiltrating into a frozen soil can be adequately described and the initial local thermal non-equilibrium can be reproduced. A difference between phase temperatures has been observed during freezing and thawing. The thawing experiments are significantly more complex to simulate because of simultaneous freezing and thawing within the soil column and the increased difference in phase temperatures due to the water infiltration. In the future, the model will provide insights into thermally controlled infiltration processes during snow melt and rainfall on frozen ground, potentially improving the understanding of rainfall triggered snow avalanches and flood events by critically evaluating the LTE assumption in such scenarios and therefore possibly improving the estimation of liquid and frozen water content in the soil. The derived model allows a consistent and physically precise formulation of water infiltration into soil without limiting assumptions regarding the thermal state of the involved phases.

Data Availability Statement
All data used for model comparison has been digitized from the publication of Watanabe and Kugisaki (2017) and Hansson et al. (2004). The Matlab code is deposited under https://doi.org/10.5281/zenodo.5542256.