The complex drift movement of snow particles during snowfall weather can last for several hours or even days, thereby causing dynamic changes in snow distribution on a roof. A sufficient change in a snow boundary will affect the local flow field. In this study, a probability-based half-steady Eulerian–Lagrangian method was proposed to simulate snow change on a roof. The snowfall process was divided into N steady phases, and the snow boundary of each phase was updated according to the snow deposition probability on the roof in the previous phase. In the numerical simulation, the motion of particles was described using the Lagrangian method. Calculation time was effectively reduced by computing the deposition probability of snow particles on the roof. This method was used to calculate snow distribution on a stepped flat roof model under different wind velocities. In addition, a wind tunnel test of snowfall simulation was conducted by spraying high-density silica sand particles. Results of the numerical simulation and wind tunnel test exhibited good agreement with the measured results, which accurately reproduced the typical snow distribution on the lower roof.