## 1 Introduction  
Gully erosion is a geomorphic process that affects landscape evolution [1, 2, 3] and can be considered the main process of sediment production within small watersheds [4, 5, 6]. Gullying phenomena are recognised as one of the leading causes of land degradation in many countries worldwide, causing massive loss of soil, severe economic devaluation, and heavy impacts on ecosystem services [3, 7, 8, 9, 10]. In the last decades, gully development was favoured by the combination of extreme rainfall events and unsuitable land use changes, such as converting forests to cultivated areas or uncontrolled deforestation and agriculture activities [10, 11].  
The gully erosion processes significantly affect several countries in the Mediterranean area, causing soil degradation and impacts on agriculture production, vital resources, and benefits of human beings [10; 11; 12; 13; 14; 15]. For instance, many zones of the Calabria region (Southern Italy) are notably affected by intense water erosion processes as gullying phenomena, giving rise to serious soil loss every year [16; 17; 18]. In this regard, detecting the likely prone areas to gullying development is fundamental in directing land conservation strategies and sustainable progress [11; 19].  
Implementing the prediction models for mapping gully erosion susceptibility is one of the most important challenges in assessing potential risks and reducing negative impacts on the land. Gullies' formation and spatial distribution are controlled by the complex interactions of a large variety of geo-environmental factors, including geology, soil, topography, rainfall, land cover, land use, and human activities [1; 20; 21; 22; 23; 24]. Thus, the core principle of gully erosion susceptibility modelling is quantifying the interplay between gullies' spatial distribution and the geo-environmental variables affecting these phenomena.  
In the last decade, many statistical and machine learning methods combined with geographic information systems (GIS) and remote sensing technologies were proposed for evaluating and mapping gully erosion susceptibility in several areas worldwide [13; 23; 25; 26]. The main bi- and multi-variate statistical methods employed in the previous research are frequency ratio [27], evidential belief function [28], certainty factors [29], information value [28], index of entropy [28], weights of evidence [30], logistic regression [25], and geographically weighted regression [13; 31]. Instead, among the most known machine learning methods, there are random forest [32], artificial neural networks [22; 33], support vector machine [33], boosted regression tree [32], and maximum entropy [20; 34]. Moreover, machine learning methods are often identified as the most reliable methods for mapping gully erosion susceptibility [26; 35].  
In this context, this work aims to test the reliability of the maximum entropy machine learning method combined with GIS-based tools for mapping gully erosion susceptibility in the Mesima catchment, Calabria region (South Italy). In the Mesima basin, geographic processes, like gullying, periodically cause considerable soil erosion. Therefore, it becomes crucial to map the future gully erosion-prone areas to mitigate the adverse effects of these processes and for suitable land management, especially for agricultural areas.