Gully erosion is an important soil degradation process, which under climate changes is projected to increase. Therefore, better understating of factors controlling gully erosion and prediction of gully headcuts' (GHs) location is still highly relevant. This study aimed to examine the spatial distribution of GHs and to assess the importance of pedological (i.e. aggregate stability, organic matter, bulk density, silt, clay, and sand content) and topographical factors (i.e. altitude, slope length, gradient, and aspect) using summary statistics and the maximum entropy (MaxEnt) model. The study was conducted in the loess-covered region of NE Iran. The highly precise data of 287 GHs locations were obtained by extensive fieldwork and the interpretation of UAV images. The spatial distribution of GHs was evaluated by univariate pair correlation function and O-ring statistics. The spatial effect of GHs density controlling factors was assessed by the cumulative density correlation function Cm,K(r). Variable importance was analyzed using the MaxEnt model, which was also for the susceptibility modelling of GHs. The results of univariate tests showed the aggregated distribution of GHs. The Cm,K(r) analyses indicated that the areas characterized by higher values of bulk density, aggregate stability, and organic matter content have lower GHs density, whereas the areas with high silt content and higher slope gradient have higher GHs density. According to the MaxEnt, there is no one single factor responsible for GHs location, but rather the combination of topographical and pedological factors with the predominance of slope gradient (0.86) and silt content (0.57). The MaxEnt modelling of GHs susceptibility has revealed that the best accuracy (0.958) is given when all pedological and topographical factors are used in the model. The susceptibility maps prepared in the study can be used for soil conversation and land use planning and, consequently, for sustainable development in the region.