Soil erosion is a serious problem affecting numerous countries, especially, gully erosion. In the current research, GIS techniques and MARS (Multivariate Adaptive Regression Splines) algorithm were considered to evaluate gully erosion susceptibility mapping among others. The study was conducted in a specific section of the Gorganroud Watershed in Golestan Province (Northern Iran), covering 2142.64 km2 which is intensely influenced by gully erosion. First, Google Earth images, field surveys, and national reports were used to provide a gully-hedcut evaluation map consisting of 307 gully-hedcut points. Eighteen gully erosion conditioning factors including significant geoenvironmental and morphometric variables were selected as predictors. To model sensitivity of gully erosion, Multivariate Adaptive Regression Splines (MARS) was used while the Area Under the Receiver Operating Characteristic (ROC) Curve (AUC), drawing ROC curves, efficiency percent, Yuden index, and kappa were used to evaluate model efficiency. We used two different scenarios of the combination of the number of replications, and sample size, including 90%/10% and 80%/20% with 10 replications, and 70%/30% with five, 10, and 15 replications for preparing gully erosion susceptibility mapping (GESM). Each one involves a various subset of both positive (presence), and negative (absence) cases. Absences were extracted as randomly distributed individual cells. Therefore, the predictive competency of the gully erosion susceptibility model and the robustness of the procedure were evaluated through these datasets. Results did not show considerable variation in the accuracy of the model, with altering the percentage of calibration to validation samples and number of model replications. Given the accuracy, the MARS algorithm performed excellently in predictive performance. The combination of 80%/20% using all statistical measures including SST (0.88), SPF (0.83), E (0.79), Kappa (0.58), Robustness (0.01), and AUC (0.84) had the highest performance compared to the other combinations. Consequently, it was found that the performance of MARS for modelling gully erosion susceptibility is quite consistent while changes in the testing and validation specimens are executed. The intense acceptable prediction capability of the MARS model verifies the reliability of the method employed for use of this model elsewhere and gully erosion studies since they are qualified to quickly generating precise and exact GESMs (gully erosion sensitivity maps) to make decisions and management edaphic and hydrologic features.