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 non-linear 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 within 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 size.

We review common modelling 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.

Model performance is demonstrated against the Neumann and Lunardini analytical solutions and by comparing results from numerical experiments with integration time steps of 1 h, 1 d, and 10 d. 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 timescale of the processes investigated.

1 Introduction

Freezing and thawing of soils affect a wide range of biogeochemical and hydrological (Walvoord and Kurylyk, 2016; Schuur et al., 2015) processes and interact with engineered structures in cold regions. Correspondingly, the simulation of freezing and thawing soil is an important an well-researched topic (Streletskiy et al., 2019; Harris et al., 2009). Climate change brings additional urgency and new phenomena of interest to these studies. It is thus not a surprise that many models of freezing and thawing soil and ice exist, some of which are reviewed in Appendix A. Here, we propose a solution to a central challenge that these models have in common.

Published models can be categorized as empirical, analytical, or numerical (Riseborough et al., 2008). Empirical methods relate ground temperature or thawing–freezing depth (TFD) to simple topoclimatic factors (Zhang et al., 2008; Riseborough et al., 2008) and are relatively simple to apply. By contrast, analytical and numerical models are based on the conservation of mass and energy and can be divided in two broad groups (Tan et al., 2011). The first group focuses primarily on freezing and thawing, commonly known as the Stefan problem. The governing equation describes energy conservation with the heat flux modelled using the Fourier law. The second group considers the coupled problem of heat transfer and water flow in soils. In this case the energy–enthalpy conservation equation includes also the advective heat flux, and it is coupled with the mass conservation equation. For both groups, the latent heat transfer during phase change of water leads to problems related to convergence, conservation, and restrictions to discretization of space and time (Bao et al., 2016).

Historically (Hu and Argyropoulos, 1996; Vuik, 1993), the first attempts to solve the problem of heat conduction considering the phenomena of solidification and melting date back to the studies by Lame and Clapeyron in 1831, as well as the analytical solutions presented by Stefan around 1890 and Neumann in 1921. Later, other analytical solutions were proposed in order to overcome some simplifications that were too restrictive (Zhang et al., 2008; Riseborough et al., 2008; Walvoord and Kurylyk, 2016). These analytical solutions, however, are limited to one-dimensional problems and constrained in their initial and boundary conditions as well as the description of soil characteristics (Kurylyk et al., 2014a).

By contrast, numerical models can accommodate complex processes or configurations, including soil heterogeneities, complicated temperature boundary conditions, intermittent freeze–thaw, and temporally variable thermal properties. Accurately representing phase transitions, however, is a non-trivial task, and several different methods have been published. They can be broadly cast into two general groups: the so-called front-tracking methods and the fixed-grid methods (Voller et al., 1990). Even though this contribution is focused on modelling heat transfer in frozen soil or ice, the following review includes, and is relevant for, other fields of research that involve phase change.

Front-tracking methods are suitable whenever the two phases are divided by a spatially smooth and continuous front, and thus the state of the system can be conveniently described by the position of this interface (Voller et al., 1990). The moving front is tracked defining a continuity (“Stefan”) condition on the heat flux across it. For example, the one-dimensional model by Goodrich (1978, 1982) uses front tracking in modelling frozen soil, and the SICOPOLIS model (Greve, 1997a, b; Greve and Blatter, 2016) uses it to model polythermal ice sheets.

In frozen soil, however, a significant proportion of water can remain liquid at temperatures well below 0 ^∘C. This depression of the melting temperature is due to the presence of solutes (Bouyoucos, 1913; Bouyoucos and McCool, 1915; Bouyoucos, 1920, 1923), surface effects in the interaction between water and soil particles as well as water and ice (Anderson and Tice, 1972; Clow, 2018), and the Gibbs–Thomson effect (Rempel et al., 2004; Watanabe and Mizoguchi, 2002). To some degree, also polycrystalline ice has a temperature-dependent liquid water content (Langham, 1974). The gradual phase change over a range of temperatures in soils is commonly described with the soil freezing characteristic curve (SFCC) (Kurylyk and Watanabe, 2013). Moreover the presence of a partially frozen region is also common in ice and snow where liquid and solid phases coexist 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 (Voller et al., 1990) 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 accurately describe 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 heat of fusion exist: using 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 nested Newton–Casulli–Zanolli (NCZ) algorithm (Casulli and Zanolli, 2010) for solving the non-linear 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.

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 non-linear numerical system. In Sect. 4, the new model is tested against analytical solutions, and in Sect. 5, its performance is compared over a range of spatial and temporal resolutions. Section 6 summarizes our findings and concludes this contribution. Appendix B contains pseudocode to facilitate the implementation of the method we describe in other models.

2 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:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial h\left(T\right)}{\partial t}}=\mathrm{\nabla }\cdot \left[\mathit{\lambda }\right(T\left)\mathrm{\nabla }T\right],\end{array}$$

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

$$\begin{array}{}\text{(2)}& C_{\mathrm{a}}{\displaystyle \frac{\partial T}{\partial t}}=\mathrm{\nabla }\cdot \left[\mathit{\lambda }\right(T\left)\mathrm{\nabla }T\right],\end{array}$$

where

$$\begin{array}{}\text{(3)}& C_{\mathrm{a}}={\displaystyle \frac{\partial h}{\partial T}}=C_T+{\mathit{\rho }}_{\mathrm{w}}{l}_{\mathrm{f}}{\displaystyle \frac{\partial {\mathit{\theta }}_{\mathrm{w}}}{\partial T}}\end{array}$$

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:

$$\begin{array}{}\text{(4)}& C_T{\displaystyle \frac{\partial T}{\partial t}}=\mathrm{\nabla }\cdot \left[\mathit{\lambda }\right(T\left)\mathrm{\nabla }T\right]-{\mathit{\rho }}_{\mathrm{w}}{l}_{\mathrm{f}}{\displaystyle \frac{\partial {\mathit{\theta }}_{\mathrm{w}}}{\partial t}},\end{array}$$

and in this equation, there are two unknowns – the temperature and the liquid fraction θ_w appearing in the source term.

The specific enthalpy per unit mass is defined as

$$\begin{array}{}\text{(5)}& h=u+pv,\end{array}$$

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 that 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 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):

$$\begin{array}{}\text{(6)}& h=h_{\mathrm{sp}}+h_{\mathrm{w}}+h_{\mathrm{i}}.\end{array}$$

Defining a reference temperature T_ref the above terms becomes

$$\begin{array}{}\text{(7)}& {\displaystyle }\begin{array}{rl}& h_{\mathrm{sp}}={\mathit{\rho }}_{\mathrm{sp}}{c}_{\mathrm{sp}}(\mathrm{1}-{\mathit{\theta }}_{\mathrm{s}})(T-T_{\mathrm{ref}})\\ & h_{\mathrm{w}}={\mathit{\rho }}_{\mathrm{w}}{c}_{\mathrm{w}}{\mathit{\theta }}_{\mathrm{w}}\left(T\right)(T-T_{\mathrm{ref}})+{\mathit{\rho }}_{\mathrm{w}}{l}_{\mathrm{f}}{\mathit{\theta }}_{\mathrm{w}}\left(T\right)\\ & h_{\mathrm{i}}={\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{i}}{\mathit{\theta }}_{\mathrm{i}}\left(T\right)(T-T_{\mathrm{ref}}),\end{array}\end{array}$$

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, and c_i are the specific 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 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

$$\begin{array}{}\text{(8)}& h=C_T(T-T_{\mathrm{ref}})+{\mathit{\rho }}_{\mathrm{w}}{l}_{\mathrm{f}}{\mathit{\theta }}_{\mathrm{w}}\left(T\right),\end{array}$$

where $C_T={\mathit{\rho }}_{\mathrm{sp}}{c}_{\mathrm{sp}}(\mathrm{1}-{\mathit{\theta }}_{\mathrm{s}})+{\mathit{\rho }}_{\mathrm{w}}{c}_{\mathrm{w}}{\mathit{\theta }}_{\mathrm{w}}\left(T\right)+{\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{i}}{\mathit{\theta }}_{\mathrm{i}}\left(T\right)$ is the bulk-heat capacity of the soil volume V_c.

SFCCs have an inflection point (Bao et al., 2016; Hansson et al., 2004) causing a sharp change in their derivative. This non-linear 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 fundamental and important issue in frozen soil models.

There is a more subtle aspect in the integration though. Analytically, Eqs. (1), (2), and (4) are equivalent because Eqs. (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 the parabolic ones.

The apparent-heat-capacity approach, Eq. (2), can suffer from large balance errors in the presence of high non-linearities 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) and Voller et al. (1990) discussed suitable techniques.

Referring to the work by Roe (1981), it can be proven that the discrete operator of C_a has to satisfy the requirement

$$\begin{array}{}\text{(9)}& {\stackrel{\mathrm{̃}}{C}}_{a_i}^{n+\mathrm{1}/\mathrm{2}}(T_i^{n+\mathrm{1}}-T_i^n)=h\left(T_i^{n+\mathrm{1}}\right)-h\left(T_i^n\right),\end{array}$$

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 obtain

$$\begin{array}{}\text{(10)}& {\stackrel{\mathrm{̃}}{C}}_{a_i}^{n+\mathrm{1}/\mathrm{2}}{\displaystyle \frac{T_i^{n+\mathrm{1}}-T_i^n}{\mathrm{\Delta }t}},\end{array}$$

and substituting the condition Eq. (9) we obtain

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{h\left(T_i^{n+\mathrm{1}}\right)-h\left(T_i^n\right)}{T_i^{n+\mathrm{1}}-T_i^n}}{\displaystyle \frac{T_i^{n+\mathrm{1}}-T_i^n}{\mathrm{\Delta }t}}={\displaystyle \frac{h_i^{n+\mathrm{1}}-h_i^n}{\mathrm{\Delta }t}}.\end{array}$$

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, 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 (Voller et al., 1990) 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 a non-convergence problem unless under-relaxation is wisely applied, and additionally it is necessary to guarantee that the liquid fraction is in the range (0,1). With DECP, the energy equation is first solved without latent heat. 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.

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 non-linear 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 non-linear 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 timescale of the process under investigation.

Table 1 Diversity of formulations and solvers in current models of heat transfer and phase change in the cryosphere. Theoretical limitations do not necessarily affect the usability of models for their intended purpose. More details are in Appendix A.

^a The governing equation is written in only in terms of temperature and the latent heat is not included. ^b Apparent-heat-capacity formulation. ^c Enthalpy formulation. ^d Source term formulation. ^e The heat flux is written in terms of enthalpy and not of temperature as in the enthalpy formulation. ^f Decoupled energy conservation parametrization. ^g Convergence of the non-linear solver can be problematic (Casulli and Zanolli, 2010, 2012). ^h Nicolsky et al. (2007b). ⁱ Voller et al. (1990).

3 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 resolve diurnal phenomena near the ground surface, and also this often arises from limitations of the numerical schemes used. Many applications related to permafrost (Teng et al., 2020), on the other hand, only require the representation of seasonal and multi-annual variation, which can be accomplished using time steps of 1 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 non-linear system of equations, and the algorithm used for this is of great importance. Existing linearization algorithms such as the Picard or the Newton algorithm require a sufficiently accurate initial guess. As reported by Casulli and 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 non-linear system of equations. Because of this the convergence of Newton-type method is often problematic (Casulli and Zanolli, 2012). Instead, the NCZ algorithm always converges in any situation, even under the use of large time steps and grid sizes when the initial guess is chosen as suggested (Casulli and Zanolli, 2010, 2012). This is an important feature of the non-linear solver, because an inexact solution of the non-linear system is not conservative (Casulli and Zanolli, 2010).

3.1 Discretization of the domain

The domain is partitioned using an unstructured orthogonal grid (Casulli and Walters, 2000), consisting of a set of non-overlapping convex volumes Ω_i, $i=\mathrm{1},\mathrm{2},\mathrm{\dots },N_v$, separated by M internal faces Γ_j, $j=\mathrm{1},\mathrm{2},\mathrm{\dots },M$. Let 𝒜_j denote the non-zero 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 non-empty intersection, are orthogonal to each other, and have distance δ_j. Each control volume Ω_i may have an arbitrary number of faces. Let ℱ_i denote the non-empty set of faces of the ith volume, with the exclusion of boundary faces. Moreover, let 𝒫(i,j) be the neighbour of volume i that shares face j with the ith control volume so that $\mathrm{1}\le \mathcal{P}(i,j)\le N_v$ for all j∈ℱ_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

$$\begin{array}{}\text{(12)}& h_i\left(T_i^{n+\mathrm{1}}\right)=h_i\left(T_i^n\right)+\mathrm{\Delta }t\left[\sum _{j\in {\mathcal{F}}_i}{\mathrm{\Lambda }}_j^{n+\mathrm{1}}{\displaystyle \frac{T_{\mathcal{P}(i,j)}^{n+\mathrm{1}}-T_i^{n+\mathrm{1}}}{{\mathit{\delta }}_j}}+S_i^n\right],\end{array}$$

where Δt is the time step size,

$$\begin{array}{}\text{(13)}& {\mathrm{\Lambda }}_j^{n+\mathrm{1}}:={\mathcal{A}}_j\text{max}\left[{\mathit{\lambda }}_{\mathrm{i}}\left(T_i^{n+\mathrm{1}}\right),{\mathit{\lambda }}_{{\mathcal{P}}_{(}i,j)}\left(T_{{\mathcal{P}}_{(}i,j)}^{n+\mathrm{1}}\right)\right],\end{array}$$

and

$$\begin{array}{}\text{(14)}& S_i=\underset{{\mathrm{\Omega }}_i}{\int }S\mathrm{d}\mathrm{\Omega }\end{array}$$

is an optional source/sink term in volume, and h_i(T) is the ith enthalpy given by

$$\begin{array}{}\text{(15)}& h_i\left(T\right)=\underset{{\mathrm{\Omega }}_i}{\int }h\left(T\right)\mathrm{d}\mathrm{\Omega }.\end{array}$$

Equation (12) can be written in matrix form as

$$\begin{array}{}\text{(16)}& \mathit{h}\left(\mathit{T}\right)+\mathbf{A}\mathit{T}=\mathit{b},\end{array}$$

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_i^{\mathrm{0}}$, at any time step $n=\mathrm{1},\mathrm{2},\mathrm{\dots }$ Eq. (12) constitutes a non-linear system for $T_i^{n+\mathrm{1}}$, with the non-linearity affecting only the diagonal of the system and being represented by the enthalpy $h_i\left(T_i^{n+\mathrm{1}}\right)$. This set of equations is a consistent and conservative discretization of Eq. (1). Therefore, regardless of the chosen spatial and temporal resolution, $T_i^{n+\mathrm{1}}$ is a conservative approximation of the new temperature.

3.2 Solution of the non-linear system

Difficulties in solving the non-linear 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 of the enthalpy is not correctly defined. The NCZ algorithm used here was discovered by Casulli and Zanolli (2010) and overcomes these difficulties applying the Newton algorithm in two subsequent 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 (Brugnano and Casulli, 2008, 2009; Casulli and Zanolli, 2010, 2012).

For each control volume the enthalpy function h_i(T) can be defined as

$$\begin{array}{}\text{(17)}& h_i\left(T\right)=\underset{T_{\mathrm{ref}}}{\overset{T}{\int }}{C}_{\mathrm{a},i}\left(\mathit{\xi }\right)\mathrm{d}\mathit{\xi },\end{array}$$

where C_a,i(T) is defined as $C_{\mathrm{a},i}\left(T\right)={\int }_{{\mathrm{\Omega }}_i}{C}_{\mathrm{a},i}\left(T\right)\mathrm{d}\mathrm{\Omega }$. 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 non-negative function with bounded variations. The second one, C2, is that there exists a $T_i^{*}$ such that C_a,i(T) is strictly positive and non-decreasing in $(\mathrm{0},T_i^{*})$ and non-increasing in $(T_i^{*},+\mathrm{\infty })$.

Since C_a,i(T) are non-negative functions with bounded variations, they are differentiable almost everywhere, admit only discontinuities of the first kind, and can be expressed by using the Jordan decomposition (Fig. 1) as the difference of two non-negative, non-decreasing, and bounded functions:

$$\begin{array}{}\text{(18)}& \begin{array}{lll}{p}_i\left(T\right)=C_{\mathrm{a},i}\left(T\right)& q_i\left(T\right)=\mathrm{0}& \text{if }T\le T_i^{*}\\ p_i\left(T\right)=C_{\mathrm{a},i}\left(T_i^{*}\right)& q_i\left(T\right)=p_i\left(T\right)-C_{\mathrm{a},i}\left(T\right)& \text{if }T>T_i^{*}.\end{array}\end{array}$$

Using Eq. (18) the enthalpy h_i(T) can be written as $h_i\left(T\right)=h_{\mathrm{1},i}\left(T\right)-h_{\mathrm{2},i}\left(T\right)$, where

$$\begin{array}{}\text{(19)}& \begin{array}{ll}{h}_{\mathrm{1},i}\left(T\right)=h_i\left(T\right)& \text{ if }T\le T_i^{*}\\ h_{\mathrm{1},i}\left(T\right)=h_i\left(T_i^{*}\right)+C_{\mathrm{a},i}\left(T_i^{*}\right)(T-T_i^{*})& \text{ if }T>T_i^{*}\\ h_{\mathrm{2},i}\left(T\right)=\mathrm{0}& \text{ if }T\le T_i^{*}\\ h_{\mathrm{2},i}\left(T\right)=h_{\mathrm{1},i}\left(T\right)-h_i\left(T\right)& \text{ if }T>T_i^{*}.\end{array}\end{array}$$

By making use of Eq. (19), the algebraic system, Eq. (16), can be written as

$$\begin{array}{}\text{(20)}& {\mathit{h}}_{\mathbf{1}}\left(\mathit{T}\right)-{\mathit{h}}_{\mathbf{2}}\left(\mathit{T}\right)+\mathbf{A}\mathit{T}=\mathit{b}.\end{array}$$

The NCZ algorithm requires first the linearization of h₂ and then the linearization of h₁. By choosing the initial guess for the NCZ algorithm such that ${\mathit{T}}^{\mathrm{0}}\le {\mathit{T}}^{*}$, a sequence of outer iterates {T^k} is obtained from Eq. (20) by linearizing h₂ as follows:

$$\begin{array}{}\text{(21)}& \begin{array}{rl}{\mathit{h}}_{\mathbf{1}}\left({\mathit{T}}^k\right)& -\left[{\mathit{h}}_{\mathbf{2}}\right({\mathit{T}}^{k-\mathrm{1}})+\mathbf{Q}({\mathit{T}}^{k-\mathrm{1}}\left)\right({\mathit{T}}^k-{\mathit{T}}^{k-\mathrm{1}}\left)\right]\\ & +\mathbf{A}{\mathit{T}}^k=\mathit{b},\end{array}\end{array}$$

so that the outer iterates are solutions of the following non-linear system:

$$\begin{array}{}\text{(22)}& {\mathit{h}}_{\mathbf{1}}\left({\mathit{T}}^k\right)+(\mathbf{A}+{\mathbf{Q}}^{k-\mathrm{1}}){\mathit{T}}^k={\mathit{d}}^{k-\mathrm{1}},k=\mathrm{1},\mathrm{2},\mathrm{\dots },\end{array}$$

where ${\mathbf{Q}}^{k-\mathrm{1}}=\mathbf{Q}\left({\mathit{T}}^{k-\mathrm{1}}\right)$ and ${\mathit{d}}^{k-\mathrm{1}}=\mathit{b}+{\mathit{h}}_{\mathbf{2}}\left({\mathit{T}}^{k-\mathrm{1}}\right)-{\mathbf{Q}}^{k-\mathrm{1}}{\mathit{T}}^{k-\mathrm{1}}$. The resulting kth (outer) residual is derived from Eq. (20) and reads as

$$\begin{array}{}\text{(23)}& {\mathit{r}}^k={\mathit{h}}_{\mathbf{1}}\left({\mathit{T}}^k\right)-{\mathit{h}}_{\mathbf{2}}\left({\mathit{T}}^k\right)+\mathbf{A}{\mathit{T}}^k-\mathit{b}.\end{array}$$

For all $k=\mathrm{1},\mathrm{2},\mathrm{\dots }$, by setting ${\mathit{T}}^{k,\mathrm{0}}={\mathit{T}}^{k-\mathrm{1}}$, a sequence of inner iterates {T^k,l} is derived from Eq. (22) by linearizing h₁(T) as

$$\begin{array}{}\text{(24)}& \begin{array}{rl}\left[{\mathit{h}}_{\mathbf{1}}\right({\mathit{T}}^{k,l-\mathrm{1}})& +\mathbf{P}\left({\mathit{T}}^{k,l-\mathrm{1}}\right)({\mathit{T}}^{k,l}-{\mathit{T}}^{k,l-\mathrm{1}})]\\ & +(\mathbf{A}+{\mathbf{Q}}^{k-\mathrm{1}}){\mathit{T}}^{k,l}={\mathit{d}}^{k-\mathrm{1}},l=\mathrm{1},\mathrm{2},\mathrm{\dots },\end{array}\end{array}$$

so that the inner iterates are determined from the following linear systems:

$$\begin{array}{}\text{(25)}& ({\mathbf{P}}^{k,l-\mathrm{1}}+\mathbf{A}+{\mathbf{Q}}^{k-\mathrm{1}}){\mathit{T}}^{k,l}={\mathit{f}}^{k,l-\mathrm{1}},l=\mathrm{1},\mathrm{2},\mathrm{\dots },\end{array}$$

where ${\mathbf{P}}^{k,l-\mathrm{1}}=\mathbf{P}\left({\mathit{T}}^{k,l-\mathrm{1}}\right)$ and ${\mathit{f}}^{k,l-\mathrm{1}}={\mathit{d}}^{k-\mathrm{1}}-{\mathit{h}}_{\mathbf{1}}\left({\mathit{T}}^{k,l-\mathrm{1}}\right)+{\mathbf{P}}^{k,l-\mathrm{1}}{\mathit{T}}^{k,l-\mathrm{1}}$. The resulting (k,l)th (inner) residual is derived from Eq. (22) and reads as

$$\begin{array}{}\text{(26)}& {\mathit{r}}^{k,l}={\mathit{h}}_{\mathbf{1}}\left({\mathit{T}}^{k,l}\right)+(\mathbf{A}+{\mathbf{Q}}^{k-\mathrm{1}}){\mathit{T}}^{k,l}-{\mathit{d}}^{k-\mathrm{1}}.\end{array}$$

The inner and the outer iterations are terminated when $\parallel {\mathit{r}}^{k,l}\parallel <\mathit{\epsilon }$ and $\parallel {\mathit{r}}^k\parallel <\mathit{\epsilon }$, respectively, with ε being a sufficiently small pre-fixed tolerance.

Figure 1 Graphical representation of the Jordan decomposition for the enthalpy of soil using the SFCC model for a silty soil (Dall'Amico, 2010). Panel (a) shows the Jordan decomposition of C_a(T), Eq. (18). For $T=T^{*}$, C_a(T) presents a maximum: for $T<T^{*}$ it is increasing, and for $T>T^{*}$ it is decreasing. This non-monotonic behaviour causes problems when solving the non-linear system. C_a(T) is thus replaced by p(T) and q(T), two monotonic functions. Consequently, in panel (b), h(T) is replaced by h₁(T) and h₂(T), Eq. (19).

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 soil water retention curve (SWRC) models and the Clapeyron equation, which in general are difficult to implement in numerical models based on the apparent heat capacity (Kurylyk and Watanabe, 2013). Functions describing the internal energy of other substances, for instance pure water (Andreas, 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, Sect. 4.1.

4 Analytical benchmarks

The numerical model is compared for the problem of a column of freezing water, i.e. the Stefan problem, with the analytical solution presented by Neumann (cited in Kurylyk et al., 2014b), and for the problem of a column of soil with three zones – frozen, partially frozen, and unfrozen – with the analytical solution presented by Lunardini (1988).

4.1 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 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 $T(z,\mathrm{0})=T_{\mathrm{0}}$, where T₀>T_m (Fig. 2). At the surface a Dirichlet boundary condition is imposed $T(z=\mathrm{0},t)=T_{\mathrm{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

$$\begin{array}{}\text{(27)}& {\displaystyle \frac{\partial h}{\partial t}}=\mathit{\lambda }{\displaystyle \frac{{\partial }^{\mathrm{2}}T}{\partial z^{\mathrm{2}}}},\text{(28)}& {\displaystyle }T(\mathit{\zeta },t)=T_{\mathrm{m}},\text{(29)}& {\displaystyle }{\mathit{\lambda }}_{\mathrm{i}}{\displaystyle \frac{\partial T}{\partial z}}{|}_{z={\mathit{\zeta }}^{-}}\mathrm{d}t={\mathit{\lambda }}_{\mathrm{w}}{\displaystyle \frac{\partial T}{\partial z}}{|}_{z={\mathit{\zeta }}^{+}}\mathrm{d}t+l_{\mathrm{f}}\mathit{\rho }\mathrm{d}\mathit{\zeta }.\end{array}$$

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 of the two substances is consumed for the phase change.

Figure 2 Scheme showing the setting of the Neumann solution for the freezing case. Initially all water is liquid, T₀>T_m. Because of the surface boundary condition, T_s<T_m, a freezing front, ζ, propagates downward.

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 Figs. 3 and 4.

Table 2 Maximum error (m) of the freezing front position from the numerical solution with the NCZ algorithm for different space and time discretizations relative to the Neumann analytical solution. For the numerical solution the position of the freezing front has been reconstructed from the linear interpolation of the temperature profile.

For comparison, Kurylyk et al. (2014b) tested the numerical model SUTRA against the Neumann analytical solution considering a soil porosity of 0.50 m³ m⁻³. For their test the time step was of 0.04–0.4 s, the vertical spatial discretization was 0.001 m, and the parameter ϵ was increased to −0.01 ^∘C to match the analytical solution. The maximum absolute error of the freezing front position was 0.00099 m.

In our model, the choice of a small melting temperature range ϵ=0.0001 ^∘C does not affect the quality of the numerical solution even at a large time step of 3600 s. Looking at Table 2 it is clear that the choice of the time step size is somehow related to the choice of the spatial discretization: using a small time step with a coarse grid does not necessarily improve the accuracy of the position of the freezing front.

Figure 3 Propagation of the freezing front compared between the Neumann analytical and the numerical solution with the NCZ algorithm. Two space discretizations are used: (a) Δz=0.005 m and (b) Δz=0.001 m. The integration time step is Δt=3600 s.

Figure 4 Panels (a) and (c) show the temperature evolution for the Neumann analytical and the numerical solution with the NCZ algorithm at various depths for a spatial discretization Δz=0.005 m and Δz=0.001 m respectively. The integration time step is Δt=3600 s. Panels (b) and (d) show the absolute error.

We use the Neumann analytical solution to assess the robustness of the NCZ algorithm in comparison with the Newton–Raphson and globally convergent Newton methods. As reported by Dall'Amico et al. (2011), Fig. 5 represents a well known case for which the Newton–Raphson algorithm cannot converge. Instead, the solution continuously cycles between the iterates T^k−1 and T^k.

Figure 5 A scheme of the problem which illustrates how the Newton–Raphson method cannot converge towards T (Dall'Amico, 2010). In this case, the Newton–Raphson method fails to converge to T since it cycles between T^k and T^k+1 values.

While the Newton–Raphson algorithm converges to the exact solution if a good initial guess for T^k exists, this represents a severe constraint for the reliable application for an iterative algorithm in a numerical model. An improvement of the Newton–Raphson algorithm can be obtained using the globally convergent Newton scheme (Dall'Amico 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). Moreover, the Newton–Raphson and the globally convergent Newton methods do not guarantee convergence for any time step size, and this represents a limiting factor whenever reaching the convergence is necessary to reduce the time step size. For example, in (Dall'Amico et al., 2011) the comparison between the Neumann solution and GEOtop has been done with a time step of 10 s.

A comparison of the numerical solutions obtained with the Newton–Raphson algorithm, the globally convergent Newton algorithm, and the NCZ algorithm shows significant differences (Fig. 6). Newton–Raphson cannot reproduce the analytical solution even if a time step of Δt=10 s is used. The globally convergent Newton is in good agreement with the analytical solution 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.

Figure 6 Comparison between the Neumann analytical solution and the numerical solution obtained with the Newton–Raphson (N. R.), globally convergent Newton (g. c. N.), and NCZ algorithms. All the numerical simulations use the same spatial discretization Δz=0.005 m.

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 algorithm, for which convergence only depends on the right definition of Eqs. (18) and (19).

Figure 7 Comparison between the Neumann analytical solution and the numerical solution obtained with the globally convergent Newton algorithm (g. c. N.). All the numerical simulations use the same spatial discretization Δz=0.005 m and a time step size of Δt=3600 s. This figure shows as the numerical solution depends on the choice of the parameter δ.

4.2 Lunardini analytical solution

Lunardini (1988) derived an analytical solution (Appendix F) for the temporal evolution of temperature during the freezing of a semi-infinite and initially unfrozen soil column. In contrast to the Neumann analytical solution, in the Lunardini analytical solution the domain is divided into three regions (Fig. 8) on the basis of temperature: unfrozen, T<T_m; partially frozen, $T_{\mathrm{m}}<T<T_{\mathrm{f}}$; and fully frozen, T>T_f. The domain is initially unfrozen with $T=T_{\mathrm{0}}=\mathrm{4}$ ^∘C. At the left boundary condition a Dirichlet boundary condition is imposed with $T(x=\mathrm{0},t)=T_{\mathrm{s}}=-\mathrm{6}$ ^∘C, and the right boundary temperature is kept equal to the initial condition, $T(z\to \mathrm{\infty },t)=T_{\mathrm{0}}$. Because the left boundary condition T_s<0 ^∘C, a freezing front propagates from left to right.

We computed benchmark T1 proposed by the InterFrost project (, ); parameters are given in Table F1. The model agrees well with the analytical solution for all the three cases of T_m in terms of both the temperature profile, Fig. 9 and Table 3, and the freezing front position, Fig. 10 and Table 4, even with an hourly time step.

For comparison, McKenzie et al. (2007) compared the numerical model SUTRA against the Lunardini analytical solution for the cases $T_{\mathrm{m}}=-\mathrm{4}$ ^∘C and $T_{\mathrm{m}}=-\mathrm{1}$ ^∘C using a time step size of 900 s and a space resolution of 0.01 m. For the first test case the maximum absolute error was 0.01 ^∘C, and for the second test case it was 0.1 ^∘C. Their parameters, however, differ from those suggested by the InterFrost consortium, making performance comparisons difficult. In particular, their porosity was 0.05 m³ m⁻³, whereas InterFrost uses 0.336 m³ m⁻³. As this determines the amount of latent heat involved in phase change, smaller errors are suspected to occur with the parameters used by McKenzie et al. (2007).

Figure 8 Scheme showing the setting of the Lunardini problem (Ruhaak et al., 2015). Initially the domain is unfrozen with T=T₀. Because of T_s<0 on the left boundary, a freezing front propagates from left to right. X₁(t) and X(t) identify respectively the isotherm corresponding to T_m and T_f.

Figure 9 Comparison between the Lunardini solution and the numerical solution with the NCZ algorithm for the three cases of the T1 benchmark: (a) $T_{\mathrm{m}}=-\mathrm{4}$ ^∘C, (b) $T_{\mathrm{m}}=-\mathrm{1}$ ^∘C, and (c) $T_{\mathrm{m}}=-\mathrm{0.1}$ ^∘C. The colours represent different time frames. The integration time step is Δt=3600 s, and the space resolution is Δx=0.01 m.

Figure 10 Propagation of the zero isotherm for the Lunardini solution and the numerical solution with the NCZ algorithm for the three cases of T1 benchmark: (a) $T_{\mathrm{m}}=-\mathrm{4}$ ^∘C, (b) $T_{\mathrm{m}}=-\mathrm{1}$ ^∘C, and (c) $T_{\mathrm{m}}=-\mathrm{0.1}$ ^∘C. The integration time step is Δt=3600 s, and the space resolution is Δx=0.01 m.

Table 3 Maximum absolute error ^∘C of the temperature after 24 h from the numerical solution with the NCZ algorithm relative to the Lunardini analytical solution. The space resolution is Δx=0.01 m.

Table 4 Maximum error (m) of the freezing front position from the numerical solution with the NCZ algorithm relative to the Lunardini analytical solution. The space resolution is Δx=0.01 m.

5 Numerical test

In the previous sections, we have demonstrated that the proposed method can reproduce the Neumann analytical solution, as 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 d time steps. The domain is a soil column of 20 m depth that is uniformly at $T=-\mathrm{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 phase change are usually greatest near the soil surface, the thickness Δz is parametrized with an exponential function (Gubler et al., 2013):

$$\begin{array}{}\text{(30)}& \mathrm{\Delta }{z}_i=\mathrm{\Delta }{z}_{\min }(\mathrm{1}+b{)}^{i-\mathrm{1}},\end{array}$$

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 years to reach a stable thermal regime. After spin-up, we performed a simulation of 100 years.

Figure 11 Comparison of the position of the zero isotherm, panel (c), after 100 years of three simulations: using an hourly boundary condition with time step of Δt=1 h, using a daily boundary condition with a time step of Δt=1 d, and using a 10 d boundary condition with a time step of Δt=10 d. Panel (a) shows the surface temperature for the hourly, the daily, and the 10 d simulations. Panel (b) shows the deviation of the position of the zero isotherm after 100 years between the hourly and the 10 d simulation, as well as between the daily and the 10 d simulation.

Figure 11 compares the zero-isotherm position computed after 100 years for the three different time steps. Interestingly, 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 (Figs. G1 and G2). These differences can be explained by the different surface temperature used to drive the simulations (see Appendix G) and by the fact that in using a larger time step we lose accuracy in the numerical solution.

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–Nicolson method, this would incur a time step restriction to guarantee the monotonicity of the solution. As this restriction is proportional to the square of the space discretization, Δz², the Crank–Nicolson method would represent a severe constrain whenever high spatial resolution is required.

Figure 12 (a) The minimum, mean, and maximum temperature profile for the hourly simulation. Panels (b), (c), and (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 d, and with a 10 d air temperature boundary condition and Δt=10 d. All three simulations last 100 years. The maximum difference of T_mean between the hourly and daily simulation is 0.04 ^∘C, while between the hourly and 10 d simulation the difference is 0.3 ^∘C.

Figure 12 compares the minimum, mean, and maximum temperature profile respectively for the three simulations. Panel (a) shows the ground temperature envelope for the hourly simulation. The maximum envelop presents an “elbow” that is due to the zero-curtain effect Fig. G3. As can be seen in panels (b) and (d), close to the soil surface the hourly simulation presents larger 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 d average.

In the mean temperature profile, the 10 d 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 reduced to 0.003 ^∘C with a maximum deviation of 0.003 ^∘C from the initial condition, which is also equal to the mean of the forcing boundary condition.

So far we have compared the results of the model with the choice of the time step size; in Appendix G we compare with the choice of the spatial discretization. The results are still in good agreement, but it is interesting to note that the zero isotherm presents some steps, independently of the size of the time step, and some details are missed, such as the joining 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.

Moreover, we checked the mean number of iterations required to solve the non-linear 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.

6 Conclusions

We have presented a new model for simulating the ground thermal regime in the presence of freezing and thawing based on 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 occurs.

Numerical experiments demonstrated the robustness of the model by comparing results at differing temporal and spatial resolutions. Results obtained with time steps of 1 h, 1 d, and 10 d are consistent. The robustness of the numerics allows 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 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 the 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 substance simply by changing the definition of the enthalpy function and the thermal conductivity function. Examples include glacier dynamics (Aschwanden et al., 2012), snowpack 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).

Appendix A: Commonly used simulation software

The heat equation can be written in different forms that are analytically equivalent but subject to differing numerical advantages and disadvantages. In the scientific literature, several simulators, i.e. software that implements a particular model (set of equations), for solving the heat equation with freezing and thawing have been presented. Here we review commonly used frozen soil models in terms of their governing equations and methods of finding numerical solutions.

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 snowpack, glaciers, and ice sheets. Our contribution does not seek to 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 non-linearity arising from phase change, and, furthermore, broadening our review to also include some snow and glacier models supports the generalization of our findings.

A1 CLM

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 non-conservative 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–Nicolson method. The temperature profile is calculated by adopting the DECP approach. This approach does not require solving a non-linear system, since the latent heat is treated in an explicit way, but Nicolsky et al. (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 solving 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 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 an 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).

A4 GEOtop

GEOtop (Rigon et al., 2006; Endrizzi et al., 2014) is a physically based distributed model of the mass and energy balance of the hydrological cycle. It includes a module for solving the energy equation in freezing soil (Dall'Amico et al., 2011); this module can also be linked with the solver for the Richards equation. The governing equation for heat transfer is written in conservative form, but when solving the equation the apparent-heat-capacity formulation is used. A globally convergent Newton algorithm is used to deal with the non-linearities arising from phase change (Dall'Amico et al., 2011). The globally convergent Newton algorithm represents an improvement over the Newton–Raphson algorithm; however, as shown in Sect. 4.1 it does not perform as well as the NCZ algorithm, and additionally the choice of the parameter δ is non-trivial.

A5 GIPL 2.0

GIPL 2.0 simulates the ground thermal regime by solving the heat equation with phase change numerically (Marchenko et al., 2008). The governing equation is written in the conservative form, and Newton's method is used to linearize the energy equation. To overcome convergence problems when solving the non-linear system, GIPL 2.0 implements a fractional time step approach – Godunov splitting. The key point of the solution regards the treatment of the enthalpy time derivative: if a non-zero gradient of temperature exists, the time derivative is approximated with a difference derivative; otherwise, it is approximated using the analytical representation.

A6 Goodrich

Goodrich (1982) presented a one-dimensional model to simulate the ground thermal regime considering the phase change of water. The governing equation is written in the non-conservative form and does not include the latent heat of fusion. Phase 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 the Picard iteration is used to linearize the algebraic non-linear system. In their paper, Hansson et al. (2004) explain that during the Picard iteration the solution can easily oscillate whenever the temperature decreases below the melting temperature. To avoid these oscillation the temperature is reset to the critical value and the iteration restarted. Hydrus 1D adopts an empirical time step adaptation criterion. It is worthwhile to note that the modified Picard iteration was proposed by Celia et al. (1990) to solve the Richards equation – a 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 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 non-linear system is solved using a Newton–Raphson method. To overcome convergence and stability problems, three modifications were introduced (Painter, 2011). The convergence of the Newton–Raphson method can be problematic (Casulli and Zanolli, 2012).

A9 NEST

Zhang et al. (2003) developed the one-dimensional physically based Northern Ecosystem Soil Temperature (NEST) model. The heat equation is written in the source term formulation and solved with the DECP approach. The numerical method is explicit in time; thus the maximum time step is 30 min to prevent oscillations in the solution.

A10 Sergueev et al.

This is a two-dimensional model, and the governing equation is written in the enthalpy form (Sergueev et al., 2003). This model 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 Newton's method. As in GIPL 2.0, the time derivative of enthalpy is computed using either the difference derivative or the analytical derivative according with the gradient of the temperature field.

A11 SoilVision

The heat equation is written using the apparent heat capacity. The equations are solved using a finite-element solver, FlexPDE suite, both explicit and implicit in time. In the 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. (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 non-linearities are solved using the 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 into account metamorphism and settling. The first versions of Crocus (Brun et al., 1989, 1992) were not enthalpy-based. The governing equation was written in terms of temperature and water content. It was solved by using the Crank–Nicolson 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 by adopting the DECP approach (Boone and Etchevers, 2001). Even though recent work Crocus is based on a simple bucket approach for liquid water percolation (Morin et al., 2012; Lafaysse et al., 2017), D'Amboise et al. (2017) implemented a routine for water flow in the snowpack based on the Richards equation, which is characterized by non-linear behaviour like the enthalpy equation. To solve it, they adopted an approach based on the Picard iteration with variable time steps (Paniconi and Putti, 1994).

A14 SNOWPACK

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 the Richards equation, and the full Richards equation (Wever et al., 2014). The full Richards equation is solved using the 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 Institut 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, 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 by 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 Max Planck Institute Earth system model (MPI-ESM) that also includes ECHAM6 for the atmosphere and the Max Planck Institute ocean model (MPI-OM) for the ocean. JSBACH simulates 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 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.

The 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 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 cannot 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 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). SICOPOLIS (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 relates to the 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 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 non-linear 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 non-linear system. To the authors' knowledge, the enthalpy gradient method has not be used in freezing soil models.

Hewitt and Schoof (2016) presented an enthalpy-based finite-volume method for polythermal ice. To solve the equation at each time step the computational domain is explicitly divided in the cold and temperate regions, and the energy equation is solved by adopting a combination of implicit and explicit methods (Hewitt and Schoof, 2016). It is worth noting 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 higher 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 symmetric and positive definite; thus the NCZ algorithm can be applied to linearize the system. To solve the resulting linear system we use 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

Using a fully implicit formulation, the discretization of Eq. (1) reads as

$$\begin{array}{}\text{(C1)}& \begin{array}{rl}{h}_i\left(T_i^{n+\mathrm{1}}\right)& =h_i\left(T_i^n\right)\\ & +\mathrm{\Delta }t\left[\sum _{j\in {\mathcal{F}}_i}{\mathrm{\Lambda }}_j^{n+\mathrm{1}}{\displaystyle \frac{T_{\mathcal{P}(i,j)}^{n+\mathrm{1}}-T_i^{n+\mathrm{1}}}{{\mathit{\delta }}_j}}+S_i^{n+\mathrm{1}}\right],\end{array}\end{array}$$

where T_i is the temperature of the ith control volume, 𝒫(i,j) is the neighbour of volume i that shares face j with the ith control volume, δ_j is the non-zero distance between the centres of two adjacent volumes which share the jth internal face, Δt is the time step size,

$$\begin{array}{}\text{(C2)}& {\mathrm{\Lambda }}_j^{n+\mathrm{1}}={\mathcal{A}}_j\text{max}\left[{\mathit{\lambda }}_{\mathrm{i}}\left(T_i^{n+\mathrm{1}}\right),{\mathit{\lambda }}_{{\mathcal{P}}_{(}i,j)}\left(T_{{\mathcal{P}}_{(}i,j)}^{n+\mathrm{1}}\right)\right],\end{array}$$

and

$$\begin{array}{}\text{(C3)}& S_i^{n+\mathrm{1}}=\underset{{\mathrm{\Omega }}_i}{\int }S\mathrm{d}\mathrm{\Omega }\end{array}$$

is an optional source/sink term in volume; h_i(T) is the ith enthalpy given by

$$\begin{array}{}\text{(C4)}& h_i\left(T\right)=\underset{{\mathrm{\Omega }}_i}{\int }h\left(T\right)\mathrm{d}\mathrm{\Omega }.\end{array}$$

For given initial condition $T_i^{\mathrm{0}}$, at every time step $n=\mathrm{1},\mathrm{2},\mathrm{\dots },$ Eq. (C1) constitutes a fully non-linear system of equation to be solved for $T_i^{n+\mathrm{1}}$. To solve Eq. (C1), one sets $T_i^{n+\mathrm{1},\mathrm{0}}=T_i^n$. Then the Picard iterations are taken to be

$$\begin{array}{}\text{(C5)}& \begin{array}{rl}& h_i\left(T_i^{n+\mathrm{1},m+\mathrm{1}}\right)=h_i\left(T_i^n\right)\\ & +\mathrm{\Delta }t\left[\sum _{j\in {\mathcal{F}}_i}{\mathrm{\Lambda }}_j^{n+\mathrm{1},m}{\displaystyle \frac{T_{\mathcal{P}(i,j)}^{n+\mathrm{1},m+\mathrm{1}}-T_i^{n+\mathrm{1},m+\mathrm{1}}}{{\mathit{\delta }}_j}}+S_i^{n+\mathrm{1},m}\right],\end{array}\end{array}$$

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), but some iterations can be used to update the thermal conductivity to the n+1th time level Casulli and Zanolli (2010).

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

$$\begin{array}{}\text{(D1)}& U_{\mathrm{c}}=U_{\mathrm{c}}(S,V,M),\end{array}$$

where S is the entropy, V is the volume, and M is the mass of the constituents. These are the independent variables and are called extensive variables since they depend linearly on the mass of the substance. The first differential of Eq. (D1) is

$$\begin{array}{}\text{(D2)}& \mathrm{d}{U}_{\mathrm{c}}=\left({\displaystyle \frac{\partial U_{\mathrm{c}}}{\partial S}}\right)\mathrm{d}S+\left({\displaystyle \frac{\partial U_{\mathrm{c}}}{\partial V}}\right)\mathrm{d}V+\left({\displaystyle \frac{\partial U_{\mathrm{c}}}{\partial M}}\right)dM.\end{array}$$

According to Callen (1985) it is possible to define

$$\begin{array}{}\text{(D3)}& {\displaystyle }\left({\displaystyle \frac{\partial U_{\mathrm{c}}}{\partial S}}\right)\equiv T\text{, the temperature},\text{(D4)}& {\displaystyle }-\left({\displaystyle \frac{\partial U_{\mathrm{c}}}{\partial V}}\right)\equiv p\text{, the pressure},\text{(D5)}& {\displaystyle }\left({\displaystyle \frac{\partial U_{\mathrm{c}}}{\partial M}}\right)\equiv \mathit{\mu }\text{, the chemical potential}.\end{array}$$

With this notation, Eq. (D2) becomes

$$\begin{array}{}\text{(D6)}& \mathrm{d}{U}_{\mathrm{c}}=T\mathrm{d}S-p\mathrm{d}V+\mathit{\mu }\mathrm{d}M.\end{array}$$

By making use of the Legendre transformation it is possible to define the enthalpy potential H_c as

$$\begin{array}{}\text{(D7)}& H_{\mathrm{c}}(S,p,M)=U_{\mathrm{c}}(S,V,M)+pV(S,p,M).\end{array}$$

The differential of the enthalpy is

$$\begin{array}{}\text{(D8)}& \begin{array}{rl}\mathrm{d}{H}_{\mathrm{c}}& =\mathrm{d}[U_{\mathrm{c}}+pV]\\ & =T\mathrm{d}S-p\mathrm{d}V+\mathit{\mu }\mathrm{d}M+V\mathrm{d}p+p\mathrm{d}V\\ & =T\mathrm{d}S+\mathit{\mu }\mathrm{d}M+V\mathrm{d}p.\end{array}\end{array}$$

If we assume that the transformation occurs at constant pressure and volume, then Eq. (D6) becomes

$$\begin{array}{}\text{(D9)}& \mathrm{d}{U}_{\mathrm{c}}=T\mathrm{d}S+\mathit{\mu }\mathrm{d}M,\end{array}$$

and Eq. (D8) becomes

$$\begin{array}{}\text{(D10)}& \mathrm{d}{H}_{\mathrm{c}}=T\mathrm{d}S+\mathit{\mu }\mathrm{d}M.\end{array}$$

Hence, from Eqs. (D9) and (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 terms 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

$$\begin{array}{}\text{(E1)}& h\left(T\right)=\left\{\begin{array}{l}{\mathit{\rho }}_{\mathrm{w}}{c}_{\mathrm{w}}(T-T_{\mathrm{ref}})+{\mathit{\rho }}_{\mathrm{w}}{l}_{\mathrm{f}}\text{if }T\ge T_{\mathrm{m}}\\ {\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{i}}(T-T_{\mathrm{ref}})\text{if }T<T_{\mathrm{m}}-\mathit{ϵ}\\ {\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{i}}(T-\mathit{ϵ}-T_{\mathrm{ref}})+h^{\prime }(T-(T_{\mathrm{m}}-\mathit{ϵ}\left)\right)\text{ otherwise},\end{array}\right.\end{array}$$

where the singularity of the enthalpy function at T=T_m has been linearized with

$$\begin{array}{}\text{(E2)}& h^{\prime }={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}{c}_{\mathrm{w}}(T_{\mathrm{m}}-T_{\mathrm{ref}})+{\mathit{\rho }}_{\mathrm{w}}{l}_{\mathrm{f}}-{\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{i}}(T_{\mathrm{m}}-\mathit{ϵ}-T_{\mathrm{ref}})}{\mathit{ϵ}}},\end{array}$$

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 Sect. 3.2. It is worth underlining that the temperature range ϵ can be chosen to be 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 grain boundaries (Langham, 1974; Nye and Frank, 1973); (b) freezing often occurs below the melting point when nucleation is relevant; and (c) in porous media such as soil, ice melts across a range of temperatures due to the Gibbs–Thomson effect in pores and surface affects at the interface 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 from a convergence problem with a time step of 3600 s.

Figure E1 (a) Comparison between the enthalpy function of pure water and the enthalpy function used in the numerical model. (b) Note that the energy jump due to the latent heat at T_m=0 ^∘C has been linearized, and the steepness is controlled by the parameter ϵ.

The thermal conductivity is defined as

$$\begin{array}{}\text{(E3)}& \mathit{\lambda }\left(T\right)=\left\{\begin{array}{ll}{\mathit{\lambda }}_{\mathrm{w}}& \text{if }T\ge T_{\mathrm{m}}\\ {\mathit{\lambda }}_{\mathrm{i}}& \text{if }T<T_{\mathrm{m}}.\end{array}\right.\end{array}$$

Defining the constant

$$\begin{array}{}\text{(E4)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{w}}={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{w}}{c}_{\mathrm{w}}}}{\mathit{\alpha }}_i={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{i}}}},\text{(E5)}& {\displaystyle }A={\displaystyle \frac{T_{\mathrm{m}}-T_{\mathrm{s}}}{\mathrm{erf}\left(\mathit{\gamma }\right)}}B={\displaystyle \frac{T_{\mathrm{m}}-T_{\mathrm{0}}}{\mathrm{erf}\left(\mathit{\gamma }\sqrt{{\displaystyle \frac{{\mathit{\alpha }}_i}{{\mathit{\alpha }}_{\mathrm{w}}}}}\right)}},\end{array}$$

the moving boundary function is

$$\begin{array}{}\text{(E6)}& \mathit{\zeta }\left(t\right)=\mathrm{2}\mathit{\gamma }\sqrt{{\mathit{\alpha }}_{i}t}\text{for }t>\mathrm{0},\end{array}$$

where the coefficient γ can be found by solving the following equation:

$$\begin{array}{}\text{(E7)}& \mathit{\gamma }\sqrt{{\mathit{\alpha }}_i}{l}_{\mathrm{f}}\mathit{\rho }-{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{i}}}{\sqrt{\mathit{\pi }{\mathit{\alpha }}_i}}}A{e}^{-{\mathit{\gamma }}^{\mathrm{2}}}-{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{w}}}{\sqrt{\mathit{\pi }{\mathit{\alpha }}_w}}}B{e}^{{\mathit{\gamma }}^{\mathrm{2}}{\displaystyle \frac{{\mathit{\alpha }}_i}{{\mathit{\alpha }}_{\mathrm{w}}}}}=\mathrm{0}.\end{array}$$

Finally the analytical solutions for problem with the Dirichlet boundary condition for the thawed and frozen zones are

$$\begin{array}{}\text{(E8)}& \left\{\begin{array}{ll}T(z,t)=T_{\mathrm{s}}+{\displaystyle \frac{T_{\mathrm{m}}-T_{\mathrm{s}}}{\mathrm{erf}\left(\mathit{\gamma }\right)}}\mathrm{erf}\left({\displaystyle \frac{z}{\mathrm{2}\sqrt{{\mathit{\alpha }}_{i}t}}}\right)& \mathrm{0}<z<\mathit{\zeta }\left(t\right)\\ T(z,t)=T_{\mathrm{0}}+{\displaystyle \frac{T_{\mathrm{m}}-T_{\mathrm{0}}}{\mathrm{erfc}\left(\mathit{\gamma }\sqrt{{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{i}}}{{\mathit{\alpha }}_{\mathrm{w}}}}}\right)}}\mathrm{erfc}\left({\displaystyle \frac{z}{\mathrm{2}\sqrt{{\mathit{\alpha }}_{\mathrm{w}}t}}}\right)& z>\mathit{\zeta }\left(t\right).\end{array}\right.\end{array}$$

Table E1 Input parameters for the comparison between Neumann analytical solution and the numerical solution with the NCZ algorithm.

Appendix F: Lunardini analytical derivation

The solution for the Lunardini problem (i.e. the Lunardini solution) as described by McKenzie et al. (2007) is given by the following set of equations:

$$\begin{array}{}\text{(F1)}& {\displaystyle }{T}_{\mathrm{1}}=(T_{\mathrm{m}}-T_{\mathrm{s}}){\displaystyle \frac{\mathrm{erf}\left({\displaystyle \frac{x}{\mathrm{2}\sqrt{{\mathit{\alpha }}_{\mathrm{1}}t}}}\right)}{\mathrm{erf}\left(\mathit{\psi }\right)}}+T_{\mathrm{s}},\text{(F2)}& {\displaystyle }{T}_{\mathrm{2}}=(T_{\mathrm{m}}-T_{\mathrm{f}}){\displaystyle \frac{\mathrm{erf}\left({\displaystyle \frac{x}{\mathrm{2}\sqrt{{\mathit{\alpha }}_{\mathrm{4}}t}}}\right)-\mathrm{erf}\left(\mathit{\gamma }\right)}{\mathrm{erf}\left(\mathit{\gamma }\right)-\mathrm{erf}\left(\mathit{\psi }\sqrt{{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{1}}}{{\mathit{\alpha }}_{\mathrm{4}}}}}\right)}}+T_{\mathrm{f}},\text{(F3)}& {\displaystyle }{T}_{\mathrm{3}}=(T_{\mathrm{0}}-T_{\mathrm{f}}){\displaystyle \frac{-\mathrm{erfc}\left({\displaystyle \frac{x}{\mathrm{2}\sqrt{{\mathit{\alpha }}_{\mathrm{3}}t}}}\right)}{\mathrm{erfc}\left(\mathit{\psi }\sqrt{{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{4}}}{{\mathit{\alpha }}_{\mathrm{3}}}}}\right)}}+T_{\mathrm{0}},\end{array}$$

where T₁, T₂, and T₃ are the temperatures at distance, x, from the temperature boundary for the frozen, partially frozen, and unfrozen zone respectively; erf and erfc are the error function and the complementary error function respectively; T₀, T_m, T_f, and T_s are the temperatures of the initial condition, the solidus, the liquidus, and the boundary temperature, respectively; and α₁ and α₃ are the thermal diffusivities for the frozen and unfrozen zone respectively, defined as ${\mathit{\lambda }}_{\mathrm{1}}/C_{\mathrm{1}}$ and ${\mathit{\lambda }}_{\mathrm{3}}/C_{\mathrm{3}}$, where C₁ and C₃ are the volumetric bulk-heat capacities of the frozen and unfrozen zones. The thermal diffusivity of the partially frozen zone is assumed to be constant across the transition region, and the thermal diffusivity with latent heat term included, α₄, is defined as

$$\begin{array}{}\text{(F4)}& {\mathit{\alpha }}_{\mathrm{4}}={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{2}}}{C_{\mathrm{2}}+{\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{d}}{l}_{\mathrm{f}}\mathrm{\Delta }\mathit{\xi }}{(T_{\mathrm{f}}-T_{\mathrm{m}})}}}},\end{array}$$

where γ_d is the dry unit of soil solids, and $\mathrm{\Delta }\mathit{\xi }={\mathit{\xi }}_{\mathrm{0}}-{\mathit{\xi }}_f$, where ξ₀ 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 $\mathrm{0}\le x\le X_{\mathrm{1}}\left(t\right)$ the temperature is T₁ and X₁(t) is given by

$$\begin{array}{}\text{(F5)}& X_{\mathrm{1}}\left(t\right)=\mathrm{2}\mathit{\psi }\sqrt{{\mathit{\alpha }}_{\mathrm{1}}t};\end{array}$$

from $X_{\mathrm{1}}\left(t\right)\le x\le X\left(t\right)$ the temperature is T₂, where X(t) is given by

$$\begin{array}{}\text{(F6)}& X\left(t\right)=\mathrm{2}\mathit{\gamma }\sqrt{{\mathit{\alpha }}_{\mathrm{4}}t};\end{array}$$

and for x≥X(t) the temperature is T₃. The unknowns, ψ and γ, are solving the set of these two equations:

$$\begin{array}{}\text{(F7)}& {\displaystyle \frac{T_{\mathrm{m}}-T_{\mathrm{s}}}{T_{\mathrm{m}}-T_{\mathrm{f}}}}{\exp }^{-{\mathit{\psi }}^{\mathrm{2}}(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{1}}/{\mathit{\alpha }}_{\mathrm{4}})}={\displaystyle \frac{{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{1}}}}\sqrt{{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{1}}}{{\mathit{\alpha }}_{\mathrm{4}}}}}\mathrm{erf}\left(\mathit{\psi }\right)}{\mathrm{erf}\left(\mathit{\gamma }\right)-\mathrm{erf}\left(\mathit{\psi }\sqrt{{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{1}}}{{\mathit{\alpha }}_{\mathrm{4}}}}}\right)}},\text{(F8)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{(T_m-T_f){\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{3}}}}}{T_m-T_f}}{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{3}}}{{\mathit{\alpha }}_{\mathrm{4}}}}{\exp }^{-{\mathit{\gamma }}^{\mathrm{2}}(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{4}}/{\mathit{\alpha }}_{\mathrm{3}})}\\ & ={\displaystyle \frac{\mathrm{erf}\left(\mathit{\gamma }\right)-\mathrm{erf}\left(\sqrt{{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{1}}}{{\mathit{\alpha }}_{\mathrm{4}}}}}\mathit{\psi }\right)}{\mathrm{erf}\left(\mathit{\gamma }\sqrt{{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{4}}}{{\mathit{\alpha }}_{\mathrm{3}}}}}\right)}}.\end{array}\end{array}$$

Table F1 Input parameters for the comparison between the Lunardini analytical solution and the numerical solution with the NCZ algorithm.

^a The first value refers to $T_{\mathrm{m}}=-\mathrm{0.1}$ ^∘C, the second value to $T_{\mathrm{m}}=-\mathrm{1}$ ^∘C, and the third value to $T_{\mathrm{m}}=-\mathrm{4}$ ^∘C.

Appendix G: Numerical test

As explained in Sect. 5, comparing the position of the zero isotherm after 100 years using three different time steps – hourly, daily, and 10 d – there is no significant deviation in the results. The larger deviation occurs when the zero isotherm is shallow: at the beginning of the thawing season as well as the freezing one.

At the beginning of the thawing season, Fig. G1, there is a time lag of about 1 month between the beginning of the thawing season for the hourly simulation and the 10 d simulation. This can be attributed to different surface temperature used to drive the simulations. In particular, in the case of the hourly simulation it is possible to see the oscillations of the position of the zero isotherm, panel (c), related to the oscillation of the surface temperature around 0 ^∘C, panel (a).

Figure G1 Detail of the beginning of the thawing season for the year 1999. Panel (a) shows the surface temperature for the hourly, the daily, and the 10 d simulations. Panel (b) shows the deviation of the position of the zero isotherm after 100 years between the hourly and the 10 d simulation, as well as between the daily and the 10 d simulation. Panel (c) shows the position of the zero isotherm after 100 years for the three simulations. In panel (b) there is a time lag of about 1 month between the beginning of thawing season for the hourly simulation and the 10 d one (dashed grey line). This can be attributed to the different surface temperature used to drive the simulations.

Figure G2 shows the detail of the freezing season. In panel (c) it is possible to note that when the zero isotherm is deep there is a good agreement between the three simulations. The main differences occur at the soil surface since with larger time steps the signal of the surface boundary condition is smoothed and does not oscillate around 0 ^∘C. Moreover, by using an hourly time step and a daily time step it is possible to capture the joining of the downward and upward freezing front, while this is not possible with the 10 d time step since the joining occurs in between two consecutive time steps.

Figure G2 Detail of the beginning of the freezing season for the year 1999. Panel (a) shows the surface temperature for the hourly, the daily, and the 10 d simulations. Panel (b) shows the deviation of the position of the zero isotherm after 100 years between the hourly and the 10 d simulation, as well as between the daily and the 10 d simulation. Panel (c) shows the position of the zero isotherm after 100 years for the three simulations. The joining of the downward and upward freezing front is captured by the hourly and the daily simulations (c). It is interesting to note that for the 10 d simulation the joining occurs in between two consecutive time steps.

As explained in Sect. 5, the maximum temperature profile, Fig. 12d, presents an “elbow” due to the so-called zero-curtain effect. The zero-curtain effect, Fig. G3, is the period of time during which the temperature remains nearly constant and very close to the freezing point because of the latent heat released during the phase change of water.

Figure G3 Hourly temperature at 1.5 m depth for the year 1999. Note the prolonged period of 43 d (11 October until 23 November) 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.

Figure G4 shows the temperature envelopes for the hourly, the daily, and the 10 d 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 G4 Temperature profile envelope considering λ_w=λ_i. (a) The minimum, mean, and maximum temperature profile for the hourly simulation. Panels (b), (c), and (d) show the comparison of the minimum, mean, and maximum temperature profile respectively for the three simulations: with an hourly surface temperature boundary condition and Δt=1 h, with a daily surface temperature boundary condition and Δt=1 d, and with a 10 d surface temperature boundary condition and Δt=10 d. 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 0.003 ^∘C.

Figure G5 shows a comparison of the zero-isotherm position by using a coarser space discretization. Again, the three simulations with the hourly time step, the daily time step, and the 10 d 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 by the hourly or by the daily simulations.

Figure G5 Comparison of the position of the zero isotherm, panel (c), after 100 years of three simulations: using an hourly boundary condition with time step of Δt=1 h, using a daily boundary condition with a time step of Δt=1 d, and using a 10 d boundary condition with a time step of Δt=10 d. Panel (a) shows the surface temperature for the hourly, the daily, and the 10 d simulations. Panel (b) shows the deviation of the position of the zero isotherm after 100 years between the hourly and the 10 d simulation, as well as between the daily and the 10 d simulation. Panel (c) shows the position of the zero isotherm. It is interesting to note that all the three curves are characterized by the presence of some steps due to the usage of a coarser grid. This trend is independent of the size of the time step. Another consequence of this is that the joining of the downward and upward freezing front is not captured by the hourly or by the daily simulations.

Table G1 Input parameters for the numerical tests.

^a We used two different space discretizations. The thickness of the ground layer is parametrized as $\mathrm{d}{z}_i=\mathrm{d}{z}_{\min }(\mathrm{1}+b{)}^{(i-\mathrm{1})}$ (Gubler et al., 2013).

For this numerical test we checked the mean number of iterations required to solve the non-linear system with the NCZ 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–Raphson nor the globally convergent Newton converges: they always reach the maximum number of iterations allowed with a consequent increase in the computational cost.

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 $\mathit{\epsilon }=\mathrm{10}\times {\mathrm{10}}^{-\mathrm{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–Raphson and the globally convergent Newton does not converge so it always reaches the maximum number of iterations allowed.