A hybrid analytical-numerical technique for solving soil temperature during the freezing process

The freeze-thaw cycle associated with climatic seasonality is a common phenomenon in cold regions affecting a wide range of subsurface processes. Due to the complex and highly nonlinear nature of the associated hydrologic processes, transient freeze-thaw dynamics are conventionally quantified in a numerical way. Here we present a hybrid analytical-numerical scheme for solving one-dimensional soil (or porous media) temperature profiles when the soil profile is subjected to unidirectional freezing (or thawing) conditions. This scheme divides the partially-frozen soil into multi-layers, each with constant thermal parameters and fixed-temperature boundaries. Temperature profiles within each layer were obtained by solving multiple moving-boundary problems. The proposed hybrid analytical-numerical scheme was tested into a freezing test of silty clay in a permafrost region on the Qinghai-Tibetan Plateau, and its solution was in good agreement with the finite element numerical solution. Results show that the proposed multi-layer method adapted well to the changes in unfrozen water content and thermal properties of soil over a wide range of subzero temperatures. By contrast, the freezing front's migration rate and penetration depth calculated by Neumann's classical solution, which only considers two zones (frozen and unfrozen), was found to be underestimated. As for our proposed multi-layer solution, by dividing the subsurface domain into many layers with smaller proportion ratios (thinner layers close to the freezing front), there was a slower penetration rate of the freezing front resulting in shallower penetration depth. The predicted profiles of temperature, thermal conductivity and diffusivity, heat flux, and dynamics of the freezing front were significantly impacted by the shape of the soil freezing curves and the magnitude of soil grain thermal conductivity, especially for the accuracy of long-term predictions.