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

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.


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.
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.

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 2 , 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 3 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.
Water 2019, x, x FOR PEER REVIEW 3 of 22 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.

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 2 , 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 3 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.

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 Tables 1 and 2. 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

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 Tables 1 and 2. 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 Water 2019, 11, 803 5 of 20 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.concentration for the grid; Sd: percent sand content (%); Si: percent silt content (%); Ci: percent clay content (%); C: percent organic carbon content of the layer (%); Lslp: slope length (m); θ: gradient of the slope (°); cv: vegetation coverage (%); rock: percent rock in the top soil layer (%).
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.Framework of revised integrated model.Ti,j is the amount of stored water or sediment in all terraces in travel time zone i at present time step j (m 3 or t); Ti,j-1 is the amount of stored water or sediment in all terraces in i at previous time j-1 (m 3 or t); and Tci is water storage capacity of all terraces in i (m 3 ).

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 ΔA1, ΔA2, •••, ΔAn, and their corresponding travel time is τ1, τ2, •••, Framework of revised integrated model.T i,j is the amount of stored water or sediment in all terraces in travel time zone i at present time step j (m 3 or t); T i,j−1 is the amount of stored water or sediment in all terraces in i at previous time j − 1 (m 3 or t); and Tc i is water storage capacity of all terraces in i (m 3 ).

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 1 , ∆A 2 , •••, ∆A n , and their corresponding travel time is τ 1 , τ 2 , •••, τ n , respectively.Then, we set ∆τ 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 If j < m, k = j; if j ≥ m, k = m, and if j > n, ∆A = 0.
For n < m: 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: where Q i,j is the outflow from travel time zone i at time step j (m 3 ); and Q i+1,j−1 is the outflow from time zone i + 1 at time step j − 1 (m 3 ).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]: 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).
For n ≥ m -r 1 1 For n < m: For travel time zone i, its runoff discharge at time step j can be calculated as follows: where Qi,j is the outflow from travel time zone i at time step j (m 3 ); and Qi+1,j−1 is the outflow from time zone i+1 at time step j−1 (m 3 ).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]: where RY is the sediment yield for the entire basin (t); Yi is the sediment yield for the sub-catchment i before routing to channel (t); B is the routing coefficient; Ti is the travel time from the sub-catchment i to the basin outlet (h); and D50i is median particle diameter of the sediment for sub-catchment i (mm).

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

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.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.

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).

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).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.

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).A terrace is the area of terrace; d is embankment height of terrace; P is rainfall; f is the infiltration; Q l is the interflow; Q j−1 is the inflow; Q d is the overflow; T j is the water stored in the terrace at the present time step j; T j−1 is the water stored in terrace at previous time j − 1; and Q j is the outflow.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: where Q i+1,j−1 is the outflow from time zone i + 1 at time step j − 1 (m 3 ); ∆T i,j is the increased water ponded in terraces that are located in time zone i at time step j (m 3 ).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: where T i,j is the amount of stored water in all terraces in time zone i at present time step j (m 3 ); T i,j−1 is the amount of stored water in all terraces in time zone i at previous time step j − 1 (m 3 ); Tc i is the maximum water storage of all terraces in the time zone i (m 3 ); 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 2 ); ∆A terrace,i is the total area of terraces in time zone i (m 2 ); 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:

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: where ∆Veg i,j is the water or sediment trapped by vegetation units in time zone i at time step j (m 3 ).∆Veg i,j is calculated as follows: Water 2019, 11, 803 9 of 20 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 2 ); 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 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.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 Tables 4 and 5, x means vegetation coverage (%) and y means runoff and sediment reduction rates.

Model Performance Evaluation Criteria
Nash-Sutcliffe efficiency (NSE) was used to evaluate the performance of the simulation: 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

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.
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 50 ) at the Pianguanhe hydrological station were acquired from the YRCC.

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 2 , 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 2 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 Pianguanhe hydrological station were acquired from the YRCC.

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 2 , 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 2 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 2 in 1980s and 1349.17km 2 in the 2010s, giving an increase of 20.47 km 2 .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.77km 2 in the 1980s and 1693.96km 2 in the 2010s, indicating that it did not change very much.The area of forestland and grassland in Pianguanhe was 1328.70 km 2 in 1980s and 1349.17km 2 in the 2010s, giving an increase of 20.47 km 2 .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.77km 2 in the 1980s and 1693.96km 2 in the 2010s, indicating that it did not change very much.

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.
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

Validation of Sediment Discharge
The hourly measured sediment concentration (kg/m 3 ) 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 Figures 8b-8e.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 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.

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 AR and RR were calculated as follows: where y is total simulated runoff of O1 (m 3 ); and y i is total simulated runoff of R1, R2 or R3 (m 3 ), A i is total area of R1 or R2 (km 2 ).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 3 /km 2 per event, and it was 15809.18m 3 /km 2 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 Water 2019, 11, 803 15 of 20 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 1 and RR 2 is inconsistent with RR 3 .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.

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 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 2 , and 2850.17t/km 2 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 1 and RR 2 was inconsistent with RR 3 .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.

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 3 /(mm•km 2 ) in the 1980s to 0.07 m 3 /(mm•km 2 ) 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 2 ) in the 1980s to 21.93 t/(mm•km 2 ) in recent years.

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.

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.

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.

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.

Figure 1 .
Figure 1.Location of the study area.

Figure 1 .
Figure 1.Location of the study area.

Figure 2 .
Figure 2.Framework of revised integrated model.Ti,j is the amount of stored water or sediment in all terraces in travel time zone i at present time step j (m 3 or t); Ti,j-1 is the amount of stored water or sediment in all terraces in i at previous time j-1 (m 3 or t); and Tci is water storage capacity of all terraces in i (m 3 ).

Figure 2 .
Figure 2.Framework of revised integrated model.T i,j is the amount of stored water or sediment in all terraces in travel time zone i at present time step j (m 3 or t); T i,j−1 is the amount of stored water or sediment in all terraces in i at previous time j − 1 (m 3 or t); and Tc i is water storage capacity of all terraces in i (m 3 ).

Figure 4 .
Figure 4. Schematic diagram of flow direction of (a) natural hillslope, and (b) the hillslope with vegetation and terrace practice.Numbers in the grid are elevations.Orange represents natural hillslope unit.Green represents vegetation unit.Yellow represents terrace unit.Arrows represent flow direction.

Figure 4 .
Figure 4. Schematic diagram of flow direction of (a) natural hillslope, and (b) the hillslope with vegetation and terrace practice.Numbers in the grid are elevations.Orange represents natural hillslope unit.Green represents vegetation unit.Yellow represents terrace unit.Arrows represent flow direction.

Figure 4 .
Figure 4. Schematic diagram of flow direction of (a) natural hillslope, and (b) the hillslope with vegetation and terrace practice.Numbers in the grid are elevations.Orange represents natural hillslope unit.Green represents vegetation unit.Yellow represents terrace unit.Arrows represent flow direction.

Figure 5 .
Figure 5. Schematic diagrams of the hydrological processes and flow distribution of the level terrace unit.(a) Section profile of a level terrace.(b) The generalization of the level terrace in (a).The flow distribution assumes that the soil profile of the terrace is deep enough for subsurface flow generation.The influence of the flow in different scenarios is listed from (c,d).A terrace is the area of terrace; d is embankment height of terrace; P is rainfall; f is the infiltration; Q l is the interflow; Q j−1 is the inflow; Q d is the overflow; T j is the water stored in the terrace at the present time step j; T j−1 is the water stored in terrace at previous time j − 1; and Q j is the outflow.

Figure 6 .
Figure 6.Extraction of terrace watershed (a) and vegetation watershed (b) of Pianguanhe in 2010s.

Figure 6 .
Figure 6.Extraction of terrace watershed (a) and vegetation watershed (b) of Pianguanhe in 2010s.

Figure 7 .
Figure 7.Comparison of observed and simulated runoff for eight events in the Pianguanhe Basin.(a) Event rainfall-runoff simulation; (b) Runoff in O1 simulation; (c) Runoff in R1 simulation; (d) Runoff in R2 simulation; (e) Runoff in R3 simulation.

Figure 7 .
Figure 7.Comparison of observed and simulated runoff for eight events in the Pianguanhe Basin.(a) Event rainfall-runoff simulation; (b) Runoff in O1 simulation; (c) Runoff in R1 simulation; (d) Runoff in R2 simulation; (e) Runoff in R3 simulation.

21 Figure 8 .
Figure 8.Comparison of observed and simulated sediment discharge for eight events in the Pianguanhe Basin.(a) Event rainfall sediment simulation; (b) Sediment in O1 simulation; (c) Sediment in R1 simulation; (d) Sediment in R2 simulation; (e) Sediment in R3 simulation.

Figure 8 .
Figure 8.Comparison of observed and simulated sediment discharge for eight events in the Pianguanhe Basin.(a) Event rainfall sediment simulation; (b) Sediment in O1 simulation; (c) Sediment in R1 simulation; (d) Sediment in R2 simulation; (e) Sediment in R3 simulation.
Water 2019, x, x FOR PEER REVIEW 17 of 21 coverage.The average sediment yield coefficient decreased from 80.68 t/(mm•km 2 ) in the 1980s to 21.93 t/(mm•km 2 ) in recent years.

Figure 9 .
Figure 9.Comparison of runoff generation coefficient (a) and sediment yield coefficient (b) in this study and in Liu's study.

Author
Contributions: J.B. conceived and designed the experiments; J.B. and Y.Z performed the experiments; Y.G. contributed analysis tools; J.B., S.Y. and X.L. analyzed the data; J.B. and Y.Z.prepared figures; J.B. wrote

Figure 9 .
Figure 9.Comparison of runoff generation coefficient (a) and sediment yield coefficient (b) in this study and in Liu's study.

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

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

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

Table 4 .
The reduction function of runoff by forestland and grassland.
O i is observed data (runoff discharge or sediment discharge); E i is simulated data (runoff discharge or sediment discharge); 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. where /215/22, 1983/235/16, 1988/199/13, 2006/224/8 and 2010/263/20 functions of normal period were chosen in Tables 4 and 5, while for Nos.1981/203/17, 1989/203/19 and 2006/195/5 functions of wet period were chosen.

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 RR 1 (%) AR 1 (m 3 /km 2 ) RR 2 (%) AR 2 (m 3 /km 2 ) RR 3 (%) RR 1 + RR 2 (%)
RR 1 and AR 1 mean the RR and AR of vegetation, respectively; RR 2 and AR 2 mean the RR and AR of terrace, respectively; RR 3 means the RR of both vegetation and terrace; Difference denotes RR 3 minus the sum of RR 1 and RR 2 . Annotation:

Table 7 .
The efficiency of sediment reduction by vegetation and terraces in the Pianguanhe Basin.RR 1 and AR 1 mean the RR and AR of vegetation, respectively; RR 2 and AR 2 mean the RR and AR of terrace, respectively; RR 3 means the RR of both vegetation and terrace; Difference denotes RR 3 minus the sum of RR 1 and RR 2 . Annotation: