Objective This study aims to construct a gully erosion susceptibility assessment model in order to provide an important scientific basis for early warning of gully erosion occurrence and its precise management in the black soil region of Northeast China. Methods Kedong County in Heilongjiang Province， a typical thick-layer black soil region， was selected as the study area. Erosion gullies were extracted based on 0.4 m high-resolution remote sensing images. Redundant variables were eliminated using the variance inflation factor （VIF）. Four machine learning models covering linear， nonlinear， and ensemble learning architectures-logistic regression （LR）， support vector machine （SVM）， random forest （RF）， and extreme gradient boosting （XGBoost）-were compared. The SHAP （Shapley Additive exPlanations） interpretability model was introduced to quantify the contribution of each influencing factor to gully erosion susceptibility. Results 1） XGBoost demonstrated the best predictive capability for gully erosion susceptibility （training set AUC=0.81， test set AUC=0.75）， followed by the SVM and RF （test set AUC=0.74）， while LR showed the weakest predictive capability （AUC=0.67）. 2） The gully erosion susceptibility predicted by the four models all exhibited a spatial distribution pattern of "low in the periphery and high in the center" at the county scale， but there were significant differences in the susceptibility levels. The RF model predicted the largest proportion of high gully erosion susceptibility， accounting for 45.0%. The gully erosion susceptibility predicted by the SVM and XGBoost models in the study area was predominantly moderate and very low levels， accounting for 56.2% and 59.2%， respectively. In contrast， the LR prediction results indicated that low gully erosion susceptibility was dominant， covering 39.4% of the study area. 3） SHAP analysis revealed that the digital elevation model （DEM） was the dominant factor influencing gully erosion susceptibility （SHAP value=0.55）， exhibiting a significant positive effect. The stream power index （SPI） had the second strongest positive driving effect on gully erosion susceptibility （SHAP value=0.30）， while the normalized difference vegetation index （NDVI） showed a clear negative effect on gully erosion susceptibility （SHAP value=0.37）. Conclusion In the thick-layer black soil region of Kedong County， significant differences are observed in the gully erosion susceptibility levels evaluated by the four machine learning models. The integration of the XGBoost model with the SHAP method enhances the interpretability of machine learning models in predicting gully erosion susceptibility. Moreover， gully erosion susceptibility is jointly driven by the coupled effects of "terrain-dominated control， vegetation regulation， and confluence hydrodynamics triggering". These findings deepen the understanding of gully erosion occurrence and provide theoretical guidance for optimizing gully erosion control measures in black soil regions.