Rill erosion represents an important transition stage from sheet erosion to gully erosion, and it is an important part of the physical process of water erosion on hillslopes. Due to its dynamic and random nonlinear properties, rill initiation and development processes on slopes are very difficult to accurately simulate and reproduce. Most of the existing erosion and sediment yield models focus on soil loss evaluation rather than to analyze the evolution of the surface. In this study, the theory and methods of Cellular Automata (CA) were applied to simulate rill erosion process while also taking into consideration soil erosion and hydraulics of the runoff system. The objective was to develop a numerical simulation model of rill erosion on the loess hillslopes in China. In this CA model, the hillslopes surface is subdivided into a series of discrete spatial cells with the same geometric features. The flow direction was determined by D8 algorithm and equal assignment rule. The amounts of interchanged water and sediment are computed using the Manning formula and the empirical sediment transport equation. The simulated results using the CA model were compared and validated with actual artificial rainfall experiments. It was confirmed that the CA model produced similar patterns of runoff and sediment as in the actual experiment. Moreover, the CA model also provided more details about the dynamic and random development of rills with spatial and temporal variability on the loess hillslopes.