Watershed Prioritization Using Morphometric Analysis & RUSLE Model for Soil Conservation Planning, in Gilgel Abay Watershed, Ethiopia

Background: Soil erosion is the predominant global land degradation process which leads to a decline in ecosystem services and functions. The Gilgel Abay watershed is one of the watersheds in Upper Blue Nile basin characterized by rapid population growth, deforestation, overgrazing, soil erosion, waterlogging and ooding. These contribute the land and water resources of the watershed and its ecosystem are under high risk. This study was aimed to prioritize high erosion risk areas of Gilgel Abay watershed for soil and water conservation planning using the Revised Universal Soil Loss Equation model and morphometric analysis. Using the Revised Universal Soil Loss Equation model, Sub watershed Prioritization was done by computing the annual soil loss rate by considering the ve model parameters. The compound factor analysis technique was used for prioritize sub-watersheds through morphometric analysis. Results: The result of estimated annual soil loss rates in Gilgel Abay watershed ranges from 0 to 781 t ha -1 year -1 , which is higher as compared to the limits of soil loss tolerance rates of the northern highlands of Ethiopia. From RUSLE model results, Sub watershed 4 was experiencing high erosion risks while, Subwatershed2 followed by sub watershed 4 are high erosion risk areas of compound factor analysis. Conclusion: The results of prioritization through both approaches show a quite alike output that is both identies sub-watershed 4 and 2 as high erosion prone area. Therefore, both the results display a good eciency in the assessment of erosion risks and giving priority in sub-watershed scale. Hence, the two approaches can be used to identify and delineate erosion-prone areas, and prioritize the areas for effective planning of sustainable land management options.


Page 3/21
The study was done in Gilgel Abay watershed which is located in West Gojjam and Awi Administrative Zones of the Amhara National Regional State (ANRS) of Ethiopia. The watershed comprises a total of 10 Woredas and 3779.16 km 2 estimated surface area. Geographically, the watershed lies between 10.95 o and 11.80 o latitudes and 36.70 o and 37.40 o longitudes (Fig. 1). The Gilgel Abay watershed falls within the cool sub humid zone with annual temperature of 17-20 0 c. The climate is generally temperate at higher elevations and tropical at lower elevations. The mean annual rainfall of the study area varies from around 1266 mm up to 2072 mm. Rugged mountainous topography characterizes most part of the watershed especially in the southern part but there is some low land within the watershed. Elevation ranges from 1805 to 3534 m a s l. Around 80% of the watershed falls in the slope range of (0-6%), 15% of the area falls in the slope range of (7-14%) and the remaining 5% is steeper than 14% (Ashena , 2007).

Data type and sources
Both satellite and meteorological data measured from different stations of the watershed were used for the study. Two satellite images, i.e. Landsat 8 OLI and Surface Radar Topographic Mission (SRTM) digital elevation model products were freely acquired from https: https://earthexplorer.usgs.gov/ website. Landsat 8 OLI image were used for land use land cover classi cation and SRTM DEM were used for hydrological analysis (delineating watershed, subwatersheds, as well as generating slope and ow accumulation). GCPs were collected from handheld GPS for creating training areas for land use classi cation and accuracy assessment.

Satellite image processing
Before the LULC classi cation, one should preprocess the distorted and degraded images to ensure the results of adequate quality with a more correct and faithful representation of the real ground scene. The Landsat 8 OLI image used for the study were preprocessed by using ERDAS imagine version 2015 software. The preprocessing activities including, atmospheric and radiometric correction to diminish the effects of clouds and the sun elevation angle of satellite images, band composition, and image enhancement operations like contrast adjustment, edge enhancement, and color composite were employed to enhance the interpretability of image data. Finally, the area of interest (Gilgel Abay watershed) was extracted using the subset tool available in ERDAS IMAGINE software.

LULC Classi cation
A computer-aided digital image classi cation procedure was used to classify satellite images to generate thematic LULC maps. The training signatures were collected for each LULC class from the eld. A supervised classi cation method based on the Maximum Likelihood Classi er (MLC) algorithm was applied to classify the Landsat image. The watershed was classi ed into a total of seven major LULC classes. The thematic layer of the classi ed LULC images was validated using in-situ data and high resolution images like Google earth. Finally, classi cation accuracy was calculated using overall accuracy and Kappa Statistics.

RUSLE model description
There are a different models used to estimate soil erosion risk over many years worldwide (Farhan & Nawaiseh, 2015;Lal, 2001). These models can be physical-based, empirical, and conceptual. The choice to apply these models depends on the availability of input data and the type of information needed. Even though there are number of soil erosion estimation models, RUSLE model is the most widely applied empirical model which was adopted in the present study. The RUSLE model is a practical tool for predicting the long-term average annual soil loss attributed to raindrop splash and runoff. It is an empirical prediction model easily applicable in different scales. To estimate annual mean soil erosion caused by rainfall, and identify the spatial pattern of the potential soil loss risks in the watershed, all RUSLE model erosion input factors are structured in raster format and co-register to a common pixel resolution and datum (Renard KG, 1997  where A = average annual soil loss per unit of area (ton ha − 1 year − 1 ); R = the rainfall erosivity factor [MJ mm, (ha − 1 year − 1 )]; derived from daily precipitation data; K = the soil erodibility factor [ton ha − 1 h MJ − 1 ha − 1 mm − 1 )]; derived from information on soil types; LS = topographic factor, i.e., length of the slope and percent of the slope steepness (dimensionless), respectively; derived from a DEM; C = the land cover and management factor (dimensionless); derived from LULC classi cation of satellite image data; and P = the conservation support factor, which accounts for soil erosion control measures (dimensionless) derived from eld observation and literature.
Once the RUSLE model is run, classi cation of soil erosion risk area was based on the estimated rate of mean annual soil loss. Thus the erosion risk area was categorized into eight classes from very low to extremely high (Uddin et al., 2016) (Table 2). Rainfall Erosivity Factor (R) Rainfall erosivity is a climatic factor, which is estimated from the rainfall data. It is a measure of the kinetic energy (E, MJ m 2 ) based on the 30 min maximum intensity of rainfall (I30, mm h-1 ) (Wischmeier & Smith, 1978). For this study, an empirical equation developed by (Kaltenrieder, 2007) (Eq. 2) were used to estimate local erosivity values from the available total annual rainfalls. The value of R in Gilgel Abbay ranges from 558 to 736 (Fig. 2).
Where, X is mean annual rainfall in mm and R is rainfall erosivity factor.

Soil Erodibility Factor (K)
The soil erodiblity (K factor) indicates soil's susceptibility to detachment and transport by agents of erosion. K factor values were determined from the available soil data. We used the monograph method (Eq. 3) developed by (Wischmeier et al., 1971) to calculate the k factor. The k value ranges from 0.0595 to 0.234 (Fig. 3). Where, fp is the particle size parameter (unitless), Pom is the percent organic matter (unitless), Sstruc is the soil structure index (unitless), fperm is the pro lepermeability class factor (unitless), Psilt is the percent silt (unitless) and Pclay is the percent clay (unitless).

Topographic Factors (LS)
Slope  Cover Management Factor (C) The cover management factor (C) represents the effects of vegetation, management and erosion control practices on soil loss rates. The value ranges from 1.0 in completely bare land (no cover) to 0.0 in water body or completely covered land surface. Land use/ cover classi cation map and normalized difference vegetation index (NDVI) are the most commonly used methods for C-value estimation. In this study, we selected land use/cover classi cation map approach since it gives comparatively precise C-value than normalized difference vegetation index (NDVI) (Lin D, 2017) and the value ranges from 0.0196 to 0.25 ( Fig. 6).

Conservation practice (P) factor
The conservation practice factor (P factor) refers to the effects of land conservation practices in minimizing the quantity and rate of rainfall-runoff and soil erosion (Wischmeier & Smith, 1978). The conservation practice factor signi es the ratio of soil erosion from land treated with a speci c conservation measure to its equivalent soil loss from up and downslope tillage. P-value can be determined by the type of conservation measure implemented. This study was employed an alternative method using a combination of slope and land use/covers for estimation of the P-value ( Fig. 7) as proposed by (Wischmeier & Smith, 1978) in (Table 5).  Compactness coe cient Cc = 0.2821 (P/A 0.5 ) (Gravelius, 1914) 2.6 Sub watershed Prioritization using RUSLE Identifying areas that are highly vulnerable to soil loss is a critical factor for designing and implementing appropriate conservation measures. Prioritization was done at sub-watershed scales considering areas with a higher soil loss and increases in erosion risk. Using RUSLE model prioritization was done by calculating the amounts of soil loss for each sub-watershed by using the Zonal Statistics tool of Spatial Analyst extension in ArcGIS 10.5 software. Finally, prioritization of sub-watersheds for soil and water conservation activity was done based on the magnitude of soil loss in each sub-watershed and subwatersheds experiencing high mean annual soil loss rate are rst prioritized for soil conservation planning.

Sub watershed Prioritization using compound factor (CF) model
Compound factor model is a very précised method to evaluate the earth surface properties particularly for the basin scale in various works of literature and comprehensively used for sustainable water resource management and scienti c planning of sub-watersheds in data scarce areas (Altaf et al., 2013). In this approach, we used the results of linear, areal or shape parameters of morphometric results for prioritization of sub-watersheds. It has been reported that linear parameters have a direct relationship with erodibility. Thus, parameters (stream length, mean stream length, bifurcation ratio, and length of overland ow) was rated as rank 1, second highest value was rated as rank 2 and so on, and the least value was rated last in rank (Biswas. et al., 1999; Ratnam et al., 2005). By contrast, the shape parameters have an inverse relation with linear parameters, so that the lower their value, the greater the erodibility (Patel & Dholakia, 2010; Patel et al., 2012). Therefore, the lowest value of shape parameters (elongation ratio, compactness factor, circularity ratio, form factor) was rated as rank 1, the next lower value was rated as rank 2 and so on and the highest value was rated last in rank. Then Compound factor is computed by summing all the values of linear parameters as well as shape parameters and then dividing by number of parameters. The subwatershed with the lowest rank was given higher priority according to (Vandam, 2013 (Fig. 9), while the statistical details of the soil loss rates are presented in (Table 7) The extent and magnitude of soil loss rate varies from one part of the country to another, one watershed to the other depending on the farming practices, population pressure, type and susceptibility of the soils to erosion, local climate, the general terrain con gurations, and variations in agro ecological setting of the area (Monsieurs et al., 2015;Tebebu et al., 2010). In this study the magnitude and extent of soil loss varies across sub watersheds. The statistical result of the study reveals that, the mean soil loss rate is relatively higher in sub watershed4, which is found in the south east direction of the watershed. From the eld observation, this sub watershed is characterized by rugged topography and gains relatively high rain fall. The least erosion rate is recorded in sub watershed1 which is found in the Northern part of the Gilgel Abbay watershed which is at in topography.

Prioritization of sub-watersheds based on RUSLE
In a watershed management program due to time and nancial limitation, it is di cult to make rehabilitation and soil and water conservation activities at one time in all places, thus it is important to study the watersheds of the area and make ordering by their risk of erosion (Tripathi et al., 2003). Proper identi cation of areas that are highly vulnerable to soil loss is a critical factor for designing and implementing appropriate SWC measures. In this study, prioritization was done at sub-watershed scales based on the magnitude of soil erosion rates in the sub-watersheds. Among the 4 sub-watersheds, Sub-watershed 4 which is found in the southeast direction of the watershed has relatively high erosion rates followed by sub watershed3 as indicated in (Table 7) and (Fig. 11). This sub-watershed also consists of the maximum erosion rate of the watershed. As a result, the priority for soil conservation activity is given for sub-watershed 4 followed by sub-watershed 3.

Morphometric analysis
The morphometric analysis results of linear and areal parameters of the 4 sub-watersheds are computed by using different formulas (Table 6) and presented in table (8 and 9).

Linear parameters
Linear aspects of the watershed are one-dimensional characteristic of morphometric analysis which closely linked with the channel patterns of the drainage network. This study accounted linear parameters like: Stream Order, bifurcation ratio, stream length, mean stream length, drainage texture, stream frequency, in ltration number.

Stream order (u)
Stream ordering is a widely applied method for stream classi cation in a river basin. It is a measure of the position of a stream in the hierarchy of tributaries (Leopold et al. 1964) and the streams of sub watersheds in Gilgel abay watershed have been demarcated according to the Strahler's system of stream ordering. The rst order stream has no tributary, and its ow depends totally on the surface overland ow to it. Similarly, the second-order stream is formed by the junction of the two rst-order streams and thus, has a higher surface ow and the third-order streams receive ow from two second-order stream (Patel et al., 2012). In the present study, sub watershed 1 and 4 are of third-order, and sub watershed 2 and 3 are fourth level stream orders.

Stream number (Nu)
Stream number is the number of the stream segments of any particular order in a sub watershed. The values were computed using Horton's law. The orderwise stream numbers and stream length of gigel abay watershed are given in (Table 8). Streams with relatively short lengths represent the areas having ner texture with steep slopes, whereas longer lengths of stream are generally indicative of gentle slopes (Strahler, 1964).

Bifurcation ratio (Rb)
Bifurcation ratio is the ratio of the number of stream segments of given order to the number of segments of next higher order (Horton, 1945). It is a dimensionless property and shows the degree of integration prevailing between streams of various orders in a drainage basin. Bifurcation ratio is an index of relief and topographic dissection and its value varies between 2 for at or rolling catchments, and 6 for watersheds distorted remarkably by geological structure. Table 8 shows a prominent variation in the bifurcation ratio (Rb) of gigel abay sub watersheds. The value of bifurcation ratio was in between 2 and 3.6.

Drainage density (Dd)
Drainage density is a measure of the closeness of spacing of channels. It reveals how well or how poorly a watershed is drained by stream channels. The measurement of Dd is a useful numerical measure of landscape analysis and runoff potential. Drainage density is affected by the initial relief, soil in ltration, and underlying rock type. A high drainage density value indicates a relatively high density of streams, high runoff, a quick stream response, and consequently a low in ltration rate, while a low drainage density means a poorly drained basin with a slow hydrologic response (Prasad et al., 2008). The Dd value in the study ranges from 0.346 to 0.42 km/km 2 (Table 8) which shows sub watersheds are less dissected and experiencing poor drainage.

Drainage texture (Dt)
Drainage texture is the total stream number of all orders per perimeter of that area (Horton, 1945

Stream frequency (Fs)
As Horton (1932) de ned stream frequency (Fs) or channel frequency is the total number of stream segments of all order per unit area. It depends upon the lithology of the basin and the texture of the drainage network. It has a positive relation with the drainage density because of the increase in stream population is connected to that of drainage density. High stream frequency means more percolation with respect to drainage density, and thus, more groundwater potential (Sreedevi et al., 2013). The value of drainage frequency ranges from 0.0662 to 0.0934 streams/km 2 (Table 8).

In ltration number (If)
It is the product of drainage density and stream frequency. It gives us an idea about the in ltration characteristics of the watershed. In ltration number is inversely proportional to the in ltration capacity of the basin. Higher in ltration number means lower in ltration and higher runoff. A greater value of 'If' reveals impermeable surface and resistance to soil loss, in the contrary lower values point towards erosive nature of the watersheds. In the present study, sub watershed 3 and 4 has the lowest 'If' value with higher risk for erosion, and sub watershed 2 has the highest value with an indicator of resistance to erosion (Table 8).

Shape parameters
The shape of a watershed is controlled by geological structure, lithology, relief and climate, and varies from narrow elongated forms to circular or semi-circular forms. The shape mainly governs the rate at which water is supplied to the main channel. In the present study, different parameters, such as: Form factor (Ff), elongation ratio (Re), circularity ratio (Rc), and compactness coe cient (Cc) were used for characterizing watershed shape, which is an important parameter from hydrological perspective.

Elongation Ratio (Re)
Elongation ratio is the ratio of diameter of a circle having the same area as of the basin and maximum basin length (Schumm, 1956). It is a very signi cant index in the analysis of basin shape which helps to give an idea about the hydrological character of a drainage basin. The lowest value of elongation ratio indicates that high relief and steep slope, while very high values indicates the plain land with low relief and low slope. The value of the elongated ratio can be classi ed in to circular (> 0.9), oval (0.9 − 0.8) and elongated (< 0.7). A circular basin is more e cient in runoff discharge than an elongated basin (Singh & Singh, 1997). In this study, the majority of the area has high relief and steep sloped. The value of the elongated ratio I the present study ranges from 0.55 to 0.61 (Table 9) which reveals the sub watersheds are slightly elongated. These brought sub watersheds for high susceptibility to erosion resulted from generation and transportation of high amount of sediment load.

Circularity ratio (Rc)
The circulatory ratio is a quantitative expression of the shape of basin, which is expressed as the ratio of basin area to the area of circle, having the same perimeter as the basin. According to miller (1953), the values of Rc is in uenced by the stream length and frequency, geology, relief, land use and land cover and climatic condition of the basin. If the value of Rc is exactly 1 the basin is set to be a perfectly circular shape, while less than 1 indicates that the watershed is elongated in shape. Analysis results indicates the sub watersheds have a circularity ratio of in between 0.17 and 0.23 (Table 9) which indicates the geologic strata of the sub watersheds are composed of relatively homogeneous and highly permeable geologic materials. In addition, slightly elongated sub watersheds experience slow runoff discharge due to runoff from various parts of the sub watershed not reaching the outlet at the same time and highly permeable subsoil.

Form factor (Ff)
Form factor is a ratio of watershed area to the square of the length of the watershed (Horton, 1932). The value of form factor would always be greater than 0.78 for perfectly circular basin (Strahler, 1964). Short duration high peak ows are common for the sub watershed with high form factors (circular), whereas elongated sub watersheds with low form factor indicate ow for longer duration (Chopra et al., 2005). In the present study, the Ff was in between 0.24 and 0.29 indicating no sudden peak discharge at the outlet in the sub watersheds.

Prioritization of sub-watersheds using morphometric analysis
Morphometric analysis is a signi cant tool for the prioritization of sub-watersheds even without considering the soil parameters (Biswas. et al., 1999). Though By prioritization of watersheds, one can conclude which watershed can lead higher amount of discharge due to excessive amounts of rainfall. In this study, Sub watersheds are prioritized for soil and water conservation planning based on morphometric analysis of the linear and shape parameters. Computing compounding factor is the technique used to prioritize sub-watersheds as indicated in (Table 10). The compound factor was computed by summing all the values of linear parameters as well as shape parameters and then dividing by the number of parameters. The sub-watershed with the lowest rank was given higher priority according to (Vandam, 2013). The Results of the analysis showed that sub-watershed2 (cp = 1.8) followed by sub watershed4 (cp = 2.3) were high erosion risk areas & prior soil control measures are required in these sub-watersheds to preserve the area from further erosion. The compounding factor value for Sub watershed1 is the least among the four sub-watersheds and gets the least rank. This study was aimed to identify the most erosion vulnerable sub-watersheds for soil erosion in Gilgel Abay watershed. The empirical erosion assessment model and morphometric analysis of the linear and shape parameters were used. The two techniques use different datasets and follow different approaches to prioritize high erosion risk areas or sub-watersheds. The RUSLE potential soil loss estimation model gives a good implication of soil erosion intensity and variability in Gilgel Abay watershed. Among the four sub-watersheds in Gilgel Abay, sub watershed4 was experiencing relatively high mean soil loss and the rst priority was given to this watershed due to a relatively high magnitude of soil erosion followed by sub watershed 3. Sub watershed1 was the least prioritized sub-watershed because of relatively less intensity of soil erosion.
The morphometric analysis of the watershed gives the potential high erosion risk areas in the watershed by considering the linear and shape aspects instead of quantifying the amount of erosion loss in each cell. The compounding factor value of sub watershed 2 is low as compared to the other sub-watersheds and it was rst prioritized for soil and conservation planning and followed by sub watershed 4.
The present study demonstrated the effectiveness of RUSLE and compound factor-based prioritization methods with the assist of GIS techniques. Both approaches can be used to identify and delineate erosion-prone areas, and prioritize the areas for effective planning of sustainable land management based on erosion severity levels in the watersheds.

Declarations
Ethics approval and consent to participate An attempt was made to conduct the research in an ethical manner. Hence, as the researcher, the Authors take full responsibility for all the contents and any mistakes here in the document.

Consent for publication
We all have agreed to submit our nal manuscript for Environmental Systems Research journal and approved the submission.

Availability of data and materials
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Competing interests
The authors declare that they have no competing interests.

Funding
No funding was received.

Authors' contributions
The corresponding author has led the overall research process starting from the design to data analysis and interpretation of results. Both the co-authors have been involved by providing critical comments and suggestions during the research process and manuscript write up. All the authors have read and approved the nal manuscript. slope length (L factor) in Gilgel Abay watershed. Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors. slope steepness (s factor) for Gilgel Abay watershed Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors. land cover management factor(C factor) Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors. Annual soil loss rates in Gilgel Abay watershed Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.

Figure 10
Soil erosion rate in the sub-watersheds. Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.