Comparison of RUSLE and MMF Soil Loss Models and Evaluation of Catchment Scale Best Management Practices for a Mountainous Watershed in India

Soil erosion from arable lands removes the top fertile soil layer (comprised of humus/organic matter) and therefore requires fertilizer application which affects the overall sustainability. Hence, determination of soil erosion from arable lands is crucial to planning conservation measures. A modeling approach is a suitable alternative to estimate soil loss in ungauged catchments. Soil erosion primarily depends on soil texture, structure, infiltration, topography, land uses, and other erosive forces like water and wind. By analyzing these parameters, coupled with geospatial tools, models can estimate storm wise and annual average soil losses. In this study, a hilly watershed called Nongpoh was considered with the objective of prioritizing critical erosion hazard areas within the micro-catchment based on average annual soil loss and land use and land cover and making appropriate management plans for the prioritized areas. Two soil erosion models namely Revised Universal Soil Loss Equation (RUSLE) and Modified Morgan–Morgan–Finney (MMF) models were used to estimate soil loss with the input parameters extracted from satellite information and automatic weather stations. The RUSLE and MMF models showed similar results in estimating soil loss, except the MMF model estimated 7.74% less soil loss than the RUSLE model from the watershed. The results also indicated that the study area is under severe erosion class, whereas agricultural land, open forest area, and scrubland were prioritized most erosion prone areas within the watershed. Based on prioritization, best management plans were developed at catchment scale for reducing soil loss. These findings and the methodology employed can be widely used in mountainous to hilly watersheds around the world for identifying best management practices (BMP).


Introduction
Soil loss resulting from erosion is a universal issue which affects agricultural production and natural resources [1,2]. Degradation of soil physical properties by soil erosion can affect the crop growth and yield by reducing the root depth, available water, and nutrient reserves, and also by affecting soil organic carbon, phosphorus, potassium, nitrogen contents, and pH [3]. It also carries soil-laden water downstream, which can produce heavy deposits of sediment that prevent the smooth flow of rivers and streams and can eventually lead to floods [2]. Soil erosion occurs when soil particles are carried off by water or wind and deposited elsewhere. Water erosion begins once rain or irrigation water detaches soil particles. Once there is an excessive amount of water on the soil surface, Given the difference in the model structures and their inherent limitations, it is critical to identify which soil loss model is more suitable for mountainous catchments in India. Therefore, in this study, a peri-urban watershed in high rainfall areas of Meghalaya, which is close to Nongpoh, district headquarters of Ri-Bhoi district, was taken for estimating erosion. The farmers are found to extend their cultivations in many other places within the study watershed. Even some crops like rice, pineapple, and tomato cultivation have been started in the foot hills areas of the watershed. As the study area is a peri-urban watershed with a high degree of slope and the site usually receives a high amount of rainfall, there is clearly pressure on the natural resources of watershed. Given the high slope aided with high rainfall combined with anthropogenic activities, the study area (and also the entire Meghalaya state) is losing the fertile topsoil layer contributing to the deterioration of soil health from an agricultural perspective [31,32]. Hence it has become paramount to quantify the soil loss and make the best management plan for maintaining the sustainability of the land and improving the agricultural production in the study area. Both RUSLE and MMF models were used for estimating soil erosion and identifying the erosion prone areas with the following objectives: (1) to prioritize the critical erosion hazard areas within the micro-catchment based on average annual soil loss, land use, and land cover, and (2) to make the appropriate management plans for the prioritized area. The findings of this study will aid in developing best management practice plans for the watershed for maintaining overall agricultural sustainability by improving soil health and reducing the application of fertilizers/manures for enhanced crop yields.

Description of the Study Area
The study area is a micro-catchment situated at Nongpoh of Ri-Bhoi district, Meghalaya. The study area lies between 25 • 54 0" to 25 • 55 05" N latitude and 91 • 52 54.7" to 91 • 54 17.7" E longitude and it covers an area of 268 ha (2.68 km 2 ). The study area falls under a humid subtropical climate. The area receives a mean annual rainfall of about 2700 mm, 77% of which is received from July to September [33]. The study area has an average elevation of 469-760 m above the mean sea level (MSL). The relief is moderate to steeply sloping and the drainage condition varies from well drained in upland areas to poorly drained in low lying areas. The location of the study area is shown in Figure 1. Additionally, Table 1 illustrates the mean monthly and average annual precipitation at four different rain gauges near the micro-catchment. The average precipitation values were derived for a period of 9 years, i.e., 2010-2019. Figure 2 showcases the severity of soil erosion in the study catchment.
The Nongpoh micro-catchment is classified into seven slope classes i.e., flat (0-3%), undulating (3-8%), moderately slope (8-15%), hilly 105 • C for 24 h and then samples were weighted, (15-30%), moderately steep slope (30-45%), steep slope (45-60%), and very steep slope (>60%) as per the recommendations of Lee and Silalahi et al. [34,35]. The slope map of Nongpoh micro-catchment is displayed in Figure 3a. Additionally, Figure 3b displays the hypsometric map of the study area. It can be seen that the contours representing higher elevation (ranging from 540 to 760) are more dominant in the edges of the catchment.        The observed soil loss data were obtained for the years 2010-2019 from an ongoing research project of College of Post Graduate Studies (CPGS), Central Agricultural University, Meghalaya, India. After each rainstorm, depth of runoff data were collected at the catchment outlet and, using a measuring cylinder of one liter of water samples, were taken for determination of soil loss. Post-collection of the samples, they were filtered through filter paper of 42 no. grade and the soil was separated from the water. Thereby, the soil was kept in an oven at 105 °C for 24 h and then samples were weighed.
The methodological flowchart is given in Figure 4. First, the digital elevation model (DEM) having a spatial resolution of 2.5 m derived from CARTOSAT-1 (downloaded from Bhuvan, https://bhuvan-app3.nrsc.gov.in/data/download/) was used in identifying the slope and catchment delineation. Additionally, the land use land covers (LULC) for the catchment was derived from the satellite image of Linear Imaging Self-Scanning  The observed soil loss data were obtained for the years 2010-2019 from an ongoing research project of College of Post Graduate Studies (CPGS), Central Agricultural University, Meghalaya, India. After each rainstorm, depth of runoff data were collected at the catchment outlet and, using a measuring cylinder of one liter of water samples, were taken for determination of soil loss. Post-collection of the samples, they were filtered through filter paper of 42 no. grade and the soil was separated from the water. Thereby, the soil was kept in an oven at 105 °C for 24 h and then samples were weighed.
The methodological flowchart is given in Figure 4. First, the digital elevation model (DEM) having a spatial resolution of 2.5 m derived from CARTOSAT-1 (downloaded from Bhuvan, https://bhuvan-app3.nrsc.gov.in/data/download/) was used in identifying the slope and catchment delineation. Additionally, the land use land covers (LULC) for the catchment was derived from the satellite image of Linear Imaging Self-Scanning The observed soil loss data were obtained for the years 2010-2019 from an ongoing research project of College of Post Graduate Studies (CPGS), Central Agricultural University, Meghalaya, India. After each rainstorm, depth of runoff data were collected at the catchment outlet and, using a measuring cylinder of one liter of water samples, were taken for determination of soil loss. Post-collection of the samples, they were filtered through filter paper of 42 no. grade and the soil was separated from the water. Thereby, the soil was kept in an oven at 105 • C for 24 h and then samples were weighed.
The methodological flowchart is given in Figure 4. First, the digital elevation model (DEM) having a spatial resolution of 2.5 m derived from CARTOSAT-1 (downloaded from Bhuvan, https://bhuvan-app3.nrsc.gov.in/data/download/) was used in identifying the slope and catchment delineation. Additionally, the land use land covers (LULC) for the catchment was derived from the satellite image of Linear Imaging Self-Scanning System IV (LISS-IV) for the year 2015. With the identified characteristics of the catchment and the Sustainability 2021, 13, 232 6 of 22 RUSLE and MMF model parameters derived within the catchment, the rate of soil loss for the years 2010 to 2019 year was calculated, which was further used in the model calibration and validation. Moreover, based on the estimated rate of soil erosion, the prioritized micro-catchments were chosen, and using ArcGIS 10.2, several best management practices for reducing soil loss at the micro-catchment were conducted.

Land Use Land Cover
Based on tonal and color variation in the satellite imagery (LISS-IV, 8-12-2015) and ground-truthing (supervised classification), the major LULC classes were identified. The spatial distribution of LULC categories, such as contour bunding, agriculture, dense forest, open forest, settlement, scrubland, terrace farming, upland paddy field, waterbody, etc., is prepared by 1:5000 scales. The land use land cover map of Nongpoh is shown in Figure 5. The data on the statistics of the LULC categories identified in Nongpoh micro-catchment are presented in Table 2.

Land Use Land Cover
Based on tonal and color variation in the satellite imagery (LISS-IV, 8-12-2015) and ground-truthing (supervised classification), the major LULC classes were identified. The spatial distribution of LULC categories, such as contour bunding, agriculture, dense forest, open forest, settlement, scrubland, terrace farming, upland paddy field, waterbody, etc., is prepared by 1:5000 scales. The land use land cover map of Nongpoh is shown in Figure 5. The data on the statistics of the LULC categories identified in Nongpoh micro-catchment are presented in Table 2.

Revised Universal Soil Loss Equation (RUSLE)
RUSLE is the method that is the most widely used approach in estimating long-term rates of inter-rill and rill erosion from field or farm size units subject to different management practices. The underlying assumption in RUSLE is that the detachment and deposition are controlled by sediment content in the runoff. The equation has following form, Equation (1).

Rainfall and Runoff Erosivity (R)
Rainfall plays a significant role in the process of soil erosion and sedimentation, and this ultimately contributes to water erosion including splash, sheet, rill, and gully erosion, triggered by water flow. The erosive force of rainfall which causes soil erosion is known as rainfall erosivity factor (R). The data of monthly rainfall (shown in Table 2) and annual rainfall are used for the calculation of R factor using Equation (2) [17,36]. The spatial variability of R is calculated in ArcGIS 10.2 domain using a spatial interpolation tool for the study area.

Revised Universal Soil Loss Equation (RUSLE)
RUSLE is the method that is the most widely used approach in estimating longterm rates of inter-rill and rill erosion from field or farm size units subject to different management practices. The underlying assumption in RUSLE is that the detachment and deposition are controlled by sediment content in the runoff. The equation has following form, Equation (1).

Rainfall and Runoff Erosivity (R)
Rainfall plays a significant role in the process of soil erosion and sedimentation, and this ultimately contributes to water erosion including splash, sheet, rill, and gully erosion, triggered by water flow. The erosive force of rainfall which causes soil erosion is known as rainfall erosivity factor (R). The data of monthly rainfall (shown in Table 2) and annual rainfall are used for the calculation of R factor using Equation (2) [17,36]. The spatial variability of R is calculated in ArcGIS 10.2 domain using a spatial interpolation tool for the study area.
where p i = monthly rainfall (mm) and P = annual rainfall (mm), and R = Rainfall and runoff erosivity (MJ mm ha −1 h −1 yr −1 ). Soil Erodibility Factor (K) Soil erodibility factor (K) is influenced by the geological characteristics, such as parent material, texture, structure, organic matter content, and porosity of the soil. As silt is the primary particle responsible for soil erosion, erodibility tends to decrease with the decrease of silt particle in the soil, irrespective of % content of sand and clay [33][34][35][36][37]. For the calculation of K factor of the study area, the soil physical properties analyses equation [38] was adopted. In this method, physical properties of soil such as soil texture, soil structure, soil permeability, and soil organic matter are used to calculate the K factor [33]. Other researchers also used this equation [33,36,37,39,40]. The soil sample was collected from the 50 different places of the watershed to calculate the physical properties of the soil. The laboratory analysis was done for all the parameters following the procedure of Deb et al. [40], and then this point value was taken in the ArcGIS domain. The K factor was calculated by using raster calculator, spatial interpolation (kriging) technique, and was calculated as in Equation (3).
where K = soil erodibility factor, M = (silt (%) + very fine sand (%) × (100-clay%)), a = organic matter (%), b = soil structure, and c = soil permeability (1-rapid, 2-moderate to rapid, 3-moderate, 4-moderate to slow, 5-slow, and 6-very slow]. Topographic Factor (LS) The LS factor represents the ratio of soil loss on a given slope length and steepness to soil loss from a slope that has a length of 22.1 m and steepness of 9%, where all other conditions are the same [35]. The LS factor is not an absolute value, instead is referenced to unity at 22.1 m slope length and 9% steepness, and hence provides significant information on the susceptibility of the terrain to soil erosion. The LS factor was computed [16,35] using the Equation (4).
where λ = slope length (m), S = angle of slope (%), m is the constant dependent on the value of the slope gradient; 0.5 if the slope angle is greater than 5%, 0.4 on slopes of 3% to 5%, 0.3 on slopes of 1 to 3%, and 0.2 on slopes less than 1%.
Crop Management Factor (C) The crop management factor (C) is used to determine the relative effectiveness of soil and crop management systems in terms of preventing soil loss and it mainly depends on vegetation canopy. The C factor was calculated using the normalized difference vegetation index (NDVI) suggested by van der Knijff et al. [41] and is given in Equation (5).
where α and β are unitless parameters that determine the shape of the curve relating to NDVI and the C factor. Van der Knijff et al. [41] found that this scaling approach gave better results than assuming a linear relationship, and the values of 2 and 1 were selected for the parameters α and β, respectively.
Conservation Practice (P) P factor reflects the effects of practices that reduce the amount and rate of water runoff and ultimately reduce the amount of erosion. The most important supporting practices are contour tillage, strip-cropping on the contour, and terrace systems. The P value ranges from 0 to 1, where 0 represents a very good manmade erosion resistance facility and 1 represents no manmade erosion resistance facility [34]. The spatial variability of input parameters of RUSLE model was calculated in ArcGIS 10.2 domain using spatial interpolation and raster calculator tool.

Modified Morgan-Morgan-Finney (MMF) Model
The Morgan-Morgan-Finney (MMF) model [26] was developed to predict soil loss from field-sized slopes. This model is flexible, simple, and also needs less data than other physically based soil erosion models [38]. The Modified MMF model separates the soil erosion process into two phases: the water phase and the sediment phase [17,28]. The water phase determines the energy of rain offered to detach soil particles from the soil mass and carried with the surface runoff. In the erosion phase, rates of soil particle detachment by rainfall and runoff were determined and also the transport capacity (TC) by runoff was calculated [28]. TC was determined by using the slope steepness, volume of overland flow, and the cover management factor (C). Improvements were made in modified MMF in the simulation of soil detachability due to precipitation, which considers leaf drainage and fall from the plant canopy.

Water Phase
Calculation of different variables within the water phase were completed using Equations (6)- (14). Equations (6)-(11) were used to calculate the kinetic energy associated with rainfall within the micro-catchment, whereas, Equations (12)- (14) were used to estimate the runoff at the micro-catchment outlet. The mean yearly rainfall was used to estimate the kinetic energy of raindrop splash and overland flow for detachment of soil particle. The total energy of rainfall was estimated by applying the intensity of rainfall and kinetic energy relationship proposed by Wischmeier and Smith (1978) [27].
where R = mean annual rainfall (mm), ER = effective rainfall (mm), and A = proportion (between 0 and 1) of the rainfall intercepted by the vegetation or crop cover.
where LD = leaf drainage (mm) and CC = proportion of canopy cover (between 0 and 1).
where DT = direct through fall (mm).
where KE = kinetic energy of the rainfall (J m −2 ). Runoff occurs when the daily precipitation exceeds the soil moisture storage capacity and indicates the moisture storage ability of the soil-crop system.
where Q = volume of overland flow (mm), R c = storage capacity of soil moisture (mm), and R o = ratio between yearly rainfall and the number of rainy days.
where R n = number of rainy days.

Sediment Phase
Soil erosion is primarily caused by the soil particle detachment from soil mass due to splashing or raindrops impact, and the soil particle being moved by the overland flow. The soil detachment caused by splash is a function of rainfall energy, soil erodibility, and the rainfall intercepted by vegetation. Rainfall energy on the surface of the ground decreases exponentially with the increase of interception rate [42,43]. The model compares the expected splash detachment rate with the overland transport efficiency, which equals the soil erosion rate to less than two values, therefore showing if the detachment is a restrictive factor. The soil transport capacity of the overland flow or runoff was established with the overland flow rate, slope, and crop management [17,28]. The calculations of different parameters in the sediment phase were completed using Equations (15)- (19).
The detachment rate by splash is mainly influenced by the soil erodibility (K) and was calculated to be the weight of soil particles detached from the soil mass per unit of rainfall energy. This is calculated by using Equation (15).
where F = annual rate of soil particle detachment by raindrop impact (kg m −2 ), K = soil erodibility factor, and KE = kinetic energy of rainfall. Transport capacity of the overland flow is influenced by the slope factor, volume of the overland flow, and the ground cover management (GC) factor. The slope factor and ground cover management factor were calculated from the satellite images in ArcGIS 10.2 domain.
where H = annual rate of soil particle detachment by runoff (kg m −2 ) and GC = ground cover (%).
where Z is constant for runoff detachment, dependent on soil cohesion, and COH = cohesion of soil in kPa.
where J = annual rate of total soil particle detachment (kg m −2 ).
where TC = annual transport capacity of overland flow (kg m -2 ), C = crop cover factor, and S = slope steepness in degree. The MMF model compared the rate of predicted splash detachment with the transport capacity by runoff or overland flow, and the lower value of both approaches was used as the rate of soil erosion [16,25]. The spatial variations over the area of all the MMF parameters were calculated with the help of ArcGIS 10.2 tools.

Calibration and Validation
The calibration and validation for both the models (RUSLE and MMF) were done for the years of 2010 to 2015 and 2016 to 2019 soil erosion data, respectively. The RUSLE calibration was done by taking the parameters like R, K, LS, C and P data from year 2010−2015 in spatial domain using ArcGIS. Values of all the parameters except C were not much changing for the calibration periods. We observed that as the land use of mainly agricultural activity changes, then the value of C factor also changed. Then, the C factor was calibrated manually by altering the α and β parameters using the trial-and-error approach with the starting values of 2 and 1 for α and β, respectively, as suggested by van der Knijff et al. [41]. Similarly, the calibration parameter for the MMF model was plant height (PH), which varies with land use type, and therefore for each land use classification the parameter value PH was optimized. The optimal value of PH was also identified in a similar approach as that of α and β parameters for RUSLE model. Following identification of the optimal values of the model parameters, the model was simulated for the validation period to evaluate the model performance relative to the observed soil loss at the microcatchment outlet.

Preparation of Best Management Plan (BMP) or Alternative Land Use Land Covers (ALULC)
The best management plan (BMP) or alternative land use land covers (ALULCs) was prepared for the prioritized areas by considering the high soil erosion rate (>40 t/ha/yr) areas within the micro-catchment, existing LULC, and slope (%). In this study, the conditions for different LULCs such as, if there is an open forest and agricultural land which have a slope of 0-8% and contribute to high soil erosion, then the suggestion was done for graded/contour bunding in the area. Similarly, for the agricultural and scrubland areas having a slope within the range of 8-33% and being observed to contribute to high soil erosion, then the recommendation was done for terracing. Likewise, for the areas having agriculture, open forest, and scrubland with a slope more than 33% and contributing to high soil erosion, afforestation was advised. Finally, for only agricultural areas having a slope of 0-8% and contributing to low soil erosion, intensive agriculture was suggested. These ALULCs were determined based on the suggestions of Srivastava et al. [44]. These suggested ALULCs were then plotted in ArcGIS environment to demonstrate a spatial map for robust planning to combat soil erosion. It is to be noted that these ALULCs were only suggested for the above-mentioned LULC classes, and the dense forest and settlement LULCs were not considered. Moreover, the output of the RUSLE approach was used in recommending the ALULC.

Results
The results obtained from field data, laboratory analyses, and analyses of models for prediction of soil loss at the micro-catchment of Nongpoh are presented.

Calibration and Validation of RUSLE and MMF Models
As mentioned earlier, the RUSLE and MMF models were calibrated and validated against six years (2010-2015) and four years (2016-2019) of observed annual soil erosion data, respectively. During the calibration of RUSLE model, the simulated and observed average annual soil erosion was found to be 59.91 and 58.22 t/ha/yr respectively. The calibration parameters α and β of the RUSLE model were found to be 2.06 and 0.97 with the standard deviation of 0.17 and 0.08, respectively. Similarly, during the validation period of the RUSLE model, the average simulated and observed annual soil erosion were noted to be 59.97 and 61.93 t/ha/yr, respectively. As the calibration parameter for the MMF model is plant height, the mean plant height results for the different land use pattern are presented in Figure 6. Clearly, the average plant height for agriculture (which is generally paddy) was 0.97 m; dense forest corresponded to the subtropical forests with an average height of 23.54 m. Mean plant height for upland paddy field was found to be 0.92 m ( Figure 6). These optimal values were identified based on literature to match the observed annual soil loss. During calibration and validation of the MMF model, the average simulated annual soil erosion was observed to be 54.05 t/ha/yr and 55.30 t/ha/yr, respectively. Figure 7 displays the scatter plot of the observed against the simulated soil loss by RUSLE and MMF models during calibration and validation. From the figure, it is clear the both models simulate the soil loss in good agreement corresponding to the observed soil loss.     The model calibration and validation statistics are provided in Table 3. The model simulation results are represented by Nash-Sutcliffe efficiency (NSE), mean absolute error (MAE), and root mean square error (RMSE). The results indicate that during both calibration and validation, the RUSLE model outperformed the MMF model. This was demonstrated by the higher values obtained for NSE, and lower values of MAE and RMSE noted for RUSLE relative to MMF model (Table 3). Nevertheless, given the values are within the ideal range of the observed values (as suggested by Moriasi et al. 2007) [44], both models are well calibrated. The following sections discuss the average RUSLE and MMF model variables during the calibration and validation process.

RUSLE Model Analysis
In this study, the RUSLE model, which has six input parameters, was used. For calculating the different input parameters, the meteorological data, soil data, satellite images (LISS-IV), and DEM were used. The values of different input parameters are discussed below.
From the analysis, it was observed that the study area comprises of high R factor values, lying between 8545.02 to 8782.42 MJ mm ha −1 h −1 yr −1 with an average of 8782.42 MJ mm ha −1 h −1 yr −1 . These high R values are due to high rainfall received in the study area which directly influences the soil erosion. The R factor was calculated from the point rainfall and then it was interpolated for the entire study watershed. The average R factor during its calibration and its spatial variability is shown in Figure 8a. Similarly, the K value is mainly influenced by the soil physical parameters and it is generally found to be high for the mountainous watersheds [33,36,37]. The present study also demonstrated the same with higher K factor values within the catchment ranging from 0.12 to 0.23 with an average value of 0.18 (Figure 8b).
The spatial distribution of LS factor is also shown in Figure 8c and it ranged between 0 and 106.24, with an average of 6.05 across the study area. Similarly, the C factor of the study area ranged from 0 to 1, whereas the dense forest area reflects much lower C values ranging from 0 to 0.12. Similarly, for agricultural land, C values ranged from 0.35 to 0.54. The highest C value was obtained in the areas of wasteland and settlement which ranged from 0.54 to 1, as shown in Figure 8d.
Soil conservation practices (P) were much lower in the study area except for some sporadic bunding and terracing in the agricultural land. The P value was 0.5 in contour bunding and upland paddy field areas and 0.6 in terrace farming areas, rest were taken as 1.
From the analysis of RUSLE model, it was found that the mean annual soil loss in the study area ranged from 0 to 1348.08 t ha −1 yr −1 with an average of 59.94 t ha −1 yr −1 , which is also falling under severe erosion risk [42,[44][45][46]. From the analysis, it was observed that the agricultural land had the highest soil erosion rate of all the portions of the study area. Figure 9a,b represents the kinetic energy dissipated by raindrop and depth of overland flow respectively. Similarly, Figure 10a,b demonstrates the total soil detachment rate and transport capacity of runoff at the study catchment respective. Additionally, the spatial distribution of soil erosion is shown in Figure 11.

MMF Model
The soil loss in MMF was determined by calculating various input parameters required in GIS environment. The input parameters of MMF models are shown in Table 4.  The spatial distribution of vegetation canopy cover and height of tree was the vital input for estimation of rainfall kinetic energy in GIS environment and it was found to have an average value of 40402.80 (J m −2 ) for the study area (Figure 9). The average overland flow (Q) for the study area was found to be 64.66 mm (Figure 8). During the calibration and validation process of the model, the annual rate of total soil particle detachment (J) ranged from 5.15-15.50 (kg/m 2 ) with an average of 5.53 (kg/m 2 ), and the annual transport capacity of overland flow (TC) lay between 0-35.72 (kg/m 2 ) with an average of 11.50 (kg/m 2 ) ( Figure 10). This model compares the rate of expected splash detachment (J) with a transport capacity of runoff or overland flow (TC), and J was noted to have lower values and it was used as the soil erosion rate.

Overall Comparison of Annual Soil Loss Estimated by RUSLE and MMF Models
The comparison of the RUSLE and MMF models showed similar trends of soil loss in the study area ( Figure 11). Both RUSLE model and MMF model results calculated for the soil erosion rate at the study area showed that the area falls under the severe erosion category. By considering the average annual soil loss from the models, it was found that the MMF model underestimated the soil erosion rate relative to the RUSLE model by 7.72%. The comparison of the RUSLE and MMF models is shown in Table 5. While considering the RUSLE model, further LULC wise soil loss was estimated ( Table 6). Both the models showed that the agricultural land contributed the maximum soil erosion rate in the study catchment. On the other hand, terracing, upland paddy fields, contour bunding, and dense forest are already protected, and they need more time to stabilize. Prioritization was done based on annual soil loss, where 40 t/ha/yr annual soil was taken as a benchmark [42,44,46,47] in the unsterilized areas of open forest, agriculture, and scrubland. This implies soil erosion values lesser and greater than this threshold are considered as a low and high erosion risk area respectively ( Figure 12). This threshold value was defined by Singh et al. 1992 [47] for the mountainous watershed of India. As mentioned earlier, the prioritization map was derived based on the soil erosion estimates of the RUSLE model. The condition for the management plan was applied in ArcGIS 10.2 environment; it is shown in Table 7. Based on the prioritization map, it can be clearly seen that the prioritization is none for open forests and settlements on the left side of the micro-catchment. An alternative LULC/suggested management plan ( Figure 13 and Table 8 Figure 13. Management plan of the study area.

Discussion
The study area is situated in mountainous regions of India which faces a major problem of soil erosion as is presented in this study. This problem occurs due to high and erratic rainfall [48,49], higher value of soil erodibility, slope length, and crop management factor (as identified in this study). Additionally, little to no conservation practices (terracing and contouring) have been employed in the study area which has led to the high soil erosion. As the watershed is mountainous with a steep slope, it inherits a high slope length factor and consequently high overland flow with high velocity. Slope length is the primary factor influencing soil erosion [50,51] and its impact is complex [52]. As the slope length increases, this leads to the high overland flow and results in the increase in soil erosion and sediment yield in the study area. High rainfall erosivity (R) was found in the study. A higher R value indicates higher kinetic energy of rainfall and surface runoff, which contributes to higher soil erosion [53]. Moreover, the high soil erodibility factor (K) directly increased the rate of soil erosion [33,37,54]. The predicted soil erosion for the study area was close to the other studies in India and other countries' hilly watersheds [54][55][56]. Soil erosion leads to loss of soil nutrients and a decrease in soil fertility that should give significant attention during the construction of the conservation management plan [57]. Based on soil loss and land use, the watershed is prioritized for conservation measures. The management plan for soil conservation and watershed development by reducing soil loss was made based on the slope, LULC, and soil loss classes. High soil erosion was found from the agricultural land and the management plans were prepared accordingly. For this, a change in the existing agricultural cultivation system should be done with intensive agriculture, contour bunding, and terracing. The conservation structures restrict soil erosion by reducing the overland flow. The contour bunding and terracing reduce the slope length and also increase the infiltration capacity of the soil.
Prior to the comparison of the RUSLE and MMF models, both models were calibrated (2010-2015) and validated (2016-2019) against observed soil loss (measured at the catchment outlet). The model parameters were adjusted manually during model calibration. Two model parameters (α and β) and one model parameter (plant height) were used in model calibration for the RUSLE and MMF models, respectively. While α and β parameters relate NDVI and the C factor, they are sensitive to land use type. Throughout the study catchment, α and β parameters varied with an average value of 2.06 and 0.97, respectively. Similarly, the plant height parameter in MMF model varied spatially due to different land use types across the catchment. As displayed in Figure 5, the calibrated values for different land use types ranged from 0 (for water body) to 23 (Table 3). These results indicate that both models are well in agreement with the observed soil loss in the study catchment with minor superior performance observed for the RUSLE model.
With regards to the BMP strategies to combat the soil erosion within the microcatchment, RUSLE method was used to identify the prioritization areas within the catchment. The prioritization plot ( Figure 12) displays the priority areas are more dominant in the right-hand side of the catchment which is predominantly dense forest, scrubland, agriculture, contour bunding, and upland paddy field types of LULC. The suggested management plan ( Figure 13) showcases that afforestation is mainly suggested for the open forest LULC, contour bunding is recommended for agricultural LULC within steep slopes, and terracing is suggested for both agricultural lands as well as scrublands (depending on slopes). Table 7  On the other hand, for scrubland LULC, 0.34 ha, 12.31 ha, and 8.34 ha land area are suggested for contour bunding, afforestation, and terracing, respectively. This indicates that a significant area is under severe erosion in case of scrubland which requires afforestation. Nevertheless, proper management practices are paramount in the study micro-catchment for reducing soil loss irrespective of the LULC.
A similar study was done using RUSLE model at West Bengal, India [58] for alike land use class of this study and found an increase in the rate of erosion with increasing areas in dense forests (134 t/ ha/yr, 55 km 2 ), degraded forests (169 t/ha/yr, 137 km 2 ), and settlement areas (30 t/ ha/yr, 105 km 2 ), which is closely in line with our results. Furthermore, several other studies at different geographic domains (mountainous and plain) of India such as in the Urmodi river watershed at Maharashtra, India with an area of 70.22 km 2 , found areas to be under very severe erosion risk (>80 t/ha/yr) [59]. In another watershed (Konar river basin), using RUSLE, the soil loss was found to be in the very high category of 5-40 t/ha/yr [60]. Synonymously, studies on similar slope and land use catchments in countries like Sri Lanka, Malaysia, and Rwanda also revealed high rate of soil erosion [55,61,62]. Moreover, in Rwanda, 91% of the study area comprised of slope more than 15 • and the average soil loss was in the range of 31-41 t/ha/yr [55], which is consistent with this study since 60% of the area is comprised of high slope (>25 • ) and falls under severe erosion classes (>40 t/ha/yr). The present study enabled us to prioritize the different segments of watershed based on soil erosion risk. This will help in proper land uses and planning conservation practices in the future extension of agricultural activities. The alternate land uses for the watershed are expected to be followed by the farming communities to achieve sustainability in agricultural productions. The results from the study will help in understanding the erosion risk of watershed vis-à-vis land uses so that the same can be extrapolated to other watersheds with similar landforms, soil, and land use. These findings from other studies are in line with the findings of this study, which indicates the wide applicability of these soil erosion models in larger catchment sizes. This study on the other hand indicates the applicability for smaller sized catchment as well for evaluating BMPs.

Conclusions
The average annual soil erosion estimated with RUSLE was found to be 59.94 t/ha/yr, whereas the MMF model estimated soil loss as 55.30 t/ha/yr. Since there was an absolute difference of only 7.74%, the estimates made by both the models can be considered at par. At this range of soil loss, the watershed can be classified as a severely eroded catchment and it needs suitable land use management plans to combat the severe soil erosion. The soil erosion was found to be highest in the agriculture dominated area (79.10 t/ha/yr), followed by scrubland (62.67 t/ha/yr). Under the alternative management plan, considering the threshold of 40 t/ha/yr, intensive agriculture is suggested only in those areas having a slope less than 8%. Within the study watershed, 68.01 ha area is prioritized for soil conservation measures to reduce the erosion potential under proposed alternative land uses. The results from the study will help in understanding the erosion risk of watershed vis-à-vis land uses so that the same can be extrapolated to other watersheds with similar landforms, soil, and land uses. These findings will aid in implementing robust BMPs in the hilly and mountainous regions in the world.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available since they are collected by the first author of this paper.