This paper describes a nonlinear elasto-plastic simulation of freezing and thawing of rock. A mathematical formulation is described in which deformation, fluid flow and heat flow are fully coupled. A non-linear elasto-plastic constitutive relationship is presented and a two dimensional (plane stress) numerical modeling is performed based on the finite element method applied to thermo-poro-elastoplasticity. It is assumed that the Mohr-Coulomb's failure criterion is valid for yield locus and plastic potential. The numerical scheme employed in the code accommodates phase change of pore-water from liquid to solid (ice). The primary aim of this paper is to compare the temperature transfer and deformation prediction obtained from the numerical code with those obtained from the freezing and thawing experiments. It is found from the numerical simulation that a relatively good prediction can be made of temperature transfer and deformation behavior. The numerical code has also been applied to a hypothetical cavern problem to demonstrate its applicability.