Simulating water and heat transport with freezing and cryosuction in unsaturated soil: Comparing an empirical, semi-empirical and physically-based approach

Freezing of unsaturated soil is an important process that influences runoff and infiltration in cold-climate regions. We used a simple numerical model to simulate water and heat transport with phase change in unsaturated soil via three different approaches: empirical, semi-empirical and physically based. We compared the performance and parameterization of each approach through testing on three experimental datasets. All approaches reproduced the observed unsaturated freezing process satisfactorily. The empirical cryosuction equation used in this study managed to capture observed cryosuction with a fixed empirical parameter value. The semi-empirical version therefore does not require calibration of a specific frozen soil related parameter. In view of simplicity, small computational demand and accurate performance, all three approaches are suitable for implementation in landuse schemes, catchment scale hydrological models, or multi-dimensional thermo-hydrological models.


Introduction
There is often a substantial flux of water in a catchment during the springtime snowmelt period in high latitude and high altitude regions ( Rango and DeWalle, 2008 ). In many cases, part of the soil has become frozen during winter, thereby reducing the soil permeability and infiltration capacity ( Stähli, 2006 ). The combination of snowmelt, rainfall and soil frost has the potential to cause flooding ( Woo, 2012 ). Frozen soil also influences snowmelt groundwater recharge and springtime contaminant transport by altering pathways and soil water retention times ( Evans et al., 2018;French et al., 2002;French and Binley, 2004;Hayashi et al., 2013;Ireson et al., 2013 ). In order to predict and understand increased flood risk, groundwater recharge and contaminant transport in areas with frozen soil, it is therefore necessary to have sufficient understanding of the mechanisms that control soil freezing, thawing and infiltration.
The part of a soil that generally undergoes seasonal freezing and thawing extends from a few centimeters to about a meter or several meters below the surface ( Loranger et al., 2017 ;Lundberg et al., 2016 ;Hayashi, 2013 ). This mostly comprises the unsaturated zone where moisture content and soil temperature respond to atmospheric dynamics licly available numerical models include a freezing module for HYDRUS-1D ( Š im ů nek et al., 1998 ), a beta version of SUTRA named SUTRA-ICE ( McKenzie et al., 2007 ) and the atmosphere-plant-soil models SHAW ( Flerchinger and Saxton, 1989 ) and COUP ( Jansson and Karlberg, 2004 ). SUTRA-ICE does not currently include the process of cryosuction and HYDRUS-1D freezing only works for unsaturated conditions. Kurylyk and Watanabe (2013) noted that the history of frozen soil model development has been characterized by inconsistencies in nomenclature and methodology, in part due to different geotechnical and hydrological backgrounds. It remains unclear how different mathematical expressions and models for unsaturated soil freezing processes compare to each other in their ability to accurately represent an unsaturated frozen soil. An exception to this is the paper by Ren et al. (2017) in which the outcomes of different SFC equations are fitted to measurements on four different frozen soils. Regarding the reduction of permeability of frozen soil, there is debate about the use of a flow impedance factor based on ice content ( Mohammed et al., 2018 ). Furthermore, previous models are often tested only on a single experimental dataset which did not include all relevant variables such as ice content and soil temperature ( Mohammed et al., 2018 ;Kurylyk and Watanabe, 2013 ). As a result, confidence in model accuracy remains limited, and it remains unclear how the model would be parameterized for a different soil than the soil used in the experiment (e.g. the impedance factor).
Another unresolved question is whether an empirical approach towards unsaturated soil freezing and cryosuction could be adequate and how such an approach would compare to more physically-based models using the phase-change temperature-pressure relationship ( Kurylyk and Watanabe, 2013 ). This question is relevant, as few multi-dimensional hydrological models, catchment-scale models or land-surface schemes have adopted approaches for unsaturated soil freezing, likely due to the complexity and often associated numerical instability of physically based simulation of unsaturated soil freezing.
In this study, we compare three different approaches for unsaturated soil freezing. These represent a fully empirical, a semi-empirical and a physically-based approach and entail different combinations of previously developed equations. The aim of this study is to compare and evaluate the performance and parameterization of each approach. Datasets from three experiments are used for testing, namely those of Mizoguchi (1990) , Watanabe et al. (2012) and Zhou et al. (2014) . The latter two datasets contain measurements of ice content, unfrozen water content and soil temperature, which allows us to test for all relevant variables. We include a discussion of current points of debate concerning the governing equations for unsaturated frozen soil dynamics, as the different approaches rely on previous insights and developments in both empirical and physically-based equations.

Soil freezing curve
Not all soil moisture freezes at the same temperature due to the dependency of the freezing point of water on pressure ( Hayashi, 2013 ). Matric pressure increases with decreasing unfrozen water content as ice replaces water. This has been found to resemble the process of soil drying ( Koopmans and Miller, 1966 ). The Soil Freezing Curve (SFC) gives the relationship between subzero temperature and unfrozen water content for a given soil; the lower the temperature, the less water remains unfrozen ( Ren et al., 2017;Spaans and Baker, 1996 ).
The following form of the Clausius-Clapeyron relationship relates pressure to temperature ( Williams and Smith, 1989 ): where P w is the pressure of liquid water (Pa), T is the temperature ( °C), L f is the latent heat of fusion (Jkg − 1 ), and V w is the specific volume of liquid water (m 3 kg − 1 ). The change in pressure of the ice phase is assumed to be zero. Furthermore, it requires the assumption of thermal, mechanic and phase equilibrium during freezing. These constraints in the context of frozen soil modeling are discussed in Ma et al. (2015) . Expressions for the SFC have been formulated that combine the Clapeyron relationship with existing equations for soil water retention during drying/wetting. An example is from Zhang et al. (2016) who combined the Clapeyron relationship with the van Genuchten equation: where uf is the unfrozen water content (m 3 m − 3 ), res is the residual total water content (m 3 m − 3 ), sat is the saturated total water content (m 3 m − 3 ), w is the density of water (kgm − 3 ), T is the temperature ( °C), T 0 the freezing point of water ( °C) and a vg , n vg and m vg are fitted model parameters.
Several empirical SFC expressions have been proposed that do not make use of the Clausius-Clapeyron relationship. The simplest approach is a linear piecewise relationship between temperature and unfrozen water which has been found to reasonably approximate measured SFC datapoints ( McKenzie et al., 2007 ;Kurylyk and Watanabe, 2013 ). Anderson and Tice (1972) developed an empirical power curve to express unfrozen water content as a function of temperature: where d is the dry density of the soil (kgm − 3 ) and and are model parameters. The power curve has been used in several models and parameter values for a wide range of soils have been established, as well as a method to derive the parameters from the specific surface area of a soil ( Kurylyk and Watanabe, 2013 ).
In the model SUTRA-ICE, an empirical exponential function for the unfrozen water content is implemented. This function has been suggested by McKenzie et al. (2007) based on work by Lunardini (1988) : where w is a fitting parameter. It is difficult, however, to know the correct value for w in Eq. (4) as very few soils have been parameterized through experimental work ( Ren et al., 2017 ). Hence, it is a parameter that often requires calibration. Ren et al. (2017) provides a verification of four different SFC expressions on experimental data (the SFC expressions of Eq. (2) , 3 and 4 included). The study shows that all three expressions work well, with some slightly better than others depending on soil type.

Reduction of hydraulic conductivity
Several studies have shown that a frozen soil has reduced permeability, and that initial water content is important among other factors (e.g. Pittman et al., 2019;Watanabe and Osada, 2017 ), but the mathematical representation of the reduction of soil hydraulic conductivity due to freezing remains a point of debate. The lack of consensus might exist due to different reasons: 1) it is difficult to accurately measure the permeability of a frozen soil, hence experimental evidence is scarce ( Azmatch et al., 2012 ;Watanabe and Osada, 2017 ); 2) ice growth in soil voids may have different effects on permeability for different soil types and different initial moisture contents and therefore, a single type of reduction may not apply to all soils and moisture conditions; 3) soil freezing rate -a factor usually not taken into account -likely plays a role as it determines how ice crystallizes in soil voids ( Azmatch et al., 2012 ); 4) soil structure can be altered by ice lens growth and it is hard to predict how this affects soil hydrological properties ( Mohammed et al., 2018 ).
Two different causes may be distinguished for a reduced permeability of a frozen soil, assuming fixed total water content and no frost heave: 1) with decreased unfrozen water content, porewater connectivity decreases, analogous to air saturation of a soil ( Stähli, 2006 ;Lundberg et al., 2016 ;Hayashi, 2013 ); 2) with an increased ice content, the pore geometry changes due to an increase in the 'solid ice' soil constituent, certain flow pathways become blocked and effective pore size is reduced ( Zhang et al., 2016 ;Azmatch et al., 2012 ). It has been argued that the second set of factors are not different from permeability reduction upon increasing air saturation and that therefore the unfrozen water content of a soil can be used to predict the permeability of frozen soil accurately (e.g. Newman and Wilson, 1997 ;Watanabe and Flury, 2008 ). A hydraulic conductivity function derived from the soil water retention curve is therefore used by most models ( Kurylyk and Watanabe, 2013 ). Some studies, however, suggest that the unfrozen water content approach overpredicts the permeability at high ice contents. These studies suggest that ice blockage of flow paths should be taken into account (e.g. Jame and Norum, 1972 ;Lundin, 1990 ;Hansson et al., 2004 ;Zhang et al., 2007 ;Shoop and Bigl, 1996). For this reason, an impedance factor has been introduced to correct for the additional reduction in permeability due to ice ( Kurylyk and Watanabe, 2013 ).
The impedance factor approach by Jame and Norum (1972) , for example, considers ice content in addition to unfrozen water content: where K f (ms − 1 ) is the hydraulic conductivity of (partly) frozen soil, K uf (ms − 1 ) is the hydraulic conductivity based on the unfrozen water content, i is the volumetric ice content (m 3 m − 3 ) and E i is a dimensionless empirical factor of impedance due to ice blockage of pores. Zhang et al. (2007) found that including an impedance factor ( E i = 17) was crucial to accurately describe observed moisture transport in the soil during freezing. Several variations of the impedance equation ( Eq. (5) ) have been developed (e.g. Taylor and Luthin, 1978 ), but without much added clarity of a most reliable expression. Shoop and Bigl (1997) suggested the following equation to determine E i based on data collected from nine different non-cohesive soils: where K sat is the saturated hydraulic conductivity (ms − 1 ). The impedance factor approach has been criticized by Newman and Wilson (1997) to be arbitrary and not physically based. It is also suggested that it is simply a parameter to correct a numerical model for overestimated cryosuction based on the Clausius-Clapeyron approach, as the extreme hydraulic gradient at the frozen fringe is difficult to simulate in a numerical model ( Miller, 1980 ). In addition, Watanabe and Flury (2008) claimed that the use of an impedance factor should be unnecessary when an accurate relative permeability function is used related to unfrozen water content, though they also stated that it is a topic of ongoing research. Several researchers have come up with alternative approaches due to the criticisms surrounding the impedance factor ( Kurylyk and Watanabe, 2013 ). An example is Watanabe et al. (2010) , who successfully employed the dual porosity model of Durner (1994) to simulate the hydraulic conductivity of a frozen unsaturated silt soil. The closure and opening of macropores has also been recognized as an important factor determining water flow in frozen soil by several studies (e.g. Mohammed et al., 2018 ;Holten, 2019 ;Pittman et al., 2020 ;Demand et al., 2019 ), which requires a dual porosity approach to incorporate in a numerical model. In our study, we will assess the need for the impedance factor in a physically-based approach as a means to represent a 'soil ice saturation heterogeneity' in an otherwise homogeneous soil (i.e., water flowing from unfrozen to frozen soil).

Cryosuction
There is no consensus on the most suitable representation of cryosuction in a soil freezing model ( Kurylyk and Watanabe, 2013 ). The process is often associated with extreme hydraulic gradients requiring high spatial discretization to successfully simulate in a numerical model. In some approaches, matric suction is increased based on ice content through an empirical equation to simulate the effect of cryosuction on water redistribution (e.g. Zhang et al., 2007 ).
The following empirical equation for cryosuction was originally developed by Kulik (1970): where t (m) is the total matric pressure of the soil (including the effect of cryosuction), uf (m) is the matric pressure due to total water content without the effect of possible ice present and C k is an empirical factor that represents the effect of ice on matric pressure. A more novel suggestion is that soil properties should be altered through a change in soil water retention parameters and soil porosity (the ice-free volume) ( Noh et al., 2012 ). In most previous numerical models however, matric suction is a function of subzero temperature based on the Clausius-Clapeyron relationship ( Eq. (1) ). Dall'Amico et al. (2011) uses the following version of the Clausius-Clapeyron relationship to determine matric pressure: where g the gravitational constant (ms − 2 ). In case of unsaturated conditions, the expression is modified to account for the matric suction resulting from air saturation in the pores. Most other cryosuction expressions using the Clausius-Clapeyron relationship do not adjust for different total water saturation levels ( Dall'Amico et al., 2011 ). In this paper, we will use both the empirical and physically-based (Clausius-Clapeyron) approaches to cryosuction for comparison.

Methods
All abbreviations or symbols used for variables and parameters are listed in Table 1 . Fixed parameters used in the model are listed in Table 2 .

Model structure and assumptions
A numerical model was used that calculates heat and water transport for a one-dimensional soil profile with any number of layers. Vertical discretization of the layers was set to one centimeter with uniform soil properties for the entire column. Heat transport occurs at the top and bottom boundary, where a fixed temperature boundary can be set. There is no water flow possible across the model boundaries. The mathematical equations are solved through explicit difference calculation for the fluxes between soil volumes. This method increases numerical stability and simplicity, but it requires high temporal discretization to maintain accuracy. Yang et al. (2009) successfully used a similar numerical method in their simulations of unsaturated flow governed by the Richards equation. In our case, it appeared to be an adequate model construct for the purpose of 1D "laboratory " type soil column simulations with freezing. The model excludes osmotic processes and density changes of unfrozen water. Also, porosity and soil structure do not change with ice saturation as ice pressure is assumed to be constant.

Unsaturated flow
Flow between soil volumes is governed by the Richard's equation, here presented in its 1-dimensional form ( Richards, 1931 ): where totw is the total water content (liquid and ice; m 3 m − 3 ), z is the elevation (m), tot is the total matric pressure (m) and K uf ( uf ) is the hydraulic conductivity (ms − 1 ) as a function of the unfrozen water content (and also ice content if an impedance factor is used).
Gravitational acceleration ms − 2 C k Empirical parameter for cryosuction -R i Empirical parameter for flow reduction - Table 2 Physical constants and parameter values used in the model.

Parameter Value Units
Cryosuction C k 1.8 -Thermal conductivity (water) c w 0.6 Wm − 1 K − 1 Thermal conductivity (ice) c i 2.14 Wm − 1 K − 1 Thermal conductivity (air) c a 0.024 The matric pressure of a soil volume is a state variable that depends on total water content, in our case given by the van Genuchten equation ( van Genuchten, 1980 ): where a ( m − 1 ) and n are model parameters. Hydraulic conductivity, the K ( uf ) term in Eq. (8) , is calculated with the following equation, which derives from the relationship between the relative permeability function of the Mualem -van Genuchten model ( Mualem, 1976 ;van Genuchten, 1980 ), the saturated hydraulic conductivity and the unfrozen water content: where K uf is the hydraulic conductivity of the soil (ms − 1 ). The effect of temperature on hydraulic conductivity by affecting the viscosity of water is neglected.

Energy exchange
There are three forms of energy exchange in the model that govern the energy balance: thermal conduction, advection and latent heat flux. These are expressed in the following general energy balance equation: , a is the volumetric air content, a is the density of air (kgm − 3 ), H a is the specific heat of air (Jkg − 1 K − 1 ), is porosity, s is the density of solid soil (kgm − 3 ), H s is the specific heat of solid soil (Jkg − 1 K − 1 ), ̄ is the average thermal conductivity of the soil (Wm − 1 K − 1 ), v is the flow velocity of unfrozen water (ms − 1 ), L f is the latent heat of fusion (Jkg − 1 ) and z is the elevation (m). Thermal dispersion is assumed to be negligible for heat transfer in small-scale unsaturated soil ( Liu et al., 2014 ;Jouybari et al., 2020 ).

Freezing approach: empirical version
The empirical approach quantifies the effect of cryosuction on flow and the effect of matric potential on freezing point depression without an underlying physical explanation. For the soil freezing curve, we use the exponential equation of McKenzie et al. (2007 ;Eq. (4) ). Only the fitting parameter w is needed to approximate the freezing curve of a soil with this equation. The rationale for an empirical approach is that the model can easily be calibrated to better fit data. In addition, it is not affected by assumptions such as thermal equilibrium phase change which is the case for the physically-based SFC. Also, an empirical SFC can more easily be applied to non-colloidal soils via calibration.
To simulate cryosuction in the empirical version, we make use of Eq. (7) . Cryosuction entails both matric potential changes and the resulting flow of water. We assume the empirical cryosuction expression represents the observed upward flow correctly, not the matric potential changes itself; hence it incorporates a possible flow impedance effect due to ice content. This approach circumvents the numerical instability associated with the extreme hydraulic gradient at the frozen fringe when using the Clausius-Clapeyron relationship. By combining the exponential SFC equation and the empirical cryosuction expression, we thus have an empirical approach requiring two parameters, w for the SFC and C k for cryosuction. The question is whether C k can be generalized for a variety of soil types, or if it should be soil type specific. We will investigate this by testing the empirical (and semi-empirical) approach on three different experimental datasets, later described. We will also include a small sensitivity analysis to the parameter C k .

Freezing approach: semi-empirical version
The semi-empirical approach that we use contains the same empirical expression for cryosuction, but it is combined with a physically-based expression for the SFC. We use Eq. (2) , shown above and originally formulated by Zhang et al. (2016) , in which the Clausius-Clapeyron relationship is incorporated into the van Genuchten soil water retention expression. The freezing curve of the soil is hence determined by the common van Genuchten soil water retention parameters a and n . The only unsaturated freezing related parameter to be calibrated for the semiempirical version is therefore C k for cryosuction.

Freezing approach: physically-based version
For the physically-based version, we use the expressions from Dall' Amico et al. (2011) to determine cryosuction. First, the freezing temperature is determined by: where T * (K) is the freezing point of water at the current matric pressure based on total water content, totw , and g is the gravitational acceleration (ms − 2 ). Subsequently, matric pressure including the effect of ice is determined by the following expression: If T ≥ T * , the equation collapses to t = totw . The van Genuchten based SFC ( Eq. (2) ) is used to determine the soil freezing curve. Hence, both matric potential and freezing point depression are based on the physical relationship between temperature and pressure. However, similar to other studies ( Kurulyk and Watanabe, 2013 ), we found that the extreme hydraulic gradient at the frozen fringe led to a strong overprediction of upward flow. Therefore, we developed a simple approach to solve this using the ice impedance factor combined with the concept of a soil structure discontinuity in the context of spatial discretization.
In any numerical model solving differential equations, soil water and heat transport need to be discretized in time and space. If water flows from a discretized location A to a location B, the hydraulic driving force between these points and the hydraulic conductivity of location A determine the flow rate. This premise would hold if the soil represented by location B has the same hydrological properties as location A. In case of ice in the soil however, the assumption of soil homogeneity cannot hold. If location B would be partly frozen, certain flow paths could be blocked as larger pores freeze first. The inflow rate is no longer dependent on the hydraulic conductivity of location A alone. Therefore, a special hydraulic conductivity reduction is needed for the frozen location that receives soil water.
We limited the flow rate of water to ice-filled soil volumes with the following formula, developed by Zhao et al., 2013 : where K m is the maximum hydraulic conductivity (ms − 1 ) for flow towards a frozen soil volume and R i is the impedance factor for flowrate reduction due to ice in the soil pores. We will assess to which extent R i varies with soil type.

Numerical method
To ensure proper functioning of the numerical method, we successfully compared the output of the model to the established models SUTRA ( Voss and Provost, 2002 ) and HYDRUS-1D ( Š im ů nek et al., 1998 ) for unsaturated flow and heat transport during nonfrozen conditions and for fully saturated frozen conditions (only SUTRA). The resulting comparisons are given in the appendix.

Experimental data
We make use of the experimental data of Mizoguchi (1990) , Watanabe et al. (2012) and Zhou et al. (2014) to test the approaches for frozen unsaturated conditions. All experimental soil columns were insulated and closed systems. Frost heave was not observed in any of the experiments. The soil properties and boundary conditions of the experiments are listed in Table 3 . The parameters w and R i were manually calibrated for each soil type. If a soil parameter was unknown, this was manually calibrated as well (mentioned in Table 3 ). For the datasets of Watanabe et al. (2012) and Zhou et al. (2014) , we compare the empirical and the van Genuchten-based SFC to the measured unfrozen water contents at several subzero temperatures. The simulated SFC's for Mizoguchi (1990) will also be shown, but without measured unfrozen water content for comparison. Regarding the C k parameter, we investigated whether a single empirical value could capture the cryosuction observed in the experiments. We tested for a range of C k values that would provide a water distribution that visually matched the experimental results; this range was between a value of 1 and 3. We include a small sensitivity test of the empirical cryosuction parameter to show how we established a single value for C k for all experiments. Mizoguchi (1990) used a 20 cm high soil column filled with Kanagawa sandy loam. It was frozen from the top with a temperature of − 6 °C, while the soil had an initial temperature of 6.7 °C. Only total water content was measured in this experiment. Several authors used the dataset of Mizoguchi for model testing, such as Dall'Amico et al. (2011) and Hansson et al. (2008). We include the model results from Dall'Amico et al. (2011) for this experiment. By comparing the three versions to their model results, we can assess how well the different mathematical expressions compare with the approach of Dall'Amico, who used the Clausius-Clapeyron relationship to simulate cryosuction ( Eq. (8) ) with an equation splitting method combined with the numerical Newton method. Furthermore, we will include the output from SUTRA-ICE for this experimental setup. SUTRA-ICE is a multidimensional saturated and unsaturated water and heat transport model that includes the depression of the freezing point of water, but not the process of cryosuction. By comparing the results of our model to SUTRA-ICE, we can assess the effect of cryosuction on simulated soil physics such as soil temperature and ice content.
Watanabe et al. (2012) did a freezing experiment on a soil column of 35 cm deep. The soil comes from the A horizon of a weeded fallow field and is characterized by a high porosity. Three different columns were prepared, each with a different initial water content: 0.31, 0.38 and 0.46. The column was brought to a homogenous temperature of 3.5 °C and subsequently frozen from the top with a temperature of − 6.2 °C. The bottom of the column was in contact with a temperature element set to 2 °C. Measurements were done after 48 h and included total water content and unfrozen water content. Zhou et al. (2014) performed a freezing soil column experiment on a sieved glacial till (sieve size 0.063 mm). The 24 cm high column initially had a temperature of 3 °C and the experiment was performed with two different initial water contents, 0.16 and 0.325. For these two different initial water content setups, freezing temperatures of − 4 and − 4.7 °C were respectively applied at the top. Total water content, unfrozen water content, ice content and soil temperature were measured 24 and 72 h after freezing started. The parameters for saturated hydraulic conductivity and thermal conductivity are unknown for this experiment and are therefore estimated based on soil type and subsequently slightly adjusted via manual calibration.

Comparison of simulated SFC to measurements
The different SFC curves corresponding to the simulations and the different experiments are shown in Fig. 1 . With the experiment of Watanabe et al. (2012) , the empirical SFC equation did not perform well without adjustment. When the measured residual water content during drying is used, the empirical SFC severely underpredicts the unfrozen water content ( Fig. 1 ). Therefore, we changed the unfrozen residual water content for the empirical SFC to 0.07 (SFC_exp2) instead of 0.006 (SFC_exp). After this adjustment, it is still apparent that the SFC underpredicts the unfrozen water content. As shown in Fig. 1 , a higher w value for the SFC does not solve the problem (SFC_exp3), as it creates a too steep decrease in unfrozen water content with decreasing temperature. In general, the empirical SFC displays a near linear relationship between unfrozen water content and temperature that quickly reaches the residual water content at relatively high subzero temperature. The van Genuchten based SFC on the contrary, displays a more gradual decrease in unfrozen water content at lower subzero temperatures, and even at − 8 °C the residual water content is not reached.
With the experiment of Zhou et al. (2014) the empirical SFC and the van Genuchten SFC represent the measured unfrozen water contents well. In the experiment of Mizoguchi (1990) no unfrozen water content measurements were performed to compare the results with, but it is also clear that the empirical SFC predicts reaching a residual water content at a much higher temperature than the van Genuchten based SFC. A crucial difference between the two SFC approaches thus seems apparent in the high matric potential range, corresponding to low temperatures (below − 1 °C).

Calibration results
The calibration results ( Table 4 ) show that w varied for the different soil types used in the experiments. The C k parameter was kept at a fixed value of 1.8. To obtain accurate results, it was important to calibrate the impedance factor, R i , for each soil type specifically as otherwise the predicted cryosuction was noticeably under-or overpredicted. Mizoguchi (1990) The comparisons with the measurements of the experiment by Mizoguchi (1990) and the simulated outputs of the three approaches in this study, SUTRA-ICE and the numerical model by Dall'Amico et al. (2011) are shown in Fig. 2 With all three approaches used in this study, the distribution of total water content after 12, 24 and 50 h is in good agreement with measurements, each performing slightly better than the model of Dall'Amico et al. (2011) . The variation amongst the empirical, semi-empirical and physically-based approaches is small. The depth to which cryosuction affects the water distribution in all simulations seems to align well with the measured water content profile.

4.3.
As can be expected, SUTRA-ICE does not reproduce the cryosuctionbased increase of total water content within the frozen zone. Also shown in Fig. 2 are the ice contents and soil temperature profiles of our model and SUTRA-ICE. The frozen zone, as well as the zero-degree temperature isotherm, is deeper in the SUTRA-ICE simulation. Although the vertical extent of the frozen zone is larger, the ice content is lower in the SUTRA-ICE simulation compared to our model simulation.  Table 4 Calibration results of frozen soil related parameters. The parameters w and C k are used for the empirical version, the semi-empirical version requires only the C k parameter and the physically-based version requires only the impedance factor, R i . Parameter Mizoguchi (1990) Watanabe et al. (2012 w (SFC) 0.5 0.1 0.5 C k (cryosuction) 1.8 1.8 1.8 R i (impedance factor) 9 11 12 4.4. Watanabe et al. (2012) The comparison between simulated output the measurements in the experiment of Watanabe et al. (2012) are shown in Fig. 3 . All versions of the model perform reasonably well to simulate the observed total and unfrozen water contents for the three different initial water content setups. Some deviation can be seen with the empirical version, as the unfrozen water content is underpredicted. The unfrozen water content drops sharply above 23 cm elevation in the empirical version, while the other versions and the measurements display a more gradual decrease in unfrozen water content above this point. Zhou et al. (2014) Figs. 4 and 5 compare the output of the simulations to the experiment of Zhou et al. (2014) for two different initial water contents and at two different measurement times. In general, model results are in good agreement with observed values. All three versions predict the increase of water in the frozen zone and other output variables with reasonable accuracy. The variation amongst the different versions is noticeable but small. The model, however, struggles to capture some of the observed cryosuction in the lower initial water content setup, as there is a strong deviation of the total water content in the lower section of the freezing front after 3 days ( Fig. 4 ). The model in general predicts a mild increase of total water content with depth in the frozen zone, while the experimental data suggests that there was a steep increase in water content at the freezing front (depth 12 -14 cm) after the first day. The physicallybased version performed slightly better in this case than the empirical and semi-empirical approach. This could imply that the empirical cryosuction expression should take total water content into account to determine the effect of ice on matric potential. For this reason, we tested with an adapted cryosuction equation for the semi-empirical approach that is dependent on total water saturation:

4.5.
where C i represents the effect of ice on matric pressure. We found a value of 0.8 for C i to match the observed cryosuction across the different experiments. The result of substituting Eq. (7) with Eq. (16) is included in Fig. 4 (CF_S2). For all other cases, results remained roughly

Sensitivity to the empirical cryosuction parameter
In addition to the calibrated value for Ck of 1.8 the following values were tested; 1, 1.5, 2.5 and 3. We used these Ck values in the simulations of several of the experimental setups with the empirical ap- Fig. 3. Measurements of total water content and unfrozen water content (Obs) after 48 h in the soil freezing column experiment by Watanabe et al. (2012) with different initial water contents (IWC) and the output of the different model approaches (CF -E, -S and -P indicating empirical, semi-empirical or physically based approaches, respectively). Porosity of the soil is 0.617. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) proach. The results are shown in Fig. 6 . It is clear that between a Ck of 1 to 3, the results are within reasonable accuracy compared to measurements, but the midpoint covering the observed cryosuction in all experiments appears to be around the calibrated Ck value of 1.8. The band of results differs for each experiment, as it is narrow below the frozen zone for Watanabe et al. (2012) and Mizoguchi (1990) , but wide for Zhou et al. (2014) .

General model performance
The model is capable of simulating three different unsaturated soil freezing experiments with reasonable accuracy. These experiments have different initial water content, freezing intensity and soil type. The three versions of the model predict the penetration depth of the freezing front with 1-centimeter accuracy in most cases. The ice content and total water contents are predicted reasonably well with an accuracy of about 0.05 (m 3 /m 3 ). It can thus be concluded that a simple freezing extension of a common soil water and heat transport model based on the Richard's equation is an adequate means of simulating freezing soil. The differences in accuracy amongst the empirical, semi-empirical and physically based approach are small, but noticeable. There is no approach that consistently performs better when considering all cases.

Empirical approach
The empirical SFC equation ( McKenzie et al., 2007 ) combined with the empirical cryosuction equation ( Kulik, 1978 ) circumvents the use of the more complex Clausius-Clapeyron relationship while the results show it can adequately capture the freezing process of unsaturated soil. A promising result is that a fixed parameter value for cryosuction (Ck), set to 1.8, simulated cryosuction well compared to measurements. Even though the experiments only represented three soil types, it suggests that in most cases -at least within the textural range of a sandy loam, silt loam and loamy silt, no soil type specific calibration would be required. This eases the applicability of this approach to a wide range of situations. Only one case showed underpredicted cryosuction, when the soil started with a low initial water content of 0.16. Adjusting the cryosuction equation ( Eq. (7) ) to include a dependency of the cryosuction effect on total water content ( Eq. (16) ), improved the fit slightly. Such a dependency on total water content could be expected because there are significant changes in matric potential with changing unfrozen water content in the low water content range of a soil water retention curve; i.e. the effect of changes in ice content on matric potential could be stronger at low total water content.
The exponential SFC requires soil type specific calibration of the empirical SFC parameter (w). It also became clear that though the empirical exponential SFC works well in most cases, it tends to underpredict the unfrozen water content in a fine soil -in this specific case, a loamy silt. The unfrozen residual water content is quickly reached at relatively high subzero temperature (between 0 and − 1 °C). The residual unfrozen water content should be set to a different value than the unsaturated residual water content to avoid a drop to residual unfrozen water content too quickly. This unfortunately limits the applicability of the exponential SFC to soils with weak soil water retention if an accurate unfrozen water content, and thus ice content, is desired. In soils with strong soil water tension, such as clays and silts, the matric potential will increase significantly at lower unfrozen water contents, leading to a strong depression of the freezing point of water. For this reason, the exponential SFC could not capture the SFC of a silt loam well.

Semi-empirical approach
The advantage of the semi-empirical approach is that it uses the empirical cryosuction equation, but it relies on the Clausius-Clapeyron relationship to determine the freezing point of water. This means the depression of the freezing point of water is thus based on the well-known physics of phase change. The van Genuchten based SFC in combination with the empirical cryosuction equation worked well in all cases considered, except for the low initial water content case. The van Genuchten based SFC captures the measured SFCs better than the empirical SFC equation and it provides a more realistic drop of unfrozen water content at low subzero temperatures. With this approach, residual unfrozen water content is reached at significantly lower subzero temperatures com-  Zhou et al. (2014) in their soil column freezing experiment and the output of the different approaches (CF -E, -S and -P indicating empirical, semi-empirical or physically based approaches, respectively) for total water content, temperature, unfrozen water content and ice content after 72 h with an initial water content of 0.325. Porosity of the soil is 0.47. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) pared to the exponential SFC. Another benefit of the semi-empirical approach is that it does not require a special parameter for the SFC. Hence, only the Ck parameter is needed. Since the empirical cryosuction equation predicts cryosuction well with a fixed value for Ck, it is implied that this approach does not need calibration of a freezing related soil parameter.

Physically-based approach
By limiting the flow of unfrozen water to the frozen zone with an impedance function, we were able to simulate cryosuction based on the Clausius-Clapeyron relationship in a simple numerical model with good accuracy. In all cases results were in good agreement with observations; only at the onset of freezing, the impedance factor tends to limit the flow towards the freezing front slightly too strongly in some cases. The main advantage of the physically-based approach is that it relies on the underlying physics of the temperature-pressure phase change relationship for both the freezing point depression and cryosuction, and it should therefore be widely applicable. The disadvantage is that, at least in our case, it requires reduction of flow to the frozen zone via an empirical impedance factor.
The impedance factor is an empirical, soil type specific parameter, and as we noted there is debate about the validity of its use. We used the impedance factor as a representation of reduced soil porosity in the frozen soil, akin to a soil heterogeneity. In our approach, the impedance parameter had to vary for the different soil types in the experiment to capture the cryosuction process well. An alternative approach could be the dual porosity model used by Watanabe et al. (2010) . A next step in a physically-based approach would likely involve changing the soil hydrological properties based on ice content, as it can be expected that the soil water retention parameters and saturated hydraulic conductivity would change with increasing ice content ( Noh et al., 2012 ), but this would require further experimental study.

The importance of simulating cryosuction
By comparing the results of our model to SUTRA-ICE, which does not simulate cryosuction, we could identify the effect it had on total  Mizoguchi (1990) , Watanabe et al. (2012) with initial water content 0.38 and Zhou et al. (2014) with initial water content 0.325. Ck was set to 1, 1.5, 1.8, 2.5 and 3. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) water content, ice content and temperature. Cryosuction logically increases the total water content in the frozen zone. Accordingly, it slows down the freezing front as the heat capacity of the upper soil is increased and more energy is used for phase change. Consequently, in SUTRA-ICE, the freezing front progresses faster and the frozen zone becomes larger, though with a lower ice content. This implies that without cryosuction, the frozen zone has a higher permeability and more space for accommodating infiltrating water. Of practical concern in flood risk assessment, cryosuction thus strongly affects the infiltration capacity of a soil. Our study however is based on medium to fine textured soils in the range of loam and silt. It can be expected that for coarse soils such as sand, in which gravitational drainage significantly precedes cryosuction, the effect of cryosuction on total water content is strongly diminished. In very fine soils such as clays, cryosuction has been found to play a limited role as well due to very low hydraulic conductivity preventing water redistribution ( Miller, 1980 ).

Conclusion
In this study, a simple 1D numerical model is used to simulate water and heat transport with phase change in unsaturated soil via three different approaches: empirical, semi-empirical and physically based. These approaches constitute new combinations of previously developed equations. The fully empirical approach uses an empirical exponential function for the soil freezing curve (SFC) and an empirical function for cryosuction. We found that the empirical SFC underpredicts unfrozen water content for fine soils at low subzero temperatures (below − 1 °C), leading to a loss of accuracy. The advantage of this approach is that it does not require accurate soil water retention parameters to work and that it does not rely on the assumptions associated with the Clausius-Clapeyron approach (such as thermal, phase and mechanic equilibrium).
The semi-empirical approach uses the van Genuchten soil water retention model combined with the Clapeyron relationship for the SFC, while cryosuction is based on the empirical equation. Since the cryosuction equation worked well with the same parameter value for a sandy loam, silt loam and loamy silt, the main advantage of this approach seems to be that calibration of a soil type related freezing parameter can be avoided. The van Genuchten.based SFC also performs better at temperatures below − 1 °C, as it more accurately links matric pressure to the freezing point depression. Therefore, if correct unfrozen water content is desired for freezing soils with significant fine particle content, the semi-empirical approach is preferred. The physically based approach was used in our numerical scheme by regarding a frozen soil volume as a soil discontinuity. Similar to other studies, it was also necessary to use an impedance function in order to not overpredict upward flow ( Kurylyk and Watanabe, 2013 ). The main advantage of this approach is that it is more physically based and therefore should be more widely applicable to different freezing circumstances, although it requires soil type calibration of the impedance factor.
The suggested approaches are useful for large-scale models in the simulation of frozen unsaturated soil. Depending on available soil data and model scale, an empirical, semi-empirical or physically-based approach could be preferred. Correct simulation of ice and water content is relevant in case of determining soil infiltration capacity and possible contaminant pathways. In addition, by simulating cryosuction correctly, it will be possible to predict zones of increased total water content which are thus susceptible to ice lensing and frost heave. Further modeling studies could investigate soil freezing and thawing dynamics in relation to actual infiltration of rain-and meltwater, which has received little attention. An important topic would for example be freezing of infiltrating water, which would lower infiltration capacity but add significant amounts of energy as latent heat.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix Comparison of model basic energy and water transport with SUTRA and HYDRUS1D
SUTRA was originally developed for saturated and unsaturated groundwater flow of variable-density with solute or energy transport, but it has been extended to include the freeze-thaw process (beta version used, called SUTRA-ICE). It uses Galerkin finite element and finite difference methods. HYDRUS-1D is a model used for water, heat and solute flow in variably saturated media. Its numerical solver is based on Galerkin type linear finite elements. Both models have been in the public domain for decades and have been applied to numerous case studies.
For the comparisons with our model, we set up a 1D model domain with a vertical extent of 1 meter, a no-flux (energy or water) bottom boundary (at a depth of 1 m) and constant top boundary conditions (at an elevation of 0 m). We chose different scenarios to test the unsaturated flow, heat conduction, latent heat flux and advection. Parameters and other conditions are listed in Table A for the six scenarios considered (parameter constants used are listed in Table 2 in the main text). Re-

Fig. C.
Comparison of simulated temperatures after 24 and 72 h of the model used in this study and SUTRA-ICE with a top boundary of 1 °C and a uniform initial temperature of − 1 °C (scenario 4). Only 0 to 50 cm depth shown of total 100 cm depth in the simulation. CF is the empirical version of the model used in this study. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. D.
Comparison of simulated water contents after 1, 2, 3, 4 and 7 days of the model used in this study, HYDRUS1D and SUTRA-ICE with a top boundary sourceflow of 2 mm per hour (scenario 5). CF is the empirical version of the model used in this study. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. E.
Comparison of simulated temperature after 12, 24 and 48 h of the model used in this study and SUTRA-ICE with a top boundary sourceflow of 5 mm per hour with a water temperature of 5 °C and no specified temperature boundaries (scenario 6). CF is the empirical version of the model used in this study. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) sults are shown in Figs. A-E . These show a nearly identical fit between output of HYDRUS1D or SUTRA-ICE for all cases considered compared to the model used in this study (empirical approach used, designated with "CF ").