Soil erosion is a major soil degradation threat to land worldwide, and gully erosion is an important soil erosion process. The lack of research on the spatial variation in gully erosion has constrained planning for gully control at a large scale. In this study, high-resolution historical and present remote sensing images were used to investigate the areal gully growth rate (RA) with systematic sampling methods on loess-covered tablelands at a regional scale. Gully growth (not disturbed by human excavation or landfill) was identified in 74 of 98 sampling watersheds in the 10,000-km(2) study area, and the average ratio of gully growth area to the gully area was 0.159 % (1969 similar to 2021), with a coefficient of variation of 0.006. Moran's I index showed a significant spatial autocorrelation with RA in the study area. Runoff depth, drainage area, and soil bulk density were significantly correlated with RA in the sampling watersheds. GeogDetector revealed that the spatial variation in RA was significantly influenced by runoff depth, precipitation, soil bulk density at the regional scale. The spatial patter n of runoff depth had the largest effect on the spatial variation in RA when comparing the individual effect of each factor, while the precipitation and soil bulk density best explained the spatial patterns of areal gully growth at the regional scale. A random forest regression model performed wel l for simulating watershed RA , with R-2 = 0.829. These results improve our understanding of the mechanisms of gully erosion in loess-covered regions.