# Numerical model of soil freezing

Following finite element framework (Zienkiewicz and Taylor, 1991), we approximate

where is a ``value'' of temperature at the -th finite element grid, and is the so-called basis function. After some standard manipulations, we derive a system of differential equations

where is the vector consisting of temperature values, and are the capacitance and stiffness matrices, respectively. A further refinement, which is often used in finite element modeling of phase change problems, is to exploit a so-called ``lumped'' formulation, i.e. the capacitance matrix is diagonal:

where is one if , or zero otherwise. Dalhuijsen and Segal (1986) provides justification for the lumped formulation on noting that it is computationally advantageous and avoids oscillations in numerical solutions when used in conjunction with the backward Euler scheme:

The main difficulty in numerical modeling of soil freezing/thawing is in consistent
calculation of the derivative in (11), where is
not a continuously differentiable function defined by (2). In many reviews,
it is proposed to employ the *enthalpy temporal averaging* to calculate . We suggest an approach that incorporates ideas of temporal averaging just
to evaluate the rapidly changing by defining as

We note that an advantage of this definition is that it does not compute temporal averaging of the heat capacity, and hence reduces numerical computations, and at the same time preserves numerical accuracy of the original idea. Even though that we avoid direct computation of the derivative by replacing it by the finite difference in (13), the calculation of is still necessary to obtain a solution of the so-called adjoint model used to construct the gradient of the cost function.

We emphasize again that is not a continuously differentiable function, and hence some regularization of is required to avoid lost of accuracy in solving the adjoint problem. We propose the following regularization of the unfrozen water content in (2):

where if then , and hence . Here, is an infinitely differentiable function, and consequently posses the same property. One of the advantages of this is that we increase a convergence rate of solving the non-linear heat equation (12) by Richardson iterations, precisely find a solution to the adjoint problem, and hence find the location of the global minimum more faster than without regularization.