Evaluation of the Impact of Land Use Changes on Soil Erosion in the Tropical Maha Oya River Basin, Sri Lanka

Abstract

Soil degradation is a serious environmental issue in many regions of the world, and Sri Lanka is not an exception. Maha Oya River Basin (MORB) is one of the major river basins in tropical Sri Lanka, which suffers from regular soil erosion and degradation. The current study was designed to estimate the soil erosion associated with land use changes of the MORB. The Revised Universal Soil Loss Equation (RUSLE) was used in calculating the annual soil erosion rates, while the Geographic Information System (GIS) was used in mapping the spatial variations of the soil erosion hazard over a 30-year period. Thereafter, soil erosion hotspots in the MORB were also identified. The results of this study revealed that the mean average soil loss from the MORB has substantially increased from 2.81 t ha⁻¹ yr⁻¹ in 1989 to 3.21 t ha⁻¹ yr⁻¹ in 2021, which is an increment of about 14.23%. An extremely critical soil erosion-prone locations (average annual soil loss > 60 t ha⁻¹ yr⁻¹) map of the MORB was developed for the year 2021. The severity classes revealed that approximately 4.61% and 6.11% of the study area were in high to extremely high erosion hazard classes in 1989 and 2021, respectively. Based on the results, it was found that the extreme soil erosion occurs when forests and vegetation land are converted into agricultural and bare land/farmland. The spatial analysis further reveals that erosion-prone soil types, steep slope areas, and reduced forest/vegetation cover in hilly mountain areas contributed to the high soil erosion risk (16.56 to 91.01 t ha⁻¹ yr⁻¹) of the MORB. These high soil erosional areas should be prioritized according to the severity classes, and appropriate land use/land cover (LU/LC) management and water conservation practices should be implemented as recommended by this study to restore degraded lands.

1. Introduction

Soil erosion has become a severe environmental issue due to the degeneration of world food production [1]. Erosion has resulted in the loss of productive land, thereby causing a decrease in agricultural output. Sediment delivery causes non-point source pollution and water quality degradation [2]. It not only harms the environment in erosion-prone areas, and also the neighborhoods [3]. The deterioration of ecosystems caused by changes in land use affects riverine environments, risking the existence of a large variety of species [4]. The two primary agents of soil erosion are water and wind. The prominent factors influencing the severity and the rate of soil erosion are topography, climate characteristics, deforestation, soil characteristics, land changes, and urbanization with high population expansion [5,6]. The Global Assessment of Human-induced Soil Degradation (GLASOD) study revealed that, globally, 56% of the land areas (excluding cold, wet, and desert regions) were susceptible to water erosion, while 28% of the land zones were vulnerable to wind erosion [7]. Water erosion is categorized into various types, such as inter-rill erosion, rill erosion, gully erosion, channel erosion, sheet erosion, and bank erosion [8,9].

Globally, the mean annual soil loss varies between 12 and 15 t ha⁻¹ yr⁻¹ [10]. The study done in Cambodia, which used the Revised Universal Soil Loss Equation (RUSLE) model, revealed that the mean annual soil erosion rate ranges between 3.1 to 7.6 t ha⁻¹ yr⁻¹ in 2002 and 2015 respectively, in agricultural land, barren land, built-up areas, paddy field, and forest [11]. According to Piyathilake et al. [12], tolerable soil erosion rate rates in the up-country wet zone, mid country wet zone, and low country dry zone area in Sri Lanka are 13.2 t ha⁻¹ yr⁻¹, 9 t ha⁻¹ yr⁻¹, and 6.7 t ha⁻¹ yr⁻¹, respectively. The mean annual soil erosion rate in 2017 in the Kalu River Basin of Sri Lanka was 0.63 t ha⁻¹ yr⁻¹ [13], which is less than the mean tolerable soil loss rate of Sri Lanka, which is 5 t ha⁻¹ yr⁻¹. The study done by Perera et al. [14] revealed that the Moneragala District in Sri Lanka had 27.8 t ha⁻¹ yr⁻¹ mean annual soil loss rate in 2020, which is about four times higher than the tolerable rate of soil erosion in the low country dry zone (6.7 t ha⁻¹ yr⁻¹). The soil loss rate of the Bibili Oya watershed was 12.9 t ha⁻¹ yr⁻¹ in 2014 [15], which is significantly higher than the soil loss tolerance limit (9 t ha⁻¹ yr⁻¹) of the mid country wet zone. The Kotmale watershed had a 9.8 t ha⁻¹ yr⁻¹ mean annual soil loss in 2018 [16], which is within the soil loss tolerance limit in the mid country wet zone. The erosion assessment done in Sabaragamuwa province in Sri Lanka revealed that estimated mean annual soil erosion was 15.53 t ha⁻¹ yr⁻¹ in 2019 [17], which is greater than the tolerable soil loss limit (9 t ha⁻¹ yr⁻¹). The study done by Fayas et al. [18] demonstrated that the estimated mean annual soil erosion (10.9 t ha⁻¹ yr⁻¹) in the Kelani River Basin is higher than the soil loss tolerance of Sri Lanka, which is 5 t ha⁻¹ yr⁻¹. De Silva et al. [19] revealed that the mean annual soil erosion rate of the Nalanda Oya catchment in Sri Lanka was 2.99 t ha⁻¹ yr⁻¹, which is far less than the soil loss tolerance limit (9 t ha⁻¹ yr⁻¹). Piyathilake et al. [20] posited that the estimated mean soil erosion of 25.6 t ha⁻¹ yr⁻¹ in Uva province, Sri Lanka was approximately three times the soil tolerance (9 t ha⁻¹ yr⁻¹) of the mid country wet zone.

In the early stages, soil erosion assessments were done by visiting the field and by collecting samples [11]. It is a complicated and time-consuming process to evaluate soil erosion in large watersheds [21,22]. With the advancement of numerical methods, different tools are available to evaluate soil erosion nowadays. They are physical-based, empirical, and conceptual models [13]. Water Erosion Prediction Project (WEPP) is utilized frequently among physical-based models [20]. Agricultural Non-Point Source Pollution (AGNPS) is one of the conceptual models developed to evaluate soil erosion from agricultural regions [23,24]. Among others, Erosion Potential Model (EPM), Coordination of Information on the Environment (CORINE), and Kinematic Runoff and Erosion Model (KINEROS) are the globally used models to estimate soil erosion [20]. The Universal Soil Loss Equation (USLE) and Revised Universal Loss Equation (RUSLE) are the frequently used empirical models to estimate soil erosion. They have been used in more than 109 countries under various environments [25,26]. In Sri Lanka, more than 50% of studies used the RUSLE empirical model, whereas 29% and 21% of studies employed the USLE and InVEST-SR/SDR approaches accordingly [20]. In addition, Ouma et al. [27] used the Soil Degradation Index (SDI) and Normalized Difference Railway Erosivity Index (NDReLI) methods to predict susceptibility of soil erosion within railway corridors in their study and their results revealed that the NDReLI method is more accurate than RUSLE and SDI models when considering the soil erosion within railway corridors.

The RUSLE model was used in this study to evaluate the soil erosion in the Maha Oya River Basin (MORB). Most of the recent studies in Sri Lanka applied the same soil erosion model in their studies: Nalanda Oya catchment [19]; Kotmale watershed [16]; Kalu Ganga River Basin [13]; Bibili Oya watershed [15]; Sabaragamuwa Province [17]; Anuradhapura, Polonnaruwa, and Vavuniya Districts [28]; and Kelani River Basin [18]. However, the RUSLE model can estimate soil erosion considering only the mass balance, which cannot evaluate the sediments deposited in the riverine and gully erosion [1,29]. The advantages of the RUSLE model are simplicity, data retrieval, efficiency, cost-effectiveness, greater reliability, and pertinency over large areas [30,31,32,33,34].

During the last three decades, the amount of forest and vegetation in Sri Lanka has decreased and the amount of land used for agriculture has increased. The studies done in the MORB mainly focused on water quality management of the river, sand mining, natural resources in the river basin, valuation of the river′s ecosystem services, and climate changes in the river basin. Seevarethnam et al. [35] found that the MORB suffers from regular soil erosion and degradation. However, even though there are some studies done for different river basins in Sri Lanka, there is no any previous research carried out for the MORB. Hence, as per the authors′ knowledge, the current study is the first which comprehensively examines the soil erosion and the possible effects of the land use/land cover (LU/LC) on the soil erosion levels in the Maha Oya basin.

The purpose of our study is to evaluate the land use/land cover (LU/LC) changes between 1990 and 2020 and how they affected soil erosion. The study′s specific objectives were to: (1) estimate the annual soil erosion in the MORB in both extent and spatial distribution; and (2) ascertain how different forms of land use and land cover contribute to soil erosion in the MORB. Estimation of soil erosion quantities in the MORB will benefit researchers, engineers, and environmentalists since the results of this study can be used to identify the erosion-prone areas and erodible LU/LC types of the basin and implement proper soil, water, and environmental conservation practices.

2. Materials and Methods

2.1. Study Area

The Maha Oya River flows through three provinces (Sabaragamuwa, Central, and Western) and four districts (Kurunegala, Gampaha, Kegalle, and Puttalam) in Sri Lanka (refer to Figure 1). The river basin is located between 7°0′0″–7°30′0″ N and 79°50′0″–80°35′0″. It has a length of about 135 km. The total drainage area is about 1521 km². The river flows from Nawalapitiya in the Aranayake region to Kochchikade, North of Negombo, where it meets the Indian Ocean. Approximately 75% of MORB is in the wet zone. The other 25% lies in the intermediate zone of Sri Lanka [36]. The MORB provides a livelihood for nearly 1 million people. The primary economic sector of the study area is agriculture [37]. The main land use types in the MORB are forest, vegetation, bare land, settlement, exposed rock, and water bodies.

2.2. Rainfall Data

The monthly rainfall data from 1990 to 2020 (30 years) were obtained from the Meteorological Department of Sri Lanka for five rain gauging stations located within the study area. These stations are Aranayake Govt. Hospital, Eraminigolla, Ambepussa Govt. Farm, Andigama Farm, and Mellawa Estate.

Table 1. Location of rain gauging stations in Maha Oya River Basin.

Rain Gauging Station	Location Coordinates
Ambepussa Govt. Farm	7°16′48″ N       80°10′12″ E
Andigama Farm	7°22′12″ N       80°7′12″ E
Aranayake Govt. Hospital	7°10′48″ N       80°28′12″ E
Eraminigolla	7°17′60″ N       80°22′48″ E
Mellawa Estate	7°19′12″ N       79°57′0″ E

shows the locations of rain gauging stations in the MORB.

Figure 2 shows the annual rainfall of 26 years for five rain gauging stations. The highest and minimum annual precipitation of 3022.7 mm and 856.5 mm were recorded in year 1993 at Eraminigolla station and in 1999 at Aranayake Govt. Hospital station, respectively.

2.3. Soil Data

The digital soil map of Sri Lanka was obtained from the Natural Resource Management Centre (NRMC) of Sri Lanka. Figure 3 shows the major soil types found in the area. The Red-Yellow Podzolic soil has the highest percentage coverage of the study area. The major soil types in MORB (Figure 3) are “Red-Yellow Podzolic soils; steeply dissected, hilly and rolling terrain″, “Red-Yellow Podzolic soils & Low Humic Gley soils″, and “Red-Yellow Podzolic soils with soft or hard laterite″ which cover 29.61%, 21.70%, and 16.50% of the total area, respectively. The texture of these soil types is sandy clay loam [38], susceptible to soil erosion [39,40].

2.4. Digital Elevation Model (DEM)

This study used Shuttle Radar Topography Mission (SRTM) DEM of 30 m spatial resolution. The DEM was obtained from United States Geological Survey (USGS) website (https://earthexplorer.usgs.gov, accessed on 12 May 2022).) and processed using Mosaic dataset tool in ArcGIS 10.5. According to the Nut et al. [11] slope classification, the study area was classified into six slope classes. The highest percentage of area (23.87%) was occupied by the low susceptibility (2–5°) sloping in the entire study region. The lowest percentage (2.63% area) has extremely high-steepness slopes (>30°).

Table 2. Slope classification of the study area with related features.

Slope Class (Degrees)	Susceptibility	Characteristics	Area (km2)	%
0–2	Very Low	Flat To Very Gently	290.62	19.12
2–5	Low	Sloping	363.03	23.87
5–10	Medium	Gently Sloping	331.24	21.78
10–15	High	Strongly Sloping	205.61	13.52
15–30	Very High	Moderately Sloping	290.22	19.08
>30	Extremely High	Steep	39.98	2.63
Total			1520.7	100

shows the characteristics of the classified slope classes in the study area.

Figure 4 shows the elevation map and slope variation of the study area in degrees (°). Higher slopes as expected can be seen upstream while the lower slopes are found downstream of the river basin.

2.5. Land Use/Land Cover Classification

The LU/LC maps were developed using Landsat images that were obtained from USGS Earth Explore website (https://earthexplorer.usgs.gov accessed on 15 November 2022). According to the study period (1990–2020), the Landsat images were obtained with a 10-year frequency. Since there were no Landsat images with low cloud cover available in 1990, 2010, and 2020 for the study area, the 1989, 2009, and 2021 years were taken as reference years, respectively, for the LU classification. The satellite image types and relevant years are shown in

Table 3. Landsat images and years of extraction.

Base Year	Landsat Type
1989	Landsat 4–5 TM C1 Level-1
2000	Landsat 7 ETM + C1 Level 1
2009	Landsat 4–5 TM C1 Level-1
2021	Landsat 8 OLI/TIRS C1 Level 1

.

The supervised maximum likelihood LU/LC classification was done using the pertinent Landsat images that were clipped to the Maha Oya River Basin after generating the composite raster. The composite raster files were created using seven bands of each Landsat image. The obtained data were then classified into seven LU/LC classes based on the visual interpretation of the Landsat images. Those classes are vegetation, waterbodies, forest, agricultural land, bare land, exposed rocks, and settlements. The features of above LU/LC categories are given in Appendix D.

Accuracy Assessment

In addition, accuracy assessment was carried out using ground truth data (ground truth data was obtained from interpreting high-resolution imagery), after finishing the LU/LC classification. To examine the accuracy, 106, 109, 113, and 131 ground truth points were used for 1989, 2000, 2009, and 2021 respectively. The accuracy was evaluated using indicators such as ‘User Accuracy′, ‘Producer Accuracy′, ‘Overall Accuracy′, and Kappa Coefficient [41,42].

2.6. Estimation of Factors in the RUSLE Model

The Revised Universal Soil Loss Equation (RUSLE) is more widely used than the Universal Soil Loss Equation (USLE) to estimate long-term average annual soil loss of river basins in Sri Lanka and beyond [18,20]. The USLE was developed originally for gently sloping agricultural circumstances, but the RUSLE model broadens the models′ use to include forest, disturbed areas, and steep-slope soil loss scenarios [15]. RUSLE represents the impact of both physical factors and anthropogenic factors on rill and sheet soil erosion caused by surface runoff and raindrops [16]. In this study, the RUSLE framework was used with the Arc GIS version 10.5 environment since it was mostly recommended by other studies [5,13,21,43,44]. The RUSLE model [45] can be mathematically denoted as Equation (1):

$$\mathrm{A} = \mathrm{R} \times \mathrm{K} \times \mathrm{LS} \times \mathrm{C} \times \mathrm{P}$$

(1)

where A is average annual soil loss (t ha⁻¹ yr⁻¹), R is rainfall-runoff erosivity factor (MJ mm ha⁻¹ h⁻¹ yr⁻¹), K is soil erodibility factor (t ha h MJ⁻¹ mm⁻¹), L is slope length factor (dimensionless), S is slope steepness factor (dimensionless), C is crop/cover management factor (dimensionless), and P is erosion control/support practice factor (dimensionless). The overall methodology is given in Figure 5 as a flowchart.

2.6.1. Rainfall-Runoff Erosivity (R) Factor

The power of rainfall that can cause erosion for a specific region can be evaluated using the erosivity factor and it depends on the intensity of the precipitation [13]. There are several methods for estimating the average annual rainfall of a river basin. The researchers tend to use several interpolating techniques including the Thiessen polygons, Inverse Distance Weighted (IDW) technique, and Kriging method. When the Thiessen polygon method is used, the average annual rainfall is uniform throughout the river basin. In practice, the accuracy of the Thiessen polygon method may be reduced because annual precipitation in a river basin may not be uniform along the catchment area. Hence, there are two types of widely used methods by researchers available to estimate the mean annual precipitation. Those are IDW method [9,16,46] and Kriging method implemented in geostatistical wizard using a semi-variogram model for spherical data [15,19,47]. The IDW method has lesser interpolation error than Kriging method [13,18]. Therefore, IDW method was used in this study to interpolate average annual rainfall data, instead of Kriging method.

The equation for estimating the R factor was originally proposed by [48]. In this study, the original equation was not used due to lack of rainfall intensity data. Therefore, the formula developed by Premalal in 1986 [49] for the Sri Lankan condition was used in this study to estimate the R Factor (Equation (2)). The same formula was used in Sabaragamuwa province [17], Kelani River basin [18], Kotmale watershed [16], Uva province [12], and for entire Sri Lanka [50,51].

$$\mathsf{R}=\frac{\left(972.75+9.95\times \mathrm{F}\right)}{100}$$

(2)

where R is rainfall-runoff erosivity factor (MJ mm ha⁻¹ h⁻¹ yr⁻¹) and F is mean annual precipitation (mm). In this study, R factor was estimated using spatially distributed average annual rainfall from 1991 to 2020 utilizing the IDW method, and generated 30 m cell resolution maps as shown in the results section. Panditharathne et al. [13] asserted IDW method as an accurate interpolation technique.

2.6.2. Soil Erodibility (K) Factor

The soil erodibility factor (t ha h MJ⁻¹ ha⁻¹ mm⁻¹) is generally used to measure the susceptibility to and involvement of soil types with soil erosion [52]. The vector/polygon version of the Sri Lanka soil map obtained from NRMC was used to mask and extract the soil type map of the MORB. Then, using relevant K factor values obtained from the literature, the 30 m resolution K factor raster map was generated in ArcGIS 10.5 environment. The obtained K factor values for soil types in the study area [12,13,18] and the K factor map of the study area are shown in the results section.

2.6.3. Slope Length and Steepness (LS) Factor

One of the most important parameters in the RUSLE model is LS factor, which can estimate the soil loss rate considering the surface runoff with gravitational forces [11,53]. According to the Chen et al. [47], slope length has a significant effect on soil erosion. The LS factor can be derived using manual methods for small water sheds. The large and various topographical watersheds/basins require some accurate methods to derive the LS factor and a lot of methods with technology based on GIS/RS were developed in the world [54]. In this study, the slope length and steepness factors were evaluated using 30 m resolution SRTM DEM obtained from USGS website (https://earthexplorer.usgs.gov, accessed on 12 May 2022). The DEM was clipped according to the MORB shape file. Thereafter, we filled the sinks of the DEM and generated flow accumulation and flow direction raster maps using tools of ArcGIS 10.5 version. The slope length factor was generated by using the equations (Equations (4)–(6)) posited by Liu et al. [55]. The formulas for LS factor estimation in the RUSLE model are as follows:

$$\mathrm{LS}= \mathrm{L}\times \mathrm{S}$$

(3)

$$\mathrm{L}={\left(\lambda /22.13\right)}^m$$

(4)

$$m= \beta /\left(1+\beta \right)$$

(5)

$$\beta =\left(\mathrm{Sin}\theta /0.0896\right)/[3\times {\left(\mathrm{Sin}\theta \right)}^{0.8}+0.56]$$

(6)

$$\mathrm{S}=10.8\times \mathrm{Sin}\theta +0.03; \theta <9\%$$

(7)

$$\mathrm{S}=16.8\times \mathrm{Sin}\theta -0.5; \theta \ge 9\%$$

(8)

where λ is horizontal slope length in meters (λ = flow accumulation × cell size) and m is slope length exponent [45]. β is the ratio of rill erosion (caused by flow) to sheet erosion (generally caused by the impact of rain drops) [46,55], and θ is slope angle in degrees [48,56]. The same equations were used by Fayas et al. [18], Panditharathne et al. [13], and Thuraisingham & Weerasinghe [15].

In this study, the slope length factor was determined by Equations (4)–(6) that were employed in the raster calculator with the input of flow accumulation raster, using ArcGIS 10.5 version. The slope steepness factor was obtained by Equations (7) and (8) using slope (in degrees) raster file as input. In ArcGIS, all trigonometric functions are in radians, so when adding equations to the raster calculator, the slope (degrees) must be converted into radians. Following that, the LS factor was determined by multiplying the L and S factors using the raster calculator.

2.6.4. Cover Management (C) Factor

The C factors represent the soil loss from an agricultural system, compared with loss from land that has been fallowed and tilled continuously [48]. Under specific conditions, the C factor estimates how much soil erosion results from cropped lands. In terms of reducing soil loss, C factor is used to examine the relative efficacy of agricultural and soil conservation approaches [57]. The vegetation cover is most important in reducing soil erosion by protecting soil from raindrops and surface runoff, while increasing the soil properties [58].

The C factor values for this study were obtained from literature performed for the Sri Lankan circumstances and assigned to all the land use maps generated for the study area. The ArcGIS 10.5 was used to input the values to the attribute table of LU/LC polygons and the C factor maps were generated. The obtained C factor values and corresponding LU/LC types are denoted in Table 4.

2.6.5. Support Practice (P) Factor

The support practice factor is typically related to the C factor since both are used to manage the negative impact on soil erosion. However, the P factor differs from the C factor in that it indicates how management methods minimize runoff and change its pattern, orientation, and velocity through the control of runoff [48]. The P factor values for this study were obtained from previous studies done in Sri Lankan conditions and the same procedure used to generate the C factor map was performed to generate the P factor map. The obtained P factor values are shown in Table 4.

Table 4. C Factor and P Factor values used for different LU/LC types in MORB. Source: [12,16,18,49,50,59,60].

LU/LC Type	C Factor	P Factor
Water Bodies	0.2	0
Agricultural Land	0.43	0.15
Forest	0.5	0.3
Bare Land	1	1
Vegetation	0.51	1
Exposed Rock	0.1	0
Settlement	0.73	0

Table 4. C Factor and P Factor values used for different LU/LC types in MORB. Source: [12,16,18,49,50,59,60].

LU/LC Type	C Factor	P Factor
Water Bodies	0.2	0
Agricultural Land	0.43	0.15
Forest	0.5	0.3
Bare Land	1	1
Vegetation	0.51	1
Exposed Rock	0.1	0
Settlement	0.73	0

2.6.6. Spatial Analysis of Land Use and Land Cover Changes

The spatial variations of land use and land cover in the MORB from 1989 to 2021 were evaluated using the intersect tool in ArcGIS 10.5 version and the LU/LC change matrix was generated in tabular format with help of Microsoft Excel. The Intersect tool can be used to determine the geometric intersection of any number of feature layers or feature classes [61].

3. Results

3.1. Accuracy Assessment of LU/LC Classification

According to the results, the overall accuracy for LU/LC classification in year 1989, 2000, 2009, and 2021 are 91.51%, 85.32%, 90.27%, and 96.95%, respectively. The Kappa coefficient for 1989, 2000, 2009, and 2021 are 89.84%, 80.72%, 88.24%, and 96.33%, respectively (

Table 5. Accuracy assessment of LU/LC classification in year 1989, 2000, 2009, and 2021.

LU/LC Types	1989		2000		2009		2021
	Users Accuracy (%)	Producers Accuracy (%)	Users Accuracy (%)	Producers Accuracy (%)	Users Accuracy (%)	Producers Accuracy (%)	Users Accuracy (%)	Producers Accuracy (%)
Vegetation	100.00	81.82	97.44	80.85	100.00	77.42	100.00	100.00
Water bodies	100.00	100.00	100.00	94.44	86.67	92.86	100.00	100.00
Forest	92.00	88.46	57.14	100.00	100.00	100.00	100.00	100.00
Agricultural Land	61.54	100.00	75.00	81.82	95.83	95.83	100.00	100.00
Bare Land	80.00	80.00	100.00	50.00	66.67	100.00	75.00	75.00
Exposed Rocks	100.00	83.33	50.00	100.00	87.50	77.78	88.89	100.00
Settlements	96.00	100.00	76.00	95.00	77.27	100.00	96.88	100.00
Overall Accuracy (%)	91.51		85.32		90.27		96.95
Kappa Coefficient (%)	89.64		80.72		88.24		96.33

). According to Congalton et al. [41], accurate classification should be equal to or greater than the Kappa coefficient percentage of 80%. The equations and obtained sample values are given in Appendix B. According to the LU/LC assessment shown in

Table 6. LU/LC types of the MORB and corresponding areas in 1989, 2000, 2009, and 2021.

LU/LC Type	1989		2000		2009		2021		Net Change
									(1989–2021)
	Area	Area (%)	Area	Area (%)	Area	Area (%)	Area	Area	Area	Area (%)
	(km2)		(km2)		(km2)		(km2)	(%)	(km2)
Water Bodies	20.5	1.35	14.9	0.98	18.4	1.22	12.1	0.80	−8.4	−0.55
Bare Land	2.7	0.18	9.8	0.64	86.1	5.66	120.3	7.91	117.6	7.73
Agricultural Land	268.6	17.66	216.7	14.25	235.3	15.47	240.5	15.82	−28.1	−1.84
Settlement	64.6	4.25	69.6	4.58	83.4	5.48	135.2	8.89	70.6	4.64
Forest	566.7	37.27	542.4	35.67	569.0	37.42	500.7	32.93	−66.0	−4.34
Exposed Rock	29.7	1.95	20.9	1.37	30.6	2.01	40.2	2.63	10.5	0.68
Vegetation	567.9	37.34	609.3	40.07	497.9	32.74	471.7	31.02	−96.2	−6.32
Clouds	0.0	0.00	37.1	2.44	0.0	0.00	0.0	0.00	0.0	0.00
Total	1521.0	100.00	1521.0	100.00	1521.0	100.00	1521.0	100.00	0.0	0.00

, most (greater than 30%) of the study area are occupied by forests in each year. The forest cover decreased from 37.27% to 32.93% between the years 1989 and 2021. Settlements increased from 4.25% to 8.89%, while vegetation cover declined from 37.34% to 31.02% between 1989 and 2021. Exposed rock increased from 1.95% to 2.63% during the period between 1989 and 2021. Agricultural land decreased from 17.66% to 15.82% and bare land increased from 0.81% to 7.91%, respectively, from the year 1989 to 2021, while water bodies decreased from 1.35% to 0.80%. The temporal variation of each LU/LC class is shown in Figure 6 as a bar chart. Figure 7 indicates the generated LU/LC maps for year 1989, 2000, 2009, and 2021, respectively.

3.2. Spatial Variations in Land Use and Land Cover

The land use variations along the MORB were evaluated from year 1989 to 2021. The generated land use and land cover (this study defines bare land as land that is prepared for cultivation) change matrix in this study is shown in

Table 7. LU/LC change matrix of MORB from 1989 to 2021 (km²).

LU/LC Type		2021							Total 1989
		Agricultural Land	Bare Land	Exposed Rock	Forest	Settlement	Vegetation	Water Bodies
1989	Agricultural land	105.5	23.3	1.7	49.7	13.6	74.1	0.6	268.6
	Bare land	0.1	0.3	0.2	0.1	1.2	0.5	0.2	2.7
	Exposed rock	0.7	3.5	5.8	9.4	2.8	7.4	0.0	29.7
	Forest	47.3	42.9	25.2	270.4	47.9	132.3	0.8	566.7
	Settlement	6.2	8.4	2.9	6.9	23.5	15.2	1.5	64.6
	Vegetation	77.3	38.6	4.4	163.5	44.3	238.9	1.0	567.9
	Water Bodies	3.3	3.4	0.0	0.7	1.8	3.4	7.8	20.5
	Total in 2021	240.5	120.3	40.2	500.7	135.2	471.7	12.1	1520.7

. Table 7 demonstrates the area of each LU/LC type change between 1989 and 2021. The diagonal values indicated in bold and italic in Table 7 represent the unchanged area (km²) of corresponding LU/LC types in the study region. The change matrix revealed that about 652.2 km² of the study area has not changed during the period between 1989 and 2021. The 42.9 km² of forest in 1989 changed into bare land, and 47.3 km² of forest changed into agricultural land in 2021. The forest, agricultural land, and vegetation in 1989 were changed into settlements by 47.9 km², 13.6 km², and 44.3 km² in 2021, respectively.

3.3. Spatial Variations in Mean Annual Rainfall

As described in the methodology, the spatial variations in average annual rainfall of the MORB were evaluated using the ArcGIS 10.5 environment, utilizing the IDW method. The estimated mean annual precipitations at each rain gauging station are shown in

Table 8. Mean annual precipitation of the rain gauging stations in Maha Oya River Basin.

Station	Longitude	Latitude	Mean Annual Precipitation (mm)
Ambepussa Govt Farm	80.17	7.28	2170.80
Andigama Farm	80.12	7.37	1820.38
Aranayake Govt. Hospital	80.47	7.18	1878.20
Eraminigolla	80.38	7.30	2038.60
Mellawa Estate	79.95	7.32	1651.17

. According to Figure 8, the average annual rainfall of the MORB varies from 1653.17 mm to 2170.80 mm.

3.4. RUSLE Factors

3.4.1. Rainfall-Runoff Erosivity (R) Factor

The values of the R factor, generated using ArcGIS, ranged from 166.32 to 217.92 MJ mm ha⁻¹ h⁻¹ yr⁻¹. The obtained R values in this study are lower than the Kelani River basins′ R values (232.37–431.20 MJ mm ha⁻¹ h⁻¹ yr⁻¹), which are situated in the wet zone in Sri Lanka [18]. As shown in Figure 8, the region that obtained higher average annual rainfall (2170.80 mm) had higher R factor values while the area with low average annual precipitation (1651.17 mm) had low R factor values. According to that, the R factor is mainly dependent on the mean annual precipitation.

3.4.2. Soil Erodibility (K) Factor

The upland areas of the MORB were occupied by soil types of “Rock knob plain″, “Red-Yellow Podzolic soils; steeply dissected″, “Red-Yellow Podzolic soils & Mountain Regosols″, “Steep rockland & Lithosols″, “Red-Yellow Podzolic soils & Low Humic Gley soils″, “Erosional remnants (Inselbergs)″, “Reddish Brown Latosolic soils″, and “Immature Brown Loams″, which covered the basin area by 0.16%, 29.61%, 4.81%, 1.05%, 21.7%, 0.04%, 8.92%, and 7.26%, respectively, while downland regions were occupied by “Alluvial soils″, “Red-Yellow Podzolic soils with soft or hard laterite″, “Regosols on Recent beach and dune sands″, and “Latosols and Regosols on old red and yellow sands″, which covered the basin area by 3.12%, 16.50%, 0.11%, and 6.72%, respectively(

Table 9. Soil erodibility factors and area percentage of soil types in Maha Oya River Basin.

Soil Type	K Factor	Area (km2)	Area (%)
Rock knob plain	0.10	2.42	0.16
Red-Yellow Podzolic soils; steeply dissected, hilly and rolling terrain	0.22	450.19	29.61
Alluvial soils of variable drainage and texture; flat terrain	0.31	47.46	3.12
Red-Yellow Podzolic soils & Mountain Regosols; mountainous terrain	0.22	73.11	4.81
Steep rockland & Lithosols	0.25	16.03	1.05
Red-Yellow Podzolic soils with soft or hard laterite; rolling and undulating terrain	0.22	250.90	16.50
Regosols on Recent beach and dune sands; flat terrain	0.48	1.60	0.11
Red-Yellow Podzolic soils & Low Humic Gley soils; rolling and undulating terrain	0.16	329.92	21.70
Latosols and Regosols on old red and yellow sands; flat terrain	0.41	102.33	6.72
Erosional remnants (Inselbergs)	0.10	0.56	0.04
Reddish Brown Latosolic soils; steeply dissected, hilly and rolling terrain	0.17	135.66	8.92
Immature Brown Loams; steeply dissected, hilly and rolling terrain	0.33	110.43	7.26
Total		1521	100

). The K values ranged between 0.1 and 0.48 t ha h MJ⁻¹ ha⁻¹ mm⁻¹ (Figure 9).

3.4.3. Slope Length and Steepness (LS) Factor

The LS factor map was generated using slope length and steepness and is shown in Figure 9. The generated L factor values ranged between 0 and 6.9 while S factor values ranged from 0.03 and 15.03. Therefore, the LS factor values vary from 0 to 23.43.

3.4.4. Cover Management (C) Factor

The values for C and P factors were obtained through the literature (Table 4). The spatially distributed C factor values ranged between 0.1 and 1 while P factor values ranged from 0 to 1. All the C factor and P factor maps that were generated in this study are shown in Appendix C.

3.5. Soil Erosion Estimation

The average annual soil loss rates (t ha⁻¹ yr⁻¹) for the considered years in this study were estimated by the Revised Universal Soil Loss Equation (RUSLE). Figure 10 shows estimated soil loss rates (t ha⁻¹ yr⁻¹) for the years 1989, 2000, 2009, and 2021.

According to Figure 10, the average annual soil loss rate of the MORB ranged between 0 to 612.96 t ha⁻¹ yr⁻¹. In 1989, the soil loss rate varies from 0 to 398.644 t ha⁻¹ yr⁻¹ with a mean erosion rate of 2.81 t ha⁻¹ yr⁻¹, while years 2000, 2009, and 2021 had the mean average annual soil erosion rates of 2.76 t ha⁻¹ yr⁻¹, 2.65 t ha⁻¹ yr⁻¹, and 3.21 t ha⁻¹ yr⁻¹, respectively. Thus, the soil erosion rate considerably increased from the year 1989 to 2021. The evaluated soil erosion rates were classified into five severity classes to identify the erosional hazard areas between 1989 and 2021 in the MORB [18]. The obtained results are denoted in

Table 10. Total soil loss and corresponding erosion hazard classes in MORB from 1989 to 2021.

Erosion Hazard Classes	Soil Erosion (t ha−1 yr−1)	Total Annual Soil Loss
		1989		2021
		(ton)	(%)	(ton)	(%)
Low	0–5	131,430.2	30.73	126,509.0	25.89
Moderate	5–12	161,406.2	37.74	158,975.3	32.54
High	12–25	96,165.2	22.48	118,956.6	24.35
Very High	25–60	32,812.8	7.67	65,162.2	13.34
Extremely High	>60	5895.3	1.38	18,950.7	3.88
Total		427,709.7	100	488,553.8	100

. The results revealed that the mean annual soil loss of the MORB which occurred due to rill and sheet erosion in 1989 ranged between 0–81.81 t ha⁻¹ yr⁻¹, while it was 0–91.01 t ha⁻¹ yr⁻¹ in 2021.

As shown in

Table 11. Distribution of mean annual soil erosion under different erosion hazard classes in MORB from 1989 to 2021.

Erosion Hazard Classes	Soil Erosion (t ha−1 yr−1)	1989			2021			Net Change (%)
		Area		Mean Soil Erosion (t ha−1 yr−1)	Area		Mean Soil Erosion (t ha−1 yr−1)
		(km2)	(%)		(km2)	(%)
Low	0–5	1236.92	81.34	1.06	1217.73	80.08	1.04	−1.26
Moderate	5–12	213.59	14.05	7.56	209.99	13.81	7.57	−0.24
High	12–25	59.72	3.93	16.10	71.85	4.73	16.56	0.8
Very High	25–60	9.76	0.63	33.63	18.90	1.24	34.48	0.61
Extremely High	>60	0.72	0.05	81.81	2.08	0.14	91.01	0.09
Total		1521	100		1521	100

, the evaluated annual soil loss rates were categorized into five erosion hazard classes: low (0–5 t ha⁻¹ yr⁻¹), moderate (5–12 t ha⁻¹ yr⁻¹), high (12–25 t ha⁻¹ yr⁻¹), very high (25–60 t ha⁻¹ yr⁻¹), and extremely high (>60 t ha⁻¹ yr⁻¹) [18]. The spatial temporal changes between 1989 and 2021 are shown in Appendix E. The estimated total annual soil losses (tons) in the MORB for 1989 and 2021 are provided in Table 11.

According to Table 11, the extremely high erosion hazard class occupied 0.72 km² (0.05%) of the study area with a mean average annual soil loss rate of 81.81 t ha−1 yr−1, and 2.08 km² (0.14%) with a mean average annual soil loss rate of 91.01 t ha⁻¹ yr⁻¹ in the years 1989 and 2021, respectively. The low-hazard-class-occupied area in 1989 and 2021 were 1236.92 km² (81.34%) with a mean soil loss of 1.06 t ha⁻¹ yr⁻¹ and 1217.73 km² (80.08%) with a mean soil loss of 1.04 t ha⁻¹ yr⁻¹, respectively. In 1989, the erosion hazard classes of moderate erosion (7.56 t ha⁻¹ yr⁻¹), high erosion (16.10 t ha⁻¹ yr⁻¹), and very high erosion (33.63 t ha⁻¹ yr⁻¹) occupied areas of 213.59 km² (14.05%), 59.72 km² (3.93%), and 9.76 km² (0.63%), respectively. In 2021, moderate erosion (7.57 t ha⁻¹ yr⁻¹), high erosion (16.56 t ha⁻¹ yr⁻¹), and very high erosion (34.48 t ha⁻¹ yr⁻¹) severity classes occupied an area of 209.99 km² (13.81%), 71.85 km² (4.73%), and 18.90 km² (1.24%), respectively. The erosion hazard maps of 1989, 2000, 2009, and 2021 are shown in Figure 11.

3.6. Influence of Slope and Elevation on Soil Erosion

The slope (degree) was categorized into six classes (

Table 12. Soil loss rates in slope classes and net changes between 1989 and 2021.

No.	Slope Class (Degree)	Area		Soil Loss (t ha−1 yr−1)		Net Change (t ha−1 yr−1)
		(km2)	(%)	1989	2021
1	0–2	290.62	19.11	0.78	1.02	0.23
2	2–5	363.03	23.87	1.55	1.90	0.36
3	5–10	331.24	21.78	2.51	2.98	0.47
4	10–15	205.61	13.52	3.81	4.07	0.26
5	15–30	290.22	19.08	5.47	5.88	0.41
6	>30	39.98	2.63	7.18	9.21	2.03

) based on FAO slope categorization according to Nut et al. (2021), and the obtained results are denoted in Table 12. The results revealed that the rate of soil erosion rapidly increases with the increase in slope. The minimum rate of 0.78 t ha⁻¹ yr⁻¹ in 1989 and 1.02 in 2021 was occupied by the slope class of less than 2° (290.62 km²) while the maximum erosion rate of 7.18 t ha⁻¹ yr⁻¹ in 1989 and 9.21 t ha⁻¹ yr⁻¹ in 2021 was occupied by the slope category that is greater than 30° (39.98 km²).

The elevation of the MORB was categorized into five different altitude classes as suggested by Nut et al. [11], and we estimated the relevant soil loss rates as shown in

Table 13. Soil loss rates in elevation classes and net changes between 1989 and 2021.

No.	Elevation (m)	Area		Soil Loss (t ha−1 yr−1)		Net Change (t ha−1 yr−1)
		(km2)	(%)	1989	2021
1	0–300	1315.12	86.48	2.57	2.80	0.26
2	300–600	128.90	8.48	4.99	5.25	2.85
3	600–900	56.97	3.75	3.24	6.09	5.74
4	900–1200	17.44	1.15	3.41	9.14	7.19
5	1200–1500	2.27	0.15	3.20	10.39	0.23

. The elevation class of 0–300 m (1315.12 km²) had a soil loss rate of 2.57 t ha⁻¹ yr⁻¹ and 2.80 t ha⁻¹ yr⁻¹ in 1989 and 2021, respectively. The 300–600 m (128.90 km²) elevation class had a soil loss rate of 4.99 t ha⁻¹ yr⁻¹ in 1989 and 5.25 t ha⁻¹ yr⁻¹ in 2021. The 1200–1500 m elevation had the maximum soil loss rate of 10.39 t ha⁻¹ yr⁻¹ in 2021 while the minimum erosion rate was occupied by the 0–300 m elevation class in both years. The maximum net erosion change (7.19 t ha⁻¹ yr⁻¹) occurred at the elevation class of 900–1200 m, between 1989 and 2021.

3.7. Impact of Land Use and Land Cover Changes on Soil Erosion

The results indicated in

Table 14. Annual soil loss for specific LU/LC types of MORB in years 1989 and 2021.

LU/LC Type	1989			2021
	Mean Soil Erosion (t ha−1 yr−1)	Annual Soil Loss (ton)	Annual Soil Loss (%)	Mean Soil Erosion (t ha−1 yr−1)	Annual Soil Loss (ton)	Annual Soil Loss (%)
Bare Land	5.93	1467.21	0.34	7.30	86,136.53	17.63
Agricultural Land	1.08	28,784.79	6.73	0.95	22,502.59	4.61
Settlement	0.42	2604.13	0.61	0.58	7619.96	1.56
Forest	2.66	152,618.82	35.69	3.26	163,204.40	33.41
Vegetation	4.21	239,739.56	56.06	4.15	200,284.89	41.00

reveal that the estimated total annual soil losses from specific LU/LC types in the MORB were approximately 425,214 tons in 1989 and 479,748 tons in 2021, with an increment of about 54,534 tons from 1989 to 2021. The soil loss of specific LU/LC types between 1989 and 2021 is denoted in Appendix F. According to Table 14, the highest mean annual soil loss occurred due to bare lands in both 1989 and 2021. Bare lands (lands that were prepared for cultivation purposes) dominated the mean annual soil loss rate of 5.93 t ha⁻¹ yr⁻¹ and 7.30 t ha⁻¹ yr⁻¹ in the years 1989 and 2021, respectively. Agricultural land and bare land acquired a total annual soil loss of 28,784.79 tons (6.73%) and 1467.21 tons (0.34%), respectively, in 1989, and soil loss from agricultural land declined to 22,502.59 tons (4.61%) in 2021 while bare land increased annual soil loss to 86,136.53 tons (17.63) in 2021. According to the obtained results, vegetation cover had the highest amount of soil loss, 239,739.56 tons (56.06%), in 1989 and it decreased to 200,284.89 tons (41%) in 2021. Forest cover had a soil loss of 152,618.82 tons (35.69%) in 1989 and it was increased to 163,204.40 tons (33.41%) in 2021. Compared to the other LU/LC types, minimum soil loss occurred in settlements (Figure 12); it was 2604.13 tons (0.61%) in 1989, and 7619.96 tons (1.56%) in 2021.

The spatially distributed LU/LC conversion from 1989 to 2021 is shown in Figure 13.

Table 15. Primary Conversions of LU/LC types in MORB from 1989 to 2021.

No.	Primary LU/LC Conversion	Max Soil Erosion Range (t ha−1 yr−1)	Mean Soil Erosion (t ha−1 yr−1)	Changed Area (km2)	Average Annual Soil Loss (ton)	Average Annual Soil Loss (%)
1	Agricultural Land—Bare Land	377.98	6.04	23.06	13,926.35	2.85
2	Agricultural Land—Forest	329.52	2.49	49.02	12,189.32	2.50
3	Agricultural Land—Settlement	121.01	0.71	13.50	953.65	0.20
4	Agricultural Land—Vegetation	252.81	3.88	73.04	28,354.61	5.80
5	Forest—Agricultural Land	414.74	1.27	46.65	5942.78	1.22
6	Forest—Bare Land	612.96	8.43	41.92	35,327.34	7.23
7	Forest—Settlement	196.22	0.79	47.23	3734.07	0.76
8	Forest—Vegetation	427.39	5.25	148.96	78,147.25	16.00
9	Settlement—Bare Land	448.78	6.85	8.19	5608.33	1.15
10	Settlement—Vegetation	256.65	4.56	14.98	6832.21	1.40
11	Vegetation—Agricultural Land	183.88	1.06	76.04	8073.42	1.65
12	Vegetation—Bare Land	210.80	5.82	38.05	22,148.77	4.53
13	Vegetation—Forest	431.03	2.72	171.81	46,760.24	9.57
14	Vegetation—Settlement	118.04	0.61	43.23	2630.94	0.54
15	Water Bodies—Bare Land	407.87	5.31	3.30	1750.52	0.36

indicates only significant LU/LC changes and their contribution to soil erosion in the MORB from 1989 to 2021. The maximum soil erosion rate of 8.43 t ha⁻¹ yr⁻¹ occurred when forest cover converted into bare land and the average annual soil loss of forest to bare land was 35327.34 tons (7.23%), while the minimum mean soil erosion of 0.61 t ha⁻¹ yr⁻¹ occurred when vegetation converted into settlements and 2630.94 tons (0.54%) of average annual soil loss was estimated. The highest soil erosion rates of 6.85 t ha⁻¹ yr⁻¹, 6.04 t ha⁻¹ yr⁻¹, 5.82 t ha⁻¹ yr⁻¹, and 5.31 t ha⁻¹ yr⁻¹ occurred when settlement, agricultural land, vegetation, and water bodies converted into bare land, respectively, and their average annual soil loss was 5608.33 tons (1.15%), 13926.35 tons (2.85%), 22148.77 tons (4.53%), and 1750.52 tons (0.36%), respectively. The highest amount of soil loss, 78147.25 tons (16%), occurred when forest cover was converted into vegetation with a 5.25 t ha⁻¹ yr⁻¹ mean erosion rate.

3.8. Extremely High Soil Erosion-Prone Locations Map

The extremely high soil erosion-prone locations of the MORB were evaluated using soil erosion-hazard maps and Google Earth Pro. The locations with average annual soil loss greater than 60 t ha⁻¹ yr⁻¹ were selected to develop the maps. The generated extremely high soil erosion-prone locations map of the MORB is given below (Figure 14), and the corresponding location coordinates are given in Appendix A.

4. Discussion

This study used the Revised Universal Soil Loss Equation (RUSLE) to estimate spatially distributed average annual soil loss in the MORB from 1989 to 2021 using remote sensing data and a GIS environment. No studies have been conducted to evaluate average annual soil loss and the effects of various ecological types on soil erosion in the MORB. In Sri Lanka, a few studies have been conducted to estimate soil loss rates in some significant basins and the literature was used to obtain the support practice factor [12,18], cover management factor [60], and soil erodibility factor [12,18] values relevant to the Sri Lankan conditions. Piyathilake et al. [20] predicted that soil erosion can be mostly influenced by the R factor. Rainfall erosivity mainly depends on the mean annual precipitation and it can particularly affect soil erosion [45]. Therefore, the rainfall erosivity factor was determined in this study using the equation that was originally developed for the Sri Lankan ecological conditions by Premalal et al. [62]. The estimated spatially distributed average annual rainfall map in the MORB is shown in Figure 8. Then, the spatiotemporal variation of R factor is predicted in Figure 8. Many researchers in Sri Lanka as well as other countries followed the same procedure (IDW interpolation method) that was used in this study to evaluate the rainfall erosivity factor [4,12,16,44]. According to the obtained rainfall data, the highest average annual precipitation was observed from the Ambepussa rain gauging station situated in the central region of the MORB (Figure 8). Thus, the evaluated highest value of R factor was obtained from Ambepussa, Alawwa, Thulhiriya, Warakapola, Mirigama, and the central region of the study area (Figure 10). The minimum average annual precipitation was obtained from the western region of the study area around Dankotuwa and Wennappuwa, while the minimum R values were obtained from the same region. The estimated R values of this study ranged between 161.315 MJ mm ha⁻¹ h⁻¹ yr⁻¹ to 217.923 MJ mm ha⁻¹ h⁻¹ yr⁻¹ and R values of the Kelani River basin ranged from 235 MJ mm ha⁻¹ h⁻¹ yr⁻¹ to 431 MJ mm ha⁻¹ h⁻¹ yr⁻¹, due to average annual rainfall that varies from 2220 mm to 4236 mm, while the MORB had a variation from 1651 mm to 2170 mm. As a result, the excessive R factor values for the study area can be obtained in the rainy season, from April/May to September, known as “Yala Season”, which occurs mainly due to the southwest monsoon for the western region and the hill country of Sri Lanka.

The K factor map of this study (Figure 9) revealed that the minimum erodibility values were found in the upper areas in the MORB. The soils with low K factor values are more resistant to soil particle separation (e.g., clay soil) while high erodibility values indicate soil with more silt and sand that can easily be susceptible to erosion due to heavy rains [20]. In this study, the erodibility values were obtained from previous literature, due to lack of data for the calculations. As predicted by Piyathilake et al. [20], K values are not changed due to the climate changes in Sri Lanka and generally remain constant within the country.

The estimated average annual soil loss for the MORB varies from zero to 612.963 t ha⁻¹ yr⁻¹, with the mean value of average annual soil erosion of 3.21 t ha⁻¹ yr⁻¹ (Figure 10), while the recent study of the Kelani River basin that is located near the MORB had average an annual soil erosion between zero and 103.7 t ha⁻¹ yr⁻¹ with a mean annual erosion rate of 10.9 t ha⁻¹ yr⁻¹ [18]. Besides that, the most recent study in the Nalanda Oya catchment, Sri Lanka revealed that the mean soil erosion of the catchment was 2.99 t ha⁻¹ yr⁻¹ [13,19]. It was evaluated that the Kalu Ganga River basin had an annual soil loss rate between zero and 134 t ha⁻¹ yr⁻¹ with a mean annual erosion rate of 0.63 t ha⁻¹ yr⁻¹. The study regarding the Bibili Oya watershed revealed that, its average annual soil loss varies from zero to 435.42 t ha⁻¹ yr⁻¹ with a mean annual soil loss of 12.96 t ha⁻¹ yr⁻¹ [20]. According to these findings, average annual soil loss can be varied from one location to another due to the LS factor and the received average annual rainfall [11].

In this study, the derived LS factor range between zero to 23.4 (Figure 9), while the Kelani, Kalu, and Bibili river basins had LS factors between zero to 111, zero to 56.7, and zero to 64.7, respectively. The LS factor can be varied due to variations in topography [11]. As the LS factor increases, the velocity of water flow rises and, consequently, the rate of erosion increases. In this way, soil loss increases proportionately with increasing slope length and steepness, and the combined effects of slope length and steepness (LS factor) lead to reasonable estimates of soil erosion rates [20,45]. The C and P factors also severely affect the soil erosion. As revealed by the results, the forest cover of 566.7 km² (27.27%) in 1989 reduced to 500.7 km² (32.93%) in 2021, and vegetation cover of 567.9 km² (37.34%) in 1989 reduced to 471.7 km², while bare lands (lands that prepared for the cultivation) increased by 117.6 km² (7.73%) from 1989 to 2021, and total cultivation lands of the MORB were increased from 271.3 km² (in 1989) to 360.6 km² (in 2021), which can be a reason for the increment of the soil erosion rate between 1989 and 2021, according to Nut et al. [11]. This study mainly focused on sheet and rill erosion that happened due to rainfall and runoff in the MORB, due to the incapability of the RUSLE model to estimate the soil loss rates of gully erosion and dispersive soils [45,48,63].

The results revealed that most of the study area (80.08%) was dominated by the low erosion hazard class (0–5 t ha⁻¹ yr⁻¹) in the considered years, due to vegetation (31.02%) and forest cover (32.93%) in the MORB, because vegetation cover and plants/residues can reduce the soil erosion potential while decelerating the surface runoff and infiltrating the excess water to protect the soil from effects of rainfall and splashes [6]. However, the low erosion hazard areas cannot be neglected, because 126508.96 tons (25.89%) of the gross soil loss estimated in 2021 were from low hazard areas (Table 11). Therefore, proper soil conservation practices need to be implemented to maintain the soil characteristics and structural stability.

The highest average annual soil loss rate of 7.30 t ha⁻¹ yr⁻¹ was obtained from bare lands (lands that prepared for the cultivation purposes) in 2021 and bare lands were increased by 117.6 km² from 1989 to 2021 while forest and vegetation cover significantly decreased (Table 14). The increment of cultivation lands can generally lead to extreme soil erosion [11]. The water bodies, settlements, and exposed rocks had the least mean annual soil loss rate in 1989 and 2021 while bare lands, agricultural land, forest, and vegetation predicted slightly high erosion rates. The cause for forest and vegetation having somewhat high mean annual soil erosion rates can be the topography of the MORB, since the slope of the study area is generally high (Figure 4) and more than 35% of the study region had strongly high slopes (Table 12), while most of the forest and vegetation cover spread in the same sloping regions (Figure 4). As Moore & Wilson [56] posited in their study, the slope considerably affects the soil loss rates. Nut et al. [11] assert that steep terrain has high rates of soil erosion. Therefore, vegetation and forest cover in this study demonstrated a significant rate of soil erosion.

The relationship between soil erosion and LU/LC changes was evaluated using LU/LC maps and average annual soil loss maps in 1989 and 2021 (Table 14). Monitoring patterns of LU/LC change and soil erosion risk can be accomplished with help from this relationship [10]. Hence, the results divulged that decreasing forest and vegetation cover with an increment of cultivation lands can lead to the proliferation of soil erosion. Since the lands prepared for cultivation were classified as bare land and that land cover type had the highest mean annual soil erosion rate among other LU/LC types (Table 14), it can be concluded that human activities affect changes in land cover, and those changes may have an impact on soil erosion. However, there are several factors that contribute to soil erosion besides human activity, such as a steep slope, heavy rainfall, particular soil types, and strong winds [6].

Furthermore, Langbein in 1958 [64] posited in their study that dry or wetter weather can reduce soil erosion, because during dry climates, the runoff will decrease, whereas in wet climates, vegetation density will increase. The high bulk density of vegetation can decrease soil erosion [64]. Besides that, the redistribution of sediments (microtopographic features) due to banned vegetation patterns can reduce soil erosion by decreasing the soil erodibility in vegetation cover [65]. As Srivastava in 2020 [66] predicted in their study, the agricultural land cover can be observed using the Variable Infiltration Capacity (VIC) approach and it can be used to evaluate the crop, environment, and water resource management strategies.

As per the evaluated results, vegetation, forest, bare lands, agricultural lands, and settlements were the erosion-prone LU/LC types in the MORB. The central, eastern, southern, and southeastern parts of the study region have been identified as erosion-prone areas (Figure 14). The extremely high erosion-hazard locations map was developed in this study using the erosion hazard classes map of 2021 (Figure 14) and Google Earth Pro software. As for the results, the most extreme erosional locations were detected from steep slope zones and near the cultivation lands. Cultivation lands (agricultural lands and bare lands) in the study area will be adversely affected by soil erosion, if forest lands are transformed into agricultural lands, which would increase the risk of soil erosion. Therefore, deforestation affected by agriculture must be diminished, particularly on steep slopes. The study done by Babel et al. [67] posited that terracing, contouring, reforestation, grassed waterways, and filter strips can mitigate the sediment yield in a watershed. Hence, to reduce soil erosion sustainably and increase agricultural yields, this study provides several recommendations: prevention of deforestation and implementation of afforestation practices; implement methods to conserve the soil; practice the contour farming method when farming on hillslope lands, considering the land′s natural contours; river sand mining should be regulated, and a healthy ecosystem should be maintained; management and conservation of water; and construction of embankments along erosion-prone areas.

Currently, one of the main focuses of the civil engineering sector is on the United Nations Sustainable Development Goals (UN-SDGs) to make the world a better place with sustainability for all life before 2030 [68]. According to SDG 2, the world strives to end hunger by ensuring food security and increasing agricultural productivity. As mentioned in the introduction to this study, soil erosion is a significant cause of reduced agricultural production and land degradation. In addition, the results of the study revealed that agricultural LU/LC types (bare land) had the highest soil loss rate (Table 14). Therefore, to reduce soil loss in agricultural land, this study recommends terrace cultivation in hilly-mountainous areas. Apart from that, hilly cultivated lands require full crop cover using intercropping methods. Furthermore, planting deep and shallow rooted crops alternately improves soil structure and minimizes soil erosion. SDG 6 talks about sustainable management of water and sanitation. One way of achieving sustainable water management is by reducing the significant delivery of sediment to surface water bodies due to sheet and rill erosion. For the most part, implementing pasture on the soil surface can reduce the impact of raindrops and control surface runoff. Therefore, proper soil and water conservation methods can improve water quality in eroding areas. According to SGD 13, the effects of climate changes on soil erosion must be mitigated. As mentioned in this study, rainfall is a major factor that can greatly increase soil erosion. Therefore, this study recommends reforestation and preventing deforestation within the erosion-susceptible regions and implementing proper soil conservation practices to reduce the impact of climate change on erosion. According to SDG 15, every life on earth should have sustainably managed lands in the future for living without any disasters. The results of this study revealed that in 2021, the total soil loss of the MORB is about 479,748 tons (Table 14), which is a significant amount of soil loss. Hence, erosion is one of the primary causes of land degradation. To prevent soil loss and rehabilitate the impacted lands, erosion-preventative measures must be adopted as mentioned in this study.

Finally, the accuracy of the evaluated soil erosion in the MORB can be validated by conducting field experiments, and future analyses of soil erosion in the study region should consider gully erosion and dispersion soils because the RUSLE model is only capable of assessing rill and sheet erosion. The high erosion-prone areas identified in this study should be protected using appropriate conservation measures. The generated maps using parameters of the RUSLE model in this study will be beneficial for addressing the erosion across the MORB and making appropriate conservation practices and investments.

5. Conclusions

The Revised Universal Soil Loss Equation (RUSLE) was utilized in this research to assess the spatially distributed average annual soil loss in the Maha Oya River Basin (MORB) from 1989 to 2021 using remotely sensed data with a GIS environment, which allowed us to evaluate the areas that were severely susceptible to soil erosion and to recommend appropriate conservation practices for mitigating the erosion sustainably. The findings of the study were used to develop maps of high erosion-prone locations in the MORB, making it possible to pinpoint high erosion-prone areas. As for the results of the study, there was a sharp increment in cultivated lands and settlements while forest and vegetation had negative growth due to environmental devastation within the study period. The study results revealed that mean annual soil loss rates of 1989, 2000, 2009, and 2021 in the MORB were 2.81 t ha⁻¹ yr⁻¹, 2.76 t ha⁻¹ yr⁻¹, 2.65 t ha⁻¹ yr⁻¹, and 3.21 t ha⁻¹ yr⁻¹, respectively. It was estimated that 6.11% of the study area would be classified in the high to extremely high erosion severity classes (>12 t ha⁻¹ yr⁻¹) in 2021, whereas 4.61% of the study area would be in these severity classes in 1989. The upland regions of the study area had the maximum rates of mean soil erosion ranging from 16.56 t ha⁻¹ yr⁻¹ to 91.01 t ha⁻¹ yr⁻¹, due to the steep slopes, vegetation/forest cover deprivation, and growth of cultivation lands. This study observed that agricultural lands, bare lands (lands that were prepared for cultivation), and lands that lay on the high steep slope regions had considerably high average annual soil loss rates. The results highlighted that the severity of soil degradation increases with the slope of the terrain. According to the LU/LC change analysis, cultivated lands was identified as the most vulnerable LU/LC type of the MORB.

The current research can also be used as an initial study in the process of achieving some of the SDGs related to food security, sustainable water and sanitation, actions on climate change impacts, and finally, achieving sustainable ecosystems (SDG 2, 6, 13, 15). When adopting proper soil conservation measures and making the necessary investments, the maps produced by this study will be helpful for addressing the erosion across the MORB. In addition, the results of this study will be beneficial for decision makers and future researchers. However, the current study is only an exploratory survey on the soil erosion rates and the spatial variation of erosion hazard levels in the Ma Oya River basin. We recommend the validation of the results of the current study by using field test data (that are not currently available for the study area) in future research to confirm the RUSLE model results obtained in our study.