1. Introduction

Influenced by the climate and land use changes, soil erosion was identified as one of the major threats to the World’s Soil Resources [1] [2] which affects agricultural productivity [3]. It is accelerated by human activities among which agriculture is the primary cause of soil degradation worldwide [4]. Soil erosion control is one option to increase crop productivity while controlling river and lake sedimentation [5].

This research intends to apply Universal Soil Loss Equation model in order to make recommendations that will build sustainability into soil erosion management in Sebeya catchment located in Western Province of Rwanda. Various researches [6] defined sustainable development as the development that meets the needs of the present without compromising the ability of future generations to meet their own needs. All USLE-type models such as USLE, Modified Universal Soil Loss Equation (MUSLE) and Revised Universal Soil Loss Equation (RUSLE) need information related with soil type, land use, landform, climate and topography to estimate soil loss. They are designed for a specific set of conditions in a particular area [2].

Presently, USLE is the most widely used model to predict soil erosion rates [7] but its applicability is limited to sheet and rill erosion without any consideration of gully erosion [8]. Soil erosion management can be facilitated by using simulation and modelling. Modelling is a useful tool for soil erosion scenario assessment that enables the adequate selection of soil erosion control measures [9] [10]. Soil erosion models can be understood as a virtual laboratory that brings together all data, observations and knowledge of different fields for sustainable environment [11].

The main objective of this study was to develop a USLE-type erosion model to be used in predicting soil loss and correlatively the crop yields at the level of parcels for sustainable agriculture management in Sebeya catchment. The specific objectives of this research were: 1) to estimate the actual soil erosion rates in Sebeya catchment; 2) to propose suitable best management practices for soil erosion control in each parcel of Sebeya catchment; 3) to correlatively establish relationships between soil erosion rates and soil productivity at the parcel level and finally; 4) to develop a soil erosion control model which will help farmers in Sebeya catchment to improve their farming system.

2. Methodology

2.1. Study Area

This study was conducted on Sebeya catchment located in the Western Province of Rwanda and shared by four administrative units namely Rubavu, Nyabihu, Rutsiro and Ngororero Districts (Figure 1). On a total area of 363.1 km², the average population density of Sebeya catchment is estimated to 644 hab/km² while the average population density of Rwanda is about 415 habitants per km² [12] [13]. The soil in this catchment favors agriculture due to its high infiltration rates and its high minerals content. Located in the high elevation region of the country with altitude varying between 1462 m to 2979 m a.b.s.l. (meters above sea level), this catchment is also characterised by steep slopes and abundant rainfall varying between 1200 mm to 1700 mm per year [13].

2.2. Data and Methods

2.2.1. Model Choice and Suitability

In choosing the soil erosion model to be applied, it is always necessary to know the availability of input data and the type of output data needed [14]. Modelling cannot be an alternative to measurement and monitoring but might be a powerful tool for simulation and prediction of the soil erosion potential [15]. Soil erosion models can be understood as a virtual laboratory which brings together all data, observations and knowledge of different fields [11].

This research aims to apply USLE-type model such as USLE, MUSLE or RUSLE. Presently, USLE is the most widely accepted method for estimating the annual soil loss because of less data demand [8]. It has been originally proposed in 1958 and actually modified to its present form in 1978 [16]. The empirical equation of this model is as follows:

$A=K\ast R\ast \left(LS\right)\ast C\ast P$ (1)

where the involved parameters are: Average annual soil loss A (t/ha * year); Soil

erodibility or K-factor (t * ha * hr/ha * MJ * mm); Rainfall erosivity or R-factor (MJ * mm/ha * hr * yr); Slope length factor or LS-factor (Dimensionless); Crop management factor or C-factor (Dimensionless); Erosion-control practice factor or P-factor (Dimensionless). As described in the next section, the designed Soil Loss and Crop Yield (SOLCY) model is based on the USLE concept.

2.2.2. SOLCY Model Design

From the Center of GIS of University of Rwanda (UR-CGIS) parcels data, a parcels map and an Excel sheet containing 259,673 parcels each of them having a specific location (District, Sector, Cell and Village) were first made available. Each parcel is characterised by a Universal Parcel Identification code abbreviated as parcel UPI. The shapes of parcels are diversified: triangular, square, rectangular and many irregular shapes.

At the large-scale application of USLE model, the catchment area has been discretized into a series of 259,673 parcels having relatively independent land use. Figure 2 shows diversified processes used by Soil and Crop Yield (SOLCY) mode designed to predict soil loss and correlatively estimated crop yields.

2.2.3. Rainfall Data

The collected data from Rwanda Meteorological Agency were used to generate rain fall distribution ranges in Sebeya catchment using ArcGIS software tool. Knowing the minimum and maximum values for each range, the reading of average annually precipitation on each parcel in Sebeya catchment was made possible by overlapping the administrative boundaries map and the parcels map of Sebeya catchment on the obtained rainfall map.

2.2.4. Soil Texture Data

In this present study, a shapefile of soil was obtained from the Center of GIS of University of Rwanda (UR-CGIS) to produce soil texture map. To determine the soil texture for each parcel within Sebeya catchment, the administrative boundaries map and Sebeya parcels map were overlapped to soil texture map.

2.2.5. Topographic Data

The Digital Elevation Model data collected from the Center of GIS of University of Rwanda (UR-CGIS) were used to generate the slope range map of Sebeya catchment. In order to detect the slope range of each parcel, the administrative boundaries map and parcels map were overlapped on the slope range map of Sebeya catchment.

2.2.6. Management Data

Terraces are more favorable in agricultural land with steep slopes ranging from 16% to 40% while progressive terraces and contour bunds are suited to flat areas of slopes less than 16% [17]. Erosion control measures in Sebeya catchment were proposed based on land slope range [17].

2.2.7. Land Cover and Land Use Data

The data collected at UR-CGIS revealed that Sebeya catchment is characterized by seven land cover types: Built-up areas, closed agriculture, forest plantation, irrigation, natural forest, open agriculture and open land. In the determination of the type of land cover and land use on each parcel within Sebeya catchment, the administrative boundaries map and Sebeya parcels map were overlapped on the land cover map.

3. Results and Discussions

3.1. Actual Status of Soil Erosion in Sebeya Catchment

Soil erosion within Sebeya catchment was categorized into 6 classes according to the estimated total soil loss in t/ha/yr which are 0 - 5 very low, 5 - 10 low, 10 - 25 moderate, 25 - 50 high, 50 - 100 very high and greater than 100 extremely high. Around 8000 ha are under high risk, around 6000 ha under very high risk while around 4000 ha are under extremely high risk of soil erosion [18]. In determining the actual soil loss in Sebeya catchment, all USLE parameters were mapped and maps integrated in GIS. The annual estimated soil loss was 130.724 t/ha/yr. In Sebeya catchment, the high risk of soil erosion results from improper management of land, heavy rainfall and human activities that disturb the soil.

3.2. Values of USLE Parameters for Each Parcel in Sebeya Catchment

3.2.1. R-Factor Determination

R-factor is the long term annual average of the product of rainfall kinetic energy (KE) in MJ∙ha⁻¹ and the maximum rainfall intensity in 30 minutes (I₃₀) in mm∙h⁻¹ [19]. The rainfall factor can vary from year to year, so an average over a number of years is usually used [20]. The average precipitation on each parcel was calculated based on the minimum and maximum precipitations found on the rainfall range map as:

$\bar{P}=\left(P_{\min }+P_{\max }\right)/2$ (2)

The observed average rainfall values were varying from 1132.5 mm to 1568 mm and were inserted in the SOLCY Excel spreadsheet with respect to each parcel UPI. Due to lack of rainfall kinetic energy and intensity data, which is the case of this study, many researchers [21] have predicted the erosivity factor R by using the following regression equation:

$R=81.5+0.38\bar{P}$ (3)

where: R = rainfall erosivity in MJ∙mm∙h⁻¹∙year⁻¹; $\bar{P}$ = average precipitation on the parcel in mm.

Applicable for areas where annual precipitation ranges from 340 and 3500 mm, Equation (3) was selected because its range of validity includes all precipitation events of this study. By inserting all calculated values for all parcels in the SOLCY Excel sheet, the R-factor in Sebeya catchment ranges from 511.85 to 677.34 MJ∙mm∙ha⁻¹∙h⁻¹∙year⁻¹.

3.2.2. K-Factor Determination

The soil erodibility factor (K) indicates the susceptibility of soil to erosion where it reflects the effect of soil properties and soil profile characteristics on soil loss [8]. The soil erodibility is the resistance of the soil to both detachment and transportation. It can be determined through the measurement of soil loss from a standard runoff plot (called “Standard USLE Plot”). A “Standard USLE Plot” is a field experimental plot having 9% slope along 22.13 m length which is kept fallow (bare soil) with periodic tillage [22]. In this case, LS = C = P = 1, the soil loss becomes a function of R and K. Then from a known R and soil loss measurement, K can be estimated. The soil erodibility factor (K) rates from 0 to 1, where 0 indicates soils with the least susceptibility to erosion and 1 for soils that are highly susceptible to soil erosion by water. From literature, different researchers have identified the variability of the soil erodibility values with the soil texture (Table 1). The soil erodibility factor is determined by soil texture and organic matter content [23]. The determination of K-factors for all parcels within Sebeya catchment was determined by reading in Table 1 and inserted in the SOLCY Excel sheet.

3.2.3. LS-Factor Determination

The topographic factor LS is defined as the ratio of soil loss under the given conditions to that at the site with the “standard” slope steepness of 9% and the slope length of 22.13 m [29]. The topographic erosivity factor LS accounts for the effect of slope length (L) and slope steepness (S) on soil erosion. The steeper and longer the slope is, the higher is the risk of soil erosion. A higher slope gradient creates a higher flow velocity which causes more detachment and transport of soil particles [30]. In 1978, Wischmeier and Smith defined the slope length (L) as: “the distance from the point of origin of the surface flow to the point where each slope gradient (S) decreases enough for the beginning of deposition or when the flow comes to concentrate in a defined channel”. Further researches extend the LS-factor to topographically complex units using a method that incorporates contributing area and flow accumulation [31]. Generally, to come up with the result of a slope length factor of each parcel, the following equation has been found appropriate to be used [32]: For the slope less than 20%:

$LS=\left(L\right)0.5\times \left(0.0138+0.00965S+0.00138S^2\right)$ (4a)

where: L = slope length (m) and S = slope gradient (%).

And when the slope is greater than 20%:

$LS=\left(L/22.2\right)0.6\times \left(S/9\right)1.4$ (4b)

Another approximation to this process was to take the maximum slope of each slope range as the slope steepness while the largest path of travel P/2 was attributed to the slope length (L) for each parcel in Sebeya catchment where P is the parcel perimeter. Finally, Equation (4) was applied to estimate the LS-factors for all parcels of Sebeya catchment and the obtained values were inserted in the SOLCY Excel sheet. For all parcels, the minimum and maximum observed slope lengths (and slope length factors) are 0.014 m (0.032) and 24707.00 m (559.798) respectively while the following Table 2 shows the percentages of areas covered in each slope range.

*: weighted average.

3.2.4. C-Factor Determination

The cropping-management factor can vary according to farming practices. In USLE-type equations, C-factor measures the combined effects of all interrelated cover and management measures. Its value includes the effects of crop cover, crop sequence, and length of growing season, tillage practices, residue management and the expected time distribution of erosive rainstorms [23]. It is expressed as the ratio of the soil loss from a cropped land under specified field conditions and the soil loss from a continuously clean-tilled bare soil surface over a given period of time. One way to reduce soil loss is to choose a crop that provides cover during the year when rainfall is most erosive.

Table 3 indicates various C-factors for 13 main crops selected with high suitability of growing in Sebeya catchment [33]. All C-factors for 13 main cultivated crops in Sebeya catchment were inserted in SOLCY Excel sheet for all 195,533 arable parcels.

3.2.5. P-Factor Determination

The P-factor represents the effect of various support practices such as contour farming, terracing and strip cropping for arresting soil erosion being taken up in the area [9] [27]. It is the ratio of soil loss from a land where conservation practices (like: contouring, strip cropping, terracing, etc.) are adopted to soil loss from a land where soil erosion control measures were not implemented. The value of P-factor ranges from 0 to 1. The value of P-factor for a water body is 0 while the value of P-factor is taken as 1 for no conservation practice [34]. From literature, Table 4 indicates typical values of P-factors of various soil erosion control measures.

Various combinations of soil erosion control measures to be implemented in Sebeya catchment have been proposed based on the land slope range [17]. As simulated by SOLCY model, Table 5 summarizes various combinations of soil erosion control measures applied on each parcel of Sebeya catchment to predict soil loss and crop yields. It is recommended that biological measures such as planting trees and grasses as well as lime and compost applications can be undertaken. The first step was to assign a slope range to each parcel before calculating P-factor corresponding to each proposed combination of the soil

erosion measures. The product of 2 or more P-factors gives a P-factor which represents the effects of the combined soil erosion control measures to be implemented. Table 4 gives all P-factors for all soil erosion control measures to be used in calculating the P-factors corresponding to various combinations of soil erosion control measures simulated in SOLCY model and presented in Table 5.

3.3. SOLCY Model Assembly in Excel Sheet

Section 2.2 of this research describes the methodology while Section 3.2 systematically shows the determination of all 5 USLE parameters for this study. At this stage, respective values of all USLE parameters have been inserted in SOLCY Excel sheet for all 259,673 parcels of Sebeya catchment.

3.3.1. Acceptable Maximum Soil Loss Tolerance Limit (T-Value)

In many developed countries, the best management practices (BMPs) for soil erosion management have been technologically implemented. These techniques helped greatly to meet the tolerable soil loss criterion [45]. T-value is a concept useful to judge if a soil has a potential risk of erosion, productivity loss and off-site damages as a river or reservoir sedimentation [46] [47]. The generally accepted maximum limit of soil loss or T-value is 11.5 t∙ha⁻¹∙year⁻¹ [16] [44] [48]. The aim of soil conservation strategies is to limit rilling using a threshold soil loss value. But the occurrence of soil loss equal to or lower than 11.5 t∙ha⁻¹∙year⁻¹ does not ensure absence of rills. The procedure used in assigning T-value has relied on multiple judgments of various researchers.

3.3.2. Predicting Soil Loss with SOLCY Model in Sebeya Catchment

In 2017, the National Institute of Statistics of Rwanda (NISR) 3 recognized agricultural seasons in Rwanda as: Season A (September-February); Season B (March-June) and Season C (July-September). Table 6 gives a typical crop rotation simulated for soil loss prediction in Sebeya catchment with SOLCY model.

SOLCY model is a USLE-type model at each parcel level:

$A=R\ast K\ast LS\ast C\ast P$

where R, K, LS are the already determined USLE parameters and now considered as physical constants in nature at each parcel. The only remaining variable is C-factor which will change according to the type of the seasonally cropping patterns. In fact, the simulated soil erosion control measures presented in Table 5 are assumed constant for any crop at any slope range.

The first part in developing this SOLCY model was to simulate how soil loss can be reduced if suitable combinations of soil erosion control measures are applied specifically to well-known site conditions at all parcels in Sebeya catchment. Comparatively by cultivating any type of the 13 main crops on any parcel, each annually average soil loss which can be predicted with SOLCY model will lay between 2 extreme values: minimum and maximum soil loss as A_min < A < A_max. The 1st step of simulation was done with the typical crop rotation described in Table 6, the 2nd step simulated 3 seasonally crops which have a minimum C-factor of C = 0.15 and the 3rd step simulated 3 seasonally crops which have a maximum C-factor of C = 0.55. Referring to the combinations of soil erosion control measures defined in Table 5, the results for these 3 steps are presented in Table 7.

3.3.3. Relationship between Soil Loss and Crops Yield

The 2nd part in developing SOLCY model was to apply successively a set of soil erosion control measures (Pi1, Pi2 and Pi3) as defined in Table 5 in order to predict the resulting soil loss and crop yield on each arable parcel in Sebeya catchment. Table 8 gives an estimate of the minimum and maximum crop yields for the 13 main cultivated crops in Sebeya catchment.

Many researchers have proposed linear relationships between soil loss and crop yields [50] [51]. In order to predict the crop yields, a linear relationship was assumed by considering the 1st point (x1 = the maximum soil loss; y1 = the minimum crop yield) and the 2nd point (x2 = the minimum soil loss; y2 = the maximum crop yield). From literature, the maximum soil loss has been fixed by previous researchers at 137 t/ha/yr [50] while the minimum soil loss can be taken as the maximum acceptable limit of soil loss or T-value of 11.5 t/ha/yr [44]. Having 2 points in plane XY Cartesian coordinates at this stage, a linear relationship becomes well defined for each crop in the form:

$Y=a\ast X+b$ (5)

where a & b are known constants for each crop yield prediction; X = A = soil loss and Y = crop yield.

3.3.4. SOLCY Model Printout

The soil loss for each of the 13 crops simulated on each arable parcel of Sebeya catchment and the value of the crop yield correlatively predicted using the above linear relationship (Equation (5)) were already simulated for each of the combination of soil erosion control measures Pi1, Pi2 or Pi3 as defined in Table 5 and now available in SOLCY Excel spreadsheet. The results as shown in Table 9 were extracted from SOLCY Excel spreadsheet by using the INDEX-MATCH type lookup formula. Taking an example of the parcel of UPI = 30,305,034,908, the user has to type only the parcel UPI and the soil loss with its corresponding crop yield are extracted for all 13 main cultivated crops in Sebeya catchment to be printed in Table 9.

3.3.5. Who Can Use SOLCY Model?

Farmers, District agriculture planners and technicians in Sebeya catchment can use SOLCY model. Most of the cases, an agriculture technician may be available at each sector and ready to facilitate farmers who will need SOLCY model guidance on how they can improve their farming systems. Also, researchers can develop similar models in other catchments based on SOLCY design concept. For 13 simulated crops, a District planner can take advantage of SOLCY to get advice in proposing a crop to be cultivated in a sector or a region based on its productivity and by suggesting BMPs for an appropriate farming system.

3.4. SOLCY Model Validation

Model validation is a process that verifies if the model is performing properly as expected according to its design objectives and intended uses. It also identifies potential limitations and assumptions, and assesses its possible impact [15].

Validation refers to the testing of the model output to confirm if the produced results reflect the reality [54].

3.4.1. Innovative Aspects of SOLCY Model

Initially, the estimation of soil loss using USLE model with GIS applications on Sebeya catchment was done at large-scale. In order to get more precise results and flexibility in assessing soil loss at parcel level, USLE model itself was limited to measure the soil loss at the outlet of the catchment but it was very difficult to print all the results of soil loss at the outlet of each of 259,673 parcels in Sebeya catchment. That’s why an improvement was proposed to apply analytically the USLE model on each parcel by simulating the cultivation of 13 main crops in Sebeya catchment with selected combinations of soil erosion control measures. A such designed SOLCY model is estimating the soil loss at each parcel and prints results of soil loss for all 259,673 parcels and results of the predicted crop yields for all 195,533 arable parcels of Sebeya catchment in an Excel sheet for each combination of soil erosion control measures.

3.4.2. Validating SOLCY Model with GIS Results

By simulating the combination Pi1 of soil erosion control measures in each parcel (Table 5), the soil loss from the entire catchment by SOLCY model was 143.27 t/ha/yr compared to 130.72 t/ha/yr of soil loss using GIS. In order to judge the significance of the difference between the means of these 2 populations, a sample of 30 parcels was selected randomly from 259,673 parcels of Sebeya catchment. Table 10 shows the estimates of soil loss using GIS in comparison with the predicted soil loss by SOLCY model.

In order to assess the significance of the difference for 2 sets of results as in the case of Table 10, One-Way ANOVA technique using F-test is usually applied [55]. Using an Excel tool for a significance level α = 5%, the F-test results obtained and shown in Table 11 allow us to draw the conclusion that the two samples have been extracted from populations having the same mean.

3.4.3. Limitations of SOLCY Model

Fertilizer has effects on soil loss and crop yield because it changes the organic matter of soil. However, its doses were not indicated in SOLCY model. Finally, SOLCY model is limited to Sebeya catchment but this study provides an appropriate methodology and procedures to be duplicated to other catchments.

4. Conclusion and Recommendations

The present study aimed to design a model that can help to improve farming system in Sebeya catchment. The 5 USLE parameters were determined on each parcel and inserted in an Excel sheet for all 259,673 parcels of Sebeya catchment. By applying a combination of 0, 1, 2 and 3 soil erosion control measures on each of 259,673 parcels, the simulated annual soil loss for the entire catchment were 849.94; 143.27; 88.64 and 28.59 t/ha/yr respectively. Correlatively on each of 195,533 arable parcels of Sebeya catchment, the crop yield was predicted from the soil loss value using an established linear relationship. The designed SOLCY model stands for Soil Loss and Crop Yield prediction for each of 13 main cultivated crops in Sebeya catchment. This research has proven high improvement in soil loss reduction if increasing the number of soil erosion control measures in the applied combinations on each of 259,673 parcels of Sebeya catchment. Farmers, District agriculture planners and technicians in Sebeya catchment can use SOLCY model to focus on the parcel soil loss and predicted crop yield for sustainable land management decisions. The development of soil erosion rates greater than the acceptable soil loss tolerance limit (T-value) should be prevented by the affected people and landholders who are capable of identifying rill erosion at its early stage through some adaptive measures based on their indigenous knowledge. Finally, SOLCY model is limited to Sebeya catchment but this study provides an appropriate methodology and procedures to be duplicated to other catchments.

Acknowledgements

The authors are grateful to the University of Rwanda (UR) for its support and wish to acknowledge the support of UR-Sweden Program in sponsoring this research and we cannot fail to commend and appreciate the work of various authors used for the present paper.