Dvid Gwpedz , Nmso Nymel , Denis Arthur Hughes ,Andrew Roert Slughter , Sukhmni Kur Mntel , Benjmin vn der Wl

a Institute for Water Research, Rhodes University P.O. Box 94, Grahamstown, Eastern Cape, South Africa

b Geography Department, Rhodes University P.O. Box 94, Grahamstown, Eastern Cape, South Africa

Keywords:MUSLE Storage component Sediment yield Erosion South Africa Inxu river

ABSTRACT The proposed site for the construction of a new dam(Laleni Dam)on the Tsitsa River in the Eastern Cape,South Africa falls within a high soil erosion zone, which is expected to affect the lifespan of the dam.Estimating the sediment yield of the catchment is important from a management perspective, and thus we designed this study to assess the applicability of the Modified Universal Soil Loss Equation (MUSLE).The original MUSLE was applied to the Inxu River sub catchment, upstream of a proposed dam, over a summer rainfall period extending from October 2016 to April 2017. A modified version of the MUSLE,incorporating a storage module,was later applied over the same period.Model calibration of the MUSLE with and without a storage model was tested against 7 months of recorded daily suspended sediment discharge data. With a storage model, we achieved a percentage bias, coefficient of efficiency and R2 of 2%, 0.83 and 0.83, respectively, while without storage the performance statistics were markedly lower.These results suggest that a calibrated MUSLE can provide estimates of sediment yield that are valuable for regional planning, and that the inclusion of a storage component further improves the MUSLE application. Given the complexity of sediment generation and transport processes, a relatively simple empirical MUSLE model is a useful tool for initial sediment yield estimation for planning purposes.

1. Introduction

Soil erosion and sediment transport are complex natural processes that are essential for ecological functioning. However,accelerated rates of soil erosion have become a global threat to hydrological and ecosystem functioning.Globally,1903×106ha of fertile soils are affected by water erosion (Arekhi et al., 2012) and 65%of Africa’s agricultural land is degraded by erosion.Much of this eroded soil is transported through river systems into reservoirs(Vanmaercke et al., 2014), and thus, sedimentation has become a major threat to reservoir lifespans and functioning(Kusimi&Attua,2015).

Nearly a quarter of all reservoirs in South Africa have lost up to 30% of their storage capacity to sedimentation (Msadala et al.,2010). The loss of reservoir storage capacity results in wide ranging problems with ecological, social and economic aspects such as a shortage of water for irrigation, domestic and industrial use, as well as a loss in reservoir flood-attenuation capacity. Apart from the consequences for water quantity,excessive sedimentation also affects water quality, such as increased turbidity, with ecological impacts on riverine communities (e.g., through disruption of filter-feeding mechanisms of some macroinvertebrates)and reduction in suitable habitat (Vanmaercke et al., 2011), as well as impacts on human society by increasing the costs of water treatment.Research-based catchment management can bring us closer to achieving the provisioning of clean water and sustainable livelihoods, in agreement with the United Nations Sustainable Development Goals (SDGs) (https://sdgs.un.org/goals). It is therefore essential to identify tools that can estimate rates of soil erosion and soil transport to guide catchment management(Gelegay&Minale,2016) and water resource planning for mitigating accelerated sedimentation of water resources infrastructure.

Erosion and sediment delivery occur at extremely variable spatial and temporal scales, making it very difficult and expensive to achieve representative observations(de Vente et al.,2007;Fryirs,2013).Numerical modelling has therefore evolved as an alternative approach for estimating sediment yields. Numerous models have been developed, including HBV-SED (Lid′en et al., 2001), Soil Loss Estimation Model for Southern Africa (SLEMSA; Elwell, 1978),Chemicals, Runoff, and Erosion, in Agricultural Management Systems (CREAMS; Kinsel, 1980), Areal Nonpoint Source Watershed Environment Response. Simulation (ANSWERS; Beasley et al.,1980), European Soil Erosion Model (EUROSEM; Morgan et al.,1998) and the Universal Soil Loss Equation (USLE; Wischmeier and Smith,1960). The USLE (and later derivatives) is probably one of the most widely used erosion model in the world(Alewell et al.,2019).

The original USLE experienced further development by Williams(1975) who replaced the rainfall factor with a runoff factor as a function of the volume of runoff(Q)and peak runoff rate(qp).The runoff factor represents the energy of flow necessary to detach and transport sediment and it eliminates the need for a delivery ratio to specify how much of the detached soil is transported to the catchment outlet.Williams and Berndt(1977)named the improved model the Modified-USLE (MUSLE). The USLE was further developed by Renard(1997)for application(in the USA)by establishing and improving parameter estimation methods, resulting in the Revised-USLE (RUSLE).

The USLE and its derivatives, RUSLE and MUSLE, have been frequently used to estimate sediment yield in South Africa.Le Roux et al.(2008)used RUSLE and GIS to estimate soil loss for the entire country. MUSLE is also used as a sub-model within established models such as the Soil and Water Assessment Tool (SWAT)(Neitsch et al., 2005) and the Agricultural Catchments Research Unit (ACRU) (Schulze, 2007). Recently, Bryson (2015) linked the MUSLE to the Pitman hydrological model (Hughes, 2013; Pitman,1973) and proposed a sediment delivery model (WQSED). The WQSED model is currently being further tested and refined at the Institute for Water Research (IWR) at Rhodes University, South Africa. The WQSED will be incorporated into the Water Quality Systems Assessment Model (WQSAM) (Slaughter et al., 2015), an existing water quality model that forms part of the freely available Spatial and Time Series Information Modelling(SPATSIM;Hughes&Forsyth,2006) framework.

The MUSLE was applied in this study because it has demonstrated flexibility and worldwide applicability for estimating event related sediment yield(Sadeghi et al,2007,2013;Sadeghi&Mizuyama,2007);however,uncertainties in the application of the model associated with the temporal and spatial scales at which it is applied have been identified (Gwapedza, Hughes, & Slaughter,2018). The MUSLE is used within the present study because of its relative simplicity and flexibility to simulate sediment yields for the Inxu River catchment using limited input data. A similar study was conducted by Le Roux (2018) using Arc-SWAT and a gully model. The interest in estimating the sediment yield of the Inxu River catchment, a tributary of the Tsitsa River,is associated with the catchment being upstream of the site of a proposed dam(Laleni Dam) and the fact that the region is known to suffer from high soil erosion rates (Le Roux, 2018), which could threaten the sustainability of the dam. The study by Le Roux (2018) also included the Inxu River catchment; however, the modelled simulations were not validated against observed data. The current study applied the MUSLE to the Inxu River catchment, including model calibration against some observed data. The previous studies (Le Roux, 2018; Le Roux et al., 2015) also provide independent simulations against which the MUSLE simulations can be compared.

The results of the present study aim to guide landscape restoration and erosion mitigation projects initiated by the South African Department of Environment, Forestry and Fisheries (DEFF;previously known as Department of Environmental Affairs, DEA)prior to dam construction.In addition,the present research aims to establish if the MUSLE can be applied as a relatively simple and effective tool to provide estimates of sediment yield from the catchment. In comparison to MUSLE, other models such as SWAT and(ACRU)(Schulze,2007),require more time,expertise and data to parameterise and run, and data-scarce African catchments generally lack such data. The key objectives focus on MUSLE application, adding a simple storage model to account for important sediment processes, and an assessment of the model’s ability to reproduce catchment erosion patterns.

2. Methods and materials

2.1. Study area

The area of interest is the Inxu River Catchment.The Inxu River is a major tributary of the Tsitsa River, itself a major tributary of the Mzimvubu River.The Inxu is situated in the northeastern part of the Eastern Cape Province of South Africa and flows in an easterly direction where it meets the Tsitsa River at 31.185°S and 28.689°E (Fig.1). The region experiences a temperate subtropical climate (Base et al., 2007) with intense summer rainfall (Van Tol et al., 2016, from November to March (Base et al., 2007). Mean annual rainfall ranges from 625 mm in the lowlands to 1327 mm in the mountains (Le Roux, 2018). Evaporation peaks in the summer season and can reach a high of 1 628 mm a year(Rutherford et al., 2006).

The Inxu River catchment includes quaternary catchments (the primary drainage basin sub-division used in South Africa for water management purposes) T35F, G, H and J and has a total area of 2062 km2. South African National Land Cover (2014) data suggest that the catchment is dominated by grassland and, to a lesser extent, plantations, croplands and forests. Numerous riverine wetlands are present in the upper catchment, where the primary land use includes forestry and commercial farming.Land use in the lower catchment is dominated by rural settlements where most of the communal land is used for subsistence farming.

The Tsitsa catchment is characterised by Triassic sediments of the Karoo sequence (Base et al., 2007) made up of sandstone and mudstone from the Molteno, Tarkastad and Elliot Formations (Le Roux & Sumner, 2012), which are also evident in the Inxu subcatchment. Though the soils in the area are highly variable in texture and sodic content(Le Roux,2018;du Plessis et al.,2020;Le Roux and van der Waal, 2020), the soils on the lower slopes are associated with highly erodible duplex soils and widespread gully erosion (Le Roux et al., 2015; du Plessis et al., 2020). Duplex soils have a notable increase in clay content from the topsoil to subsoil and an abrupt transition in structure,texture and consistency.Due to the low permeability of the soils, water does not quickly percolate to the subsurface, thus increases the frequency of subsurface flow.The dispersive nature of the soil at this horizon results in pipe and subsequent gully erosion(Le Roux,2018;Le Roux and van der Waal,2020).According to Le Roux et al.(2008),the lower regions of the Inxu River catchment have some of the highest erosion potential in the country.

2.2. Streamflow measurements

The streamflow record was based on an instrumented and rated natural channel cross-section near the outlet of the Inxu River.Stage readings were barometrically corrected.The rating curve was based on five discharge measurements and modelled data using the Manning’s equation. Discharge measurements were based on the velocity-area method during the study period (October 2016 to April 2017). For low-flow measurements (water depth ＜0.7 m),depth-velocity measurements were taken across the channel with an electromagnetic current meter (Gordon et al., 2004). During high flows(river depth ＞0.7 m),floats(oranges)were timed over a 50 m reach at the instrumented site to calculate flow velocity.Depth and cross-sectional area for high flow measurements were based on the depth data and cross-sectional survey. The measurements of total streamflow were separated into surface runoff and baseflow using a statistical baseflow separation method by Hughes et al. (2003).

Fig.1. Map showing the Inxu River catchment upstream of the proposed site of the Laleni Dam in the Eastern Cape, South Africa.

2.3. Suspended sediment measurements

Suspended sediment samples were collected at the channel cross-section using a citizen-based monitoring method, as described by Bannatyne et al. (2017). A local resident collected suspended sediment samples using a depth-integrating sampling pole (Gordon et al., 2004) with a 400 ml jar attached and the nozzle facing the flow. Samples were collected sub-daily (early morning and late afternoon) during moderate to low flow conditions (Q ＜20 m3s-1) and at 30-min intervals during high flows(Q ≥20 m3s-1)when the majority of the sediment is transported,in order to capture the rising and falling limb of the stream hydrograph.Regular(under moderate to low flow conditions)and event-based (high flow conditions) sampling strategies were applied to improve the representativeness of the collected data.Three repeat samples were collected on Thursday mornings to test for variability in measurement accuracy. Suspended sediment concentration (SSC; mg l-1) was calculated as the ratio of the weight of sediment after the samples were dried(at 60°C)and the volume of water;SSC was then converted to sediment flux(ton d-1)based on a relationship between SSC and discharge(Bannatyne et al., 2017). Turbidity, measured with a UTech TN-100 Turbidity Meter, was used as a surrogate for sediment concentration when the samples had very low SSC.The SSC-turbidity relationship was established using 100 measurements ranging from 10 to 1000 NTU.

2.4. Application of MUSLE in this study

The MUSLE model(Eq.(1))was used to estimate daily sediment yields at the outlet of the catchment. Estimates of sediment yield(Eq. (1)) are given as a product of six factors which include vegetation cover,soil erodibility,erosion control practices,slope length and steepness and erosivity calculated form the runoff volume and peak (Williams,1975):

whereSyis sediment yield(tons)for the entire catchment,Qis the volume of runoff(m3),qpis the peak flow rate(m3s-1),andaandbare location coefficients.Theaandbcoefficients are typically fixed at 11.8 and 0.56, respectively (Williams & Berndt,1977).K,L,S,C,andPare the soil erodibility, slope length, slope steepness, cover management,and soil erosion control practice factors,respectively.

Data were gathered and prepared mainly using readily available GIS data sets from national databases for estimating the MUSLE factors. Gwapedza, Slaughter, et al. (2018) found that the use of readily available spatial datasets simplifies parameter estimation and yields reasonable results. The flow data used to drive the MUSLE comprised of a discharge record extending from October 01,2016 to April 28,2017.A distribution approach based on slope zones was implemented in the model to account for at least some of the spatial variability in the erosion estimates and reduce the scale problems associated with MUSLE application (Gwapedza, Hughes,& Slaughter, 2018). The Inxu River catchment was distributed into three slope zones (High, Moderate and Low [HML]) using the Shuttle Radar Topography Mission (SRTM) 30m digital elevation model (DEM) (https://earthexplorer.usgs.gov/Accessed 01-02-2017) and ArcGIS 10.3.1. The DEM was classified into three( HML)slope zones using a manual classification where 0 to ＜7,≥7 to ＜20 and ≥20°slopes represented the low, moderate and high slopes,respectively (FAO, 2006). The sediment yield was calculated for each slope zone and aggregated to provide the total estimate for the catchment.The model also includes a storage sub-model(described below) that allows the erosion estimates to be routed through several storage units before aggregation.

2.4.1. Runoff factor

The runoff equation applied(Eq.(3))is a modification of Eq.(2)and it estimates sediment yield in t ha-1(Williams&Berndt,1977)using values of total runoff depth and peak runoff in mm h-1:

whereQDis runoff depth in mm,qdpis peak runoff in mm h-1,andAreais the total catchment area in ha.

Runoff depth (QD) and peak discharge (qdp) are calculated separately for each HML slope zone. The calculation procedure is briefly presented in this paper but is described in detail in Bryson(2015).

The model is run on a daily time step using the surface runoff estimates derived from the separated total flow data(Qobs)(in m3s-1), which is first converted to runoff depths (QD):

The total catchment is divided into three runoff zones in an attempt to account for spatial variability in both hydrological response as well as the other components of the MUSLE.The parts of the catchment allocated to each runoff zone were based on the identification of slope zones reported in Section 2.4. The total surface runoff depth (QD) is used to calculate runoff depths for the three zones based on the assumption that the high runoff zone generates 75%more runoff than the moderate runoff zone,which in turn is assumed to generate 75% more runoff than the low runoff zone (Bryson, 2015):

where the subscriptsH,M,andLrefer to the high, moderate, and low runoff zones, QD?Is the assumed runoff depth from the different zones and PA?refers to the proportion of the total catchment area assigned to each zone.

To obtain an estimate of the peak runoff (qdp), the storm duration associated with each day’s total runoff depth was calculated using Eq. (8):

whereDis the duration (hours),DP,andDconare power and constant parameters, respectively. All three runoff zones are assumed to have the same duration of runoff,and the peak runoffs for each runoff zone are calculated assuming a simple symmetrical triangle(QDH, QDM,and QDL) by the storm duration (D) using Eq. (9):

where the qdpH, M orLvariables are the peak runoff flows for the high, medium, and low runoff zones, respectively(in mm h-1).

2.4.2. Soil erodibility (K) factor

The soil erodibility (K) factor is the amount of soil erosion per unit erosive force (Wang et al., 2016), and determines the susceptibility of the soil to erosional forces. MostKfactor computations are based on the USLE method(Wischmeier&Smith,1978).Since empirical observations are generally not possible for catchment scale studies, soil erodibility calculations are typically based on readily available soil data. Example calculation approaches can be found in the literature for the USLE(Wischmeier&Smith,1978),RUSLE (Renard et al., 1997), and the Erosion Productivity Impact Calculator (EPIC: Sharply & Williams,1990).

In South Africa, the erodibility of various soil forms and series have been rated by the Agricultural Research Council-Institute for Soil, Climate and Water (ARC-ISCW 2005) based on the Binomial Classification System.K-values for South African soils were derived by Schulze and Horan (2007), as part of the South African Atlas of Climatology and Agrohydrology(Schulze et al.,2007),based on the soil form and series information for each land type. Schulze and Horan (2007) applied the RUSLEK(Equation (10)) (Renard et al.,1997) to calculateK-values for soils with limited available information;other variants of RUSLEKwere used where more soil data were available.

Dg is mean particle size diameter, fiis the primary particle size fraction in per cent, and miis the arithmetic mean of the particle size limits of that size. The updated version of the AUTOSOILS model of Pike and Schulze(1995)was used to computeK-values for the entire country(see Schulze and Horan,2007).TheKfactor data was accessed as a 1: 250 000 scale polygon map. These data were processed in ArcGIS following the method of Gwapedza,Slaughter,et al. (2018) and used to parametrise the MUSLEK-values for the Inxu River catchment.

2.4.3. Topographic factor (LS)

The topographic(LS)factor consists of a combined slope length and slope steepness factor. TheLSfactor for this study was determined using a DEM in a GIS environment(ArcMap10.3.1).A review of MUSLE applications by Sadeghi et al.(2013)highlighted that the use of a DEM in GIS greatly improvedLSestimation and higher DEM resolution improved accuracy. The SRTM 30m DEM was obtained from USGS Explorer (https://earthexplorer.usgs.gov/Accessed 01-02-2017). The raw data were processed to be depressionless and the slope and flow accumulation values obtained from the processed DEM were used within the GIS raster calculator to compute theLSfor all grids lying inside the boundary of each HML slope zone using the following equation:

Soon he agreed to “dangle,” sit on the side of the bed to build up sitting tolerance6, energy and balance. He agreed to “work” with physical therapy if I d return “to talk.”

whereFAis the flow accumulation, α is the slope gradient in degrees,and cell size is the DEM resolution.The meanLSof each HML zone was used as an input in the MUSLE. Equation (12) has been widely used to calculateLSby others who reported obtaining reasonable results,including an application in a nearby Umzintlava catchment(Phinzi et al.,2020)and for the entire South Africa by Le Roux et al. (2008).

2.4.4. Vegetation cover (C) factor

TheCfactor is one of the most important factors controlling soil detachment and transport.An increase in vegetation cover reduces the erosion capacity of runoff(Wang et al.,2016),and also reduces sediment concentrations.An experiment conducted by Hudson and Jackson(1959)showed that the difference in erosion rates between bare and vegetated soils can be greater than 90% (Morgan, 2005),suggesting that vegetation cover can drastically reduce erosion rates by reducing the velocity of runoff and trapping sediment particles (Qiang et al., 2016). TheCfactor for this study was estimated using the South African National Land cover (NLC, 2014)database(https://www.sanbi.org/biodiversity/foundations/national-vegetation-map/accessed 21-07-2016) that was developed by the Department of Environmental Affairs (DEA). The data was processed using the Spatial Analyst tools within ArcGIS 10.3.1,and catchment specific land cover was derived. TheCfactor was calculated by extracting land cover attributes, assigningCfactor values based on look-up tables (Wischmeier & Smith, 1978), and weighting theCfactor values for each land cover category against the total area of each HML slope zone. A summation of the weightedCfactor values for each HML zone is used as a mean vegetation cover value for the Inxu River catchment.

2.4.5. Erosion control practice (P) factor

ThePfactor is associated with conservation practices on cultivated land. Although no active database has information on the practice factor for the catchment, field observations indicated that contour ploughing is a common practice. ThePvalue is between 0 and 1;zero indicates best practice,while 1 represents no erosion control practices.ThePvalue is set to 1 for all non-cultivated land uses such as rangelands and forests. Best practises reduce the impact of runoff velocity that drives the erosion processes (Jang et al., 2015). The catchmentPfactor was derived by using a DEM to calculate catchment slope categories (Gwapedza, Slaughter,et al., 2018) and then estimating thePfactor using the slope categories and a map of cultivation (extracted from the NLC, 2014) by overlaying the cultivation map on the slope categories map (Jang et al., 2015). This process results in a layer that specifies the percentage of cultivated land in each slope category and Wischmeier and Smith (1978) guidelines were used to assign thePfactor values for cultivated land in each slope category, considering that contouring is practised in the study area.The resultingPfactor(for each HML zone)was weighted against the total area of that zone in order to derive the meanPfactor that was input into the MUSLE.

2.4.6. Sensitivity analysis

Environmental models should be subjected to sensitivity analysis to determine how sensitive the model results are to changes in the parameter values. For theK,LS,C,andPfactors in the MUSLE,this analysis is trivial as they all represent linear scaling factors in the equation (Eq. (1)). The important issue for these factors is therefore the likely range over which they could vary depending on the methods for estimating them. The runoff factorRis, however,calculated as a non-linear relationship with runoff depth and peak flow, while peak flow is based on estimates of expected storm duration, which are highly uncertain and require some form of calibration. To assess this uncertainty, simulations derived from a range of possible storm durations were compared with the calibrated duration based on the best fit to the observed sediment yield data. The plausible range was based on assumptions that a maximum storm duration would be about 18 h(experienced during frontal weather systems), while a minimum of 4 h might be expected for intense thunderstorms.

A simple uncertainty analysis was conducted on the erodibilityCandKparameters to analyse the potential equifinality associated with the MUSLE model structure. TheCandKparameters were considered most uncertain becauseCcan vary temporally (largely in response to seasonal or interannual changes in soil moisture availability), and this is not currently considered in the method used to calculateC,whereasKwas estimated from poor resolution soil data.The One factor At a Time(OAT)method(see Pianosi et al.,2016) was used by setting the M zone to values representing feasible mean,minimum and maximumCandKvalues.With the M zone parameters held constant, the H and L zones were systematically changed based on a range of realistic maximum and minimum parameter values (see Table 3). Different parameter combinations that yielded NSE and PBIAS values of ＞0.6 and ＜30%,respectively were noted, from which the best-fit parameter combinations were identified.

2.5. Storage model

A storage function was added as an optional extra component of the MUSLE applied in the current study, based on the assumption that sediment movement through a 2000 km2catchment will inevitably be subject to delays and attenuation.The first step was to estimate the proportion of the daily MUSLE sediment output(MSi)that contributes directly to sediment delivery at the catchment outlet (SDi) on the same day (i). This is based on a power function(parameter SPOW1)of the average of the available sediment at the start and end of the day relative to a maximum possible storage parameter(SSmax).The initial estimate of the available sediment at the end of the day(TSi) is given by:

where Qfactiis the streamflow factor, Qi is the current day streamflow(m3s-1),Qmax is the maximum daily streamflow in the input data, and QPOW is a power parameter. The 1.2 multiplier is included on the assumption that the input data do not represent the highest possible flows and could be considered a parameter,even though it was fixed in this study.Total flow is used to drive the movement of sediment through the storages of the storage model whereas only surface runoff was used to drive the MUSLE as subsurface flow is unlikely to impact on slope erosion. However, the baseflow separation procedure(Hughes et al.,2003)remains highly uncertain as very little data are available for validating outcomes.

The sediment storage factor (SSfacti) is given by a relationship similar to equation(12)using a further power parameter(SPOW2):

The final downstream output from storage (SSDi) is given by equation(18),after which the storage is updated(equation(19))for the following day,and the total sediment delivery becomes the sum of SSDiand SDi,

The main objective of using average sediment storage values in equations (14) and (17) is to smooth the effect of very large sediment inputs from the MUSLE in any one day.

2.6. Assessing model performance

To assess the goodness of fit between observed and simulated data,various objective functions were used.The choice of objective functions may also be guided by the variable being simulated or the observed data that are available.In this study,the statistics include:the coefficient of determination R2,the Nash Sutcliffe coefficient of efficiency (NSE) (Nash & Sutcliffe, 1970), and the percentage bias(PBIAS).

A PBIAS of ＜±15% is considered to be the target for successful calibrations and model performances. Moriasi et al. (2007) highlights that the level of acceptability based on the objective functions varies with the type of variable being simulated.An NSE value of ＞0.50 and a PBIAS ＜±30% for sediment can be judged as good calibrations(Moriasi et al.,2007),and an NSE of greater than 0.75 is rated as very good.

3. Results

3.1. Stream and sediment discharge

The river discharge was measured over a short period, but the discharge patterns show a large variability in daily flows.The river experienced one large discharge between the end of January 2017 and the start of February 2017 (Fig. 2). Other notable discharges occur in mid-November 2016, January 2017 and March 2017. The river exhibited a flashy hydrological discharge pattern that characterises semi-arid catchments in South Africa.Notably,more than 50%of all the discharge is recorded between January and February.The suspended sediment concentrations (SSC) shown in Fig. 2 do not always follow the same pattern as the discharge. The mean daily SSC is 860 mg l-1and the maximum lies above 12×103mgl-1.The largest peak in SSC coincides with the largest peak discharge.However, a significant portion of the SSC peaks do not coincide with high discharges. Fig. 2 illustrates that there is no clear and discernible relationship between the amplitudes of discharge and SCC within the times series.A linear regression indicated a positive but weak relationship between SSC and discharge (R2= 0.39).

3.2. MUSLE model with estimated parameters

The result of the initial application of the MUSLE using parameters estimated from GIS coverages(Table 1)is illustrated in Fig.3a.A strong linear regression was exhibited between the observed and simulated sediment. However, there is a large systematic error between simulated and observed sediment,resulting in a negative NSE and large PBIAS(NSE=-36;PBIAS=1011%).Fig.3b illustrates the effects of forcing MUSLE using total flow instead of separated surface runoff, which produced even worse statistics. However,applying separated surface runoff resulted in days with zero sediment yields, although low sediment outputs were observed in the channel.

The over-estimation noted within the initial application is most influenced by the H slope zone which produced 75% of the cumulative erosion output. The high contribution from the H zone is related to a combination of highLSand runoff factors within this zone. A lowLSfactor for the L slope zone resulted in a meagre 2%contribution from the L zone to the cumulative sediment yield.

Table 1 Estimated and calibrated MUSLE parameters for the Inxu River catchment.

Fig. 2. Time series of observed river discharge and suspended sediment concentration for the Inxu River catchment.

3.3. Calibrated MUSLE model

The results of the sediment yield estimation (Fig. 3c), after calibrating the MUSLE parameters (but without any structural modifications),showed high temporal variability in the catchment sediment yield output. Sediment output was generally lower than 5×104tons d-1for most of the period.However,a large sediment yield event of up to 35 × 104tons d-1occurred from 16 to 19 February 2017.Comparisons with the observed data indicated that the calibrated MUSLE was outside the reasonable range, as exhibited by the high PBIAS (PBIAS = 51%; NSE = 0.76). Though unsatisfactory,the calibration presented a large improvement from the original results (PBIAS=1011% and NSE= -36) obtained when uncalibrated parameters were used. The calibrated MUSLE performance is illustrated by the sediment yield time series (Fig. 3c),which shows that the majority of recessions over-simulate the sediment load.

The table of calibrated parameters(Table 1)shows the variation of MUSLE factors across the slope zones. The calibrated MUSLE model was focused on reducing the over-estimation noted in Section 3.2 by reducing MUSLE scaling parameters. TheCfactor was the primary focus of the calibration because theCfactor estimation method (Section 2.4.4) yielded a highCfactor that is believed to have caused a large over-simulation when estimated parameters were applied in Section 3.2. The high estimatedCfactor is influenced by the land use map that broadly classified a large proportion of the H zone as grassland. Unlike degraded grasslands in the L zone, H zone vegetation is generally pristine grasslands and covered with sparse bush.

Additionally, theCfactor estimation method failed to sufficiently account for spatial and temporal changes in vegetation cover.A separate study conducted by Phinzi et al.(2020)in a nearby Umzintlava catchment(T32E)reported that lowCvalues of 0.01 are generally found in elevated areas, whereas highCvalue are more prominent in the lower areas where settlements are prevalent.The authors applied a spatially distributed Normalized Vegetation Difference Index(NDVI)converted toCfactor.The calibration reduced theCfactors in the H and M zones, and the outcome was a significantly reduced sediment output in the two zones. TheCfactor in the H zone(Table 1)was set to 0.01 because the vegetation composes mainly of undisturbed grass,forest and bush.The L zone had the highestCfactor because anthropogenic activities and grazing can degrade the dominant grass vegetation in that zone.Adjusted factors remained within the minimum and maximum ranges reported in Table 4.

ThePfactor was not adjusted because there is no knowledge beyond the available GIS data that can justify the reduction of the parameter.TheLSfactor is generally expected to remain unchanged(because it is based on analysis of the DEM data); however, it was slightly reduced because of previous analysis (see Gwapedza,Slaughter, et al., 2018) that indicated that the lumped method ofLScalculation might be affected by spatial scale, making the calculatedLSslightly higher.The availableKfactor data used within the model is of low resolution and based on physical knowledge of the catchment, theKwas reduced in the H zone and increased in the L zone in an attempt to account for sediment from gully erosion which is not modelled by MUSLE. The durationsDPandDconparameters were set to 0.8 and 12, respectively; resulting in a maximum storm duration of 13.5 h. The adjusted storm durations were considered reasonable given the size of the catchment and the local climate characteristics. The above calibration approach yielded a relatively good statistical fit when compared to observed data(NSE=0.76; PBIAS=51%). Although the PBIAS remained weak, the result was achieved by adjusting the parameters to match the physical attributes of the catchment as much as possible.

Fig. 3. Graphs (a) and (b) represent the uncalibrated MUSLE forced using separated surface runoff and total flow, respectively. Graph (c) represents the calibrated MUSLE forced using separated surface runoff, and (d) is M-MUSLE, the storage was forced using total flow and MUSLE using separated surface runoff. The y-axis uses a Log scale.

3.4. Modified MUSLE (M-MUSLE) model incorporating storage model

The third sediment yield simulation was conducted using the MUSLE with a storage model. The simulation using the M-MUSLE achieves an improved fit between the simulated and observed data(Figs.3d and 4).The calibration statistics improve significantly and show that the introduction of the storage model improves sediment yield estimation. A notable improvement in the statistics(PBIAS=2%;NSE=0.83 and R2=0.83)versus the original MUSLE calibration statistics (PBIAS = 51%; NSE = 0.76 and R2= 0.77)indicated that the calibrated M-MUSLE produces results are acceptably close to the observed data. The observed cumulative sediment yield for the total period was 1072 × 103tons, while simulated yield was 1622×103tons using the original MUSLE and 1090×103using MUSLE with a storage model.The highest yields of sediment were achieved in January and February.

The application of M-MUSLE yielded a significant improvement in the recession of sediment loads (Fig. 3). The M-MUSLE simulation still reflects the typical variability exhibited in the first simulation using the original MUSLE. The same calibration parameters used in calibrated MUSLE were maintained as these were represented the catchment conditions best and did not need to be changed in order for the storage model to report a good result.The core inputs into the storage model were Qmaxand SSmax.While the Qmaxwas derived from the total flow time-series, the SSmax was set to 900×103based on a calculation of area,bulk density and soil depth.TheQPOWandSPOW1andSPOW2were calibrated to values 1.8, 2.7 and 1, respectively.

3.5. Sensitivity analysis

3.5.1. Storm duration

The results of the sensitivity analysis (Table 2) indicated that model outputs are very sensitive to estimates of storm duration and the way in which these affect the calculations of peak discharge and, therefore, theRfactor in the MUSLE. Shorter duration thunderstorms have higher peak intensity compared to longer duration frontal rainfall, and these inevitably generate more sediment yield(Fig. 5 and Table 2). The calibrated minimum 12-h average storm duration paired with representative catchment parameters yielded the best model fit to observed data. The calibrated duration was considered to be acceptable given the size of the study area and represented a balance between short duration and long duration rainfall.

3.5.2. MUSLE erosivity (R factor)

The results of the analysis of the effect of theRfactor on sediment output are shown in Table 3.When the model is applied with average estimated MUSLE parameters (Table 1), the H zone accounts for 75%of the total sediment,although this zone constitutes only 12%of the total catchment area.When the effect of the linear scaling parameters (LS,K,CandP) was eliminated, the H zone sediment output(Re1;Table 3)was 41%indicating that theRfactor had a significantly strong influence on the sediment output.The H zoneLSfactor (Table 1) was likely the primary influence for the other 34%sediment output driven by the linear scaling parameters.After model calibration, the effects ofR(Re2; Table 3) reduced to 29%in the H zone,and the L zone sediment contribution rose from 2% to 14%. TheRfactor maintains a considerable influence on sediment output as the model is driven by runoff.

Fig.4. Linear correlation of simulated and observed sediment yield(Log X-and Y-axes).Graph(a)represents the original MUSLE output and(b)represents the outputs of M-MUSLE incorporating the storage model.

3.5.3. MUSLE erodibility (C, K factors)

The result of uncertainty analysis indicated that of the ＞3×103model runs performed,more than 500 runs had an NSE ＞0.6 and a PBIAS of ＜30%. The R2remained consistently ＞0.75 across all the runs. Notably, the parameter values used for all the runs were within the maximum and minimum expected ranges (Table 4)based on the physical catchment attributes. Further analysis indicated that it is possible to attain an identical statistical fit using various combinations of theCandKparameters while still maintaining the range of possible maximum and minimum values.Table 5 displays a sample of 15 parameter sets that achieve the same result of NSE 0.7 and PBIAS 1%.

The uncertainty analysis revealed that all the combinations that met the set statistical criteria had theCfactor set to the minimum(Cmin) in the high runoff (H) zone. To be more specific, no combination of parameters was able to achieve the minimum statistical criteria of NSE = 0.5 or PBIAS ＜30 when theCfactor was set to＞0.05 in the H zone, indicating the higher model sensitivity to changes inCin the zone. The high runoff factors could be influencing the sensitivity of the scaling parameters in the H zone.

4. Discussion

4.1. MUSLE factors and sediment yield in the Inxu River catchment

The HML parameters indicate that the variations in MUSLE factors across the slope zones can be quite significant. TheLSfollows the variation of the HML zones as both are associated with the slope. Consequently, the high values of LS (Table 1) are associated with high rates of erosion and sediment yield(Oliveira et al.,2013),although this is dependent on the distribution of other factors.ThePfactor exhibited the same pattern asLS,and this is attributed to the spatial distribution of agricultural activity within the catchment,which is related to topography.Notably,agricultural activity is not prevalent within the higher slope zones;hence thePfactor is lower in theHandMslope zones.TheKandCfactors gradually get higher with a decrease in slope because the lower catchment areas are characterised by degraded grasslands (Gwapedza, Slaughter,et al., 2018) and highly erodible duplex soils (Le Roux & van der Waal, 2020; Le Roux, 2018). Field studies conducted by the authors revealed that the grasslands that dominate the catchment are degraded by overgrazing, further worsening the erosion risk.

Table 2 Storm duration changes and sediment yield outputs compared to the calibrated minimum12 h duration.

Fig. 5. The sensitivity of changing storm duration on sediment yield simulations.

Table 3 A summary of the distribution of sediment output,the effects of the R factor(Re1 and Re2) on sediment output and proportion of catchment area in each HML zone.

Table 4 Maximum and minimum ranges of the C and K factor calculated for the High(H)and low(L) slope zone of the Inxu River catchment.

Table 5 Selected combinations C and K of parameters that gave identical statistics(NSE = 0.7; PBIAS = 1%). C and K factors were held constant for the M slope zone.

Past erosion risk studies have highlighted that the Inxu River catchment is a high erosion risk area in South Africa(Borrelli et al.,2017;Le Roux et al.,2008,2015).Resultant land degradation caused by erosion threatens the viability of agricultural systems (van Tol et al., 2016). Van Tol et al. (2016) further suggests that human activities that include farming could be contributing to high erosion rates because of rapid erosion on abandoned farmland. Water quality problems associated with high sediment discharge on the riverine ecosystem have been identified in the area by Akamagwuna et al. (2019). Additionally, high sediment loads are expected to result in reservoir storage loss as estimated by Le Roux(2018). The problem of erosion and resultant sediment yield on reservoir storage has been identified as a problem in many parts of the world(see Amore et al.,2004;Ezzaouini et al.,2020;Smetanov′a et al.,2020;Yang and Lu,2014)and is likely to impact the planned reservoir that will depend on the Inxu River catchment for inflows.Restoration and rehabilitation of abandoned agricultural lands(see Nearing et al.,2017)could significantly reduce the impact of erosion and sedimentation in the catchment.

Current measurements indicate an observed sediment yield rate of 5 t ha-1.The sediment yield is lower than the African average of 6.35 t ha-1year-1 reported in Vanmaercke et al. (2014). However,the sediment yield is still high and indicative of a rapidly eroding catchment. The sediment yield estimates from the uncalibrated MUSLE yields a 50 t ha-1yield rate within the catchment,which is double compared to the 0.1-25 t ha-1estimated by Le Roux et al.(2015). Both the current and previous estimates are significantly higher than the 5 t ha-1observed yield rate.However,the observed data record was collected over a short period,and perhaps a longer record would provide a different outcome, or perhaps this result suggests that both the current model and previous models of Le Roux et al. (2015) overestimate soil erosion. Within the current study,it was noted that the over-estimation was linked to how theCfactor was estimated.The traditional method(used in this study)of assigningCfactor values from literature to classified land cover images gave results that are constant over large areas (de Asis &Omasa, 2007) and fail to represent the spatial variability of vegetation that is common over large areas(Wang et al.,2002).Errors in classification(Sepuru and Dube,2018)are also introduced into the resultantCfactor map (de Asis & Omasa, 2007). Therefore, it is essential to acquire high-resolution imagery, where possible, for land use classification and account for temporal/seasonal changes in the vegetation(C) factor.

4.2. Evaluating MUSLE applicability in the Inxu River Catchment

The model evaluation showed that M-MUSLE performed better compared to the original MUSLE based on typical model evaluation guidelines. An NSE of 0.83 and a PBIAS of 2%for the M-MUSLE are very good when compared to the 0.5 (NSE) and ±30% (PBIAS)standard for evaluating sediment yield model calibrations(Moriasi et al., 2007). However, the uncertainty analysis (Section 3.5) indicates that it is possible to get a good fit with many different sets of parameters for the MUSLE.The typical equifinality(Beven&Freer,2001)exhibited by the sample parameters in Table 4 indicates that the ’multiplier’ effect of the MUSLE may lead to a good result, but not necessarily for the right reason.Therefore,a modelling exercise involving the MUSLE should be guided by ensuring that good results are associated with parameter values that are a realistic representation of the catchment physical properties and their variation across the slope zones.However, it is also acknowledged that the manner in which the total discharge is distributed across the slope zones has a considerable influence on the model results and is also highly uncertain in the absence of empirical data. The current study indicates that the MUSLE model can be successfully calibrated. While a calibration has been achieved, we need to gather more data to validate the result and assess the stability of the calibrated parameters.Sadeghi and Mizuyama(2007)indicates that the MUSLE would perform well if longer-term observed data are available for calibration.However,observed daily sediment are generally scarce because of financial constraints associated with collecting and processing sediment data;many studies(e.g.Arekhi et al., 2012; Sadeghi & Mizuyama, 2007) resorted to assessing MUSLE performance using a few event based sediment samples.

This study demonstrated that the MUSLE outputs could be further improved by the addition of even a simple storage model.The storage model introduces a storage component that absorbs sediment during low flow events and flushes it out when high flow events occur;a phenomenon that is consistent with highly variable hydrology (in time and space) of most South African catchments(Rowntree et al., 2016). The M-MUSLE results in a significant improvement of the recession sediment loads and enhances model performance.

The applicability of the MUSLE model in the Inxu River catchment shows that a relatively simple method can be used to estimate sediment yield. The result is vital given that large complex models, such as SWAT, require intensive data to operate and these data are usually not available in data-scarce catchments of South Africa. This leads to a great deal of uncertainty in providing more complex models with the required parameters. Given its relative simplicity, the application of the MUSLE requires less modelling expertise. The current study also shows that while the MUSLE calibration performs relatively well, the addition of a storage component greatly improves its performance. However, it is essential to highlight that the MUSLE was applied with observed flow data, whereas SWAT applications involve a hydrological modelling component to drive the sediment processes. The next step is to evaluate how the MUSLE performs with flow data simulated by a hydrological model, and evaluate the possibility of resolving the combined uncertainties emanating from both the hydrological and sediment models.The result reported in this study indicates that simple models such as MUSLE may be calibrated using a relatively short observed data record.However,an extended observed data record will be valuable for validating the calibration.Furthermore, in the context of observed sediment data scarcity, it would be essential to assess if the use of an even shorter record would still achieve a good calibration and a successful validation.

5. Conclusions

This study evaluated the applicability of the MUSLE model in estimating sediment yield in a catchment with high soil erosion rates within the Eastern Cape region of South Africa. The study found that an un-calibrated MUSLE overestimates sediment yield.Overestimation was mainly linked to theCfactor estimation method that fails to account for changing vegetation patterns over time. The calibration effort performed relatively well, yielding some satisfactory performance statistics.However,the addition of a simple storage model significantly improved sediment yield prediction and statistical performance.There remains a critical need to collect more observed data in order to extend model validation.The study also found that the sediment yield of the Inxu catchment was highly variable, and the majority of the total sediment was transported through a few high flow events.However,it should be noted that the success of the MUSLE application in this study may be attributed to the use of observed flow data to drive the MUSLE.Nevertheless,the results indicate that the MUSLE may be applied as a first alternative model to estimate sediment yield, rather than expending more resources on a more complex model. This is a significant result, as simple models are required in most African catchments where both data(to parameterise complex models)are scarce, and modelling expertise is limited.

Declarations of interest

The authors declare no conflict of interests.

Acknowledgements

We are grateful that part of the project has been supported by the Water Research Commission of South Africa under project K5/2448, which also provided partial support for the post-graduate bursary for the first author. Additional post-graduate financial support was provided by the Carnegie Corporation of New York under the Regional Initiative for Science Education (RISE)programme of the Carnegie Corporation of New York and the Oppenheimer Memorial Trust (OMT). The Department of Environment, Forestry and Fisheries (DEFF), Chief Directorate: Natural Resource Management Programmes (NRM), Directorate: Operational Support and Planning provided financial support to conduct the sediment sampling programme and post-graduate funding for the second author.