Mapping of Soil Erosion Risk in Bukavu (Democratic Republic of Congo): Using RUSLE, Remote Sensing and GIS ()

1. Introduction
Soil erosion is one of the biggest global environmental problems resulting in both on-site and off-site effects. In developing countries, the economic implications of soil erosion are more serious because of lack of capacity to cope with it and also to replace lost nutrients. These developing countries have also high population growth which leads to intensified use of already stressed resources and expansion of production to marginal and fragile lands. Indeed, such processes aggravate erosion and productivity declines, resulting in a population-poverty-land degradation cycle. Worldwide, soil erosion associated with land degradation is spatio-temporal phenomena that are growing [1] [2] [3] [4]. This phenomenon of soil erosion is an unavoidable natural phenomenon that becomes a serious environmental and economic problem when it is accelerated by human activities [5] [6]. According to [7] [8], water erosion phenomena have accelerated all over the world.
However, water erosion is a natural geological phenomenon that has shaped the earth’s surface over geological ages [9]. Thus, climate change and landscape under the influence of demographic pressure and the extension of farming crops, have contributed to the increase in the exposure of the land to the runoff process, and consequently, to the degradation soils by erosion [10] [11]. Hence, various kinds of human activities such as agricultural practices, logging, grazing, construction of roads and buildings tend to modify the phenomena of erosion, often accelerating it considerably [12] [13]. Then, water erosion affects the quality of soils by inducing the deterioration of their surface layers rich in organic matter. Soil erosion leads to loss of soil productivity and land degradation.
However, the city of Bukavu is subject to strong constraints which negatively affect their hydrological functioning and make this topography a natural environment particularly sensitive to anthropogenic actions. Among these constraints, we particularly cite a very rugged topography, soils poor in organic matter and weakly thick, impermeable source rocks, friable substrates, anarchic construction, a harsh and sometimes brutal climate, and an often sparse vegetation cover, which makes it very exposed to the phenomenon of erosion [14] [15] [16].
The quantification of water erosion is normally carried out using models adapted to the intended purpose [17]. The mathematical models used in this quantification are continuously improved by researchers in many countries by adapting them to local conditions. According to [18], several models are currently developed and used in several countries and in different regions. Erosion prediction models help decision-makers in planning the management of natural and agricultural environments. Although it is difficult to find a model that considers all forms of erosion, some models have been developed to help water and soil managers identify the area’s most vulnerable to erosion. The RUSLE model [19] [20] is widely used in the Mediterranean region and in East Africa. It is an empirical model based on the USLE model [21]. The use of remote sensing and geographic information systems has made it possible to estimate and spatialize soil erosion at a reasonable cost [22] [23] [24] [25]. The principle consists in integrating in a geographical information system the factors appearing in the empirical model RUSLE
The objective of this work is to analyze and spatialize the various factors involved in the phenomenon of erosion and to arrive at the map of the risks of erosion and soil loss in the city of Bukavu. Indeed, the mapping of soil erosion factors and the identification of vulnerable areas would make it possible to assess the risks of erosion with a view to considering soil conservation and water resource management measures.
2. Methodology
2.1. Description of the Study Area
The area of the present study is located in the East of the Democratic Republic of Congo, Province of South Kivu in the city of Bukavu. The latter is located on a fracture point of the earth’s crust and between 2˚27'31'' and 2˚33'17.8'' South latitude and 28˚48'3.75'' and 28˚53'37'' East longitude. It is bounded to the north by Lake Kivu, to the west by the Nyamuhinga river which constitutes its border with Kabare territory, and to the south by the Mudusa groupement in Kabare territory and to the east by the Ruzizi river. The city is separated from Rwanda by Lake Kivu, as well as by the Ruzizi River [26] [27]. According to Chamaa et al. (1981) cited by [27], it has an area of 63 km2 including 43 km2 of dry land and 20 km2 occupied by the lake.
2.2. RUSLE Calculation
The factors causing erosion such as climate, soil properties, vegetation cover and management practices are considered for estimating soil loss. The RUSLE equation is a multiplicative function of five factors controlling the rill and inter-rill erosion [19] and can be expressed as:
(1)
Where:
A is the mean annual soil loss expressed in ton\ha*yr
R is rainfall and runoff erosivity index (in MJ*mm\ha*yr)
K is soil erodibility factor (in ton*ha*h/ha*MJ*mm)
LS is slope Steepness and slope Length factor (dimensionless)
C is the cover factor (dimensionless)
P is the conservation practice factor (dimensionless).
To achieve our objectives, different types of data were collected on the city of Bukavu: The Digital Terrain Model (DTM) will be useful to us in the characterization of flows (the length and the values of slopes, the accumulation of flows). These elements (Flow accumulation and Flow Direction) will be used in the calculation of the LS factor, which is one of the factors in the Soil loss equation. Landsats images: these images will be used in the characterization of the land cover. This will be useful to us in determining the C factor, given that each type of coverage corresponds to a corresponding C value. We will proceed by the classification of LandSat images to discriminate the different types of land use in the city of Bukavu. Data on the type of soil: each soil texture corresponds to a factor K (soil erodibility factor). Rainfall data: these data will be taken at points around the city of Bukavu. At each point considered, the corresponding R factor will be determined according to the following equation:
(2)
where:
R: erosivity of rain (Jm-2) and P: amount of annual precipitation (mm).
Being quantitative data on points around the city of Bukavu, the R factors corresponding to each point will be interpolated to transform them into Raster by the IDWI method thanks to “Geostatistical Analyst” of the ArcGIS software. The “GA Layer to Grid” conversion tool will be used for the conversion of the interpolated layer.
Before the interpolation, to have the dimensions of the pixels in meters, we will first project the point layer in WGS 84 UTM Zone 35 S through the “Projection and Transformation” tool. The K factor will be determined based on the particle size composition of the soil, the proportions of three elements of which will be taken into account: %
, %
and %,
. By adding to these three elements the proportion of Organic Matter (OM), the K Factor will then be determined by the equation of [28]:
where:
.
K is expressed in ha/MJmm, OM is the percentage of organic matter in the soil,
is the fraction of sand in the soil,
is the fraction of silt in the soil and fclay is thefraction of clay in the soil. The determination of the Factor P: This factor depends on land use practices and the implementation of any protection measures against soil erosion in the watersheds (such as reforestation, stabilization of slopes, vegetative erosion traps, etc.). Typical P-Factor values based on type of land use and erosion protection measures [21] are shown in Table 1.
To find the raster corresponding to the factor P, we will proceed as follows: Reclassify (using the Reclassify tool) the slope raster according to the 7 classes as listed in the table. Re-reclassify the previously created raster by assigning each slope class its P value.
The data will be processed in the ArcGIS Software using the following tools: mainly the Spatial Analysis tools will be taken into account in the various calculations: the conversion tools, thanks to “Polygon to Raster” will allow us to convert the soil layer into a Raster layer with the K Factor as an attribute. Thanks to
![]()
Table 1. Typical values of the P factor.
the “Raster to Polygon” tool, we will be able to convert the Land use raster into polygon (Vector). After adding the “C” field to this layer and completing each type of coverage with its corresponding value, the “Polygon to Raster” tool will help us convert the vector layer back to a raster. But first of all we will use the “Dissolve” tool to group the polygons of each category into a single group.
・ Hydrology tools:
○ Fill: to eliminate any gaps in the Digital Terrain Model
○ Flow Direction: to model the direction of flows on slopes
○ Flow Accumulation: to model the accumulation of flows from flows
・ The Surface tools: here we will use the “Slope” tool to generate the values of the slopes in degrees. This slope will then be transformed into radians by applying the equation:
(3)
With Sr: Slope in radians and Sl: slope in degrees
・ The Extraction tools: the “Extract by Mask” tool will be used to cut out the different images. The layer of Mask will be that of the study area, that is to say that of the administrative limit of the city of Bukavu;
・ “Map Algebra” tools: The Raster Calculator tool will be used to apply different equations:
○
(1)
○
(4)
○
(5)
○
(6)
Equation (4) remains valid for slopes less than 21%. Equation (5) will be applied to slopes with values greater than 21%. Equation (6) will be applied to combine the two equations by applying the “if” function.
Since 21 corresponds to the slope of 21%, we will use the corresponding of this slope in degrees (the tan of 0.21). This corresponds to 11.86˚.
The general calculation procedure in GIS is illustrated in the diagram below (Figure 1).
Table 2 lists the different data sources used in this study.
![]()
Figure 1. General calculation procedure in GIS.
![]()
Table 2. Data sources used in this work.
3. Results
3.1. R Factor
Rainfall data for the 34-year period were collected from eight points around and over the city of Bukavu as shown in Table 3.
In view of this table, it appears that the average annual rainfall is between 1465 and 1680 mm over the 35 years studied.
The R-factor interpolation results are shown in Figure 2.
3.2. K-Factor
According to SOTER soil data, the city of Bukavu represents only one type of soil (clay soil). This is how a single value of K will be determined for the entire map of Bukavu. According to the results of research on the soil of Bukavu, it would be made up of fine particles (below 75 microns) which would vary from 60-85%. Until now, no research has been able to discriminate between the different grain sizes in fine particles. By applying the equation of [28], we find the following as shown in Table 4.
![]()
Table 3. Rainfall data for the 35-year period (1986 to 2020) from eight points around and over the city of Bukavu.
![]()
Table 4. K-factor of fine particles.
The two proportions attributed to fsil and fclay come from the average between 60 and 85% of the proportions of the fine fractions of the soil of Bukavu as published in scientific works. Some publications attribute up to 5% the proportion of sand (fsand). Organic Matter can exist up to 0.105%.
3.3. K-Factor Mapping
Having the soil layer in vector format, we will create a field of K and assign its calculated value to it, then we will transform the vector layer into a Raster layer. With the “Polygon to Raster” conversion tool, the input layer will be the soil vector and the output layer will be the Raster of K which will be used in the calculation of the universal soil loss equation.
3.4. Characterization of Land Use/Land Cover
Figure 3 presents the characterization of Land Use and Land cover
Referring to the work of [29] assigning the different types of occupation the value of C, we also assigned to the five types of occupations of Bukavu the values as follows in Table 5.
![]()
Figure 3. Characterization of land user/land cover.
![]()
Table 5. Five types of occupations in Bukavu.
3.5. Factor C Mapping
We will then transform the vector layer of the Land Use back into a Raster by taking the Factor C as “Value” (using the “Polygon To Raster” Tool). Figure 4 presents the result of the characterization of factor C
On this map, the gaps (for example at Labotte) are due to the fact that the ground data (shapefile) comes from cartographic data at very small scales, while our work is done on a large scale. The C Factor The determination of Factor C is based on the characterization of Land use/Land Cover. This is how we will gradually proceed with the characterization of Land Use/Land Cover and then the mapping of Factor C. For each type of land cover, there is a value of C.
3.6. The P Factor
Figure 5 presents the result of the slope reclassification
The percentage and pixel count of slope classes presented in Table 6 below.
With regard to Table 6, the zone is dominated by slopes >5˚, i.e. more than 70% of proportion. After re-classify of P [21] for each class of slope, the result follows (Figure 6).
3.7. The LS Factor
3.7.1. Direction of Flows
Figure 7 shows the Flow Direction result. The torrential rivers crossing the city of Bukavu and the Ruzizi river.
3.7.2. Accumulation of Flows
Figure 8 presents the result of the Accumulation of Flows.
3.7.3. Slope
Figure 9 presents the result of the slope and applying Equation (3), we obtained the slope gradient.
3.7.4. Calculation of LS
The influence of the length of the slope and its inclination on the erosive processes expressed by the LS factor. For slopes whose value is less than 21%: in raster calculator we will enter
where, FA: Flow accumulation raster and Sr: slope Raster in Radians. The output Raster (output) we call it LS1. Figure 10(a) presents the map of the LS factor whose slope is less than 21% and greater than 21% (Figure 10(b)).
For slopes with a value greater than 21%: in raster calculator we will enter Power
where, FA: Flow accumulation raster and Sr: slope Raster in Radian. The output Raster (output) we call it LS2. We will now combine the two previously created rasters and apply the condition on the degree slopes raster. In the Raster Calculator, we’ll enter Con (“Slope” < 11.86, “LS1”, “LS2”). 11.86˚ corresponds to 21% slope. The Output Raster (output) we call it LS end. This is the final LS raster to be used in the soil loss calculation equation.
3.8. Universal Soil Loss Equation Mapping
The rasters of 5 factors involved in the universal soil loss equation having already been found, we will now calculate the final raster by calculating the product of these five factors with the “Raster Calculator” tool. The erosivity map thus produced will be classified into 5 risk levels as follows (Table 7).
We will reclassify the map to be able to find the proportions of each level of risk.
Figure 11 presents the erosivity of the city of Bukavu.
The pixels are very coarse (331 m) since the erosivity raster is derived from the product of several rasters (of 5 factors) with different spatial resolution levels (Table 8).
(a)
(b)
Figure 10. (a) Map of the LS factor whose slope is less than 21%; (b) Map of the LS factor whose slope is greater than 21%.
This previous table represents the distribution of pixels according to the different levels of risk over the entire extent of the city of Bukavu. In the calculations, although some parts of the outskirts of the city are mapped, we have taken into account only the part strictly covering the city of Bukavu. On this we have cut the erosivity raster taking as Mask the administrative limit of the city of Bukavu. However, we can analyze the levels of risk according to the municipalities (Table 9)
4. Discussion
4.1. The Topographic Factor (LS)
The slopes and topographic factor determine the area from the 21% slope value. Below this slope, the erosion is very weak, essentially characterized by the form of diffuse erosion or in weak incisions. On the other hand, beyond 21% slope, erosion takes the form of immense ravines. In the details, it has been recognized, as [30] [31], zones of low erosion on slopes less than 8%, that of 8% to 12% as a zone of diffuse runoff and rills, and the zone of slope greater than 12% as an area of heavy erosion, essentially gully erosion from which the city of Bukavu suffers gully erosion.
4.2. Rain Erosivity (R)
The highest R values of 256.730943 MJ∙mm/ha∙h∙an, are located in the city of Bukavu. Overall, the erosivity decreases going up from this zone towards the highest relief of the catchment area. It is necessary to see in this variation the impact of the Foehn effect well known to climatologists by which the rains decrease while rising on the relief because of the loss of humidity. Indeed, the rains that fall in this part seem to indicate a north-east direction on the plain towards south-west in the hilly area. The combination of this energy with the various degradations of the environment is at the origin of the strong erosion in these zones. On the other hand, [32] found higher values in the Lake Kivu basin (average R = 4623 MJ∙mm/ha∙h∙year) where the average annual rainfall is 1285 mm.
4.3. Soil Erodibility Index (K)
This index is low on the Labotte. It does not reflect the susceptibility of clay soils in the region. However, these are only values from the first estimate based solely on soil texture. Further calculation integrating soil texture with soil organic matter content, soil structure and permeability, not included here, could yield arguably higher values. Erodibility varies according to soil types. [32] found K values between 0.009 and 0.11 t∙ha∙h/ha∙MJ∙mm in the Lake Kivu basin.
4.4. The Cultural Factor (C)
The basin includes areas of urban habitat, natural forests and areas of rural exploitation. [32] found 0.15 in the Lake Kivu basin. Estimating the C factor is not easy. In the absence of real agricultural development, the calculation of this factor is based on land occupation.
4.5. The P Factor
There are no anti-erosion devices voluntarily installed in the city of Bukavu. This factor can only be estimated indirectly. However, some traditional agricultural practices such as ridges or flat cultivation are used by the population in relation to the slope of the land. These cultural methods could be considered compatible with soil conservation. Areas with a high P value are those where anti-erosion devices are essential. The data obtained indicate that these areas coincide with spaces with high erosion sensitivity.
5. Conclusion
In conclusion, it can be said that soil erosion risk map has been created from the combination of many parameters interacting each other in a complex way generating the final quantitative RUSLE values. The data showed that factors affecting in a stronger way the erosion process are the one taking into account the topography (LS factor), and support practices (P factor), followed by the Cover parameter (C factor). This study is carried out in the city of Bukavu, with the aim of analyzing and spatializing the various factors involved in the phenomenon of erosion and of arriving at the map of the risks of erosion and soil loss in the city of Bukavu for sustainable land management. Concretely, it is a question of using the universal land loss equation (USLE), coupled with GIS and remote sensing. The results give the orders of magnitude of the parameters generating erosivity of 256.730943 MJ∙mm/ha/h/year while the erodibility (K) goes up to about 0.2. With regard to potential erosion, the calculated land losses are between 2 and 50 t/ha/year. Data analysis shows that more than 70% of the surface of the city of Bukavu is within the soil loss tolerance threshold in the intertropical region, i.e. 20 t/ha/year. However, it appears that 1% of the watershed is gullied in the parts less covered by vegetation, especially in areas occupied by uncontrolled urbanization. Thus, at the end of this study, some erosion control methods are recommended in order to conserve the soil.