Gully erosion susceptibility (GES) mapping is crucial for controlling gully erosion hazards and has become a significant focus of global research and management efforts. Machine learning models have proven effective in this field. However, in areas with different terrain complexity, the model shows significant variation in optimal resolution and algorithms, factor importance and spatial distribution of the model results, which limits their broader application. This study compares GES mapping in two small watersheds: one located in the complex terrain of the Loess Plateau and the other in the relatively flat terrain of the Northeast China Mollisol region. The model predictive accuracy was evaluated using 30% of the datasets that were excluded from model training. The results revealed that: 1) significant differences in optimal resolution of GES mapping in the two regions, which were 1-2.5 m for the Mollisol region, and 2.5-5 m for the Loess Plateau. The extreme boosting tree (XGBoost) algorithm achieved the best simulation results compared to random forest (RF) and gradient boosting decision tree (GBDT) in both regions. 2) Slope gradient and contributing area influenced gully distribution in both watersheds, with land use being critical in the Loess Plateau and distance from streams more important in the Mollisol region. 3) In the Loess Plateau watershed, 25% of the area was highly susceptible to gully erosion, while only 1% of the Mollisol watershed was highly susceptible. This research compared GES mapping in two watersheds with different terrain complexity, which would be beneficial for better use of machine learning in gully research.