A robust and energy-conserving model of freezing variably-saturated soil

. Phenomena involving frozen soil or rock are important in many natural systems and, as a consequence, there is a great interest in the modeling of their behavior. Few models exist that describe this process for both saturated and unsaturated soil and in conditions of freezing and thawing, as the energy equation shows strongly non-linear characteristics and is often difﬁcult to handle with normal methods of iterative integration. Therefore in this paper we propose a method for solving the energy equation in freezing soil. The solver is linked with the solution of Richards equation, and is able to approximate water movement in unsaturated soils and near the liquid-solid phase transition. A globally-convergent Newton method has been implemented to achieve robust convergence of this scheme. The method is tested by comparison with an analytical solution to the Stefan problem and by comparison with experimental data derived from the literature.


Introduction
The analysis of freezing/thawing processes and phenomena in the ground is important for hydrological and other land surface and climate model simulations (e.g.Viterbo et al., 1999;Smirnova et al., 2000).For example, comparisons of results from the Project for Intercomparison of Land Surface Parameterization Schemes have shown that the models with an explicit frozen soil scheme provide more realistic soil temperature simulation during winter than those without (Luo et al., 2003).Freezing soil models may be divided into three categories: empirical and semiempirical, analytical, Correspondence to: M. Dall'Amico (matteo@mountain-eering.com) and numerical physically-based (Zhang et al., 2008).Empirical and semiempirical algorithms relate ground thawingfreezing depth to some aspect of surface forcing by one or more experimentally established coefficients (e.g.Anisimov et al., 2002).Analytical algorithms are specific solutions to heat conduction problems under certain assumptions.The most widely applied analytical solution is Stefan's formulation, which simulates the freezing/thawing front using accumulated ground surface degree-days (either a freezing or thawing index) (Lunardini, 1981).Numerical physicallybased algorithms simulate ground freezing by numerically solving the complete energy equation, and in natural conditions they are expected to provide the best accuracy in simulating ground thawing and freezing (Zhang et al., 2008).However, this approach has difficulties, especially regarding the treatment of phase change, which is strongest in a narrow range of temperatures near the melting point, and thus represents a discontinuity that may create numerical oscillations (Hansson et al., 2004).Furthermore, the freezing process has a profound effect also to the water fluxes in the soil, as it changes the soil hydraulic conductivity and induces pressure gradients driving water movements.Therefore, a coupled mass and energy system is needed to simulate both the thermal and hydraulic characteristics of the soil.
The objectives of the paper are: (1) to revisit the theory of the freezing soil in order to provide the formulation for the unfrozen water pressure, which can accomodate variablysaturated soils; (2) to outline and describe a numerical approach for solving coupled mass and energy balance equations in variably-saturated freezing soils, based on the splitting method; (3) to provide an improved numerical scheme that: (i) is written in conservative way, (ii) is based on the globally convergent Newton scheme, and (iii) can handle the high non-linearities typical of the freezing/thawing processes.
Published by Copernicus Publications on behalf of the European Geosciences Union.
The algorithm is tested against the analytical solution of unilateral freezing of a semi-infinite region given by Neumann (Carslaw and Jaeger, 1959;Nakano and Brown, 1971), in order to test convergence under extreme conditions, and against the experimental results published by Hansson et al. (2004), in order to test the coupled water and heat flow.

Freezing-soil models
One of the first attempts to include soil freeze-thaw into a numerical model is the work of Nakano and Brown (1971), who assumed the advance of the freezing front as a moving boundary.They applied the analytical solution of Carslaw and Jaeger (1959) to a porous soil, and introduced the effect of an artificial freezing zone of finite width between the frozen and the unfrozen parts; this avoids the problem of shock wave propagation in the transition between the frozen and thawed state, typical of the "freezing front" assumption.No water flow, i.e. no mass balance equation, is taken into account and the energy balance is expressed through the definition of an apparent heat capacity during phase change, as proposed by Lukyanov and Golovko (1957).Harlan (1973) is probably the first to attempt to solve coupled mass and energy balance equations for freezing soil by making an analogy between the mechanism of water transport in partially frozen soils and those in unsaturated soils.He solves the system on a homogeneous rigid porous medium through a fully implicit finite difference scheme, where the unknowns are temperature and soil water potential; phase change in the water balance appears in the source/sink term.He also uses the apparent heat capacity formulation in the energy balance.The results show that the freezing process induces the movement of both heat and mass from warm to cold regions, causing the moisture content in the unfrozen soil zone to decrease sharply towards the freezing front.Soil texture and initial moisture conditions seem to be crucial in affecting the availability and mobility of water.Fuchs et al. (1978) developed a theory of soil freezing applicable to unsaturated conditions with solute presence in the soil.They considered that solutes tend to depress the freezing point temperature and modified the relationship between temperature, moisture content and apparent thermal properties of the soil.Phase change was taken into account with the apparent heat capacity theory, and the water flow contribution is accounted for in the apparent thermal conductivity.Therefore, the simultaneous heat and water transport equations result in a merged single differential equation for heat.Jame and Norum (1980) further developed the model of Harlan (1973), and highlighted that the effect of mass transfer on the thermal state of soil is an important factor to be considered.Newman and Wilson (1997) proposed a theoretical formulation for unsaturated soils using soil-freezing and soil-water characteristic curve data to combine the heat and mass transfer relation-ships into a single equation for freezing or frozen regions of the soil.Christoffersen and Tulaczyk (2003) constructed a high-resolution numerical model of heat, water, and solute flows in sub-ice stream till subjected to basal freeze-on.They proposed a formulation of the equilibrium relation without assuming zero ice pressure through the full version of the generalized Clapeyron equation (Mikkola and Hartikainen, 2001), which enables them to model segregation ice onto the freezing interface and so develop stratified basal ice layers.Hansson et al. (2004) introduced a new method for coupled heat transport and variably saturated water flow using the Richards' equation.They accounted for water flow due to gravity, pressure gradient and temperature gradient, both for liquid and vapor phase.McKenzie et al. (2007) proposed the freezing module of SUTRA for saturated conditions.Daanen et al. (2007) developed a 3-D model of coupled energy and Richards equation to analyze the effects that lead to the formation of non-sorted circles, an example of a relatively stable patterned-ground system.Their model is very similar to Hansson et al. (2004), but it differs in the linkage between ice content and temperature.Watanabe (2008) also used a similar model to Hansson et al. (2004).He reproduced directional freezing experiments on silty soil and compared this to experimental data.
Basing on these studies, some important aspects in freezing soil algorithms may be highlighted: 1. the degree of soil saturation is often not specified, and thus soil is tacitly assumed fully saturated; 2. the relation between liquid water content and temperature, usually derived by the "freezing = drying" assumption (Miller, 1965), is not linked to the soil saturation degree; 3. the method used to solve the mass and energy balance is seldom accurately explained, especially regarding the approach to the coupled differential equations.

Mass and energy in the soil
A soil is usually composed of soil particles (subscript "sp") and water.The content and phase (liquid or solid) of water vary with soil temperature and saturation degree: in unsaturated-unfrozen conditions liquid water (subscript "w") and air are present, whereas in unsaturated-frozen conditions ice (subscript "i") also occurs (Fig. 1).
Let us assume the following hypothesis: (1) rigid soil scheme, i.e. the volume V c is constant and no volume expansion during freezing is allowed.Therefore the density of water ρ w and of ice ρ i (kg m −3 ) must be considered equal, because otherwise the change of the water density upon freezing would lead to unrealistically large gauge pressures that cannot be converted into an expansion of the soil matrix, due to the lack of a mechanical model; (2) the "Freezing = drying" assumption (Miller, 1965;Spaans and Baker, The Cryosphere, 5, 469-484, 2011 www.the-cryosphere.Basing on these studies, some important aspects in freezing soil algorithms may be highlighted: 1. the degree of soil saturation is often not specified, and thus soil is tacitly assumed fully saturated; 2. the relation between liquid water content and temperature, usually derived by the "freezing=drying" assumption (Miller, 1965), is not linked to the soil saturation degree; 3. the method used to solve the mass and energy balance is seldom accurately explained, especially regarding the approach to the coupled differential equations.1996), which implicitly implies that: (i) the freezing (thawing) water is like evaporating (condensing) water; (ii) the ice pressure is equal to the air pressure; (iii) the water and ice content in the soil may be related to the soil matric potential according to the soil water retention curve, e.g.Brooks and Corey (1964), Clapp and Hornberger (1978), Gardner (1958) or Van Genuchten (1980); (iii) the volume V c is representative of soil volumes above the water table or at shallow depths below the water table, where the water pressure may be considered similar to the atmospheric pressure.
Let us define ψ w0 (m) as the matric potential corresponding to the total water content (liquid and ice), and ψ (m) as the matric potential corresponding to the liquid water content only.Indicating with f () the water retention curve function relating water content to matric potential, it is: θ w = f (ψ) and v = f (ψ w0 ), where θ w and θ i (−) are respectively the volumetric water and ice contents.v := θ w + θ i (−) is the total volumetric content of ice and water, as explicated in Appendix A. Since the liquid water content cannot be larger than the total water content and f is normally represented as a monotone function, it follows that: where θ s (−) the soil saturated water content.When the soil volume V c is saturated, ψ w0 = 0 and v = θ s .When the soil is unfrozen ψ = ψ w0 and θ w = v .
The energy content of the soil volume V c is represented by the internal energy U (J m −3 ), which, neglecting the energy of air and excluding the work of volume expansion passing from the liquid to the frozen state, can be calculated as the sum of the internal energy of the soil particles, ice and liquid water: In particular, with regard to a reference temperature T ref (K), each of the above terms becomes: where T (K) is the temperature, L f is the latent heat of fusion (J kg −1 ), ρ sp (kg m −3 ) is the density of soil particles, and c sp , c i and c w (J kg −1 K −1 ) are the specific thermal capacity of soil particle, ice and liquid water, respectively.Usually T ref is set to the melting temperature T m at atmospheric pressure, equal to 273.15 K. Summing the components in Eq. ( 3) in a more compact form gives: where ) is the total thermal capacity of the soil volume.

The conservation law
Water and energy fluxes through the boundaries of V c cause variations of water content and internal energy.Neglecting the ice flux within the soil, the mass conservation in a soil volume may be written as (Harlan, 1973): where M w and M i (kg) are the mass of water and ice respectively present in V c , the symbol "∇•" refers to the divergence operator and S w (s −1 ) is a sink term due to evapotranspiration.After some calculations one obtains: where m (−) is the total water content as derived by the fraction of the total mass of water and ice M t = M w + M i (kg) and the maximum mass of water that would fit in V c , as explained in Appendix A. J w (m s −1 ) is the water flux within the soil following the Darcy-Buckingham formulation: "∇" is the gradient operator, K H (m s −1 ) is the hydraulic conductivity, p (m) is the liquid water pressure head (equal to the soil matric potential when soil is unsaturated), and z f (m) is the elevation with respect to a reference.Analogous to the case of water content controlled only by drying processes, the hydraulic conductivity is dependent on the soil matric potential ψ associated with liquid water (Mualem, 1976), with its maximum value K H,sat when water is unfrozen and soil saturated.However, the presence of ice may significantly reduce the hydraulic conductivity due to the apparent pore blockage effect exerted by ice.This is accounted for by further www.the-cryosphere.net/5/469/2011/The Cryosphere, 5, 469-484, 2011 reducing the hydraulic conductivity by an impedance factor smaller than 1 and equal to 10 −ωq (Zhao et al., 1997;Hansson et al., 2004), where ω is a coefficient and q is the ice fractional content given by θ i θ s −θ r , where θ r is the residual water content.Zhao et al. (1997) and Hansson et al. (2004) both set ω equal to 7.
The energy conservation in a soil volume V c may be written as (Fuchs et al., 1978): where: 1. S en (W m −3 ) represents a sink term due to energy losses; 2. G (W m −2 ) is the conduction flux through the volume boundaries.According to the Fourier law, G is defined as: where is the thermal conductivity of the soil within V c that depends on the thermal conductivity of the soil particles (λ sp ) and on the proportions of ice and water and the temperature T (Johansen, 1975;Farouki, 1981).
3. J (W m −2 ) is the heat advected by flowing water: Eventually the equations to be solved are the mass conservation (Eq.6) and energy conservation (Eq.8), which represent a system of two differential equations.If the soil is unfrozen, the unknowns are ψ and T , and so the system is solvable; if the soil is frozen, the unknowns are: ψ, ψ w0 and T , and so, in order to solve the system, it is necessary to provide a further closure relation relating the temperature and soil water pressure head.

Pressure and temperature under freezing condition
When the soil is unfrozen, the liquid water content (θ w ) coincides with the total water content.At freezing conditions, however, it depends also on temperature.According to the "freezing = drying" assumption, the ice pressure is set to null, and so the generalized Clapeyron equation becomes: where p w (Pa) is the water pressure.This equation states that the variation of water pressure head during phase change is only dependent on water temperature.
If the soil volume V c is saturated, we can assume that the water pressure at the soil surface equals the atmospheric pressure p a .In order to derive the pressure p w at freezing conditions, where the temperature is lower than the melting temperature at atmospheric pressure (T m ), one should integrate Eq. ( 11) as follows: The left hand term of Eq. ( 12) may be approximated as: Usually the atmospheric pressure is set to zero (p a = 0).Combining Eqs. ( 12) and ( 13) one obtains: Considering that p w ≤ 0 and ρ w • g • ψ w = p w , where g (m s −2 ) is the acceleration due to gravity, one can find the liquid water matric potential (ψ) subject to freezing conditions: If the soil is unsaturated, the surface tension at the water-air interface decreases the water melting temperature to a value T * < T m .Let us call p w0 the pressure of water corresponding to the pressure head ψ w0 .The value of T * may be found integrating Eq. ( 11) in temperature from T m to the unknown T * , which means integrating in pressure head from the value p a = 0 to the degree of saturation given by p w0 .One obtains: which, considering the approximation of Eq. ( 13), provides the melting temperature T * at unsaturated conditions: When T ≥ T * , the soil is unfrozen whereas when T < T * , the soil is under freezing conditions.The above relation, applying the Laplace equation for determining the radius r from ψ w0 , becomes similar to Eq. ( 5) in Watanabe and Mizoguchi (2002).
The liquid water pressure p w depends on the intensity of freezing condition provided by T , and may be found integrating Eq. ( 11) in pressure from p w0 to p w and in temperature from T * to T .One obtains: and consequently: The  Freezing curve for pure water and various soil textures, according to the Van Genuchten parameters given in Table 1.
5 The decoupled solution: splitting method The final system of equations is given by the equations of mass conservation (Eq.6) and energy conservation (Eq.8): The previous system is a function of T and ψ w0 and can be solved by the splitting method, as explained in Appendix B.
In the first half time step, the Richards' equation is solved and the internal energy is updated with only the advection contribution.In the second half, no water flux is allowed, which makes the volume a closed system, and the internal energy is updated with the conduction flux in order to find the new temperature and the new combination of water and ice contents.

Fig. 2.
Freezing curve for pure water and various soil textures, according to the Van Genuchten parameters given in Table 1.
The above equation is valid for T < T * : in fact, when T ≥ T * , the freezing process is not activated and the liquid water pressure head is equal to the ψ w0 .Equations ( 17) and ( 19) collapse in Eq. ( 15) for a saturated soil (i.e.ψ w0 = 0).Thus the formulation of the liquid water pressure head ψ(T ) under freezing conditions, valid both for saturated and unsaturated soils, becomes: which can be summarized using the Heaviside function H ( ) as: If the soil water retention curve is modeled according to the Van Genuchten (1980) model, the total water content becomes: where θ r (−) is the residual water content.The liquid water content θ w becomes: Equation ( 23) gives the liquid water content at sub-zero temperature and is usually called "freezing-point depression equation" (e.g.Zhang et al., 2007 andZhao et al., 1997).Differently from Zhao et al. (1997), it takes into account not only the temperature under freezing conditions but also the water 1.0 0.0 4 × 10 −1 2.50 sand 0.3 0.0 4.06 × 10 −3 2.03 silt 0.49 0.05 6.5 × 10 −4 1.67 (Schaap et al., 2001) clay 0.46 0.1 1.49 × 10 −3 1.25 (Schaap et al., 2001) depressed melting temperature T * , which depends on ψ w0 .It comes as a consequence that the ice fraction is the difference between v and θ w : It results that, under freezing conditions (T < T * ), θ w and θ i are function of ψ w0 , which dictates the saturation degree, and T , that dictates the freezing degree.Equation ( 24), usually called "freezing-point depression equation", relates the maximum unfrozen water content allowed at a given temperature in a soil.Figure 2 reports the freezing-point depression equation for pure water and the different soil types according to the Van Genuchten parameters given in Table 1.Equations ( 21) and ( 17) represent the closure relations sought for the differential equations of mass conservation (Eq.6) and energy conservation (Eq.8).

The decoupled solution: splitting method
The final system of equations is given by the equations of mass conservation (Eq.6) and energy conservation (Eq.8): The previous system is a function of T and ψ w0 and can be solved by the splitting method, as explained in Appendix B.
In the first half time step, the Richards' equation is solved and the internal energy is updated with only the advection contribution.In the second half, no water flux is allowed, which makes the volume a closed system, and the internal energy is updated with the conduction flux in order to find the new temperature and the new combination of water and ice contents.

Step 1: water and advection flux
Let us indicate with the superscript "n" the quantities at the time step n, with "n + 1" the quantities at the time step n+1: then t n+1 = t n + t ( t being the integration interval), and with "n + 1/2" the quantities at the end of the first step (temporary quantities).In the first step of the splitting method, www.the-cryosphere.net/5/469/2011/The Cryosphere, 5, 469-484, 2011 the temperature T is kept constant and phase change is not allowed, i.e. the ice content remains constant.The resulting equations to be solved are Eqs.( B3) and (B4), whose solution are represented by the flux J w , the water pressure ψ n+1/2 = ψ n + ψ, which results in ψ n+1 w0 = ψ n w0 + ψ (see Eq. 21), and eventually in the new total water content Eq. 22).The variation of liquid water content θ fl w is subject to the following limitation given by volume conservation: The internal energy of the System is changed by the quantity:

Step 2: energy flux and phase change
In the second step of the splitting method, water flux is not allowed and v = θ n w + θ n i + θ fl w is kept constant; for mass conservation this means that the mass of water subject to phase change is equal to the mass of ice subject to phase change.The resulting equations to be solved are Eqs.( B5) and (B6), whose solution are represented by the new temperature T n+1 and the new combination of ice and liquid water content.The variation of liquid water due to phase change θ ph w is subject to the following limitation due to the volume conservation: The internal energy of the system is eventually subject to the following variation: 6 The numerical implementation Equations (B3) and (B5) share the common characteristics of the nonlinear diffusion-advection equations.Let us now propose a numerical scheme to solve this type of equations in one dimension, based on the Newton method.The notation is based on Eq. ( B5), but the same scheme applies also to Eq. (B3).
Let us consider a 1-D domain composed by a soil column with uniform area (m 2 ) and let us divide the column into N cells.Let us indicate with the subscript "l" the cell number l (1 ..., N): according to this notation, z l represents the position of the cell l, z l+1 the position of the cell l + 1, located below the cell l, and z l+1/2 identifies the lower border of the cell l, that coincides with the upper border of the cell l + 1.It results that: z l+1/2 = z l + z l /2 where z l is the depth of the cell l.
In one dimension the lateral energy fluxes are neglected, therefore the energy conservation equation becomes: where, for simplicity of notation, S is intended as S en .Integrating the above equation in the cell l of volume V c from the upper border (z l−1/2 ) to the lower border (z l+1/2 ) one gets: where the quantities U * (J) and S * (W) are to be considered as the integral in the cell: Applying the divergence theorem yields: where: After discretization, the equation is still exact and no approximation of any of its quantities has been made.The time derivative in Eq. ( 31) can be estimated, by using a finite difference scheme between t n and t n+1 .This yields: where energy is integrated over the cell volumes.The discretized equation can be written as: R l (T ) l = 1..., N (J ) is a component of an array of N functions, said residuals; f l is the sum of the fluxes and the source term at time step n + 1: ) is a non linear function of the temperature in the l-th volume at time t n+1 , and finding its root (for any l) is equivalent to find the solution of Eq. ( 30).This problem can The Cryosphere, 5, 469-484, 2011 www.the-cryosphere.net/5/469/2011/be solved iteratively through the Newton method (e.g.Kelley, 2003), which consists of approximating the non linear functions as: where "l,k" represent a layer index, "m" the Newton iteration number, "n" the time step, (J n+1,m R ) l,k is the Jacobian matrix of R n,m k .The approximate solution is obtained by solving iteratively the linear system Eq.( 39), until the following condition is met: where (W) is the tolerance on the energy balance that can be set by the user.Eventually the new temperature becomes: with: The Jacobian is defined as: Considering a constant cell area in the same soil column, the above matrix becomes tridiagonal and is composed by: where all the symbols are defined in Table C1.At each of the Newton iterations, say m + 1, finding the approximate roots means solving the linear system Eq.( 39).In this 1-D case the matrix is tridiagonal and the linear system may be solved by the Thomas algorithm (Conte and De Boor, 1980).The details of the numerical procedure may be found in Dall'Amico (2010).

The boundary conditions
The boundary conditions may be of Neumann or Dirichlet type.The Neumann boundary conditions are represented by the water and energy fluxes at the boundaries.For the mass conservation, the flux is represented by rain or snow melting (mm s −1 ) at the top, and the water leakage (mm s −1 ) at the bottom.For the energy conservation, at the top it is represented by the soil heat flux (W m −2 ), as a result of the surface energy balance; at the bottom, it represents the geothermal heat flux (W m −2 ).The Dirichlet boundary conditions, both at the top and at the bottom, are given by values of soil water pressure head for the mass conservation and temperature for the energy conservation.

The non-linearity in the energy conservation equation
The derivative dU * dT n+1,m l in Eq. ( 44) deserves special attention.From the definition (Eq.32) one gets: Deriving Eq. ( 4) with respect to T and applying the first Equation of the System (Eq.B6), after some intermediate steps, one gets: The variation of liquid water content in temperature may be calculated through the derivative chain rule: The term ∂θ w /∂ψ := C H (m −1 ) is defined as the hydraulic capacity and is the slope of the soil water retention curve.
The second term ∂ψ/∂T represents the slope of the pressuretemperature relation (Eq.21) derived by the Clausius-Clapeyron equation.Eventually inserting Eq. ( 47) into Eq.( 45) yields: where C a (J m −3 K −1 ) is the so-called apparent heat capacity (Williams and Smith, 1989): and is the sum of the sensible heat transmitted to the soil matrix and the latent released by phase change.The apparent heat capacity formulation is one of the approaches usually used to deal with phase change, as it has the advantage of relating the latent heat term of the equation to temperature, without the need to split the pure conduction and phase www.the-cryosphere.net/5/469/2011/The Cryosphere, 5, 469-484, 2011  (Shoop and Bigl, 1997).However, with this approach the formulation of the internal energy becomes highly non-linear near T * , where it increases by several orders of magnitude, often inducing numerical oscillations and convergence problems (cf.Hansson et al., 2004).The following section proposes an improved scheme to deal with the high non-linearities typical of phase change problems.

Globally convergent iteration
The Newton method transforms the initial non-linear problem into a sequence of linear problems and, according to a tolerance , preserves the conserved variable.However, it only works well if the initial guess is close enough to the true solution; typically, a region that is well-behaved is first located with some other method and the Newton method is then used to refine the solution (the root finding).It may happen that the residual array of function is not convex or concave close to the solution (always decreasing or increasing with increasing or decreasing derivative, respectively), and the method may pass from T 0 to T 1 without converging, as shown in Fig. 3.In the energy equation solution this happens during the phase transition, when temperature passes from positive to negative values or vice-versa.At positive temperatures the heat capacity is ≈2 MJ m −3 K −1 and at -0.1 • C it assumes more or less the same value.All the latent term of the equation, in fact, is comprised in very small temperature intervals, where the peak of the apparent heat capacity is positioned, and may increase by three orders of magnitude depending on the Van Genuchten parameters.Hansson et al. (2004) recommend, in order to converge, to set the value of the heat capacity to its maximum value when passing from positive to negative temperature (let us hereafter indicate this method as Newton C-max).However, for very steep freezing characteristic curves, Newton C-max is not converging and therefore precludes model verification by comparison with the analytical solution (see Fig. 4 and section 7.1 for the test explanation).
A considerable improvement was obtained using the socalled globally convergent Newton scheme, hereafter Newton Global (e.g.Kelley 2003).This is based on the fact that the direction of the tangent given by the Newton scheme is always good in the sense that it points to the direction of decreasing residuals.Yet the final point may be too far from the solution causing the scheme to oscillate.In order to avoid this, the globally convergent Newton scheme tests the residual: This test implies that, far from the solution, the increment is multiplied by a reduction factor δ with 0 ≤ δ ≤ 1.If δ = 1 the scheme coincides with the normal Newton-Raphson scheme.This method, together with the maximum heat capacity imposition, allows the scheme to converge.This scheme is also applied by Tomita (2009) to solve surface energy balance equation, when the surface temperature shows oscillations caused by the exclusion or poor consideration of the surface temperature dependence of the turbulent transfer coefficient at the surface.

Model testing
This section deals with the verification of the model and the numerical algorithm.The first test regards the application of the model to a special case of freezing, know as the Stefan problem, of which an analytical solution exists for special cases.Neglecting water movement and energy advection, the objective of this comparison is to test the proper implementation and accuracy of the numerical algorithm applied to Eq. (B5).The second test is a comparison against experimental data provided by Hansson et al. (2004), where a soil column, under controlled boundary conditions, is forced to freezing.The objective is to test the splitting method for solving both mass and energy conservation (Eq.25) in a case where both water and heat fluxes are allowed.In a final experiment, the solution of infiltration into frozen soil is demonstrated.

Phase change without water flux: analytical solution
The model and the numerical algorithm are compared against the analytical solution of unilateral freezing of a change terms in the equation (Shoop and Bigl, 1997).However, with this the formulation of the internal energy becomes highly non-linear near T * , where it increases by several orders of magnitude, often inducing numerical oscillations and convergence problems (cf.Hansson et al., 2004).The following section proposes an improved scheme to deal with the high non-linearities typical of phase change problems.

Globally convergent iteration
The Newton method transforms the initial non-linear problem into a sequence of linear problems and, according to a tolerance , preserves the conserved variable.However, it only works well if the initial guess is close enough to the true solution; typically, a region that is well-behaved is first located with some other method and the Newton method is then used to refine the solution (the root finding).It may happen that the residual array of function is not convex or concave close to the solution (always decreasing or increasing with increasing or decreasing derivative, respectively), and the method may pass from T 0 to T 1 without converging, as shown in Fig. 3.In the energy equation solution this happens during the phase transition, when temperature passes from positive to negative values or vice-versa.At positive temperatures the heat capacity is ≈2 MJ m −3 K −1 and at −0.1 • C it assumes more or less the same value.All the latent term of the equation, in fact, is comprised in very small temperature intervals, where the peak of the apparent heat capacity is positioned, and may increase by three orders of magnitude depending on the Van Genuchten parameters.Hansson et al. (2004) recommend, in order to converge, to set the value of the heat capacity to its maximum value when passing from positive to negative temperature (let us hereafter indicate this method as Newton C-max).However, for very steep freezing characteristic curves, Newton C-max is not converging and therefore precludes model verification by comparison with the analytical solution (see Fig. 4 and Sect.7.1 for the test explanation).
A considerable improvement was obtained using the socalled globally convergent Newton scheme, hereafter Newton Global (e.g.Kelley 2003).This is based on the fact that the direction of the tangent given by the Newton scheme is always good in the sense that it points to the direction of decreasing residuals.Yet the final point may be too far from the solution causing the scheme to oscillate.In order to avoid this, the globally convergent Newton scheme tests the residual: if || > ||R(T ) n+1,m || (50) then This test implies that, far from the solution, the increment is multiplied by a reduction factor δ with 0 ≤ δ ≤ 1.If δ = 1 the scheme coincides with the normal Newton-Raphson scheme.This method, together with the maximum heat capacity imposition, allows the scheme to converge.This scheme is also applied by Tomita (2009) to solve surface energy balance equation, when the surface temperature shows oscillations caused by the exclusion or poor consideration of the surface temperature dependence of the turbulent transfer coefficient at the surface.

Model testing
This section deals with the verification of the model and the numerical algorithm.The first test regards the application of the model to a special case of freezing, know as the Stefan problem, of which an analytical solution exists for special cases.Neglecting water movement and energy advection, the objective of this comparison is to test the proper implementation and accuracy of the numerical algorithm applied to Eq. (B5).The second test is a comparison against experimental data provided by Hansson et al. (2004), where a soil column, under controlled boundary conditions, is forced to freezing.The objective is to test the splitting method for solving both mass and energy conservation (Eq.25) in a case where both water and heat fluxes are allowed.In a final experiment, the solution of infiltration into frozen soil is demonstrated.

Phase change without water flux: analytical solution
The model and the numerical algorithm are compared against the analytical solution of unilateral freezing of a semi-infinite region given by Neumann.The features of this problem are the existence of a moving interface between the two phases, in correspondence of which heat is liberated or absorbed, and the discontinuity on the thermal properties of the two phases (Carslaw and Jaeger, 1959).The assumptions are: (  Fig. 6.Cumulative error associated with the the globally convergent Newton method.Solid line: cumulative error (J), dotted line: cumulative error (%) as the ratio between the error and the total energy of the soil in the time step.was set to 1E-8.
the two phases, in correspondence of which heat is liberated , respectively.The initial conditions are: T i (t = 0,z > 0) = +2 • C and a substance completely unfrozen: θ w (t = 0,z) = 1 and θ i (t = 0,z) = 0.The boundary conditions of Dirichlet type: for the top and bottom boundary, respectively.As Nakano and Brown (1971) did for the case of an initially frozen soil, in Appendix C we reported the complete derivation of the solution both for freezing and thawing processes.
As the analytical solution considers the freezing of pure water, in the numerical scheme we have considered a soil with porosity θ s = 1 characterized by a very steep soil water retention curve with no residual water content, approaching a step function (Table 1).
The domain is composed of 500 cells characterized by a uniform depth z = 10 mm; the integration time t = 10 s.We have already shown in Fig. 4a the results of the test without the globally convergent method.Figure 5 shows the comparison between the numerical and the analytical solutions of the soil temperature profile using the globally convergent Newton's method.The analytical solution is represented by the dotted line and the simulation according the numerical model by the solid line.The results are much improved with the globally convergent method, as the simulated temperature follows the analytical solution very well.The temperature evolution shows a change in the slope that coincides with the separation point between the upper frozen and the lower thawed part.Figure 4b reports the comparison on the time axis (days) at different depths (m).The numerical simulation result shows oscillations, which begin at the time of phase change and then dampen with time.In the numerical solution the temperature starts decreasing only when all the water in the grid cell has been frozen.Furthermore, T l is influenced by the phase change of T l+1 by the release of latent heat and thus the temperature oscillation continues also in the frozen state.Therefore, oscillation amplitude is both linked to the grid size and to the time: increasing the grid size, the oscillation amplitude increases, as the mass of water to freeze increases before the temperature may decrease.The oscillation amplitude dampens with time as the freezing front moves away from z l ; it may be reduced but not eliminated, as it is embedded with the fixed-grid Eulerian method, where the freezing front may move in a discrete way and not in a continuum as in the reality.In order to test the energy conservation capabilities of the algorithm, the error, defined as the difference between the analytical and the simulated solution, was calculated at each time step as the p-norm (p = 1) of all the components and was cumulatively summed for the duration of the simulation.The tolerance on the energy www.the-cryosphere.net/5/469/2011/The Cryosphere, 5, 469-484, 2011 Fig. 6.Cumulative error associated with the the globally convergent Newton method.Solid line: cumulative error (J), dotted line: cumulative error (%) as the ratio between the error and the total energy of the soil in the time step.was set to 1E-8.
the two phases, in correspondence of which heat is liberated balance was set to 1 × 10 −8 .Figure 6 shows the cumulated error in J (solid line) and in percentage as the ratio between the error and the total energy of the soil in the time step.With set to 1×10 −8 , after 75 days of simulation, the error in percentage remains very low (< 1 × 10 −10 ), suggesting a good energy conservation capability of the algorithm.

Coupled water and energy flux: experimental data
In order to test the splitting time method for solving the coupled water and energy conservation, as done by Daanen et al. (2007), the model was tested against the experiment of Hansson et al. (2004).The soil considered represents a Kanagawa sandy loam, with the following parameters: and saturated hydraulic conductivity K H (sat)= 0.0032 mm s −1 .The column was considered initially unfrozen, with a uniform total water content v = 0.33 which, given the parameters of the soil retention curve, corresponds to ψ w0 = −2466.75mm, and initial temperature T = 6.7 • C uniform.The boundary conditions are of Neumann type: for the energy balance, at the top a flux F = −28 • (T 1 + 6) was considered, where T 1 is the temperature of the first layer, and a zero flux condition at the bottom.For the mass balance equation, a zero flux at both top and bottom boundaries were used.
Figure 7 shows the comparison of the profile of the total water content v .Starting from a thawed condition and a uniform water content v = 0.33, the liquid water content Fig. 6.Cumulative error associated with the the globally convergent Newton method.Solid line: cumulative error (J), dotted line: cumulative error (%) as the ratio between the error and the total energy of the soil in the time step.was set to 1E-8.
the two phases, in correspondence of which heat is liberated Fig. 6.Cumulative error associated with the the globally convergent Newton method.Solid line: cumulative error (J), dotted line: cumulative error (%) as the ratio between the error and the total energy of the soil in the time step.was set to 1 × 10 −8 .
decreases from above due to the increase of ice content.It is visible that the freezing of the soil sucks water from below.The increase in total water content reveals the position of the freezing front: after 12 h it is located about 40 mm from the soil surface, after 24 h at 80 mm and finally after 50 h at 140 mm.Similar to Hansson et al. (2004), the results were improved by multiplying the hydraulic conductivity by an impedance factor, as described is Sect.3.1.It was found that the value of ω that best resembles the results is 7.
Figure 8 shows the cumulative number of iterations required by the Newton C-max and the Newton global schemes to converge.It is clear that the number of iterations of the new method is much lower than the previous, indicating that this method provides improved performance on the total simulation time.

Infiltration into frozen soil
The coupled mass and energy conservation algorithm was finally tested with simulated rain (infiltration) during the thawing of a frozen soil.The soil geometry is a 20 cm depth column discretized in 800 layers of 0.25 mm depth; the top boundary condition for the energy equation is of Neumann type, with 10 W m −2 constant incoming flux (no daily cycle) and the bottom boundary condition is a zero energy flux.The initial conditions on the temperature are: T i (t = 0,z > 0) = −10 • C; the soil is considered initially unsaturated, with water pressure head ψ w0 given by a hydrostatic profile based on a water table at 5 m depth.This profile The Cryosphere, 5, 469-484, 2011 www.the-cryosphere.net/5/469/2011/after 12 hours q q q q q q q q q q q q q q q q q q q q Sim q Meas water content after 24 hours q q q q q q q q q q q q q q q q q q q q Sim q Meas water content after 50 hours q q q q q q q q q q q q q q q q q q q q Sim q Meas Fig. 7. Comparison between the numerical (solid line) and the experimental results (points) obtained by Hansson et al. (2004) of the total water (liquid plus ice) after 12 (left), 24 (center) and 50 hours (right).

Infiltration into frozen soil
The coupled mass and energy conservation algorithm was finally tested with simulated rain (infiltration) during the thawing of a frozen soil.The soil geometry is a 20 cm depth column discretized in 800 layers of 0.25 mm depth; the top boundary condition for the energy equation is of Neumann type, with 10 W m −2 constant incoming flux (no daily cycle) and the bottom boundary condition is a zero energy flux.The initial conditions on the temperature are: T i (t = 0,z > 0) = −10 • C; the soil is considered initially unsaturated, with water pressure head ψ w0 given by a hydrostatic profile based on a water table at 5 m depth.This profile corresponds, according to Eq. ( 23) and ( 24), to the water and ice contents at each level.The soil texture and thermal parameters are as in the experiment described in Section 7.2, whereas the saturated hydraulic conductivity is 0.3 mm s −1 .As far as the boundary condition on the mass is concerned, the bottom boundary is characterized by a no flux condition, whereas the top boundary condition varies along two simulations: zero flux (without rain) and constant 10 mm h −1 flux, resembling a constant precipitation (with rain).As can be seen in Fig. 9, in the first 6-7 days the behaviors with rain (solid line) and without precipitation (dotted line) are almost equal, because hydraulic conductivity is so low due to the fast ice-saturation of the first layer during cold conditions neat the beginning of the experiment.Then, when some ice is melted, hydraulic after 12 hours q q q q q q q q q q q q q q q q q q q q Sim q Meas water content after 24 hours q q q q q q q q q q q q q q q q q q q q Sim q Meas water content after 50 hours q q q q q q q q q q q q q q q q q q q q Sim q Meas Fig. 7. Comparison between the numerical (solid line) and the experimental results (points) obtained by Hansson et al. (2004) of the total water (liquid plus ice) after 12 (left), 24 (center) and 50 hours (right).

Infiltration into frozen soil
The coupled mass and energy conservation algorithm was finally tested with simulated rain (infiltration) during the thawing of a frozen soil.The soil geometry is a 20 cm depth column discretized in 800 layers of 0.25 mm depth; the top boundary condition for the energy equation is of Neumann type, with 10 W m −2 constant incoming flux (no daily cycle) and the bottom boundary condition is a zero energy flux.The initial conditions on the temperature are: T i (t = 0,z > 0) = −10 • C; the soil is considered initially unsaturated, with water pressure head ψ w0 given by a hydrostatic profile based on a water table at 5 m depth.This profile corresponds, according to Eq. ( 23) and ( 24), to the water and ice contents at each level.The soil texture and thermal parameters are as in the experiment described in Section 7.2, whereas the saturated hydraulic conductivity is 0.3 mm s −1 .As far as the boundary condition on the mass is concerned, the bottom boundary is characterized by a no flux condition, whereas the top boundary condition varies along two simulations: zero flux (without rain) and constant 10 mm h −1 flux, resembling a constant precipitation (with rain).As can be seen in Fig. 9, in the first 6-7 days the behaviors with rain (solid line) and without precipitation (dotted line) are almost equal, because hydraulic conductivity is so low due to the fast ice-saturation of the first layer during cold conditions neat the beginning of the experiment.Then, when some ice is melted, hydraulic corresponds, according to Eqs. ( 23) and ( 24), to the water and ice contents at each level.The soil texture and thermal parameters are as in the experiment described in Sect.7.2, whereas the saturated hydraulic conductivity is 0.3 mm s −1 .As far as the boundary condition on the mass is concerned, the bottom boundary is characterized by a no flux condition, whereas the top boundary condition varies along two simulations: zero flux (without rain) and constant 10 mm h −1 flux, resembling a constant precipitation (with rain).As can be seen in Fig. 9, in the first 6-7 days the behaviors with rain (solid line) and without precipitation (dotted line) are almost equal, because hydraulic conductivity is so low due to the fast ice-saturation of the first layer during cold conditions neat the beginning of the experiment.Then, when some ice is melted, hydraulic conductivity increases so that some water can infiltrate.At this point some incoming water may freeze because the soil is still cold, the ice content is increased (not shown here) and the zero-curtain effect is prolonged.As a result, infiltrating water provides energy (latent), so that temperature rises above 0 • C earlier than in the case of without rain.This earlier complete thaw is more evident at greater depths: in the case without rain, one has complete thawing at ≈ −0.1 • C (change of slope of dashed curve) because the soil remains unsaturated and is characterized by T * < T f , while in the case with rain complete thawing occurs at 0  conductivity increases so that some water can infiltrate.At this point some incoming water may freeze because the soil is still cold, the ice content is increased (not shown here) and the zero-curtain effect is prolonged.As a result, infiltrating water provides energy (latent), so that temperature rises above 0 • C earlier than in the case of without rain.This earlier complete thaw is more evident at greater depths: in the case without rain, one has complete thawing at ≈ −0.1 • C (change of slope of dashed curve) because the soil remains unsaturated and is characterized by T * < T f , while in the case with rain complete thawing occurs at 0 • C because there is saturation.When soil is thawed, in the with rain case soil temperature rises more slowly due to the lower thermal diffusivity of the soil.It is interesting to notice that water only partially infiltrates into the frozen soil (see that dotted line is equal to the solid line) because, as the soil is very cold, it freezes and so the first layers become saturated with ice thus forming an impervious barrier for the rain.The simulation results are based on the assumption of infiltration in thermodynamical equilibrium.

Conclusions
In this paper we have proposed a new method for robustly solving energy and mass balance equations in simulations dealing with phase change under in variably-saturated soil.
The differences with respect to previous methods are: (1) a very robust solution based on the globally convergent New-ton scheme is used; (2) the notation of the energy equation is based on the internal energy, which makes it possible to generalize the problem to a diffusion-advection equation, similar to the mass balance equation (Richard's equation), and thus the numerical method used in the energy equation can be further used to the mass balance equation; (3) the soil freezing curve implies that the ice content depends not only on temperature but also on the total water content, making this scheme useable in non-saturated conditions; and (4) the detailed explanation of the splitting method to decouple the system of equations is included.The test against analytical solutions shows good agreement, improving previous results obtained with methods based just on the maximum heat capacity formulation.The model was applied to simulate temperature-driven water flow in freezing soil, and the results were compared to the experimental findings of Hansson et al. (2004) with good agreement.Similar to previous studies, a high sensitivity to the value of the impedance factor ω of the hydraulic conductivity was found.This algorithm can now be applied in more realistic configurations, with complex boundary conditions accounting for the soil-atmosphere energy exchange, and thus contribute to improve our understanding of the factors controlling the soil freezing and thawing processes in the Alpine and Arctic cryosphere.Endrizzi et al. (submitted), who included this algorithm in the open-source hydrological model GEOtop (Rigon et al., 2006), present a first application.
one obtains m which is defined as the total equivalent liquid water content present in V c .The relation between v and m is: The splitting method The system represented by Eq. ( 25) can be arranged as: where the superscript "fl" refers to the changes due to the "flux" of water, and the superscript "ph" refers to the changes due to the "phase change" of water.The above system is equivalent to: The last two Equations are only function of T and can be solved considering ψ w0 = cost: Eq. (B5) represents the energy conservation equation under a "no-flux" condition, that gives the new temperature and the amount of mass that undergoes phase change; Eq. (B6) can be used to update the ice and liquid water content with the mass of water that undergoes phase change.In the thawing case the frozen and thawed parts are below and above z = Z(t) respectively, and the analytical solution of v 1 and v 2 becomes: where ζ is given by the solution of:

Fig. 1 .
Fig. 1.Frozen soil constituents and schematization of the control volume V c

Fig. 1 .
Fig. 1.Frozen soil constituents and schematization of the control volume V c .
Fig.2.Freezing curve for pure water and various soil textures, according to the Van Genuchten parameters given in Table1.

Fig. 4 .
Fig. 4. Comparison between the analytical solution (dotted line) and the simulated numerical (solid line) at various depths (m).The numerical model in panel (A) uses Newton C-max the one in panel (B) uses Newton Global.Both have a grid spacing of 10 mm and 500 cells.Oscillations are present in panel (B) but not in panel (A) where no convergence is reached.

Fig. 5 .
Fig. 5. Comparison between the simulated numerical and the analytical solution.Soil profile temperature at different days.Grid size=10 mm, N=500 cells

Fig. 4 .
Fig. 4. Comparison between the analytical solution (dotted line) and the simulated numerical (solid line) at various depths (m).The numerical model in panel (A) uses Newton C-max the one in panel (B) uses Newton Global.Both have a grid spacing of 10 mm and 500 cells.Oscillations are present in panel (B) but not in panel (A) where no convergence is reached.

Fig. 4 .Fig. 5 .
Fig. 4. Comparison between the analytical solution (dotted line) and the simulated numerical (solid line) at various depths (m).The numerical model in panel (A) uses Newton C-max the one in panel (B) uses Newton Global.Both have a grid spacing of 10 mm and 500 cells.Oscillations are present in panel (B) but not in panel (A) where no convergence is reached.

Fig. 5 .
Fig. 5. Comparison between the simulated numerical and the analytical solution.Soil profile temperature at different days.Grid size = 10 mm, N = 500 cells.

Fig. 4 .Fig. 5 .
Fig. 4. Comparison between the analytical solution (dotted line) and the simulated numerical (solid line) at various depths (m).The numerical model in panel (A) uses Newton C-max the one in panel (B) uses Newton Global.Both have a grid spacing of 10 mm and 500 cells.Oscillations are present in panel (B) but not in panel (A) where no convergence is reached.

Fig. 8 .
Fig. 8. Cumulative number of iterations of the Newton C-max and the Newton global methods on the simulation test based on the experimental results obtained byHansson et al. (2004)

Fig. 8 .
Fig. 8. Cumulative number of iterations of the Newton C-max and the Newton global methods on the simulation test based on the experimental results obtained byHansson et al. (2004)

Fig. 8 .
Fig. 8. Cumulative number of iterations of the Newton C-max and the Newton global methods on the simulation test based on the experimental results obtained by Hansson et al. (2004).

Table 1 .
Porosity and Van Genuchten parameters for water and different soil types as visualized in Fig.2.

Table C1 .
Table of symbols used.