Summary: We consider coupled models of flow and energy with phase change in thawing soils in permafrost at the pore- and Darcy scales. We develop the underlying models from first principles and define robust and conservative discretization algorithms based on cell-centered finite volume discretization on Cartesian grids and implicit time stepping for individual models. To couple the models we consider sequential approaches as well as time-lagging of some of the nonlinear properties. We also derive practical new surrogate models to obtain Darcy scale data from the pore-scale simulations. Finally, we study the sensitivity of models to the data in realistic scenarios. In particular we find that for the dynamics at the time scale of days, the nonlinear flow parameters have most significant impact on pressure, while temperature is less sensitive to data.