Assessing the Impact of Terraces and Vegetation on Runoff and Sediment Routing Using the Time-Area Method in the Chinese Loess Plateau

Abstract

Terracing and vegetation are an effective practice for soil and water conservation on sloped terrain. They can significantly reduce the sediment yield from the surface area, as well as intercept the sediment yield from upstream. However, most hydrological models mainly simulate the effect of the terraces and vegetation on water and sediment reduction from themselves, without considering their roles in the routing process, and thus likely underestimate their runoff and sediment reduction effect. This study added the impact of terraces and vegetation practice on water and sediment routing using the time-area method. The outflow in each travel time zone was revised in each time step by extracting the watershed of the terrace units and the vegetation units and calculating the water or sediment stored by the terraces or held by the vegetation. The revised time-area method was integrated into the Land change Model-Modified Universal Soil Loss Equation (LCM-MUSLE) model. Pianguanhe Basin, in the Chinese Loess Plateau, was chosen as the study area and eight storms in the 1980s and 2010s were selected to calibrate and verify the original LCM-MUSLE model and its revised version. The results showed that the original model was not applicable in more recent years, since the surface was changed significantly as a result of revegetation and slope terracing, while the accuracy improved significantly when using the revised version. For the three events in the 2010s, the average runoff reduction rate in routing process was 51.02% for vegetation, 26.65% for terraces, and 71.86% for both terraces and vegetation. The average sediment reduction rate in routing process was 32.22% for vegetation, 24.52% for terraces, and 53.85% for both terraces and vegetation. This study provides a generalized method to quantitatively assess the impact of terraces and vegetation practice on runoff and sediment reduction at the catchment scale.

1. Introduction

Soil erosion and water loss is a serious environmental problem in the Chinese Loess Plateau [1,2,3]. It can cause soil deterioration and loss of sustainable productivity in croplands [2,4,5]. In addition, sediment yields and chemical loadings associated with soil erosion can cause severe degradation of surface water quality [6,7,8]. To control soil and water losses, terrace engineering was implemented since the 1980s and the Grain for Green (GFG) project was launched in 1999 [9,10,11,12]. Therefore, it is necessary to quantitatively analyze the effect of terraces and vegetation on runoff and sediment in the Loess Plateau.

Slope terracing and vegetation planting are common practices for soil and water conservation on sloped terrain susceptible to water erosion, and have been proven to be effective at retaining water and soil [6,13,14,15,16,17,18,19,20]. Terraces and vegetation can reduce the sediment yield from themselves significantly. Terraces reduce the peak runoff rate by reducing the slope gradient and slope length of the hillside [21,22]. Owing to the topographic slope and the embankment, terraced fields have a certain storage capacity similar to a reservoir, which can intercept surface runoff and sediment and promote runoff infiltration, evaporation, and sediment deposition [23,24,25]. The vegetation canopy can intercept rainfall directly and influence the rainfall kinetic energy and erosion rates. The stems, roots, and litter layer of vegetation can reduce runoff discharge by promoting infiltration, increasing surface roughness, and slowing down the overland flow and peak runoff. Vegetation also can reduce soil erosion by reducing surface flow volume and increasing sediment trapping through reducing flow velocity [21,26,27]. In addition to this, both terraces and vegetation can intercept the sediment yield from upstream and achieve sediment reduction in the valley by reducing runoff flowing from the slope into the valley [16,24,28,29].

Numerous studies have provided many insights into how terraces and vegetation control water erosion at local scales using observational experiments [20,27]. Terraces could be classified into level terraces, slope terraces, slope-separated terraces, and zig terraces according to their structure [23]. In the Loess Plateau of China, the main terrace type is level terraces. Yao [30] found a terraced field could reduce soil erosion by 92–100% compared with sloped farmland, while Wu [31] found the average benefit of level terraces on soil and water conservation were 86.7% and 87.7%, respectively. Huo and Zhu [32] combined soil water moisture and 137Cs content analysis and found the average soil and water conservation benefit of level terraces was 53%_ENREF_18, while Pan and Shangguan [27] reported that grassplots had 14–25% less runoff and 81–95% less sediment yield compared to a bare soil plot. Meng et al. [33] conducted a series of laboratory flume simulation experiments and the results showed that vegetation could reduce the mean velocity by 31–65%_ENREF_36. As the effectiveness of terraces and vegetation is limited by many factors, such as climate, soil properties, topography, land use, vegetation type, and spatial patterns, the diversity and natural variability of previously conducted erosion studies limit their potential extrapolation to the catchment scale [23,34,35,36]. Thus, how to assess the effects of terraces and vegetation on water erosion control at the catchment scale remains a crucial issue.

Estimation of rainstorm-generated sediment yield by means of a hydrological model is an important way to quantitatively evaluate the effect of soil and water conservation measures, such as by using the Soil and Water Assessment Tool (SWAT) model [37,38,39], the Agricultural Non-Point Source pollution (AGNPs) model [40], and the Agricultural Policy/ Environmental eXtender (APEX) model [41]. In these models, accounting for the impact of terraces and vegetation on runoff and sediment yields has focused on reduction from themselves through adjusting the key input variables, such as the Soil Conservation Service Curve Number (SCS-CN), slope gradient, slope length, Universal Soil Loss Equation (USLE) support practice factor (P-factor), and cover and management factor (C-factor) [11,21,22,25,29,42,43], without considering the roles of water and sediment reduction in the routing process. Then, the runoff and sediment reduction effect of terraces and vegetation is likely to have been underestimated. Therefore, it is necessary to take into account the water and sediment reduction effect of terraces and vegetation in the model’s routing process.

Explicitly simulating the interaction between conservation practice and watershed response is difficult. The time-area method, with its inherently distributed concept, explains which parts of a watershed contribute to runoff during a specific period [44]. It provides a simple and useful tool to understand runoff mechanisms and is widely used at the catchment scale [44]. It has also been indicated that the routing of sediment through time-area segments in a catchment produces better results than the conventional routing through a network of cells [45,46]. Her and Heatwole [47] revised the time-area method and considered the upstream contribution for routing sediment, such that the new method provided detailed spatial representation_ENREF_15. These studies mainly focused on the effect of surface heterogeneity on routing time or flow velocity [48,49]. Topographic features are delineated through a Digital Elevation Model (DEM), slope gradient, flow direction, and so on. However, the lack of representation of a specific terrace and vegetation process makes it difficult to quantitatively distinguish vegetation and terraced fields and their integrated ability for water and sediment reduction.

Thus, the purpose of this study was to incorporate the effect of terraces and vegetation on runoff and sediment routing in the time-area method and assess their impact in water and sediment reduction. The outflow in each travel time zone was revised in each time step by extracting the watershed of the terrace units and the vegetation units, and calculating the water or sediment stored by the terraces or retained by the vegetation based on their properties.

2. Materials and Methods

2.1. Study Area

In this study, we selected the Pianguanhe River in the key soil erosion region of the Loess Plateau as our case study. The Pianguanhe River is located in the hilly and gully region of the Loess Plateau (111°26′ E–112°17′ E, 39°15′ N–39°42′ N) (Figure 1). It is a tributary of the Yellow River with a watershed area of 2089 km², with altitude ranging from 984 to 2162 m. The study area features a semi-arid continental climate, with an average annual rainfall of 429 mm. The uneven seasonal distribution of precipitation results in more than 80.9% of the annual precipitation occurring from May to September [50]. The average annual runoff and sediment discharge is 39.48 million m³ and 12.58 million t, respectively. The sediment modulus is 65.23 t/ha/yr [50]. The overall loss of soil and water is a serious environmental issue.

2.2. Methods

The serious soil and water loss in the Loess Plateau has always occurred because of the surface runoff caused by heavy rains. In this paper, runoff generation was simulated based on the Land Catchment Model (LCM), which is a flood event forecasting model that was developed through more than 300 artificial rainfall experiments in the Loess Plateau [51,52]. The model has been modified into a distributed model [53,54,55]. The Modified Universal Soil Loss Equation (MUSLE) was chosen to estimate the sediment yield of individual heavy rainfall events [56], and numerous studies had proven its applicability to estimate the sediment yield in the Loess Plateau [57,58,59]. The main functions of LCM-MUSLE model are listed in

Table 1. Functions and parameters related to runoff generation in the Land Catchment Model-Modified Universal Soil Loss Equation (LCM-MUSLE) model.

No.	Module Name	Equations Reference	Equations Reference
1	Canopy interception	$S_v=V_f\times S_{\max }\times [1-{{\displaystyle e}}^{-\eta \frac{P}{S_{\max }}}]$ $S_{\max }=0.935+0.498\times LAI-0.00575\times LA{I}^2$ $\eta =0.046\times LAI$ $P^{\prime }=P-S_v$	[62]
2	Surface runoff	$Q_d=P^{\prime }-f=P^{\prime }-R\times P^{\prime r}$	[51]
3	Interflow	$Q_l=L_a(W_s/W_{sm})f$	[51]
4	Base flow	$Q_b=K_b(G{W}_s+REC)$ $REC=R_c(W_s/W_{sm})(f-Q_l)$	[51]

Annotation: S_v: interception (mm); V_f: vegetation coverage (%); S_max: canopy storage capacity (mm); η: a correction coefficient; P: precipitation (mm); LAI: leaf area index; Q_d: surface runoff (mm); P′: effective rainfall (mm); f: infiltration (mm); R and r: both infiltration coefficient; L_a: interflow coefficient; W_s: unsaturated soil moisture storage (mm); W_sm: maximum soil moisture storage capacity of the soil layer (mm); Q_l: interflow (mm); Q_b: base flow (mm); K_b: base flow coefficient; GW_s: groundwater storage (mm); REC: groundwater recharge (mm); R_c: groundwater recharge coefficient.

and

Table 2. Functions and parameters related to sediment yield in the Land Catchment Model-Modified Universal Soil Loss Equation (LCM-MUSLE) model.

No.	Module Name	Equations Reference	Equations Reference
1	Sediment yield	$Sed=11.8\times {(Rs\times q_{peak}\times A_{pixel})}^{0.56}\times K\times C\times P\times LS\times CFRG$	[63]
2	Runoff factor	$q_{peak}=\frac{{\alpha }_{tc}\times Rs\times A_{pixel}}{3.6\times t_{conc}}$	[37]
3	Soil erodibility factor	$\begin{array}{l}K=\left\{0.2+0.3\exp \left[-0.0256S_d(1-\frac{S_i}{100})\right]\right\}{(\frac{S_i}{C_i+S_i})}^{0.3}\left[1-\frac{0.25C}{C+\exp (3.72-2.95C)}\right]\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\left[1-\frac{0.7(1-S_d)}{1-S_d+\exp (-5.51+22.9(1-S_d))}\right]\end{array}$	[63]
4	Topographic factor	$L={(L_{slp}/22.1)}^{0.44}$ $S=\{\begin{array}{l}10.8\sin \theta +0.03,\hspace{1em}\theta <5°\\ 16.8\sin \theta -0.5,\hspace{1em}5°\le \theta <10°\\ 21.9\sin \theta -0.96,\hspace{1em}\theta >10°\end{array}$	[64,65]
5	Cover and management factor	$C=\{\begin{array}{l}1, c_v=0\\ 0.6508-0.3436\lg c_v, 0<c_v\le 78.3\%\\ 0, c_v>78.3\%\end{array}$	[66]
6	Coarse fragment factor	$CFRG=\exp (-0.053\times rock)$	[63]

Annotation: Sed: sediment yield (t); Rs: surface runoff (mm); q_peak: peak runoff rate (m³/s); A_pixel: area of the grids (ha); K: soil erodibility factor (0.013·t·m²·h/(m³·t·cm)); C: cover and management factor; P: support practice factor; LS: topographic factor; CFRG: coarse fragment factor; α_tc: fraction of rainfall that occurs during the time of concentration (for event modeling the value of α_tc is 1); t_conc: the time of concentration for the grid; S_d: percent sand content (%); S_i: percent silt content (%); C_i: percent clay content (%); C: percent organic carbon content of the layer (%); L_slp: slope length (m); θ: gradient of the slope (°); c_v: vegetation coverage (%); rock: percent rock in the top soil layer (%).

. The reader is referred to Luo et al. [59] for a detailed description of the LCM-MUSLE model, which is not described here. Runoff and sediment generation from sub-catchments was routed to the outlet by considering overland flow and stream channels. Time-area method was used for overland routing [49] and Muskingum method was used for channel routing [60,61].

With the purpose of simulating the influence of terrace and vegetation units on water and sediment reduction at confluences, as the critical hydrologic process, both terrace and vegetation modules were added to the time-area method. In the time-area method, the catchment is divided into a number of travel time zones via isochrones [67]. By extracting the watershed of terrace units and vegetation units, and calculating the water stored by the terraces or intercepted by vegetation, the runoff yield and the outflow in each travel time zone are revised. This represents the spatiotemporally varied flow in the routing simulation. Sediment reduction was achieved in a similar way.

The integrated structure of the model is shown in Figure 2. Specific equations and methods are discussed in the following sections. The integrated model can represent landscape heterogeneity in detail, if a higher spatial resolution DEM is used.

2.2.1. Time-Area Method for Overland Routing

Isochrones are defined as the contours of equal travel time to the catchment outlet, and the travel time zone is the area between two adjacent isochrones [67]. The isochrones of the whole basin are calculated, and then these are adjusted within each sub-catchment, as shown in Figure 3. It is assumed that the area of each time zone is ΔA₁, ΔA₂, ···, ΔA_n, and their corresponding travel time is τ₁, τ₂, ···, τ_n, respectively. Then, we set Δτ = τ_i − τ_i−1 (i = 1, 2, ···, n). It is assumed that m is the time of runoff generation in hours, ΔR_r (r = 1, 2, …, m) is the runoff depth in m, and the runoff discharge of the outlet at time step j is:

For n ≥ m

$$Q_j=\frac{1}{\Delta \tau }{\displaystyle \sum _{r=1}^k\Delta R_r\Delta A_{j-\mathrm{r}+1}},(j=1,2,3,\dots ,n+m-1)$$

(1)

If j < m, k = j; if j ≥ m, k = m, and if j > n, ΔA = 0.

For n < m:

$$Q_j=\frac{1}{\Delta \tau }{\displaystyle \sum _{r=1}^k\Delta A_r\Delta R_{j-r+1}},(j=1,2,3,\dots ,n+m-1)$$

(2)

If j < n, k = j; if j ≥ n, k = n, and if j > m, ΔR = 0.

For travel time zone i, its runoff discharge at time step j can be calculated as follows:

$$Q_{i,j}=Q_{i+1,j-1}+\Delta R_j\Delta A_i$$

(3)

where Q_i,j is the outflow from travel time zone i at time step j (m³); and Q_i+1,j−1 is the outflow from time zone i + 1 at time step j − 1 (m³).

Unlike runoff routing, sediment transport simulation gave consideration to the sediment particle size and the hydraulic characteristic of the basin. It was calculated using Equation (4), as suggested by Williams [68]:

$$RY={\displaystyle \sum _{i=1}^n{Y}_i{e}^{-B{T}_i\sqrt{D{50}_{\mathrm{i}}}}}$$

(4)

where RY is the sediment yield for the entire basin (t); Y_i is the sediment yield for the sub-catchment i before routing to channel (t); B is the routing coefficient; T_i is the travel time from the sub-catchment i to the basin outlet (h); and D50_i is median particle diameter of the sediment for sub-catchment i (mm).

2.2.2. Extraction of the Terrace Units and Vegetation Units Watershed

The impact of terrace and vegetation units on water and sediment flow depends upon their size and position in the catchment. The flow direction of the natural hillslope is shown in Figure 4a. When there is a terrace in the hillslope, the connectivity of its original flow pathways is broken, and the whole slope could be divided into three segments: upstream section, terrace section, and downstream section (Figure 4b). Here, we define the upstream section as the watershed of terrace (or runoff contributing area). The area impacted by the terrace units in terms of water and sediment reduction is the total of the terrace and its watershed.

The terrace watershed was extracted by searching upward from the terrace cell based on deterministic-8 (D8) flow direction method [69]. For each cell, once there was a route that water could follow to reach the terrace cells, the cell was considered to belong to the terrace watershed. The extraction of vegetation watershed was the same as that for the terrace watershed.

2.2.3. Consideration of Terrace Units in the Time-Area Method

The revised time-area method simplifies a level terrace as a dynamic water tank, and its storage capacity is the product of the terrace area and the embankment height (Figure 5a,b). In this study, we made the assumption that the stream flow in a terrace should only be considered as the overflow, regardless of the drainage discharge. Namely, when the water trapped by the terrace exceeds the terrace storage capacity during heavy rainfall events, the surplus water will overflow (Figure 5c,d).

Referring to Equation (3), when there are terrace units occurring in time zone i, the outflow Q_i,j from time zone i at time step j is revised, as follows:

$$Q_{i,j}=Q_{i+1,j-1}+\Delta R_j\Delta A_i-\Delta T_{\mathrm{i},\mathrm{j}}$$

(5)

where Q_i+1,j−1 is the outflow from time zone i + 1 at time step j − 1 (m³); ΔT_i,j is the increased water ponded in terraces that are located in time zone i at time step j (m³). When ΔT_i,j equals 0, the storage capacity of the terraces is filling up and they can no longer have a retention function.

Before filling up, the current terrace storage volume is the sum of its previous storage volume and the current terrace watershed inflows. Due to the possible distribution of several terraces in a time zone and due to the fact that each terrace may correspond to multiple cells, the statistics of the water interception of each terrace requires a large amount of calculations. In order to simplify the complex problem, we estimated the terrace watershed inflows as a percentage of the runoff generation in the time zone, and the percentage equals the area ratio of all terraces control area to the time zone. Then ΔT_i,j is calculated as follows:

$$\Delta T_{i,j}=T_{i,j}-T_{i,j-1}$$

(6)

$$T_{i,j}=\{\begin{array}{c}T{c}_i\\ T_{i,j-1}+\Delta R_j\Delta A_i\times T{u}_i\end{array}\begin{array}{l}(\mathrm{Filling}\mathrm{up})\\ (\mathrm{Unfilling}\mathrm{up})\end{array}$$

(7)

$$T{u}_i=\Delta A_{Tcontrol,i}/\Delta A_i$$

(8)

$$T{c}_i=\Delta A_{terrace,i}\times d_i$$

(9)

where T_i,j is the amount of stored water in all terraces in time zone i at present time step j (m³); T_i,j−1 is the amount of stored water in all terraces in time zone i at previous time step j − 1 (m³); Tc_i is the maximum water storage of all terraces in the time zone i (m³); Tu_i is the ratio of all terrace control area in the time zone i; ΔA_Tcontrol,i is the area of all terrace control area in time zone i (m²); ΔA_terrace,i is the total area of terraces in time zone i (m²); and d_i is embankment height of terrace (m). As runoff is the main carrier of sediment, the increased sediment stored by terraces in time zone i at time step j equals the sediment yield in time zone i at time step j multiplied by the runoff trapped rate of terraces in time zone i at time step j (namely, the ratio of ΔT_i,j to ΔR_jΔA_i).

Taking Equation (6) and Equation (7) into Equation (5), we can get the outflow Q_i,j as follows:

$$Q_{i,j}=\{\begin{array}{c}{Q}_{i+1,j-1}+\Delta R_j\Delta A_i-T{c}_i+T_{i,j-1}\\ Q_{i+1,j-1}+\Delta R_j\Delta A_i\times (1-T{u}_i)\end{array}\begin{array}{l}(\mathrm{F}\mathrm{i}\mathrm{l}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{u}\mathrm{p})\\ (\mathrm{U}\mathrm{n}\mathrm{f}\mathrm{i}\mathrm{l}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{u}\mathrm{p})\end{array}$$

(10)

2.2.4. Consideration of Vegetation Units in the Time-area Method

Unlike terraces, which can act as a water tank, vegetation does not have a direct storage volume, but it can resist the confluence process of runoff and sediment. Vegetation coverage and canopy density are the main indicators influencing vegetation’s impact on water and soil conservation [70,71,72]. In this study, we chose vegetation cover that was easily derived to incorporate the vegetation module into the time-area method.

Similar to Equation (5), when there are vegetation units occurring in time zone i, the outflow Q_i,j from time zone i at time step j is as follows:

$$Q_{i,j}=Q_{i+1,j-1}+\Delta R_j\Delta A_i-\Delta Ve{g}_{\mathrm{i},\mathrm{j}}$$

(11)

where ΔVeg_i,j is the water or sediment trapped by vegetation units in time zone i at time step j (m³).

ΔVeg_i,j is calculated as follows:

$$\Delta Ve{g}_{\mathrm{i},\mathrm{j}}=(Q_{i+1,j-1}+\Delta R_j\Delta A_i)\times V{u}_i\times \overline{f(Veg)}$$

(12)

$$V{u}_i=\Delta A_{Vcontrol,i}/\Delta A_i$$

(13)

where Vu_i is the ratio of all the vegetation control area in the time zone i; ΔA_Vcontrol,i is the total area of vegetation in time zone i (m²); Veg is vegetation coverage (%); $f(Veg)$ is function of Veg, it refers to runoff or sediment retention rate of the forestland or grassland with certain vegetation coverage, and $\overline{f(Veg)}$ is the average value of runoff or sediment retention rates in the time zone i.

Most flume test researches in the Loess Plateau lack the complete information about the runoff and sediment retention rate of grassland and forestland with different cover [73,74,75,76]. Xiong et al. [77] systematically deconstructed the experimental data from different slope runoff plots in the Loess Plateau, and summarized benefit indices of runoff and sediment reduction by forestland and grassland of different qualities in years with different runoff and sediment levels [77], as shown in

Table 3. The reduction percentage in runoff and sediment generation of forestland and grassland of different quality.

Vegetation Coverage (%)		Dry Year		Normal Year		Wet Year
		Runoff Reduction (%)	Sediment Reduction (%)	Runoff Reduction (%)	Sediment Reduction (%)	Runoff Reduction (%)	Sediment Reduction (%)
Forest	70	100	100	100	98	76.5	57.7
	60	100	100	96.5	92.9	72.2	51
	50	99	99	90.1	86.9	64.2	46.2
	40	94	96	73.2	69.8	48.8	33.3
	30	80	89	52	48.2	28.4	19.2
	20	55	73	26.7	20.2	11.1	6.4
Grass	70	100	100	96.3	94.4	64.8	50
	60	100	100	92.6	89.9	59.3	45.1
	50	98	99	83.7	82.5	51.2	40
	40	86	95	67.8	66.5	37.7	30
	30	72	85	42.7	41.8	22.1	16.9
	20	45	69	19.5	18.6	8.2	5.9

. This paper refers to the study of Xiong et al. [77]. For the convenience of distributed calculation, we fitted the data in Table 3 to get the runoff and sediment reduction functions of forestland and grassland under different conditions, as shown in

Table 4. The reduction function of runoff by forestland and grassland.

Land	Dry Year	Normal Year	Wet Year
Forest	y = 0.3619ln(x) − 0.468	y = 0.6145ln(x) − 1.557	y = 0.5551ln(x) − 1.565
Grass	y = 0.4537ln(x) − 0.854	y = 0.6498ln(x) − 1.748	y = 0.4733ln(x) − 1.357

and

Table 5. The reduction function of sediment by forestland and grassland.

Land	Dry Year	Normal Year	Wet Year
Forest	y = 0.2136ln(x) + 0.133	y = 0.6429ln(x) − 1.701	y = 0.4239ln(x) − 1.222
Grass	y = 0.2527ln(x) − 0.028	y = 0.6384ln(x) − 1.721	y = 0.3669ln(x) − 1.053

, x means vegetation coverage (%) and y means runoff and sediment reduction rates.

2.2.5. Model Performance Evaluation Criteria

Nash–Sutcliffe efficiency (NSE) was used to evaluate the performance of the simulation:

$$\mathrm{Nash}=1-{\displaystyle \sum _{i=1}^N{{(\mathrm{O}}_i-{\mathrm{E}}_i)}^2}/{\displaystyle \sum _{i=1}^N{{(\mathrm{O}}_i-\overline{\mathrm{O}})}^2}$$

(14)

where O_i is observed data (runoff discharge or sediment discharge); E_i is simulated data (runoff discharge or sediment discharge); $\overline{\mathrm{O}}$ is average observed data; N is number of values. NSE varies from negative infinity to 1, where NSE closer to 1 indicates a better simulation.

A total of eight isolated storms with observed runoff and sediment yield were selected to calibrate and verify the model. There were five in the 1980s, and three in the 2010s, and each storm was encoded with its start time as Nos. year/day/hour. Among them, Nos. 1981/203/17, 1983/215/22 and 1983/235/16 were used for calibration, while Nos. 1988/199/13, 1989/203/19, 2006/195/5, 2006/224/8 and 2010/263/20 were used for validation.

The runoff and sediment simulation was implemented in four cases: O1—original simulation without considering terrace and vegetation practice; R1—revised with vegetation module; R2—revised with terrace module; R3—revised with vegetation and terrace modules. For the five events in the 1980s, as the terrace data was unavailable and terraces just accounted for a small proportion of the area, only O1 and R1 were simulated. For three events in the 2010s, all four cases were simulated. The terraces in the study area are all level terraces with good quality [78], and we made an assumption that all terraces had an embankment height of 20 cm. For vegetation, the retention functions of runoff and sediment were chosen according to the rainfall amount of each event and the rainfall condition of the two days before each event. Here, for Nos. 1983/215/22, 1983/235/16, 1988/199/13, 2006/224/8 and 2010/263/20 functions of normal period were chosen in Table 4 and Table 5, while for Nos. 1981/203/17, 1989/203/19 and 2006/195/5 functions of wet period were chosen.

2.3. Data Source

Input data for the model included hourly precipitation, DEM, land use, vegetation coverage, soil data, and particle size of the sediment. In addition, observed runoff and sediment discharge of hydrological stations were used in the model’s calibration and validation.

(1)

Hourly precipitation was collected from 10 rain-gauge stations (Figure 1). The gauge data were interpolated by the inverse distance weighted (IDW) method to acquire the spatial data.

(2)

Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM) was selected to extract topographical information, river network, sub-basins and isochrones (http://www.gscloud.cn/).

(3)

Land use and vegetation coverage were derived from Landsat images (http://www.gscloud.cn/). Land use was classified into 11 types, including cropland (slope < 6°), cropland (6° < slope < 25°), cropland (slope ≥ 25°), forest, shrub, open woodland, immature forest land and orchard, grassland, water area, developed land, and other land. Vegetation coverage was derived based on its relation with normalized difference vegetation index (NDVI) [79], while the vegetation coverage of cropland was not included in routing process. Land use and vegetation coverage data from 1978 and 2010 were used to represent data for 1980s and 2010s, respectively.

(4)

Terrace data in 2012 was acquired from the Yellow River Conservancy Commission (YRCC), which was interpreted from ZY-3 images with a spatial resolution of 2.5 m and an accuracy of 94% [78].

(5)

The soil types, sand content, silt content, clay content, organic carbon, and gravel content were derived from the Harmonized World Soil Database (HWSD), and the soil saturated moisture of each soil type were determined from the software Soil-Plant-Atmosphere-Water Field & Pond Hydrology (SPAW).

(6)

The measured discharge, sediment concentration and median of sediment particle size (D₅₀) at the Pianguanhe hydrological station were acquired from the YRCC.

All data were reproduced at 30 × 30 m spatial resolution and projected to Albers, using the World Geodetic System-84 (WGS84) datum. Data were processed by ENVI5.1 (Harris Corporation, Melbourne, FL, USA), GisNet [80], ArcGIS10.4.1 (Environmental Systems Research Institute, Redlands, CA, USA), or programs written in IDL (Interactive Data Language).

3. Results

3.1. Watershed of Terrace Unit and Vegetation Unit

The terrace ratio is the proportion of terrace area to water and soil loss area [78]. In 2012, the total terrace area of Pianguanhe was 189.32 km², accounting for 9.87% of the whole basin and giving a terrace ratio of 10.18%. The terrace ratio of Pianguanhe was 1.80% in 1979 [50]. In view of data availability, we chose terrace data from 2012 (to represent the 2010s) and extracted its terrace watershed. The terrace watershed was 169.90 km² in 2012. The terrace control area is the sum of terrace area and terrace watershed and accounted for 1.89 times the area of the terrace, as shown in Figure 6a.

The area of forestland and grassland in Pianguanhe was 1328.70 km² in 1980s and 1349.17 km² in the 2010s, giving an increase of 20.47 km². The average vegetation coverage was 22.27% in 1980s and 62.74% in the 2010s, giving an increase of 40.47%. Thus, the vegetation quality demonstrated a significant increase compared to the area. We extracted the vegetation watersheds in the 1980s and 2010s. Figure 6b shows the vegetation watershed of Pianguanhe in the 2010s. Vegetation control area was the sum of vegetation units and vegetation watershed. It was 1675.77 km² in the 1980s and 1693.96 km² in the 2010s, indicating that it did not change very much.

3.2. Validation of Runoff Discharge

Figure 7a presents the event-based comparison between the measured runoff discharge and the runoff discharge predicted under different cases. The simulated runoff discharge in most events were comparable to their measured values. For the five events in the 1980s, the peak values of R1 were lower than those of O1. For the three events in the 2010s, the peak values rank from lowest to highest: R3 < R1 < R2 < O1 (Figure 7a). This shows that the revised simulation with vegetation can significantly reduce the runoff peak compared to the simulation with terrace, and that the simulation considering both terrace and vegetation was closest to the estimated value.

Comparison of the simulated runoff with the measured runoff is shown in Figure 7b–e. Figure 7b shows the validation of O1, in which the five events in the 1980s achieved an average NSE of 0.46, while the three events in the 2010s achieved an average NSE of −15.29. The data points of the 1980s are distributed close to the 1:1 line, while the data points of the 2010s are distributed farther away from the 1:1 line. It shows that the original model did not achieve good performance for recent years, since the surface had undergone significant change. Figure 7c shows that the average NSE of the 1980s and 2010s both increased in R1, with the average NSE of the 2010s showing a greater increase (value of −1.32), and the data points of the 2010s distributed closer to the 1:1 line than in the O1 validation shown in Figure 7b. Figure 7d shows the R2 validation, and indicates that the average NSE of the 2010s has also increased compared with that of O1, but that the extent of the increase is less than for R1. Figure 7e shows the R3 validation and indicates that the average NSE of the 2010s is 0.39, which for the first time is positive, and that the data points of the 2010s are closely distributed around the 1:1 line. Overall, the simulation accuracy of the 2010s is highest in R3, followed by R1. These figures show that the revised model was better at simulating the runoff in recent years and can reflect the effect of significant surface change (i.e., slope terracing and revegetation) on runoff.

3.3. Validation of Sediment Discharge

The hourly measured sediment concentration (kg/m³) along with the measured runoff discharge were converted to get sediment discharge (t/h). Figure 8a shows the event-based comparison between the measured and predicted sediment discharge under different cases. The simulations of sediment discharge were in good agreement with their measured values in most events. The difference between the sediment discharge of O1 and measured values is smaller than the runoff discharge between O1 and measured values in recent years; it was because the sediment reduction was partly achieved in the sediment yield process through the C and P factors in MUSLE. For the five events in the 1980s, the peak values of R1 were lower than those of O1. For the three events in the 2010s, the peak values rank from lowest to highest were: R3 < R1 < R2 < O1 (Figure 8a). This shows that the revised simulation with vegetation can significantly reduce the sediment peak compared to the simulation with terrace, and that the simulation considering both terrace and vegetation was closest to the estimated value.

NSE was also used to evaluate the performance of the simulated sediment discharge, as shown in Figure 8b–e. Figure 8b shows the O1 validation, in which the five events in the 1980s achieved an average NSE of 0.08, while the three events in the 2010s achieved an average NSE of −2.47. It shows that the original model did not achieve good performance in the recent years. Figure 8c shows the R1 validation, in which the average NSE of the 1980s and 2010s both increased. The data points of the 1980s and 2010s were distributed closer to the 1:1 line than in Figure 8b. Figure 8d shows the R2 validation, in which the average NSE of the 2010s have also increased, and greater than that of R1, which indicates that terrace can significantly reduce sediment. Figure 8e shows the R3 validation, in which the average NSE of the 2010s is 0.37, and the data points of the 2010s are distributed closely around the 1:1 line. Overall, the accuracy of the simulation in the 2010s is highest in R3, followed by R2. These figures show that the revised model performs better in recent years, and can reflect the effect of significant surface change (i.e., slope terracing and revegetation) on sediment yield.

4. Discussion

4.1. Effect of Terraces and Vegetation on Runoff Reduction

To quantify the reduction effects of vegetation and terracing during runoff and sediment routing process, the absolute reduction in unit area (AR) and the reduction rate (RR) were calculated. In the case of runoff, for example, we first summed the hourly runoff discharge to obtain the total runoff of each rainfall, and counted the total simulated runoff of each event per case. We then calculated the AR of R1 and R2, and the RR of R1, R2, and R3, as shown in

Table 6. The efficiency of runoff reduction by vegetation and terraces in Pianguanhe Basin.

No. of Storm (year/day/hour)	Class	Practice
		Vegetation		Terrace		Vegetation and Terrace
		RR1 (%)	AR1 (m3/km2)	RR2 (%)	AR2 (m3/km2)	RR3 (%)	RR1 + RR2 (%)	Difference (%)
1981/203/17	Calibration	35.95	4096.15	-	-	-	-	-
1983/215/22		30.42	1642.02	-	-	-	-	-
1983/235/16		39.60	642.21	-	-	-	-	-
1988/199/13	Validation	35.77	3959.85	-	-	-	-	-
1989/203/19		26.18	494.47	-	-	-	-	-
2006/195/5		48.22	3318.46	26.13	9122.72	69.02	74.35	−5.33
2006/224/8		53.17	5034.97	26.26	18,060.69	73.26	79.43	−6.17
2010/263/20		51.66	5223.21	27.57	20,244.11	73.30	79.23	−5.93

Annotation: RR₁ and AR₁ mean the RR and AR of vegetation, respectively; RR₂ and AR₂ mean the RR and AR of terrace, respectively; RR₃ means the RR of both vegetation and terrace; Difference denotes RR₃ minus the sum of RR₁ and RR₂.

. The AR and RR were calculated as follows:

$$A{R}_i=\frac{y-y_i}{A_i},(i=1,2)$$

(15)

$$R{R}_i=\frac{y-y_i}{y}\times 100\%,(i=1,2,3)$$

(16)

where y is total simulated runoff of O1 (m³); and y_i is total simulated runoff of R1, R2 or R3 (m³), A_i is total area of R1 or R2 (km²).

As shown in Table 6, the average RR of R1 on runoff was 33.58% in the 1980s. For the three events in the 2010s, the average RR of R1 was 51.02%, the average RR of R2 was 26.65%, and the average RR of R3 was 71.86%. The average AR of R1 was 3051.42 m³/km² per event, and it was 15809.18 m³/km² for R2. The results show that the runoff reduction rate of vegetation was significantly higher than that of terraces, as the area of vegetation is seven times larger than that of terraces. However, terraces could reduce more runoff per unit area. Influenced by the revegetation and increase in vegetation coverage, the average RR of R1 in recent years increased by 17.44% over that in the 1980s. In general, the effectiveness of runoff reduction was highest under R3. Besides, the sum of RR₁ and RR₂ is inconsistent with RR₃. The reason for the latter might be that when the terraces had played the role of water reduction, the water reducing capability of vegetation in the same time zone would be decreased. The assumption that all terraces in the Pianguanhe Basin are of good quality and have the same embankment height may be an oversimplification that overestimated the runoff reduction efficiency of the terraces.

4.2. Effect of Terraces and Vegetation on Sediment Reduction

We counted the total simulated sediment yield of each event per case, and also calculated the AR of vegetation and terraces, and the RR of R1, R2, and R3, as shown in

Table 7. The efficiency of sediment reduction by vegetation and terraces in the Pianguanhe Basin.

No. of Storm (year/day/hour)	Class	Practice
		Vegetation		Terrace		Vegetation and Terrace
		RR1 (%)	AR1 (t/km2)	RR2 (%)	AR2 (t/km2)	RR3 (%)	RR1 + RR2 (%)	Difference (%)
1981/203/17	Calibration	13.66	796.05	-	-	-	-	-
1983/215/22		11.52	299.66	-	-	-	-	-
1983/235/16		17.31	96.91	-	-	-	-	-
1988/199/13		15.37	748.46	-	-	-	-	-
1989/203/19		8.71	129.34	-	-	-	-	-
2006/224/8	Validation	33.22	341.94	23.62	1765.63	54.02	56.84	−2.83
2006/195/5		30.36	528.25	24.08	3041.80	51.75	54.44	−2.69
2010/263/20		33.07	659.02	25.87	3743.10	55.80	58.94	−3.14

Annotation: RR₁ and AR₁ mean the RR and AR of vegetation, respectively; RR₂ and AR₂ mean the RR and AR of terrace, respectively; RR₃ means the RR of both vegetation and terrace; Difference denotes RR₃ minus the sum of RR₁ and RR₂.

. The average RR of R1 for sediment was 13.31% in the 1980s. For the three events in the 2010s, the average RR of R1 was 32.22%, the average RR of R2 was 24.52%, and the average RR of R3 was 53.85%. The average RR of R1 in recent years increased by 18.91% over that in the 1980s. The average AR of R1 was 449.95 t/km², and 2850.17 t/km² for R2. The results also show that the sediment reduction rate of vegetation was higher than that of terraces, but that terraces could reduce more sediment per unit area. In general, the effectiveness of sediment reduction of R3 was highest. Besides, the sum of RR₁ and RR₂ was inconsistent with RR₃. As mentioned previously, the reason might be that when terraces had played the role of sediment reduction, the sediment reducing efficiency of vegetation in the same time zone would be decreased. The assumption that all terraces in the Pianguanhe Basin are of good quality and have the same embankment height may have overestimated the sediment reduction efficiency of terraces.

4.3. Reliability Analysis of the Water and Sediment Reduction Efficiency of Terraces and Vegetation

Liu et al. [78,81] introduced a runoff coefficient and sediment yield coefficient to discuss the flood and sediment yield for different vegetation conditions at the catchment scale in the Loess Plateau. To analyze the reliability of the simulation results, we also calculated the runoff coefficient and sediment yield coefficient of each event in the Pianguanhe Basin, and compared with the results of Liu et al. [78], as shown in Figure 9a,b. The runoff coefficient is the runoff yield per unit of precipitation per unit area, while the sediment yield coefficient is the sediment yield per unit of precipitation per unit area. In Liu et al. [78], the runoff coefficient was calculated based on annual runoff and precipitation data, and the latter only considered rainfall greater than 25 mm. In this study, the runoff coefficient was calculated based solely on storm event discharge and rainfall data. Thus, there is a disparity between the calculation results of these two studies, but yet the results are comparable to those of Liu et al. [78]. The runoff coefficient decreased along with the increase of the vegetation coverage. The vegetation coverage of the 1980s and 2010s was 16.17% and 46.29%, respectively, and the average runoff coefficient decreased from 0.12 m³/(mm·km²) in the 1980s to 0.07 m³/(mm·km²) in recent years.

Similarly, the sediment yield coefficients calculated from the event data in this study are different from those of Liu’s results based on calculations using annual data, and yet the results are comparable. The sediment yield coefficients also decreased along with the increase of vegetation coverage. The average sediment yield coefficient decreased from 80.68 t/(mm·km²) in the 1980s to 21.93 t/(mm·km²) in recent years.

4.4. Limitations and Potential Improvements

While this study assessed the impact of terraces and vegetation practice on runoff and sediment routing process using the revised time-area method, it still has several assumptions and limitations in the method that need to be clarified and studied further. First, embankment height is a key indicator to evaluate the quality and storage volume of terraces. The assumption that all terraces in the Pianguanhe Basin are of good quality and set the embankment height with a fixed value may overestimate the runoff and sediment reduction efficiency of the terraces. This paper referred to the study of Xiong et al. [77] as the limitation of the flume test data. While the runoff and sediment retention function of vegetation is not independent of slope steepness or vegetation structure, these factors have not been considered yet. In addition, the effect of engineering measures such as check dam and human factors such as mining and road construction should also be considered to further improve simulation efficiency.

5. Conclusions

In this study, we added the impact of terrace and vegetation practice on runoff and sediment routing in the time-area method. The revised time-area method was integrated into the LCM-MUSLE model which is suitable to estimate the water and sediment yield in the Loess Plateau. Eight isolated storm events in the 1980s and 2010s in Pianguanhe Basin were selected to calibrate and verify the original LCM-MUSLE model and its revised version. It is shown that the original model was not applicable in the more recent years, since the surface had changed significantly as a result of revegetation and terrace engineering. The revised model considered the impact of vegetation and terracing on runoff and sediment routing and its accuracy had been improved significantly.

The effect of the level terraces and vegetation was parameterized effectively according to their location, size, embankment height, and vegetation coverage. These parameters could be easily obtained and were used to represent the landscape heterogeneity at the catchment scale. Besides, the revised time-area method was loosely coupled with the LCM-MUSLE model. Therefore, the method could be readily applied in other regions and integrated into other hydrological models and erosion models. Consequently, this study provides a generalized method to quantitatively assess the impact of terrace and vegetation practice on runoff and sediment reduction at the catchment scale, which has great significance in runoff change analysis and implementation of soil and water conservation.