Impact of Land Use Changes on the Erosion Processes of a Degraded Rural Landscape: An Analysis Based on High-Resolution DEMs, Historical Images, and Soil Erosion Models

: Soil erosion is one of the major natural risk factors for developing high-value crops and an accurate estimation of spatial distribution and rates of soil degradation can be crucial to prevent crop degradation. In this paper, we use comparisons between high-resolution DEMs and soil erosion models to uncover the short-term landscape evolution of hazelnut crop yields, which are affected by incipient processes of rill development. Maps of rill initiation and evolution were extracted from the analysis of UAV-based multitemporal DEMs and the application of soil erosion models. A comparison between such a short-term analysis and historical orthophotos was carried out. Such a comparison shows how the USPED model predicts, very reliably, where linear erosion occurred. In fact, a reliable overlay between the linear erosive forms predicted by the USPED model and those captured by the UAV images can be observed. Furthermore, land use changes from 1974 to 2020 are characterized by a transition from abandoned areas (1974) to areas with high-value cultivation (2020), which has a strong impact on the spatial distribution of erosion processes and landslide occurrence. Such data represent a key tool for both the investigation of the spatial distribution of hot-spots of soil degradation and the identiﬁcation of effective mitigation practices of soil conservation.


Introduction
Soil erosion is one of the most prominent natural risks in the rural landscapes of Mediterranean areas due to peculiar geological, geomorphological, and climate characteristics such as high relief, lithological features, and the rainfall regime [1]. The abandonment of rural areas has a strong influence on the spatial distribution of soil erosion processes since it can drive the relevant changes in several factors controlling geomorphological processes such as land cover, vegetation recovery, human-induced topographic modifications, and water runoff and infiltration [2]. Many works dealing with the impact of land abandonment have provided contrasting results, since the response to land abandonment in agricultural areas can increase or decrease the soil loss depending on several local factors such as climate setting, topographic features, the rainfall regime, and the type and magnitude of land cover changes/vegetation recovery [2][3][4]. In the arid and semi-arid landscapes of the Mediterranean region, a general trend of increase in soil loss can be observed [5][6][7]. Such an increase in soil loss has been frequently ascribed to natural and anthropic factors such as land abandonment and climate change [6,7]. An increase in the soil loss rate represents a critical issue for landscape management, the mitigation of geological risks, and land degradation. Soil loss rates are higher than the soil formation rates in the Mediterranean basin [8]. Mediterranean soils are easily eroded due to different characteristics: marked relief, 45% of the area having a slope greater than 8%; a high frequency of heavy rains in autumn and winter; poor, shallow, and skeletal soils; and sparse natural vegetation linked to severe summer droughts [9]. Several works have documented the strong impact of the abandonment of farming areas in the Mediterranean region on the hillslope degradation and the increase in the rates of erosion and landslide processes [10][11][12][13][14][15]. Land degradation and badland development is a common problem in the foredeep area of the southern Italian Apennine. This is as a consequence of human and natural factors such as climate setting, deforestation, widespread outcrops of silty clay, and high relief [10,[16][17][18][19].
This work is part of a wider project named "Corylus, Sustainable hazelnut cultivation in Basilicata region" funded by the European Union as part of the program "Rural development 2014-2020". The main objective of the project is to investigate the impact and sustainability of the intensive cultivation of hazelnuts on the rural landscapes of the Basilicata (Figure 1). The hazelnut crop was planted in 2018 and its possible impact on erosion processes has been verified through the visual inspection of historical images and UAV surveys. in the Mediterranean basin [8]. Mediterranean soils are easily eroded due to different characteristics: marked relief, 45% of the area having a slope greater than 8%; a high frequency of heavy rains in autumn and winter; poor, shallow, and skeletal soils; and sparse natural vegetation linked to severe summer droughts [9]. Several works have documented the strong impact of the abandonment of farming areas in the Mediterranean region on the hillslope degradation and the increase in the rates of erosion and landslide processes [10][11][12][13][14][15]. Land degradation and badland development is a common problem in the foredeep area of the southern Italian Apennine. This is as a consequence of human and natural factors such as climate setting, deforestation, widespread outcrops of silty clay, and high relief [10,[16][17][18][19]. This work is part of a wider project named "Corylus, Sustainable hazelnut cultivation in Basilicata region" funded by the European Union as part of the program "Rural development 2014-2020". The main objective of the project is to investigate the impact and sustainability of the intensive cultivation of hazelnuts on the rural landscapes of the Basilicata (Figure 1). The hazelnut crop was planted in 2018 and its possible impact on erosion processes has been verified through the visual inspection of historical images and UAV surveys. The paths of the hazelnut tree in historical times were varied and do not allow a certain attribution of its origins. The movement of the hazelnut species from east to west, in the presence of a plant which has existed since prehistoric times, was widespread in classical times, as confirmed by the scientific term, which is composed of a Latinized Greek flanked by a Latin term (Corylus avellana). The coexistence of two names with different roots is a precise metaphor for the widespread diffusion of the hazelnut tree in the classical world, as evidenced by both historical sources [20,21] and recent studies [22]. Today, the hazelnut is a very ancient plant with a growing market demand and a high degree of The paths of the hazelnut tree in historical times were varied and do not allow a certain attribution of its origins. The movement of the hazelnut species from east to west, in the presence of a plant which has existed since prehistoric times, was widespread in classical times, as confirmed by the scientific term, which is composed of a Latinized Greek flanked by a Latin term (Corylus avellana). The coexistence of two names with different roots is a precise metaphor for the widespread diffusion of the hazelnut tree in the classical world, as evidenced by both historical sources [20,21] and recent studies [22]. Today, the hazelnut is a very ancient plant with a growing market demand and a high degree of adaption to pedological and morpho-climate conditions. Moreover, an analysis of historical sources highlights that hazelnut cultivation is widespread in southern Italy, mainly in landscape sectors affected by land abandonment and degradation. For these reasons, the local government authority and an international food company started an extensive project of planting hazelnut crops in different rural and degraded areas of southern Italy with Land 2021, 10, 673 3 of 18 the aim of: (i) increasing the social and economic growth of internal areas; (ii) promoting good practices of land management and mitigating the impact of land abandonment on the hillslope degradation of steep landscapes.
Extensive farming of internal areas with high-value cropland can be an effective solution to land abandonment and the related deterioration of the rural landscape [23] but these mitigation actions should be integrated by detailed monitoring and investigation of soil loss and land vulnerability estimations. Soil erosion models have been widely used as effective tools to assess the spatial distribution and the rates of erosion and deposition processes, but their calibration/validation is a critical issue [24,25]. The source to sink approach is one of the most used to evaluate the prediction ability of soil erosion models and it has been largely applied worldwide [25][26][27]. Alternative methods of evaluating the robustness of the prediction ability of soil erosion models can take advantage of the rapid increase in tools that are useful for the generation of multitemporal high-resolution DEMs, such as lidar surveys, UAV photogrammetry, and GPS surveys [28][29][30]. Among them, unmanned aerial vehicles (UAVs) have recently become an effective tool to obtain detailed information on the spatial and temporal distribution of different types of geomorphological processes such as gully erosion [31][32][33], fluvial dynamics [34], shoreline retreat [35], and mass wasting [36]. UAV-based photogrammetry, which generates digital elevation models (DEMs) from photo alignments, is a remarkable alternative to soil change measurements [33,37]. It offers many benefits, including a high data accuracy and the avoidance of any impact to the investigated surface. The use of UAV photogrammetry was made possible thanks to digital photogrammetry and computer vision that are part of the Structure from Motion (SfM) technique. Through the use of software based on SfM, it is possible to make 3D reconstructions of models starting from multiple images [38].
In this work, we estimated the spatial distribution of erosion and deposition through the application of the Unit Stream Power Erosion Deposition (USPED) model [39] in a catchment area in the eastern sector of the Basilicata region, southern Italy, (Figure 1), where a high-value crop has recently been planted. Validation of the model results has been carried out by comparing the USPED map with high-resolution DEMs deriving from the UAV workflow. After the preliminary evaluation of the robustness of the model results, we investigated the impact of the land use change (i.e., planting of hazelnut cropland) on the distribution and rates of the soil erosion processes.

Methods
The workflow of the study is summarized in Figure 2. After a detailed analysis of the geological, geomorphological, climate, and land use features of the study area based on classical and consolidated approaches such as field surveys, photo-interpretations, and DEM analysis, we prepared the basic parameter maps of the USPED model. More specifically, the following data were used to calculate the model parameters:  Finally, the results derived from the application of the USPED erosion model were compared with the photo-interpretations of the historical images and the short-term analysis of the erosion/deposition processes derived by the UAV-based high-resolution DEM. The UAV DEM was reconstructed as follows: DJI Phantom 3 Std UAV (Dà Jiāng Innovations, Shenzhen, China) was applied to take the photos, Leica Viva GS14 GNSS (Leica Geosystems, Heerbrugg, Switzerland) was applied to acquire the GCPs and a workstation was used for data processing.

USPED Model
The Unit Stream Power Erosion and Deposition (USPED) model estimates the average soil loss (A) using the basic structure of the RUSLE empirical equation [42]: where, A (Mg ha −1 yr −1 ) is the annual average soil loss, R (MJ mm h −1 ha −1 yr −1 ) is the rainfall intensity factor, K (Mg h MJ −1 mm −1 ) is the soil erodibility factor, LS is the topographic factor, C (dimensionless) is the land cover factor and P (dimensionless) is the soil conservation or prevention practice factors [42]. Finally, the results derived from the application of the USPED erosion model were compared with the photo-interpretations of the historical images and the short-term analysis of the erosion/deposition processes derived by the UAV-based high-resolution DEM. The UAV DEM was reconstructed as follows: DJI Phantom 3 Std UAV (Dà Jiāng Innovations, Shenzhen, China) was applied to take the photos, Leica Viva GS14 GNSS (Leica Geosystems, Heerbrugg, Switzerland) was applied to acquire the GCPs and a workstation was used for data processing.

USPED Model
The Unit Stream Power Erosion and Deposition (USPED) model estimates the average soil loss (A) using the basic structure of the RUSLE empirical equation [42]: is the rainfall intensity factor, K (Mg h MJ −1 mm −1 ) is the soil erodibility factor, LS is the topographic factor, C (dimensionless) is the land cover factor and P (dimensionless) is the soil conservation or prevention practice factors [42].
The USPED model adopts an empirical equation for net erosion and deposition in a transport capacity limited regime, which differs from the RUSLE for the computation of the LS-factor. The equation, based on the upslope contributing area, has been developed by the authors of [39] to include topographic complexity by considering both the profile (in the downhill direction) and the tangential (perpendicular to the downhill direction) curvature [43,44]. Net erosion or deposition within a grid cell (ED) are computed as the divergence of sediment flow (change in sediment transport capacity) using the following equation [39,43]: where, α is the aspect of the terrain (in degrees). The physical basic assumption of the model is that net erosion pixels coincide with areas of profile convexity and tangential concavity (flow acceleration and convergence) while net deposition areas correspond to areas of profile concavity (decreasing flow velocity). The widespread use of semi-empirical erosion models such as RUSLE and USPED on a large scale is due to their relatively good degree of prediction ability with easily accessible data. As a matter of fact, many applications in different landscapes showing accurate results adopt simplified input parameters, which can be extracted by regional-scale and easily available data ( [1,17,26,36,45], among others). A similar approach was used in this work, where detailed soil samples, small-scale plots of vegetation/soil features, and rainfall data at a high temporal resolution were not available. For these reasons, the spatial variations of the C-and K-factors were derived from land cover and lithological maps at a 1:5000 scale, whereas the R-factor was extracted from an empirical equation that correlated daily rainfall data with rainfall erosivity [46]. The R-factor depends on rainfall-runoff characteristics, which in turn are influenced by geographic location and altitude. In the original formulation of the USLE, the R-factor was calculated as the product of the total kinetic energy of the storm (E) and its maximum 30-min intensity (I30). Since this information was not available in the study area, we used the power-law equation proposed by the authors of [46] for the Basilicata region. This empirical equation was statistically derived from 20-min and hourly precipitation data higher than 10.0 mm. The formula is: where, EI30 is the rainfall erosivity in MJ mm h −1 ha −1 yr −1 and P 24 the daily rainfall amount in mm. The annual averaged R-factor value was estimated for the period from January 2010 to December 2020.
The soil erodibility factor (K) is the rate of soil loss per rainfall erosion index unit, Mg h MJ −1 mm −1 , as estimated on a standard Wischmeier erosion plot [47]. The K-factor represents the soil profile response to the erosive power of rainfall events and was estimated using the following equation that includes soil texture and permeability: where, M is the organic matter content (%), Si is the silt fraction (%), 2 to 50 µm, fS is the fine sand content (%) 50 µ to 100 µm, C is the clay content (%) less than 2 µm, S is the sand content (%) 50 µm to 2 mm, A is the structure, and P is the permeability class (within the top 0.60-0.70 m) and the factor 0.001317 is derived from the division by 100 of the conversion value (0.1317) to the Si.
Texture, structure, and permeability data for each soil unit of the study area were extracted from previous works that have dealt with the physical and chemical characterization of deposits outcropping in the foredeep area of the southern Apennine chain [17,[48][49][50]. The K-factor map was drawn by combining lithological and land use maps with the abovementioned literature data about the texture and permeability of the main lithological units of the study area.
The topographic parameters, slope, and specific catchment area were derived from a 1 m DEM, deriving from a LIDAR survey acquired in 2013.
The modified LS-factor is given by the following relation: where, Ac (m) is the upslope contributing area per unit of width; β is the slope angle in degree; α0 is the length (72.6 ft, equal to 22.13 m) and β0 is the angle of the standard terrain  [39,43]. The C-factor (dimensionless), which reflects the effects of cropping and management practices on soil erosion rates, was calculated according to the literature data [51,52], with values assigned on the basis of the vegetative cover, its density, and the monthly rainfall runoff erosivity.
The information on land use was extracted from a detailed (i.e., 1:5000 scale) inventory of land use/vegetation cover provided by the Regional Spatial Data Infrastructure of the Basilicata Regional Authority (http://rsdi.regione.basilicata.it, accessed on 25 June 2021). The classification follows the procedures and recommendations of the Corine Land Cover project [41]. The C-values were assigned to the resulting 7 land use classes according to a revision of the C-values proposed in the literature (see for example [53]).

Photogrammetric Data of UAV-Based DEMs
In this work, the UAV DJI Phantom 3 Standard (quadricopter) was used to acquire the photos ( Table 1). The drone was equipped with a stabilized camera, with a sensor able to compensate for involuntary movements due to the wind. In this way, the acquired images were guaranteed to have the correct orientation with respect to the ground. To ensure a high image resolution, the drone flew at a low elevation above the ground level (about 50 m). The drone recorded the coordinates for each photo through an internal GPS. The features of the camera (Table 1) represent the essential factors for a good flight plan in addition to the resolution and the size of the sensor, especially the camera viewing angle (FOV-Field of View). In the flight plan, we selected the FOV value, the flight altitude, and the overlap percentage between adjacent swipes (side overlap). We calculated the shutter speed value to have a good percentage of longitudinal overlap (frontal overlap). We set a 60% side overlap and a frontal overlap of 80%. With these characteristics, the study area was covered with about 1100 photos. Agisoft Metashape (Agisoft LLC., St. Peterburg, Russia) was used to generate 3D models after the acquisition of the images ("Software environment" in Step 5 UAV Photogrammetry- Figure 2). The SfM algorithm calculated the position of the photos and the geometry of the scene, comparing all the images and finding the common points [54]. This process involved different steps. In the first step, homologous image points (Tie points) between the images were matched. In this way, the images were oriented in the space and a cloud of low-density points, namely, the dense sparse cloud, was created. Starting from this sparse cloud, a dense point cloud was built in the next step. In the last step, the DEM and orthophoto of the study area were built starting from the dense point cloud. The orientation of the 3D model and the generation of the DEM with a spatial resolution of 3.5 cm were drawn using 13 ground control points (GCPs) acquired by a RTK GPS.

Geological and Climate Setting
The study area was a small catchment area located in the foredeep of the southern Italian Apennine (Figure 3).
images and finding the common points [54]. This process involved different steps. In the first step, homologous image points (Tie points) between the images were matched. In this way, the images were oriented in the space and a cloud of low-density points, namely, the dense sparse cloud, was created. Starting from this sparse cloud, a dense point cloud was built in the next step. In the last step, the DEM and orthophoto of the study area were built starting from the dense point cloud. The orientation of the 3D model and the generation of the DEM with a spatial resolution of 3.5 cm were drawn using 13 ground control points (GCPs) acquired by a RTK GPS.

Geological and Climate Setting
The study area was a small catchment area located in the foredeep of the southern Italian Apennine (Figure 3).  This sector of the chain is a hilly landscape characterized by widespread outcrops of a thick (i.e., several tens of meters) succession of marine grey-blue silty clays Early Pleistocene in age [55,56]. These deposits crop out in the central and northern sectors of the study area and are unconformably covered by marine to transitional coarser-grained clastic wedges, Early-Middle Pleistocene in age. They are composed of sand and gravel deposits and record the progressive emersion of the foredeep basin in response to the regional uplift of the southern Apennine chain [55][56][57]. Coarse-grain deposits are located in the headwater sectors of the catchment (Figure 3).
From a climate viewpoint, the study area exhibits a typical Mediterranean climate.
Daily rainfall records from the Irsina station in the last two decades highlight a mean annual value of total rainfall of 578 mm, with a strong seasonal variability and a prevalent concentration of rainfall in autumn and spring (Figure 4). Campomaggiore unit; (8) Mesozoic-Cenozoic shallow-water carbonates of the Apulian platform; (9) thrust front of the chain; (10) volcanoes. This sector of the chain is a hilly landscape characterized by widespread outcrops of a thick (i.e., several tens of meters) succession of marine grey-blue silty clays Early Pleistocene in age [55,56]. These deposits crop out in the central and northern sectors of the study area and are unconformably covered by marine to transitional coarser-grained clastic wedges, Early-Middle Pleistocene in age. They are composed of sand and gravel deposits and record the progressive emersion of the foredeep basin in response to the regional uplift of the southern Apennine chain [55][56][57]. Coarse-grain deposits are located in the headwater sectors of the catchment (Figure 3).
From a climate viewpoint, the study area exhibits a typical Mediterranean climate. Daily rainfall records from the Irsina station in the last two decades highlight a mean annual value of total rainfall of 578 mm, with a strong seasonal variability and a prevalent concentration of rainfall in autumn and spring (Figure 4).
In general, the seasonal distribution of rainfall is shaped by dry summers and cold winters, and rainfall peaks occur in autumn. In this season, we observed heavy rainfall events related to one-day or multi-days of severe precipitation. An analysis of the recent trends in rainfall data suggest a general decrease in the total annual rainfall and an increase in both dry periods and short-term extreme events (Figure 4). In general, the seasonal distribution of rainfall is shaped by dry summers and cold winters, and rainfall peaks occur in autumn. In this season, we observed heavy rainfall events related to one-day or multi-days of severe precipitation. An analysis of the recent trends in rainfall data suggest a general decrease in the total annual rainfall and an increase in both dry periods and short-term extreme events (Figure 4).

Geomorphological Analysis
The study area included two catchment areas with elongated shapes, ranging in altitude from 500 to 210 m a.s.l. It covers a total area of about 7 kmq and is located along a Land 2021, 10, 673 9 of 18 SN-trending slope developing along the right-orographic side of the Basentello River. The studied catchment area is mainly carved in grey-blue silty clays and shows a north-dipping gentle slope (i.e., a mean slope of about 10 • ), which is transversally incised by a dendritic drainage network (Figure 1). Higher-altitude sectors of the study area are mainly carved in gravel deposits, thus exhibiting steeper slopes and higher relief. Channel incision, rills, and gullies developments are the main erosion processes of the study area, although small and shallow earthflows can be observed in the central and northern sectors of the catchment areas, where clay-rich lithological units crop out.
The present-day land use of the study area was derived from a digital inventory of the Basilicata Region (http://rsdi.regione.basilicata.it, accessed on 25 June 2021). It roughly follows the III level of the Corine Land Cover project [56] and highlights the prevalence of arable lands, which cover about 80% of the total catchment area. Natural grasslands and sclerophyllous vegetation represent the other relevant land use classes of the study area. They show a high level of degradation and are frequently affected by mass movement processes and accelerated linear erosion. Hazelnut cropland was also found in the catchment area.

USPED Model
The mean annual average erosivity factor computed for the study area using the hourly rainfall data of the Irsina rain gauge is equal to 664 MJ mm ha −1 h −1 yr −1 for the 2000-2020 period.
The land use and land management factor (C-factor) in the drainage basin ( Figure 5A) varies between a maximum value of 0.35 assigned to bare or degraded areas and a minimum value of 0.001 related to urban areas, roads, and small areas of broad-leaved forest. Intermediate values ranging from 0.1 to 0.25 were assigned to permanent crops (olive groves and hazelnut croplands), pastures, and arable lands. Lower values of the C-factor ranging from 0.025 to 0.035 were attributed to schlerophyllous vegetation and transitional woodland/sparse scrubs, respectively ( Figure 5A).
The soil erodibility factor (K-factor) ( Figure 5B) followed the stratigraphic contacts of the lithological map and ranged in the drainage basin between 0.02 Mg h MJ −1 mm −1 and 0.033 Mg h MJ −1 mm −1 . The fine-grained clastic deposits outcropping in the intermediate and lower sectors of the catchments are inorganic silty clay, typically defined by a percentage of silt higher than 55-65% [50]. Coarser-deposits of the study area were composed of more sandy fractions and by a higher permeability, thus exhibiting higher values of the K-factor according to the Wischmeier equation. The assigned values of the K-factor are in accordance with similar estimations in adjacent areas of southern Italy (see for example [17]).
The terrain analysis showed a prevalence of areas with moderate slopes characterized by an average inclination lower than 12 • . Accordingly, the average topographic LS-factor of the watershed ( Figure 5C) is low, with an average of 1.74 and a standard deviation of σ = 22.6. Minimum LS-factor values occur in flat areas such as drainage divides.
The results of the application of the USPED model are shown in Figure 6. Each pixel has been classified using the following eleven classes of erosion/deposition: Land 2021, 10, 673 10 of 19   The frequency distribution of the USPED model classes highlights that stable and low erosion/deposition areas of the study area cover about 84% of the entire catchment area, corresponding to a surface of about 6.1 kmq (Figure 7). The model results also indicate the diffuse presence of areas affected by extreme erosion and deposition, mainly located along the main channels of the drainage network. For example, the erosion surface area with values higher than 10 Mg ha −1 yr −1 is 0.23 kmq, which corresponds to about 3.1% of the whole catchment. The frequency distribution of the USPED model classes highlights that stable and low erosion/deposition areas of the study area cover about 84% of the entire catchment area, corresponding to a surface of about 6.1 kmq (Figure 7). The model results also indicate the diffuse presence of areas affected by extreme erosion and deposition, mainly located along the main channels of the drainage network. For example, the erosion surface area with values higher than 10 Mg ha −1 yr −1 is 0.23 kmq, which corresponds to about 3.1% of the whole catchment.

UAV Data and Landscape Changes
About 1060 photos were taken by UAVs at an altitude of 50 m. The results and timing of data processing, through the use of SfM, are summarized in Table 2. The resolution of the DEM and of the orthophoto obtained (Table 2) underlines the high accuracy that can be obtained using this methodology. In fact, the errors of the 3D

UAV Data and Landscape Changes
About 1060 photos were taken by UAVs at an altitude of 50 m. The results and timing of data processing, through the use of SfM, are summarized in Table 2. The resolution of the DEM and of the orthophoto obtained (Table 2) underlines the high accuracy that can be obtained using this methodology. In fact, the errors of the 3D model on the Ground Control Points (GCP) resulted in: East 0.061 m; North 0.047 m; Altitude 0.092 m.
In the cultivated hazelnut area, linear erosive forms were observed ( Figure 8A), which can be also observed in the USPED model map ( Figure 8B). In the cultivated hazelnut area, linear erosive forms were observed ( Figure 8A), which can be also observed in the USPED model map ( Figure 8B). Within the study area, two changes were observed in the investigated time period. In fact, when both changes are compared to 1974, it can be seen that the area occupied by the cultivation of hazelnuts in 2020 was abandoned in 1974. The abandonment is noticeable by the pronounced linear erosive forms and by the presence of small gravitational movements.
The other element that has undergone changes compared to 1974 is the position of the road to reach the hazelnut area ( Figure 10).  Within the study area, two changes were observed in the investigated time period. In fact, when both changes are compared to 1974, it can be seen that the area occupied by the The other element that has undergone changes compared to 1974 is the position of the road to reach the hazelnut area ( Figure 10).  Before 2017 (Figure 9), the study area was mainly used as arable land; the only change that can be seen compared to 2017 is the presence of hazelnut cultivation. Before 2017 (Figure 9), the study area was mainly used as arable land; the only change that can be seen compared to 2017 is the presence of hazelnut cultivation.

Discussion and Concluding Remarks
Many works have independently used soil erosion models and a multi-temporal analysis of aerial photos to investigate the spatial distribution and magnitude of geomorphological processes and the related control factors such as external forcing or anthropogenic impact [24,28,33,36,[58][59][60][61][62][63]. Although land degradation and rural abandonment are frequently ascribed as relevant factors related to an increase in the occurrence and rates of linear and slope erosion processes [10,11,64,65], few studies have documented the role of cropland on the mitigation of geomorphological processes. In this work, we combined the USPED erosion model, a DEM derived from UAV photos, and an analysis of land use change from 1974 to 2020 to investigate the long-and short-term geomorphological processes of a rural landscape. These data were used to evaluate historical land use changes and their impact on the short-and long-term slope and fluvial processes of the study area. More specifically, we investigated the possible role of localized farming in the catchment area with high-value cropland as a sustainable action to prevent issues of land degradation and deterioration.
In foredeep areas of Basilicata, landscape remodeling activities and the use of machinery are typical agricultural practices in rural landscapes. Such practices are mainly related to cereal cultivation and several works have demonstrated their negative effect on soil erosion, landslide occurrence, badland development and, more generally, land degradation [17,23,50]. Land use changes from an abandoned area to widespread cultivation of Land 2021, 10, 673 15 of 18 cereals can then promote an increase in soil erosion and land degradation. On the contrary, the planting of more sustainable cropland such as hazelnut plants can represent a key strategy to prevent soil erosion and land degradation [23].
Our data are substantially based on a visual inspection of the multitemporal dataset and suggest a positive effect of farming practices (i.e., hazelnut cultivation) on the mitigation of slope processes and geomorphological processes. As a matter of fact, land use of the studied catchment area underwent several changes between 1974 and 2020. It was observed that in 1974 ( Figure 9A), the central area was mostly disturbed by gravitational processes, which were not present in the images of the years 2013, 2017, and 2020 ( Figure 9B-D). This difference is probably due to two factors that have changed over time: (i) change in land use; and (ii) change of the position of the road present in the western part of the images (Figure 10).
A short-term analysis of linear erosion processes based on the USPED model and its comparison with UAV images highlights that rill initiation and development is the main geomorphological process acting on the hazelnut cropland. In general, the comparison between the USPED erosion map and the DEM obtained from the SfM workflow demonstrates the reliability of the USPED in predicting the formation of linear erosion forms. Our results show that UAV images and the related data (i.e., ortophotos and high-resolution DEMs) can contribute to the identification of soil degradation with a similar ability to the soil erosion model. Such an integrated approach can represent a key-point in combining the validation and calibration of the two methods. The USPED model shows a diffuse and fairly superficial linear erosion ( Figure 6). During the generation of the model parameters, the presence of the hazelnut area was inserted in the C-factor. As shown in Figure 8, the presence of the cultivated hazelnut favored the surface runoff following the rows of the hazelnut grove. The model predicted rills between the alignment of the hazelnut plant, which can be also observed by analyzing the UAV-based DEM. A relevant control factor for the spatial distribution of erosion processes is the position of the road (Figure 10): the modification of its location after the land use change from natural areas to agricultural land promoted the modification of the hydraulic conditions of this sector of the catchment area. In fact, linear elements such as trails and roads act as preferential routes along a slope for water flow [66]. In this way, the quantity of water received by the slope by runoff is reduced ( Figure 10). The energy with which the water erodes is also reduced. This aspect is also highlighted by the USPED model ( Figure 10D) and can be useful to drive the possible mitigation actions of erosion processes. Another important factor when calculating plantation erosion rates is the age of cultivation. In our case, the hazelnut plants were three years old, with an average height of 1.5 m and a poorly developed root system. The erosion rate of young plantations is higher than for older plantations, as demonstrated by Rodrigo Comino et al. (2017) [66] in vineyard plantations.
The study of land use dynamics and its impact on changes in soil properties and surface runoff will help toward the proper planning, management, and sustainable development of land resources [67]. The information can be employed to predict potential changes that would likely exist in the future. Furthermore, it can contribute information to allow good practices promoting soil conservation and mitigating land degradation to be undertaken. Our approach has a practical relevance since it can be useful for predicting localized areas of higher vulnerability or accelerated erosion.
In conclusion, the results show how the USPED erosion prediction model works well in areas characterized by steep slopes and poor soil cover. It also highlights how the cultivation of hazelnuts (high value), especially in the first years, can be affected by soil erosion and prevent the development of plants. This factor can have several negative effects, including slowing down the future harvest of hazelnuts with negative social and economic consequences. The cultivation of hazelnuts can also be an effective solution to land abandonment and the related deterioration of rural areas [23]. These mitigation actions should be integrated by detailed monitoring and an investigation of soil loss and land vulnerability estimations at a basin scale and can be applied to other areas with similar types of soil, climate, and socio-economic conditions.