Assessment of Soil Erosion Using the RUSLE Model for the Epworth District of the Harare Metropolitan Province, Zimbabwe

Urban development without adequate soil erosion control measures is becoming a major environmental concern in developing urban areas across Africa. These environmental disturbances encompass rampart Land Use and Land Cover changes (LULC) due to a high population growth rate and increased economic activities. To understand the influence of accelerated LULC changes and urban expansion as major drivers in landscape degradation in the Epworth district of the Harare Metropolitan Province, the RUSLE model was employed. This considers land use, soil, climate and topography as input parameters in the assessment of the extent and impact of these drivers on soil erosion. The Revised Universal Soil Loss Equation (RUSLE) was used to predict the potential erosion between 1984 and 2018 and soil erosion risk for the years 2000 and 2018. The mean rate of the predicted potential soil erosion was 13.2 t ha−1 yr−1 (1984–2018); areas especially vulnerable to erosion were predicted for foot slope areas with direct tributaries to the major streams and steep sloping zones. The average soil erosion risk was estimated at 1.31 t ha−1 yr−1 for the year 2000 and 1.12 t ha−1 yr−1 for 2018. While the overall potential soil loss decreased between 2000 and 2018, the potential soil loss was observed to increase tremendously in residential areas, which doubled in extent between 2000 and 2018. The findings reveal that about 40% of the Epworth district was threatened by unsustainable soil loss resulting from increased soil erosion risk within the built-up areas.


Introduction
Urbanization is a continuous process that has boldly accelerated with population increase, expansion and spread of built-up structures in a designated urban area [1]. Urbanization in Africa has been growing at alarming rates with an anticipated annual growth rate of approximately 3.9% [2]. Population growth and increasing economic activities have been linked to aggravate Land Use and Land Cover (LULC) changes [3][4][5]. As such, urban development inevitably involves construction and sealing activities that alter natural landscapes [6] resulting in an increase in impervious surfaces, which replace natural vegetation and reduce the capacity for water infiltration.
This in-turn results in surface runoff that substantially threatens soil loss in vulnerable landscapes through erosion processes [7]. While multiple studies on urbanization processes focus on the social and planning aspect and the assessment of the effects of urbanization processes on the (quasi-) natural environment, in the presented study we assess the effects of urbanization on surface processes by water in the strongly urbanizing area of Harare in tropical Africa; as case study we selected the 35 km 2 largest growing informal urban district of the Harare Metropolitan Province, which showed urbanization Settlement activities such as urban agriculture, construction, illegal sand mining, brick molding and effluent discharge from industries pose water quality threats to the Lake Chivero drainage basin [8,58]. Epworth district is largely dominated by high density residential areas. These are characterized by overcrowdings and concentrated housing residential structures due to densification; only little vegetation occurs across the settlements [8].
The general topography of Epworth district is characterized by undulating and slightly rolling terrain in the upland areas, being part of the southern Africa Highveld. Elevations range from 1455 m to 1556 m a.s.l. For Epworth district, clayey Fersiallistic soils (moderately leached soils of the kaolinitic order) occur in contact zones and Paraferrallitic soils (comprises of highly leached soils) are widely distributed across the entire district [59,60]. The sub-tropical climate of Harare Metropolitan Province is dominated by four distinct seasons: the cool-dry season (mid-May to August), hot-dry season (September to mid-November), rain-wet season (mid-November to mid-March) and the post rainy season (mid-March to mid-May) [11]. Harare Metropolitan Province receives 470-1350 mm of precipitation annually; rainfall predominantly occurs during the four months of rainy season. Average temperatures range from 7 °C to 20 °C during dry periods and from 13 °C to 28 °C in hotdry periods [11].

Parameter estimation for Soil erosion risk assessment using RUSLE
The RUSLE model is simple and the mostly used computerized version of the Universal Soil Loss Equation (USLE), a statistical model developed to estimate the annual soil loss per unit area based on erosion factors [43,61]. The RUSLE model has been widely implemented for the prediction of average annual soil losses caused by sheet and rill erosion and to display the spatial distribution of potential erosion risk [23,55,[61][62][63][64][65]. The application of the RUSLE model for soil erosion risk

Parameter Estimation for Soil Erosion Risk Assessment Using RUSLE
The RUSLE model is simple and the mostly used computerized version of the Universal Soil Loss Equation (USLE), a statistical model developed to estimate the annual soil loss per unit area based on erosion factors [43,61]. The RUSLE model has been widely implemented for the prediction of average annual soil losses caused by sheet and rill erosion and to display the spatial distribution of potential erosion risk [23,55,[61][62][63][64][65]. The application of the RUSLE model for soil erosion risk considers the rainfall erosivity factor (R), soil erodibility factor (K), slope length and steepness factor (LS), land cover and management factor (C) and the support practice factor (P) [43]. In the current study, the RUSLE model was adapted for mapping potential erosion using C and P factors as identity elements (C and P = 1) and for the spatial distribution of soil erosion risk.
According to [43], the Revised Universal Soil Loss Equation (RUSLE) states that: where: A is the annual average of soil erosion rate factor (t ha −1 yr −1 ); R is the rainfall erosivity factor (MJ mm ha −1 h −1 yr −1 ); K is the soil erodibility factor (t hMJ −1 mm −1 ); LS is the dimensionless slope length and steepness factor; C is the dimensionless crop management factor (ranging between 0 and 1) and P is the dimensionless conservation support practice factor (ranging between 0 and 1). The calculation of the potential erosion is based on the same formula while adjusting factors C and P to one [3,65]. Spatial analyses in the study were performed using ArcGIS 10.2 software to assess the dynamics of soil erosion over time in the heterogenous urban landscapes for the years 2000, 2018 and the overall, 1984-2018 long-term rainfall data were used for rainfall erosivity factor computation ( Figure 2). The LS factor RUSLE geospatial input factor was computed using a hydrology module (LS-factor field base) in SAGA GIS [66]. The acquired geospatial input parameters for the RUSLE model (Table 1) were used to produce thematic maps for the estimation of potential erosion and soil erosion risk generated within every cell grid [67][68][69]. For data harmonization, we resampled all data sources to determine 30 m × 30 m grid cell size using nearest neighborhood technique so as to retain original pixel value before carrying out grid cell calculations. This was performed to enhance data compatibility from varying data sources used for modelling [70].
According to [43], the Revised Universal Soil Loss Equation (RUSLE) states that: where: A is the annual average of soil erosion rate factor (t ha −1 yr −1 ); R is the rainfall erosivity factor (MJ mm ha −1 h −1 yr −1 ); K is the soil erodibility factor (t hMJ −1 mm −1 ); LS is the dimensionless slope length and steepness factor; C is the dimensionless crop management factor (ranging between 0 and 1) and P is the dimensionless conservation support practice factor (ranging between 0 and 1). The calculation of the potential erosion is based on the same formula while adjusting factors C and P to one [3,65]. Spatial analyses in the study were performed using ArcGIS 10.2 software to assess the dynamics of soil erosion over time in the heterogenous urban landscapes for the years 2000, 2018 and the overall, 1984-2018 long-term rainfall data were used for rainfall erosivity factor computation ( Figure  2). The LS factor RUSLE geospatial input factor was computed using a hydrology module (LS-factor field base) in SAGA GIS [66]. The acquired geospatial input parameters for the RUSLE model (Table  1) were used to produce thematic maps for the estimation of potential erosion and soil erosion risk generated within every cell grid [67][68][69]. For data harmonization, we resampled all data sources to determine 30 m × 30 m grid cell size using nearest neighborhood technique so as to retain original pixel value before carrying out grid cell calculations. This was performed to enhance data compatibility from varying data sources used for modelling [70].   Rainfall erosivity factor (R) describes the ability of rainfall to trigger soil erosion [69,73,74]. The RUSLE model makes use of this erosivity factor (R [MJ mm ha −1 h −1 yr −1 ]) in integrating the effects of raindrop impact, rainfall duration and resulting runoff rates, which are coupled with the amount and the energy within each recorded rainfall pattern [43,69]. Rainfall erosivity was calculated using mean annual rainfall data collected from three gauging stations in Harare Metropolitan Province (Harare Belvedere, Harare airport and Harare Kutsaga; Table 2) following Equation (2) [40,64]: where R = Rainfall erosivity factor, M = Mean annual rainfall. To analyze possibly changing rainfall erosivity since 1984, annual rainfall averages were calculated for the time steps: 1984-2000, 2000-2018 and overarching, 1984-2018. The mean annual precipitation data were interpolated over entire district by applying the inverse distance weighting (IDW) interpolation technique and converted to rainfall erosivity by applying Equation (2).

Soil Erodibility Factor (K)
The responsive effect of a particular soil in a given location to the erosive power of rainfall and runoff impacts is referred to as the soil erodibility factor (K) [16,75]. Soil erodibility is regarded as a function of the soil texture, organic matter content, soil structure and the degree of permeability [41,76]. Soils being highly susceptible to erosion have soil erodibility values close to 1, whereas corresponding values close to 0 indicate a resistive nature of the soil [51,77]. In the current study, information on soil structure and profile permeability was not available. Therefore, the K factor was estimated using the ISRIC (International Soil Reference Information Centre)-World Soil Information data [78], following the equation by [79].

Topographic Factor (LS)
The LS factor summarizes the effects of topography on soil erosion and combines the influence of slope length and slope angle on soil loss; while the S-factor measures the effect of slope steepness, the L-factor defines the impact of slope length. The slope length L is defined as the distance between the upslope starting point of a slope segment to the downslope point where deposition begins [33,66]. Increasing slope length and steepness per unit area results in increased runoff and flow velocity and consequently in increased soil loss exposure [42,75,80]. The Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) with a spatial resolution of 30 m was used for the LS factor computation applying SAGA software 2.3/Hydrology module [66,80]. For data pre-processing, fill sink algorithm was employed using the spill elevation method for filling sinks on the DEM [81]. Multiple flow direction tool (MFD) incorporated in the Hydrology module was applied to the DEM to assign flow directions and flow accumulation [82,83]. The Land cover and management factor C's pivotal role is to capture differences in soil loss in vegetated areas by dissipating raindrop impact on the soil surface compared to bare areas [43,84,85]. The C factor decreases from 1 to 0 depending on vegetation cover and cropping management systems implemented to mitigate soil erosion [41,80]. In order to determine the C factor, land use maps generated from the supervised classification of satellite images were adapted for use [9]. The C factor values for each land use are a result of weighted average of the soil loss ratio deduced from a reference plot (bare) with a C factor of 1 [43,66]. Henceforth, the C factor values of each land use type were evaluated from literature [3,[86][87][88]; weighting of the data was performed following field observations and the biophysical characteristics per sampling plot (canopy cover, prior land uses, and surface cover) [43,89]. The weighted C factor evaluation considered plant growth, height and the extent of canopy cover in-situ [89,90], in relation to bare and sealed area. Further, previous urban farming practices and residues of the plant material influence were majorly factored in during surface cover subfactor assessment per sampling plot. For croplands the C factor was estimated by averaging the values of predominant crops within the study area's croplands [91]. The availability of remote sensing data, especially the supervised classification maps, appropriately aided the evaluation of spatial variability of the C factor [89]. Overall, the weighted C factors were estimated as a result of multiplying the scaled percentages of the evaluated C subfactors (Table 3), reviewed RUSLE C factors according to the literature [86][87][88] and the ratio of sealed area proportion [9] in relation to the reference plot (bare).

Support Practice Factor (P)
The Support practice factor P expresses the effects of surface management practices that are applied to reduce soil loss through erosion processes [30,43]. These practices include among others terracing, strip cropping and contour ploughing [43]. The P factor value ranges between 0 and 1, where 0 shows the highest effectiveness of the conservation practice and 1 indicates that there are no support practices or measures implemented [43,92]. In this study a P factor of 1 has been all-over applied due to the area wide absence of support or management practices.

Mapping and Surveying Soil Erosion Dynamics
The complexity of the setup of urban built-up areas and the distribution of the different land use require soil erosion field survey mapping to receive reference data [23]. In this regard, a simple snap-shot sampling procedure was implemented during the field survey in December 2019 to estimate the spatial extent of eroded areas on plots 40 m × 80 m in size [23]. In total, 49 sites were randomly surveyed, covering varying land use. Mapping of spatial soil erosion assessment was geocoded using a hand-held GPS (Garmin 60Cx); parallel on-site soil erosion features were measured, considering the erosion features (inter-rills and rills) individually calculating their area and volumes [55,93]. Rills were Sustainability 2020, 12, 8531 8 of 24 defined as linear erosion channels not more than 0.5 m deep and with a cross-sectional area <929 cm 2 to make them distinct from ephemeral gullies and deep incised gully features [94,95].
Overall total percentages of disturbed surface area data were registered for the sample plots in the Epworth district to estimate the spatial extent of perturbed regions of the district. The spatial extent of eroded areas was estimated as percentage per plot [23] and classified into five severity classes. Mapping results are displayed in a positional diagram map.
The complexity of the setup of urban built-up areas and the distribution of the different land use require soil erosion field survey mapping to receive reference data [23]. In this regard, a simple snapshot sampling procedure was implemented during the field survey in December 2019 to estimate the spatial extent of eroded areas on plots 40 m x 80 m in size [23]. In total, 49 sites were randomly surveyed, covering varying land use. Mapping of spatial soil erosion assessment was geocoded using a hand-held GPS (Garmin 60Cx); parallel on-site soil erosion features were measured, considering the erosion features (inter-rills and rills) individually calculating their area and volumes [55,93]. Rills were defined as linear erosion channels not more than 0.5 m deep and with a cross-sectional area < 929 cm 2 to make them distinct from ephemeral gullies and deep incised gully features [94,95].
Overall total percentages of disturbed surface area data were registered for the sample plots in the Epworth district to estimate the spatial extent of perturbed regions of the district. The spatial extent of eroded areas was estimated as percentage per plot [23] and classified into five severity classes. Mapping results are displayed in a positional diagram map.

Potential Erosion Risk Analysis
The potential erosion risk map was derived from the application of the "natural" RUSLE factors for soil characters, rainfall erosivity and topography (K, R and LS) for the Epworth district. Resulting data are classified into five potential erosion risk classes showing how erosion varies in the Epworth district ( Figure S2). The potential erosion risk map indicates the vulnerability of the landscape independent from vegetation cover and crop management. The findings reveal that very high to extreme erosion risk areas occur in areas of steep slopes ( Figure S2); the only spatially slight variations of R and K factors cause a strong control of potential erosion risk by topographic factors' LS. Due to the general orientation of the drainage network to the south and downstream with increasing inclination of the slopes along the valley flanks, the areas with high erosion risk expand from the north to south. Very high to extreme erosion risk areas are also observed in the eastern parts of the Epworth district resulting from locally occurring steep slopes towards the middle course of the Jacha river and tributaries. In contrast, low to moderate potential erosion risk areas occur in the plateau areas of the Highveld. Extreme potential erosion risk zones occur in the immediate vicinity of streams and at steep slopes. The area wide average potential erosion risk rate in the Epworth district was 13.2 t ha −1 yr −1 since 1984, referring to the precipitation data 1984-2018 ( Figure S2, Table 5).

Potential Erosion Risk Analysis
The potential erosion risk map was derived from the application of the "natural" RUSLE factors for soil characters, rainfall erosivity and topography (K, R and LS) for the Epworth district. Resulting data are classified into five potential erosion risk classes showing how erosion varies in the Epworth district ( Figure S2). The potential erosion risk map indicates the vulnerability of the landscape independent from vegetation cover and crop management. The findings reveal that very high to extreme erosion risk areas occur in areas of steep slopes ( Figure S2); the only spatially slight variations of R and K factors cause a strong control of potential erosion risk by topographic factors' LS. Due to the general orientation of the drainage network to the south and downstream with increasing inclination of the slopes along the valley flanks, the areas with high erosion risk expand from the north to south. Very high to extreme erosion risk areas are also observed in the eastern parts of the

Soil Erosion Risk
The estimated soil erosion risk maps were generated for 2000 and 2018. The estimated soil erosion risk averaged 1.31 t ha −1 yr −1 in 2000 and 1.12 t ha −1 yr −1 in 2018; with highest total soil loss rate for the Epworth district amounting to 92.79 t ha −1 yr −1 in 2000. The spatial patterns of the estimated soil erosion risk indicate areas with high soil erosion loss predominantly along the river courses ( Figure 5).
Correspondingly, areas of high soil erosion risk can also be found in the southwestern and southeastern parts of the Epworth district with highest soil erosion risk close to the river courses.

Soil Erosion Risk
The estimated soil erosion risk maps were generated for 2000 and 2018. The estimated soil erosion risk averaged 1.31 t ha -1 yr -1 in 2000 and 1.12 t ha -1 yr -1 in 2018; with highest total soil loss rate for the Epworth district amounting to 92.79 t ha -1 yr -1 in 2000. The spatial patterns of the estimated soil erosion risk indicate areas with high soil erosion loss predominantly along the river courses ( Figure 5). Correspondingly, areas of high soil erosion risk can also be found in the southwestern and southeastern parts of the Epworth district with highest soil erosion risk close to the river courses.  The soil erosion damage map shows plot-wise surface area damages in the Epworth district calculated from the soil erosion feature dimensions, expressed in percentages of the total surface area for each plot ( Figure 5). According to this soil erosion damage map, the southeastern and southwestern plots of the Epworth district experience higher soil erosion compared to the other areas investigated. The spatial extent of the soil eroded area ranges from 0% to 1.4% for the surveyed plots in the Epworth district. The magnitude of soil erosion observed during field mapping in 2019 indicated that slope, high proportion of sealed and impervious surfaces attributed to increased soil erosion damages in the Epworth district ( Figure 6). calculated from the soil erosion feature dimensions, expressed in percentages of the total surface area for each plot ( Figure 5). According to this soil erosion damage map, the southeastern and southwestern plots of the Epworth district experience higher soil erosion compared to the other areas investigated. The spatial extent of the soil eroded area ranges from 0% to 1.4% for the surveyed plots in the Epworth district. The magnitude of soil erosion observed during field mapping in 2019 indicated that slope, high proportion of sealed and impervious surfaces attributed to increased soil erosion damages in the Epworth district ( Figure 6). The estimated soil erosion risk for the year 2000 highlights that 56.3% of the Epworth district was exposed to low soil erosion risk and 25.9% to moderate soil erosion risk, while 15% of the Epworth district was exposed to high, and 2.8% to very high and extreme soil erosion ( Figure 6, Table  6). For 2018 modelling of soil erosion risk displays a slight decline of risk of exposure with 59.5% of the area being exposed to low soil erosion risk and 29.3% to moderate soil erosion risk; the spatial extent of areas exposed to high soil erosion risk declined to 10% and areas exposed to very high to extreme soil erosion risk covered 1.2% of the Epworth district (Table 6).  The estimated soil erosion risk for the year 2000 highlights that 56.3% of the Epworth district was exposed to low soil erosion risk and 25.9% to moderate soil erosion risk, while 15% of the Epworth district was exposed to high, and 2.8% to very high and extreme soil erosion ( Figure 6, Table 6). For 2018 modelling of soil erosion risk displays a slight decline of risk of exposure with 59.5% of the area being exposed to low soil erosion risk and 29.3% to moderate soil erosion risk; the spatial extent of areas exposed to high soil erosion risk declined to 10% and areas exposed to very high to extreme soil erosion risk covered 1.2% of the Epworth district (Table 6).

Magnitude of Soil Erosion in Epworth District
A spatial extent of about 765 m 2 was eroded with an average area of 31 m 2 affected by soil erosion as calculated from the 49 randomly selected sample plots in Epworth district during the field survey in 2019 ( Table 7). The soil erosion damage measured approximately 0.5% of the total area mapped (15.7 ha). The occurrence of soil erosion features varied in the surveyed plots corresponding to vegetation cover, slope characteristics and human activities. Model validation was done using the empirical RUSLE model data in comparison with on-site field measurements. The results showed good RUSLE model performance as there was satisfactory moderate positive correlation between field measurements and model results for sample areas (r = 0.76 and R 2 = 0.581, p <0.05) (Figure 7). This provides confidence in the application of the model for sustainable land use planning and decision-making processes.

Magnitude of soil erosion in Epworth district
A spatial extent of about 765 m 2 was eroded with an average area of 31 m 2 affected by soil erosion as calculated from the 49 randomly selected sample plots in Epworth district during the field survey in 2019 ( Table 7). The soil erosion damage measured approximately 0.5% of the total area mapped (15.7 ha). The occurrence of soil erosion features varied in the surveyed plots corresponding to vegetation cover, slope characteristics and human activities. Model validation was done using the empirical RUSLE model data in comparison with on-site field measurements. The results showed good RUSLE model performance as there was satisfactory moderate positive correlation between field measurements and model results for sample areas (r = 0.76 and R 2 = 0.581, p <0.05) (Figure 7). This provides confidence in the application of the model for sustainable land use planning and decision-making processes.

Land Use and Soil Loss Analysis
The results show that about 50,408 tons of soil were estimated to be lost under 2000 LULC conditions, while an estimated total soil loss of 42,934 tons was calculated for 2018 (Table 8). For the industrial areas of Epworth district, approximately 40 tons of soil loss were estimated for 2000, while an increase of up to 47 tons of soil loss was estimated for the same land use type for 2018. For 2000, for the land use type "less concentrated residential area" (15.5% of the Epworth district in 2000) 6218 tons of soil loss were estimated while, for the land use type "concentrated residential areas" (38% of the Epworth district in 2000) about 14,018 tons total soil loss were estimated. An increase in soil erosion risk for less concentrated and concentrated residential areas were estimated to amount 12,203

Land Use and Soil Loss Analysis
The results show that about 50,408 tons of soil were estimated to be lost under 2000 LULC conditions, while an estimated total soil loss of 42,934 tons was calculated for 2018 (Table 8). For the industrial areas of Epworth district, approximately 40 tons of soil loss were estimated for 2000, while an increase of up to 47 tons of soil loss was estimated for the same land use type for 2018. For 2000, for the land use type "less concentrated residential area" (15.5% of the Epworth district in 2000) 6218 tons of soil loss were estimated while, for the land use type "concentrated residential areas" (38% of the Epworth district in 2000) about 14,018 tons total soil loss were estimated. An increase in soil erosion risk for less concentrated and concentrated residential areas were estimated to amount 12,203 tons for the "less concentrated residential areas "(31.5% of the Epworth district in 2018) and 19,858 tons for the "concentrated residential areas" (52.6% of the Epworth district in 2018). A decline in the estimated soil loss was observed for land use types either of agricultural use or covered by green spaces (undifferentiated) between 2000 and 2018, decreasing proportional to the reduction of the areas of these land use types (Table 8).

Discussion
For the tropics, studies reported average soil loss rates of 5 t ha −1 yr −1 [96,97], while [98] highlights that a soil loss limit of 11 t ha −1 yr −1 may be accepted as reasonable mean annual loss due to soil erosion. However, [99] argues that for sensitive and fragile land areas average soil loss tolerance rate of 2 t ha −1 yr −1 could be recommended. In contrast, considering the slow rate of soil formation and spatio-temporal effects of soil loss on water quality and productivity, tolerance limits for soil erosion loss are set for the tropics at 1 t ha −1 yr −1 [3,65,100]. The occurrence of soil loss exceeding 1 t ha −1 yr −1 was considered as the critical rate for the Epworth district due to the low rate of soil formation as typical for the tropics [3]. This agrees with the recommendation of [98] that the rate of soil loss through soil erosion should be relatively balanced with the soil formation rate to minimize excessive environmental damage. As such, average soil loss rates exceeding the suggested 1 t ha −1 yr −1 are classified as unsustainable to continue supporting land use [101,102].
The empirical RUSLE model implemented in this study was applied to predict potential erosion risk and soil erosion risk for the Epworth district of the Harare Metropolitan Province ( Figure S2 and Figure 5). The calculation of the potential erosion risk is based on the assumption that there is no land use and no land management as well as no support practice; potential erosion is understood as the erosion processes only controlled by physical factors. Consequently, potential erosion risk depicts areas vulnerable to erosion even without considering land use [65]. For the Epworth district the potential erosion risk was averaged at 13.2 t ha −1 yr −1 between 1984 and 2018, significantly exceeding the soil loss tolerance limit of 1 t ha −1 yr −1 for the tropics [3,100].
Estimated soil losses due to soil erosion risk averages amounted to 1.31 t ha −1 yr −1 and 1.12 t ha −1 yr −1 in the years 2000 and 2018, respectively. Correspondingly, the revealed soil loss due to soil erosion risk for Epworth district was slightly above the recommended tolerable limits of 1 t ha −1 yr −1 [3,65,100]. Considering the proposed range of tolerable maximum annual soil loss in tropical regions, it can be deduced that the mean estimated soil loss during all study periods slightly causes irreversible soil erosion. The occurrence of high potential erosion risk compared to low soil erosion risk is due to the assumption that for the calculation of potential erosion risk, the factors land cover and management (C) and support practice (P) are not considered, and consequently, these factors are mathematically handled as identity elements in order to assess the impact of RUSLE "natural factors" on the study area. In the application of the RUSLE this corresponds to dealing with C and P factors as bare ground [3]. Such conditions earmark the impact of soil erosion on cleared land area and also reveal the significance of vegetation cover in dissipating raindrop energy impact on the bare ground. Nevertheless, such scenarios merely occur in urban areas as a result of built-up densification and spread of impervious surfaces unless croplands and disturbed green spaces exist at larger scale.
The decrease of soil loss due soil erosion risk for the Epworth district calculated for the period 2000-2018 probably resulted from the expansion of built-up areas at the expense of green spaces and cropland areas and thus with sealing the underground impeding surface exposure to erosion [5,103]. However, soil disturbance has been triggered by human activities in the district due to high population pressure and demand for shelter [104,105]. This contributes to increased soil loss within the urban built-up areas compared to the previous years ( Table 8). The replacement of green spaces with impervious surfaces causes a reduction in surface area for water infiltration [106,107], causing increased overland flow either by sheet flow or by concentrated surface runoff. This substantially threatens soil loss through erosion where the overland flow reaches areas where soils are exposed to the surface making them highly vulnerable landscapes with mostly bare grounds [5,7].
Construction activities create artificial slopes and reduce vegetation canopy cover exposing soil surfaces to raindrop impact thereby exacerbating rates of soil detachment and transportation during rainstorm events [30]. Beyond, construction activities repeated earth movements affect the stability of the soil structure and increase soil erodibility [74,108] and soil compaction resulting in the reduction of infiltration capacity and increase surface runoff generation [107,108]. Therefore, continued LULC changes dominated by the expansion of built-up areas largely perpetuate soil loss within the built-up areas [103]. To ensure that soil erosion risk thresholds remain sustainable, land suitability analysis should be considered to enhance land use allocation and proper land servicing require to be implemented by responsible authorities in order to meet competing demands for land. The findings of this study highlight that zones with high potential erosion risk and soil erosion risk correspond to zones of strong relief and high slope lengths (high LS factors).
The weighted C factor values derived based on the biophysical properties observed and measured in the field [42,66]. The approach integrated LULC maps derived from satellite images. However, there are limitations on the use of LULC data resulting from misclassifications, heterogeneity and spatial distribution of vegetation densities across the entire district. The resulting approach anticipates that the same LULC class poses the same land cover and management factor C value [47,89,109,110]. However, uncertainties and limitations influencing the results could have also emanated from omitting other biophysical characteristics during estimations of the C factor including surface roughness and below-ground biomass [42,111,112].
An analysis of the relationship between LULC and estimated soil erosion risk was performed by overlaying LULC and soil erosion risk maps for the time slices 2000 and 2018 (Table 8). This relationship has been observed as a useful tool monitoring patterns of LULC change against soil erosion risk [3,100]; the analysis reveal spatial patterns and changes for each LULC class majorly as influenced by human activities in relation to soil erosion risk. Soil erosion risk was extremely high in rain fed croplands in the year 2000, highlighting their vulnerability to water induced erosion. This might have been propelled by leaving little to no crop residues as soil cover on the fields, exposing bare soil to rainfall at the onset of rainy season [22,24,113]. In addition, pre-season land preparation exposes fine tilled land to raindrop impact, exacerbating soil loss due to reduced surface and ground cover in dissipating and scattering raindrop energy [114], while concurrently ploughing increases surface roughness and pore volume, both fostering infiltration [115,116].
The estimated soil loss in areas covered by green spaces initially was predicted to be high considering the vegetation's ability to intercept raindrop impacts [66]. However, it has to be clearly differentiated what character the vegetation cover has and whether the area is disturbed or undisturbed by human impact. Especially areas with sparse disturbed green spaces are exposed to erosion processes by the first rains coming after a dry period due to the hydrophobic character of dried out soils while parallel the soil stabilization by roots is insufficient [22,23]; footpaths spreading across areas increase soil sealing along the paths and thus foster generation of concentrated runoff [24,117]. Other local human activities including harvesting of firewood for domestic use, burning of bricks during their processing and fencing of gardens have negative impacts on vegetation cover. The more intense land use gets, especially transferring fallow land covered by sparse green spaces into built-up areas, the higher the soil erosion risk in this area, which is highly affected by the increasing areas characterized by sealed and impervious surfaces [118]. Soil erosion risk investigations on the Epworth district show that, despite the decreasing of the overall estimated soil loss from 50,408 tons to 42,934 tons between 2000 and 2018, the expansion of urban built-up areas at the expense of croplands and green spaces has locally increased soil loss risk.
The expansion and spread of residential areas have been linked with increases in soil erosion risk within areas of intense human activities ( Table 8). The relationship between LULC and soil erosion risk analysis reveals that in concentrated residential areas in 2000 about 14,018 tons of estimated soil loss occurred, while about 19,858 tons of estimated soil loss occurred in 2018. The estimated increase of soil loss in concentrated residential areas corresponds to the massive growth of built-up areas and the coinciding increase of impervious surfaces [119]. The observed changes in built-up areas are significantly attributed to population growth in the Epworth district. The population was estimated to have increased from approximately 114,067 in 1992 to about 167,462 in 2012 [57]. The resulting reduced infiltration rates contributing to high rates of surface runoff [120], result in soil damage down the slope and parallel to roads [117]. The increase in estimated soil loss in the residential areas of the Epworth district are comparable with those in Kinshasa/DRC, where it is observed that highest soil erosion risk rates spatially correlate with steep slopes along river flanks and increasing density of informal settlements [38]. Estimated soil loss in less concentrated residential areas increased from 6218 tons in 2000 to approximately 12,203 tons in 2018. This corresponds to a doubling in size of the land use class "less concentrated residential areas" in the Epworth district from 15.5% in 2000 to 31.5% in 2018 (Table 8). The increase in soil loss for less concentrated residential areas may have resulted from the registered decrease in green spaces (Table 8). This is possibly attributed to high population pressure and demand for shelter, hence propelling landowners to informally construct shacks on their backyards to curtail housing demand. These activities inevitably occur at the expense of green spaces and ground cover resulting from clearing of land for construction [6,121].
Steep areas have high potential erosion risk as well as estimated soil erosion risk, especially on the flanks of the streams [122,123]. This corresponds to the high loading of slope steepness in the RUSLE model [124,125]. The combination of slope inclination with slope length contributes to the cumulative effect of increasing surface flow with increasing the drainage basin resulting in increased soil erosion risk [30,38]. For the Epworth district, wall around homesteads and industrial areas most likely acted as physical barriers for surface runoff, reducing slope length and affecting flow velocity and flow direction of surface runoff.
The spatial pattern of current soil erosion damage as documented in a diagram map is based of field survey in December 2019 ( Figures 5 and 6). The field measurements data serve as a validation tool for the estimated soil erosion risk modelling for the Epworth district. The plot-based field measurements recorded reflect the extent of erosion prevalent within the mapped area. The RUSLE model is presumed to predict the amount of soil moved on the field [126,127]; henceforth, the spatial damage data from the field measurements represented the proportions of soil moved. These mapped sites of the spatial extent damage map spatially concur with soil loss estimated using the RUSLE model. This, however, improves the model evaluation despite the lack of sheet erosion assessment. Even though, [128,129] reiterated that empirical modelling requires long-term field measurements and the analysis of sedimentation rates for validation purposes. The utilization of point-like plot-based data from field surveys for validation improves the evaluation of the model outputs and its understanding [55,130,131]. Validation of soil erosion risk modelling by comparison with outcomes of current damage mapping was further coupled with an analysis of the relationship between LULC and estimated soil erosion risk [62,100].
A total area of 765 m 2 was subjected to soil erosion in 2019 recorded for the 49 randomly selected sampling plots covering an area of 15.7 ha in the Epworth district. The randomized locations of the plot measurements indicate that eroded areas occurred in high frequency in the southwestern and southeastern parts of the Epworth district. Evidently, the spatial extent soil erosion damage in the diagram map indicates altogether low percentages of disturbance with a range of 0%-1.4% damaged area in the surveyed plots. However, during the field survey, the observed damage resulting from soil erosion ( Figure 6) appears to be greater in extent and magnitude than the depicted damage illustrated on the soil erosion damage map ( Figure 5). In comparison to the maps displaying the soil erosion risk based on the application of the RUSLE model, the spatial extent of soil erosion damage has been observed to spatially concur majorly in the southwestern and southeastern parts of the district where LS factor is presumably high. Areas of high surface damage could be identified in the southeast of the district and were predominantly observed in croplands and areas undergoing construction.
The field survey positively contributes to the study through the identification and registration of areas vulnerable to soil erosion ( Figure 5), which therefore lessens the burden of the resource strained land managers and local boards in developing conservation strategies direct on hotspots rather than concentrating on the entire district [5,132]. Areas of high soil erosion risk which occur in zones with high topographic (LS) factors can be confirmed by strong surface damages ( Figure 5). Nevertheless, the exclusion of sheet erosion recording and quantification during conducting damage mapping reduces the usability of the soil erosion damage map for the validation of the soil erosion risk mapping applying the RUSLE model [133]; this is because sheet erosion has a major contribution to erosion damage and is included in the RUSLE model [64,86,125]. However, the conclusions were drawn on the basis of field observations due to the heterogeneity of the urban set-up and the widespread nature of impervious surfaces.
The modelling findings reveal that topographic characteristics (LS factor) significantly influence potential erosion and soil erosion risk in the Epworth district, which concurs with the findings by [23] that model-based distribution patterns of soil erosion risk in the area of Windhoek, Namibia were mainly defined basing on the spatial structure of slope. The high soil erosion risk observed on sloppy areas in Epworth district corresponds to areas with convex to straight profile curvature and to the occurrence of ridges [134,135]. Furthermore, along the channel flanks modelled soil erosion risk was in general high, predominantly controlled by relief [122]. Due to the southward drainage of the stream network in Epworth district and thus, southward increasing incision of the streams, increased soil erosion risk along the river flanks can be predominantly observed and also on the southeast areas ( Figure 5). In contrast, the analysis of the effects of LULC change on soil erosion pointed out that increasing distribution of built-up areas as a result of high population pressure and demand for shelter substantially propels soil erosion risk within residential areas.

Conclusions
Soil erosion is a global environmental concern impacting negatively on agricultural productivity, accessibility to properties and also posing flooding risk in urban areas. The empirical RUSLE model was implemented to assess vulnerable areas and the computation of soil erosion risk through the integration of GIS and remote sensing techniques. The quantitative assessment of average annual soil loss for the Epworth district using the RUSLE model considers climate, soil, land use and topographic datasets as input parameters. Areas with high soil erosion risk were found to spatially correlate with topographic characteristics, especially slope length and slope steepness. The unrestricted LULC changes resulting from rampart informal settlements growth have accentuated soil erosion risk in the Epworth district. The analysis of LULC and estimated soil erosion risk improves the understanding of the spatial distribution patterns of soil loss for the different land uses in the years 2000 and 2018. The predicted soil erosion loss in the Epworth district amounted to 50,408 tons in 2000 while, 42,934 tons were estimated for 2018. Thus, the findings reveal that estimated soil erosion risk in total decreased over the study period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). This is attributed to the reduction of croplands and areas covered by green spaces at the expense of built-up areas. Soil loss massively increased in the residential areas from 20,236 tons in 2000 to 32,061 tons 2018, regardless of the concentration of built-up areas (concentrated and less concentrated residential areas); in total the area covered by residential areas almost doubled between 2000 and 2018. Increasing impervious surfaces, sealed areas and avoidance of paved areas during high traffic flow have been observed as contributing factors towards increased generation of surface runoff and hereby affected soil erosion risk in the growing residential areas of the Epworth district.
The soil erosion damage map generated from the field measurements served as a validation tool for the study as it revealed areas vulnerable to soil erosion within the Epworth district that concur with the results of the application of soil erosion risk models. The area affected by soil erosion in the surveyed plots showed damages of 0%-1.4% in the spatial extent of the mapped plot area. Therefore, field mapping data have been observed as necessary in ascertaining and improving an understanding of quantitative soil erosion modelling. Model validation demonstrated that the RUSLE model performance was good due to positive correlation between field measurements and model results basing on sample areas (r = 0.76, R 2 = 0.58, p < 0.05). It can be concluded that the spread of urban built-up areas without implementation of sound conservation practices, such as proper land suitability analysis and the construction of runoff drainage canals, will increase soil erosion damage by water. Although, the research has predicted potential erosion risk and soil erosion risk, it is important to outline that there are uncertainties with the modelled data provided by the RUSLE model arising from the lack of site-specific parameterization. Limited studies on water-induced soil erosion for urban areas within the region reduce options for data comparison. However, the computed soil erosion risk maps may assist environmental managers and land and policymakers on planning mitigatory measures for the study area.