Land degradation from gully erosion poses a significant threat to the Erer watershed in Eastern Ethiopia, particularly due to agricultural activities and resource exploitation. Identifying erosion-prone areas and underlying factors using advanced machine learning algorithms (MLAs) and geospatial analysis is crucial for addressing this problem and prioritizing adaptive and mitigating strategies. However, previous studies have not leveraged machine learning (ML) and GIS-based approaches to generate susceptibility maps identifying these areas and conditioning factors, hindering sustainable watershed management solutions. This study aimed to predict gully erosion susceptibility (GES) and identify underlying areas and factors in the Erer watershed. Four ML models, namely, XGBoost, random forest (RF), support vector machine (SVM), and artificial neural network (ANN), were integrated with geospatial analysis using 22 geoenvironmental predictors and 1,200 inventory points (70% used for training and 30% for testing). Model performance and robustness were validated through the area under the curve (AUC), accuracy, precision, sensitivity, specificity, kappa coefficient, F1 score, and logarithmic loss. The relative slope position is most influential, with 100% importance in SVM and RF and 95% importance in XGBoost, while annual rainfall (AR) dominated ANN (100% importance). Notably, XGBoost demonstrated robustness and superior prediction/mapping, achieving an AUC of 0.97, 91% accuracy, 92% precision, and 81% kappa while maintaining a low logloss (0.0394). However, SVM excelled in classifying gully resistant/susceptible areas (97% sensitivity, 98% specificity, and 91% F1 score). The ANN model predicted the most areas with very high gully susceptibility (13.74%), followed by the SVM (11.69%), XGBoost (10.65%), and RF (7.85%) models, while XGBoost identified the most areas with very low susceptibility (70.19%). The ensemble technique was employed to further enhance GES modeling, and it outperformed the individual models, achieving an AUC of 0.99, 93.5% accuracy, 92.5% precision, 97.5% sensitivity, 95.4% specificity, 85.8% kappa, and 94.9% F1 score. This technique also classified the GES of the watershed as 36.48% very low, 26.51% low, 16.24% moderate, 11.55% high, and 9.22% very high. Furthermore, district-level analyses revealed the most susceptible areas, including the Babile, Fedis, Harar, and Meyumuluke districts, with high GES areas of 32.4%, 21.3%, 14.3%, and 13.6%, respectively. This study offers robust and flexible ML models with comprehensive validation metrics to enhance GES modeling and identify gully prone areas and factors, thereby supporting decision-making for sustainable watershed conservation and land degradation prevention.