Gully erosion is a significant threat to the sustainability of soil in Mediterranean basins. Despite its impact, there is a lack of research providing accurate regional-scale cartography of complete gully network. This study aims to automatically map the gully network in the olive-growing landscapes of the Guadalquivir basin (Spain) using Machine Learning (ML) algorithms: Random Forest (RF), Support Vector Machine (SVM), Decision Tree (DT), and Logistic Regression (LR). We integrated these models with 17 predictive variables (including hydrotopographic, climatic, and edaphic factors) and the Gully Head Initiation (GHI) index. RF was the most suitable model, achieving an Area Under the Curve (AUC) of 0.91 and an F1-score of 0.83 and enabled the delineation of a gully network totalling 8439.05 km. Variable importance analysis revealed that flow accumulation (17.33 %) and the GHI index (nearly 30%) were the primary predictors, with the Rainy Day Normal (RDN)-based formulation outperforming the maximum daily precipitation (Pmax)-based one. Spatially, countryside hills landscapes exhibited the highest gully densities (42.50 m/ha). The results demonstrate the effectiveness of combining ML with physically-based indices to generate high-resolution gully cartography for soil conservation planning in Mediterranean olive groves.