Investigating Behavior of Six Methods for Sediment Transport Capacity Estimation of Spatial-Temporal Soil Erosion

: Estimation of sediment transport capacity (STC) plays a crucial role in simulating soil erosion using any physics-based models. In this research, we aim to investigate the pros and cons of six popular STC methods (namely, Shear velocity, Kilinc-Richardson (KR), Effective stream power, Slope and unit discharge, Englund-Hansen (EH), and Unit stream power) for soil erosion/deposition simulation at watershed scales. An in-depth analysis was performed using the selected STC methods integrated into the Grid Surface Subsurface Hydrologic Analysis model for investigating the changes in morphology at spatial-temporal scales at the Cheoncheon watershed, South Korea, over three storm events. Conclusions were drawn as follows. (1) Due to the ability of the KR and EH methods to include an additional parameter (i.e., erodibility coefﬁcient), they outperformed oth-ers by producing more accurate simulation results of sediment concentration predictions. The KR method also proved to be superior to the EH method when it showed a more suitable for sediment concentration simulations with a wide range of sediment size and forcing magnitude. (2) We further selected 2 STC methods among the 6 methods to deeply explore the spatial distribution of erosion/deposition. The overall results were more agreeable. For instance, the phenomenon of erosion mainly occurred upstream of watersheds with steep slopes and unbalanced initial sediment concentrations, whereas deposition typically appeared at locations with ﬂat terrain (or along the mainstream). The EH method demonstrated the inﬂuence of topography (e.g., gradient slope) on accretionary erosion/deposition results more signiﬁcantly than the KR method. The obtained results contribute a new understanding of rainfall-sediment-runoff processes and provide fundamental plans for soil conservation in watersheds.


Introduction
Soil loss is a severe problem worldwide causing poor water quality, ecosystem destruction, reduced reservoir storage, and decreased agricultural productivity. The process of soil loss comprises the correlated subprocesses of sediment detachment, transportation, and erosion/deposition caused mainly by rainfall impacts and surface flow [1]. Detachment estimation and net deposition are decided by comparing sediment load with sediment transport capacity (STC). Sufficiently high energy must be available to evacuate a soil particle since it has been detached; otherwise, the particle deposits. As a pivotal input function of physics-based soil erosion models, the STC of overland flow is the maximal equilibrium sediment load that surface runoff can convey for specific hydraulic conditions [2,3]. The above-mentioned environmental impacts of soil erosion have persuaded researchers to develop physics-based and computational models that can simulate the effects of sediment movement. Various physics-based models are accessible for soil loss The schematic diagram of the study is illustrated in Figure 1 and the structure of the paper is designed as follows. Section 2 provides general information about the study site and input data. Methodologies of 6 STC methods and model description/setup are presented in Section 3. Rainfall-runoff-sediment evaluations, the spatial distribution of soil erosion/deposition, and the effect of topography on spatial distribution are visualized and analyzed in Section 4, followed by summarized conclusions in Section 5.

Study Area and Data Availability
The Cheoncheon basin is located in the upper part of the Yongdam dam basin ( Figure 2). With the area corresponding to 290.1 km 2 the sub-basin contains approximately 31% of the Yongdam dam catchment. The Cheoncheon basin is the main inflow supply into the Yongdam Dam which contributes up to 74% of the total inflow [42]. Korean water resources company (K-water) has selected the Cheoncheon basin as an experimental site to evaluate water quantity and quality. The watershed is mainly characterized by a hilly topography with elevation ranging from 277 to 1452 m (Figure 3c). The parent materials are mainly dominated by acidic rocks and metamorphic rocks. The main stream is approximately 28 km long with an average slope of 0.246 m/m. The annual average temperature and humidity are 14 • C and 74%, respectively. Rainfall seasons with high intensity occur at the site between July to August, accounting for approximately half of the total annual rainfall.
Land use/land cover (LULC) information about the study area was extracted from cloud-free Landsat ETM+ image on November 22, 2002 (Path 115, Row 035) and processed via the ArcGIS 10.5 software using a Maximum-likelihood-classifier. As shown in Figure 3a, the classification output depicted the distribution of six categories, including the preeminence of mixed forest (63% of the basin) and mixed field (22%). The soil type data is extracted from National Soil Survey Projects [43]. The soil survey project had started in 1964 with a 1:250,000 Korean soil map as a result. After various phases, the final highly detailed digital soil maps (1:5000) product has been published online in http://soil.rda.go.kr/ (accessed on 8 August 2021). Two main techniques had been applied to determine soil texture including textural analysis in the laboratory and feel method in the field. The soil properties are illustrated in Table 1. Soil types in this study area are dominated by Oesan (19%) and Samgag (39%) (Figure 3b). Based on the percentages of sand, silt, and clay in the soil and the USDA soil taxonomy [44], the prevailing soil types in this study are classified into silty loam and sandy loam, respectively. In the study of Sastre et al. [45], the authors stated that silty loam texture was highly susceptible to causing erosion. The DEM displays cartographic information. Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) with 30-m resolution was used to draw the basin, extract information about the Cheoncheon topography through WMS watershed delineator tools.
Hourly rainfall data for three storm events was obtained from Water Resources Management Information System, South Korea. Those events were historical events that happened in 2002, 2003, and 2007. The characteristics of three events are shown in Table 2. Four meteorological stations and a hydrological station were utilized in this study. Two meteorological stations are located outside the Cheoncheon basin. For model calibration/validation, the observation of flow and sediment discharge at the outlet were extracted from regression equation curves (Table A1) [46][47][48]. Table 3 summarizes information about the data presented in this study.

Six Methods for Sediment Transport Capacity Estimation
Firstly, Kilinc and Richardson [33] investigated soil erosion behavior from surface flow by using an artificial rainfall system at the Engineering Research Center of Colorado State University. The initial KR equation was developed for sand-size particles and then was modified by Julien [49] in 1995 to simulate smaller particles. In 2001, Orden and Heilig [50] advocated for the final version (Equation (1)) to include a reduction parameter for use at event scales.
T c = 25, 500q 2.035 S 1.664 f K 0.15 (1) where T c represents the sediment unit discharge (ton m −1 s −1 ), q represents the unit discharge (m −1 s −1 ), S f represents the friction slope (unitless), and K represents the erodibility factor. The K ranges from 0 to 1 and is used as a calibration parameter. Secondly, Equation (2) was developed by Englund and Hansen and based on flume data and sediment size of bed material as input variables [28]. Equation (2) is to calculate the transport rate in total by multiplying the percentage of each soil particle size fraction (e.g., sand, silt, clay): where P i represents the volumetric sediment transport rate of i-th size fraction (m 3 s −1 ), F j represents the proportion of j-th faction in the layer (0-1), B represents the flow width (m), V represents the mean water velocity (m s −1 ), h represents the flow depth (m), s represents the specific gravity of j-th fraction (unitless), g represents the gravitational acceleration (m s −1 ), and D j represents the mean size of the j-th fraction (m). Thirdly, four STC methods (Equations (3)-(6)) were established by Everaert [17] in relation the of sediments to uncover empirical interrelationships between sediment discharge, slope, velocity, and shear velocity using flumes experiment.
SUD Equation: SV Equation:  (6) where T c represents the sediment unit discharge (g cm −1 s −1 ); q represents unit discharge (cm 2 s −1 ); u * cr represents the critical shear velocity (cm s −1 ); u * represents the actual shear velocity (cm s −1 ); Ω represents the effective stream power; a, b, c, d, and e represent the empirical parameters.

Physics-Based Model: GSSHA
The GSSHA model is a fully distributed, physics-based, and hydrologic model capable of simulating hydrologic processes, water quality analysis, and sediment transport on either event-based or continuous configuration [51]. The GSSHA model is implemented into Watershed Modeling System (WMS) version 11.0 software by Aquaveo, which allows modelers to visualize the model [52]. The model is based on grid-cell structure and uses two-dimensional (2D) diffusive wave equations to calculate surface runoff (e.g., Alternative Direction Explicit (ADE), ADE Predictor-Corrector, and Explicit). The up-gradient explicit method is adopted to simulate one-dimensional channel flow. For infiltration simulation, several approaches (namely, Green and Ampt with soil moisture redistribution, Green and Ampt multilayer, and Richard's infiltration) are applied and integrated into the GSSHA model. Inverse Distance Weighted and Thiessen polygons methods are utilized for spatial distribution of precipitation. Lateral groundwater flow and evapotranspiration (ET) are simulated using 2D vertically averaged and the Penman-Monteith or Deardorff method, respectively. Various auto-calibration methods have been established into the GSSHA model (e.g., Levenberg-Marquardt/Secant LM, Multi-start, Trajectory Repulsion, Multilevel Single Linkage, and Shuffled Complex Evolution (SCE)). The GSSHA model is also a powerful tool for simulating soil erosion/deposition, which involves complicated processes such as detachment by raindrop impact and overland runoff.
Raindrop impact plays an initial role in the detachment and transport of particles. The equation below shows the inter-relationship between rainfall momentum, surface water depth, vegetable cover, and canopy interception [53]: where D R represents the detachment capacity rate (kg m −2 s −1 ), K R represents the soil erodibility factor by raindrop detachment (J −1 ), H w represents the water depth correction factor (unitless), C C represents the canopy cover factor (unitless), C M represents the covermanagement factor (unitless), and M R represents the momentum squared for rainfall Soil detachment by surface flow plays a paramount role in soil loss estimation [54]. The philosophy of the detachment depends on a specific threshold using shear stress that destroys the links between soil particles, as follows: where D F represents the detachment capacity rate (kg m -2 s -1 ), α and β represent the empirical coefficients, Ʈ represents the shear stress of flow (Pa), Ʈ c represents the critica shear stress (Pa), G represents the sediment load (kg m -2 s -1 ), and T c represents the STC of overland flow (kg m -2 s -1 ).

Model Setup and Evaluation
The ArcGIS 10.5 (ESRI) platform [55] was deployed to prepare the necessary input data (e.g., Digital Elevation Model (DEM), LULC, and soil type). To avoid adverse effects resulting from inaccuracies of DEMs, all streams were smoothed and assumed trapezoidal in the section. Topographic information is critical for accurate model predictions. A DEM with a resolution of 250-m is a suitable size for distributed rainfall-runoff modeling [56] thus, the basin was discretized into 5155 computational grids at 250 × 250-m resolution The computation time step in the GSSHA model was selected as 30 s. The ADE method was employed for simulating 2D surface flow, whereas the Green and Ampt with soil where D F represents the detachment capacity rate (kg m -2 s -1 ), α and β represent the empirical coefficients, Ʈ represents the shear stress of flow (Pa), Ʈ c represents the critical shear stress (Pa), G represents the sediment load (kg m -2 s -1 ), and T c represents the STC of overland flow (kg m -2 s -1 ).

Model Setup and Evaluation
The ArcGIS 10.5 (ESRI) platform [55] was deployed to prepare the necessary input data (e.g., Digital Elevation Model (DEM), LULC, and soil type). To avoid adverse effects resulting from inaccuracies of DEMs, all streams were smoothed and assumed trapezoidal in the section. Topographic information is critical for accurate model predictions. A DEM with a resolution of 250-m is a suitable size for distributed rainfall-runoff modeling [56]; thus, the basin was discretized into 5155 computational grids at 250 × 250-m resolution.
where D F represents the detachment capacity rate (kg m -2 s -1 ), α and β represent the empirical coefficients, Ʈ represents the shear stress of flow (Pa), Ʈ c represents the critical shear stress (Pa), G represents the sediment load (kg m -2 s -1 ), and T c represents the STC of overland flow (kg m -2 s -1 ).

Model Setup and Evaluation
The ArcGIS 10.5 (ESRI) platform [55] was deployed to prepare the necessary input data (e.g., Digital Elevation Model (DEM), LULC, and soil type). To avoid adverse effects resulting from inaccuracies of DEMs, all streams were smoothed and assumed trapezoidal in the section. Topographic information is critical for accurate model predictions. A DEM with a resolution of 250-m is a suitable size for distributed rainfall-runoff modeling [56]; thus, the basin was discretized into 5155 computational grids at 250 × 250-m resolution. The computation time step in the GSSHA model was selected as 30 s. The ADE method was employed for simulating 2D surface flow, whereas the Green and Ampt with soil represents the critical shear stress (Pa), G represents the sediment load (kg m −2 s −1 ), and T c represents the STC of overland flow (kg m −2 s −1 ).

Model Setup and Evaluation
The ArcGIS 10.5 (ESRI) platform [55] was deployed to prepare the necessary input data (e.g., Digital Elevation Model (DEM), LULC, and soil type). To avoid adverse effects resulting from inaccuracies of DEMs, all streams were smoothed and assumed trapezoidal in the section. Topographic information is critical for accurate model predictions. A DEM with a resolution of 250-m is a suitable size for distributed rainfall-runoff modeling [56]; thus, the basin was discretized into 5155 computational grids at 250 × 250-m resolution. The computation time step in the GSSHA model was selected as 30 s. The ADE method was employed for simulating 2D surface flow, whereas the Green and Ampt with soil moisture redistribution was used to describe infiltration. The rainfall pattern was spatially distributed using the Thiessen polygons method application. To optimize nine model parameters reported in Table 4, the SCE method was selected [57]. The initial values of these parameters were referred to literature and GSSHA User's Manual [58]. ET was neglected in storm events for the sake of simplicity. To evaluate the model's performance, four statistical indicators have been employed (Equations (9)-(12)). Each criterion has its merits and weaknesses. The coefficient of determination R 2 provides information about the linear relationship between simulated and predicted values. Ranging between 0.0 and 1.0, a higher value of R 2 shows less bias-variance. A value higher than 0.5 is acknowledged to be acceptable in hydrological simulations. The Root Mean Square Error (RMSE) ranges from 0.0 to +∞, and in contrast to R 2 , small RMSE corresponds to a better model. The optimal value of Percent Bias (PBIAS) is 0.0, ranging from −∞ to +∞. Negative values of PBIAS indicate the overestimation of model performance, whereas positive values imply model underestimation bias [59]. The total bias of the volume is illustrated by Volume Conversation Index (VCI) with the best value corresponding to 1.0.
where O k and S k denote the observed and predicted values, respectively; n represents the total number of paired values, and O and S denote the average observed and simulated values, respectively.

Rainfall-Runoff Calibration
The reliability of sediment transport prediction depends on the accuracy of the runoff process. Downer and Ogden [58] suggested that the retention depth, surface roughness, river roughness, soil hydraulic conductivity, soil moisture depth, and the top layer depth are sensitive parameters that should be calibrated in the GSSHA model. As a process-based model, the number of calibrated parameters should be strictly selected to a minimum for robust simulation [60]. The calibrated parameters for estimating flow discharge are listed in Table 4.
Overall, there was a consistent trend of simulated and observed lines of flow discharge in the three events ( Figure 4). The coefficient of determination suggested very high correlations for the model. In the calibration procedure, four values of the statistical indicators illustrated high agreement between the observed and simulated flow discharge (R 2 = 0.90, RMSE = 116.79 m 3 /s, PBIAS = 20.17%, and VCI = 0.998). The compromise between the observation and simulation values gradually decreased during the validation experiments, but was still satisfactory and within the acceptable range (Table 5). Without considering the small rainfall peaks, Event 2 has a more similar rainfall pattern with Event 1 compared to Event 3. This explains why the result of R 2 from Even 2 (0.92) is higher than that of Event 3 (0.73). It is also noted that the RMSE result of Event 1 is the highest among the three events. The observation value of peak flow in Event 1 (1385 m 3 /s) is higher than in Event 2 and 3 (932 and 1131 m 3 /s, respectively); and the duration of Event 1 is also the shortest (Table 2). This can explain the reason why the RMSE value of the calibration event is higher than that of the validation ones since the RMSE detects the mean of the square bias. The PBIAS results ranging from 20.17% to 61.47% in three events indicate that the model tends to underestimate the rainfall-runoff predictions. As a result, the GSSHA model behaved as a potentially robust predictor for flow discharge.  Table 6 reports the results of sediment concentrations produced by six STC methods. Overall, the sediment concentration predictions were less accurate than the flow discharge simulations in both calibration and validation tests. During the calibration and validation processes, the sediment concentration prediction was determined to be applicable only with the KR and EH methods. The former and the latter simulated sediment concentrations with acceptable results in the calibration task ( Figure 5). The values of R 2 , RMSE, PBIAS, and VCI obtained by the KR method were 0.82, 532.65 m 3 /s, 20.17%, and 0.798, respectively. In sequence order, these values were 0.77, 599.65 m 3 /s, 11.55%, and 0.885 for the EH method. The results of R 2 imply that the linear relationship between simulated and predicted sediment concentrations simulated by the KR method is higher than that of the KR method. However, the PBIAS and VCI indexes show that the latter has less bias in terms of total sediment volume predictions than the former. During the validation experiments, the KR method outperformed the EH method. Detailed values are listed in Table 6. According to Johnson et al. [61], these values are among an acceptable range for sedimentation simulations. For the Everaert methods, minimal values of R 2 denoted by "-" pinpointed no linear relationship between simulated and observed values of sediment concentrations. The remaining statistical indicators also illustrated poor performance.    The main reasons for Everaert methods' failure to predict sediment concentration are related to particle size and topography characteristics. On the one hand, the median diameter of grain sizes from Everaert's study varied from 0.033 to 0.39 mm, whereas there is a wider range of grain sizes in the Cheoncheon basin (Table 1). This presents the hypothesis that the Everaert methods could not simulate very fine and coarse textures accurately. The success of stream power theories in STC estimation depends on the grain size of sediment as pivotal input [12]. On the other hand, previous studies had the same agreement that sediment loss and transport capacity are sensitive to the slope gradient [21,62,63]. The slopes (0-42.6 • ) in the Cheoncheon basin are significantly steeper than gentle slopes (<10 • ) in Everaert's study, and approximately 73% of the basin has a slope value higher than 10 • (Figure 3e). Zhang et al. [64] indicate that the significant implications of slope gradient trigger the increment in sediment concentrations. This explains the failure estimation of sediment concentrations using the Everaert methods. In the GSSHA model, these formulas had no parameters for calibration tasks and were likely to vanish in the next release version.

Sediment Transport Evaluation
The EH method is the only method that originated from stream flow conditions in the six mentioned STC methods. In the calibration task, the performance of the EH equation showed a reasonable agreement between sediment concentration values of observation and simulation. The reason probably is that the movement of sand particles in surface runoff situations includes rolling, creeping, sliding, and saltating processes [65], which resemble bedload transport of sand in rivers. The validation periods show a decreasing performance of the EH equation in sediment concentration prediction missions due to the discrepancies in hydraulic conditions between overland flow and stream flow. Surface flow is significantly shallower than river flow, and flow conditions in shallow flows vary temporally owing to surface roughness [66], which depends on the LULC change. The limitation of the study is that uniform LULC has been used in the rainfall-runoff-sediment simulation for three events. This contributes to the declining agreement between observed and simulated values of sediment concentrations. The KR method is robust in predicting sediment concentrations, which agrees with the conclusion of Downer et al. [37]. Experiments conducted by Kilinc and Richardson were similar to the situation in the study watershed with slope gradients ranging up to 22 • .

Spatial Distribution of Soil Erosion/Deposition
Owing to the striking performance in sediment concentration predictions, the KR and the EH methods were spatially analyzed in more detail. For the sake of visual simplicity, spatial information about erosion and deposition was detailed by valued-colored grid cells using the ArcGIS platform. Figure 6 displays spatial-temporal variations in erosion/deposition maps predicted by the KR and EH methods. Overall, the two STC methods agreed on the spatial and temporal dissemination of erosion/deposition. Spatially, the erosion pattern was spotted over the watershed, whereas the deposition pattern was mainly distributed along the river networks. The sediment budget results revealed severe sediment detachments form a rill erosion type occurring around the basin boundaries in steep slope zones, particularly in the southeast and northwest of the watershed. Overall, erosion is more prevalent than deposition, which only accounted for approximately one-third of the basin area shown in the stacked column charts (Figure 7). In Event 1, the percentage of soil erosion simulated by the KR method was 66.0% (approximately 174.1 km 2 ). After 5 years, this percentage increased to 71.8% (approximately 208.9 km 2 ). This assumedly climate-change-related trend [67,68] has also been validated by results produced by the EH method. The values are 64.1% and 67.3%, respectively. The percentage of basin area that impacts soil erosion with magnitude ranging from 0.0 to 0.05 cm is approximately 30-60%. The values for soil deposition are lower, approximately 15-20%. Both soil erosion and deposition with a magnitude higher than 0.5 cm only account for less than 10% of the basin area. This concluded that the soil loss in the Cheoncheon basin is mild and in compromise with previous results of Yu et al. [69]. Since the Cheoncheon watershed was mainly covered by forests (approximately 63%), the rainfall kinetic energy rainfall impacts could be partly precluded by the canopy from directly contacting the land surface [70], and this may reduce the magnitude of soil erosion. The proportions of erosion in 2002 and 2007 were 64.1% and 67.3%, respectively.

Effect of Topography on Spatial Distribution
The LULC, soil type, topography, rainfall, and lithology of the bedrock are the main components affecting soil erosion. However, the topography is of significant influence on soil loss [71]. In Figure 8, slope gradients greater than 15 • were prone to cause erosion. The augment in slope gradient triggered overland flow velocity and accelerated the transport capacity at these sites [72]. The eroded particles carried by surface runoff were getting weaker and then deposited on evener and flatter sites with slope gradients ranging from 0-5 • . The moderate-slope areas (5-15 • ) considering both erosion and deposition values corresponded to transition sites. These findings are consistent with the previous study of Busacca et al. [73]. Figure 9 reveals the relationship between land-use type and slope gradient. Locations of forest and pasture higher than the remaining land-use type were more strongly related to soil erosion. The box-whisker plots also show that the EH method had a higher influence on topography than the KR method in soil loss spatial distribution, even though both methods had the same results in sediment concentration prediction.  Besides the effect of topography, crop and vegetation stems also affected deposition, notably in agricultural areas. In this watershed, cultivated areas are predominantly located on two sides of rivers. Vegetation stems are instigators that intercept runoff and trap sediment particles [3,27]. Thus, crop stems could contribute a significant influence on the performance of the STC in surface conditions and increase the deposition status of the soil. The spatial distribution patterns of soil organic carbon (SOC) in the watershed (Figure 3d) also indicated that dense SOC content was distinguishable from where the deposition was located. This result is consistent with the study of Zhu et al. [74].

Conclusions
In this research, the GSSHA model was deployed to evaluate the rainfall-sedimentrunoff process in the Cheoncheon river basin. It is demonstrated that it is possible to reproduce flow discharge over three storm events. Six STC methods were analyzed for simulating sediment concentrations. The principal outcomes are compiled as follows: (1) From this study, none of the Everaert methods achieved satisfactory results to predict sediment concentrations in three storm events. By allowing the adjustment of an additional parameter, the KR and EH methods surpassed the Everaert methods by producing more accurate results. The KR method also proved to be superior to the EH method when it showed a better performance for sediment concentration simulations with a wide range of sediment size and forcing magnitude.
(2) The spatial distribution of erosion/deposition outputs predicted by the KR and EH methods had similar patterns. While eroded sites were distributed over the watershed, deposited areas were spotted around streams and water bodies. The EH method illustrated the influence of topography on the distribution of erosion/deposition more significantly than the KR method. The results also indicated that the magnitude of soil loss is mild. Slope gradient was the factor spatially contributing to the soil loss distribution. Results depicted that areas with slope gradients ranging from 0 • -5 • were prone to be deposited. By contrast, erosion was likely to happen with a slope gradient higher than 15 • .
The KR method has a simple structure, allowing developers to apply it to their models, whereas the EH method is the only available option for sediment concentration simulations with specific gravity values of particles different from 2.65. However, it is worth noting that applying the EH equation to STC in specific overland-flow areas should be scrutinized conscientiously owing to the divergence in the original establishment of the equation, such as slope gradients and surface conditions. The constraint in the range of sediment sizes and slope gradients contributed to poor performances of the Everaert methods. These formulas could be suited for mild slope areas.
There were some limitations in this study. Firstly, sampling sediment directly in rivers or reservoirs is the most trendy and preferred approach to measure sediment loads with high accuracy. It is not always feasible during extreme storm events due to safety policy and the difficulty of executing such a task. The absence of sediment data has motivated the deployment of sediment rating curves. Nevertheless, such a method is associated with a significant lack of accuracy in extracting sediment loads from flow discharge [75]. Secondly, initial soil moisture values have a significant impact on the results of rainfallrunoff simulation [76]. The wet conditions of soil trigger the overestimation of soil loss, whereas the soil's dry status leads to the underestimation of sediment concentration for physics-based models [77]. The initial soil moisture data partly affected the results since soil moisture data in this study was provided from the literature. Continuous monitoring of sediment and soil moisture installation is thus essential to reduce data uncertainty and enhance proper evaluation of STC methods since observed datasets of spatial erosion/deposition distribution are excluded in this study.
Soil loss is a worldwide concern contributing to the deterioration of agricultural productivity and water quality. The comprehensive analyses presented above could provide pieces of useful knowledge for choosing an appropriate STC approach to estimate soil erosion/deposition at watershed scales. This can be used as a reference to land-use decision-making and to the design of hydraulic installations for soil conservation.