Application of the Adapted Approach for Crop Management Factor to Assess Soil Erosion Risk in an Agricultural Area of Rwanda

: Land use and land cover (LULC) management inﬂuences the severity of soil erosion risk. However, crop management (C) is one factor of the Revised Universal Soil Loss Equation (RUSLE) model that should be taken into account in its determination, as it inﬂuences soil loss rate estimations. Thus, the present study applied an adapted C-factor estimation approach (CvkA) modiﬁed from the former approach (Cvk) to assess the impact of LULC dynamics on soil erosion risk in an agricultural area of Rwanda taking the western province as a case study. The results disclosed that the formerly used Cvk was not suitable, as it tended to overestimate C-factor values compared with the values obtained from t CvkA. An approximated mean soil loss of 15.1 t ha − 1 yr − 1 , 47.4 t ha − 1 yr − 1 , 16.3 t ha − 1 yr − 1 , 66.8 t ha − 1 yr − 1 and 15.3 t ha − 1 yr − 1 in 2000, 2005, 2010, 2015 and 2018, respectively, was found. The results also indicated that there was a small increase in mean annual soil loss from 15.1 t ha − 1 yr − 1 in 2000 to 15.3 t ha − 1 yr − 1 in 2018 (1.3%). Moreover, the soil erosion risk categories indicated that about 57.5%, 21.8%, 64.9%, 15.5% and 73.8% had a sustainable soil erosion rate tolerance ( ≤ 10 t ha − 1 yr − 1 ), while about 42.5%, 78.2%, 35.1%, 84.5% and 16.8% had an unsustainable mean soil erosion rate (>10 t ha − 1 yr − 1 ) in 2000, 2005, 2010, 2015 and 2018, respectively. A major portion of the area fell under the high and very high probability zones, whereas only a small portion fell under the very low, low, moderate and extremely high probability zones. Therefore, the CvkA approach presents the most suitable alternative to estimate soil loss in the western province of Rwanda with reasonable soil loss prediction results. The study area needs urgent intervention for soil conservation planning, taking into account the implementation of effective conservation practices such as terracing for soil erosion control.


Introduction
Globally, soil erosion incidences enhanced by human activities have been a serious environmental problem since the last century [1]. They are a known risk factor to an ecosystem's integrity and cause major environmental impacts, including a reduction in Land 2021, 10, 1056 2 of 24 agricultural productivity; they can even affect the freshwater ecology of different regions. Previous studies have reported a soil loss of approximately 73.5 Pg yr −1 in agricultural land [2,3]. Topography, soil characteristics, local climatic conditions and land use are considered as the main influencing factors of increasing soil erosion risks [3][4][5].
However, it is not easy to estimate soil loss without relating it to various human, natural or physical processes. Quantitative and qualitative descriptions of soil erosion have been executed by field observations and measurements and via meta-analyses of soil erosion levels [6]. Owing to the scarcity of in situ or field data, as well as their expensiveness and time-consuming aspects, researchers have established quantitative soilerosion modeling approaches that apply physical factors including topography, climate, soil features and vegetation type, with the aim of mapping spatial distribution rates to better understand the mechanisms of soil erosion. Moreover, these erosions provide a clear understanding of natural phenomena such as the transportation and deposit of sediment, including soil erosion prediction [7,8].
It is difficult to find a model that matches all types of erosion. This situation contributes to the intricacy in identifying areas prone to practicing soil conservation measures in terms of a reduction in soil loss. Therefore, the development and planification of numerous alternative land-use scenarios are of great significance as far as the management of soil erosion risk is concerned [9].
The primary challenge in modeling soil erosion risk is mostly associated with accuracy assessments, especially in data-scarce regions such as Rwanda [4]. Therefore, it is challenging to determine the degree of risk without the existence of numerous soil erosion models. For the prediction of soil erosion, different approaches have been proposed by different authors depending on study objectives, time, data availability and spatial scales [10][11][12]. The approaches for modeling soil erosion are easily understandable because of their powerful reliance on geospatial techniques that include remote sensing (RS) and geographic information system (GIS) platforms [11]. To this end, the USLE (Universal Soil Loss Equation) and its modified form, RUSLE (the Revised Universal Soil Loss Equation), proposed by [13], are largely applied due to their ability to integrate the aforementioned geospatial techniques. The integrated application of RUSLE and GIS combined with RS technology has been verified as an effective approach for assessing the spatial variability of soil erosion occurrence and its drivers [11,[14][15][16][17].
In Rwanda, studies have applied the RUSLE model to simulate soil erosion at both national and local scales. However, those considering agricultural areas are still understudied. Among the few existing studies, Nyesheja and Chen [18] modeled soil erosion in the Congo Nile Ridge region of Rwanda by estimating serious soil loss from different land use/land cover (LULC) types such as grasslands, forestlands and farmlands that expressed impaired ecosystems. Additionally, Byizigiro et al. [19] estimated soil erosion in the Ngororero district taking the Satinskyi Catchment as a case and disclosed the geomorphological features and processes of the area to be the major natural factors that aggravated increased soil erosion rates pronounced in the area. In conjunction with its seven-year development strategies on the sensitization of the population regarding aspects related to soil conservation, Niyonsenga et al. [20] considered the Nyamasheke district to spatially analyze the unknown sensitivity of soil erosion. They reported that soil erosion and the associated high level of land degradation caused by vegetation cover change are leading to a loss in fertile soil in the area.
Nevertheless, keeping the basic form of the RUSLE equation, these studies presented the overestimation of the crop management (C) factor, which seemed to provide unsatisfactory results in the considered areas [21,22]. In addition, all the above-mentioned studies were limited to only the district and catchment scales of the province opposite this study, which considered the agricultural areas of the entire western province. Useful information related to LULC change, erosional features such as rainfall interception and gullies has been delivered by multi-temporal satellite images to analyze the LULC change in relation to soil erosion hazards at a regional scale [23]. LULC change is central to the study of global environmental change, including land degradation [24]. It is the most significant indicator between a natural environment and human activity. Rwanda is considered one of the regions highly vulnerable to water-induced soil erosion [25]. This is because of its rapid population growth with a low economic base, poor agricultural practices, fragile soil [26] and steep topography, which contributes to heavy rainfall [27]. Similarly, the effects of LULC change studies have been conducted at different small watersheds at district, national, regional and global levels [11,[28][29][30]. All these studies agree on the strong influence of land-use change on soil erosion.
Rwanda has experienced a significant regional change in its land-use patterns, especially in the western province, with limited scientific studies tackling this issue. In line with this, Nambajimana and He [21] used the RUSLE model for Land-Use Change Impacts on water erosion in Rwanda and identified the driving forces of soil erosion; however, this study was limited by its inability to capture degraded areas in its impact assessment. This can have some unforeseen significant impacts on land area and the wellbeing of local residents. Based on this background and the context of evolving land-use changes in Rwanda, this present study intends to address all of the above-mentioned shortcomings in the previous studies, and establish an adapted approach that is suitable for the estimation of C-factor in the RUSLE model that exhibits overestimated values. Thus, the objectives of this study are to: (a) examine the adapted approach for the estimation of the crop management (C) factor in the western province of Rwanda; (b) analyze the LULC changes and their impacts on soil erosion risk; (c) estimate the potential annual soil loss rate for the study area; and (d) delineate the spatial distribution of soil erosion probability zones. Addressing these objectives will fill the knowledge gap in identifying the most degraded parts of agricultural lands in the study area, which is crucial for effective soil conservation planning and environmental management.

Study Area
The study covers the entire western part of Rwanda with −2 • 22'" and 17.99" south latitude and 29 • 12" and 21 • 00" east longitude spatial coordinates. It is made of seven districts, namely Ngororero, Karongi, Rubavu, Nyabihu, Rutsiro, Rusizi and Nyamasheke ( Figure 1). The area is mountainous and covers 5883 km 2 , with the highest altitude in the country. It has an elevation ranging between 921 m and 4492 m, with an annual average temperature of 18 • C. The country's climate is made up of four seasons, namely the long dry (June-September), short dry (mid-December-mid-February), long rainy (late February-late May) and short rainy (late September-early December) seasons [31]. The area receives up to 1800 mm annual rainfall, which usually contributes to high average runoff [18]. High annual precipitation and steep topography (> 26 • ) constitute major natural factors that influence morphological developments in the area. The latter is aggravated by agriculture and mining operations, which contribute to increased soil erosion rates observed in the area.

Data and Preprocessing
For the successful application of the RUSLE model in this study, GIS and RS techniques were used for satellite data processing and analysis. The Landsat 7 ETM+ and Landsat 8 OLI-TIRS images (30m ×30 m, WRS path 172 and rows 61) products were acquired from the United States Geological Survey (USGS) EROS data center. These images were pre-processed for radiometric and atmospheric rectification [32] and then used for calculating the vegetation index, as well as LULC classification using a supervised maxi-

Data and Preprocessing
For the successful application of the RUSLE model in this study, GIS and RS techniques were used for satellite data processing and analysis. The Landsat 7 ETM+ and Landsat 8 OLI-TIRS images (30 m × 30 m, WRS path 172 and rows 61) products were acquired from the United States Geological Survey (USGS) EROS data center. These images were pre-processed for radiometric and atmospheric rectification [32] and then used for calculating the vegetation index, as well as LULC classification using a supervised maximum classification algorithm Land 2021, 10, 1056 5 of 24 supported by expert judgment and ground truth observation. The main land-cover classes were water bodies, forest, grassland, cropland, wetland and settlements. In addition, the Digital Elevation Model (SRTM-DEM) at 30 m resolution (http://www.earthexplorer.usgs.gov, accessed on 4 October 2021) was acquired to calculate the flow accumulation and the slope gradient (degrees), which were both used to generate the factors representing the length of the slope and its steepness (LS). Moreover, a soil-type map from the Food and Agriculture Organization of the United Nations (FAO) was used in addition to meteorological data (rainfall) acquired from the Rwanda Meteorological Agency. Further details of the datasets are illustrated in Table 1.

Estimations Based of RUSLE Model
RUSLE is an erosion empirical model aiming to predict the long-term averages of soil loss per year caused by runoff from specific slopes area in specified cropping and management systems. The RUSLE model collects many factors on the process of soil erosion into different groups (relief, climate, vegetation, land use, soil profile and the practice of land management). These five groups are subdivided as soil erosion factors, including the rainfall erosivity (R) factor, the soil erodibility (K) factor, the slope length and slope steepness (LS) factor, the crop management (C) factor and the erosion management practice (P) factor. The expected soil loss in t/ha/yr was estimated by multiplying these factors' values based on the dimensions used in the climate and soil factor.
The RUSLE equation is given below: where A is the estimated spatial average soil loss and temporal average soil loss per unit area, expressed in t/ha/yr.

Rainfall Erosivity (R) Factor
Rainfall is a key factor triggering soil erosion in the tropic region. The R factor defines the erosive power of precise rainfall or the intensity of rainfall as the factor causing soil erosion [33]. The understanding of this factor relies on the kinetic force of raindrops on the soil's surface [34]. In the RUSLE model, R factor estimation depends on the product given by the storm power based on a rainfall intensity of 30 min, expressed as R = EI30 [34]. The R factor (Figure 2a) was estimated using Equation (2) proposed by [35]: where R represents the rainfall erosivity factor while the mean annual rainfall (mm) is represented by P.

Crop Management (C) Factor
C is a factor representing crop management as a fraction of soil loss from land collected under specified conditions to the soil loss from the same area in tilled unceasing fallow land. It represents the effects of vegetation and other ground covers [38]. The traditional approach for the calculation of the C-factor in space is allocating values to LULC classes. The effects of vegetation are generally expressed using different vegetation indices, whereby each index has its specific expression of green vegetation and its particular suitability for specific applications. Therefore, the choice of a specific vegetation index or soil biophysical indices needs to be made cautiously by widely considering and analyzing their prevailing advantages and shortcomings. For this specific study, the Normalized Difference Vegetation Index (NDVI) was deemed worthy of being used for the estimation of the RUSLE-based C-factor due to its important advantages such as its ability to indicate green vegetation potency [39], show fluctuations in LULC [40] and estimate values that represent the utmost heterogeneity in the stage of crop growth [41]. Different researchers developed many techniques to calculate the C-factor using the NDVI for assessing soil loss with the USLE/RUSLE model. This index can be derived from the spectral reflectance difference between near infrared (NIR) and red. It ranges between −1 and +1, with the highest values accredited to areas with higher vegetation cover [42,43], and it is estimated using the equation below: We used NDVI and two approaches for determining the soil cover management factor. The first (Equation (6)) was proposed by Van der Knijff and Jones [22], and the second

Soil Erodibility (K) Factor
The K-factor expresses the susceptibility of soil-to-soil erodibility. K-factor ( Figure 2b) was computed based on the soil properties, including a fraction of sand, clay, silt and organic carbon using the Equation (3): where A represents the soils with high-and rough-sand contents and those with less sand (Equation (3a)), B specifies the soils with high clay-to-silt fractions (Equation (3b)), C is soils rich in organic carbon (Equation (3c)) and D defines soils very rich in sand (Equation (3d)). Practically, A, B, C and D were multiplied with 0.1317 values to convert the K-factor from the American system to the metric system unity/International System of Units (SI). where SAN is the sand content percentage (0.05-2.00 mm diameter); SIL specifies the silt content percentage (0.002-0.05 mm diameter); CLA specifies the clay content percentage (<0.002 mm diameter); C is the organic carbon content percentage; and SN1 = 1 − (SAN/100).

Slope Length and Its Steepness (LS) Factor
LS is a factor derived from DEM obtained by multiplying the length (L) and steepness (S) of a slope. L is the fraction of soil loss from the ground to soil loss from 22.1 meters under similar conditions, while S is the fraction of soil loss from the ground slope gradient to the loss of soil from 9% under similar conditions. L was derived from an algorithm (Equation (4)) developed by Desmet and Govers [36], while the S factor was computed using the McCool and Brown [37] method (Equation (5)). Finally, the LS (Figure 2c) was calculated using the equations below: where L and S indicate the slope length and steepness factor (dimensionless); λ is the slope length (in m); m denotes the variable power based on E; E denotes the rill to interrill erosion fraction; and θ is the slope angle (in degrees). S1 and S2 represent the slope percentage <9% and slope percentage ≥ 9%, respectively.

Crop Management (C) Factor
C is a factor representing crop management as a fraction of soil loss from land collected under specified conditions to the soil loss from the same area in tilled unceasing fallow land. It represents the effects of vegetation and other ground covers [38]. The traditional approach for the calculation of the C-factor in space is allocating values to LULC classes. The effects of vegetation are generally expressed using different vegetation indices, whereby each index has its specific expression of green vegetation and its particular suitability for specific applications. Therefore, the choice of a specific vegetation index or soil biophysical indices needs to be made cautiously by widely considering and analyzing their prevailing advantages and shortcomings. For this specific study, the Normalized Difference Vegetation Index (NDVI) was deemed worthy of being used for the estimation of the RUSLE-based C-factor due to its important advantages such as its ability to indicate green vegetation potency [39], show fluctuations in LULC [40] and estimate values that represent the utmost heterogeneity in the stage of crop growth [41]. Different researchers developed many techniques to calculate the C-factor using the NDVI for assessing soil loss with the USLE/RUSLE model. This index can be derived from the spectral reflectance difference between near infrared (NIR) and red. It ranges between −1 and +1, with the highest values accredited to areas with higher vegetation cover [42,43], and it is estimated using the equation below: We used NDVI and two approaches for determining the soil cover management factor. The first (Equation (6)) was proposed by Van der Knijff and Jones [22], and the second was suggested as an adapted version of the first approach (Equation (7)). In the approach proposed by Van der Knijff and Jones [44] (Cvk), the C-factor is calculated as: where α and β are the parameters that determine the shape of the curve that links NDVI with the C-factor. Van der Knijff and Jones [22] used values 2 and 1 for European climate conditions as the best representatives of the equation parameters α and β, respectively. Nevertheless, for tropical climate conditions with high intensities of rainfall, the C-factor tends to be overestimated compared to that calculated using the approach by Van der Knijff and Jones [22] for the same vegetation cover. Hence, a newly adapted approach for calculating the RUSLE C-factor was presented in this study in comparison to the former one.
Due to the bias identified in the previous studies, this study recommends an adjustment factor of 2, which has been divided by the aforementioned equation (Equation (6)) to minimize the bias for the tropical climate conditions. Therefore, the adapted C-factor (CvkA) was proposed as follows: To establish the adequacy of the used index (NDVI) in the estimation of the RUSLEbased C-factor, a linear regression analysis was then executed using the extracted values against the related NDVI-based C-factor values [42]. Additionally, the shape of the mean cover management factor values calculated for all images was compared to determine the relationship between the two approaches (Cvk and CvkA). To achieve this, a sample of spatially distributed random points within the study area was produced in the ArcGIS 10.8 platform, and then its tool was applied to extract multi-values to points from both variables (NDVI and C-factor). To quantify the sensitivity of NDVI-derived CvkA values to biophysical variables, the influence of biophysical variables ( Table 2) on CvkA (adapted C-factor) in Equation (7) with regard to this NDVI was carried out.

Support Practice (P) Factor
P is a factor that represents the support practice to express the rate of soil loss with erosion control practices such as strip-cropping, terracing and contouring [18]. The P factor (Figure 2d) can be computed for the croplands of different slope classes. Depending on the slope grades, either erosion control practice can be used. However, the values of the P factor were not estimated as part of the equation due to the unavailability of data. Therefore, this study adopted the P factor values (Table 3) suggested by Shin and Noh [45].

Soil Erosion Probability Zones Delineation
For the identification and generation of the areas with a probability of soil erosion, maps were prepared using ArcGIS 10.8 version. Key factors including slope, rainfall intensity, LULC and soil properties were utilized as the most influential factors on soil erosion. Weights were assigned to each theme taking into consideration its role in soil erosion. A raster overlay analysis process (weighted index overlay) was used to generate the zones with the highest probability of soil erosion occurrence. Therefore, the feature with the highest susceptibility was given the maximum value, while the minimum value was assigned to the feature with the lowest susceptibility.

Adapted Approach (C vk A) for C-Factor Estimation
By comparing the C-factor values from both methods (Cvk and its adapted version CvkA), the results exhibited Cvk with higher overestimations for the mean value in different LULCs (forestland, cropland, grassland and wetland) than CvkA in all the periods studied ( Figure 3). The results estimated from the adapted C-factor (CvkA) presented greater consistent values compared with those of the Cvk approach, which makes the estimation of soil loss more reliable. Figure 4 shows the spatial distribution of C-factor values estimated from two approaches for the LULC in the western province of Rwanda.    Figure 4). Cvk was found not to be suitable for estimating the C-factor in the study area, as it presented overestimated values compared to those from its adapted equation. CvkA was the most suitable based on the spatial variability of vegetation coverage of the area and the obtained C-factor values.
In general, NDVI is one of the several numerical combinations of satellite bands reported to be sensitive indicators in the existence and state of green vegetation considering the LULC of an area. For this, based on the spatial variability of soil coverage in the study area, the method using NDVI to obtain the C-factor was the most suitable. It was applied with an assumption that there was a linear relationship between NDVI and the C-factor linear model ( Figure 5). The two variables have shown a tendency to move in the opposite direction from one another, such that when NDVI increased, C-factor decreased, and vice versa. The study achieved the highest correlation coefficient between the NDVI and Cfactors in all the considered periods, with an R 2 higher than 85% ( Figure 5). Thus, NDVI was judged to be relevant in the estimation of the C-factor in the adapted equation.  Figure 4). Cvk was found not to be suitable for estimating the C-factor in the study area, as it presented overestimated values compared to those from its adapted equation. CvkA was the most suitable based on the spatial variability of vegetation coverage of the area and the obtained C-factor values.
In general, NDVI is one of the several numerical combinations of satellite bands reported to be sensitive indicators in the existence and state of green vegetation considering the LULC of an area. For this, based on the spatial variability of soil coverage in the study area, the method using NDVI to obtain the C-factor was the most suitable. It was applied with an assumption that there was a linear relationship between NDVI and the C-factor linear model ( Figure 5). The two variables have shown a tendency to move in the opposite direction from one another, such that when NDVI increased, C-factor decreased, and vice versa. The study achieved the highest correlation coefficient between the NDVI and C-factors in all the considered periods, with an R 2 higher than 85% ( Figure 5). Thus, NDVI was judged to be relevant in the estimation of the C-factor in the adapted equation.

Influence of Soil Erodibility on NDVI-Derived C-Factor (CvkA)
The executed regression analysis revealed that the vegetation index-derived C-factor values were affected by the soil erodibility (K) values ( Figure 6). The correlation analysis varied from 0.77 to 0.86. The sensitivity of the adapted C-factor estimation (CvkA) to the variation of soil background can be clarified through the spatial variability of soil erodibility (K) values in the study area. In the current study, an increase in the value of soil erodibility led to an increase in the values of NDVI-derived C-factor (Figure 6), while the magnitude probably changed following different seasons of a year. A higher value in the K-factor resulted in a higher near-infrared reflectance (NIR) band, while the lower the NIR reflectance, the lower NDVI values. This fact can explain the possible overestimation of the values in the prediction of soil erosion risk using the Cvk method under certain climate conditions, which may generate different K-factor values that compound the impact of the K-factor and NDVI-derived C-factor values in the RUSLE model ( Figure 8).

Influence of Soil Erodibility on NDVI-Derived C-Factor (CvkA)
The executed regression analysis revealed that the vegetation index-derived C-factor values were affected by the soil erodibility (K) values ( Figure 6). The correlation analysis varied from 0.77 to 0.86. The sensitivity of the adapted C-factor estimation (CvkA) to the variation of soil background can be clarified through the spatial variability of soil erodibility (K) values in the study area. In the current study, an increase in the value of soil erodibility led to an increase in the values of NDVI-derived C-factor (Figure 6), while the magnitude probably changed following different seasons of a year. A higher value in the K-factor resulted in a higher near-infrared reflectance (NIR) band, while the lower the NIR reflectance, the lower NDVI values. This fact can explain the possible overestimation of the values in the prediction of soil erosion risk using the Cvk method under certain climate conditions, which may generate different K-factor values that compound the impact of the K-factor and NDVI-derived C-factor values in the RUSLE model ( Figure 8). Land 2021, 10, x FOR PEER REVIEW 13 of 25

Influence of Rainfall on NDVI-Derived C-Factor (CvkA)
The correlation between annual mean NDVI and rainfall variable is shown in Figure  7. The studied period (2000-2018) clearly showed the geographical distribution of the relationship between NDVI and rainfall. The results revealed NDVI to be closely interacting with rainfall erosivity (Figure 7 (Figure 7). When the annual rainfall erosivity increased, the values of NDVI also increased at the same time. The analyses proved that more rainfall would definitely lead to higher NDVI value, which might affect the estimation of C-factor for accurate soil erosion predictions. The variation in annual NDVI was determined by climate variable including rainfall. The correlation between NDVI and rainfall is stronger due to the fact that rainfall is the dominant factor influencing vegetation growth in the area.

Influence of Rainfall on NDVI-Derived C-Factor (CvkA)
The correlation between annual mean NDVI and rainfall variable is shown in Figure 7. The studied period (2000-2018) clearly showed the geographical distribution of the relationship between NDVI and rainfall. The results revealed NDVI to be closely interacting with rainfall erosivity (Figure 7 (Figure 7). When the annual rainfall erosivity increased, the values of NDVI also increased at the same time. The analyses proved that more rainfall would definitely lead to higher NDVI value, which might affect the estimation of C-factor for accurate soil erosion predictions. The variation in annual NDVI was determined by climate variable including rainfall. The correlation between NDVI and rainfall is stronger due to the fact that rainfall is the dominant factor influencing vegetation growth in the area.

Influence of Topographic Features on NDVI-derived C-Factor (CvkA)
The results showed that NDVI-derived CvkA values had a significant impact on the varying elevation of the area (Figure 8)

Influence of Topographic Features on NDVI-Derived C-Factor (C vk A)
The results showed that NDVI-derived C vk A values had a significant impact on the varying elevation of the area (Figure 8

Soil Loss Rate Estimation
To facilitate the analysis and identify soil erosion severity classes that are of high priority for conservation practices, the estimated soil erosion maps ( Figure 9) were classified into six categories. From total surface area (5,883 km 2 ), erosion-prone areas were estimated at 4,670.8 km 2 (79.4%). Based on the LULC, the remaining 20.6% in each was occupied by non-erosive lands. The results revealed a mean soil loss of 15. An estimated soil loss rate of ≤ 10 t ha −1 yr −1 was considered as a threshold for soil erosion rate tolerance; however, an area that experiences a soil erosion rate of > 10 t ha −1 yr −1 occupied by the area with high erosion susceptibility requires urgent soil erosion control measures (Figure 9). Analyses showed that an area of about 2,684. 16 (Figure 9 and Table 4). However, our study found a high increase in annual mean soil loss by 68.

Soil Loss Rate Estimation
To facilitate the analysis and identify soil erosion severity classes that are of high priority for conservation practices, the estimated soil erosion maps ( Figure 9) were classified into six categories. From total surface area (5883 km 2 ), erosion-prone areas were estimated at 4670.8 km 2 (79.4%). Based on the LULC, the remaining 20.6% in each was occupied by non-erosive lands.   (Table 3). High, very high and extremely high An estimated soil loss rate of ≤10 t ha −1 yr −1 was considered as a threshold for soil erosion rate tolerance; however, an area that experiences a soil erosion rate of >10 t ha −1 yr −1 occupied by the area with high erosion susceptibility requires urgent soil erosion control measures (Figure 9). Analyses showed that an area of about 2684. 16 (Figure 9 and Table 4). However, our study found a high increase in annual mean soil loss by 68.  A soil erosion severity map was developed from the estimated soil erosion of the western province of Rwanda ( Figure 9). The result showed that about 40.7% (1900.12 km 2 ), 45.4% (2119.3 km 2 ) and 59.2% (2764.44 km 2 ) of the total erosive land highly was in the very low category (0-5 t ha −1 yr −1 ) of soil erosion risk in 2000, 2010 and 2018, which contributed to about 0.64%, 0.58% and 0.45% of the mean annual estimated soil loss, respectively. These results also show that only about 30% (1397.3 km 2 ) and 29% (1350.06 km 2 ) of the total erosive land was mostly covered by the category of high soil erosion rate (20-50 t ha −1 yr −1 ) in 2005 and 2015, respectively (Table 3). High, very high and extremely high soil loss categories experienced increases of 13.8%, 12.4% and 12.42% in terms of total erosive land from 2000 to 2005 (Figure 9a,b), respectively, while increases of about 16.77%, 15.7 % and 21.06% were recorded from 2010 to 2015 (Figure 9c,d), respectively. Therefore, the high (20-50 t ha −1 yr −1 ), very high (50-100 t ha −1 yr −1 ) and extremely high (>100 t ha −1 yr −1 ) erosion categories were found mainly in a small part of the Rusizi district (Bweyeye and Butare sectors), Nyamasheke district (Karambi, Mahembe, Cyato, Rangiro and Bushekeri sectors), Ngororero district (Muhanda sector) and in Nyabihu district (in Muringa Sector and north-east part of Kabatwa sector (Figure 9b,d,e). Commonly, this was mainly influenced by steep slopes in the area, rainfall intensity and land-cover management factors.

LULC Change and Its Impacts on Soil Erosion Rate
Classification and statistical results are presented in Figure 10 and Table 4, while estimated soil erosion rates per each LULC are shown in Table 5 these values were generated at the expense of a high decrease in forest cover (Table 5).

Estimated Mean Soil Loss for Each LULC Category
Comparing the soil erosion based on land-use types, the results showed a significant increase in the estimated soil erosion rates in cropland in 2005 and 2015 (Table 6) , the predicted mean soil erosion in cropland is more than two times (15.6 t/ha/yr to 34.3 t/ha/yr) and more than four times (10.4 t/ha/yr to 45.41 t/ha/yr) higher than mean soil erosion rates, respectively (Table 6). Water bodies are considered as non-erosive land. The  Considering the entire surface area, the net annual forest loss was 53% from 2000 to 2005, and 43% from 2010 to 2015, respectively ( Table 5). The built-up area had less significant land cover, but steadily increased from 0.43% to 0.6 and 0.57 to 2.39% between  (Table 6). Therefore, there was a mean increase of 18.7 t/ha/yr (about 120%) and 35.01 t/ha/yr (about 337%) from 2000 to 2005 and 2010 to 2015 in cropland, respectively; these values were generated at the expense of a high decrease in forest cover (Table 5).

Estimated Mean Soil Loss for Each LULC Category
Comparing the soil erosion based on land-use types, the results showed a significant increase in the estimated soil erosion rates in cropland in 2005 and 2015 (Table 6) , the predicted mean soil erosion in cropland is more than two times (15.6 t/ha/yr to 34.3 t/ha/yr) and more than four times (10.4 t/ha/yr to 45.41 t/ha/yr) higher than mean soil erosion rates, respectively (Table 6). Water bodies are considered as non-erosive land. The forestland and grassland are expected to have low erosion rates, but the results show that meant soil loss increased. Therefore, the poor coverage of trees due to a decline in forestland in the steep slopes accelerates soil erosion rates in the area.

Soil Erosion Probability Zones
The zones with a probability of soil erosion have been classified into six categories, including very low, low, moderate, high, very high and extremely high risk. It was observed that approximately 44.9% of the total eroded land produced a very high mean soil erosion of 697.6 t ha −1 yr −1 (Figure 11), while zones with high probability occupied 42.9% of the total eroded land, with 315 t ha −1 yr −1 mean soil loss. That said, zones of extremely high probability were found to cover 5.03% of the total eroded land, producing a mean soil loss of 1164.6 t ha −1 yr −1 . Therefore, the three smallest probability zones were very low, low and moderate, which covered about 0.46%, 2.10% and 4.5% of the total eroded area and produced a mean soil loss of 3.47 t ha −1 yr −1 , 32.79 t ha −1 yr −1 and 75.8 t ha −1 yr −1 , respectively. erosion of 697.6 t ha -1 yr -1 (Figure 11), while zones with high probability occupied 42.9% of the total eroded land, with 315 t ha -1 yr -1 mean soil loss. That said, zones of extremely high probability were found to cover 5.03% of the total eroded land, producing a mean soil loss of 1164.6 t ha -1 yr -1 . Therefore, the three smallest probability zones were very low, low and moderate, which covered about 0.46%, 2.10% and 4.5% of the total eroded area and produced a mean soil loss of 3.47 t ha -1 yr -1 , 32.79 t ha -1 yr -1 and 75.8 t ha -1 yr -1 , respectively. Figure 11. Zones with probability of soil erosion.

Cover Management (C) Factor and Biophysical Variables
The change in soil loss prediction can be reflected by NDVI with some variation in the values of cover management practices. The C-factor is one of the decisive parameters of the RUSLE, preventing bias from any C-factor estimation approach, which affects the prediction of soil erosion rate estimations [46]. For this, estimating the C-factor using NDVI by considering biophysical variables should be taken into account, as the applied equation has an influence on the estimated values. For this reason, the NDVI used to estimate C-factor was analyzed by checking its sensitivity to some biophysical factors such as rainfall, elevation and soil erodibility. The results revealed a good agreement between both variables (NDVI and biophysical factors) with acceptable and good correlation coefficients ( Figures 5-8). This is in agreement with a previous study [42] reporting that the NDVI-derived C-factor used to estimate soil erosion can be sensitive to different biophysical variables, such as soil condition, topographic features and vegetation phenology. In this sense, Jin and Zhang [47] argued that the spatial distribution of vegetation cover through NDVI was usually affected by elevation, concluding that elevation was the main controlling factor in vegetation growth in their study area. The result of this study revealed that variation in annual NDVI was determined by climate variable including rainfall; this is aligned with an earlier study [48] that reported that various climate factors, including rainfall, played a vital role in vegetation growth. Moreover, topographic reasons also conformed with an earlier study [19] that argued that highly mountainous and elevated sites are susceptible to soil erosion due to the domination of high inclination slopes and heavy rainfall distribution. Ayalew and Deumlich [42] revealed that a higher erodibility condition was associated with higher C-factor values which, in turn, might affect the prediction of soil erosion. To this end, this study found that the Van der Knijff and Jones [44] approach (Cvk) was not fit to estimate the C-factor in the study area, as it overestimated the values [18]. This can be clarified by the fact that this approach assessed the monthly C-factor based on NDVI and then calibrated the factor in a European study area that was totally different from the current study area (tropical) in terms of climate conditions and vegetation features. This is in agreement with Almagro and Thomé [43], who found underestimated Cvk values and concluded that this method was not suitable for tropical regions such as Rwanda; thus, this approach should be changed to CvkA. Generally, this study adapted the aforementioned equation by also performing a correlation analysis between the used NDVI and the produced C-factor ( Figure 5), leading to a good result in the prediction of soil loss in the study area without overestimation.

Variation in Land Use Land Cover and Soil Erosion
This study found that the main reason behind an accelerated soil erosion rate is attributable to changes in LULC ( Figure 10). The changes were depicted by the expansion of agricultural areas resulting from the decrease in forestland, steep slopes and high annual rainfall, which are also considered to be core soil erosion risk factors. Compared with different parts of the world, studies [49][50][51] have disclosed different rates of mean soil erosion under different land uses. The same conclusion was reached for the average annual soil losses following the fluctuation of land-use types triggered by changing rainfall patterns, while the highest estimated average annual soil loss per LULC has been accredited to seasonal agricultural land [52]. Moreover, based on global circulation models (GCMs) data, projected future storm rainfall and soil erosion rates in a monsoon-dominated region of eastern India showed that these regions were challenged by relatively greater erosion problems triggered by intense precipitation. They concluded that the loss of soils during erosion processes leads to decline in soil fertility, and thus the deterioration of soil resources in the region; this consequence should be expected in the study area examined in this paper. All the aforementioned drivers could be accelerated by LULC change, as the area showed dynamism. In the study area, the LULC change affecting soil erosion can be depicted clearly in the areas where a decrease in forestland and expansion of cropland areas resulted in an increase in the mean soil erosion of cropland cover from 2000 to 2005 and from 2010 to 2015 (Table 6). Previous studies found that the LULC change affects soil erosion rates by producing higher mean soil loss [18,25,30]. In the temperate forest zone of Russia characterized by multifaceted levels of land-use history, Zhidkin and Fomicheva [53] compared the simulated rates of soil erosion given by two models based on the type of precipitation (rainfall and snow erosion). Concerning rainfall acting as a trigger and the same types of land-use cover as examined in this present study, they found that variations in arable land (cropland) area and cropping structure were factors that determined changes in soil erosion to a greater extent, while other factors had no directional trend. This latter point is in accordance with the findings of this study, where the risk of erosion increases when land use changes from forestland in 2000 and 2010 to cropland in 2005 and 2015, respectively ( Figure 10). The features of each LULC category could be considered in changing the average soil loss estimated from one land-use class to another [54,55]. Similar to this study, some LULC types such as cropland and built-up land increased, while others, including forestland, decreased ( Figure 10). The change in forestland to uncontrolled agricultural lands is the main cause blamed for the increase in soil erosion, with triggers such as heavy rainfall and the length and steepness of slopes ( Figure 9); these results are in accordance with a study [52] conducted in Kenya located in the same region as our case study disclosing that the increment of slope steepness and its length (LS-factor) increases the level of soil erosion due to the high velocity and erosivity of runoff. Generally, in rural areas of the western province, a large number of the residents depends on agriculture as a means of livelihood [56]. Land scarcity for agriculture and the continuous increase in demand for agricultural commodities are the main reasons pushing communities to convert forest and grasslands into cultivated lands [56]. Unfortunately, unprecedented climate changes are expected in the region, and the most prominent climatic change impact on rain-fed agriculture results in a reduction in crop productivity [57,58]. Thus, areas with a considerable conversion from forest to less vegetated areas, coupled with steep slopes and rainfall intensity, in the study area should expect higher soil erosion potential [25,[59][60][61]. Table 4 shows the estimated soil erosion intensity classes and the mean soil erosion rate of a certain magnitude of soil loss and its area coverage over the previous two decades. In the entire study area, soil erosion severity was classified from very low to extremely high. These classes were concentrated in cultivated land and steep slope areas, which indicated that areas in these categories were highly susceptible to soil erosion. These findings are in accordance with a previous study [62] providing details on the intensity of soil erosion in different slope gradients, whereby the analysis of soil erosion intensity exhibited a slightly intense erosion on largely gentle (<5 • ) slope terrains and steep terrains (15-25 • ), while severe soil erosion was found to be on moderate (5-15 • ) slopes. However, the very low and low soil erosion classes fell under a <10 t ha −1 yr −1 severity class, making them less susceptible to soil erosion based on slope gradient [25,30]. This study also indicates a significant increase in high, very high and extremely high soil erosion severity classes of about 13.85, 12.3% and 12.42%, respectively, between 2000 and 2005. Increases of about 16.77%, 15.7% and 21.06%, were found in high, very high and extremely high soil erosion severity class, respectively, from 2010 to 2015, accelerated by the expansion of cropland in steep slopes. Nevertheless, a small decline in three severity classes (very low, low, moderate) from 2005 to 2010 and from 2015 to 2018 resulted from various conservation activities implemented by local leaders and sector land officers in the study area for sustainable land-use management [25,29]. The activities include the implementation of bench terraces on steep cultivated land and afforestation campaigns. While this optimistic initiative to lessen the soil loss caused by erosion has proven popular, problems related to clearing natural forests and unsuitable cultivation on steep slopes still exist. Therefore, rainfall intensity and forest clearance become core factors for increasing zones prone to soil erosion within the study area. It was also demonstrated that about 44.9 % of the study area fell under very high soil erosion risk, while high-probability zones covered 42.9% of the study areas ( Figure 11). Compared with an earlier study [63] in the Urmodi river watershed (India), a significant output regarding increased conversion was obtained from very slight to very severe risk zone; research concluded the existence of extensive variation in very severe soil vulnerable risk zones over the study period. Overall, the findings of this study and other studies carried out in Rwanda and east African regions similar to results reported by Karamage [30,62], which seem to overestimate soil loss rates.

Estimation of Soil Loss and Probability Zones
The conversion of lower to higher risk category areas of soil erosion was detected all over the study area. Hence, croplands and grasslands, with their related conversions as the main LULC classes, should lead the immediate execution of more suitable agricultural practices; terracing integrated with surface residue mulching and cover-cropping and crossslope farming should be suggested for reducing soil erosion potential within agricultural areas in the western province of Rwanda.

Conclusions
In this study, an adapted approach (CvkA) for the estimation of cover management (C) factor from RUSLE was explored in comparison with the formerly used approach (Cvk) based on NDVI. We also analyzed the correlation between the adapted C-factor (CvkA), NDVI and different biophysical variables such as the rainfall, elevation and soil erodibility. The Cvk approach presented an overestimation of the C-factor values, making it unsuitable for the study area. Therefore, using the adapted C-factor, which was found to be suitable for this study, we analyzed the impacts of LULC dynamics on soil erosion risk in agricultural areas of the western province of Rwanda. The results showed significant correlations between the NDVI and the biophysical variables, which implied that using remote sensing data to capture the spatial and temporal changes in C-factor estimation can serve as an input factor for soil erosion modeling procedures that can themselves be enhanced by considering the quantified sensitivity of NDVI-based C-factor estimations.
Moreover, the results revealed that the area experienced estimated annual mean soil losses of about 15. . It was also observed that the rate of erosion varied based on steep slopes and LULC changes. The noted differences in the studied period indicated that LULC changes affected the rate of soil erosion significantly. The application of the weighted index overlay (WIO) method supported the categorizing of areas into different probable zones of soil erosion risk, which was useful for developing sustainable soil conservation measures. The estimated amount of soil loss suggests the urgent implementation of conservation practices such as terracing, integrated with surface residue mulching and cover cropping and cross-slope farming in the study area. The results of this study may be useful for soil and water conservation planning and soil loss estimations. However, the application of the adapted C-factor (CvkA) approach can be used in similar tropical regions after evaluating its performance. Future research should attempt to compare NDVI with other spectral indices such as the soil adjusted vegetation index and enhanced vegetation index in the context of C-factor estimation.