A Monte Carlo algorithm to simulate the isothermal recrystallization process of snow is presented. The snow metamorphism is approximated by two mass redistribution processes: surface diffusion and sublimation-deposition. The algorithm is justified and its parametrization is determined. The simulation results are compared to experimental data, in particular, the temporal evolution of the specific surface area and the ice thickness. We find that the two effects of surface diffusion and sublimation-deposition can accurately model many aspects of the isothermal metamorphism of snow. Furthermore, it is shown that sublimation-deposition is the dominant contribution for temperatures close to the melting point, whereas surface diffusion dominates at temperatures far below the melting point. A simple approximation of gravitational compaction is implemented to simulate density change.