Using soil erosion as an indicator for integrated water resources management: a case study of Ruiru drinking water reservoir, Kenya

Functions and services provided by soils play an important role for numerous sustainable development goals involving mainly food supply and environmental health. In many regions of the Earth, water erosion is a major threat to soil functions and is mostly related to land-use change or poor agricultural management. Selecting proper soil management practices requires site-specific indicators such as water erosion, which follow a spatio-temporal variation. The aim of this study was to develop monthly soil erosion risk maps for the data-scarce catchment of Ruiru drinking water reservoir located in Kenya. Therefore, the Revised Universal Soil Loss Equation complemented with the cubist–kriging interpolation method was applied. The erodibility map created with digital soil mapping methods (R2 = 0.63) revealed that 46% of the soils in the catchment have medium to high erodibility. The monthly erosion rates showed two distinct potential peaks of soil loss over the course of the year, which are consistent with the bimodal rainy season experienced in central Kenya. A higher soil loss of 2.24 t/ha was estimated for long rains (March–May) as compared to 1.68 t/ha for short rains (October–December). Bare land and cropland are the major contributors to soil loss. Furthermore, spatial maps reveal that areas around the indigenous forest on the western and southern parts of the catchment have the highest erosion risk. These detected erosion risks give the potential to develop efficient and timely soil management strategies, thus allowing continued multi-functional use of land within the soil–food–water nexus.


Introduction
Soils form a critical component influencing the hydrological cycle and are crucial in supporting and protecting food and energy production. Several sustainable development goals (SDGs) defined by the United Nations (Goals 1-3, 6-7, 11 and 13) are linked to soil health (United Nations 2015). Soil health must therefore be maintained as they are critical elements of the water-energy-food nexus (Lal et al. 2017). However, the ability of soils to provide their functions is under threat due to degradation mainly through accelerated water erosion. Most soils today can be regarded as being in a state of fair, poor or very poor condition (FAO ITPS 2015). This is the result of human activities such as deforestation, overgrazing, intensive tillage on steep terrain, over-cropping and poor soil management practices (Borrelli et al. 2017;Ebabu et al. 2019). While onsite effects of erosion include land degradation and the loss of soil fertility (Kidane et al. 2019;Lambin et al. 2003), off-site effects include siltation, sedimentation, eutrophication of surface water bodies and even induced flooding (Borrelli et al. 2017;Ozsoy et al. 2012). Borrelli et al. (2017) rated African countries around the equator as erosion hotspots. Due to high population growth compounded by an agriculturebased economy, existing forests are subjected to increasing pressure. This has resulted in deforestation to pave way for other competing land uses such as settlement and agriculture (Carr 2004;Mulinge et al. 2016). Deforestation coupled with high rainfall and fragile terrains with exposed soils has led to severe soil erosion. Maloi et al. (2016) showed that a shift in land use, primarily from forests to cultivation, has increased siltation in the Ruiru Reservoir (Kenya), thereby decreasing its storage capacity by 11-14% in 65 years. There is an urgent need to identify areas with high susceptibility to erosion to implement strategies to mitigate erosion to maintain soil health and prevent impairment of water quality of surface water bodies such as rivers, lakes, and reservoirs. Patil (2018) suggested that the assessment of the soil erosion risks can assist in watershed assessment in areas where erosion is a major threat. Moreover, the processes involved in soil loss by water underlie a high spatio-temporal variability which should be taken into account in integrated water resource management (IWRM). IWRM refers to "a process which promotes the coordinated development and management of water, land and land related resources, in order to maximise the resultant economic and social welfare in an equitable manner" (GWP 2000). This denotes that water and land/soil are interrelated (Calder 2006). Problems concerning water quality and quantity cannot be treated in isolation and a "system thinking approach" between the water-soil nexus needs to be adopted. Using average erosion soil rates can camouflage erosion occurring at smaller scales (Hatfield et al. 2017). Furthermore, extreme events in precipitation as a result of climate change have increased the potential for soil erosion in ways that are not well understood and assessed. Hence, proper spatio-temporal analysis could provide a better understanding of the climate variability and soil ecosystem services intersection. Based on such information, decision-makers and land resource managers can develop more timely and targeted cost-effective measures in the frame of IWRM (Mwangi et al. 2016a, b).
Different studies have been undertaken to quantify soil erosion. The most accurate and ideal methods involve conducting field experiments using erosion plots. Here, erosion measurements are carried out under natural conditions (Ampofo et al. 2002;Boardman and Evans 2019). Alternatively, the rain events can be simulated, where parameters such as rainfall intensity, drop size and spatial variability can be adjusted (Duiker et al. 2001;Stroosnijder 2005;Ries et al. 2013;Iserloh et al. 2013). Other methods include the use of radionuclide tracers (Maina et al. 2018). Nevertheless, direct erosion measurements are expensive, unstandardised, time-consuming and limited in terms of spatial and temporal variation (Lal 2001;Stroosnijder 2005;Ries et al. 2013). The setbacks faced by direct measurements have resulted in the use of models to predict soil erosion. These models are primarily classified into three (physical, conceptual and empirical) depending on model input and the extent of the underlying theoretical principles (Igwe et al. 2017;Patil 2018). Physical models take the individual components affecting soil erosion into account, e.g., spatial and temporal variability. Such models include the European Soil Erosion Model (EUROSEM) (Morgan et al. 1998). Conceptual models like the Soil and Water Integrated Model (SWIM) (Krysanova et al. 1998) represent processes in terms of fluxes at different spatial and temporal resolutions. Empirical models are the least data intense and focus primarily on observed data and their responses. These include the Universal Soil Loss Equation (USLE) (Wischmeier and Smith 1978) and its revised version, the Revised Universal Soil Loss Equation (RUSLE) (Renard et al. 2017), which are conventionally applied to quantify long-term average soil loss . They are used to calculate the annual erosion rate using factors which include rainfall erosivity (R-factor), soil erodibility (K-factor), slope steepness and length (LS-factor), cover management (C-factor) and support practices (P-factor). The R-factor highly correlates with rainfall amount and intensity (Schmidt et al. 2016). Hence, it is expected that there will be an inter-annual and seasonal variation of this factor wherefrom dynamic soil erosion risks can be identified. Recent studies on the temporal variability of the R-factor revealed a distinct seasonality influenced by intense rainfall events in Switzerland (Schmidt et al. 2016). The K-factor, defined as soil loss rate per erosion index, expresses the susceptibility of a soil to be detached and transported by rainfall and runoff (Renard et al. 2017). Direct measurements of this factor require long-term soil erosion studies or the use of soil property data (Wischmeier and Smith 1958). Although there exists a dataset on soil erodibility at higher resolution (up to 500 m) for some developed regions such as Europe (Panagos et al. 2014), only coarse estimates are available for Africa (Borrelli et al. 2017). Since the K-factor is a significant parameter in the soil erosion process, its spatial variability across different landscapes should be considered (Zhu et al. 2010). The spatial variability of the K-factor has been mapped through interpolation methods (Addis and Klik 2015;Avalos et al. 2018) and machine learning algorithms such as the cubist method (Panagos et al. 2014). The LSfactor accounts for the effects of slope steepness and length on erosion (Alexandridis et al. 2015). The C-factor represents the effect of plants as well as soil cover, biomass and disturbing activities on erosion. Schmidt et al. (2018) identified crop cover as one of the main triggers of soil erosion as it is highly dependent on their seasonal dynamics and growth curve. The P-factor includes erosion control measures such as contour cropping, terracing, grass strips and all other enhancements that reduce slope length, thereby reducing the amount and rate at which the soil is lost (Renard et al. 2017).
The RUSLE has been reported to generate quantitative estimates of soil loss in ungauged catchments that have been applied in designing sound conservation measures (Asis and Omasa 2007). Moreover, recent applications of soil loss models such as RUSLE have been complemented by the application of GIS and remote sensing technologies. Such technologies include the use of satellite images and geospatial algorithms such as interpolation distance weighing (IDW) (Angulo-Martínez et al. 2009;Gaubi et al. 2017), kriging (Wang et al. 2002;Kouli et al. 2009), regression equations (Panagos et al. 2014) and machine learning available through digital soil mapping (DSM) (Avalos et al. 2018;Taghizadeh-Mehrjardi et al. 2019). Jenny (1941) established that soil physical and chemical properties are influenced by five soil-forming factors, namely, climate, organisms, topography, parent material and time. This can be exploited with DSM which assumes that a relationship exists between environmental conditions and soil properties (McBratney et al. 2003). By using prediction algorithms, a suite of environmental covariates and measured soil properties, this relationship can be obtained and extended to areas without observations. This means that soil properties such as soil erodibility can be mapped at much lower costs and effort. With access to both geospatial algorithms and higher quality input data, models such as RUSLE can be applied at higher temporal and spatial scales (Uddin et al. 2016;Schmidt et al. 2019;Wang et al. 2019). Schmidt et al. (2016) demonstrated that it is possible to predict soil loss dynamics by modelling the intra-annual variability of the R-factor and C-factor within the RUSLE. By applying the RUSLE equation on Swiss grassland with a sub-annual resolution at a national scale, the authors showed that the mean monthly soil loss by water was 48 times higher in summer as compared to winter. However, the effectiveness of the application of RUSLE depends on the availability of datasets. For developing countries, global datasets are sometimes the only data available. Although providing valuable data, this may still be considered too coarse for effective application. But by undertaking purposive soil sampling and applying digital soil mapping, time and spatial variability of soil properties and vulnerabilities can effectively be determined in data-scarce areas such as Africa (Kamamia et al. 2021). In the presented study, we applied a combination of DSM with the RUSLE to estimate erosion loss on monthly time steps for a data-scarce watershed in Kenya, East Africa. We additionally compared the predicted monthly soil loss with estimates from a global soil loss dataset for this region. With the application of the temporally and spatially variable RUSLE prediction, we wanted to address the following questions. (1) What is the impact of rainy periods on erosion and soil loss in the study area? (2) What is the contribution of the different land-use classes (LULC) to sediment yields in the catchment? (3) Can the assessment of spatio-dynamic erosion inform the planning of soil conservation measures to mitigate soil erosion?
This paper is organised into five sections. "Ruiru Reservoir catchment" (following the introduction) describes the study area. "Methodology" summarises the site selection procedure, data collection phase and laboratory analysis. It further presents a step-by-step methodological approach adopted including the plausibility check. "Results and discussion" discusses the spatio-temporal analysis and links this to the major land uses. Final section concludes the impact of the major findings on the soil-water nexus.

Ruiru Reservoir catchment
The Ruiru Reservoir catchment (Fig. 1a) is located in Kiambu County, Kenya. The catchment covers 51 km 2 from the uplands close to the Rift Valley escarpments to the Ruiru Reservoir at the catchment's outlet. A humid highland subtropical climate with wet and dry seasons characterises the local climate. An average annual rainfall of 1300-1500 mm is received in the catchment. Long rains are experienced between March and May, while short rains are experienced between October and December. Daily temperatures range between 13.0 and 24.9 °C. Temperatures are highest from January to March and lowest in July-August (Nyakundi et al. 2017). Nitisols are the dominating soils in the catchment, whereas a small portion of Andosols can be found in the upper part of the catchment (Fig. 1a). These soils are influenced by pyroclastic and igneous volcanic parent material. The region belongs to the tea-dairy zone and subsistence farming is characterised by low-input low-output production (Kamamia et al. 2021). The Ruiru Reservoir in which the Ruiru River drains, located in the lowest part of the catchment, was designed to supply 23,000 m 3 /day of water to the residents of Nairobi, the capital city of Kenya. Maloi et al. (2016) reported that land-use change, majorly from forest to agriculture, has increased sediments in runoff into the reservoir especially during the rainy seasons, thus increasing the treatment costs. Hence, land-use management must address the water insecurities in the catchment to ensure sustainable supply of good quality water (Calder 2006 Figure 2 presents a flowchart of the approach followed in the spatio-temporal analysis. The first step involved the collection and synthesis of raw data. In the second step, the individual RUSLE factors were determined. These factors were then integrated and compared to the Global Soil Erosion Model (GloSEM) in the third step. Step 1: creation of land-use land cover map, data collection and synthesis

Supervised classification of the major land uses in the Ruiru catchment
A Landsat 8 image Operational Land Imager (OLI) for August was downloaded from the Earth explorer (EarthExplorer, USGS 2019) and was used in combination with the semi-automatic classification plugin in QGIS 3.4 (QGIS.org 2019). All the bands were then clipped to the study area size and converted to reflectance (Young et al. 2017). A band set with Band5, Band4 and Band3 (Bands 5-4-3), which represents a temporary virtual raster that allows for the display of composite colours, was created (NASA 2013). Next, the training sites were defined by creating regions of interest (ROI). In addition to georeferenced points obtained during the field campaign, Google Earth was used to increase the number of ROIs. The main defined land-use classes were water (WTR), bare land (BARE), cropland (AG), tea (TEA), grassland and shrubland (GR/SH), built-up (BLD) and forests (FOR). The classification was carried out by using the maximum likelihood algorithm (MLA) (Benediktsson et al. 1990). The MLA is a rule-based algorithm that is based on Fig. 2 Flowchart of the application of RUSLE equation and plausibility check; digital elevation model (DEM), rainfall erosivity (R-factor), soil erodibility (K-factor), slope steepness and length factor (LS-factor), cover management factor (C-factor), support practices factor (P-factor) and Global Soil Erosion Model (GloSEM) (Borrelli et al. 2017) the probability that a given pixel belongs to a particular class. MLA was applied iteratively and with the Band 5-4-3 (see description in Table 1) band combination.
Each time the classification was done, the spectral signatures (reflectance as a function of the shortwave and different objects have unique signatures which can be used for classification) (NASA 1999) were assessed. However, during classification similar spectral signatures may be recorded for different materials which could lead to misclassification. To overcome this, more training sites were delineated to allow MLA to discriminate between the various vegetation cover and between bare land and built-up areas. The post-processing step included the removal of raster polygons smaller than the minimum mapping unit (MMU), which was set at 25. As the GR/SH vegetation was limited, a new class "Natural Vegetation (NV)" was created by aggregating the grasslands, shrublands and forests. A final land-use map with an accuracy of 71% obtained from the confusion matrix was created (Fig. 3). The distribution of the main LULC for determining soil loss was BARE (4.8%), AG (38.3%), TEA (34.8%), NV (16.6%). The remaining distribution represented built-up areas and the Ruiru Reservoir.

Selection of sampling sites
A 30 m sink-filled digital elevation model (DEM) (Earth Explorer, USGS 2019) was used to extract derivatives using the System for Automated Geoscientific Analyses (SAGA-GIS) (Conrad et al. 2015). From the DEM, slope, aspect, plan curvature, profile curvature, relief, elevation, Topographic Wetness Index (TWI) and hill shading derivatives were extracted. Additional data acquired include a Soil map at a scale 1:250,000 from the Soil and Terrain (SOTER) database of the International Soil Reference and Information Centre (ISCRIC) (Batjes 2008) and a geology map from the Survey of Kenya (SOK). All DEM derivatives and additional data were stacked together and split into 30-m segments in ArcGIS (ArcMap version 10.2). This means that each segment was represented by a series of a different set of covariates which included the DEM derivatives, geology, soil, and land-use values. The sampling distribution was determined by using the conditioned Latin hypercube sampling algorithm (cLHS) (Minasny and McBratney 2006) included in the cLHS package (Roudier et al. 2012) of the R statistical software, version 3.5.1 (R Core Team 2019). 90 soil sampling sites whose distribution is presented in the supplementary information (Online Resource 1) were selected for this study. Similarity sites within a 500-m radius for each site were determined. These were set as alternative sampling sites to the originally selected sampling sites. For example, the sites would act as substitute sites where permission to access private land was not granted or where accessibility was constrained due to geographical barriers or safety reasons.

Fieldwork and laboratory analysis
The soil sampling sites were located in the field with the aid of a GPS. Caution was exercised to ensure that samples were collected from the exact coordinates obtained from the selection process. At each site, undisturbed soil samples from the top 0-30 cm were collected by using soil cores for the determination of bulk density, texture and organic carbon content (Kamamia et al. 2021).
For soil texture analysis, the combined wet sieving and sedimentation method was applied. The sand fraction was determined using wet sieving, while the silt and clay fractions were analysed using sedimentation analysis using the Sedimat 4-12 equipment (Umwelt-Geräte-Technik-UGT GmbH, Müncheberg, Germany). The organic carbon content (C org ) was determined by dry combustion using the Vario TOC Cube (Elementar Analysensysteme GmbH, Langenselbold, Germany). A summary of this data is presented as supplementary information (Online Resource 2).

STEP 2: calculating the individual RUSLE factors
Soil loss was calculated with RUSLE (Renard et al. 2017) by multiplying all factors represented below: The data sources and derivation of the RUSLE factors are given in Table 2.

R-factor
The required daily rainfall data for the determination of R month for the years between 2011 and 2017 was obtained from the Kenya Forest Service (KFS) for the Upland Station located in the western part of the Ruiru Reservoir catchment. This is the only existing station in the catchment. These data were complemented with gridded daily rainfall of the Climate Hazard Group Infrared Precipitation (CHIRPS) (Climate Hazards Center-UC Santa Barbara 2020) (Funk et al. 2015). This freely available high-resolution data combines 0.05° resolution satellite imagery with in situ station data to create gridded time series data starting 1981 to near-present. CHIRPS daily precipitation data for the period corresponding to that of the observed daily rainfall was extracted and used to fill any gaps present in the observed data. The extracted Modified Fournier Index (MFI) (Arnoldus 1977) and consequently the R month factor were determined using Eq. 3 (Table 2). For each month, R month was assumed to be spatially static. It was deemed adequate, as the area of the river catchment is 51 km 2 and does not experience a large monthly rainfall variation.

K-factor
Soil texture, permeability and organic matter content of the topsoil were used to determine the K-factor for the 90 sampled points using the Schwertmann et al. (1987) approach (Eq. 4). The cubist method (Quinlan 1992)

combined with
where A is the annual soil loss in t ha −1 yr −1 R is the rainfall erosivity MJ mm ha −1 h year K is the soil erodibility factor t ha −1 R unit −1 LS is the topographic factor (dimensionless) C is the cropping management factor (dimensionless) P is the practice support factor (dimensionless) where A month is the monthly soil loss in t ha −1 month −1 regression kriging (Malone et al. 2017) and selected covariates (Table 3) was applied as DSM to create a continuous map of soil erodibility.
The cubist model divides data into partitions based on rules associated with covariates and fits a regression equation to each subset. Predictions are then determined based on the relative importance of the covariates. Moreover, three parameters must be established: rules-maximum number of partitions allowed, committees-maximum number of boosting iterations, and extrapolation-model constraints (Malone et al. 2017). On the other hand, regression kriging is a spatial interpolation technique that uses a semi-variogram to quantify the spatial structure of residuals (difference between the predicted and observed values) (Ma et al. 2017). In this study, soil erodibility data from the 90 sampling sites were split into two segments representing the training and testing data (∼ 70% and ∼ 30%, respectively). First, the Cubist model was applied where the data was partitioned based on the most relevant covariates present. This represented the deterministic component of the predictions. Regression kriging, representing the stochastic component, was then undertaken. Finally, these two components were added together to arrive at the final prediction. A leave-one-out cross-validation scheme was applied to assess the accuracy of the predictions, which were represented using the coefficient of determination R 2 . A detailed description can be found in Malone et al. (2017).

C-factor
The monthly C-factors were adjusted based on the normalised difference vegetation index (NDVI) (Jong 1994), which was calculated from Landsat 8 images for the year 2017. The NDVI values range between −1, for almost bare surfaces and water bodies, to 1 for densely vegetated surfaces. The values for NDVIs were afterward converted to the C-factor using Eq. 6 (Almagro et al. 2019), by applying the Raster Calculator tool in ArcGIS (ArcMap version 10.2).  (Arnoldus 1977) pi is the monthly rainfall P is the annual rainfall Soil erodibility, K-factor Measured soil texture, soil organic matter, soil permeability K = 2.77 * 10 −6 * M 1.14 * (12 − OM) + 0.043 * (A − 2) + 0.033 * (4 − D) (4) with: M = (%silt + %f ine sand) * (%silt + %sand(f ine sand excluded)) OM = %Organic matter A = aggregate stability class D = permeability class Cover and management (monthly), C month factor Landsat 8 + (monthly) NDVI = NIR−IR NIR+IR (5) where NIR is the near infrared IR is the reflection in the visible spectrum Slope length and slope steepness, LS-factor Digital elevation model LS = X 22.13 n (0.0065 + 0.045s + 0.0065s 2 ) (7) where X is slope length in meters s is the slope gradient% n is the exponent according to the angle of slope In the months where cloud cover was ≥ 10% (April, May, November), the C-factor values were obtained from spline interpolation. The spline interpolation method is a minimum curvature function that passes through the data with the accuracy of their mean errors. In this method, all the data points influence the value of the interpolated point, with those closest to the main station having the greatest impact on the value of the interpolated point (Niedzielski 2015). Using the polynomials, the first and second derivatives can easily be derived, making them applicable in biological modelling such as developing plant growth curves (Quero et al. 2015).

LS-factor and P-factor
Using the DEM as an input, the LS-factor was determined by using Eq. 7 in Table 2 in ArcGIS (ArcMap version 10.2). Only sporadic conservation practices were observed in the Ruiru Reservoir catchment. This included 'fanya juu' terraces (Mati 2005) on steep agricultural lands and grass strips along the riparian region. To account for this traditional land management practice, a threshold of 25% slope was set for agricultural lands. This means that this conservation measure would be implemented in areas with a slope above 25%. As not all farmers have adopted this conservation measure, 10% of all the possible pixels were randomly selected. Using the dply package in R statistical software version 3.5.1 (R core Team 2019), all pixels were loaded and the first filter according to the slope was implemented. The 10% of remaining pixels were selected using random sampling and a value of 0.7 was assigned according to Angima et al. (2003). Furthermore, the P-factor value for all pixels within 30 m of the Ruiru River was adjusted to 0.9 following Mwangi et al. (2015) to account for grass strips. These have been developed along the Ruiru River as part of an integrated water resource protection initiative by Water Rural User Association (WRUA) (Kamamia et al. 2021). For all other LULC, the P-Factor was set to 1.

STEP 3: plausibility check with GloSEM data
As a plausibility check, the monthly soil loss values (REG_SL) were compared with annual soil loss estimates from GloSEM (Borrelli et al. 2017). As presented in the flowchart (Fig. 2), the individual RUSLE factors were overlaid and multiplied with each other using the Raster Calculator tool in ArcGIS (ArcMap version 10.2). The K-factor, LS-factor and P-factors remained static, while the R-factor and C-factor were substituted for each month. Monthly soil loss was then compared for the different LULC by using the land-use map created from Landsat 8 + data.

Map of soil erodibility
A spatial soil erodibility map with an R 2 of 0.63 for the validation dataset was obtained from the DSM analysis (Fig. 4). The most important predictors were slope, elevation, MrRTF, Band5d, ndvid and GVI. Given that no comprehensive model evaluation guideline was found on the model output, these results were compared with other published works. Avalos et al. (2018) reported R 2 values of 0.31 for DSM by using only terrain attributes and obtained a significant linear relationship between slope and erodibility. Taghizadeh-Mehrjardi et al. (2019) reported R 2 values of 0.71. Panagos et al. (2014) reported an R 2 of 0.74 on cross-validation when applying the cubist method with the multi-level B splines interpolation technique with both spectral and terrain attributes. Although the combined use of these covariates improved the predictions, latitude and elevation covariates were ranked as the most important predictors.
The importance of terrain factors may be attributed to the information they contain on landscape positioning. These covariates can discriminate between hillslopes which are dominated by erosion and transport processes and consequently high erodibility (Gallant and Dowling 2003;Jones et al. 2020). Noteworthy is the preference for MrRTF over MrVTF as a predictor. In general, the RUSLE equation does not account for depositional areas such as valley bottoms represented using the MrVFT . The importance of spectral data may be attributed to the information they contain on vegetation development, which is an input source of soil organic matter. Taghizadeh-Mehrjardi et al. (2019) found a strong negative correlation between the K-factor and soil organic matter. Zhao et al. (2018) and Addis and Klik (2015) reported that soil erodibility has an indirect relationship with vegetation type, which influences the soil organic matter and soil particle distribution. This means that soil erodibility is a dynamic property that can be highly influenced by anthropogenic activities such as changes in land-use or land management (Lal 2001).
According to Schwertmann et al. (1987), erodibility can be classified as follows: K < 0.1-'very low', 0.1 < K < 0.2-'low', 0.2 < K < 0.3-'medium', 0.3 < K < 0.5-'high' and K > 0.5-'very high'. Hence, 46% of the catchment can be ranked as either having medium or high erodibility (Fig. 4). The highest erodibility is observed in areas surrounding the indigenous forest on the western part of the catchment and in the northern and southern extremes (see Fig. 3). Moreover, steeper slopes are classified as having higher erodibility as compared to lower slopes. For the different LULCs (Fig. 5), the highest K-factors were recorded under cropland and bare land. The lowest K-factors were predicted under some tea plantations and not under the indigenous forest as expected. The K-factor is affected by the complex interaction between the different environmental factors, which vary within the different LULC. For the silt-dominated soils of the Loess Plateau of China, Zhao et al. (2018) reported that for the native vegetation, soil properties and topography were the dominant factors which influence soil erodibility. For managed and restored vegetation, soil organic matter highly influenced erodibility (Zhao et al. 2018). Terrain factors such as slope and elevation influence the physical and chemical properties of soils leading to changes in soil particle composition and soil erodibility. Higher elevations were associated with higher erosion of silt material, which is then deposited in the lower elevation (Zhao et al. 2018). Against this background, the presence of andosol (high content of silt and organic matter) in the high elevation could have recorded lower K-factor values due to the effect of long-term erosion. Conversely, the use of high organic amendments such as mulch from pruned residues could explain the low erodibility observed under some tea plantations. The use of mulch protects the soil against Fig. 6 a Water erosion risk maps for the Ruiru Reservoir catchment for January-June. b Water erosion risk maps for the Ruiru Reservoir catchment for July-December the impact of raindrops and increases soil organic matter, which stabilises the soil aggregates and making them less prone to erosion (Ni et al. 2016;Xianchen et al. 2020).

Temporal and spatial soil erosion dynamics
An overview of the spatial monthly soil erosion ( Fig. 6a  and b) shows that most of the soil loss occurs within two distinct periods. The first is between March and June and the second is between October and December. Furthermore, they reveal that the areas in the vicinity of the indigenous forest on the western and southern parts of the catchment have at the highest erosion risk. These are mostly deforested areas that are now either bare land or have sparse shrub cover. Consequently, increased rainfall further increases the risk of erosion along the slopes in the Ruiru catchment. These results can be corroborated by Kamamia et al. (2021) who found that these areas are characterised by low aggregate stability and are therefore highly susceptible to erosion.
From Fig. 7 a cumulative monthly soil loss of 2.23 t/ ha and 1.68 t/ha was estimated from the long rainy season (March 0.11 t/ha/month, April 1.70 t/ha/month, May 0.43 t/ha/month) and short rainy season (October 0.56 t/ha/ month, November 0.89 t/ha/month, December 0.22 t/ha/ month). The highest soil loss was observed in April, while the lowest soil loss of 0.003 t/ha/month was observed in July. Ongoma (2019) reported that precipitation in Kenya reaches its peak in April and is lowest in July. A comparison with the monthly average of 0.27 t/ha/month estimated from GloSEM shows that the average soil loss in April predicted by REG_SL is 6.3 times higher than the average cumulative yearly soil loss from GloSEM. This illustrates that the use of such averages for the design of erosion control measures could lead to the under-designing and construction of ineffective control measures. This situation could further be exacerbated where such measures fail, resulting in damage to crops and even loss of lives. Soil loss can be described with two main mechanisms: erosion by raindrop impact and erosion by surface runoff. In reality, both mechanisms occur simultaneously. Thus, soil erosion mechanisms should be classified based on the degree of susceptibility of one mechanism relative to another (Kinnell 2005). There are two cropping seasons in the Ruiru Reservoir catchment that take advantage of the bimodal rainfall in the catchment. Before the onset of the rainy seasons, the soils are dry and possess high infiltration rates. Despite this, the soils are loose and bare due to tilling. At the start of the rains, the soils are more susceptible to erosion by raindrop impact. Over time, the soil pores gradually get saturated with water and therefore increased precipitation causes more raindrops to penetrate through the flow detaching soil particles which are then lifted and transported (Kinnell 2005). Higher soil loss recorded during the long rains is due to increased precipitation events. Wei et al. (2009) concluded that rainfall events with long durations but relatively low intensities play important roles in inducing severe erosion. Therefore, it is imperative to note that measures protecting the soil against rainfall drops (such as increased plant cover) would markedly reduce soil losses in the Ruiru catchment at the onset of rainfall.

Soil loss dynamics under different LULC
The monthly soil loss trend observed in Fig. 8 for all months was BARE > AG > TEA > NV. Despite occupying the smallest area (4.8%), BARE contributes largest to soil erosion. Online Resource 2 indicates that these areas have the lowest clay percentages making them highly prone to erosion and probably the largest contributor to siltation in the Ruiru Reservoir catchment. This LULC includes some deforested areas and areas around homesteads and roadsides. These roads recorded high soil losses as observed from the maps ( Fig. 6a and b). Most of the feeder roads within the Ruiru catchment occur on steep slopes and are made of earthen material. They are constantly used by heavy vehicles that carry milk and tea from the smallholder farmers to the various processing units. This impact compresses the roads, forcing the edges to break away leaving the bare soils susceptible to water erosion. Cerdan et al. (2010) found a positive relationship between erosion rates and slope length on bare soils. As this LULC records the highest soil loss, vegetative or structural stabilisation of roadsides for example, by use of gabions (Laflen et al. 1985) could offer an all-year solution to mitigate erosion in these areas. The establishment of grasses/turfs around homesteads and reforestation of areas surrounding the forests (see Fig. 3) also recording the high soil loss in March and November could protect the bare soils against erosion.
Cropland with an aerial coverage of 38.3% is the dominating land-use in the Ruiru catchment and recorded the second-highest soil loss throughout the year. Most of the farmers in the Ruiru catchment practice rain-fed mixed annual agriculture. At the onset of the long rains in March and short rains in November, the soils are usually tilled and sowed. Tilling activities destroy soil aggregates and lack of vegetation exposes the soil to the direct impact of raindrops.
Exposure of the already loose soil to rain leads to severe denudation of the topsoil (Feng et al. 2016;Muoni et al. 2020). This variation also changes with the rainfall intensity with higher erosion rates being experienced during higher precipitation (Fig. 8). Adoption of more appropriate crop and tillage activities such as conservation tillage and strip cropping could be used as a strategy to reduce the soil loss recorded at the onset of both the long and short rainfall. Mulching activities can further protect the soil surface from the erosive forces of raindrop impact and overland flow especially during times where crops are less developed. Finally, the adoption of agroforestry can serve as a more permanent solution to reducing erosion by water.
Tea plantations in the Ruiru catchment, covering an area of 34.8%, occur mostly along steep slopes. High soil loss may occur due to the general topography and poor management strategies (Krishnarajah 1985;Mupenzi et al. 2011). In consequence, soil erosion on top of slopes could lead to loss of productivity, which may be irreversible as the rate of soil loss greatly supersedes its formation rate. Soil loss occurs during (1) establishment of new tea plantations, which are undertaken without any conservation measures and (2) management activities such as pruning and weeding often lead to trampling, which loosens the soil, making it more susceptible to erosion. Soil loss in the tea plantations is particularly high within the first two months of the long (April and May) and short (October and November) rains. Before these months, productivity is usually low, there is largescale pruning and weeding in the catchment, as a result, the soil is loosened and exposed. This soil is then lost during the next rainy season. Contrary to this, some studies have concluded that well-managed tea bushes or developed tea bushes have erosion rates comparable to natural forests (Krishnarajah 1985;Allaway and Cox 1989). For Sri Lanka, Krishnarajah (1985) reported that conservation measures such as mulching with pruning residues and grass cuttings as well as the use of terraces resulted in almost complete elimination of soil erosion in tea plantations along steep slopes. The author further observed that mulching during replanting reduced soil erosion. Reduction of water erosion risk for this LULC can be attained by applying mulch during the establishment of new tea plantations. Moreover, mulching before pruning and weeding activities could protect the soil against trampling which induces erosion.
The NV LULC with a share of 16.6%, recorded the lowest soil loss throughout the year. NV included indigenous forest located on the western part of the catchment, afforested land as well as grasslands and shrublands. Afforested lands are part of an initiative by the Kenya Forestry Service (KFS) to restore degraded or deforested areas using mixed species of indigenous trees. Some of the afforested areas contain young tree stands. In addition to possessing underdeveloped tree canopies, the spacing between planted trees is usually wider during the establishment phase (Oliveira et al. 2013). This leaves a lot of open spaces where the impact of raindrops is much higher due to the lack of canopy, litter and underdeveloped tree roots (Drzewiecki et al. 2014). Over time, the trees' canopies develop, and the leaf area index (LAI) increases. The increased litter forms a mulch layer protecting the soil surface (Oliveira et al. 2013). Only then, trees are efficiently capable of controlling soil erosion while increasing soil biomass as in the case of indigenous ones. Nevertheless, the major advantages associated with trees, especially indigenous, should be exploited through the adoption of agroforestry systems. During the formative years, the crops can provide ground cover for the young trees. The trees would then later on control soil erosion and increase base flow through increased infiltration (Nair 2008). This was also observed by Ligonja and Shrestha (2015), who concluded that on-farm tree planting contributed significantly to a 7% reduction of areas under very high erosion between 1986 and 2008 in addition to increasing flow during the dry periods in Kondoa area, Tanzania. Likewise, Nambajimana et al. (2020) recommended a reforestation scheme of rapidly growing tree species as an important feature for erosion control in eroded areas of Rwanda. Waithaka et al. (2020) reported that most grasslands and shrublands (now occupying ~ 6% of the catchment) have been converted into settlement areas in the Ruiru catchment. The remaining scattered portions serve as alternate pasture for the dairy cattle in the catchment. Grass and shrublands have shown the potential to reduce soil erosion.
Thus, it is clear that vegetation coverage is an important factor for soil erosion, as it reduces the direct impact of rain drops and alleviates runoff through increased infiltration which finally affects the quality of surface and river runoff (Cerdan et al. 2010;Feng et al. 2016). It should be however noted that extreme rainfall, such as that received in April, may at times increase the uncertainty of the effectiveness of vegetation in reducing erosion (Wei et al. 2009). Finally, the spatial distribution of vegetation on the slopes is a key factor for the Ruiru Reservoir catchment. This affects the source-sink spatial landscape patterns which influence the runoff generation and sediment transport (Liu et al. 2018). For instance, if done properly, adjusting vegetation patterns along slopes could greatly reduce soil losses into streams and the Ruiru Reservoir and therefore improve water quality and preserve the water quantity.

Annual regional soil loss vs GloSEM
An average cumulative yearly soil loss of 3.24 t/ha/yr (range: 2.33-17.52 t/ha/yr) was obtained for the study area while using the GloSEM (Online Resource 3a). Although the REG_SL (Online Resource 3b) recorded a much lower average cumulative soil loss of 1.72 t/ha/yr, the soil loss range was much wider (0.0003-38.29 t/ha/yr). Online Resource 3a-b show that there is a spatial agreement between the two estimates of high erosion risk on the farthest western part of the catchment. However, the differential map further shows that the GloSEM overpredicted soil loss for most of the catchment and underpredicted for the areas where REG_SL recorded high erosion rates. The cumulative REG_SL unlike the GloSEM displayed the small-scale spatial heterogeneity, which exposed the various erosion hotspots. Borrelli et al. (2020) argued that although the GloSEM model provides pioneering assessment to determine potential erosion at global scales, it is heavily data dependent and is not able to capture all the varying conditions to which it is applied given its global scale. Therefore, the coarse nature of the GloSEM grids could have resulted in spatial aggregation /generalisation and the loss of small soil loss inclusions.

Conclusion
Sustainable soil management, in the context of integrated water resource management, improves soil functions which impact positively on the SDGs. Suitable indicators such as risk of soil erosion is needed to transfer the complexity of the soil-water nexus into a format that can be assessed and measured. More often these indicators have not been well explored and little or no data is available in developing countries. Here, we established that through the combination of efficient soil sampling, analysis of remote sensing data and DSM useful information such as maps of monthly soil loss can be developed in data-scarce areas to aid in developing IWRM measures. Moreover, we demonstrated that the static annual erosion could misrepresent the true dimensions of soil loss with averages disguising areas of low and high erosion potential. The monthly erosion risk maps reveal that a monthly average obtained from average cumulative yearly soil loss from the GloSEM average soil loss is 6.3 times higher than the average soil loss obtained from April. Furthermore, the bare land LULC (occupying the smallest areal coverage) was the largest contributor to erosion indicating that at times a small portion of the catchment is responsible for a large proportion of the total erosion.
As the RUSLE does not account for all complex soil erosion exchanges, the results obtained should be taken as estimates and not absolute values (Benavidez et al. 2018). But as highlighted in this study, it is ideal for determining soil loss at landscape areas especially in data-scarce areas. Varying both the R-factor and C-factor greatly impacted the dynamics of soil loss in the Ruiru catchment. While the spatio-temporal distribution of rainfall is largely uncontrollable, the crop management factor is greatly influenced by the type and intensity of human intervention. Thus, it is only natural that erosion control focuses on improving catchment management practices "where" and "when" erosion is highest. This means that the spatio-temporal maps can be used by different stakeholders during the development of watershed management plans (Mwangi et al. 2016b, c). Within these plans, landscape-scale measures including timely allocation of scarce erosion mitigation and protection measures, proper crop selection that reduces erosion, as well as time-dependent planting and harvesting techniques for agriculture, can be purposefullly incorporated. For Ruiru drinking water reservoir catchment, a successful watershed management plan will require joint effort from the different stakeholders: small-scale farmers/communities, large-scale tea estate farmers, The Nairobi Water and Sewerage Company and the Kenya Forest Service. Only then will the sustainable use of soil and water resources and related ecosystem services in the catchment be achieved.