A method for solving heat transfer with phase change in ice or soil that allows for large time steps while guaranteeing energy conservation

The accurate simulation of heat transfer with phase change is a central problem in cryosphere studies. This is because the nonlinear behaviour of enthalpy as function of temperature can prevent thermal models of snow, ice and frozen soil from converging to the correct solution. Existing numerical techniques rely on increased temporal resolution in trying to keep corresponding errors withing acceptable bounds. Here, we propose an algorithm, originally applied to solve water flow in soils, as a method to solve these integration issues with guaranteed convergence and conservation of energy for any time step 5 size. We review common modeling approaches, focusing on the fixed-grid method and on frozen soil. Based on this, we develop a conservative formulation of the governing equation and outline problems of alternative formulations in discretized form. Then, we apply the nested Newton-Casulli-Zanolli (NCZ) algorithm to a one-dimensional finite-volume discretization of the energy-enthalpy formulation. 10 Model performance is demonstrated against the Neumann and Lunardini analytical solutions and by comparing results from numerical experiments with integration time steps of one hour, one day, and ten days. Using our formulation and the NCZ algorithm, the convergence of the solver is guaranteed for any time step size. With this approach, the integration time step can be chosen to match the time scale of the processes investigated. Copyright statement. © Author(s) 2020. This work is distributed under the Creative Commons Attribution 4.0 License. 15

in thick isotherm layers.
With phase change occurring over a range of temperatures, rather than at one specific temperature, front-tracking methods become computationally expensive  and conceptually ambiguous. This is the case in many industrial (Voller and Cross, 1981) and environmental problems. Additionally, front tracking is complicated because it requires either a deforming grid or a transformed coordinate system (Aschwanden and Blatter, 2009). By contrast, fixed-grid methods can accu-60 rately describing the thermodynamics of the problem without requiring additional complications in handling the computational domain. For these reasons, fixed-grid methods are generally preferable to front-tracking methods when simulating frozen soil.
Fixed-grid methods include the latent heat of fusion in their governing equation, avoiding the necessity to define a continuity condition across the moving boundary and related implementation problems. All contemporary fixed-grid methods we reviewed aim to solve the numerical integration using globally convergent algorithms. Three differing approaches for treating the latent 65 heat of fusion exists: the enthalpy method, using a source term, and using apparent heat capacity. As analytical expressions, these methods look the same because their governing equations can be obtained from each other by the chain rule of derivation.
As we will illustrate in the next section, problems can arise in the discrete domain where this rule is not always valid.
Here we present a numerical model of heat conduction with freezing and thawing in soils without water flow that guarantees exact energy conservation for any time step size and for a wide range of soil freezing characteristics. It is novel in using the 70 nested Newton-Casulli-Zanolli (NCZ) algorithm (Casulli and Zanolli, 2010) for solving the nonlinear system obtained from discretizing the governing equation, written in terms of the specific enthalpy, using an implicit finite volume scheme. The NCZ algorithm has previously been applied to solving water flow in soils and to our knowledge this is the first application for solving the heat equation. Long time steps are desirable in several applications including permafrost thaw, or surface components of climate models, and models dedicated to avalanche prediction. 75 The remainder of the paper is organized as follows. Section 2 reviews established approaches to study freezing and thawing phenomena in soils and points to relevant issues. Section 3 describes the new approach we propose. It details the discretization of the governing equation and the NCZ algorithm used to solve the resulting nonlinear numerical system. In Section 4, the new model is tested against analytical solutions and in Section 5, its performance is compared over a range of spatial and temporal resolutions. Section 6 summarises our findings and concludes this contribution. The Appendix B contains pseudocode to 80 facilitate the implementation of the method we describe in other models.

The governing equation and their numerical issues
The governing equation of the problem in the first of the three approaches is written in terms of both the total enthalpy and temperature ∂h(T ) ∂t = ∇ · [λ(T )∇T ] (1) 85 where h(T ) is the specific enthalpy, T is temperature, λ(T ) is the thermal conductivity, and t is the time.
In the approach relying on apparent heat capacity, the governing equation is where is the apparent heat capacity that is the sum of the actual heat capacity C T and a term representing the additional thermal capacity arising from phase change with the local derivative of the SFCC (Dall'Amico, 2010).
In the approach using a source term for latent heat, it is considered as a heat source and in this equation, there are two unknowns: the temperature, and the liquid fraction θ w appearing in the source term.

95
The specific enthalpy per unit mass is defined as where u is the specific internal energy, p is pressure, and v is the specific volume, the inverse of density. Assuming that the heat transfer occurs at constant pressure and volume the differential of the specific energy and of the specific enthalpy are equal (Appendix D). However, since the term enthalpy method is commonly used in the literature, we will refer to enthalpy instead 100 of internal energy.
The specific enthalpy of a control volume of soil V c can be calculated as the sum of the enthalpy of the soil particles, liquid water and ice (Dall'Amico et al., 2011): Defining a reference temperature T ref the above terms becomes where l f is the specific latent heat of fusion, ρ sp , ρ w and ρ i are the densities of the soil particles, water, and ice, c sp , c w , c i are the specif heat capacity of the soil particles, water, and ice, θ w (T ) is the unfrozen water content, and θ i (T ) is the ice content.
The liquid water content and the ice content are evaluated using SFCCs (Dall'Amico et al., 2011) which are dependent on 110 temperature and, in the general case, on temperature and water saturation. Usually the reference temperature, T ref , is set to 273.15 K, the melting temperature of pure water at standard atmospheric pressure. By using Eq. (7) the enthalpy Eq. (6) can be rewritten as 4 where C T = ρ sp c sp (1 − θ s ) + ρ w c w θ w (T ) + ρ i c i θ i (T ) is the bulk heat capacity of the soil volume V c .

115
SFCCs have an inflection point (Bao et al., 2016;Hansson et al., 2004) causing a sharp change in their derivative. This nonlinear behaviour gives rise to convergence problems during the solution of the system of equations resulting from the numerical approximation of the governing equation (Voller, 1990;Casulli and Zanolli, 2010). This is true for any method used such as finite differences (Westermann et al., 2016;Bao et al., 2016;Sergueev et al., 2003), finite elements (McKenzie et al., 2007), and finite volumes (Dall'Amico et al., 2011). As a consequence, the robustness (stability) of the numerics used is a 120 fundamental and important issue in frozen soil models.
There is a more subtle aspect in the integration though. Analytically, Eq. (1, 2, and 4) are equivalent because Eq. (2 and 4) are derived from Eq. (1) by applying the chain rule of derivative on the enthalpy under the general assumption that the enthalpy is a differentiable variable. However, this is not necessarily so in the discrete domain where the derivative chain rule is not always valid. This is a known issue when dealing with hyperbolic equations (Roe, 1981), but often overlooked when treating 125 the parabolic ones.
The apparent heat-capacity approach, Eq. (2) can suffer from large balance errors in the presence of high nonlinearities and strong gradients (Casulli and Zanolli, 2010). The key to deriving a conservative numerical method here concerns the discretization of the apparent heat capacity, and Nicolsky et al. (2007b) as well as Voller et al. (1990) discussed suitable techniques.

130
Referring to the work by Roe (Roe, 1981), it can be proven that, the discrete operator of C a has to satisfy the requirement ensuring preservation of the chain rule at the discrete level. In Eq. (9) i refers to the spatial discretization index and n to the time discretization index. Approximating the time derivative in Eq.
(2) using a backward Euler scheme we obtaiñ 135 and substituting the condition Eq. (9) This shows that solving Eq. (2) with a conservative numerical method, i.e. by making use of Eq. (9), is equivalent to solving the enthalpy formulation, Eq. (1). Roe's condition, however, is often not checked in numerical models.
The source-term approach presents problems analogous to those of the apparent heat-capacity formulation. Specifically,

140
Eq. (4) is derived from Eq. (2) by moving the latent heat term to the right-hand-side of the equation. Equation (4) can be solved numerically using an iterative procedure  or the Decoupled Energy Conservation Parametrization method (DECP) (Zhang et al., 2008). As pointed out by Voller et al. 1990, the numerical solutions based on an iterative procedure may suffer from non-convergence problem unless under-relaxation is wisely applied, and additionally, it necessary to guarantee that the liquid fraction is in the range (0, 1). With DECP, the energy equation is first solved without latent heat.

145
Then, soil temperature and the liquid and solid fractions are readjusted to ensure energy conservation during phase change.
This method is mainly used in land-surface models (LSMs) (Dai et al., 2003;Foley et al., 1996;Verseghy, 1991). In this case, Nicolsky et al. (2007a) showed that it results in an artificial stretch of the phase change region, with consequent inaccuracies in the simulation of active-layer thickness. For comparison, Table 1 shows the diversity of formulations and solvers used in current models representing heat transfer and phase change in the cryosphere.

150
In summary, the governing equation can be written using three different approaches that are equivalent analytically, but not in their discrete formulation. Of the three, the enthalpy approach remains conservative, even when discretized, and should be preferred. An additional fundamental problem is the solution of the nonlinear system of equations. Current algorithms either require time step adaptation or may fail to converge, leading to unstable simulations and reduced computational efficiency (Casulli and Zanolli, 2010). Here we address this fundamental challenge by using the NCZ algorithm to solve the nonlinear 155 system of equations. Compared to other algorithms, it guarantees convergence of the solution for any integration time step.
When the time step is not constrained by numerical issues, it can be chosen to better match the time scale of the process under investigation.

A soil heat-transfer model using the NCZ algorithm
Frozen soil models are typically solved with time steps between seconds and hours. This may be motivated by the desire to 160 resolve diurnal phenomena near the ground surface, and also, this often arises from limitations of the numerical schemes used.
Many applications related to permafrost (Erum et al., 2019), on the other hand, only require the representation of seasonal and multi-annual variation, which can be accomplished using time steps of one or more days if permitted by the numerical schemes employed. In order to have a numerical scheme that does not suffer from time step restriction due to a stability condition, an implicit time integration is required. An implicit formulation includes the necessity of solving a nonlinear system 165 of equations and the algorithm used for this is of great importance. Existing linearization algorithms such as the Picard or the Newton, require a sufficiently accurate initial guess. As reported by Zanolli (2010, 2012) this can be obtained by using small time steps, often requiring an empirical criterion for time-step adaptation. Therefore, even if the numerical scheme does not require a time step restriction, one may still be required to solve the nonlinear system of equations. Because of this the convergence of Newton-type method is often problematic (Casulli and Zanolli, 2012). Instead, the NCZ algorithm 170 always converge in any situation, even under the use of large time steps and grid sizes when the initial guess is chosen as suggested Zanolli (2010, 2012). This is an important feature of the nonlinear solver, because an inexact solution of the nonlinear system is not conservative (Casulli and Zanolli, 2010).

Discretization of the domain
The domain is partitioned using an unstructured orthogonal grid (Casulli and Walters, 2000), consisting of a set of nonoverlap-175 ping convex volumes Ω i , i = 1, 2, .., N v , separated by M internal faces Γ j , j = 1, 2, ..., M . Let A j denote the nonzero jth face area. Within each control volume a centre must be identified in such a way that the segment joining the centres of two adjacent volumes and the face shared by the two volumes have a nonempty intersection, are orthogonal to each other and have distance δ j . Each control volume Ω i may have an arbitrary number of faces. Let F i denote the nonempty set of faces of the ith volume, with the exclusion of boundary faces. Moreover, let P(i, j) be the neighbour of volume i that shares face j with the ith control 180 volume so that 1 ≤ P(i, j) ≤ N v for all j ∈ F i . The discrete variables h i and T i are located at centre of each element Ω i . Using an implicit finite volume method, the discretization of Eq.
(1) reads as where ∆t is the time step size, is an optional source/sink term in volume, and h i (T ) is the ith enthalpy given by Eq. (12) can be written in matrix form as where T = {T i } is the vector of unknowns, h(T ) = h i (T i ) is a vectorial function representing the discrete enthalpy, A is the energy flux matrix, and b is the right-hand-side vector of Eq. (12), which is properly augmented by the known Dirichlet boundary condition when necessary. For a given initial condition T 0 i , at any time step n = 1, 2, . . . Eq. (12) constitutes a nonlinear system for T n+1 i , with the nonlinearity affecting only the diagonal of the system and being represented by the 195 enthalpy h i (T n+1 i ). This set of equations is a consistent and conservative discretization of Eq. (1). Therefore, regardless of the chosen spatial and temporal resolution, T n+1 i is a conservative approximation of the new temperature.

Solution of the nonlinear system
Difficulties in solving the nonlinear system of Eq. (16) arise from the non-monotonic behaviour of the derivative of the enthalpy, h(T ), with respect to temperature, and because in some parametrizations used for substances such as water, the derivative 200 of the enthalpy is not correctly defined. The NCZ algorithm used here was discovered by Casulli and Zanolli (2010) and overcomes these difficulties with a nested Newton algorithm, two subsequent Newton-type iterations. NCZ is based on Jordan decomposition (Chistyakov, 1997) of the enthalpy function, rewriting it as the difference of two monotonic functions on which the Newton algorithm can be applied separately in a nested iteration, as explained below. A mathematical proof of convergence exists for NCZ Casulli, 2008, 2009;Zanolli, 2010, 2012).

205
For each control volume the enthalpy function h i (T ) can be defined as where C a,i (T ) is is defined as C a,i (T ) = Ωi C a,i (T )dΩ. C a,i (T ) has to fulfil two requirements. The first one, C1, is that C a,i (T ) is defined for every T and that it is a nonnegative function function with bounded variations. The second one, C2, is that there exist a T * i such that C a,i (T ) is strictly positive and non decreasing in (0, T * i ) and nonincreasing in (T * i , +∞).

210
Since C a,i (T ) are nonnegative functions with bounded variations, they are almost everywhere differentiable, admit only discontinuities of the first kind, and can be expressed by using the Jordan decomposition ( Fig.1) as the difference of two nonnegative, nondecreasing, and bounded functions: By making use of Eq. (19) the algebraic system, Eq. (16), can be written as The NCZ algorithm requires first the linearization of h 2 and then the linearization of h 1 . By choosing the initial guess for the NCZ algorithm such that T 0 ≤ T * , a sequence of outer iterates {T k } is obtained from Eq. (20) by linearizing h 2 as follows: so that the outer iterates are solutions of the following nonlinear system: where The resulting kth (outer) residual is derived form Eq. 20 and reads as follow  For all k = 1, 2, . . ., by setting T k,0 = T k−1 , a sequence of inner iterates {T k,l } is derived from Eq. (22) by linearizing h 1 (T ) as follows so that the inner iterates are determined from the following linear systems where P k,l−1 = P (T k,l−1 ) and f k, form Eq. 22 and reads as follow The inner and the outer iterations are terminated when r k,l < ε, and r k < ε, respectively, with ε being a sufficiently small 240 prefixed tolerance.
The most commonly used constitutive SFCCs (McKenzie et al., 2007;Kozlowski, 2007;Dall'Amico et al., 2011;Sheshukov and Nieber, 2011;Watanabe et al., 2011), used to define the enthalpy of frozen soil, satisfy the assumptions C1 and C2. In particular, the NCZ approach can be successfully applied to SFCCs models derived from the combination of existing SWRC models and the Clapeyron equation that in general are difficult to implement in numerical models based on the apparent heat 245 capacity (Kurylyk and Watanabe, 2013). Functions describing the internal energy of other substances, for instance pure water (Andreas et al., 2005), satisfy the assumptions C1 and C2 and, therefore, the NCZ method can be successfully used to model phase change problems as it will be shown for the original Stefan problem, Sections (4.1).

Analytical Benchmarks
The numerical model is compared for the problem of a column of freezing water, i.e. the Stefan problem, with the analytical 250 solution presented by Neumann (cited in Kurylyk et al. (2014b)), and for the problem of a column of soil with the three-zone with the analytical solution presented by Lunardini (1988).

Neumann analytical solution
The Neumann analytical solution gives the solution of unilateral freezing of a semi-infinite domain for both the temperature profile and the position of the moving boundary. Kurylyk et al. (2014b) recommended the Neumann solution due to its ability 255 to represent differences between the thermal diffusivities of the thawed and frozen zones. Here we consider the freezing of pure water instead of soil since it is more numerically demanding. Consider a semi-infinite domain of pure water at temperature Fig. (2). At the surface a Dirichlet boundary condition is imposed T (z = 0, t) = T s , with T s < T m . As a consequence a freezing front ζ propagates downward separating the solid and the liquid phase. The governing equations are At the moving boundary ζ(t), the temperature is equal to the melting temperature of water, and the time evolution of ζ(t) is described by the third equation, the Stefan condition. This condition states that the difference of the heat fluxes at the interface 265 of the two substances is consumed for the phase change. The parameters used in the comparison are given in Table E1. The numerical model is able to simulate the freezing problem of water well as seen in Fig et al., 2011). It uses the Newton-Raphson algorithm to provide the right search direction and, in order to avoid overshooting, a reduction factor δ is used to find the new estimate. This represents an improvement over the Newton-Raphson method, but its ability to converge depends on the choice of the parameter δ and on the treatment of the apparent heat capacity (Hansson et al., 2004;Nicolsky et al., 2007b;Dall'Amico et al., 2011). As such, this algorithm does not guarantees to converge for any time     if ∆t = 10 s. With an hourly time step, however, the example with the globally convergent Newton method is not able to reproduce the position of the freezing front over longer periods of time. By contrast, the NCZ algorithm reproduces the analytical solution well using ∆t = 3600 s. The quality of the solution obtained with the globally convergent Newton algorithm depends not only on the time step duration but also on the definition of the parameter δ (Fig. 7). The additional necessity for an arbitrarily chosen parameter in the globally convergent Newton algorithm further underscores the robustness of the NCZ 295 algorithm, for which convergence only depends on the right definition of Eq. (18) and Eq. (19).

Lunardini analytical solution
Lunardini (1988)  We computed benchmark T1 proposed by the InterFrost project (InterFrost Project), parameters are given in Table ( T Freezing front Frozen Partially frozen Unfrozen Figure 8. Scheme showing the setting of Lunardini problem (Ruhaak et al., 2015). Initially the domain is unfrozen with T = T0. Because of Ts < 0 on the left boundary, a freezing front propagates from left to right. X1(t) and X(t) identify respectively the isotherm corresponding to Tm, and T f .

Numerical test
In the previous sections, we have demonstrated that the proposed method can reproduce the Neumann analytical solution, as 315 well as the Lunardini analytical solution, even when using larger time steps than other numerical models.
After comparing simulation results with analytical solutions, we now analyse the difference between solutions using hourly, daily, and 10-day time steps. The domain is a soil column of 20 m depth that is uniformly at T = −3°C, initially. The bottom  boundary condition is adiabatic and at the surface, we use a Dirichlet boundary condition. The original forcing has hourly resolution and for longer time steps, corresponding averages are computed. As temperature gradients and the influence of 320 phase change are usually greatest near the soil surface, the thickness ∆z is parameterized with an exponential function (Gubler et al., 2013) where ∆z min is the thickness of the first layer, b is the growth rate and i is the layer index, being one at the ground surface and increasing downward. The parameters used are reported in Table G1. All three simulations were spun-up for a period of 1400 325 years to reach a stable thermal regime. After spin-up, we performed a simulation of 100 years. there are no significant deviations in the results. The larger deviations occur when the zero-isotherm is shallow: at the beginning of the thawing season as well as the freezing one, Fig. (G1, G2). This can be attributed on one side to the diurnal cycles of surface boundary condition, and on the other side that using a larger time step we lose accuracy in capturing the timing of 330 thawing/freezing even if we use the same boundary condition.
With larger time steps, we lose some of the information of the boundary conditions and the accuracy of the numerical model decreases because it is first-order accurate in time. The overall performance relative to simulations with smaller time steps, however, is largely preserved. While the order of accuracy can be increased to second order in time using the Crank-Nicholson method, this would incur a time step restriction to guarantee the monotonicity of the solution. As this restriction is proportional 335 to the square of the space discretization, ∆z 2 , the Crank-Nicholson method would represent a severe constrain whenever high spatial resolution is required. shows the ground temperature envelope for the hourly simulation. The maximum envelop presents an 'elbow' that is due to the phase-change effects Fig. (G3). As can be seen in (b) and (d), close to the soil surface the hourly simulation presents larger 340 values for both the minimum and maximum temperature due the fact that the hourly boundary condition presents a greater amplitude that is smoothed computing the daily and 10-day average.
In the mean temperature profile, the 10-day simulation presents a larger deviation from the hourly simulation than the daily simulation. The large deviation can be explained with the interaction of the time-step size with the thermal offset effect (Fig.   G4). If the thermal conductivity of water is set equal to that of ice, the maximum difference between the three profiles is 345 reduced to 0.003°C with a maximum deviation of 0.003°C from the initial condition, that is also equal to the mean of the forcing boundary condition.
Regarding the spatial discretization Fig. (G5) reports a comparison of the zero-isotherm position obtained using an hourly time step, a daily time step, and a 10 day time step. The results are still in good agreement, but is it interesting to note that the zero-isotherm presents some steps, independently on the size of the time step, and some details are missed, such as the joining 350 of the downward and upward freezing fronts captured with the finer grid. These steps are caused by the greater thickness of the grid elements. Because temperature is computed in the middle of each control volume, more time is required to achieve complete phase change of water, resulting in slower variation of the zero-isotherm position.
These synthetic experiments demonstrate that spatial and temporal discretization can be chosen accordingly to the aim the study without any constrains due to the convergence and stability issues of the numerical scheme. shows the deviation of the position of the zero-isotherm after 100 years between the hourly and the 10-days simulation, and between the daily and the 10-days simulation.
Moreover, we checked the mean number of iterations required to solve the nonlinear system for a simulation lasting 1 year with a time step ∆t = 1 h and for different spatial discretizations. We compared the NCZ against the Newton-Raphson algorithm and the globally convergent Newton. Results are reported in Table G2.

Conclusions
We have presented a new model for simulating the ground thermal regime in the presence of freezing and thawing based on 360 the heat-transfer equation and the application of the NCZ algorithm. To our knowledge, this is the only method that guarantees convergence while also permitting large time steps. The numerical model was implemented and verified against the Neumann and Lunardini analytical solutions. In both cases, the results were in good agreement even with an hourly integration time step.
For the Neumann solution, we considered pure water instead of saturated soil since it is more numerically demanding, and no convergence problems were encountered despite choosing a narrow temperature range (0.0001°C) over which phase change 1 h 1 d 10 d Figure 12. (a) The minimum, mean, and maximum temperature profile for the hourly simulation. (b), (c), (d) show the comparison of the minimum, mean, and maximum temperature profile respectively for the three simulations: with an hourly air temperature boundary condition and ∆t = 1 h, with a daily air temperature boundary condition and ∆t = 1 day, with a ten day air temperature boundary condition and ∆t = 10 day. All three simulations last 100 years. The maximum difference of Tmean between the hourly, and daily simulation is of 0.04°C, while between the hourly, and ten-days simulation is of 0.3°C.
the user to choose both the space and time discretization without any restriction due to stability and convergence issues. As a consequence, this method is effective for simulating permafrost thaw, a phenomenon that occurs at depth, in response to 370 seasonal and multi-annual cycles, and often over tens, hundreds or even thousands of years. Furthermore, phenomena like hysteresis or the variation of solute concentration upon freezing (Clow, 2018) can be included in the numerical model if the enthalpy function (i.e. its parameters) does not change within the current time step of integration.
While we presented a finite volume method, the NCZ algorithm can be also used with finite difference and finite element method. Beyond applications to frozen soil, it can be used to study other geophysical phenomena that involve phase change of a 375 substance simply by changing the definition of the enthalpy function and the thermal conductivity function. Examples include, glacier dynamics (Aschwanden et al., 2012), snow pack evolution (Brun et al., 1992;Lehning et al., 1999), and magma bodies (De Lorenzo et al., 2006). This may be even further expanded to industrial problems involving phase change materials used in energy recovery systems (Mongibello et al., 2018;Nazzi Ehms et al., 2019) or casting problems of pure metals and alloys (Lewis and Ravindran, 2000).

codeavailability
The source code is written in Java using the object-oriented programming paradigm. It can be found at https://github.com/ geoframecomponents/FreeThaw1D (Tubini, c). The OMS3 project can be found at https://github.com/GEOframeOMSProjects/ OMS_FreeThaw1D (Tubini, b). FreeThaw1D is deployed as an open source code to work alone or inside the Object Modelling System version 3 framework (David et al., 2013). In the latter case it can be connected at run time with the many 385 20 other components developed along with the GEOframe system (Formetta et al., 2014;Bancheri, 2017) for providing hydrometeorological forcings and other fluxes, like the evapotranspiration. The simulations presented here can be found at http: //dx.doi.org/10.5281/zenodo.4017668 (Tubini, a). The code must be run inside the OMS3 Console or using the Dockerized version of OMS3. For setting up the environment please follow the steps described in the README file present it he Github repository https://github.com/GEOframeOMSProjects/OMS_FreeThaw1D and in the GEOframe pages at https://geoframe.blogspot. Heat transfer with phase change of water is a cross-cutting problem existing in many geophysical phenomena other than frozen soil. This includes, for example, the seasonal snow pack, glaciers, and ice-sheets. Our contribution does not seek to 400 present an improvement in the description of these problems and we ignore typical processes such as metamorphism and settling in seasonal snow or strain heating and deformation in glaciers and ice sheets. Nevertheless, corresponding models may benefit from the NCZ algorithm in the treatment of the nonlinearity arising from phase change and, furthermore, broadening our review to also include some snow and glacier models supports the generalisation of our findings.

405
The Community Land Model (CLM) is the LSM for the Community Earth System Model (Oleson et al., 2004). It includes a module to simulate the ground temperature considering freezing and thawing. The governing equation is written in the nonconservative form and does not include the latent heat term (Oleson et al., 2004) (Lawrence et al., 2019). The heat conduction equation is solved using a Crank-Nicholson method. The temperature profile is calculated adopting the DECP approach. This approach does not require to solve a nonlinear system, since the latent heat is treated in an explicit way, but Nicolsky et al.

410
(2007a) have pointed out that this two-step procedure can overestimate the region where the phase change occurs, resulting in inaccuracies in the simulation of active-layer thickness.

A2 CoupModel
The CoupModel (Jansson and Karlberg, 2011) is a one-dimensional numerical model to simulate the heat and water flow as well as carbon and nitrogen budgets in a soil-plant-atmosphere system (Hollesen et al., 2011). The governing equation for heat flow in the soil is defined using the apparent heat capacity, and solved with an explicit numerical method. This does not require to solve a non-linear system but sets a time step restriction to avoid numerical oscillation.

A3 CryoGrid
CryoGrid 2 simulates the ground thermal regime based on conductive heat transfer in the soil and in the snowpack (Westermann et al., 2013). The heat equation is written using the apparent heat capacity and solved using the method of lines (Westermann 420 et al., 2013). The resulting system of ordinary differential equations is solved numerically with the package CVODE of Sundials that implements a modified Newton method, and Inexact Newton method, or a fixed-point solver to linearize the algebraic system resulting from the discretization of the heat transfer equation. The convergence of the Newton-type methods can be problematic (Casulli and Zanolli, 2012).  change is treated with the front tracking method, which offers good accuracy for problems in which phase change occurs at a fixed temperature (Goodrich, 1982). This model does not use a SFCC, and instead, the soil is represented as homogeneous layers with distinct frozen and thawed thermal properties.

A7 Hydrus 1D
Hydrus 1D includes a module to simulate water flow and heat transport in frozen soil. The governing equation is written using the apparent heat capacity formulation and Picard iteration is used to linearize the algebraic nonlinear system. In their paper, Hansson et al. (2004) explain that during the Picard iteration the solution can easily oscillate whenever the temperature decrease below the melting temperature. To avoid these oscillation the temperature is reset to the critical value and iteration restarted.

450
Hydrus 1D adopts an empirical time-step adaptation criterion. It is worthwhile to notice that the modified Picard iteration was proposed by Celia et al. (1990) to solve the Richards equation -problem for which the NCZ algorithm was originally proposed (Casulli and Zanolli, 2010).

A8 MarsFlo
MarsFlo is a three-phase numerical model to simulate the heat transfer and water flow in partially frozen, partially saturated 455 porous media (Painter, 2011). The heat equation is written in the conservative form. The equation is solved using an implicit finite difference method, and the resulting nonlinear system is solved using a Newton-Raphson method. To overcome convergence and stability problems, three modification were introduced (Painter, 2011). The convergence of the Newton-Raphson method can be problematic (Casulli and Zanolli, 2012).

A10 Sergueev et al.
This is a two dimensional model and the governing equation is written in the enthaply form (Sergueev et al., 2003). This model 465 implements a fractional time step approach (Godunov splitting): each time step is divided into two steps and at each step, a different dimension is treated implicitly. The system of finite difference equations is non-linear and is solved with the Newton's method. As in GIPL-2.0, the time derivative of enthalpy is computed either using the difference derivative or the analytical derivative according with the gradient of the temperature field.

470
The heat equation is written using the apparent heat capacity. The equation are solved using a finite element solver, FlexPDE suite, both explicit and implicit in time. In case of implicit methods, the resulting non-linear system is solved using the Newton-Raphson method. The convergence of the Newton-Raphson method can be problematic (Casulli and Zanolli, 2012).

A12 SUTRA
SUTRA is an established USGS groundwater flow and coupled transport model (Voss and Provost, 2002). McKenzie et al. 475 (2007) and McKenzie and Voss (2013) have extended the model to simulate freezing and thawing processes in the soil. The heat equation is written using the apparent heat capacity formulation and nonlinearities are solved using Picard iteration. The convergence of the Newton-type method can be problematic (Casulli and Zanolli, 2012).

A13 Crocus
Crocus is a one-dimensional finite difference model that solves the mass and energy balance within the snowpack taking 480 into account metamorphism and settling. The first versions of Crocus (Brun et al., 1989(Brun et al., , 1992 were not enthalpy-based. The governing equation was written in terms of temperature and water content. It was solved by using the Crank-Nicholson method, and the phase change is treated by using the DECP approach (Brun et al., 1992). After the integration within SURFEX (Vionnet et al., 2012), Crocus uses the enthalpy formulation and the numerical scheme is fully implicit, based on the numerics of ISBA-ES (Boone and Etchevers, 2001). Similarly to the previous version, the heat balance equation is solved adopting the DECP it, they adopted an approach based on Picard iteration with variable time steps (Paniconi and Putti, 1994).

490
SNOWPACK (Lehning et al., 1999) solves the heat transfer and creep/settlement equations using a Lagrangian finite element method. The governing equation is written using the source/sink formulation and it is solved using the DECP approach (Bartelt and Lehning, 2002;Lehning et al., 1999). Regarding the water flow, SNOWPACK implements three different schemes: a simple bucket-type approach, an approximation of Richards equation, and the full Richards equation (Wever et al., 2014). The full Richards equation is solved using Picard iteration with variable time steps (Paniconi and Putti, 1994).

A15 ORCHIDEE
ORCHIDEE is terrestrial biosphere model and it is part of the IPSL-CM4 Earth system model developed by the Institute Pierre Simon Laplace (IPLS) (Krinner et al., 2005). In the version 1.9.6 the snow is described with a single layer of constant density (Wang et al., 2013). Because of the limitations of the this approach, Wang et al. introduce a three-layer snow model, ORCHIDEE-ES, largely inspired from ISBA-ES (Boone and Etchevers, 2001) to consider snow settling, water percolation, 500 and water refreezing. The governing equation is written in the non-conservative form and does not include the latent heat term.
The temperature profile is calculated adopting the DECP approach.

A16 JSBACH
JSBACH is the land surface model developed by the Max Plank Institute (Ekici et al., 2014). It is a component of the Earth System Model (MPI-ESM) that also include ECHAM6 for the atmosphere and MPI-OM for the ocean. JSBACH simulates 505 both the frozen soil and the snowpack. In both cases the heat conduction is assumed to be the dominant method of heat transfer.
The governing equation is written in the source term formulation and solved with the DECP approach (Ekici et al., 2014).

A17 Ice-sheet models
For glacier and ice-sheet models it is necessary to distinguish between cold and temperate ice. Following Aschwanden and Blatter (2005), "ice is treated as temperate if a change in heat content leads to a change in liquid water content alone, and is 510 considered cold if a change in heat content leads to a temperature change alone." This means that cold ice is always below the melting temperature and thus the phase change does not occur. As result, present-day ice sheet models can be classified into: 'cold-ice method' models and polythermal models.
'Cold-ice method' does not consider the phase change of ice. Because of this the heat capacity can be assumed to be constant and therefore the governing equation can be written in terms of only temperature. These models are easy to implement, but their 515 applicability is restricted since in general temperate zones can be present (Aschwanden and Blatter, 2009). In fact, since the phase change of ice is overlooked, locally, the 'cold-ice method' violates the energy conservation, overestimates the temperate region (Aschwanden and Blatter, 2009), and can not quantify the liquid water content that affects viscosity in temperate ice (Lliboutry and Duval, 1985).
By contrast, polythermal ice-sheet models consider the phase change of ice. Similar to freezing soil models, the polythermal 520 ice-sheet models can be classified in two groups on the base of the treatment of the phase change: front tracking method and enthalpy method (Nedjar, 2002). SICOPLOIS (Greve, 1997a, b;Greve and Blatter, 2016) is the only 'truly' polythermal ice sheet model. It employs the polythermal two-layer scheme (Greve, 1997b): the temperature field and the water content field are computed separately for the ice and temperate domain and a Stefan-type condition is applied at the cold-temperate surface (CTS). This model defines the CTS for both energy flux and mass flux. The drawback of this method relate to the 525 implementation and restriction on the geometry and topology of the CTS (Aschwanden et al., 2012). Aschwanden and Blatter (2009) presented an enthalpy gradient method. This is a fixed-grid method that differs from the enthalpy method commonly used for freezing soil in its definition of the energy flux. In the enthalpy method, the heat flux is expressed in terms of the temperature gradient, whereas in the enthalpy gradient method it is expressed in terms of enthalpy, assuming that the heat capacity is constant (Aschwanden and Blatter, 2009). The enthalpy approach combines the advantage of 530 solving one equation for the entire domain, cold-ice models, and the correct description of the thermodynamics of temperate ice (front tracking model). This model is implemented in COMSOL Multiphysics (Aschwanden and Blatter, 2009), where nonlinear problems are solved using either a Newton algorithm or a damped Newton algorithm. Also in this case the NCZ may represent a valid option to solve the nonlinear system. To the authors' knowledge, the enthalpy gradient method has not be used in freezing soil models.

535
at each time step the computational domain is explicitly divided in the cold and temperate regions, and the energy equation is solved adopting a combination of implicit and explicit methods (Hewitt and Schoof, 2016). It is worth to note that in the temperate region, temperature is set equal to the melting temperature of the ice. This limits the application of this model to simulate freezing soil, where temperature can be larger than the melting temperature of water.

Appendix B: Pseudocode
We present the pseudocode for a one-dimensional implementation of the NCZ algorithm. Since the matrix A in Eq. (16) is tri-diagonal we can efficiently compute only the non-zero diagonal: the upper diagonal, the main diagonal, and the lower diagonal. We use the generic expression Discretize the governing equation since here, we can choose to adopt either a finite volume method, as presented in this paper, a finite element method, or a finite difference method. Moreover, the matrix A is 545 symmetric and positive definite thus within the nested Newton algorithm the linearized algebraic system can be easily solved with the Thomas algorithm. Here, it is worthwhile to point out that when we move to the two-dimensional or three-dimensional problem, the linearized algebraic system cannot be solved with the Thomas algorithm as the matrix is no longer tri-diagonal.
In these cases, iterative schemes such as the Conjugate Gradient Method need to be used (Shewchuk, 1994).
Appendix C: Fully implicit formulation and the problem of the thermal conductivity 550 Using a fully implicit formulation, the discretization of Eq.
(1) reads as where T i is the temperature of the ith control volume, P(i, j) is the neighbour of volume i that shares face j with the ith control volume, δ j is the nonzero distance between the centers of two adjacent volumes which share the jth internal face, ∆t is the time step size, and is an optional source/sink term in volume, and h i (T ) is the ith enthalpy given by where the index m refers to the Picard iteration. By using the NCZ algorithm local and global energy conservation is enforced at each Picard iteration. However, convergence of the Picard iterations is not essential to get a conservative solution for Eq. (1), 565 but some iterations can be used to update the thermal conductivity to the n + 1th time level Casulli and Zanolli (2010).

27
Appendix D: Enthalpy and internal energy Following the work by Dall'Amico (2010), the internal energy in its canonical form, U c , can be written as where S is the entropy, V is the volume, and M the mass of the constituents. These are the independent variables and are called 570 extensive variables since they depend linearly on the mass of the substance. The first differential of Eq. (D1) is According to Callen (1985) it is possible to define The differential of the enthalpy is If we assume that the transformation occurs at constant pressure and volume then Eq. (D6) becomes and Eq. (D8) Hence, from Eq. (D9) and Eq. (D10) the differential of the internal energy and the differential of enthalpy are equal. Therefore the governing equation, Eq.
(1), can be equivalently written in term of either the specific enthalpy or the specific internal energy.

Appendix E: Neumann analytical derivation
In this section we report the derivation of the Neumann analytical. The enthalpy is defined as where the singularity of the enthalpy function at T = T m has been linearized with and is a parameter defining the temperature range over which the phase change of water occurs, Fig. (E1). In the following tests is set to be equal to 0.0001°C. The introduction of this linearization is necessary since the enthalpy function needs to be continuously differentiable according to assumption C1 in Section 3.2. It is worth to underline that the temperature range can be chosen sufficiently small in order to make this approximation negligible when compared to the physical behaviour of water, considering that: (a) The melting of water in temperate ice is known to actually occur progressively below 0°C along 600 grain boundaries (Langham, 1974;Nye and Frank, 1973). (b) Freezing often occurs below the melting point when nucleation is relevant. (c) In porous media such as soil, ice melts across a range of temperatures due to the Gibbs-Thompson effect in pores and surface affects at the interfaces between ice and particles (Rempel et al., 2004;Watanabe and Mizoguchi, 2002).
Even though the internal energy function is very steep, the code used does not suffer of convergence problem with a time step of 3600 s. The thermal conductivity is defined as: Defining the following constant: the moving boundary function is where the coefficient γ can be found solving the following equation Finally the analytical solution for problem with Dirichlet boundary condition for the thawed and frozen zones are:

Appendix F: Lunardini analytical derivation
The solution fo the Lunardini problem (i.e. the Lunardini solution) as described by McKenzie (McKenzie et al., 2007) is given by the following set of equations: where T 1 , T 2 , and T 3 are the temperatures at distance, x, from the temperature boundary for the frozen, partially frozen, and unfrozen zone respectively; erf and erf c are the error function, and the complementary error function respectively; T 0 , T m , T f , and T s are the temperatures of the initial condition; the solidus, the liquidus, and the boundary temperature, respectively; α 1 and α 3 are the thermal diffusivities for the frozen, and unfrozen zone respectively, defined as λ 1 /C 1 and λ 3 /C 3 where C 1 and C 3 are the volumetric bulk-heat capacities of the frozen and unfrozen zones. The thermal diffusivity of the partially frozen 630 zone is assumed to be constant across the transition region, and the thermal diffusivity with latent heat term included, α 4 , is defined as: where γ d is the dry unit of soil solids, and ∆ξ = ξ 0 − ξ f where ξ 0 and ξ f are the ratio of unfrozen water to soil solids for the fully thawed and frozen conditions respectively. For a time, t, in the region 0 ≤ x ≤ X 1 (t) the temperature is T 1 and X 1 (t) is given by and from X 1 (t) ≤ x ≤ X(t) the temperature is T 2 , where X(t) is given by and for x ≥ X(t) the temperature is T 3 . The unknowns, ψ and γ, are solving the set of these two equations:  Appendix G: Numerical test 645 As explained in Section 5, comparing the position of the zero-isotherm after 100 years using three different time step, hourly, daily and 10-days, there are no significant deviation in the results. The larger deviation occur when the zero-isotherm is shallow: at the beginning of the thawing season as well as the freezing one.
At the begging of the thawing season, Fig. (G1), there is a time lag of about one month between the beginning of the thawing season for the hourly simulation and the 10-days simulation. This can be attributed to different surface temperature used to    660 Figure (G4) shows the temperature envelopes for the hourly, the daily, and the 10-days simulations setting the thermal conductivity of water equal to the thermal conductivity of ice, λ w = λ i . It is interesting to note that mean temperature, panel (c), is constant throughout the soil column. The mean temperature is very close to the initial temperature profile that is also equal to the mean surface boundary condition.
Figure (G5) shows a comparison of the the zero-isotherm position by using a coarser space discretization. Again, the three 665 simulation with the hourly time step, the daily time step, and the 10-days time step are still in good agreement. By using a coarser spatial discretization, the zero-isotherm, panel (c), presents some 'steps' that are not present when using a finer grid, Fig. (11). Moreover, the joining of the downward and upward freezing front is not captured neither by the hourly nor by the daily simulations.
For this numerical test we checked the mean number of iterations required to solve the nonlinear system with the NCZ 670 algorithm, the Newton-Raphson algorithm, and the globally convergent Newton algorithm. We performed a simulation lasting 1 year with a time step ∆t = 1 h and for different spatial discretizations. As can be seen in Table G2, neither the Newton- when temperature remained within ±0.1°C. This is the so-called zero-curtain effect, and it is due to latent heat of fusion that is continually released during the freezing of soil moisture. with an hourly surface temperature boundary condition and ∆t = 1 h, with a daily surface temperature boundary condition and ∆t = 1 day, with a ten day surface temperature boundary condition and ∆t = 10 day. All three simulations last 100 years. Because λw = λi the mean temperature, panel (c), is constant throughout the soil column and it is not possible to appreciate the thermal offset. The mean temperature is very close to the initial temperature profile, the maximum error is of 0.003°C.
Raphson nor the globally convergent Newton converge: they always reach the maximum number of iterations allowed with a consequent increase of the computational cost. shows the deviation of the position of the zero-isotherm after 100 years between the hourly and the 10-days simulation, and between the daily and the 10-days simulation. By using a coarser spatial discretization, the zero-isotherm presents some 'steps', panel (c), independently on the size of the time step. Another consequence of this is that the joining of the downward and upward freezing front is not captured neither by the hourly nor by the daily simulations.  (Gubler et al., 2013). Table G2. Summary of the mean number of iterations for the NCZ algorithm, the Newton-Raphson algorithm (N. R.), and the globally convergent Newton algorithm (g. c. N.). The simulation lasts 1 year with a time step ∆t = 1 h. We considered different spatial discretizations.
The tolerance ε = 10e − 11 has been rescaled with the water latent heat of fusion and the water density. The maximum number of iteration for each time step is 40. As can be seen the Newton-Rapshon and the globally convergent Newton does not converge so it always reaches the maximum number of iteration allowed.