Abstract. Heat and water movement in variably saturated freezing soils is a tightly coupled phenomenon. Strong coupling of water and heat movement in frozen soils moves considerable amounts of water from warmer to colder zones. The coupling is a result of effects of sub-zero temperature on soil water potential, heat carried by water moving under pressure gradients, and dependency of soil thermal and hydraulic properties on soil water content. This makes water and heat movement in variably saturated soils a highly non-linear process in mathematical terms. This study presents a one-dimensional cellular automata (direct solving) model to simulate coupled heat and water transport with phase change in variably saturated soils. The model is based on first order mass and energy conservation principles. The water and energy fluxes are calculated using first order empirical forms of Buckingham–Darcy's law and Fourier's heat law, respectively. The water-ice phase change is handled by integrating along experimentally determined soil freezing curve (unfrozen water content and temperature relationship) obviating the use of apparent heat capacity term. This approach highlights a further subtle form of coupling one in which heat carried by water perturbs the water content – temperature equilibrium and exchange energy flux is used to maintain the equilibrium rather than affect temperature change. The model is successfully tested against analytical and experimental solutions. Setting up a highly non-linear coupled soil physics problem with a physically based approach provides intuitive insights into an otherwise complex phenomenon.