Impact of Conservation Agriculture on Soil Erosion in the Annual Cropland of the Apulia Region (Southern Italy) Based on the RUSLE-GIS-GEE Framework

The processes of soil erosion and land degradation are more rapid in the case of inappropriate agricultural management, which leads to increased soil loss rates. Moreover, climatic conditions are one of the most important determining factors affecting agriculture, especially in the Mediterranean areas featuring irregular rainfall and high summer temperatures. Conservation agriculture (CA) can make a significant contribution to reducing soil erosion risk on the annual cropland (ACL) of the Mediterranean region in comparison with conventional management (CM). The objective of this study is to provide soil loss rate maps and calculate the values for each altitude and slope class and their combination for the Apulia region in four annual production cycles for the scenarios CM and CA. The present study estimates the significance of the adoption of CA on soil erosion assessment at regional scale based on the Revised Universal Soil Loss Equation (RUSLE) model. The parameters of the RUSLE model were estimated by using remote sensing (RS) data. The erosion probability zones were determined through a Geographic Information System (GIS) and Google Earth Engine (GEE) approach. Digital terrain model (DTM) at 8 m, ACL maps of the Apulia region, and rainfall and soil data were used as an input to identify the most erosion-prone areas. Our results show a 7.5% average decrease of soil loss rate during the first period of adoption of the four-year crop cycle—from 2.3 t ha−1 y−1 with CM to 2.1 t ha−1 y−1 with the CA system. CA reduced soil loss rate compared to CM in all classes, from 10.1% in hill class to 14.1% for hill + low slope class. These results can therefore assist in the implementation of effective soil management systems and conservation practices to reduce soil erosion risk in the Apulia region and in the Mediterranean basin


Introduction
Soil erosion is one of the main parameters for assessing soil quality [1][2][3]. It is defined as "the movement and transport of soil by various agents, particularly water, wind, and mass movement" [4]. In addition, soil quality is defined "the capacity of a soil to function, within natural or managed ecosystem boundaries, to sustain plant and animal productivity, to maintain or enhance water and air quality, and support human health and habitation" [5]. Soil erosion and soil quality are strongly correlated phenomena. Soil quality affects the rate of soil loss and is in tum affected by it. Erosion effects on soil quality are determined so far to assess soil erosion, for example, the Pan-European Soil Erosion Risk Assessment (PESERA) model and erosion rates based on runoff plot data [32,33]. However, as pointed out by Panagos et al. [32], such models do not capture the effects of conservation practices and their potential to mitigate soil erosion risk. Thus, even the effectiveness of policies that promote such practices cannot be addressed adequately. On the other hand, this becomes possible by adopting other models, such as the RUSLE, which has been used in the past years to monitor erosion, and which is specific to the EU context [34]. Moreover, soil erosion is one of the agro-environmental indicators adopted by the European Commission services for monitoring agricultural and environmental policies. This model can capture the impact of land-use changes, and thus highlight the efficacy that the European agri-environmental policies would have on restoring soil health [32,34]. Such a model has been used to estimate soil loss rate in several Mediterranean areas [35].
The objective of this study is to provide soil loss rate maps and calculate the averaged values for each altitude and slope class and their combination for the ACL of Apulia region for both the scenarios. This part of Italy suffers from heavy rainfall (fall/winter period) as well as decreasing precipitation. Moreover, increase of temperatures and consequently drier conditions (summer) caused by more and more evident effects of climate change and CM pose an additional threat to Apulian cropland [36,37]. Our approach is based on the RUSLE-GIS-GEE framework [3,38], using more suited databases at regional scale. This could provide greater detail and accuracy in calculating soil loss rate for ACL, separately for the two management systems: CM and CA for the period 2016-2020, following the introduction in 2016 of the specific sub-measure M10.1 "Conservation Agriculture", in which only a part of the farmers participated, being a voluntary measure.

Study Area
The Apulia region ( Figure 1) is situated in the southeastern part of the boot-shaped Italian Peninsula bordering the Adriatic and Ionian seas along the east and southeast coasts, respectively. The region is divided into 257 municipalities grouped into five provinces and covers a surface area of approximately 19,500 km 2 . The region is characterized by low mountains located in the Gargano promontory and in the Daunian Sub-Apennine, respectively, in the north and east of the Foggia province; the Tavoliere plain (the second largest plain in Italy), which extends for 3000 km 2 in the central and southern part of the Foggia province; and the Murgia plateau, which covers a surface of 4000 km 2 between the provinces of Barletta-Andria-Trani and Bari [39]. The Apulia region features hot and dry summer seasons, as well as mild and rainy winter seasons typical of semiarid Mediterranean climate. Annual precipitation varies between 450 and 550 mm in much of the region (two-thirds concentrated from autumn to winter). The highest values of precipitation, with more than 900 mm y −1 , are observed in the Gargano area, in the province of Foggia, whereas the lowest values, around 400 mm y −1 , are observed in the Tavoliere plain. The hydrological regimes are irregular, of torrential type, with high stream and river flow rates during the rainy season and practically no water flow in summer [37]. Agriculture plays a very important role in the economic context of the Italian territory; this region is second in terms of production of olive oil, wine, and fresh vegetables. Particularly important is the production of durum wheat in the Tavoliere plain, orchards (cherries and figs) near Bari, and tobacco in the province of Lecce [40]. Durum wheat and, in general, cereals (barley, rye, oats) are the most important typical ACL in Apulia region. For ACL in Southern Italy, the period between 15 October to 1 August corresponds approximately to the period from sowing to harvesting. in Southern Italy, the period between 15 October to 1 August corresponds approximately to the period from sowing to harvesting.

Soil Loss Modeling: RUSLE Factors
The RUSLE model has been widely used for both agricultural and natural land to estimate annual soil loss rate and to evaluate soil erosion risk [41,42]. This model is accurate, easy to apply, and needs a moderate amount of data. Its usage has increased over the past few decades, particularly with the increase of RS and GIS applications [43,44]. The advent of RS and GIS application has increased the interest in developing new methods of calculation and sharing data, using cloud-computing platforms. GEE was developed as an open-source platform for analyzing geospatial data. GEE has been used worldwide for retrieving and processing many Earth observation data, which nowadays cover all geospatial data needed to build the RUSLE-GIS-GEE framework in a comprehensive and robust cloud-based environment. GEE's capabilities can be used to process large amounts of geospatial data: especially, with improvements in these data's availability and processing time, for this reason, it is successfully used in several fields on both regional and global scales [45][46][47]. In the current study, soil loss rate estimation based on RUSLE was implemented in the GEE environment to increase the ability to determine susceptibility to erosion risk.
The RUSLE model [48,49] was used to estimate the soil loss rate for the scenarios for the CM and CA system in the Apulia region with limited annual cropland for the first period of adoption (2016-2020). RUSLE provides an ideal framework to assess soil loss rate and its factors for both the scenarios. Specifically, RUSLE considers support practices (P), rainfall (R), soil erodibility (K), topography (LS), and cover management (C) as important factors affecting soil loss rate. RUSLE can be mathematically expressed as is the rainfall erosivity factor, K (Mg h MJ −1 mm −1 ) is the soil erodibility factor, LS (unitless)

Soil Loss Modeling: RUSLE Factors
The RUSLE model has been widely used for both agricultural and natural land to estimate annual soil loss rate and to evaluate soil erosion risk [41,42]. This model is accurate, easy to apply, and needs a moderate amount of data. Its usage has increased over the past few decades, particularly with the increase of RS and GIS applications [43,44]. The advent of RS and GIS application has increased the interest in developing new methods of calculation and sharing data, using cloud-computing platforms. GEE was developed as an open-source platform for analyzing geospatial data. GEE has been used worldwide for retrieving and processing many Earth observation data, which nowadays cover all geospatial data needed to build the RUSLE-GIS-GEE framework in a comprehensive and robust cloud-based environment. GEE's capabilities can be used to process large amounts of geospatial data: especially, with improvements in these data's availability and processing time, for this reason, it is successfully used in several fields on both regional and global scales [45][46][47]. In the current study, soil loss rate estimation based on RUSLE was implemented in the GEE environment to increase the ability to determine susceptibility to erosion risk.
The RUSLE model [48,49] was used to estimate the soil loss rate for the scenarios for the CM and CA system in the Apulia region with limited annual cropland for the first period of adoption (2016-2020). RUSLE provides an ideal framework to assess soil loss rate and its factors for both the scenarios. Specifically, RUSLE considers support practices (P), rainfall (R), soil erodibility (K), topography (LS), and cover management (C) as important factors affecting soil loss rate. RUSLE can be mathematically expressed as where A (Mg ha −1 y −1 ) is the longtime average annual soil loss rate, R (MJ mm h −1 ha −1 y −1 ) is the rainfall erosivity factor, K (Mg h MJ −1 mm −1 ) is the soil erodibility factor, LS (unitless) is the slope length factor and slope steepness factor, C (unitless) is the land cover and management factor, and P (unitless) is the soil conservation or prevention practices factor [49].

R-Factor
R-factor is an index of rainfall erosivity that quantifies the potential capacity of rain to cause erosion [49]. The factors that are affected by rainfall erosivity are amount, intensity, terminal velocity, drop size, and drop size distribution of rain [50]. For a given location, it is the long-term average of the annual R aj values which, in turn, are given by the sum of all the erosion index (EI) single-storm EI 30 values for year j [51].
The annual rainfall-runoff erosivity (R-factor) values (MJ mm h −1 ha −1 y −1 ) for Italy were computed by the following equation [51]: (2) where N is an N-year period, and EI is the rainfall erosivity index (MJ mm ha −1 h −1 ), expressed as EI 30−annual = 12.142(abc) 0.6446 (3) where the variables a, b, and c are the annual precipitation, the maximum annual daily precipitation, and the maximum annual hourly precipitation, respectively-all expressed in centimeters. Variable a represents less erosive precipitations, with a cumulative effect over a long period. Variables b and c describe very erosive effects due to extreme rainfalls in storms and heavy showers. This study estimated the R-factor based on ERA5-Land (E5L) gridded weather data, freely available as product of the Copernicus Climate Change Service (C3S). The E5L is a high-resolution reanalysis dataset which covers the period from 1981 to present on a regular grid with a spatial resolution of 0.1 • × 0.1 • latitude-longitude referred to as geographic coordinate system WGS84 (EPSG:4326), corresponding to a horizontal resolution of 9 km. The data are provided with an hourly time-step and released with a delay of 2-3 months from present.
All grid cells (for a total of 231) located in the Apulia region were considered for retrieving total hourly precipitation (estimated in millimeters, mm) for the period from January 1981 to December 2019. The accumulated precipitation values were processed to derive the total hourly precipitation (mm) from 00 UTC to the hour ending at the forecast step.
The R-factor was calculated as average of EI yearly values regarding four N-year periods: 1981-2016, 1981-2017, 1981-2018, and 1981-2019. R-factor descriptive statistics, such as minimum, maximum, standard deviation, and weighted average values, for the Apulia region were estimated, and the results are listed in Table S1. Data processing was performed through R (R Core Team, 2020). For RUSLE calculation (see Section 2.4), the results obtained for R-factor were projected into the EPSG: 3035 reference system and then masked for the ACL 2015 (1981-2016) and for ACL 2018 (1981-2017, 1981-2018, and 1981-2019) in GEE (Table S1).

K-Factor
The K-factor is the soil erodibility factor (t ha h ha −1 MJ −1 mm −1 ), an empirical parameter based on the measurements of specific soil erodibility [52]. This parameter is measured based on these soil properties: texture, organic matter, structure, and permeability of the topsoil profile [53]. In this study, the reference value of K-factor is the one obtained from the "Soil Erodibility in Europe High Resolution dataset" [53] provided by the Joint Research Center (JRC) of European Soil Data Centre (ESDAC) and clipped for the Apulia region by using QGIS. The K-factor is estimated for the 20,000 field sampling points (133 for Apulia region) included in the Land Use/Cover Area frame (Land Use and Coverage Area, LUCAS) survey [54] and then interpolated with a Cubist regression model [55] using spatial covariates such as remotely sensed data and terrain features to produce a 500 m resolution K-factor map of Europe [53].

P-Factor
Support practice (P-factor) is an expression of the effects of agricultural management practices that reduce the erosion potential of runoff by influencing drainage patterns and the concentration and velocity of runoff [52]. The adoption of supporting conservation practices decreases the p value, which ranges between 0.2 (terraces with reverse slope) and 1.0 (no erosion control practices). The average value is estimated at around 0.95 in agricultural land, and for the EU it is estimated at 0.97. This effect is considerably greater in sensitive regions such as the Mediterranean area, although the reduction rate can generally be relatively small when adopting supportive practices [56]. In the present study, the reference value of P-factor is the one obtained from the EU datasets [57] provided by the JRC's ESDAC with 1 km resolution, and then clipped for the Apulia region by using QGIS.

LS-Factor
In the RUSLE model, the L and S factors represent the influence of the terrain topography on the sediment transport capacity of the overland flow [52]. Slope length (L) is defined as the point where overland flow starts to the point in which deposition occurs or runoff waters are channelized [58]. Slope steepness (S) describes how erosion increases with slope angle [58]. The combined LS-factor (dimensionless) describes the potential of surface runoff in accelerating soil loss rate, and, in most studies, determines the spatial resolution (cell size) of the modeled soil loss estimates. The topographic LS-factor was calculated by using the 8 m high-resolution DTM provided by the Apulia region (available on http://www.sit.puglia.it, accessed on 1 November 2021).
The LS-calculation was performed by using the equation proposed by Desmet and Govers [59]: where A i,j−in is the contributing area at the inlet of grid cell (i,j), measured in m 2 . D is the grid cell size (meters), X i,j = sinα i,j + cosα i,j , the α i,j is the aspect direction of the grid cell (i,j). This equation was implemented by using the System for Automated Geoscientific Analyses (SAGA) software in QGIS, which incorporates a multiple flow algorithm and contributes to a precise estimation of flow accumulation [58]. It provides a comprehensive set of modules for data analysis, focusing on DTMs and terrain analysis [60,61]. A multiple flow algorithm present in SAGA (LS-factor field based) allows the calculation of LS-factor [58].

C-Factor
Among the inputs of RUSLE, the cover and management factor (C-factor) is the one most sensitive factor [62] that ranges between 0 and 1. The C-factor follows plant growth and rainfall dynamics [52,63] and can be managed by farmers and managers to control soil erosion in agricultural activities [34]. The C-factor represents the effect of cropping and management practices on soil erosion by water [49]. The decrease of the C-factor can be promoted by changing the amount of vegetation cover and tillage practices and soil management measures (e.g., reduced or no tillage and cover crop residues) that protect the soil surface, disperse raindrop energy, and reduce surface runoff [34,49]. Land-use types, crop rotation, and cultivation and management practices show obvious spatial and temporal variations that affect the accuracy of the C-factor estimate, ultimately affecting soil loss rate estimated by RUSLE [64,65]. Therefore, it is necessary to improve the ability to capture the space-time dynamics of the C-factor. In our study, to estimate the C-factor within the ACL of the Apulia region, we started from the equation proposed by Panagos et al. [34]: where C crop is the C-factor based on the crop composition of an agricultural area, and C management quantifies the influence of management practices (reduced tillage, cover crop, and crop residues) on soil loss rate reduction. With regard to C management , the combined effect of tillage practice (C tillage ), plant residues (C residues ), and cover crops (C cover ) was taken into account for the estimation of management factor [34]: where a value of 0.176 was used for CA. This value was derived from the multiplication of the three factors adopted for CA (C tillage = 0.25; C residues = 0.88; C cover = 0.80) [34], while a value of 1 was used for CM, in accordance with the multiplication of three factors adopted for CM (C tillage = 1; C residues = 1; C cover = 1) [34]. Instead, C crop was estimated by taking into account the Normalized Difference Vegetation Index (NDVI), as proposed by van der Knijff et al. [12], for regional-scale applications: where α and β are parameters of the NDVI-C correlation. An α-value of 2 and a β-value of 1 seem to give reasonable results, because these values permit achieving a linear relationship, according to van der Knijff et al. [66]. This method has been employed by several studies worldwide [67][68][69][70]. In the present work, the C-factor was calculated for Apulia ACL as follows: For the calculation of C-factor's NDVI, GEE provides Sentinel-2 images that have a resolution of 10 m and are available for the two layers: bottom of atmosphere (BOA), level 2A) and top of atmosphere (TOA), level 1C). The TOA level is not provided with atmospheric correction, while the BOA level has atmospheric correction due to Sen2Cor [71]. The noncorrect atmospheric images have been continuously available since the launch of the satellite (23 May 2015), while the corrected ones have been available since 28 March 2017. As shown in the literature [72,73], there is a correlation between vegetational indices calculated at the two different layers (NDVI BOA > NDVI TOA). Subsequently, the correction factor between NDVI is calculated at the two levels, without which the C-factor would have been underestimated. We proceeded by calculating the NDVI Sentinel-2 at level 2a over the whole region for the interval from 1 April 2017 to 28 January 2021, by taking into consideration 2033 images in which NDVI was calculated both at the BOA and TOA level. The correction factor shows an average increase between TOA and BOA of 27% in the whole Apulia region and an increase of 29% in ACL. For the period 2016-2017, NDVI is calculated by using NDVI level 1c added to correction factor. The annual interval for each NDVI calculation is considered for the period between 15 October to 1 August of the following year, as it corresponds approximately to the period from sowing to harvesting in Southern Italy. Table 1 reports the number of images that have a cloudy pixel percentage lower than 20% and are processed for the calculation of each NDVI. Finally, the calculation of the two factors, C crop [38,74] and C management, and their multiplication, is possible in GEE environment.

Identification of ACL
ACL areas are calculated using LUCAS points, which is a network that gives accurate and detailed information, it aims at computing statistical estimates at EU level with fine scale [75], and it collects information on European land. This examination provides spatial information that can be used for agricultural goals and to define the impact on the environment and natural resources. Each LUCAS point collects information including land cover, land use, and environmental parameters, identified as microdata [76]. For the calculation of crop land, LUCAS microdata was used (from class B10 to class B50). These data are available online and both are downloadable; for the year 2015 (https://ec. europa.eu/eurostat/web/lucas/data/primary-data/2015, accessed on 1 November 2021), for the year 2018 (https://ec.europa.eu/eurostat/web/lucas/data/primary-data/2018, accessed on 1 November 2021). The areas were calculated in QGIS (QGIS.org, 2021, QGIS Geographic Information System, QGIS Association. http://www.qgis.org, accessed on 1 November 2021) following Gallego and Bamps [75]. ACL, mapped by LUCAS in the Apulia region, resulted in 733,801 ha and 773,828 ha for 2015 and 2018, respectively. Subsequently, areas detected by LUCAS were mapped by using Global Land Cover [77,78] in Google Earth Engine [79]. In Global Land Cover, cultivated areas are aggregated into crop land. The separation between ACL and crop land is performed through a multitemporal analysis, by using the LANDSAT8 Collection 1 Tier 1 composite dataset (https://developers. google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C01_T1_8DAY_NDVI, accessed on 1 November 2021) which allows the calculation of the averaged NDVI-using six images for 2015 and six images for 2018. This allows the creation of a mask to distinguish ACL from croplands, setting the threshold between 0.08 and 0.30 during summer periods (1 July-15 August), in which differences between these two categories are more evident in Mediterranean areas. This threshold guarantees significant difference to areas detected by LUCAS. This analysis makes it possible to map ACL areas-732,052 ha for 2015 and 772,654 ha for 2018, respectively (

RUSLE Factors Multiplication
The multiplication of all the RUSLE factors was carried out in GEE, carrying all the factors to the same resolution scale, 10 m resampling, and reducing resolution for each factor by using bilinear interpolation, and with the same EPSG:3035 reference system ( Figure 4). RUSLE is calculated for four consecutive annual production cycles, as shown in Table S2.
For an optimal understanding of RUSLE distribution in the Apulia region, altitudes and slopes were calculated by using DTM in GEE. Altimetry was classified in three categories (according to the Italian Institute of Statistics-ISTAT, https://www.istat. it (accessed on 1 November 2021), based on their altitude above sea level: plain (0 ≤ altitude ≤ 300 m a.s.l.), hilly (300 < altitude < 800 m a.s.l.), and mountain (altitude ≥ 800 m a.s.l.). Each of these three categories are designated to a numerical class: 100 for the plain, 200 for the hilly, and 300 for the mountain. The slope's division in categories is carried out by dividing the slope range into three quantiles: low slope (under 1.8%), medium slope (1.8-3.7%), and high slope (over 3.3%). Each class was renominated by using a number: 1 for low slope, 2 for medium slope, and 3 for high slope. Moreover, the two categories were combined to obtain nine classes for both the scenarios (Table S3). As a final step, the statistics for each category were calculated.

RUSLE Factors Multiplication
The multiplication of all the RUSLE factors was carried out in GEE, carrying all the factors to the same resolution scale, 10 m resampling, and reducing resolution for each factor by using bilinear interpolation, and with the same EPSG:3035 reference system (Figure 4). RUSLE is calculated for four consecutive annual production cycles, as shown in Table S2. For an optimal understanding of RUSLE distribution in the Apulia region, altitudes and slopes were calculated by using DTM in GEE. Altimetry was classified in three categories (according to the Italian Institute of Statistics-ISTAT, https://www.istat.it (accessed on 1 November 2021), based on their altitude above sea level: plain (0 ≤ altitude ≤ 300 m a.s.l.), hilly (300 < altitude < 800 m a.s.l.), and mountain (altitude ≥ 800 m a.s.l.). Each of these three categories are designated to a numerical class: 100 for the plain, 200 for the hilly, and 300 for the mountain. The slope's division in categories is carried out by dividing the slope range into three quantiles: low slope (under 1.8%), medium slope (1.8-3.7%), and high slope (over 3.3%). Each class was renominated by using a number: 1 for low slope, 2 for medium slope, and 3 for high slope. Moreover, the two categories were combined to obtain nine classes for both the scenarios (Table S3). As a final step, the statistics for each category were calculated.

Rainfall Erosivity (R-Factor)
The average annual R-factor for the Apulia region in the 1981-2019 year-period totals 318.3 MJ mm h −1 ha −1 y −1 , with a standard deviation of 57.4 MJ mm h −1 ha −1 y −1 (Table S1) (Table S1; Figure 5). The high R-factor corresponds to mountain areas ( Figure 5) proximity to the Daunian Sub-Apennine, in the province of Foggia in the northwest, and to the Murgia plateau in the provinces of Barletta-Andria-Trani and Bari.
Erosivity Database on the European Scale (REDES) (Borrelli et al. [80]; Panagos et al. [81]). Considering the ACL of the Apulia region, R-factor ranged from 280.4 to 284.4 MJ mm h −1 ha −1 y −1 (1981-2016 and 1981-2019, respectively) (Table S1; Figure 5). The high R-factor corresponds to mountain areas ( Figure 5) proximity to the Daunian Sub-Apennine, in the province of Foggia in the northwest, and to the Murgia plateau in the provinces of Barletta-Andria-Trani and Bari.

Soil Erodibility (K) and Support Practice Factor (P)
The spatial variation of K-factor values from the ESDAC dataset in the Apulia region depend on the variation of ACL area for both years. The results showed a mean value of 0.0331 t ha h ha −1 MJ −1 mm −1 for 2015 ( Figure 6a

Soil Erodibility (K) and Support Practice Factor (P)
The spatial variation of K-factor values from the ESDAC dataset in the Apulia region depend on the variation of ACL area for both years. The results showed a mean value of 0.0331 t ha h ha −1 MJ −1 mm −1 for 2015 ( Figure 6a In both years, the highest value of K-factor is distributed in the areas of Tavoliere and in the Daunian Sub-Apennine, in the province of Foggia in the northwest, and in the Murgia plateau in the provinces of Barletta-Andria-Trani and Bari. From the ESDAC dataset for the P-factor, it was possible to extrapolate factors for the ACL area for both years. The P-mean for 2015 is 0.8736, and for 2018, 0.8682. In both maps (Figure 7), the higher values are found in the mountainous regions where the CA system is mainly applied. In both years, the highest value of K-factor is distributed in the areas of Tavoliere and in the Daunian Sub-Apennine, in the province of Foggia in the northwest, and in the Murgia plateau in the provinces of Barletta-Andria-Trani and Bari. From the ESDAC dataset for the P-factor, it was possible to extrapolate factors for the ACL area for both years. The P-mean for 2015 is 0.8736, and for 2018, 0.8682. In both maps (Figure 7), the higher values are found in the mountainous regions where the CA system is mainly applied. In both years, the highest value of K-factor is distributed in the areas of Tavoliere and in the Daunian Sub-Apennine, in the province of Foggia in the northwest, and in the Murgia plateau in the provinces of Barletta-Andria-Trani and Bari. From the ESDAC dataset for the P-factor, it was possible to extrapolate factors for the ACL area for both years. The P-mean for 2015 is 0.8736, and for 2018, 0.8682. In both maps (Figure 7), the higher values are found in the mountainous regions where the CA system is mainly applied.

Topographic Factor (LS)
LS-factor is calculated for the ACL of the Apulia region. The averaged values ( Figure  8  The high LS values are found in the mountainous area with steep topography, especially in the west where mountains are prevalent (Daunian Sub-Apennine and Murgia), while the lowest values are distributed along the Adriatic and Ionic coasts. To determine the spatial resolution (cell size) of the soil erosion model results, and, therefore, to incorporate the soil erosion potential due to surface runoff, we used a high-resolution (8 m) DTM of the Apulia region to calculate LS-factor. This resolution of Apulia region DTM is the best available at regional scale. Generally, the DTM resolutions used on a European scale are in the range of 25 to 100 m [56], while at regional scale, the resolutions range from 5 m to 40 m [82,83]. As reported by Bircher et al. [84], several authors have mentioned the importance of assessing the risk of soil erosion based on the size of the cells and the accuracy of the DTM. The risk that can occur is that a low DTM resolution (large cells) is not able to map relevant topographic details [85][86][87]. The high LS values are found in the mountainous area with steep topography, especially in the west where mountains are prevalent (Daunian Sub-Apennine and Murgia), while the lowest values are distributed along the Adriatic and Ionic coasts. To determine the spatial resolution (cell size) of the soil erosion model results, and, therefore, to incorporate the soil erosion potential due to surface runoff, we used a high-resolution (8 m) DTM of the Apulia region to calculate LS-factor. This resolution of Apulia region DTM is the best available at regional scale. Generally, the DTM resolutions used on a European scale are in the range of 25 to 100 m [56], while at regional scale, the resolutions range from 5 m to 40 m [82,83]. As reported by Bircher et al. [84], several authors have mentioned the importance of assessing the risk of soil erosion based on the size of the cells and the accuracy of the DTM. The risk that can occur is that a low DTM resolution (large cells) is not able to map relevant topographic details [85][86][87].

Cover-Management (C-Factor)
The four-year mean for the C-factor values of the two management systems is shown separately for each year and crop stage in Table 1. The four-year mean C-factor for CM in ACL of the Apulia region is 0.2867, while for the CA it is 0.2745, with an average percentage reduction of the C-factor of 4.2% according to the cropping system adopted. The lowest C-factor values for CM and CA were registered both in the year 2019/2020 (0.2656 and 0.2552, respectively), while the highest were in the year 2016/2017 (0.3268 and 0.3122, respectively; Supplementary Figures S1 and S2). The four-cropping cycle average of NDVI for the study site varied from 0.3653 to 0.4078 (in 2016/2017 and 2019/2020, respectively) with an average over 0.3928. The ground-cover percentage was directly measured in all ACL with 10 m resolution. Based on the two management systems in the study area, higher C-values were observed in the CM management system, thus indicating higher potentiality for soil erosion risk in these areas ( Figure 9). Spatiotemporal dynamics are required to understand impacts and risks at regional scale. For this reason, it is necessary to identify the agricultural area that can be most affected by soil erosion by water and identify, where possible, which CA system will reduce these processes. Our approach was focused on the use of multitemporal satellite images to capture the temporal variability of the determination of factor C. This approach is used for ecosystem modeling studies on a large scale [34,88] and places the focus on the importance of explicit consideration of temporal variability on soil management systems to protect agricultural land from the impact of soil erosion by water [89]. Generally, in large-scale modeling applications, the estimation of multiple sub-factor C parameters is derived from preexisting literature [34,89,90] and does not use spatially explicit land-cover data for the determination of the C-factor; it is, rather, based on statistical data. Current approaches on a European scale adopt the CORINE Land Cover Database for the calculation of the C-factor [12]. This data is inadequate in its spatial and thematic resolution compared to RS sensors with up-to-date information on proximal measures of vegetation, which is a key aspect for a correct and accurate calculation of C-factor [22]. This aspect becomes even more important for areas where the spatial and temporal dynamics of the vegetation cover are provided by the cultivation of crops [83,91]. On the other hand, data obtained from field experiments carried out to measure the values of the C-factor take a long time and are rarely available [22].
Within this work, by using a spatial resolution of 10 m (obtained through Sentinel a b Our approach was focused on the use of multitemporal satellite images to capture the temporal variability of the determination of factor C. This approach is used for ecosystem modeling studies on a large scale [34,88] and places the focus on the importance of explicit consideration of temporal variability on soil management systems to protect agricultural land from the impact of soil erosion by water [89]. Generally, in large-scale modeling applications, the estimation of multiple sub-factor C parameters is derived from preexisting literature [34,89,90] and does not use spatially explicit land-cover data for the determination of the C-factor; it is, rather, based on statistical data. Current approaches on a European scale adopt the CORINE Land Cover Database for the calculation of the C-factor [12]. This data is inadequate in its spatial and thematic resolution compared to RS sensors with up-to-date information on proximal measures of vegetation, which is a key aspect for a correct and accurate calculation of C-factor [22]. This aspect becomes even more important for areas where the spatial and temporal dynamics of the vegetation cover are provided by the cultivation of crops [83,91]. On the other hand, data obtained from field experiments carried out to measure the values of the C-factor take a long time and are rarely available [22]. Within this work, by using a spatial resolution of 10 m (obtained through Sentinel 2A) and applying the formula proposed by Van der Knijff et al. [12] integrated with that of Panagos et al. [34], we obtained a more detailed resolution on the ACL (previously mapped) to multitemporal NDVI calculation. This NDVI-derived method can variably capture the actual soil cover status [82] rather than using aggregate and fixed data over time. This approach can be useful in increasing the accuracy of the calculation of the C-factor [67,92,93]. Our results show how, for the estimation of the C-factor in the Apulia region, the values are comparable to those present in the literature for ACL [67,74,94,95], which fluctuate between 0.01 and 0.44 for ACL, but with a much more accurate regional scale of detail. Concerning the difference between the CM and CA systems, our results demonstrate how, on average, the adoption of CA reduces the C-factor by 4.2% in line with other analyzed scenarios [34]. CA has been reported in many studies as an effective strategy to control erosion processes, maintain soil fertility, increase soil carbon sequestration, and improve cropping system sustainability [96][97][98].

Loss Rates
The loss rates (A), in t ha −1 y −1 , for the Apulia region in this study are generated by using the RUSLE model to calculate the mean for the 2016-2020 period. The results are shown in Table 2. The four-year agricultural annual crop cycle average for CM is 2.28 t ha −1 y −1 , while for CA it is 2.11 t ha −1 y −1 (Figure 10).

Loss Rates
The loss rates (A), in t ha −1 y −1 , for the Apulia region in this study are generated by using the RUSLE model to calculate the mean for the 2016-2020 period. The results are shown in Table 2. The four-year agricultural annual crop cycle average for CM is 2.28 t ha −1 y −1 , while for CA it is 2.11 t ha −1 y −1 (Figure 10). Table 2. RUSLE values for the 2016-2020 period in the Apulia region for each management system.   Figures S3 and S4 and Table S4.

RUSLE
Interestingly, the areas featuring high risk of soil erosion are in the Tavoliere and in the Daunian Sub-Apennine, in the northwest, as well as in the Murgia plateau in the center of the Apulia region. As for the soils under CA, there is instead a reduction in soil erosion risk in the same areas as well as in the south of the Apulia region (Salento, province of Lecce).
Current results from the RUSLE model reveal a serious soil erosion risk where the CM system is adopted, while the CA system showed a trend over the years to contain the rate of soil loss for ACL in the Apulia region, representative of pedoclimatic conditions widely present in large agricultural areas of the Mediterranean basin. As suggested by Wischmeier and Smith [52], practices of improved tillage, such as no-till and cover crops, were considered as conservation cropping and management practices and implemented in the C-factor. A recent scenario analysis carried out on 54 countries, that reported information on the implementation of CA to the Food and Agriculture Organization of United Nations (FAO), assumes a 45% reduction of soil erosion risk in CA compared to CM [99]. Furthermore, a previous study by Borrelli et al. [100] shows how in Italy, when good agricultural and environmental conditions (GAEC) including CA are adopted, there is a potential to reduce erosion by 8.5%, while for the Apulia region, the potential rate of erosion is 5%. The results obtained ( Table 2) at this territorial scale are in line with this scenario and, moreover, by improving the spatial resolution of DTM (8 m), the erosion rate calculated over the four years (2016-2020) shows an average of 7%, in line with the scenario assumed for the Apulia region.

Mean Loss Rates for Altitude and Slope Classes
RUSLE values were calculated for altimetric and slope classes for each agricultural annual crop cycle ( Table 3).
The erosion rate is variable across the regional territory because the pedological and climatic variability is high, but even more so is the morphology of the Apulia region, which is also found in the discontinuous adoption of CA. In the Apulia region, the ACL lying on plain (0 ≤ 300 m a.s.l.) and hilly (300 m < 800 m a.s.l.) areas covers 99% of the surfaces, and about 36% of the study areas is located on slopes greater than 3.7%. In our study, with the introduction of the CA management system, the soil loss rate in these areas ranged from −1.7% for the plain + medium slope terrain to −18.8% for the mountain + low slope terrain, with an overall average of all terrain classes considered, during the four years, of −8.5%.
In the plains, 68% of the total ACL is present, and in four years, the contribution of CA can reduce the soil loss rate by −3.7%, compared to the CM system. For the hilly areas (32% of the total ACL), the erosion rate decreases by −10.1% when CA is adopted. ACL distributed in the slope classes are very similar, and the contribution of CA is higher on the high slopes, decreasing the erosion rate by −7.6% compared with the CM system, according to the scenario analysis for Italy by Borrelli et al. [100].

Combination of Altitude and Slope Classes
In the combination of altitude and slope classes, the highest effect is on the hilly + low slope and hilly + medium slope. The CA management system contributes, respectively, −14.1% and −12.2% to the erosion rate over the four years, compared with CM. These values are in line with the scenario proposed by Panagos et al. [101], in which a projection of soil loss by water erosion in Europe by 2050 is calculated, and it is reported that there will be a potential of reduction of soil losses of 17-22% with the contribution of reduced tillage combined with cover crops. This scenario is confirmed by the results of Table 3, in which the highest absolute values of soil loss rate containment are obtained for the altitude and slope classes in which CA is currently marginally adopted.

Agricultural Seasons
lated over the four years (2016-2020) shows an average of 7%, in line with the scenario assumed for the Apulia region.

Mean Loss Rates for Altitude and Slope Classes
RUSLE values were calculated for altimetric and slope classes for each agricultural annual crop cycle (Table 3). The erosion rate is variable across the regional territory because the pedological and climatic variability is high, but even more so is the morphology of the Apulia region, which is also found in the discontinuous adoption of CA. In the Apulia region, the ACL lying on plain (0 ≤ 300 m a.s.l.) and hilly (300 m < 800 m a.s.l.) areas covers 99% of the surfaces, and about 36% of the study areas is located on slopes greater than 3.7%. In our study, with the introduction of the CA management system, the soil loss rate in these areas ranged from −1.7% for the plain + medium slope terrain to −18.8% for the mountain + low slope terrain, with an overall average of all terrain classes considered, during the four years, of −8.5%.

Numerical Class
Altitude

Conclusions
Soil erosion is among the most critical environmental hazards of modern times. Vast areas of Mediterranean land now being cultivated with CM may be rendered unproductive, or at least economically unproductive, if erosion continues unabated. Soils cultivated with annual crops in Mediterranean climatic conditions under conservative agriculture can benefit from a permanent cover for the direct increase of surface water infiltration, significantly reducing surface runoff and therefore soil erosion risk.
In order to estimate soil loss rate at regional scale, empirical models are accurate and easy to interpret and require modest resources. They can be processed with readily available inputs to identify areas exposed to high risk of erosion. In this study, RUSLE models integrated with GEE and QGIS were used to estimate soil loss rate on the ACL of the Apulia region for a period of four annual crop cycles-from 2016 to 2020-for both the scenarios.
Results show that where the CA system is applied continuously in ACL, there is an annual average reduction of soil loss rate over 7% compared with CM; furthermore, it is significatively different for altimetric, slope, and combined classes, showing that the important contribution of the CA system can reduce soil loss rate in hilly areas by 10.1% and in hilly + low slope terrain by 14.1%. These results represent a baseline to estimate the effects of the adoption of the specific agro-environment-climate measure on soil erosion risk during the first phase of transition from conventional to conservation management systems. Consequently, the results of this study can represent an objective target baseline for the planning of the new CAP 2023-2027, which provides for the selection of reliable concrete and achievable result indicators, including erosion by water, whose values can be monitored and verified periodically. The goal is to increase the agricultural area under CA in the Apulia region and in those areas with semiarid Mediterranean climate where there is greater loss of soil due to water erosion. Such data, which should be monitored periodically, could be used to evaluate soil conservation management planning processes, and help determine the dimension and duration of a transition phase to support farmers in providing ecosystem restoration services, reduction of erosion by water, and improvement of soil healthy and agricultural productivity.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/agronomy12020281/s1, Figure S1: Annual average of cover and management (C) factor of Apulia region under conventional management for annual cropland; Figure S2: Annual average of cover and management (C) factor of Apulia region under conservation agriculture for annual cropland; Figure S3: Annual soil loss rate using RUSLE model for Apulia region under conventional management for annual cropland; Figure S4: Soil loss rate using RUSLE model for Apulia region under conservation agriculture for annual cropland; Table S1: Descriptive statistical values of the R factor (MJ mm h −1 ha −1 y −1 ) for each year-periods of RUSLE calculation; Table S2: Calculation of annual cropland areas in Apulia region in four different years; Table S3: Nine new classes obtained by combined altitude class with three quantiles; Table S4: Intermediate RUSLE values in the Apulia region for each management system.