In this joint work with many students and collaborators we describe robust and conservative computational schemes for coupled process of flow, deformation, and energy with phase change (TpHM) in the soils in permafrost regions. The models present challenges due to the free boundary of freezing/thawing, strong dependence of constitutive parameters on the micro-physics of TpHM, disparate time scales, and micro- and macro heterogeneity. We also discuss how to get the data for the Darcy scale [TpH] ad [HM] models from the models at the pore-scale by computational upscaling.