A method of modeling freeze-thaw cycles of naturally deposited snowpacks is presented. The model involves the Stefan condition as an independent governing equation on the exterior moving boundary to calculate snowpack thinning, flow of water through a variably saturated layered porous medium as described by the Richards equation, and heat conduction with a phase change. The heat conduction problem was treated in two ways. Local heat conduction between a snow grain and its surrounding water film was treated by using a simple energy balance. Global heat conduction with a phase change (the Stefan problem) was introduced to calculate the space-time temperature distribution. In order to handle multiple interior moving boundaries, a specific form of the enthalpy formulation was used for heat conduction with a phase change. Changing material properties were considered according to the calculated meltwater refreezing. ¿ American Geophysical Union 1994