Gully erosion is a problem; therefore, it must be predicted using highly accurate predictive models to avoid losses caused by gully development and to guarantee sustainable development. This research investigates the predictive performance of seven multiple-criteria decision-making (MCDM), statistical, and machine learning (ML)-based models and their ensembles for gully erosion susceptibility mapping (GESM). A case study of the Dasjard River watershed, Iran uses a database of 306 gully head cuts and 15 conditioning factors. The database was divided 70:30 to train and verify the models. Their performance was assessed with the area under prediction rate curve (AUPRC), the area under success rate curve (AUSRC), accuracy, and kappa. Results show that slope is key to gully formation. The maximum entropy (ME) ML model has the best performance (AUSRC = 0.947, AUPRC = 0.948, accuracy = 0.849 and kappa = 0.699). The second best is the random forest (RF) model (AUSRC = 0.965, AUPRC = 0.932, accuracy = 0.812 and kappa = 0.624). By contrast, the TOPSIS (Technique for Order Preference by Similarity to Ideal Solution) model was the least effective (AUSRC = 0.871, AUPRC = 0.867, accuracy = 0.758 and kappa = 0.516). RF increased the performance of statistical index (SI) and frequency ratio (FR) statistical models. Furthermore, the combination of a generalized linear model (GLM), and functional data analysis (FDA) improved their performances. The results demonstrate that a combination of geographic information systems (GIS) with remote sensing (RS)-based ML models can successfully map gully erosion susceptibility, particularly in low-income and developing regions. This method can aid the analyses and decisions of natural resources managers and local planners to reduce damages by focusing attention and resources on areas prone to the worst and most damaging gully erosion.