The investigation into the freeze-thaw (F-T) damage holds paramount theoretical significance in comprehending the mechanisms underlying F-T disasters, predicting catastrophes, and devising engineering protection systems. F-T damage in rock evolves progressively from the surface towards the interior, however, current models exhibit conspicuous deficiencies in simulating such progressive damage. Drawing upon experimental findings concerning the evolution of F-T damage in water-saturated rock, an improved discrete element model (DEM) for rock F-T damage is developed. In this model, thermal conduction model and volume expansion theory are integrated to describe the volumetric expansion characteristics during the freezing of internal pore water, thus effectively simulating the initiation and progression characteristics of frost heave cracks. We validated this model through comparisons with theoretical analyses and experiments, exploring the effects of model parameters. Employing moment tensor, the spatial gradient distribution patterns of F-T damage in water-saturated rock from the vantage point of fracture energy was analyzed. The results reveal that at −20 °C, 20 F-T cycles are requisite for the propagation of damage from the outer layer to the center of the Sichuan sandstone samples. The extent of damage in the outer layer exceeds that of the center, with the maximum fracture energy density of the outer layer being approximately 38.7 times that of the center. This model aptly delineates the progressive development of F-T damage characteristics and holds promise for further applications in identifying and delineating F-T damage zones in cold region rock engineering.