Is Overgrazing Really Influencing Soil Erosion ?

Soil erosion is a serious problem spread over a variety of climatic areas around the world. The main purpose of this paper is to produce gully erosion susceptibility maps using different statistical models, such as frequency ratio (FR) and information value (IV), in a catchment from the northeastern part of Romania, covering a surface of 550 km2. In order to do so, a total number of 677 gullies were identified and randomly divided into training (80%) and validation (20%) datasets. In total, 10 conditioning factors were used to assess the gully susceptibility index (GSI); namely, elevation, precipitations, slope angle, curvature, lithology, drainage density, topographic wetness index, landforms, aspect, and distance from rivers. As a novelty, overgrazing was added as a conditioning factor. The final GSI maps were classified into four susceptibility classes: low, medium, high, and very high. In order to evaluate the two models prediction rate, the AUC (area under the curve) method was used. It has been observed that adding overgrazing as a contributing factor in calculating GSI does not considerably change the final output. Better predictability (0.87) and success rate (0.89) curves were obtained with the IV method, which proved to be more robust, unlike FR method, with 0.79 value for both predictability and success rate curves. When using sheepfolds, the value decreases by 0.01 in the case of the FR method, and by 0.02 in the case of the success rate curve for the IV method. However, this does not prove the fact that overgrazing is not influencing or accelerating soil erosion. A multi-temporal analysis of soil erosion is needed; this represents a future working hypothesis.


Introduction
Soil erosion became a worldwide issue with regard to the material losses associated (reduction of soil fertility and cultivated areas) with this phenomenon [1].Among the most common types of erosional processes, we can distinguish erosion by water [2], wind erosion [3], and harvest erosion [4].At a European level, significant research is conducted to assess [5] and reduce [6] the negative effects of soil erosion.One of the most erosive processes acknowledged for water erosion is gullying [7].
In this study, we will refer to the gully erosion susceptibility mapping based on a GIS approach and using two statistical models in a catchment with an area of 550 km 2 from the northeastern part of Romania; besides the normal conditioning factors, grazing was added as a conditioning factor.Gully erosion can be defined as an erosion process in which deep channels are generated by runoff water removing topsoil to a certain depth [8].Once the gully channels are formed, the process of growing is very fast; this is also induced by the local environmental factors, climate, and land use changes [9].
Gully erosion represents one of the most erosive processes worldwide and leads to significant land degradation for agricultural lands, sedimentation of adjacent rivers and lakes, loss of livelihoods, Over the last few years, statistical modelling for gully erosion has encountered an upward trend; whether employing one or more statistical models to determine gully erosion susceptibility, many authors are evaluating gully erosion susceptibility based on some local environmental factors and the weights derived from each method: the information value (IV) method [22], logistic regression [23], multivariate adaptive regression splines (MARS) [24], weights of evidence (WOE); frequency ratio (FR) [25], and index of entropy (IOE) [8].
Located in the eastern part of Europe, Romania is well known for its degraded lands [26] and high gully erosion potential; a few studies refer to the conditions of gully initiation, development, and monitoring [27][28][29][30][31]. Gullying has been identified as the process with the highest sediment source, and summarises 26.4% of the total erosion processes in Romania [28].In order to correlate gully erosion initiation, it is needed to envisage areas where further gully development will occur or new gullies might initiate [23], adding a factor that has been neglected until now; overgrazing, this being one of the main purposes of this paper.
However, there is only one study that employs statistical modelling to gully erosion in Romania; the study was made to predict gully erosion in the southern part of Romania [32].Taking into consideration the high potential of gully erosion in the northeastern part of Romania, this study comes to fill the gap, by employing statistical modelling to gully erosion susceptibility.High and very high susceptible areas will be correlated to the negative effects on the cultural heritage of the area (107 Neolithic sites).At an international level, there is a limited number of studies that tackle this issue [31,33].This study comes as a necessity in an area with a very high potential to soil erosion and Over the last few years, statistical modelling for gully erosion has encountered an upward trend; whether employing one or more statistical models to determine gully erosion susceptibility, many authors are evaluating gully erosion susceptibility based on some local environmental factors and the weights derived from each method: the information value (IV) method [22], logistic regression [23], multivariate adaptive regression splines (MARS) [24], weights of evidence (WOE); frequency ratio (FR) [25], and index of entropy (IOE) [8].
Located in the eastern part of Europe, Romania is well known for its degraded lands [26] and high gully erosion potential; a few studies refer to the conditions of gully initiation, development, and monitoring [27][28][29][30][31]. Gullying has been identified as the process with the highest sediment source, and summarises 26.4% of the total erosion processes in Romania [28].In order to correlate gully erosion initiation, it is needed to envisage areas where further gully development will occur or new gullies might initiate [23], adding a factor that has been neglected until now; overgrazing, this being one of the main purposes of this paper.
However, there is only one study that employs statistical modelling to gully erosion in Romania; the study was made to predict gully erosion in the southern part of Romania [32].Taking into consideration the high potential of gully erosion in the northeastern part of Romania, this study comes to fill the gap, by employing statistical modelling to gully erosion susceptibility.High and very high susceptible areas will be correlated to the negative effects on the cultural heritage of the area (107 Neolithic sites).At an international level, there is a limited number of studies that tackle this issue [31,33].This study comes as a necessity in an area with a very high potential to soil erosion and with a tradition of sheep breeding; obtaining quantitative data and locating the precise areas that have a high and very high susceptibility to gully erosion will help local authorities and stakeholders in hazard mitigation, a proper management of soil and water resources, and mitigation towards the rich cultural heritage of the area.Therefore, the specific objectives of this study are (1) spatially predict the gully erosion in the Bahluiet , river basin; (2) employing two statistical models, FR and IV; (3) to evaluate the final susceptibility maps performance by using ROC curves; (4) to analyse the final susceptibility maps with and without the overgrazing as a conditioning factor; (5) to determine the number of Neolithic sites located in areas with high and very high susceptibility to gully erosion; and (6) to investigate the final susceptibility maps in the framework of big infrastructure projects (A8 motorway).
Europe's cultural heritage is in danger from natural hazards [34] and anthropic interventions [35]; the cultural heritage from the northeastern part of Romania is no exception [19,[36][37][38][39].However, significant efforts are made every year in order to highlight the importance of Neolithic sites in establishing the Romanian identity and to promote them from a tourism point of view [40].Besides the approach towards cultural heritage, the final maps will be analysed as a possible working direction in order to reduce time and money for one of the biggest infrastructure projects from the northeastern part of Romania-the A8 motorway.
Two working scenarios will be made; in the first one, the final gully susceptibility index (GSI) will be assessed using all the conditioning factors, and the second scenario will involve realising the final GSI by adding the sheepfold density as a conditioning factor.

Study Area and Archaeological Background
Bahluiet , catchment is located in the northeastern part of Romania (Figure 2a) and has a surface of 550 km 2 .It is a right side tributary of Bahlui River (catchment surface 2023 km 2 ), one of the most important rivers in this part of the country.Bahluiet , River has a length of 50 km, a sinuosity coefficient of 1.23, an average slope of 13%, and average catchment altitude of 163 m a.s.l.[41].
Water 2018, 10, x FOR PEER REVIEW 3 of 16 with a tradition of sheep breeding; obtaining quantitative data and locating the precise areas that have a high and very high susceptibility to gully erosion will help local authorities and stakeholders in hazard mitigation, a proper management of soil and water resources, and mitigation towards the rich cultural heritage of the area.Therefore, the specific objectives of this study are (1) spatially predict the gully erosion in the Bahluieț river basin; (2) employing two statistical models, FR and IV; (3) to evaluate the final susceptibility maps performance by using ROC curves; (4) to analyse the final susceptibility maps with and without the overgrazing as a conditioning factor; (5) to determine the number of Neolithic sites located in areas with high and very high susceptibility to gully erosion; and (6) to investigate the final susceptibility maps in the framework of big infrastructure projects (A8 motorway).Europe's cultural heritage is in danger from natural hazards [34] and anthropic interventions [35]; the cultural heritage from the northeastern part of Romania is no exception [19,[36][37][38][39].However, significant efforts are made every year in order to highlight the importance of Neolithic sites in establishing the Romanian identity and to promote them from a tourism point of view [40].Besides the approach towards cultural heritage, the final maps will be analysed as a possible working direction in order to reduce time and money for one of the biggest infrastructure projects from the northeastern part of Romania-the A8 motorway.
Two working scenarios will be made; in the first one, the final gully susceptibility index (GSI) will be assessed using all the conditioning factors, and the second scenario will involve realising the final GSI by adding the sheepfold density as a conditioning factor.

Study Area and Archaeological Background
Bahluieț catchment is located in the northeastern part of Romania (Figure 2a) and has a surface of 550 km 2 .It is a right side tributary of Bahlui River (catchment surface 2023 km 2 ), one of the most important rivers in this part of the country.Bahluieț River has a length of 50 km, a sinuosity coefficient of 1.23, an average slope of 13%, and average catchment altitude of 163 m a.s.l.[41].The lithological setting of the area is characterised by non-cohesive deposits of Bessarabian age (an alternation of marls, clays with intrusions of sand, and oolithic limestone); the fine granulation of these lithological deposits led to the triggering and developing of geomorphological processes, The lithological setting of the area is characterised by non-cohesive deposits of Bessarabian age (an alternation of marls, clays with intrusions of sand, and oolithic limestone); the fine granulation of these lithological deposits led to the triggering and developing of geomorphological processes, especially gullying (Figure 3).More details regarding the geomorphic characteristics of this area can be found in the literature [35][36][37][38][39].
Water 2018, 10, x FOR PEER REVIEW 4 of 16 especially gullying (Figure 3).More details regarding the geomorphic characteristics of this area can be found in the literature [35][36][37][38][39]. Cucuteni culture is one of the most representative prehistoric cultures from Eastern Europe and is part of the well-known Cucuteni-Ariușd-Trypillia Cultural Complex (approximately 350,000 km 2 ).Compelling all the Neolithic sites from the databases of the National Archaeological Registry (RAN), the National Heritage Institute (INP), the Institute of Cultural Memory (cIMEC), and field expeditions with archaeologists, resulted in a total number of 107 Neolithic sites (Figure 2b); our efforts were to index as many sites as possible in order to have a complete archaeological inventory.Neolithic sites from this part of the country have a high density and represent one of the most iconic archaeological assets from the northeastern part of Romania.

FR Method
FR method represents a quantitative and recognised method to generate gully erosion susceptibility maps with high accuracy [25].The main principle of this method is that future gullies will occur in places with the same geographical conditions as those in the past; the method is based Cucuteni culture is one of the most representative prehistoric cultures from Eastern Europe and is part of the well-known Cucuteni-Arius , d-Trypillia Cultural Complex (approximately 350,000 km 2 ).Compelling all the Neolithic sites from the databases of the National Archaeological Registry (RAN), the National Heritage Institute (INP), the Institute of Cultural Memory (cIMEC), and field expeditions with archaeologists, resulted in a total number of 107 Neolithic sites (Figure 2b); our efforts were to index as many sites as possible in order to have a complete archaeological inventory.Neolithic sites from this part of the country have a high density and represent one of the most iconic archaeological assets from the northeastern part of Romania.

FR Method
FR method represents a quantitative and recognised method to generate gully erosion susceptibility maps with high accuracy [25].The main principle of this method is that future gullies will occur in places with the same geographical conditions as those in the past; the method is based on the observed relationship between the distribution of gullies and individual gully-induced factors.In order to determine the frequency ratio for each factor class (with the help of Equation ( 1)), the ratio between gully occurrence and non-occurrence was calculated.
where E is the number of pixels with the landslide for each factor, F is the total number of gullies, M is the number of pixels in the class area, and L is the total number of pixels [25].The weights of each factor were obtained and by summarising their weights, the gully susceptibility index (GSI) was calculated with the following equation (Equation ( 2)): where FR i represents the frequency ratio for each factor, and FR represents the area where the gullies occurred.

IV Method
IV method is a bivariate statistic approach that was introduced for the first time by van Westen [42].It is widely used in studies to calculate gully erosion susceptibility; the method is recognised for its simplicity and robustness is recommended for large areas and where the spatial data is limited [43].The basic principle of this method is the distribution of gullies in each factor class and is calculated with the help of the following equation: where IV is the weight of a certain class i of factor j; F ij represents the gully density within class i of the factor j; F represents the total gully density in the area of interest; L ij is the number of gullies in a certain class i of the factor j; L T represents the total number of gullies; P ij is the number of pixels in a certain class i of factor j; and P L represents the total number of pixels in the study area.Each conditioning factor was combined with the gully map, in this way the weighting value (IV) was obtained for each parameter class (Table 1, column 7).Positive values show a high incidence of gully occurrence, while negative values are associated with a very low probability of gully occurrence.
In order to follow the evolution of the Băiceni-Cucuteni gully, a series of old maps [35,36] were used, along with a Leica GPS 1200 System in RTK (Real-Time Kinematic) mode; the surveys were made to determine whether the mitigation measures had a positive effect on the gully evolution.The last surveys were conducted in May 2018.

Gully Inventory Map and Conditioning Factors
A complete and precise gully inventory map is the key to produce high-quality gully erosion susceptibility maps, as well as to calibrate and validate the final susceptibility models.As there are no gully inventories available in Romania, or for the study area, a complete inventory was made using topographic maps, Google Earth images, and field investigations.Gullies were manually digitised as polygons; it resulted that only 0.9% of the study area is affected by landslide processes.From the total number of 677 gullies, 542 gullies (80%) were randomly selected for model training, and the remaining 135 gullies (20%) were used to test the model's performance (Figure 2b).The total number of 107 Neolithic sites was compiled from the available databases, and field trips with archaeologists.
During the digitisation process, it was observed that gully erosion is closely related to another process-landslides.Many gullies are initiated by the unconsolidated deposits of the landslides; therefore, there are higher chances of gully triggering around landslides.Out of the total number of 677 gullies, 457 gullies intersect with a landslide (there is a total of 764 landsides) [39].In order to realize the final susceptibility maps, a set of ten environmental factors were used; namely, elevation, precipitations, slope angle, curvature, lithology, drainage density, topographic wetness index, landforms, aspect, and distance from rivers.To these ten conditioning factors, we have added sheepfolds as a conditioning factor.We try to figure out if overgrazing is influencing the final GSI and if better validation results are obtained in this way.
Precipitations-the distribution of precipitation regulates the soil's water content and is a significant factor in calculating gully erosion susceptibility.In the study area, the values of this parameter have been concluded according to [41], and distributed into four classes (mm/year) (Figure 4b); 500-550, 550-600, 600-650, and 650-700.The precipitation class with the highest probability of contributing to gully erosion is 600-650 mm/year.
Slope angle is considered to be one of the most important factors when it comes to slope stability, being frequently used in preparing gully erosion susceptibility maps [25].This parameter is derived from DEM (Digital Elevation Model) with the help of ArcGIS Spatial Analyst Tools.In the study area, the slope was divided into four classes 0-3, 3-7, 7-14, and >14 (Figure 4c).
The curvature is defined as the rate of change of slope gradient in a particular direction; the values represent the morphology of the topography [44].Negative curvature represents concave, zero curvature represents flat, and positive curvature represents convex areas; therefore, the map was classified into three classes (Figure 4d).
Topographic wetness index (TWI) represents the effect of topography on the size and location of saturated source areas of runoff triggering; this parameter was developed in a rainfall-runoff model named TOPMODEL [45].The soil moisture affects the material on the slopes, thus diminishing soil stability.It is a common factor used in gully erosion susceptibility studies [8].Within the study area, the topographic wetness index was calculated and classified into four classes −7-−2, −1-−0.1,0-30, and >30 (Figure 4e).
The lithology factor is a significant variable in the analysis of natural hazards, but especially for gully erosion; the characteristics of lithological units have distinct gully erosion susceptibility values, and thus represent the base on which gullies are occurring [25].The lithological units were extracted from the soil maps, scale 1:10.000, and comprises four main lithological classes: sandstone (20.42%), limestone (4.11%), clay and sand (60.06%), and clay (15.41%) (Figure 4f).
The landforms parameter was derived with the help of ArcGIS-Topography Tools [46].The main landforms calculated with the help of ArcGIS are plains (52%), followed by open slopes (45%), U-shaped valleys (1.3%), upper slopes (1.2%), deeply incised streams (0.06%), high ridges (0.04%), shallow valleys, and mid-slope ridges (0.001%) (Figure 4g).The landforms with the highest probability of contributing to gully erosion are the U-shaped valleys, mid-slope ridges, deeply incised streams, and open slopes.and if better validation results are obtained in this way.
Elevation represents a widely used conditioning factor in gully susceptibility analysis [25,43].Within the study area, this parameter was obtained from the digital elevation model (5 × 5 m/pixel) and reclassified into eight classes using natural breaks (Jenks), as follows: 50.2-100, 100-150, 150-200, 200-250, 250-300, 300-350, 350-400, and >400 m (Figure 4a).A high incidence of gullies is encountered in the elevation class 200-250 m.Slope aspect represents another important factor in estimating gully erosion susceptibility; the slope aspect has a direct influence on gully erosion because it controls the vegetation cover, sunlight exposition, moisture, and evapotranspiration [25,47].Dominant in the area are north-facing slopes and flat areas (Figure 4h).Drainage density is another important factor in assessing gully erosion susceptibility [8,25]; high values of this parameter are associated with large surface runoff ratio.This parameter was calculated using the line density tool in ArcGIS and reclassified into four classes as follows <0.64, 0.65-1.2,1.3-1.8,and 1.9-3.2km/km 2 (Figure 4i).
Distance from rivers parameter is a significant conditioning factor because streams decrease slope stability by erosion of slopes and facilitate the evacuation of the eroded material from upslope areas [25,48]; in our study area, five different buffer classes were created as follows 0-150, 150-300, 300-450, 450-600, and >600 m (Figure 4j).
Another significant factor in triggering and acceleration of soil erosion processes is overgrazing [49]; there are many studies worldwide, from different geographic areas, highlighting the negative effects of overgrazing.One of the approaches is referring to the desertification [50,51], erosion of mountain pastures [52], degradation of forests [53], and so on.Many studies have discussed and documented the influence of overgrazing on the development of gully erosion [54], but none have managed to "quantify" it; an area with extensive studies is the southern part of Africa [55][56][57].
Transhumance is a very old habit of the Romanian sheep breeding community; many of them came from Transylvania during the autumn towards the rich area in grass of Moldavian Plain [58], especially within our study area.The area is known in the literature as <<Poarta Târgului Frumos>> (Târgul Frumos Gate) [59].The number of sheepfolds has increased from 119 to 175 from 1984 until 2012.For this study, the sheepfolds were digitized as a point from the aerial images of 2012; the next step was to make a kernel point density analysis.The data were then reclassified into four classes of density as follows 0-0.3, 0.3-0.9,0.9-1.5, and 1.5-3 sheepfolds per km 2 (Figure 4k).
The A8 Motorway Târgu Mures , -Ias , i-Ungheni (Figure 4l) represents one of the biggest infrastructure projects aimed to connect Moldavia and Transylvania as part of the Pan-European Transport Corridors no.IV and IX; the motorway will have a length of approximately 310 km, according to the National Company for Road Infrastructure Administration (CNAIR).A significant part will transit our study area and this study may be of great importance; in Romania it is well known that many infrastructure projects are delayed because of the discovery of significant archaeological remains and of areas with a high potential of natural hazards (e.g., landslides and gully erosion).

Results and Discussion
FR and IV were used to calculate the GSI.The surface and distribution (%) of each susceptibility class were calculated for each method.The classes of the GSI are as follows: low, medium, high, and very high susceptibility.The method with the highest predictability rate, which is IV, was used to evaluate the number of Neolithic sites that are in danger, as well as how the future A8 motorway will further develop under the high and very high gully erosion susceptibility areas.

GSI Using FR
This method selected curvature, lithology, precipitations, landforms, and distance from rivers as being the most important when it comes to gully erosion initiation and development.In the case of curvature, high values of frequency ratio are connected to concave surfaces (0.64).When it comes to lithology, the second most important factor related to gully erosion, clay and sand are the most susceptible (0.41), followed by a clay lithology (0.29).Regarding the quantity of precipitation, high values of the frequency (0.47) are related to precipitations ranging between 600-650 mm/year class; higher values of precipitations are not necessarily related to a higher incidence of gully initiation because in this area, specifically torrential summer rains occur.Landform classes four and nine, which are U-shaped valleys and mid-slope ridges, respectively, have higher frequency ratio weights (0.33 and 0.24, respectively).The parameter distance from river shows that the highest number of gullies occur in class 0-150 m (0.40).These five factors are followed by aspect, TWI, elevation, slope angle, and drainage density.The factor sheepfold density has the same weight as the factor drainage density.High values of frequency (0.28) are in the class with 1.5-3 sheepfolds/km 2 .
Water 2018, 10, 1077 10 of 16 The final GSI is calculated and shown in Figure 5a.According to the gully susceptibility map developed with the help of the FR method, 18.71% (102.72 km 2 ) of the total catchment is found under a very high susceptibility, 52.37% (287.56 km 2 ) high susceptibility; medium and low susceptibility represent 24.03% (131.93 km 2 ) and 4.89% (26.80 km 2 ) from the total area, respectively.High and very high classes occupy more than half of the basin's surface (over 70%), which is an alarming fact for the cultural heritage sites of the area.When calculating the GSI including the sheepfold density factor (Figure 5b), 19.81% (108.73 km 2 ) of the total catchment is found under a very high susceptibility, 51.45% (282.5 km 2 ) with high susceptibility; however, medium and low susceptible classes represent 23.87% (131.06 km 2 ) and 4.87% (26.72 km 2 ) from the total area, respectively.As can be observed, there are no significant differences between the maps and percentages when we added sheepfold density as a conditioning factor.
Water 2018, 10, x FOR PEER REVIEW 10 of 16 Figure 5. Gully susceptibility index produced using (a) frequency ratio (FR) method, (b) FR method including sheepfold density as a contributing factor, (c) information value (IV) method, and (d) IV method including sheepfold density as a contributing factor.

GSI Using IV
The relation of each gully erosion conditioning factor and its spatial relationship with gully triggering is shown in Table 1, column 7. A positive value indicates a high correlation between the class factor and gully erosion phenomena, while negative values highlight a very low correlation with gully erosion phenomena.High values are related to the following conditioning factors: landforms (U-shaped valleys-0.753,mid-slope ridges-0.573,canyons and deeply incised streams-0.516),curvature (concave-0.516),aspect (NW-0.511,W-0.288), distance from river (0-150 m class-0.255),and sheepfold density (class 1.54-3.01-0.078).The final GSI map is shown in Figure 5c.The four classes of susceptibility are distributed as follows: very high susceptibility class occupies 20.93% (45.97 km 2 ) of the total area, high susceptibility with 46.50% (102.12 km 2 ), medium susceptibility occupies 26.41% (58 km 2 ), while low susceptibility areas occupy 6.16% (13.52 km 2 ).High and very high classes occupy more than a half of the basin's surface (70.93%).When it comes to the IV method, there are even less significant changes in the final GSI when including sheepfolds as a contributing factor (Figure 5d).

GSI Using IV
The relation of each gully erosion conditioning factor and its spatial relationship with gully triggering is shown in Table 1, column 7. A positive value indicates a high correlation between the class factor and gully erosion phenomena, while negative values highlight a very low correlation with gully erosion phenomena.High values are related to the following conditioning factors: landforms (U-shaped valleys-0.753,mid-slope ridges-0.573,canyons and deeply incised streams-0.516),curvature (concave-0.516),aspect (NW-0.511,W-0.288), distance from river (0-150 m class-0.255),and sheepfold density (class 1.54-3.01-0.078).The final GSI map is shown in Figure 5c.The four classes of susceptibility are distributed as follows: very high susceptibility class occupies 20.93% (45.97 km 2 ) of the total area, high susceptibility with 46.50% (102.12 km 2 ), medium susceptibility occupies 26.41% (58 km 2 ), while low susceptibility areas occupy 6.16% (13.52 km 2 ).High and very high classes occupy more than a half of the basin's surface (70.93%).When it comes to the IV method, there are even less significant changes in the final GSI when including sheepfolds as a contributing factor (Figure 5d).

Validation and Comparison of the GSI Maps
AUC method-for the study area, validation and comparison of the landslide susceptibility maps produced with FR and IV models were checked using the area under curvature (AUC).This method is widely used in studies dealing with gully susceptibility evaluation and creates success rate and prediction rate curves.The success rate curves for the two methods are shown as follows FR (Figure 6a) and IV (Figure 6c); it can be observed that the IV model has a higher area under the curve (AUC) value than FR, which shows the fact that this method is more robust.Regarding the prediction rate of the two models, the IV model has a higher value (AUC = 0.89) (Figure 6d) than the FR model (AUC = 0.87) (Figure 6c).The final GSI maps were classified into four classes: low, medium, high, and very high.

Cultural Heritage Affected by Gully Erosion
According to previous studies from the area [35], for the 107 Neolithic sites, a buffer zone of 100 m was performed; this was done in order to obtain the average surface of a Neolithic site-3 ha.Out of the total number of 107 Neolithic sites, 14 sites are found in areas with very high susceptibility, 58 sites in areas with high susceptibility, 28 sites in medium susceptibility areas, and 7 in low susceptibility areas; this indicates the fact that more than a half of the sites (67%) are located in areas with high and very high susceptibility to gully erosion.Therefore, they are in danger in the future.12 Neolithic sites (Table 2) are part of the List of Historical Monuments (LMI) and National Archaeological Registry (RAN); this means that they have a significant historical value and

Cultural Heritage Affected by Gully Erosion
According to previous studies from the area [35], for the 107 Neolithic sites, a buffer zone of 100 m was performed; this was done in order to obtain the average surface of a Neolithic site-3 ha.Out of the total number of 107 Neolithic sites, 14 sites are found in areas with very high susceptibility, 58 sites in Water 2018, 10, 1077 areas with high susceptibility, 28 sites in medium susceptibility areas, and 7 in low susceptibility areas; this indicates the fact that more than a half of the sites (67%) are located in areas with high and very high susceptibility to gully erosion.Therefore, they are in danger in the future.12 Neolithic sites (Table 2) are part of the List of Historical Monuments (LMI) and National Archaeological Registry (RAN); this means that they have a significant historical value and systematic archaeological excavations have been undertaken.Following the excavations, not only very important data was discovered regarding the periodisation of Cucuteni culture; it offered many insights over the way the prehistoric people were building their settlements, what they were eating, how they were keeping the seeds, and so on.However, if no immediate mitigation measures are being taken, a very significant part of the Neolithic sites (some of them of European significance) from the northeastern part of Romania will disappear.Some mitigation measures have been undertaken [14] within our study area by installing concrete thresholds and planting trees in 2014.Following our last field surveys from May 2018, the gully head advance has been considerably reduced (Figure 7a), and presently, the trees formed a well-developed young forest (Figure 7b); all these measures, along with local authorities understanding of the erosion process, have improved and reduced the erosion process.However, if no immediate mitigation measures are being taken, a very significant part of the Neolithic sites (some of them of European significance) from the northeastern part of Romania will disappear.Some mitigation measures have been undertaken [14] within our study area by installing concrete thresholds and planting trees in 2014.Following our last field surveys from May 2018, the gully head advance has been considerably reduced (Figure 7a), and presently, the trees formed a well-developed young forest (Figure 7b); all these measures, along with local authorities understanding of the erosion process, have improved and reduced the erosion process.

Discussion
When it comes to analysing the data with the sheepfolds included as a contributing factor, there is no big difference between the curves.Whether or not overgrazing is added as a contributing factor in a GSI calculation does not guarantee the fact that the results will be better.However, this does not

Discussion
When it comes to analysing the data with the sheepfolds included as a contributing factor, there is no big difference between the curves.Whether or not overgrazing is added as a contributing factor in a GSI calculation does not guarantee the fact that the results will be better.However, this does not prove the fact that overgrazing is not influencing or accelerating soil erosion.A multi-temporal analysis of soil erosion is needed; this represents a future working hypothesis.Besides this, in our study, we took into account only the sheepfolds, as they have more or less a permanent status.We could not evaluate the overgrazing from cows, for example, because cows are being taken each morning and brought back to their owner each afternoon during the grazing season (which stretches from April-May until October), which could vary depending on the meteorological conditions.

Conclusions
Soil erosion by water, especially gully erosion, is a serious geo-environmental problem at a global scale causing degradation and damages to agricultural lands, and sedimentation of water reservoirs; Romania makes no exception, particularly in the northeastern part of the country.In order to mitigate the effects of soil erosion on Neolithic cultural heritage sites from Bahluiet , catchment, mapping gully erosion susceptibility is a crucial step.Two different statistical models (FR and IV) were employed to create GSI maps for Bahluiet , river basin, from the northeastern part of Romania; the model was prepared using a total of ten conditioning factors and a complete gully inventory.The areas affected by gully erosion were randomly selected into training data (542 gullies) to build the model and testing data (135 gullies) in order to test the predictive model's performance.The weight of each factor was calculated for the two selected statistical models and summed into the final GSI maps.The validation results by AUC method shows that the area under the curve for FR and IV models are 0.79 and 0.78 (without using sheepfolds and using sheepfolds, respectively) and 0.87 and 0.86 (without using sheepfolds and using sheepfolds, respectively), respectively, with a prediction accuracy of 0.79 and 0.78 (without using sheepfolds and using sheepfolds, respectively) and 0.89 and 0.87 (without using sheepfolds and using sheepfolds), respectively.Analysing the GSI maps, it can be observed that approximately 70% of the sites are located in areas with high and very high susceptibility to gully erosion; this means that the Neolithic sites from the study area are in real danger.Another use of the GSI maps is for local authorities and stakeholders to plan the economic activities (e.g., future infrastructure projects), minimise damages costs, environmental and cultural heritage protection, better management of water resources, and hazard mitigation.
Funding: This research received no external funding.

Figure 1 .
Figure 1.Sand columns created inside the Băiceni-Cucuteni gully (Romania), after the collapse from May 2013 (nowadays they do not exist anymore).

Figure 1 .
Figure 1.Sand columns created inside the Băiceni-Cucuteni gully (Romania), after the collapse from May 2013 (nowadays they do not exist anymore).

Figure 2 .
Figure 2. (a) Location of the study area in Romania; (b) Location of training and testing gully dataset and of Neolithic sites.

Figure 2 .
Figure 2. (a) Location of the study area in Romania; (b) Location of training and testing gully dataset and of Neolithic sites.

Figure 3 .
Figure 3. Examples of gully erosion from the study area.

Figure 3 .
Figure 3. Examples of gully erosion from the study area.

Figure 5 .
Figure 5. Gully susceptibility index produced using (a) frequency ratio (FR) method, (b) FR method including sheepfold density as a contributing factor, (c) information value (IV) method, and (d) IV method including sheepfold density as a contributing factor.

Water 2018 , 16 Figure 6 .
Figure 6.(a) Success rate curves produced with FR method, (b) prediction rate curves produced with FR method, (c) success rate curves produced with IV method, and (d) prediction rate curves produced with IV method.AUC-area under the curve.

Figure 6 .
Figure 6.(a) Success rate curves produced with FR method, (b) prediction rate curves produced with FR method, (c) success rate curves produced with IV method, and (d) prediction rate curves produced with IV method.AUC-area under the curve.

Figure 7 .
Figure 7. (a) Băiceni-Cucuteni gully evolution over the last 39 years, (b) the forest developed on the upper part of the gully.

Figure 7 .
Figure 7. (a) Băiceni-Cucuteni gully evolution over the last 39 years, (b) the forest developed on the upper part of the gully.

Table 1 .
Statistical index values for the conditioning factors.

Table 2 .
Neolithic sites listed in the List of Historical Monuments (LMI) and National Archaeological Registry (RAN).GSI-gully susceptibility index; IV-information value.

Table 2 .
Neolithic sites listed in the List of Historical Monuments (LMI) and National Archaeological Registry (RAN).GSI-gully susceptibility index; IV-information value.