A new method is presented to account for phase changes in a fully implicit numerical model for coupled heat transport and variably saturated water flow involving conditions both above and below zero temperature. The method is based on a mixed formulation for both water flow and heat transport similar to the approach commonly used for the Richards equation. The approach enabled numerically stable, energy- and mass-conservative solutions. The model was evaluated by comparing predictions with data from laboratory column freezing experiments. These experiments involved 20-cm long soil columns with an internal diameter of 8 cm that were exposed at the top to a circulating fluid with a temperature of −6°C. Water and soil in the columns froze from the top down during the experiment, with the freezing process inducing significant water redistribution within the soil. A new function is proposed to better describe the dependency of the thermal conductivity on the ice and water contents of frozen soils. Predicted values of the total water content compared well with measured values. The model proved to be numerically stable also for a hypothetical road problem involving simultaneous heat transport and water flow. The problem was simulated using measured values of the surface temperature for the duration of almost 1 yr. Since the road was snow-plowed during winter, surface temperatures varied more rapidly, and reached much lower values, than would have been the case under a natural snow cover. The numerical experiments demonstrate the ability of the code to cope with rapidly changing boundary conditions and very nonlinear water content and pressure head distributions in the soil profile.