Soil Erosion Susceptibility Prediction in Railway Corridors Using RUSLE, Soil Degradation Index and the New Normalized Difference Railway Erosivity Index (NDReLI)

: This study presents a remote sensing-based index for the prediction of soil erosion susceptibility within railway corridors. The empirically derived index, Normalized Difference Railway Erosivity Index (NDReLI), is based on the Landsat-8 SWIR spectral reﬂectances and takes into account the bare soil and vegetation reﬂectances especially in semi-arid environments. For the case study of the Botswana Railway Corridor (BRC), the NDReLI results are compared with the RUSLE and the Soil Degradation Index (SDI). The RUSLE model showed that within the BRC, the mean annual soil loss index was at 0.139 ton ha − 1 year − 1 , and only about 1% of the corridor area is susceptible to high (1.423–3.053 ton ha − 1 year − 1 ) and very high (3.053–5.854 ton ha − 1 year − 1 ) soil loss, while SDI estimated 19.4% of the railway corridor as vulnerable to soil degradation. NDReLI results based on SWIR1 (1.57–1.65 µ m) predicted the most vulnerable areas, with a very high erosivity index (0.36–0.95), while SWIR2 (2.11–2.29 µ m) predicted the same regions at a high erosivity index (0.13–0.36). From empirical validation using previous soil erosion events within the BRC, the proposed NDReLI performed better that the RUSLE and SDI models in the prediction of the spatial locations and extents of susceptibility to soil erosion within the BRC.


Introduction
Soil erosion refers to the physical degradation of topsoil through the detachment and removal by flowing waters resulting from rainfall and runoff or by blowing winds. The process of water-or wind-induced soil erosion occurs naturally and is accelerated by anthropogenic activities such as deforestation, urbanization, agricultural activities and climate change [1,2]. The resulting movement and deposition of the soil and rock particles contributes to overland flow as well as surface runoff, which can affect infrastructural assets and the environment [3]. The factors and mechanisms that influence the process of water-or wind-induced erosion such as ground cover, topographic structure, rainfall intensity and soil type are geographically dynamic and vary in spatiotemporal scales [4][5][6][7].
Globally, it has been estimated that the total land area damaged by soil erosion is about 1.1 billion hectares, which results in the annual transportation of approximately 2-2.5 × 10 10 Mg of soils [8]. In Africa, approximately half of the 494 million hectares of land that are subjected to different types of degradations are influenced by water [9]. Further, using the RUSLE model, [10] estimated that approximately 6.1% of the global land area experienced very high soil erosion rates, which exceeded the 10 Mg ha − 1 yr −1 estimated in 2012 [11]. To this extent and to assess soil erosion by water and predict its intensity and risks, several approaches and models have been developed. The models Remote Sens. 2022, 14, 348 2 of 27 obviously have different demands for data and complexity in terms of the processes and calibration requirements [12,13].
Models such as the Universal Soil Loss Equation (USLE) by [14] Wischmeier and Smithy (1978), the Revised USLE (RUSLE) by [15] Renard et al. (1997), the Water Erosion Prediction Project (WEPP) [16] (Laflen et al. 1991) and the Soil and Water Assessment Tool (SWAT) [17] (Arnold et al. 1998) have popularly been used in the estimation of soil erosion. Notably, the choice of the most suitable model is a logical process influenced by several factors including land-use, the characteristics of the catchment and, most importantly, the data available [2,13].
Physical and conceptual models such as the European Soil Erosion Model (EUROSEM) and WEPP are considered superior to the RUSLE since they incorporate runoff in determining the erosion processes. The models are, however, complex and require substantial data, which usually preclude their broad application in most areas. The empirical models on the other hand require specific parametrizations. For example, the Soil Loss Estimation Model for Southern Africa (SLEMSA) is site specific, while the Modified USLE (MUSLE) requires monitored runoff and peak flow inputs [18]. Comparatively, the RUSLE model is considered to be simple, flexible and time-and cost-effective. Though limited in terms of reliability and spatial coverage for large areas, the RUSLE model has the advantage of being integratable with GIS and remote sensing data [8,19].
Most soil erosion models have been used extensively in the estimation of soil loss, mostly at catchment scales [19][20][21][22][23]. However, no investigations have focused on the prediction of soil erosion and vulnerability mapping within road and railway corridors. The road and rail network corridors are critical regions within the national planning and economic systems, as they act as the engines for productivity and improved quality of life.
For a given region, the causes of landscape degradation that triggers soil erosion can be categorized as on-site and off-site impacts [24]. Off-site impacts are characterized by widespread damages such as the destruction of roads and railway lines, clogging of drainage systems and siltation of water reservoirs [7,24]. The off-site impacts on road and railways are also attributed to the long-term effects of anthropogenic activities that affect river systems by altering the channel flows. For road and railway infrastructural assets, the steep and lengthy slopes have additional contributions to sediment loss as the embankments accentuate the potential for soil erosion [25].
For sustainability, road and rail transport network operators should continuously monitor the vulnerabilities of the network corridors with the aim of understanding and mitigating the conditions contributing to potential infrastructural failures. However, no robust models have been developed to monitor the vulnerability of rail and road corridors to soil erosion. Reference [26] proposed the Normalized Difference Road Landslide Index (NDRLI) to recognize and classify road-related landslides. The method takes advantage of the unique spectral responses of bare soil, land covers and topographic characteristics. However, the NDRLI approach is unsuitable for flat areas without potentials for landslides and works under the assumption that farmland and bare soils have the same and higher reflectance than actual landslides. Since most road-induced landslides normally comprise of rocks, soils and vegetation, this assumption may not be applicable especially in flat and semi-arid areas with clay-sandy soils. Thus, the NDRLI may not be able for effectively separate barren and bare soils from the other soil erosion objects or features. Further, while the widely used RUSLE model is suitable for the quantification of the relative soil loss, it does not directly predict the spatial location and extent of soil loss. As such, the use of spectral index-based methods have been proposed [27,28].
The objective of this study is to develop and validate a remote sensing-based spectral index for the prediction and mapping of soil erosion vulnerability within railway corridors in semi-arid environments. The proposed approach is termed as the Normalized Difference Railway Erosivity Index (NDReLI). The NDReLI is a remote sensing reflectance index derived at medium spatial resolution data and is modelled based on: (1) the maximization of the reflectance of railway-induced landslides using the SWIR1 (1.57-1.65 µm) and SWIR2 (2.11-2.29 µm) spectral wavelengths; (2) minimization of the reflectance by erosion causing waters and (3) taking advantage of the high spectral reflectance of bare soil in the SWIR bandwidths. The rationale in using the SWIR wavelengths to detect soil erosivity during rainy events is in the ability of the SWIR signals to penetrate thin clouds and discriminate between vegetation and soil in terms of their water content and the resulting reflectances [29,30]. The novelty of the proposed index is in the ability to automate the process of predicting the spatial location and extent of soil erosivity within railway corridors. In addition, the approach does not require data related to climate and hydrology. The NDReLI results are compared with the RUSLE soil loss model and the Soil Degradation Index (SDI) for the case study of the Botswana Railway Corridor (BRC).

Study Area
The BRC is situated within Botswana's Limpopo River Basin (BLRB). The railway line is about 886 km in length running from Lobatse and Gaborone City in the south to Francistown in the north ( Figure 1). In the recent past, there have been reported incidents of derailments of trains along the railway line. The main cause has been attributed to soil erosion within the railway line corridor. This is based on the fact that the basin is semi-arid, and the vast bare-soils along the railway corridor are susceptible to flooding and subsequent soil erosion during rainy seasons. This coupled with the flooding of the rivers across the railway line increases the flooding threats during rainy seasons. In this study, a 200 m buffer corridor was delineated to represent the case study of the BRC. The 200 m buffer was considered optimal since it included the width of the railway reserve and captured the contributions and impacts of the areas around the railway line in terms of the terrain features and the influence of elevation and slope.

Data
The assessment of soil erosion within a catchment area is dependent on the regional climate, soil type, land-use and the topographic characteristics of the area under study. Table 1 presents a summary of the data used in this study. The study data were acquired in 2019 and harmonized to 12.5 m × 12.5 m grid spatial resolution using nearest-neighbor resampling. The Landsat images were calibrated for surface reflectance using the FLAASH (Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes) atmospheric correction model, which incorporates the MODTRAN (Moderate Resolution Transmittance) radiation transfer code to remove the effects of atmospheric gases and aerosols to produce a surface reflectance image. For maximum surface reflectance, the Landsat digital numbers (DNs) were converted to the Top Of Atmosphere (TOA) reflectance using the respective Landsat radiometric rescaling coefficients [31].

Revised Universal Soil Loss Equation (RUSLE)
The RUSLE model formulated by [32] represents the influence of topography, soils, climate and land use on rill and inter-rill soil erosions as triggered by the impacting raindrops and the resulting surface runoff [15]. The RUSLE model for soil loss index (SLI) estimation is given as in Equation (1): where A is the estimated mean annual soil loss (ton h −1 y −1 ); R represents the rainfall erosivity factor (MJ mm ha −1 h −1 y −1 ); K is the soil erodibility factor (ton h MJ −1 mm −1 ); L is the slope length (m); S is the slope steepness (%); C is the crop/land-cover management factor and P represents the conservation practice factor, with C and P being dimensionless factors.

Rainfall Erosivity (R-Factor)
The RUSLE R-factor is traditionally calculated as the product of the intensity of 30-min rainfall and the kinetic energy of rainfall [14]. Due to the lack of continuous records on rainfall data and high-resolution pluviograph data, this approach is unsuitable in most parts of the world [33,34]. This drawback is overcome by establishing the relationship between mean monthly rainfall data and its erosivity [35]. Using the [14] empirical formulation, the mean monthly rainfall data for 2019 from 13 meteorological stations within the BLRB were used to calculate the R-factor (Equation (2)): where R represents the rainfall erosivity factor (MJ mm ha −1 h −1 yr −1 ); p i is the total monthly precipitation (mm) and p is the mean annual precipitation (mm).

Soil Erodibility (K-Factor)
Remote Sens. 2022, 14, 348 5 of 27 K refers to the effect of a given soil type on the erosivity impact of rainfall and runoff. The factor is influenced by the soil's organic matter content, texture, structure and the degree of permeability [36]. Soils that are highly susceptible to erosion have K-values close to 1, and K-values of close to 0 represent soils with a high resistance to erosion. In this study the K-factor was derived from the Digital Soil Map of the World (DSMW) at a scale of 1:5 million. For each soil type, the percentage content of clay, silt, and sand were used to estimate the value of the K-factor. Since the study comprises mostly of sandy soil, Equation (3) [37] was adopted: where K is the soil erodibility factor (thMJ −1 mm −1 ); f csand is a factor that lowers the K indicators in soils with high coarse-sand content and is higher for soils with less sand content; f cl−si gives low soil erodibility factors for soils with high clay-to-silt ratios and f orgc reduces K values in soils with high organic carbon content, while f hisand lowers K values for soils with extremely high sand content. The f-parameters are derived as in Equations (4)- (7): where R represents the rainfall erosivity factor (MJmm. ha −1 h −1 yr −1 ); pi is the total monthly precipitation (mm) and p is the mean annual precipitation (mm).

Soil Erodibility (K-Factor)
K refers to the effect of a given soil type on the erosivity impact of rainfall and runoff. The factor is influenced by the soil's organic matter content, texture, structure and the degree of permeability [36]. Soils that are highly susceptible to erosion have K-values close to 1, and K-values of close to 0 represent soils with a high resistance to erosion. In this study the K-factor was derived from the Digital Soil Map of the World (DSMW) at a scale of 1: 5 million. For each soil type, the percentage content of clay, silt, and sand were used to estimate the value of the K-factor. Since the study comprises mostly of sandy soil, Equation (3) [37]was adopted: where K is the soil erodibility factor (thMJ −1 mm −1 ); csand f is a factor that lowers the K indicators in soils with high coarse-sand content and is higher for soils with less sand content; fhisand ≟ 0.7 1 100 1 1 exp 5.51 22.9 1 100 100

Slope Length and Steepness (LS) Factor
The topographic effect on the degree of erosivity is determined by the influence of the length of the slope (L), the steepness of the slope (S) and the morphology of the slope. The soil loss per unit area increases with increase in slope length and as the runoff accumulates downslope. Increased slope steepness results in increased soil erosion as the velocity and erosivity of the accumulated runoff increases. The L-parameter defines the ratio of rill erosion flow to the impact of the raindrop resulting into inter-rill erosion for a plot length of 22.13 m. The S-factor represents the horizontal movement from the beginning of the overland flow to the deposition point [15], and the S-parameter defines the influence of the slope gradient on erosivity as compared to the plot with a slope steepness of 5.16°. The combined LS-factor is then defined as the ratio of soil loss per unit area from a field 22.1 m long, with a slope of 5.16°. The LS-factor is calculated from the 12.5 m × 12.5 m ALOS-PALSAR digital elevation model (DEM) as in Equations (8)

Slope Length and Steepness (LS) Factor
The topographic effect on the degree of erosivity is determined by the influence of the length of the slope (L), the steepness of the slope (S) and the morphology of the slope. The soil loss per unit area increases with increase in slope length and as the runoff accumulates downslope. Increased slope steepness results in increased soil erosion as the velocity and erosivity of the accumulated runoff increases. The L-parameter defines the ratio of rill erosion flow to the impact of the raindrop resulting into inter-rill erosion for a plot length of 22.13 m. The S-factor represents the horizontal movement from the beginning of the overland flow to the deposition point [15], and the S-parameter defines the influence of the slope gradient on erosivity as compared to the plot with a slope steepness of 5.16 • . The combined LS-factor is then defined as the ratio of soil loss per unit area from a field 22.1 m long, with a slope of 5.16 • . The LS-factor is calculated from the 12.5 m × 12.5 m ALOS-PALSAR digital elevation model (DEM) as in Equations (8) and (9): where LS is a dimensionless factor representing the slope length and slope steepness; A is the flow accumulation x cell-value of the upslope contributing area per unit cell (m); m (=0.4) is a variable slope length exponent and n (=1.3) is a slope steepness exponent [28].

Conservation Practice (P-Factor)
The land conservation practice factor is an indicator of the influence of the conservation and conservation practices on the impact of runoff on soil erosivity [15]. The factor represents the ratio of soil loss after intervention through a specific conservation practice to the determined soil loss after the up-down cultivations. The range for P-factor is [0, 1], with 0 indicating strong protective measures and 1 representing areas with no cultivation practice support [34]. Within the Botswana's LRB, there are seasonal agricultural activities due to the arid to semi-arid climatic conditions. However, since there are non evident conservation measures within the study area, the P-factor was assumed as equivalent to 1 as proposed by [14,38]. Furthermore, for the case study the slope factor in the support practice factor is insignificant since more than 90% of the corridor has slopes varying in the range of 0 • -5 • , which are classified as very low slopes.

5.
Land Cover and Management (C-Factor) The C-factor is considered as second most significant factor that influences soil erosion, after the topographic factors. C-factor is a measure of the runoff and soil erosion rate as controlled by the cropping and management practices within a region [14,15]. Since the kinetic energy of the raindrops are dissipated by the presence of vegetation on the soil surface, it implies that vegetation cover and cropping systems (plant residue and tillage) are the key indicators of the degree of the C-factor [39]. For the current application, the Cfactor was derived from the Normalized Difference Vegetation Index (NDVI), as calculated from Landsat-8 OLI from which the C-factor is derived (Equations (10) and (11)). NDVI determines the health and distribution of vegetation and hence takes into consideration the degree of ground cover by green vegetation.
where ρ i indicates the spectral band of Landsat-8 OLI in NIR and red wavelengths; C is the C-factor and α and β are parameters that determine the shape of the NDVI curve. The values α = 2 and β = 1 were adopted in this study, as they have been proven to produce reasonable results [28,40].

Soil Degradation Index
Being a semi-arid region, the change and degradation of the natural environment can be detected by changes in the color and mineralogy of the soil and also by the vegetation variations in structure and spatial distribution [41]. To determine the degree of degradation of soil within the railway corridor, the following soil reflectance indices were derived: NDVI, brightness index (BI), coloration index (CI) and form index (FI) (Equations (10) and (12)- (14)).
BI is a useful indicator in determining the soil degradation state, as it can be used to distinguish between plant cover and bare soils since low BI implies presence of vegetation cover while high BI is attributed to bare soil. CI, developed by [42], indicates that degraded Remote Sens. 2022, 14, 348 7 of 27 soils are generally brighter and clearer, as they contain low organic matter and hence low CI, and vice versa for less degraded soils. Additionally developed by [42], FI in Equation (14) characterizes the soil types and their levels of degradation using the reflectance curves derived from the visible bands, especially in arid and semi-arid regions. The NDVI is computed as in Equation (10) and is useful in soil degradation assessment as it maps the inverse correlations between the decrease in vegetation cover and increase in soil loss rate [43]. The spectral indices (BI, CI, FI, NDVI) are iteratively combined to determine the linear combination that best characterizes the soil degradation and bare soil occupation within the study area.

Normalized Difference Railway Erosion Index (NDReLI)
The spectral characteristics of soils depend on their composition and moisture content. Eroded soils tend to have a lower humus content, lighter shade and low coloration index. As shown in Figure 2a, bare soil absorbs most energy in the visible but has a higher spectral reflectance in the SWIR1 than in SWIR2. At the SWIR wavelengths, urban built-up areas have higher absorption in the SWIR1 and higher reflectance in the SWIR2 (Figure 2a,b. On the other hand, water reflects most of the energy from the visible to the infrared spectrum but attracts energy in the SWIR1 and SWIR2 wavelengths (Figure 2a). Most vegetation, as expected, absorbs the SWIR wavelengths, whereas the NIR energy is mostly reflected (Figure 2a). In order to detect bare soil, the proposed NDReLI spectral index utilizes the dissimilarities among bare soil, urban built-up and vegetation features in NIR and in the SWIR bands to isolate the bare soil areas that are the most susceptible to soil erosion. This enables NDReLI to enhance bare soils and spectrally isolate the similarly reflecting urban areas as depicted in Figure 2b. The advantage of using SWIR in detecting the erosion susceptibility during rainy seasons is the ability of the SWIR signals to penetrate thin clouds and discriminate vegetation and soil [29,30]. Particularly, SWIR2 band is sensitive to soil and vegetation water content [44], thus providing information on typical conditions that occur after the rainy events.
Conceptually, within the railway corridors the water-induced soil losses are enhanced by mapping the bare soils, while water is suppressed. Compared to the remote sensing index for road-induced landslides (NDRLI) by [26], which is based on the 20 m spatial resolution of the Sentinel-2 SWIR1 wavelength, NDReLI can be derived from both SWIR1 and SWIR2 spectral bands of the Landsat data (Equations (15)- (17)). The NDReLI index takes into account the existence of bare soil by using the Dry Bare-Soil Index (DBSI) [45]. The inclusion of DBSI improves the suitability of the NDReLI for the detection of bare soil areas in arid and semi-arid areas [46].
where ρ SW IRi , ρ blue , ρ green , ρ N IR and ρ red are, respectively, the spectral reflectance in the shortwave infrared, blue, green, near-infrared and red wavelengths from Landsat satellite data; DBSI is the dry bare-soil index and NDReLI index is in the range [−1, +1]. From its formulation in Equations (15)- (17), the NDReLI index considers the following factors in its modelling: (1) the maximization of the reflectance of railway-induced landslides using SWIR1 and SWIR2 spectral reflectances and the presence of bare soil; (2) minimization of the reflectance by erosion causing waters and (3) taking into consideration the high reflectance of the bare soil features in the SWIR spectral bandwidths to detect potential erosion events within railway corridors. For the case study, the bare soils include barren land with soil cover as well as fallow agricultural lands which are highly vulnerable to sediment transport and soil erosion from heavy rains.
While the widely used RUSLE model is suitable for the quantification of the relative soil loss, it does not directly predict the spatial location and extent of soil loss, and the use of remote sensing indices has been proposed [27,28]. NDReLI is considered a rapid approach for the effective prediction of areas with a high to very high susceptibility to soil erosions along railway corridors and is capable of distinguishing the noise from bare soils, water, vegetation and farmlands. Further, compared to NDRLI, NDReLI does not rely on trial-and-error and thresholding for the determination of the optimum results. The implementation steps for the derivation and evaluation of NDReLI is presented in Figure 3.

Soil Erosion Factors for RUSLE Model
For this study, the cell-by-cell calculation of the mean annual soil loss rate (ton ha −1 year −1 ) was derived from the RUSLE model integrated with GIS and remote sensing data to predict the areas susceptible to soil erosion along the Botswana railway corridor.

R-Factor
For the entire basin, the average annual rainfall varied from 25.83 mm to 1852.85 mm (Figure 4a), and the corresponding R-factor values obtained for the region ranged between 6.35 MJ mm ha −1 h −1 year −1 and 86.79 MJ mm ha −1 h −1 year −1 , as shown in Figure 4b. For the 200 m corridor (Figure 4c), the R-factor varied from 18.82 MJ mm ha −1 h −1 year −1 to 86.79 MJ mm ha −1 h −1 year −1 with the highest R-factor observed around the middle of the corridor. The rainfall intensities were spatially interpolated from 13 gauging stations using the inverse distance weighting (IDW) as applied in [47]. From the rainfall distribution within the corridor, there is observed high intensity especially in the lower, middle and upper sections of the corridor. These sections are expected to be at higher erosion risks due to the higher likelihood of runoff with high rainfall intensity and less vegetation cover.

K-Factor
The soil maps for the BLRB and BRC are, respectively, presented in Figure 5a,b, and the corridor has four main soil types. The derived soil erodibility factor from the soils is presented in Figure 5c, with the results showing the K-factor ranging from 0.068 to 0.150. The area has a mixture of clay and sandy coarse textured soils characterized by the K-values of about 0.071-0.122 and 0.138-0.150, respectively. Most of the soils within the corridor are not easily detached due to the predominantly low K-values of between 0.068 to 0.122. Few sections of the railway line are likely to be prone to soil erosion, as the K-values are slightly higher, at 0.138-0.150. Most of the areas with higher K-factor values coincidentally receive higher precipitation and could be more vulnerable to soil erosion due to the loose soils.

Slope Length (L) and Slope Steepness (S) Factor
The LS-factor is determined from the elevation, slope and the flow accumulation terrain parameters. While the elevation for the LRB ranges from 734 m to 1519 m above mean seal level (inset in Figure 5a), the 200 m corridor is elevated between 852 m and 1510 m (Figure 6a). The corresponding slope derived for the area is presented in Figure 6b and indicates that most of the corridor is gentle and flat at 0-5.8 • of slope. The LS-factor (Figure 6c) varied from 0 to 1.538, with a mean value of 0.044 over the corridor. The LS-factor was distributed as follows within the basin: 0-0.042 (65.6%); 0.042-0.139 (27%); 0.139-0.332 (6.1%); 0.332-0.845 (1.2%) and 0.845-1.538 (0.1%), implying that more than 90% of the corridor is relatively flat with a low C-factor. This implies that it would take time and a certain intensity of precipitation before runoff is generated.

C-Factor
The C-factor was computed using the NDVI, since the degree of soil erosivity is mostly influenced by the degree of vegetation cover. The NDVI was preferred to the land-use land-cover (Figure 7a), since the study area is semi-arid and sparsely settled with scarce vegetation cover. The NDVI within the corridor varied between −0.25 and +0.61 (Figure 7b) with most of the areas covered with bare soils and sparse vegetation. The corresponding C-factor values ranged from 0.043 to 1.488 (Figure 7c). From the results, the low C-values represents strong vegetation cover characterized by closed shrubs, which is predominant in the southern parts of the railway corridor and areas near Mahalapye in the middle of the corridor. The higher C-factor values imply scarce or no vegetation cover. The LULC map (Figure 7a) was used to verify the correspondence of the C-factor to the land-use type.

P-Factor
The P-factor for the study area, as discussed in Section 2.3.1. above, was assigned the value 1 since most of the land cover is covered with bare soil, with no erosion control practices and the seasonal agricultural practices mainly comprising of crop cultivation and pastoral farming.

Mean Annual Soil Loss Using RUSLE
As depicted in the schematic flow in Figure 3, the mean annual soil loss index within the BRC was determined through the integration of the five RUSLE factors with the results in Figure 8. It is observed that approximately 90% of the railway corridor is prone to very low and low soils rates with only 1% prone to high to very high rates of soil loss. The implication from the RUSLE results is that most of the railway corridor is not susceptible to significant soil erosion. Based on the results in Figure 8, the main areas affected by erosion-induced derailments are located in zones L2 and L3. These areas characteristically receive higher precipitation during the rainy seasons, which increases the localized R-factor. The increased R-factor, especially in areas A1, B1 and C1 (Figure 4), will trigger the high soil K-factor to detach, hence causing soil erosion. For areas with higher slopes in Figure 6, the LS-factor will further exacerbate the erosion events, and if the ground surfaces are not covered (low C-factor), the bare soils will be subjected to soil erosion mainly in areas L2 and L3 in Figure 8.
The RUSLE model was initially developed for small agricultural catchments in temperate regions characterized by gentle slopes and thick soil types with the objective of mapping the degrees of rill and sheet erosion [15]. Clearly, from the current findings, the RUSLE approach may not be relied upon for the detection of the degree of soil erosion at least within the BRC. This further motivates investigations into the use of remote sensing spectral reflectance for the delineation soil erosion susceptibility in such case studies and environments.

Soil Degradation Index (SDI) Using Spectral Indices
In the spectral analysis for SDI derivation, the BI, CI and FI indices were found to be the most significant. The respective results for the individual SDIs are presented in Figure 9a-c, and the combined SDI results are presented in Figure 9d. The BI map (Figure 9a) shows a good distinction between plant cover and bare soil, with built-up areas showing higher BI values as compared to areas with vegetation covers. The CI results in Figure 9b depicts characteristically low Cl values towards the southern regions of the railway line, from Mahalapye to Gaborone and to Lobatse. The low CI is attributed to degraded soils with low organic matter content. In correspondence to the CI, the FI results (Figure 9c) indicate that there is a higher degradation of soils in the southern part of the corridor. This is expected, as the southern sections of the railway line are more populated and developed. The northern sections of the corridor indicate higher-form index values due to more agricultural activities in comparison to building developments in the southern sections. The results indicate that most of the corridor is not highly degraded, as only 19.4% of the corridor is susceptible to soil degradation and hence has higher chances of being eroded depending on the amount and intensity of received rainfall. Compared to the RUSLE-SLI results in Figure 8, the SDI (Figure 9d) is able to predict the locations of the areas vulnerable to soil erosion within the BRC in the outlined regions L1, L2 and L3. The lower and middle sections of the corridor are observed to be susceptible to moderate and high erosivity, while the upper regions around Francistown are prone to moderate to very high erosivity according to the SDI index values.
The region bounded by L1 is slighted elevated, has a low K-factor, a low R-factor and low to very low slopes with moderate vegetation cover. As such, the L1 area is expected to have low soil degradation. Thus, as compared to the RUSLE results above, the SDI based on BI, FI and CI overestimated the soil loss especially in the region bounded by L1.

NDReLI
Using Landsat-8 OLI SWIR spectral wavelengths for the BRC case study, the NDReLI results represent the degree and extent of soil erosivity in terms of the spectral index for the SWIR1 (NDReLI1) and SWIR2 (NDReLI2). Figure 10a,b respectively shows the NDReLI1 and NDReLI2 soil erosivity results for the case study. For both NDReLI1 and NDReLI2, there are four regions within the corridor with the highest index values, namely the regions near Francistown (L1), Palapye Town (L2), Mahalapye/Dikgatlhong areas near Bonwapitse River (L3) and Lobatse Town (L4). Using the same index range for the two spectral indices, it is observed that while the same regions are predicted to be susceptible to soil erosion, NDReLI1 based on SWIR1 has higher index values with more areas under high vulnerability as compared to the NDReLI2 (SWIR2). This implies that the replacement of SWIR1 with SWIR2 reduces the influence of built-up and impervious surfaces that have similar spectral reflectance with bare soil.
The NDReLI results in Figure 10a,b are compared with the Normalized Difference Road Landslide Index (NDRLI) results in Figure 11. It is observed that from the NDRLI road index, the SWIR1 mapped most of the corridor as subject to moderate to very high susceptibility, while SWIR2 only marginally reduced the very-high and high vulnerabilities to moderate and high degrees of susceptibilities. The results show that for the study area, NDRLI is not able to automatically detect the low and very low susceptibility areas. The main difference between the two indices is that NDReLI considers and enhances the presence of bare and vulnerable soils using the DBSI index.
The NDReLI results shows that within the study area, not only is the spectral reflectivity of soil erosion varying to a high degree, it also tends to be similar to non-erosion features. Based on this observation, the susceptibility of the railway line to erosion is detected and determined based on the spatial distance of the predicted erosion event from the railway line and is considered positively if the event horizontally traverses the railway line. This ensures that the non-erosion features are not included as soil erosion features.

Empirical Validation and Comparisons of RUSLE, SDI and NDReLI
It is notable that the results of the compared models are predictive indicators and not the absolute rates of soil erosion. Since there were no measured in-situ soil loss estimates to validate the models, the accuracy of the susceptibility to erosion was empirically evaluated in terms of the absence and presence of soil erosion. The empirical evaluation was carried out using historical reports of the soil erosion events within the case study area and also tested for two railway corridors in Namibia.
For the BRC, the reported erosion-induced derailments in December 2019 and April 2020 near Bonwapitse River in the Mahalapye/Dikgatlong and Lobatse areas were used to compare the RUSLE soil loss index, SDI and NDReLI results. In both incidences, it was established that the causes of the derailments were attributed to the track structures being flooded by eroded soils. From the RUSLE soil loss prediction results (Figure 12a), it is observed that for both case study areas, the high to very high vulnerabilities are marginally detected with most of the areas within the corridor subjected to low to very low vulnerabilities.
The results imply that while the RUSLE soil loss index may estimate the relative amount of soil loss per year, it is not capable of accurately predicting the spatial locations and extents of the areas susceptible to erosivity especially in semi-arid and narrow environments. This observation is supported by [28] in their case study in Umzintlava catchment (T32E) in South Africa, where they concluded that the RUSLE method was not able to delineate all erosion in the study area and recommended the use of satellite-based estimates. This could be attributed to the fact that RUSLE and the related soil loss models are more suited for the modelling of the absolute soil loss values in areas affected by rill and inter-rill erosion and not for mapping and predicting the spatial extents and distributions of the soil erosion phenomena.  The results in Figure 12b show that for the case of the Lobatse area (L4), the SDI accurately showed that the line before and after the Lobatse town is susceptible to moderate to very high soil degradations. For the case of the Mahalapye area (L3), moderate to high degrees of potential degradations are predicted within the Mahalapye town and not near the river crossing which is the likely hotspot. The results from the NDReLI1 and NDReLI2 are, respectively, presented in Figure 12c,d. The NDReLI1 results (Figure 12c), which are based on SWIR1, are observed to overestimate the spatial extent and degree of susceptibility to soil erosivity, with most areas within the railway segments having high to very high vulnerabilities. This may not be the actual case, since the normal amount of rainfall received within the region may not trigger such high degrees of vulnerabilities. Compared to the NDReLI2 results in Figure 12d, the SWIR2 is observed to predict the vulnerabilities to soil erosion with high to moderate intensities and is more accurate in locating the hotspots than the NDReLI1, RUSLE and SDI.
The NDReLI results for the case study of the BRC were further validated with field observations to empirically determine the validity of the new spectral index. The validation results are presented in Figure 13a-d for sections L1-L4, as indicated in Figure 10. In Figure 13a-c, the circled areas indicate the detected hotspot area mapped as either high or very high in terms of vulnerability to erosion along the railway corridor. The results in Figure 13a-d show the NDReLI1 and NDReLI2 results for ground truth points where high and very-high susceptibilities are detected along the railway line. Alongside the index images, the transverse horizontal profiles indicating the variability of the index values across the railway line are presented. In Figure 13a-c, the NDReLI1 shows very high disturbance along the cross-section as compared to the NDReLI2 results. NDReLI2 indicates no significant vulnerabilities along the hotspot cross-section. The same empirical results are observed for the three hotspot regions L1-L3, except for section L4 located in the southernmost part of the railway line where the railway line passes through towns including Mochudi, Gaborone, Ramotswa and Lobatse.
In Figure 13d, despite NDReLI1 mapping the river crossing as having a very high vulnerability, the horizontal profile of the index pixels does not reflect any significant variabilities that can result in a very high susceptibility. On the other hand, NDReLI2 predicts the possibility of minimal susceptibility as the rail line crosses the river. These results show that both indices can be significant in the prediction and identification of the presence and susceptibility to erosion. That is, NDReLI1 can be used to point to the general location and degree of susceptibility, while NDReLI2 can be used to infer the spatial locations and extents of the actual risk areas.
For the section L1 (Figure 13a), the transverse horizontal profile from NDReLI1 shows a spike at a specific section on the railway line. The ground-truth shows that along this section, the surrounding area has a higher slope with a bare soil cover and thus a higher susceptibility to erosion effects. For the same section, NDReLI2 results shows that the area has a high susceptibility to erosion, and the corresponding horizontal profile depicts that the section should be marked for inspection especially during high precipitations. Similar observations in section L1 are made from the ground-truth for sections L2 and L3, where the bare farmlands and bare soils with minimal grass/shrubland covers within proximity to the railway line are identified as either subject to high or very high susceptibility to erosion events. These observations further point to the fact that NDReLI1 should be used as an indicator to susceptibility detection, and NDReLI2 together with field validations applied to determine the actual degree of susceptibility in a given section or area.
There are also false alarm results from both the NDReLI1 and NDReLI2 results, as visually observed in section L4, where the railway crosses the river. In this section, the susceptibility to erosion is detected; however, the horizontal profile plots show that there is insignificant susceptibility due to the vegetation cover along the river and in the surrounding areas. Thus, the application and inference from the indices should incorporate both visual image analyses, statistical evaluation using the transverse horizontal profiles and ground-truth inspections.
Overall, from the case study validations, the sections detected as having a high susceptibility using NDReLI2 were confirmed to be more susceptible to soil erosion based on their high degree of bare soil cover, having slightly higher degrees of slope and being closer to water bodies, especially rivers and their close proximity to built-up urban areas.
The NDReLI was further tested for previously reported erosion-induced railway derailments in Namibia (Figure 14a,b). For the case study of Namibia in Figure 14a, the blowing away of dunes near the Walvis Bay was reported to have triggered the derailment in July 2015. This was captured by the SWIR2 in NDReLI2 results and not the SWIR1 in the NDReLI1 results (Figure 14a). The inability of SWIR1 to detect the wind induced erosion confirms the higher sensitivity of SWIR2 to soil and vegetation in semi-arid conditions. In the second case study in Namibia, there was a reported derailment due to erosion between Asab and Tses in December 2020 (Figure 14b). In mapping the areas vulnerable to erosion along the Asab and Tses railway line, both NDReLI1 and NDReLI2 mapped the areas around the river crossings as the most vulnerable to erosion with NDReLI1 having higher index values as compared to NDReLI2 with high to moderate vulnerabilities along the railway line.

Discussions
Previous studies have recommended the use of remote sensing methods in identifying eroded areas and to monitor erosion processes at regional and local levels [27,48,49]. These recommendations have been based on the use of spectral data, vegetation indices and combinations of remote sensing and morphological data. From the results in this study, the proposed spectral index-based NDReLI can predict the spatial location and extent of soil erosivity within the railway corridor. It is however noted that while in few cases the NDReLI indices yield the same results, in most cases the NDReLI1 is observed to overestimate the probability of occurrence of soil erosion as compared to NDReLI2. In overall, the NDReLI performed better than the RUSLE model and the soil degradation index.
Because arid and semi-arid areas are directly exposed to the sun and the field of view of remote sensors, it is possible to monitor and detect the susceptibility to soil erosion, since it is easy to detect the low vegetation cover due to low cloud cover. Additionally in these areas, the reflectance signals from eroded or bare soil surfaces are characterized by their physicochemical properties, which are easily detectable due to their variations in color brightness values and textures and the structures of the soils [50]. Thus, the correlation between the composition of the eroded soils and their recorded surface reflectances can be used to derive the spectral indices that characterize the spatial level of soil erosion [41,43,[50][51][52][53].
In semi-arid areas, previous studies have also demonstrated the significance of spectral indices such as the correlation index, hematite index, vegetation index and shape index in addition to other indices in the characterization of the soil surface conditions using different optical sensors [43,51,54]. While vegetation-based indices have been particularly used in soil erosion studies [23,[55][56][57], they have been considered inadequate, as they are often mixed with the soil background and interfered with by atmospheric conditions [58,59]. Further, the popularly used RUSLE method has also been found to be unable to delineate all erosion, and the use of satellite-based estimations have been recommended as suggested by [27,28].
While various pixel-based and object-oriented remote sensing techniques have been developed for mapping soil erosion, they tend to rely on finer spatial resolution imagery such as QuickBird, WorldView and IKONOS [60,61]. These high spatial resolution imageries are, however, costly and therefore the associated remote sensing techniques may not be implementable in most case studies. In addition, the high spatial resolution images tend to have lower spectral resolutions, rendering them limited in the detection of ground information in the SWIR [61]. Therefore, in the absence of high-spatial resolution data, the use of medium spatial resolution images such as Sentinel-2 MSI and Landsat is still suitable for the detection and mapping of soil erosion phenomena [55,62].
The empirical validation shows positive results in the mapping of the spatial distribution of the erosion-susceptible areas in the different case studies. For soil erosion management by land managers and decision makers who are mostly interested in the quantification of the spatial extent of the soil erosion, hence inferring the degree of susceptibility and not the absolute soil loss values [22], the proposed NDReLI approach provides a fast and useful framework. However, as already noted, the detection of some small erosion features and events may still be a challenge, especially at medium spatial resolutions such as 12.5 m × 12.5 m. The challenge is likely to be compounded by the spectral homogeneity of the non-erosion features and the soil erosion. To resolve this problem and be able to map small and sensitive erosion risks, the enhancement of the spatial resolutions of the medium resolution sensors is recommended.
The study proposes that the NDReLI can be integrated with RUSLE to improve on the prediction and estimation of soil erosion and soil loss in areas with similar climatic and topographic conditions. In terms of practicability, the proposed index provides an easy-to-implement spectral index approach for the mapping of the spatial distribution of soil erosion more so in data scarce regions where physical and conceptual models may not be usable.

Conclusions
This study presents an approach for the detection of soil erosion susceptibility in railway corridors using a new Normalized Difference Railway Erosion Index (NDReLI). For arid and semi-arid environments, the index is based on the maximization of the reflectance of railway-induced landslides using the Landsat SWIR spectral reflectance enhancement of bare soil cover while minimizing the reflectance of blue light by erosion-causing waters. Using ground-truth validations, the results shows that the SWIR1-based NDReLI1 tends to overestimate the degree of susceptibility to erosions as very high, while NDReLI2 results predicts the same sections as being prone to high erosivity events with moderate degrees of susceptibility as mapped using the transverse horizontal profiles. The results shows that the replacement of the SWIR1 with the SWIR2 band reduces the influence of built-up and impervious surfaces on the detection of bare soils, which are susceptible to soil erosions. Compared to the RULSE soil loss index and the soil degradation index (SDI), the NDReLI based on SWIR2 (NDReLI2) can accurately predict the regions within railway corridors that are susceptible to soil erosion. With different degrees of errors of commission from the ground-truth investigations, the empirical results from the NDReLI index illustrates the ability to effectively isolate bare soils from other features and to predict areas that are vulnerable to soil erosion. With positive predictive results, the NDReLI is considered applicable for finer wavelengths from satellite sensors with similar spectral wavelengths and can also be applied in the assessments of the vulnerabilities of main truck roads and to map the extents of bare-soil cover in arid and semi-arid regions. In conclusion, the study provides a quick assessment approach that can be used to monitor and mitigate the impacts of soil erosion on infrastructure networks in a cost-effective manner. For optimal application of the NDReLI, the study recommends applications especially in semi-arid environments with gentle-flat terrains. (3) soil data from the FAO Map Catalog (https://data.apps.fao.org/map/catalog/srv/eng/catalog. search#/home) (accessed on 16 August 2021) and (4) precipitation data from the Department of Meteorological Services (Botswana). The rest of the data used in the study are as presented in this paper.

Conflicts of Interest:
The authors declare no conflict of interest.