Forestation does not necessarily reduce soil erosion in a karst watershed in southwestern China

I Introduction

Forests play a key role in preventing soil erosion by moderating the soil and hydrological ecosystem and water balance. The water and soil conservation functions of forests have long been recognized (Altieri et al., 2018; Dong et al., 2020; Zhang et al., 2015). In areas dominated by hydraulic erosion, natural forestland reduces soil erosion mainly in the following ways: a forest canopy reduces rainfall kinetic energy and sediment concentration; surface litter absorbs rainfall and increases soil infiltration by providing organic matter; and roots can consolidate and hold soil (Borrelli et al., 2017; Zhang et al., 2011, 2015). However, it has also been indicated that there is a negative correlation between soil erosion and forest cover in a certain range on a hillslope scale (Miura et al., 2015; Sun et al., 2010). The effect of forestation on preventing soil erosion could significantly vary in different areas (Altieri et al., 2018; Li et al., 2011; Miura et al., 2015). Therefore, more research on the effect of forestation on soil erosion control and its sustainability on different areas with varied climate, soil and vegetation are urgently required.

Soil erosion is known to be extremely serious in karst regions around the world – for example, in western Texas in the US (Cao et al., 2009; Rodda et al., 2010), New Zealand (Wilcox et al., 2002) and in southwestern China (Li et al., 2011). Located in the center of the karst region in southwestern China, Guizhou is largely covered by carbonate rocks with karst landscapes (Jie et al., 2012). The carbonate rocks in Guizhou Province cover an area of 130,000 km², accounting for 73% of the province’s land area (Wan, 2003). As a result of long-term unsustainable human activities, such as extensive logging and steep-slope farming from the 1950s to 1990s (Chicas and Omine, 2015; Silva-Sánchez et al., 2014; Wan, 2003; Zha et al., 1992), the forest fractional cover of the province decreased from 45% before the 1950s to 12.6% in the 1990s (Han, 2006). Correspondingly, the area of soil erosion increased from approximately 14% to 50% (Wan, 2003). The soil erosion area in Guizhou Province reached 76,700 km², and the sediment flowing into the Yangtze River and Pearl River was 2.7×10⁸ t in 1995, eroding the topsoil layer by 0.3 cm per year by the end of last century (Wan, 2003). Substantial topsoil loss was caused by erosion, leading to rocky desertification (Cao et al., 2005; Tong et al., 2018; Zhang et al., 2006). The Returning Farmland to Forest Project (RFFP) was launched in 1999 across China, with one of the aims being to reduce soil erosion by increasing forest cover percentages (Wang et al., 2011). Since 2000, the implementation of the RFFP has significantly increased forest cover in southwestern China’s karst areas (Chen et al., 2019; Tong et al., 2018) and has led to improved ecological restoration in the karst rocky desertification areas (Tong et al., 2014). Recently, the great potential influence of bedrock geochemistry on plant growth in karst areas by controlling the regolith water-holding capacity was proposed (Jiang et al., 2020), implying that particular attention should be given to the lithological impacts on sustainable growth of forests under the local climate after forestation, as well as its effect on preventing soil erosion in karst areas. Therefore, Guizhou Province is an ideal place to explore whether forestation can prevent soil erosion in karst areas.

Numerous studies have investigated the soil erodibility and erosivity, and different models have been developed with varied research objects and purposes (Batista et al., 2019; Dutta, 2016), including empirical prediction models such as the Revised Universal Soil Loss Equation (RUSLE) (Renard and Ferreira, 1993) and physical process models such as the Water Erosion Prediction Model (Nearing et al., 1989). Process-based models require calibration of different parameters, such as runoff and sediment delivery, which is often limited by the scarcity of data resources. In contrast, empirical prediction models are relatively intuitive and simple, and can be adjusted according to the specific region of the factor (Gao et al., 2014). A typical example of the empirical prediction models is the RUSLE, which has been widely used to depict the soil erosion rates and grades of erosion risk at various spatial scales (e.g. runoff plot scale, watershed scale and basin scale) in the karst region of southwestern China (Feng et al., 2016; Xu et al., 2008; Zeng et al., 2017). Based on the calculation of annual soil erosion rates through RUSLE, a spatial distributed model called the Water and Tillage Erosion Model and Sediment Delivery Model (WaTEM/SEDEM) was developed (Oost et al., 2000; Rompaey et al., 2001). It considered the routing of the eroded sediment to the river channel network and the transport capacity of each spatial unit, and, thus, can simulate the soil redistribution process (Rompaey et al., 2001). With this model, the influence and contribution of each factor on the soil erodibility and erosivity can be evaluated within the parameters of the model. Moreover, it has few parameters to adjust, a relatively simple structure and high simulation efficiency, which have led to its wide application in soil erosion simulations (Alatorre et al., 2010; Haregeweyn et al., 2012; Sheng and Fang, 2014). To date, the applicability of the WaTEM/SEDEM model has been assessed in several different areas and has further been used to analyze the driving forces of soil erosion, sedimentological connectivity and other attributes for drainage basins (Alatorre et al., 2010; Liu and Fu, 2016; Sheng and Fang, 2014; Sheng et al., 2015). Moreover, it has been applied in a karst watershed in Guizhou by evaluating the model with the dating data of the reservoir sediments, showing the suitability of the model in karst areas (Qiu et al., 2019). The driving forces of soil erosion in different catchments have been widely studied with WaTEM/SEDEM simulation (Fang, 2017; Fang and Sun, 2017; Feng et al., 2010), but few were conducted in karst regions. In the simulation of WaTEM/SEDEM, the driving forces of soil erosion include dynamic factors (rainfall and land use) and static factors (terrain and soil). For each grid, the impact of rainfall on soil erosion was expressed as R factor (the rainfall erosion factor), and the influence of land use on soil erosion was assigned into C factor (the vegetation cover and crop management factor) and P factor (the erosion control practice factor mainly influenced by the tillage measures); thus, the change of rainfall and land use and its effect on soil erosion can be quantitatively captured and analyzed.

In this study, we took the Chong’an River, a small watershed in Guizhou Province of China, as a case study. Based on the observed annual sediment yield at the outlet of the basin and WaTEM/SEDEM simulation, we estimated the soil erosion distribution of the basin for the period 1986–2018. To clarify the causes of the temporal fluctuation of the sediment yield and effectiveness of large-scale forestation on the reduction of soil erosion, the spatially distributed soil erosion results were used to do the temporal correlation analysis with the dynamic factors of soil erosion on independent pixels. Changes in vegetation cover were interpreted from satellite image. We aimed to answer the following scientific questions: (i) Can forestation significantly reduce soil erosion in this karst basin? (ii) What is the most related dynamic factor influencing erosion changes in the basin during the study period? (iii) Can lithology affect vegetation restoration and further soil erosion?

II Materials and methods

1 Study area

The Chong’an River is a major tributary flowing into the Qingshui River in eastern Guizhou. The Qingshui River drains to the Yuanjiang, which is a major tributary of the Yangtze River. The Chong’an River is a rain-fed river in a mountainous area. The river length is 141 km, the elevation is 1190 m for the source and 546 m for the water outlet, and the drop is 644 m (Song and Luo, 2015). Originating from the hills, it flows from southwest to northeast in the drainage basin, and pours into Qingshui River in Guizhou, with a drainage basin spanning 107°35’–108°27’E, 26°40’–27°47’N (Figure 1).

OPEN IN VIEWER

Figure 1.

Location and topography of the Chong’an River Basin and the distribution of the interpretation training samples. The topography data of the study area were obtained from the Advanced Spaceborne Thermal Emission and Refection Radiometer Global Digital Elevation Model (http://earthexplorer.usgs.gov) (Tachikawa et al., 2011), with a spatial resolution of 30 m.

In the basin, the majority type of physiognomy is dissolution residual hills and low hills composed of carbonate rocks, followed by the middle and low mountain landforms of erosion and denudation with Proterozoic rocks (Song and Luo, 2015). The basin is located in the subtropical temperate monsoon climate zone. The average annual temperature is 14.7°C, the annual precipitation is 1278 mm and precipitation is greater in mountainous areas than in hilly areas (Song and Luo, 2015). The major soil exposed in the basin is lime soil. The basin is mainly covered by mixed deciduous forest with subtropical evergreen broad-leaved forests and Masson pine forest (Song and Luo, 2015). According to the observed sediment yield of Wanshui hydrological station, which is located at the outlet of the Chong’an River Basin, the average annual sediment yield was 627,011 t. Based on the location of the Wanshui hydrological station and the terrain data, the Chong’an River basin area was extracted, which was approximately 2562.1 km² (Figure 1), and the average slope of the basin is 16.7°.

3 Random forests classification

Random forests classification is a multi-decision tree classification method to integrate multiple trees through the ideal of integrated learning (Breiman, 2001). It has been widely used for its robust performance and generalization ability; in particular, it is highly efficient in classification-based multi-feature (Feng et al., 2019). The Landsat images used for classification were preprocessed by radiation calibration, atmospheric correction of the fast line-of-sight atmospheric analysis of spectral hypercubes method, terrain correction and cloud removal (Zhu and Woodcock, 2012). The regions of interest of bare land, cropland, forestland, grassland, residential land, shrubland and waterbody made on the Google Earth images in 2018 were taken as training samples (Figure 1). A total of 23 classification features (Table 1) were extracted based on the preprocessed Landsat images and the topographic data; thus, the classification model was constructed based on these 23 features of the training samples. To eliminate the influence of sensor differences (between Landsat 5 TM and Landsat 8 OLI) on model accuracy, the classification features extracted from the Landsat images were normalized. Additionally, the test samples were developed based on the Google Earth images in the study area of 2002 and 2018, respectively, so as to detect the classification accuracy of the model for different sensors (TM and OLI). Thus, the land use maps of the basin in 1986, 1994, 2002, 2007, 2014 and 2018 were made. These maps were used to reconstruct the change in forestland in the past 33 years and to derive the input factors of the WaTEM/SEDEM model.

OPEN IN VIEWER

Table 1.

Features used for random forests classification.

Feature type	Features
Spectral features	Band_B
	Band_G
	Band_R
	Band_NIR
	Band_SWIR1
	Band_SWIR2
Index features	NDVI ( ρ N I R − ρ R E D ) / ( ρ N I R + ρ R E D )
	NDBI ( ρ M I R − ρ N I R ) / ( ρ M I R + ρ N I R )
	NDWI ( ρ G R E E N − ρ N I R ) / ( ρ G R E E N + ρ N I R )
Topographic feature	DEM
	Slope
Texture features	Mean_G
	Mean_R
	Mean_NIR
	Variance_G
	Variance_R
	Variance_NIR
	Homogeneity_G
	Homogeneity _R
	Homogeneity _NIR
	Contrast_G
	Contrast_R
	Contrast_NIR

4 WaTEM/SEDEM model

4.1 Model description

The spatially distributed WaTEM/SEDEM model is developed on the basis of RUSLE, WaTEM and SEDEM (Renard and Ferreira, 1993). The WaTEM/SEDEM model takes into account the interception of soil erosion by land use patterns and sediment transport processes and can be used to simulate the average annual soil erosion rate and sediment intensity at the pixel scale (Sheng and Fang, 2014; Vente et al., 2008). The model contains three parts: an annual soil erosion module and an annual sediment transport capacity module based on the pixel scale and sediment flow algorithm module.

The formula for calculating the annual soil erosion was adopted from the RUSLE model (Renard and Ferreira, 1993):

A = R * K * L S 2 D * C * P

1

where A is the average amount of soil loss caused by gully erosion (t·ha^-1·yr^-1), R is the rainfall erosion factor (MJ·mm·ha^-1·yr^-1), K is the soil erosion factor (t·h·MJ^-1·mm^-1) and LS_2D is the topographical slope and length factor. C is the vegetation cover and crop management factor (dimensionless), and P is the erosion control factor (dimensionless). The specific calculation methods for each factor are shown in the supplementary materials.

WaTEM/SEDEM deals with the sediment transport capacity (TC) by calculating the maximum sediment amount that can pass a pixel. We calculated TC in this study according to the following equation (Oost et al., 2000; Desmet and Govers, 1995):

T C = K tc * R * K * ( L S 2 D − 4.12 * s * 0.8 )

2

Where s represents the slope, and Ktc is the transport capacity coefficient that needs to be calibrated. Based on the soil erosion and sediment transport capacity of the pixel unit, the multiple confluence algorithm (Desmet and Govers, 1995; Oost et al., 2000) was used to simulate the process of sediment transported into the river from high to low according to the change in confluence direction.

4.2 Calibration and validation

Model calibration aims to make the model suitable for the study area through adjusting the key parameters. In WaTEM/SEDEM, by adjusting Ktc_min and Ktc_max (two transport capacity coefficient parameters) (Alatorre and Beguería, 2009; Alatorre et al., 2010), and comparing the sediment yield data observed at the outlet with simulation results of the model, it is possible to optimize accuracy and stability of model simulation. Following the different abilities of land use types to transport sediment by overland flow, Ktc_min was used for forestland, grassland and shrub land, while Ktc_max was used for cropland and bare land. Due to the regional scale difference and climate difference, the optimal Ktc parameter differed greatly in different regions; for example, the Ktc_min varied from 0.1 to 70 (Alatorre and Beguería, 2009; Fang, 2017; Sheng and Fang, 2014), but the ratio between the Ktc_max and Ktc_min was mostly between 1.0 and 3.89 (Keesstra et al., 2009; Verstraeten and Poesen, 2001; Verstraeten et al., 2001). In the calibration results, the ratio was 3.80, which fell within the valid range. The indexes used to evaluate the model simulation efficiency and accuracy were the Nash–Sutcliffe efficiency coefficient (NSE) (Nash and Sutcliffe, 1970) and relative root mean square error (RRMSE), respectively:

N S E = 1 − ∑ i = 1 n ( O i − P i ) 2 ∑ i = 1 n ( O i − O m e a n ) 2

3

R R M S E = 1 n ∑ i = 1 n ( O i − P i ) 2 1 n ∑ i = 1 n O i

4

where O_i is the observed value, P_i is the predicted value, O_mean is the mean observed value and n is the number of observations. NSE can range from –∞ to 0. The closer the value of NSE is to 1, the better the simulation performance of the model. Since the observed sediment transport rate data of Wanshui station only included five years among the simulated years, the annual sediment transport rate data observed at the Wanshui station of three years (2007, 2014, 2018) were selected for calibration, and the other two years (1986, 2002) were selected for validation.

III Results

1 Accuracy test and validation

1.1 Accuracy test of interpretation results based on random forests

The land use maps of the seven classes (supplementary Table S1) were obtained. Since the training features were extracted from the images of OLI (2014 and 2018) and TM (1986, 1994, 2002, 2007) that are two different sensors, we used samples of 2018 to test the interpretation results of OLI and used samples of 2002 to test the interpretation results of TM. Finally, the overall accuracies of 2002 and 2018 were 71% and 73%, respectively. All interpretation results were converted into the input factors required by the WaTEM/SEDEM and used for the sediment yield simulation.

1.2 Evaluation of WaTEM/SEDEM

The model calibration results are shown in Figure 2(a), with a comparison between the observed and simulated sediment transport rates at Wanshui station in 2007, 2014 and 2018. The results showed that when the Ktc_min value was 10, the Ktc_max value was 38, the NSE reached a maximum (0.84) and the RRMSE reached a minimum (0.16), which suggested that the model simulation had obtained the best simulation performance and the best accuracy. With these two parameters set, the comparison of the observed values with the simulated values of all five years is shown in Figure 2(b). Through the model with the two optimal parameter combinations, the sediment yield of the basin in 1994, which lacks observation data, was simulated.

OPEN IN VIEWER

Figure 2.

(a) Convergence graphs of NSE and RRMSE of the model according the observed sediment yield data of 2007, 2014 and 2018. (b) Comparison of observed sediment yield at the outlet of the basin and simulated sediment yield – the year with “*” represents the data used for calibration, and the remaining years without “*” represent the data used for validation.

2 Relationship between forestland and sediment yield in the basin

During 1986–2018, the most related dynamic factor for the change in the sediment yield on the pixel scale was mainly the C factor. Based on the temporal correlation analysis between factor C (the vegetation cover and crop management factor), factor P (the soil erosion control practice factor) and factor R (the rainfall erosion factor) with the soil erosion simulation results of 1986–2018, the most related dynamic factor of each pixel is shown in Figure 3. The “NaN” value in the figure represents the absence of data or the absence of factors that had a significant correlation with variation of sediment yield on the pixel scale. According to the statistical proportions of dynamic factors most correlated with sediment yield, factor C was the factor with the strongest correlation with sediment yield change among the three dynamic factors of the model, accounting for 43% of pixels of the whole basin. The pixels of with the R factor or P factor as the most related dynamic factor to erosion accounted for 9% and 7%, respectively. For the whole basin, pixels with significant temporal correlation between the variation of sediment yield and the dynamic factors of the model accounted for 59%. Moreover, by comparing forestland proportional changes in the basin with the sediment yield data observed in the basin (Figure 4), it was found that the varied proportions of forest cover in the basin showed a significantly negative correlation with annual sediment yield.

OPEN IN VIEWER

Figure 3.

The distribution of dynamic factors that were mostly related with the variation of sediment yield in the Chong’an River Basin from 1986 to 2018, where “NaN” means the absence of data or there was no significant correlation between the variation of sediment yield with factor C, factor R or factor P.

OPEN IN VIEWER

Figure 4.

Annual proportion of forestland (according to the interpretation results) and the change in the annual sediment yield of the basin (observed at the outlet of the basin). The blue dot represents the year without observed data, represented by the simulation result of WaTEM/SEDEM.

3 Change in forest cover on different lithologies

A significant difference occurred with the reduction in forestland on dolomite and limestone from 2007 to 2014. According to the land cover interpretation results from satellite imagery, the forest cover in the basin was divided into four significant stages based on the change (Figure 5(a)). From 1986 to 2002, the forest cover decreased by 16% after 16 years. From 2002 to 2007, it increased sharply by 12%. However, it decreased by another 15% from 2007 to 2014 and gradually recovered after 2014. To further investigate the growth differences in trees on dolomite and limestone in the karst area, the areas where limestone and dolomite were separately distributed (Ye et al., 2017), which accounted for 45% of the whole basin (Figure 5(b)), were used for the following analysis. The annual change rates of the forest cover at different stages in the separate distribution areas of dolomite and limestone were calculated and are shown in Figure 5(c). There was no significant difference between the annual change in forest cover in the separate dolomite and limestone distributed areas between 1986 and 2007 (p > 0.05). In 2002–2014, the change in forest cover showed a significant difference in the separate dolomite and limestone areas (p < 0.05). It indicated a significant difference between the forest covers in the dolomite and limestone areas during 2007–2014. Specifically, the annual decrease rate of forestland on dolomite was 0.77% higher than that on limestone in 2007–2014 (Figure 5(c), (d)).

OPEN IN VIEWER

Figure 5.

(a) Changes in forestland at different stages in the Chong’an River Basin from 1986 to 2018. (b) Annual change rates of forestland in the dolomite distribution area and limestone distribution area at different stages in 1986–2014. (c) The lithology distribution map of the Chong’an River Basin (Ye et al., 2017). (d) Tree change mechanism on different lithologies after severe drought in karst areas.

IV Discussion

Our results showed that the change in forest cover influenced the variation in the sediment yield of the basin most. Factor C reflects the comprehensive inhibition of soil erosion by vegetation and crop management factors and the protection effect of the surface. The value of the factor is mainly determined by the specific vegetation cover and crop management measures (Lu et al., 2011). During 1986–2018, the proportions of grassland and shrubland in the basin were less than 3%; therefore, the impact on vegetation cover can be ignored. Previous studies in karst areas of other parts of the world have illustrated the effects of forest cover changes on runoff soil erosion through the yield observation and model simulation. In New Zealand, the effect of increase of forest cover on reducing runoff and soil erosion was illustrated by the soil erosion models under different land use scenarios in the catchments (Cao et al., 2009; Rodda et al., 2010). In karst regions of Texas, the juniper removal could increase streamflow through reducing rainfall interception and infiltration (Wilcox, 2002). Moreover, several runoff models have been implemented to illustrate the increase in surface runoff from controlling or clearing shrubs, and thus increased the risk of soil erosion in karst regions of Edwards Plateau, Texas (Afinowicz et al., 2005; Brown and Raines, 2002; Wu et al., 2001). The effect of reduction in forest cover on increase of sediment yield was also confirmed by field observation in karst catchments in western Texas (Huang, 2006).

Forestation and severe drought resulted in obvious stage characteristics for the change in forest cover in the basin. Since 2002, when the implementation of the RFFP was launched, large-scale forestation activities, such as returning farmland to forests and grasslands, have started, leading to an increase in foreland coverage (Chen et al., 2019; Tong et al., 2018; Zhang et al., 2017b). However, from 2007 to 2014, Guizhou experienced frequent and severe droughts: severe autumn–winter–spring droughts from autumn 2009 to spring 2010 (Wang et al., 2010; Zhang and Kovacs, 2012a) and extreme summer droughts in 2011 and 2013 (Hu et al., 2019; Zhang et al., 2012b), resulting in a sharp reduction in forestland in a short time. It had been reported that drought could hamper karst vegetation recovery of ecological projects (Zhang et al., 2017a). Since drought could lead to decreased vegetation productivity, which, in turn, could increase sediment yield on hillslope-scale in karst region (Weltz and Spaeth, 2012). Therefore, the sediment yield increased during the period with the forestland reduction due to severe drought after forestation in our study area. Our study implied that the effect of forestation on preventing soil erosion is strongly influenced by the sustainability of the forest after forestation, while drought is an important external factor for the growth of the forest after forestation.

Large-scale forestation has led to significantly different variations in the reduction rates of forestland on dolomite and limestone under drought conditions. It was found that the main factor that affected the vegetation change in eastern Guizhou was water availability from 1982 to 2011 (Zhang et al., 2017b). The importance of water availability is closely related to the thin soil layer, uneven soil distribution and low water-holding capacity in this area. In the study area, limestone and dolomite are widely distributed as the main lithologies, interbedded with sandstone, shale, mudstone and other rocks. The rock crevice structure in the limestone area was better developed than that in the dolomite area, which was more suitable for tree growth through a more effective water supply in karst region (Liu et al., 2019). The influence of bedrock geochemistry on vegetation change in karst areas was mainly through controlling the regolith water-holding capacity on vegetation growth (Jiang et al., 2020). Therefore, lithologies with different crevice developments and geochemical elements are particularly important to vegetation growth in this region in that they determine the structure and distribution of the soil. It suggested that the risk of tree mortality during severe drought is greater in the dolomite area than that in the limestone area due to the lack of an effective water supply in the root system (Figure 5(d)). We suggested that forestation should not only consider climate condition but also the plant growth potential affected by the lithology in karst areas. Also, more quantitative analysis is in need about the specific tree growth influenced by the climate and lithology in karst areas.

V Conclusions

Through the remote sensing interpretation by random forests and the simulations of soil erosion by WaTEM/SEDEM in the Chong’an River basin, we obtained the distribution of forest cover and soil erosion in six time nodes of 1986–2018. According to the temporal correlation analysis between the dynamic factors about the rainfall and land use change and simulated erosion results on the pixel scale, we found that in the past 33 years, the change in forest cover was the most related factor affecting the drastic fluctuation in sediment yield. Due to large-scale forestation and severe droughts, the forestland varied sharply in 2002–2014. Consequently, soil erosion has fluctuated correspondingly and reversely in the basin with the forest cover change due to forestation and severe drought. Since rock crevices are more developed in limestone than in dolomite, when serious drought occurs after large-scale forestation, the annual rate of decrease in forest cover in dolomite is significantly higher than that in limestone, weakening the role of forestation in soil erosion control. Our study suggested that the effect of forestation on preventing soil erosion were influenced by the sustainability of forests influenced by both lithology and climatic drought in karst areas.

Supplemental material