• A novel analytical-numerical scheme for calculating temperature profiles in porous media with temperature-dependent thermal properties during the freezing process; • The hybrid analytical-numerical method can deal with different types of nonlinear soil freezing functions; • Neumann's two-layer solution underestimates the penetration rate and depth of the freezing front; • The profiles of temperature, equivalent thermal conductivity and diffusivity, conductive heat flux, and dynamics of the freezing front were significantly impacted by the shape of the unfrozen water content curve and the magnitude of soil grain thermal conductivity. 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.