Abstract
The Guizhou Plateau, SW China is largely underlain by carbonate rocks. Because soils are thin, soil loss remains a serious problem despite low erosion rates. Further understanding the impacts of changes in rainfall, land use and differences in topography on sediment yield and delivery may assist in the development of suitable policies to reduce soil erosion on the plateau. A spatially distributed soil erosion and sediment delivery model (WaTEM/SEDEM) was applied to investigate temporal–spatial changes in soil erosion between 1985 and 2014 in three watersheds (Dadukou (DDK), Caopingtou (CPT) and Gaoche (GC)) located in the southwest Guizhou Plateau. The WaTEM/SEDEM model was calibrated and validated using data on sediment yields measured at the watershed scale. The total sediment yield (SY) and soil erosion modulus (SEM) firstly decreased followed by an increase, whereas the sediment delivery ratio (SDR) remained almost unchanged over the 30-year period. The major sediment source was dry farmlands. SY was the highest in the largest DDK watershed. The highest SEM occurred in the CPT watershed due to steep terrain and high ratio of dry farmland areas on steeper slopes. SEM was the lowest in the GC watershed where slope gradient and ratio of dry farmland on steeper slopes are low. SDR was the highest in the GC watershed because of its topographic characteristics. SEM was sensitive to precipitation fluctuations in the GC, DDK and particularly in the steep and intensively eroded CPT watershed, while changes in dry farmland ratio influenced the SEM in the CPT and DDK watersheds but not in the gentle and mildly eroded GC watershed. Changes in forest ratio had significant impacts on SEM only in the GC watershed. Since responses of soil erosion to variations or differences in the main impact factors differ in the different watersheds, soil conservation strategies should be watershed specific.
Keywords
I Introduction
Soil erosion has become an important obstacle to the long-term sustainability of agriculture and ecosystem services (e.g. food productivity, climate regulation, water quality and flood prevention) (Brevik et al., 2015; Guerra et al., 2016). To identify hotspots of soil erosion in the world, special attention should be given to regions where erosion rates are low but soils are thin (Boardman, 2006). Guizhou Plateau (Figure 1a) is such a place. Guizhou Province, SW China (24o30’–29o13’N, 103o31’–109o30’E), with an area of 176,128 km2, is mostly a plateau at 1000–2800 m above sea level (a.s.l.) (Figure 1a). The plateau has a humid, sub-tropical monsoon climate, with mean annual temperature of 15.77±0.43oC and mean annual precipitation of 1184.95±188.12 mm according to the data for 30 years (1981–2010) from meteorological stations in Guizhou Province (http://data.cma.cn). However, only 37.1% of it is covered by forest (NBSC, 2016), mainly because of human-mediated deforestation. The Guizhou Plateau is largely underlain by limestone, dolomite and other carbonate rocks (Yuan, 1997). Pedogenic processes are generally slow on these rocks and soil thickness is usually < 50 cm, which is typically only 30–50% of the soil thickness as observed on non-carbonate rocks under similar climatic conditions (Wan et al., 1995). Moreover, mountains occupy 97% of the plateau (An et al., 1999; Zhou and An, 2000). In many parts of the plateau, thin soil patches on the sparsely vegetated, steep and rugged slopes are rather quickly washed out by seasonal and abundant sub-tropical rainfall, leaving the bedrock exposed and leading to “rock desertification” (An et al., 1999; Wang et al., 2008; Yuan, 1997). Karst landforms are the most common landscape on the Guizhou Plateau. Through fissures, sinkholes and/or other tunnels, varying proportions of the eroded soils may leak into the caves and underground streams (e.g. Ma and Zhang, 2018; Xiong et al., 2012). These soil particles are at least partly trapped in the caves and underground channels and fail to arrive at outlets of watersheds. Therefore, such so-called underground or subterranean soil loss has somehow complicated estimations of total amounts of soil loss by monitoring sediments at watershed outlets. Furthermore, erosion processes are amplified by the implementation of inappropriate agricultural land practices (e.g. over-cultivation and over-deforestation) (An et al., 1999; Xiong et al., 2012). Hence, there exists an urgent need to obtain a thorough understanding of the long-term impact of rainfall, land use and topography on the spatial and temporal dynamics of soil loss on the plateau (Zhang et al., 2018). Furthermore, quantification and mapping of erosion rates at the watershed scale also allows for specific soil conservation policies to be implemented in places where the benefit will be greatest (de Vente et al., 2008; Luo and Wang, 2019).

Location of (a) the Guizhou Plateau in China, (b) the Beipanjiang Basin, (c) the Dadukou (DDK), Caopingtou (CPT) and Gaoche (GC) watersheds within the Beipanjiang Basin.
In recent years, a wide range of studies have evaluated the combined impacts of control measures on sediment yield and soil redistribution by using field surveys and experimental methods (Cao et al., 2013; Evans et al., 2017; Parsons and Foster, 2011; Prosdocimi et al., 2016). However, although traditional field-based studies have provided insight into the response of soil erosion to changes in land use and/or climate, such approaches can be viewed as cost inefficient, because they require multiple years of monitoring in order to obtain reliable soil erosion reference values. In addition, the results are, at best, semi-quantitative, allowing trends to be identified, but the precise relationship and complex interactions between the influencing factors of soil erosion and wide range of geomorphic processes cannot be quantified (López-Vicente and Navas, 2010). To address this knowledge gap, various soil erosion models have been developed and applied to estimate sediment yield, with the majority using the Universal Soil Loss Equation (USLE)/Revised Universal Soil Loss Equation (RUSLE) approach (Pandey et al., 2016; Zhang et al., 2019). A clear advantage of spatially distributed models is that they can be used for implementing measures to prevent soil erosion and sediment generation, since they assess the impacts of changes in land use or climate in a spatial, explicit way. However, such model applications may over- or under-estimate local soil erosion rates, as these models do not consider the details of hydrological and geomorphologic processes determining sediment transport across large spatial scales (i.e. from the within-field to the watershed scale) and across the terrestrial–aquatic continuum (e.g. Fu et al., 2008; Glendell et al., 2018). It is also difficult to assess the accuracy and precision of these models because of the lack of detailed observations restricting the ability to make associated data for reliable estimations of model error, especially conducting long-term soil erosion predictions. Hence, recent attempts to improve the soil erosion models have been focusing on developing a factor database and reducing these sources of uncertainty (Zhang et al., 2019). A remarkable example of such an improved soil erosion model is the WaTEM/SEDEM model, which is a spatially distributed erosion and sediment transport model based on the RUSLE equation (Van Oost et al., 2000). The model predicts sediment delivery to streams using a sediment transport capacity equation and a cascading transport model (Van Rompaey et al., 2005). The WaTEM/SEDEM model is calibrated by comparing the differences between predicted sediment yields and observed sediment data at the outlet of a watershed. This model has been used to simulate soil erosion and redistribution under a wide range of topographic, climatic and land use conditions, including Mediterranean agro-ecosystems (e.g. Quijano et al., 2016; Vente et al., 2008), large agricultural catchments in tropical regions (e.g. Didoné et al., 2017) and a small catchment located in a national park at Lake Balaton, Hungary (e.g. Jordan et al., 2005). Such studies have often generated satisfactory estimations of soil loss and sediment yields. For example, Alatorre et al. (2010) estimated the long-term average annual sediment yield in the watershed of a reservoir with the WaTEM/SEDEM model, and the simulated values of sediment yields were very close to the observed ones; Fang and Sun (2017) simulated soil erosion and its response to the soil conservation measures in a catchment dominated by mollisols in northeastern China with this model, obtaining a satisfactory result for model calibration and valuation. These and other applications suggest that WaTEM/SEDEM is a valid and widely applicable tool for estimating soil erosion and sediment yields. However, applications of the WaTEM/SEDEM model to investigate soil erosion in the carbonate-rock-dominated Guizhou Plateau have, to date, never been reported.
Therefore, we have applied this model to the present study on soil erosion in the southwestern part of the plateau. This study is the first of its kind to calibrate and apply the WaTEM-SEDEM in southwest China, allowing us to determine the relative importance of the specific regional factors influencing sediment production and deposition such as land use, rainfall and topography. The main objective of this study is to (a) calibrate and validate the WaTEM/SEDEM model for its application in three watersheds under varying topography and land-use conditions in the southwestern Guizhou Plateau using annual suspended sediment records; (b) reconstruct sediment production and deposition history of the watersheds during the past 30 years (1984–2015) when rainfall and land use have considerably changed and (c) evaluate the impacts of rainfall, land use and topography on soil erosion by comparing differences of sediment yield, soil erosion modulus, sediment delivery ratio and sediment sources among the watersheds.
II Materials and methods
2.1 Study area
The southwestern part of the Guizhou Plateau has been categorized as a “strongly eroded” region within the plateau based on the second, most recent, province-wide soil erosion survey, with 1:500,000 remote-sensing images completed in 1997 (An et al., 1999). The Beipanjiang Basin (24.85–26.78oN, 103.78–106.23oE) is located in this region (Figure 1b). The area of the basin is 25,958 km2 and has an average slope of 15.4°, with an altitudinal range of 298–2852 m a.s.l. The Beipanjiang River flows southeastwards across the basin and drains into the Pearl River, the third longest river in China (Figure 1b). This basin has a typical humid sub-tropical monsoon climate, with a mean annual temperature of 16.2°C and a total annual precipitation amount ranging between 886 and 1526 mm during the past 50 years, with an overall average value of 1255 mm (Compilation Committee of Guizhou Annals, 2011). The rainy season typically lasts for three months from June to August, during which 80% of the annual rain falls. The annual rainfall decreases from the southeast to the northwest across the basin (Compilation Committee of Guizhou Province Annals, 2011). The basin is generally characterized by heterogeneous relief, vegetation and soils.
The Dadukou (DDK) watershed (25.60–26.78oN, 103.78–104.72oE) is located in the northern part of the basin (Figure 1c). The headwater of the Beipanjiang River is located within this watershed. DDK has an area of 7687 km2 and an average slope gradient of 9.04° (Table 1). This watershed is mountainous, with an altitudinal range of 977 to 2841 m a.s.l. The main land use types are forestland and dry farmland, which account for 44.82% and 24.72% of the total area respectively. The main soil types are yellow brown earth, yellow earth and calcareous soil. Other soil types include red earth, paddy soil and purple soil.
Characteristics of Dadukou (DDK) watershed, Caopingtou (CPT) watershed and Gaoche (GC) watershed in the Beipanjiang Basin.
The Caopingtou (CPT) watershed (25.63–25.98oN, 104.51–105.03oE) is located to the southeast of the DDK watershed (Figure 1c). This watershed has an area of 1249 km2 and an average slope gradient of 9.07° (Table 1), lying at altitudes of 936–2313 m a.s.l. Most of the watershed is characterized by a rolling hill topography, with forestland accounting for ∼40.77% and dry farmland 32.03%. The soil types in the CPT watershed are the same as those in the DDK watershed but are dominated by yellow earth and calcareous soils.
The Gaoche (GC) watershed (25.85–26.31oN, 105.25–106.05oE) is located in the northeastern part of the Beipanjiang Basin (Figure 1c). The area of the watershed is 1788 km2 and average slope gradient is 6.24° (Table 1). The watershed is at 644–1808 m a.s.l. Although a hilly landscape occupies most of the watershed, the overall topography within the watershed is much flatter compared to the other watersheds. Forestland and dry farmland cover 26.6% and 24.2% of the total area of the watershed, respectively. The watershed is overlaid by only three soil types: yellow earth, calcareous soil and paddy soil.
Areas and other characteristics of the three watersheds are listed in Table 1, which highlights the topographical complexity of the three watersheds. There are three hydrological stations, located at each outlet of these three watersheds (Figure 1c).
2.2 Model description
WaTEM/SEDEM is a grid-based model used to estimate long-term mean annual soil erosion. The main aim of the model is to obtain spatial patterns of erosion, transport and deposition within drainage basins of interest and to predict sediment delivery to river channels. It can be applied at different scales from the catchment scale to the continental scale, under a wide range of conditions (Alatorre et al., 2010; Borrelli et al., 2018; Lieskovský and Kenderessy, 2014; Verstraeten and Prosser, 2008). A detailed description of the WaTEM/SEDEM model can be found in Van Oost et al. (2000) (WaTEM), and Van Rompaey et al. (2001) (SEDEM). A brief description of two basic processes of sediment production and sediment deposition is provided below.
The erosion module is used to calculate the annual soil erosion per grid cell, based on a RUSLE approach but working with a 2D application of the length and slope (LS) factor (Desmet and Govers, 1996):
where E is the mean annual soil loss (kg m–2 y–1), R is the rainfall erosivity factor (MJ mm m–2 h–1 y–1), K is the soil erodibility factor (kg h MJ–1 mm–1), LS2D is the dimensional topographic factor (dimensionless), C is the crop management factor (dimensionless) and P is the erosion control practice factor (dimensionless).
The sediment delivery module is used to calculate the sediment yields based on the assumption that sediment on the hillslopes is transported by overland flow, according to the transport capacity equation (Van Rompaey et al., 2001):
where TC is the transport capacity (kg·m–1·y–1), Ktc (m) is an empirical transport capacity coefficient that depends on the land cover and S is the slope gradient (mm–1).
The model uses a mass balance approach to determine the net sediment yield of each grid cell. The amount of sediment transported from a specific grid cell into a downslope adjacent cell is calculated by adding the amount of sediment produced in the grid cell itself to the sediment transported into it from the upslope adjacent cell (Alatorre et al., 2010). The sediment transport capacity (Eq. (2)) is determined by the topographical factors of the considered grid cell. Sediment deposition occurs when the sediment transport capacity is lower than the incoming sediment flux, otherwise the sediment will be transported by overland flow to the downslope adjacent grid cells (Alatorre et al., 2010). However, as the models estimate mainly soil loss resulting from inter-rill and rill erosion processes, other geomorphological processes contributing to the sediment yield at the catchment level, such as permanent gully erosion, bank and channel erosion, and the re-entrainment of landslide sediments, are not explicitly considered (Borrelli et al., 2018; Verstraeten and Prosser, 2008).
2.3 Data
The input to the WaTEM/SEDEM model includes spatial data on topography, extracted from digital elevation model (DEM), land use, climate and soil management in order to create corresponding factor maps of the RUSLE (Renard et al., 1991). These need to be supplied in the form of IDRISI GIS (Clark Labs, Inc., USA) raster layers, and with model parameters dealing with high and low transport capacities (i.e. Ktcmax and Ktcmin model coefficients, respectively).
Short-term records of sediment yields are available in the hydrological stations at the outlets of the three watersheds (Figure 1c), which allow estimations of mean annual sediment volumes used for the WaTEM/SEDEM model calibration and validation. The sediment measurements have been made at the three stations according to the Standards of Measurements for Suspended Sediments of Rivers (GB/T501598-2015) provided by the Minister of Water Resources (MWR) of the People’s Republic of China (PRC). Floods mainly occur (and thus the majority of the sediments are produced) during June, July and August in the study area as well as other parts of the Guizhou Plateau. For each flood, sediments have been measured at least five times, depending on its discharge and duration. Between five and 10 measurements have been made monthly during low-flow seasons. More specifically, sediments have been measured 70–130 times per year at the outlet of the DDK watershed and 33–169 times per year at the outlet of the GC watershed, respectively, since the hydrological stations were set up at the outlets of the catchments. Annual data of the sediment yield based on such measurements should be able to reflect the overall regimes of sediment production in the three catchments.
Moreover, the availability of spatial data on land use, soils, topography and climate for the study area has enabled WaTEM/SEDEM model simulation to be completed for the past 30 years (1985–2014), which was further divided into six periods: P1 (1985–1989), P2 (1990–1994), P3 (1995–1999), P4 (2000–2004), P5 (2005–2009) and P6 (2010–2014). The DEM for the watersheds (30 m×30 m) (Figure 2a, 2d and 2g) were obtained from advanced spaceborne thermal emission and reflection radiometer global digital elevation model (ASTER GDEM) released by NASA and METI (USGS – http://earthexplorer.usgs.gov). The DEM were used to derive information on the spatial variability of topography, associated drainage and river network for the three watersheds. In particular, the 2D slope LS factor of the WaTEM/SEDEM model was calculated (Desmet and Govers, 1996). Specifically, the river network used in this study was generated by topological analysis, with the topographic data derived with the aid of an outline map of the fluvial system framework for the whole Guizhou Plateau. Maps of soil distribution (Figure 2b, 2e and 2h) and data on the properties of these soils (texture and organic matter content) were obtained from the Institute of Soil Science, Chinese Academy of Sciences and were used to calculate and map the K-factor for the study area with the erosion productivity impact calculator (EPIC) method (Sharpley and Williams, 1990).

Maps of Digital Elevation Model (DEM) (a, d, g), soil distribution (b, e, h) and land use (c, f, i) of the Dadukou (DDK), Caopingtou (CPT) and Gaoche (GC) watersheds.
The land use data for the three watersheds were obtained by interpreting Landsat imageries (USGS – http://earthexplorer.usgs.gov) for the midpoint year or so of each of the six time periods (1987, 1992, 1997, 2002, 2007 and 2013). The data may approximately represent the five-year land use conditions of P1, P2, P3, P4, P5 and P6, based on the assumption that the land use has remained unchanged in the short term. Almost all of the visited and consulted local seniors, soil conservation workers and officials working in the governmental agencies of agriculture, land resources and regional planning thought that such a simplification on land use and its changes was acceptable or fairly reasonable. Using ENVI 5.3 software, seven land use categories were identified: (a) water-body, (b) barren and construction land, (c) dry farmland, (d) paddy land, (e) grassland, (f) shrub forest and (g) forestland (Figure 2c, 2f and 2i, and Table 2). Temporal changes in percentages of dry farmland and of forestland, the two types of land use usually strongly affecting soil erosion, were somehow different in the three watersheds (Figure 3a and b). The percentage of the dry farmland was low in P1, increased in P2, was high in P3 and P4, dropped to the lowest values in P5 and increased again in P6 in the DDK watershed (Figure 3a). In the CPT watershed, this percentage was low in P1 and P2, further decreased and was thus the lowest in P3 and P4, and increased distinctly in P5 and P6 (Figure 3a). The percentage of the dry farmland tended to decrease persistently during the six periods in the GC watershed (Figure 3a). The percentage of the forestland was the highest in P1, dropped in P2, further decreased and was hence the lowest in P3 and P4, and increased in P5 and P6 in the DDK watershed (Figure 3b). In the CPT watershed, this percentage was high in P1, low in P2, P3, P4 and P5, and increased and was thus high again in P6 (Figure 3b). The percentage of the forestland was low in P1 and P2, increased sharply and was the highest in P3, and remained high in P4, P5 and P6 in the GC watershed. (Figure 3b).

Temporal changes in the percentages of dry farmland and forestland and annual precipitation in the Dadukou (DDK), Caopingtou (CPT) and Gaoche (GC) watersheds.
Percentages (%) of areas of different types of land use in Dadukou (DDK) watershed, Caopingtou (CPT) watershed and Gaoche (GC) watershed.
The rain erosivity factor (R-factor) was calculated using the relationship between the R-factor and the meteorological parameters as given by Arnoldus (1977), with selected monthly rainfall series from 12 meteorological stations in and around the Beipanjiang Basin (Figure 1c). For each of the six periods, the rainfall data used in calculating the R-factor was the average for the five years of the period. Thus, the calculated R-factor may have largely reflected the average conditions for the whole period. Maps of R-factor values were then generated for the three watersheds based on the Kriging interpolation method as demonstrated in other studies (Oliver and Webster, 1990). The overall trends of the changes in annual precipitation were also explored with the rainfall data for the three watersheds (Figure 3c). The temporal changes in the annual precipitation were the same or very similar across the three watersheds (Figure 3c). The annual rainfall was high in early and mid P1 but sharply decreased and was low in late P1 and early P2. It increased during mid and late P2 but decreased and was low during early P3. In mid and late P3, P4 and P5, the annual rainfall was generally high though it did fluctuate. However, it dropped to the lowest values in early and mid P6. Afterwards, a remarkable increase of the annual rainfall happened in late P6.
The crop management factor (C-factor) is used to define the susceptibility of various land uses and covers to erosion by water. Maps of C-factor are usually generated based on land cover classification maps by assuming that the C-factor value is the same for a certain land cover (Verstraeten, 2006). For the study area, C-factor maps of each category of land use were determined according to relevant reports and local conditions. The C-value was determined as 0 for water-body, as 1 for barren and construction land, as 0.3 for dry farmland, as 0.2 for paddy land, as 0.08 for grassland, as 0.02 for shrub forest and as 0.01 for forestland (Sun et al., 2016; Wang et al., 2014; Xu and Peng, 2008).
WaTEM/SEDEM version 2006 was used in order to estimate the annual sediment yield (SY). Soil erosion modulus (SEM) was calculated by dividing the SY by the watershed area. Sediments produced from soil erosion within a given catchment may be either trapped within the catchment itself or moved away from the catchment by the stream at its outlet. As such, the sediment delivery ratio (SDR) is defined as the proportion of sediments carried by the stream at the catchment outlet to the gross soil erosion across the entire catchment, and hence was calculated by dividing the SY by the total sediment production in the catchment.
The classification criterion of SL190-2007 (MWR of PRC, 2007) divides soil erosion into six classes in terms of erosion modulus, i.e. very slight erosion (0–2 t ha–1 y–1), slight erosion (2–25 t ha–1 y–1), moderate erosion (25–50 t ha–1 y–1), strong erosion (50–80 t ha–1 y–1), very strong erosion (80–150 t ha–1 y–1) and severe erosion (>150 t ha–1 y–1). Since this classification system is the most widely adopted one in soil erosion studies and soil conservation in China, we have decided to use it in this study. This would facilitate comparing the results of this study with those of others in the future. Based on this system, soil erosion was classified and its spatial variability was mapped in each of the three watersheds in terms of the sediment yields predicted with the WaTEM/SEDEM model.
2.4 Model calibration and validation
Two control parameters, high (Ktcmax) and low (Ktcmin) transport capacity coefficient, need to be calibrated in order to ensure the precision of the model. The Ktcmax and Ktcmin represent the slope length under farmland and forestland conditions, respectively. In fact, the predicted amounts of soil erosion by WaTEM/SEDEM are very sensitive to these two parameters (Ktcmax and Ktcmin), which are usually calibrated using sediment yield data (e.g. Alatorre et al., 2010) or using radionuclide proxies such as 137Cs (e.g. Quijano et al., 2016). However, these methods are not always applicable because of cost effectiveness (i.e. establishment of long-term and/or spatially detailed erosion monitoring networks or soil sampling campaigns, e.g. Alatorre et al., 2010). In this study, this kind of detailed measurement was lacking. The model calibration for the three watersheds was performed by using annual data on sediment yield measured at the catchment outlets during 1990–1992 for the DDK watershed, during 1995–1997 for the GC watershed and during 1993–1995 for the CPT watershed. The calibrated model was then validated in order to evaluate its performance with the annual data on sediment yield measured during 1993–1995 for the DDK watershed, during 1998–2000 for the GC watershed and during 1996–1998 for the CPT watershed. A trial-and-error solution was employed for manually calibrating the WaTEM/SEDEM model at discrete steps within a predefined range. Two evaluation criteria were used to assess the simulated sediment yields: (a) relative root mean square error (RRMSE), and (b) Nash–Sutcliffe (NS).
RRMSE is a common tool for assessing model accuracy; the calculation is shown in Eq (3):
where Oi is the observed value and Pi the predicted value.
NS is widely used as a likelihood measure to evaluate model efficiency; it is calculated as shown in Eq (4):
where Omean is the mean observed value and n is the number of observations. NS can range from –∞ to 0. The closer the value of NS is to 1, the more efficient the model is. Thus, the optimum values of Ktcmax and Ktcmin were determined as those that maximized the NS or minimized the RRMSE.
III Results
3.1 Model calibration and validation
The results of the model calibration for the three watersheds using the data on sediment yield measured at the outlets during calibration periods are displayed in Figure 4 and Table 3. The goodness-of-fit plots indicate that there are remarkable trends in the NS and RRMSE distributions (Figure 4). A good convergence of the model to a global optimum point coincides with the maximum NS and the minimum RRMSE, so the model calibration was optimized with this particular parameter setting (Figure 4). For the DDK watershed, when the Ktcmin and Ktcmax parameter of 3 m and 10 m was selected, the model efficiency statistics of NS and RRMSE were 0.59 and 0.07, respectively (Figure 4a and 4d). For the GC and CPT watersheds, the best parameter set was 5 m and 13 m, respectively (Figure 4b, 4c, 4e and 4f). In both cases the NS reached the maximum (0.66 and 0.84, respectively) while the RRMSE reached the minimum (0.13 and 0.41, respectively), which are at least considered satisfactory.

Goodness-of-fit plots of calibration of the transport capacity parameters Ktcmax and Ktcmin for the Dadukou (DDK), Caopingtou (CPT) and Gaoche (GC) watersheds, with annotation of NS and RRMSE trends across parameter settings.
Annual sediment yields predicted by the WATEM/SEDEM model based on the best parameterization of Ktcmax and Ktcmin and their comparisons with observed values for the model’s calibration and validation.
Ktcmax: high transport capacity coefficient; Ktcmin: low transport capacity coefficient; CPT: Caopingtou watershed; GC: Gaoche watershed; NS: Nash–Sutcliffe; RRMSE: relative root mean square error; SD: standard deviation.
The optimal Ktcmin and Ktcmax parameter settings for the CPT and GC watersheds are slightly higher than those for the DDK watershed. However, they are generally lower than the parameter values as calibrated in other studies, such as 7 m and 23 m in a watershed in the central part of the Spanish Pyrenees (Alatorre et al., 2010), 20 m and 50 m in semi-natural catchments in northern Italy (Van Rompaey et al., 2005), and 75 m and 250 m in 21 catchments in southern Flanders, Belgium (Verstraeten, 2006).
According to Van Rompaey et al. (2001), it is not surprising that Ktc values are so variable across different systems, because different routing methods and/or grid sizes may have been considered in the calibration processes. On the other hand, the ratio of Ktcmax to Ktcmin was affected by differences in sediment yields caused by different vegetation cover types (Alatorre et al., 2010). This ratio was 3.33 in the DDK watershed, while the ratio in the GC and CPT watersheds was 2.6. Since the three watersheds are representative of the whole Beibanjiang Basin, the parameter ratio of the model may be also within the range of 2.6–3.33 for the whole basin. Both values fall in the reported range of 1.44 to 3.8 (Alatorre et al., 2010; Fang and Sun, 2017; Van Rompaey et al., 2003; Van Rompaey et al., 2005; Verstraeten et al., 2006, 2007), which also indicates a reliable calibration of the model.
The calibrated model reasonably predicted sediment yields with an NS >0.5 and an RRMSE <0.48. Moreover, there was a very good agreement with Pearson’s R of 0.99 (p < 0.001) between the predicted total annual sediment yields and measured total annual sediment yields, indicating a satisfactory model validation (Table 3). However, it should be borne in mind that the sediment data used in this study to calibrate and validate the model are to a considerable extent offering us a generalized expression of the total soil loss for the entire watershed. Due to the high topographical complexity and heterogeneity in other erosion-determining factors within each of the watersheds, the simulation results as regards the spatial distribution of soil erosion and sediment contribution are apparently less accurate than those of the overall soil loss for the whole watershed. This highlights the need to further improve the WaTEM/SEDEM model in future studies, in particular by taking into account further the intra-watershed variability. As such, it would be more appropriate to consider the simulation results presented in this study as a preliminary estimate of soil erosion and sediment contribution in this type of highly diverse karst landscape.
3.2 Sediment yield, soil erosion modulus and sediment delivery ratio simulated for three watersheds
Total annual SY, SEM and SDR were estimated using the WaTEM/SEDEM model with the calibrated Ktcmax and Ktcmin for the three watersheds (Figure 5). Analysis of variance (ANOVA) and multiple comparison analysis by Fisher’s least significant difference (LSD) were conducted to assess the possible differences among the DDK, CPT and GC watersheds with respect to the mean SY, SEM and SDR (Table 4). The two analyses indicate that SY, SEM and SDR were temporally significantly different (p < 0.001) across the three watersheds.

Sediment yield (SY), soil erosion modulus (SEM) and sediment delivery ratio (SDR) over six periods in the Dadukou (DDK), Caopingtou (CPT) and Gaoche (GC) watersheds. t: metric tons; t/ha: metric tons per hectare.
Analysis of variance (ANOVA) and least significant difference (LSD) statistics of sediment yield (SY), soil erosion modulus (SEM) and sediment delivery ratio (SDR) for the Dadukou (DDK), Caopingtou (CPT) and Gaoche (GC) watersheds.
Among the three watersheds, SY was the highest in the DDK watershed (Figure 5a), which has the largest area (Table 1). During the six temporal periods, SY ranged from 14.89 ×106 to 20.41 ×106 t y–1 with the mean value of 17.42 ×106 t·y–1. The maximum value occurred in P3, and the minimum value occurred in P5. The CPT watershed had a moderate SY ranging from 3.00 ×106 t y – 1 in P4 to 4.72 ×106 t y – 1 in P5. The mean SY for this watershed was 4.15 ×106 t y – 1, which is 76% lower than for the DDK watershed. The total sediment yields were lower in the GC watershed compared to the two other watersheds for all of the six temporal periods. The lowest value of 0.94 ×106 t y – 1 occurred during P3, while the highest values of 2.64 ×106 t y – 1 occurred in P2. Moreover, the SY displayed a higher temporal variability in the DDK watershed, while fluctuations in SY were relatively low at the CPT and GC watersheds.
The temporal trend of variations in SEM was similar to that in SY across the three watersheds (Figure 5b). Moreover, the maximum and minimum values of SEM occurred during the periods of the maximum and minimum SY across the three watersheds. However, the highest SEM occurred in the CPT watershed, ranging from 24.95 to 39.22 t ha – 1, although the values of SY were moderate for this watershed. In contrast, the SEM was medium, with a mean value of 22.24 t ha – 1 in the DDK watershed, where the average SY is the highest. The average SEM was the lowest and ranged from 5.22 to 14.62 t ha – 1 in the GC watershed.
Unlike changes in SY and SEM, the values of SDR for the three watersheds displayed low temporal variability and had changed little over the 30-year period (Figure 5c). In fact, the SDR is a parameter that has been related to the area of the watershed, topographical conditions and relief energy (Chow, 1964). The relatively stable SDRs in the three watersheds imply that the topographic factor played an important role in determining the sediment delivery to the river channels. There are remarkable differences in topography among the three watersheds (Table 1). Furthermore, the highest average SDR-value of 72% was found in the GC watershed. Although the average SDR-value for the CPT and DDK watersheds was much lower than for the GC watershed (i.e. SDR values of ca. 55% and 47%, respectively), the values for the two watersheds (CPT and DDK) were still higher compared to other studies, e.g. 10% in the Barasona Reservoir watershed in the central part of the Spanish Pyrenees (Alatorre et al., 2010), 28% in a catchment in the Czech Republic (Van Rompaey et al., 2005) and 20–39% in the catchments of 164–2173 km2 in Australia (Verstraeten et al., 2007). This means that a large fraction of the soils eroded in these three watersheds were exported by overland flow/streams and not redeposited elsewhere within the watersheds themselves, indicating a high degree of connectivity between the terrestrial and aquatic system as demonstrated for other catchments in previous studies (e.g. Glendell et al., 2018).
Based on the classification criterion of SL190-2007 (MWR of PRC, 2007) and SEM simulated, seven erosion classes were identified and their distributions were mapped in each of the three watersheds (Figure 6). In the DDK watershed, 14.5% (average of six periods) of the area was characterized by strong, very strong and severe soil erosion, which was mostly situated in the southeast and northeast of the watershed (Figure 6a–6f). These parts of the watershed are relatively low in altitude but were mainly occupied by dry farmlands and paddy lands (Figure 2c). The southwestern part of the watershed was characterized by relatively low soil erosion, despite the vast areas of farmlands in this part of the watershed (Figure 2c). In the CPT watershed, soil erosion was generally intensive, with 17.6% (average of six periods) of the area of the watershed dominated by strong, very strong and severe soil erosion (Figure 6g-6l). The areas with relatively low erosion were at the southern and northern part of the watershed (Figure 6g-6l), which was largely covered by forests and shrub forests (Figure 2f). Soil erosion rates were the lowest over the past decades in the GC watershed. Still, there was a small part, lower than 4.23% (average of six periods), of strong, very strong and severe soil erosion located in the southwestern part of the watershed and another small part of intensive soil erosion in northeast of the area (Figure 6m–6r).

Spatial patterns of erosion classified based on the classification criterion of SL190-2007 (MWR of PRC, 2007) and modeled soil erosion modulus over six periods in the Dadukou (DDK), Caopingtou (CPT) and Gaoche (GC) watersheds.
3.3 Major sediment sources
The average soil erosion modulus for each of the five land use types and the ratio of the soil loss from each of the five land use types to total soil loss were calculated for each of the three watersheds to explore the sediment sources and assess contributions from land of different types (Figure 7). The principal sediment source was dry farmland in the DDK watershed, ranging from 43.5% (P3) to 60.1% (P1) as well as in the CPT watershed, ranging from 28.7% (P3) to 49.9% (P1) (Figure 7).

Sediment source ratio and average soil erosion modulus (SEM) metric tons per hectare (t/ha) of five land use types (dry farmland, paddy land, grassland, shrub forest and forestland) over six periods in the Dadukou (DDK), Caopingtou (CPT) and Gaoche (GC) watersheds.
Although the soil loss from the dry farmlands has also made important contributions to the total soil loss over the past decades in the GC watershed, these ratios were considerably lower (i.e. ranging from 27.1% in P4 to 34.3% in P1) and much closer to the soil loss ratio of other land use types (e.g. forestland and grassland). The highest soil loss ratio of dry farmland was found in P1 across all of the three watersheds, whereas the ratio was the lowest during P3 in the DDK watershed, during P3 and P5 in the CPT watershed and during P4 in the GC watershed. Furthermore, the forestland has also made a relatively high contribution to the total soil loss, and therefore was the second largest source of sediments in all of the three watersheds. More specifically, the average soil loss from the forestland in the CPT and GC watersheds accounted for 25.2% and 25.8% of total soil loss, whereas in the DDK watershed this was 17.8%. In contrast, the shrub forest in the DDK watershed accounted for the smallest contribution to the total soil loss (>5.93%). This is the same in the CPT watershed with the exception of P1. The inter-land use distributions of soil loss ratio in the GC watershed are much more complex in comparison to the other two catchments. The soil loss ratio for the paddy land was the lowest in P1, P2, P3 and P5, while the ratio for the shrub forest and grassland was the lowest in P4 and P6 respectively.
The differences in SEM between the five different land use types were notable in the CPT and DDK watershed, whereas this was clearly smaller in the GC watershed (Figure 7). In the CPT watershed, the highest SEM, with a mean value of 62.4 t/ha over the six periods, was found in dry farmland, while the lowest SEM (mean value: 21.5 t/ha) was found in forestland. For the DDK watershed, dry farmland (mean value: 36.8 t/ha) and forestland (mean value: 14.2 t/ha) also had the highest and lowest SEM, but the difference was less evident than the predicted SEM value for the CPT watershed. As the SEM in the GC watershed did not change much, the maximum difference of erosion modulus was only 3.51 t/ha. The soil loss ratio did not match well with soil erosion modulus in these three watersheds. For example, the paddy land, characterized by relatively high SEM, did account for a rather low ratio of soil loss in both the DDK and GC watersheds, while the forestland, characterized by the lowest SEM, did account for a relatively high ratio of soil loss. This may imply that the land use structure is a more important factor for determining the total amount of soil erosion within the watershed, rather than the SEM value of individual land use type.
IV Discussion
4.1 Impacts of precipitation and land use on soil erosion
Changes in precipitation and land use have various impacts on soil erosion (Ochoa-Cueva et al., 2015). Regression analyses between SEM and various influencing factors (i.e. precipitation, dry farmland ratio and forest ratio) were made for the three watersheds (Figure 8). This analysis suggests that the effects of these factors on soil erosion are different in the three watersheds. There is a strong positive relationship between SEM and precipitation in the GC (r = 0.81, p < 0.01), the CPT (r = 0.78, p < 0.05) and DDK (r = 0.64, p < 0.1) watershed (Figure 8a). In fact, soil losses from watersheds are directly proportional to rainfall erosivity (Renard et al., 1997) because the R-factor is a numerical descriptor of the ability of rainfall to erode soil (Wischmeier and Smith, 1978). Therefore, rainfall could be considered as the main factor resulting in the temporal patterns of soil erosion within three watersheds (Figure 5a and 5b). The regression slope identified for the CPT watershed was higher than that for the DDK and GC watersheds (Figure 8a). Thus, soil erosion in the CPT watershed was more sensitive to rainfall and a small increase in annual rainfall could trigger a significant increase in SEM in the watershed. Furthermore, since the soil erosion was mostly sensitive to fluctuations in rainfall in the watershed with the highest SEM, more attention should be given to reducing the impact of rainfall on highly erodible parts in such watersheds. It should be pointed out that, due to the lack of more detailed data on rainfall (particularly those on heavy or torrential rainfall), we have adopted the general R-factor definition and the value used in the simulation is the result of a multi-year average. Hence, the response of soil erosion to rainfall across the different watersheds in the present study reflects a comprehensive response over multiple years. However, in future research it could be of interest to investigate the response considering a much shorter time window if more detailed meteorological and hydrological data is available.

Regression relationships between soil erosion modulus (SEM) and precipitation, metric tons per hectare (t/ha), dry farmland and forest in the Dadukou (DDK), Caopingtou (CPT) and Gaoche (GC) watersheds.
Wan et al. (2004) reported that serious land degradation took place on the Guizhou Plateau as a consequence of intensive land use. It has been found that rapid changes in land use distribution for cultivation purposes could increase the extent of soil erosion by altering the sediment source distribution and enlarging the spatial connectivity (García-Ruiz and Lana-Renault, 2011). In the CPT and DDK watersheds, SEM and dry farmland ratio displayed a general positive correlation, with Pearson correlation coefficient (r) values of 0.71 and 0.47, respectively (Figure 8b). In particular, the regression slope identified for the CPT watershed was higher than that for the DDK watershed (Figure 8b). However, there was no clear correlation between SEM and dry farmland ratio in the GC watershed (r = 0.12) (Figure 8b). Therefore, increases of dry farmland would particularly intensify soil erosion in the steeper CPT watershed that has a higher average slope gradient than the DDK and GC watersheds (Table 1). In contrast, impacts of changes in the percentages of dry farmland could be minimal on soil erosion in the relatively flat GC watershed with its lower average slope gradient (Table 1). In essence, these results underline the relatively high importance of dry farmland coverage as a driver for erosion as compared to other land use types. There was only a significant negative (p < 0.05) relationship between SEM and forest ratio (r = –0.78) for the GC watershed, whereas no clear correlation was found in the CPT and DDK watersheds (Figure 8c). This seems to suggest that changes in forest areas had little impact on the watershed in which soil erosion was intensive and relatively intensive (i.e. the CPT and DDK watersheds). Therefore, the more effective measure to control soil erosion is to decrease the area of dry farmland rather than increase the forest area. Hence, limiting the expansion of farmland and particularly dry farmland should be encouraged to reduce sediment yields in these watersheds (Luo et al., 2018). Ouyang et al. (2018) also suggested that controlling farmland extent could remarkably weaken soil erosion and growing trees might be not quite effective in reducing erosion in highly erodible regions. For the GC watershed where SEM has been persistently low, increases of forest area could further reduce soil erosion. However, it should be pointed out that rill and inter-rill are not the predominant erosion processes on the steeply sloped forested areas, and therefore a USLE-based model such as the WaTEM/SEDEM may not be the most appropriate model in order to precisely map erosion patterns within this particular environment. Nevertheless, increasing the proportion of forest land has been proven to be a very effective large-scale policy measure in order to decrease soil erosion, such as in the context of the Grain for Green program.
4.2 Impacts of topographic factor
Topography has been regarded as another important factor influencing soil erosion (e.g. Baggaley and Potts, 2017; Fu et al., 2011; Prosdocimi et al., 2016; Sun et al., 2014). The inter-watershed differences in soil erosion are at least partly attributable to the topographic differences across the DDK, CPT and GC watershed. As shown by the highest average SEM, soil erosion was the most intense in the CTP watershed (Figure 5b), which was probably due to the steepest topographies (i.e. lowest percentage of the slopes <5º, highest percentage of the slopes between 5º and 15º and highest average slope gradient) in this watershed (Table 1). In contrast, the percentage of the slopes <5º is the highest and the percentages of the slopes between 5º and 15º, the percentages of the slopes >15º and average slope gradient are the lowest in the GC watershed (Table 1). As a result, soil erosion was the slightest and SEM was the lowest in this watershed (Figure 5b). SDR was higher for the GC watershed than those for the CPT and DDK watershed (Figure 5c), which could be explained by the highest average profile curvature, average plan curvature and drainage density, and lowest average surface roughness, in the GC watershed (Table 1). The overall good connectivity with the stream network may have facilitated transportations of the sediments, leading to the higher SDR for the GC watershed as for several other agricultural catchments (Blake et al., 2012).
Figure 9 shows the SEM and dry farmland ratio for three slope classes (i.e. <5°, 5–15° and >15°). These figures indicate clear differences among the three watersheds, both in SEM and dry farmland distribution across the different slope classes. Only 32% of the area of dry farmland in the CPT watershed (Figure 9b) was located on the slopes < 5°, which is 15.4% and 49.2% lower than that in the DDK and GC watershed (Figure 9a and c). Moreover, the ratio of dry farmland area on the slopes between 5° and 15° in the CPT watershed was higher than those in the DDK and GC watersheds (Figure 9). Although the ratio of dry farmland area with slopes >15° in the DDK and CPT watershed was approximately equal (circa 14%), the SEM in the CPT watershed is much higher compared to the DDK watershed, i.e. 83.65 t/ha versus ca. 49.74 t/ha (Figure 9a and 9b). This is most probably because most parts of the CPT watershed have better connectivity than in the DDK watershed, as indicated by the higher SDR for the CPT watershed (Figure 5c). Furthermore, this may demonstrate that a more inappropriate spatial pattern of dry farmland distribution could trigger more serious soil erosion in the CPT watershed. Thus, agricultural activity, and in particular the spatial distribution of related land use patterns, should be strictly controlled in such hilly watersheds. In contrast, the GC watershed had the lowest SEM across all slope classes, as most of the dry farmland is located on relatively flat landscape positions (Figure 9c). However, it should be noted that the average profile curvature, average plan curvature and particularly drainage density in the GC watershed are the highest and average surface roughness is the lowest among these watersheds (Table 1), which results in an overall good connectivity with the stream network. Consequently, the soils eroded by rainfall and tillage were more likely to be transported outside of the watershed, as described by high SDR values in the GC watershed (Figure 5). Therefore, local authorities should recognize and apply the most appropriate and effective soil conservation measures despite the fact that the total sediment yields are relatively low. For example, increasing the diversity of land use and/or landscape could be an effective measure (Hou et al., 2016; Zhang et al., 2017). Some studies have reported that land use diversity has a clear impact on sediment connectivity by affecting the link between the sediments produced on the hillslopes and the material transiting into the river (Prosdocimi et al., 2016; Samec et al., 2018; Zhang et al., 2017). In addition, other measures aimed at reducing sediment connectivity (e.g. the increase of land use and/or landscape diversity and construction of water conservancy facilities such as reservoirs and dams) will also act as a physical barrier and prevent the sediments from reaching the water bodies when they are located on the main runoff/sediment flow pathways.

Soil erosion modulus (SEM) (box plot) metric tons per hectare (t/ha) and dry farmland ratio (line graph) for slopes of different gradients in the (a) Dadukou (DDK), (b) Caopingtou (CPT) and (c) Gaoche (GC) watersheds.
4.3 Possible impacts of subterranean soil loss
On the Guizhou Plateau where karst landscapes have widely developed, some of the eroded soils may be discharged into and partly or totally trapped in caves and underground channels. As a result, these sediments will not reach outlets of watersheds promptly and not be measured at the associated hydrological stations. Therefore, the sediment yield data used in calibrating and validating the WaTEM/SEDEM model may not have been complete with regard to the subterranean sediment transport pathways. This kind of subterranean soil loss has been often mentioned and regarded as a common phenomenon on the Guizhou Plateau, so it is likely that the WaTEM/SEDEM model calibrated has underestimated the total soil loss from each of the three watersheds. However, it is difficult to have an idea of the extent of this underestimation, as specific percentages of the subterranean loss in the total soil loss have been rarely reported so far (e.g. Ma and Zhang, 2018; Xiong et al., 2012). The percentage of the subterranean loss may differ greatly from basin to basin, depending on the levels of karst developments. Geissen et al. (2007) related this to land use types, soil types, geological formations and particularly the presence/absence of sinkholes when considering a study site in tropical Mexico. Field observations carried out by Li et al. (2012) in a small watershed in the central part of the Guizhou Plateau indicated that the percentage of the subterranean loss was extremely low, i.e. accounting for only 0.81% of the total soil loss. Peng et al. (2016) related the extents of karst developments to percentages of areas of exposed bedrocks and suggested that percentages of the underground loss will be low if the percentages of the areas of exposed bedrocks are below 20%. There were tiny patches of exposed bedrocks scattered among the barren and construction land in the three watersheds. The percentage of the barren and construction land ranged from 0.53 to 2.55% with an average of 1.11% in the DDK watershed and from 0.25 to 0.51% with an average of 0.41% in the CPT watershed (Table 2). In the GC watershed, the percentage was slightly higher but still within the range of 2.55–8.89%, with an average of 5.17% (Table 2). The percentages of the barren and construction land were generally low and those of the areas of the exposed bedrocks incorporated in them must be even lower in the three watersheds. Therefore, percentages of subterranean loss could be low overall in these watersheds and the associated underestimation resulting from the subterranean loss is believed to be minimal.
Subterranean soil losses are not considered in any predictive model for soil erosion, for example RUSEL, USLE, WEPP, SIDASS, etc. (De la Rosa et al., 2005; Geissen et al., 2007; Nearing et al., 1989; Wischmeier and Smith, 1978); therefore, there is still the need for more research in this field, with a particular focus on developing erosion models that allow the researcher to take into account site specific subterranean sediment pathways within a given karst-catchment. Nevertheless, the primary task in such studies would be the spatial identification of these pathways, which remains a key methodological hurdle until present. Hence, although we recognize the topographical complexity of the karst regions and the implications that this may have on spatial erosion patterns, we are convinced that WaTEM/SEDEM is the preferred model to be used in this particular context (given the current state of knowledge).
V Conclusions
The WaTEM/SEDEM model was fairly well calibrated and validated for the three watersheds. Good fits between the predicted and observed sediment yield suggest that the WaTEM/SEDEM model is an effective tool to estimate soil erosion for such watersheds on the carbonate-rock dominated Guizhou Plateau. The trends of temporal changes in soil erosion were the same in the three watersheds from 1985–2014, i.e. SY and SEM decreased firstly and then increased, and SDR were almost unchanged. Dry farmland was the most important source of the sediments. Increased rainfall played an important role in intensifying soil erosion in each of the watersheds. SY was the highest from the largest watershed (DDK). SEM was the highest in the CPT watershed due to steep topographies and a high ratio of dry farmland area on relatively steep slopes. In contrast, SEM was the lowest in the GC watershed where terrain was generally gentle and dry farmland was mainly distributed on the relatively gentle slopes. SDR was the highest for the GC watershed, in which the average profile curvature, average plan curvature and drainage density are the highest and average surface roughness is the lowest. Variations in precipitation have influenced SEM in the three watersheds. In particular, SEM was remarkably sensitive to variations in precipitation in the CPT watershed where erosion has been persistently intense. Changes in ratio of dry farmland have also had impacts on SEM in the DDK watershed and particularly in the CPT watershed, but not in the GC watershed, in which SEM was the lowest. Changes in ratio of forestland have had no noticeable impacts on SEM in the CPT and DDK watersheds, but strongly influenced the SEM in the GC watershed. Since the responses of erosion in the different watersheds to changes in the areas of dry farmland and forest areas were different, soil conservation measures adopted should suit local conditions in each of the watersheds.
Due to the deficiency of the sediment data used in calibrating and validating the model and high topographical complexity and heterogeneity in other erosion-determining factors within each of the watersheds, it would be more appropriate to consider the simulation results presented in this study as a preliminary estimate of soil erosion and sediment contribution in this type of highly diverse karst landscape.
Footnotes
Author’s Note
Jeroen Meersmans is now affiliated to TERRA Teaching and Research Centre, Gembloux Agro-Bio Tech, University of Liège, Gembloux, Belgium.
Acknowledgements
Yao Luo is grateful to University of Exeter and NERC for offering him a NERC-NEWTON Research Scholarship to support his stay in UK. We would like to thank three anonymous reviewers whose comments have improved the manuscript.
Declaration of conflicting interests
The authors declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.
Funding
The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This study is granted by the National Natural Science Foundation of China (No. 41571130044) and the Natural Environmental Research Council of the UK and the Newton Foundation (No. NE/N007603/1, NE/S009175/1 and NE/S009116/1).
