Assuming you have got the physics right and you wrote the right equations, which is not given for granted when you deal with freezing soils, you have to solve the equations. This paper, deals with this last topic: given the freezing soil equation it implements a new algorithm to solve it. This algorithm was invented by Casulli and Zanolli in their 2010 paper. They called it nested Newton, we renamed it NCZ from Newton-Casulli-Zanolli. It was implemented for solving Richards equation which present a spiky term called hydraulic capacity, that poses serous challenges to the convergence of the solver. We extended here to a new equation with the same type of terms. In fact all equations that involve phase transitions have terms of this type.
To someone a new algorithm for integration of some equation can seem a minor achievement but, while in some type of simulation, the numerical errors of traditional methods can be somewhat constrained, in most of the simulation they do not and tend to increase up to a point that any prediction either quantitative or qualitative becomes useless. Obviously our case has an enormous effect when dealing with simulations of permafrost areas under the threat of climate change. If this introduction makes you curious, you can find the the preprint at The Cryosphere page, by clicking on the Figure above.
Casulli, Vincenzo, and ZANOLLI. 2010. “A Nested Newton-Type Algorithm for Finite Colume Methods Solving Richards’ Equation in Mixed Form.” SIAM Journal of Scientific Computing 32 (4): 2225–73.
Tubini, Niccolò, Stephan Gruber, and Riccardo Rigon. “A Method for Solving Heat Transfer with Phase Change in Ice or Soil That Allows for Large Time Steps While Guaranteeing Energy Conservation.”