Physics-Based Simulation of Hydrologic Response and Sediment Transport in a Hilly-Gully Catchment with a Check Dam System on the Loess Plateau, China

: Check dams are among of the most widespread and e ﬀ ective engineering structures for conserving water and soil in the Loess Plateau since the 1950s, and have signiﬁcantly modiﬁed the local hydrologic responses and landforms. A representative small catchment was chosen as an example to study the inﬂuences of check dams. A physics-based distributed model, the Integrated Hydrology Model (InHM), was employed to simulate the impacts of check dam systems considering four scenarios (pre-dam, single-dam, early dam-system, current dam-system). The results showed that check dams signiﬁcantly alter the water redistribution in the catchment and inﬂuence the groundwater table in di ﬀ erent periods. It was also shown that gully erosion can be alleviated indirectly due to the formation of the expanding sedimentary areas. The simulated residual deposition heights ( ∆ h) matched reasonably well with the observed values, demonstrating that physics-based simulation can help to better understand the hydrologic impacts as well as predicting changes in sediment transport caused by check dams in the Loess Plateau.


Introduction
The Chinese Loess Plateau has suffered from severe water and soil loss for decades [1,2]. Many measures, including artificial forestation, terraced farming, and check dam construction, have been implemented to conserve soil and water since the 1950s. Since the first check dam appeared in the Loess Plateau 400 years ago, the effective engineering structure has prevailed all over the Loess Plateau, especially in the Loess Mesa Ravine Region and the Loess Hill Ravine Region, to create productive farmlands and conserve soil and water [3]. There were 122,028 check dams in the Loess Plateau at the end of 2005, which held 2.1 × 10 10 m 3 of sediments and formed 3340 km 2 of dam farmlands [4,5]. The number of check dams and the area of dam-farmlands is expected to double by 2020, with the completion of check dam systems for all the main tributaries of the Yellow River in the Loess Plateau.
Check dams have been shown to be an effective engineering structure to reduce water discharge [6,7] and sediment yields [8][9][10] at the basin scale. For example, Xu et al. [6] applied Soil and Water Assessment Tool (SWAT) in a 7725 km 2 watershed with 6572 check dams and found that the annual runoff was reduced by 14.3%, comparing to when there were very few dams.
Ran et al. [9] compared the sediment retention by check dams in five typical watersheds in the Hekou-Longmen section of the Yellow River and found that the average sediment reduction ratio can reach 60% when the percentage of the basin area above check dams in the catchments reached 3.0%. Previous numerical simulations studying the impacts of check dams mainly focused on water Water 2019, 11, 1161 3 of 22 check dams. These check dams have been in operation for more than 50 years under various current conditions (i.e., some have been fully deposited, some have been partially destroyed, the others still remain in good condition). The 42 check dams effectively prevented sediments being transported downstream into Wudinghe River and further into the Yellow River by surface runoff. Other water and soil conservation measures aiming at slope erosion control such as terraced farmlands and woods planting were started in 1962. Figure 2 compares the land use conditions of 1962 and 2010. The primary land use type 50 years ago was slope farmlands (Figure 2a), which was a major erosion source of the catchment. After 50 years, 60% of the slope farmlands were turned into terraced farmlands (Figure 2b), which have a higher water erosion resisting performance. Large areas of grasslands were transformed into sparse woods (Robinta pseudoacacia and Platycladus orientalis) or orchards (apple tree and peach tree) by local farmers. Table 1 summarizes the data set used in this study. Detailed field measurements were conducted from April to September 2014. The soil sample tests provide supplemental soil information including saturated hydraulic conductivities (K sat ) and soil water characteristic curves of various land use types at different sampling sites (Figure 1b). Two wells were monitored for water table data set, which was used for calibrating the initial groundwater table. The residual deposition heights (∆h) behind each dam, defined as the height between dam crest and the surface elevation behind the dam, were also carefully measured to present the sedimentary processes. The field data provides a reference for how the catchment functions currently, and serves as validation for the simulated responses. Ten rainfall-runoff events recorded by the local Suide Soil and Water Conservation Station from 1962 to 1964 were used for model calibration. Two topography maps of WMG catchment were used to construct the computational meshes. Engineering measures focusing on gully erosion control (i.e., check dams) were carried out in the catchment to alleviate soil and water loss since 1953, before which the average erosion rate was 18000 t km −2 per year. 42 check dams were constructed from 1953 to 1959, but most of the dams were destroyed by heavy storms due to low design criteria. Figure 1b shows the 22 currently existing check dams. These check dams have been in operation for more than 50 years under various current conditions (i.e., some have been fully deposited, some have been partially destroyed, the others still remain in good condition). The 42 check dams effectively prevented sediments being transported downstream into Wudinghe River and further into the Yellow River by surface runoff. Other water and soil conservation measures aiming at slope erosion control such as terraced farmlands and woods planting were started in 1962. Figure 2 compares the land use conditions of 1962 and 2010. The primary land use type 50 years ago was slope farmlands (Figure 2a), which was a major erosion source of the catchment. After 50 years, 60% of the slope farmlands were turned into terraced farmlands (Figure 2b), which have a higher water erosion resisting performance. Large areas of grasslands were transformed into sparse woods (Robinta pseudoacacia and Platycladus orientalis) or orchards (apple tree and peach tree) by local farmers. Table 1 summarizes the data set used in this study. Detailed field measurements were conducted from April to September 2014. The soil sample tests provide supplemental soil information including saturated hydraulic conductivities (Ksat) and soil water characteristic curves of various land use types at different sampling sites (Figure 1b). Two wells were monitored for water table data set, which was used for calibrating the initial groundwater table. The residual deposition heights (Δh) behind each dam, defined as the height between dam crest and the surface elevation behind the dam, were also carefully measured to present the sedimentary processes. The field data provides a reference for how the catchment functions currently, and serves as validation for the simulated responses. Ten rainfall-runoff events recorded by the local Suide Soil and Water Conservation Station from 1962 to 1964 were used for model calibration. Two topography maps of WMG catchment were used to construct the computational meshes. Guangdigou catchment (GDG) is a headwater catchment locating in the south of WMG catchment, with a relatively constant dam-system and less land use changes from the 1950s to now. There are five check dams, and Table 2 shows the characteristics of them. Guangdigou catchment (GDG) is a headwater catchment locating in the south of WMG catchment, with a relatively constant dam-system and less land use changes from the 1950s to now. There are five check dams, and Table 2 shows the characteristics of them.
As shown in Figure 1b, GDG1 is located near the outlet of the main gully and installed in series with GDG2 and GDG4 along the main gully. GDG3 and BTG are located on the outlet of two tributary gullies. The five dams are still in good condition. Compared to other sub-catchments in WMG, the As shown in Figure 1b, GDG1 is located near the outlet of the main gully and installed in series with GDG2 and GDG4 along the main gully. GDG3 and BTG are located on the outlet of two tributary gullies. The five dams are still in good condition. Compared to other sub-catchments in WMG, the landuse change of GDG was relatively simple during the 50 years, with a large area of slope farmlands being replaced by terraced farmlands.

The Integrated Hydrology Model (InHM)
Based on the blueprint of the distributed physics-based hydrological model proposed by Freeze and Harlan [17], InHM was developed to quantitatively simulate, via a fully coupled approach, 3-dimensional variably saturated flow in soil and 2-D flow and sediment transport across land surface [18][19][20][21]. The 3-D Richards' equation was implemented to describe variably saturated flow in soil, while the 2-D diffusion-wave equation coupled with depth-integrated multiple-species sediment transport was applied to describe surface flow movement and sediment transportation. Those surface and subsurface governing equations were discretized in space using the control volume finite element method and coupled in one coherent framework using physics-based first-order flux relationship driven by pressure head gradients. Newton iteration was used to implicitly solve each coupled system of nonlinear equations. More details of InHM can be found in Appendix A.
With the most important and innovative characteristic (i.e., no priori assumption of specific runoff-generation mechanisms), InHM has been successfully employed in many different catchments across the world for event-based or continuous hydrological-response simulations [22,23] as well as hydrologically-driven sediment transport simulations [24,25]. In the Loess Plateau catchments of Wangmaogou and its sub-catchment Guangdigou, InHM is capable of simulating rainfall-runoff processes dominated by the infiltration-excess surface flow mechanism and also rainsplash and hydraulic erosion processes in flood events. Spatially distributed information (e.g., surface flow velocity and sediment flux) of these processes can be provided by the model.

Scenario Setting and Modelling Procedure
Two simulation stages were designed in this study. The first stage is calibration and validation simulations, to obtain the actual values of important and sensitive parameters. In the second stage, which was the focus of the study, annual continuous simulations were conducted to evaluate the hydrologic effects and sedimentary processes of a check dam system. The effects of check dam operation in the early period (1955)(1956)(1957)(1958)(1959)(1960)(1961)(1962) and in the current period (2010-2013) were both studied. In the second stage, the model took us to "revisit" what had happened in the first few years after the construction of check dams and compared the water table changes and sedimentary processes of the current dam-system with the early dam-system.

Calibration and Validation Simulations in WMG Catchment
Calibration and validation simulations were conducted in WMG catchment, using the data of ten rainfall-runoff events recorded at the discharge station ( Figure 1b) from 1961 to 1965, which were the only available observation data (Table 1). Four events were used for calibration and the other six events were used for validation. For a specified flood event, check dams were manually added/removed from the mesh according to their status (already existing or destroyed) during the event.
Simulated water discharges and sediment discharges were compared with observed values and the preliminary results were improved by adjusting surface and subsurface parameters of InHM. Parameters, related to infiltration and runoff-generation (i.e., saturated hydraulic conductivity and Manning's roughness coefficient) and erosion (i.e., rainsplash coefficient and surface erodibility), were carefully adjusted to improve the simulated results in the course of calibration. The Nash and Sutcliffe modelling efficiency (EF) [26], given by Equation (1), was employed here to evaluate the model performance.
where O i was the observed value (water discharge or sediment discharge), O was the average value of the observed values, and S i was the simulated value and n was the number of samples. Considering the lack of available measurements to provide land surface information of WMG catchment from 1955 to 1962, these calibration and validation efforts helped to obtain important parameters (e.g., calibrated K sat for various soil types) to construct a more realistic boundary-value problem for the following simulations.

Annual Continuous Simulations in GDG Catchment
Using the calibrated parameters but different meteorological and land use data, annual continuous simulations were carried out for research period 1953-1962 and 2010-2013. It should be noticed that we had to use event-based calibrated parameters to conduct annual continuous simulations because of the lack of long-term observation data. Considering the fact that most of the observed runoff events in a year only occurred in several rainfall storms in the catchment and InHM performed reasonably well in the ten rainfall-runoff events with various rainfall characteristics (rainfall intensity, rainfall amount and duration, see details in Section 3.1), it was technically sound to conduct the following annual continuous simulations using event-based calibrated parameters.
A WMG-catchment boundary-value problem (BVP) was first built to conduct the calibration simulations in the first stage because the discharge station was located at the outlet of WMG catchment. Another BVP was needed to simulate the influence of check dams because: (1) it was difficult to evaluate the hydrologic effects and sedimentary processes of a complicated check dam system which contained 42 check dams 50 years ago and only 22 now; (2) different land use types of WMG catchment after the construction of check dams also increased the degree of complexity. To simplify the problem, the GDG sub-catchment, within WMG with a relatively constant dam-system and less land use changes, was chosen as an example to study the effects of check dams. As mentioned above, the geological structure of WMG was relatively simple. The soil types and features of GDG gully were the same as those of the rest part of WMG catchment. The parameters calibrated in the WMG BVP were mainly related to soil characteristics (e.g., saturated hydraulic conductivity, rainsplash coefficient) and could be directly used in GDG BVP. It was technically reasonable to conduct the GDG BVP simulations using parameters from WMG BVP simulations.
Four scenarios were designed in this study to evaluate the impacts of check dams on hydrologic response and the sedimentary processes (Table 3). Pre-dam scenario (PD) represented the situation before the construction of the first check dam (i.e., base case), with a two-year simulation because the available precipitation data started in 1953. Single-dam scenario (SD) represented the situation after the construction of GDG1 dam in the downstream. Early dam-system scenario (EDS) represented the conditions after four extra dams (i.e., GDG2-4, BTG) were constructed in the upstream of the gully. Current dam-system scenario (CDS) represented the current conditions including large areas of terraced farmlands and a more than 50-year-old check dam system. Table 3. Description of the four scenarios in the study.

Scenario Name Simulation Time Range Check Dam Involved
The impacts of check dams on water redistribution were compared among the four scenarios, based on water balance calculation. Changes in groundwater table of channel reach A0-A2 ( Figure 3) were analyzed. A series of observation nodes, located on/near the A0-A2 profile at a 5-m (near check dam) and 20 m interval, were set. Pressure head and soil water content were calculated in InHM for each observation nodes at every time step. The surface and subsurface nodes with zero pressure head values together outlined the groundwater table profile of A0-A2 reach. Soil erosion of the GDG catchment in the four scenarios was also evaluated by calculating the eroded sediment mass of different surface zones (Table 4), which was calculated in InHM by the integration of sediment flux through the boundary we set for each surface zone.
Water 2019, 11, x; doi: www.mdpi.com/journal/water (near check dam) and 20 m interval, were set. Pressure head and soil water content were calculated in InHM for each observation nodes at every time step. The surface and subsurface nodes with zero pressure head values together outlined the groundwater table profile of A0-A2 reach. Soil erosion of the GDG catchment in the four scenarios was also evaluated by calculating the eroded sediment mass of different surface zones (Table 4), which was calculated in InHM by the integration of sediment flux through the boundary we set for each surface zone.      [20], calibrated from the first stage.; c Immobile water depth [20], identical for all surface zones.; d Surface residual saturation [20], identical for all surface zones.; e Average height of non-discretized micro-topography [20], identical for all surface zones.; f Rainsplash coefficient [18], calibrated from the first stage.; g Rainsplash depth dampening factor [27], identical for all surface zones.; h Rain intensity exponent [27], identical for all surface zones.; i Rain-induced turbulence coefficient [27], identical for all surface zones.; j Surface erodibility coefficient [27], calibrated from the first stage. The values for CDS scenario were derived from Gao et al. [15].; k Sedimentary fields are usually used as productive farmlands after averagely 2-year deposition, increasing manning's roughness coefficient to 0.12.

Boundary Conditions and Initial Conditions
The 3D meshes for the four scenarios were all constructed by adding layers to 2D surface meshes. The 2D surface meshes of GDG gully for the first three scenarios were constructed based upon the topographic map surveyed in 1953. A surface mesh with 2165 nodes and 4252 triangular elements was first generated for PD scenario. GDG1 dam was then added to the first surface mesh according to characteristics of GDG1 in Table 2 to generate the second surface mesh for the SD scenario, which has 2243 nodes and 4408 triangular elements. The other four dams were then added on the second mesh to generate the third surface mesh for EDS scenario (Figure 3), which consists of 2499 nodes and 4920 triangular elements. The fourth surface mesh for CDS scenario, which has 2451 nodes and 4799 triangular elements, was generated using the DEM in 2010 with 5 m horizontal resolution. The discretization of all the surface meshes varied from 50 m along the boundary to 20 m along the gullies and 5 m on the five check dams. 30 subsurface layers were added below the surface mesh. Considering that the 0.5 m deep soil near surface is the most sensitive to land use changes and important for surface hydrology, a constant thickness of 0.05 m was assigned to the first ten sublayers. The second ten sublayers (layer 11 to layer 20) had a uniform thickness of 0.5 m. A base elevation of 976.5 m (the elevation of GDG gully ranges from 992 to 1188 m a.s.l) was set as the bottom of the 3D mesh, creating variable thickness for the third ten layers (layer 21 to layer 30).
The vertical discretization for layer 21 to layer 30 ranged from 1.0 m near gully's outlet to 20.5 m at the boundary upland of the gully. This discretization method both ensured a fine resolution in the hydrologically active areas (i.e., the channel, sedimentary land, and the near-surface soil) and simultaneously saved computation resources in the relatively inactive areas (i.e., the headwater regions and deep spaces).
Using a similar method as Heppner and Loague [24], the initial 3D pressure-head distribution of GDG-gully for PD scenario was generated by conducting a one-year quick-drainage simulation starting from the following condition: where ψ t=0 [L] was the initial pressure head of the simulation for a certain node, z sur f [L] referred to the surface elevation directly above the node and z[L] was the elevation of the node. 1060 m-contour was the typical dividing line of gully and slope. The choice of the initial condition for the one-year simulation was motivated by the fact that the groundwater table of the north-western Loess Plateau is around 5~10 m below surface in gullies and nearly 50~150 m below surface for slopes. This quick-drainage simulation used a synthetic rainfall time-series which only contained several small rainfall events to represent a dry year before 1953. A self-consistent head distribution, obtained from the quick-drainage simulation, was assigned to the first GDG-gully BVP (i.e., PD scenario). The initial conditions of SD scenario and EDS scenario were gleaned from the simulation results of PD scenario and SD scenario, respectively. The initial pressure head values for newly added nodes in SD scenario and EDS scenario were determined as the weighted average values of the nearest eight nodes.
To generate the initial 3D pressure-head distribution for CDS scenario, another quick-drainage simulation for the whole WMG catchment was conducted starting from the following condition: The choice of 0.999 in the first part of Equation (3) generated a high water-table shape along gullies at the beginning of the simulation. Then, the quick-drainage simulation started with no flux applied at the surface and ended when the simulated water table depths matched with observed average water table depths at the well-1 and well-2 (i.e., 8.0 m and 5.0 m, respectively) in Figure 1b. The pressure head values of all subsurface nodes for GDG catchment were extracted from the drainage simulation and assigned to the fourth GDG gully-BVP (i.e., CDS scenario) as initial subsurface conditions. Three subsurface boundary conditions were assigned to the 3D boundary-value problems: (1) impermeable for each lateral face; (2) leaking at the saturated hydraulic conductivity for the basal boundary; (3) a local sink (i.e., head-dependent flux) at the down-gradient face. Specified fluxes for precipitation and evapotranspiration and a critical depth (d = 0 m) at the gully outlet were the three surface boundary conditions. The rainfall time-series spanning from 1953 to 1962 for the former 3 scenarios and from 2010 to 2013 for the CDS scenario were obtained from the precipitation station ( Figure 1b). Figure 4 shows the cumulative rainfall and annual potential evapotranspiration (ET 0 ) for all the simulation years. The FAO56 recommended revised-Penman-Monteith method was used to estimate the daily potential evapotranspiration, using meteorological data such as daily temperature, daily vapor pressure, daily atmospheric pressure, daily solar radiation and wind speed from the meteorological station. The calculated daily potential evapotranspiration was then incorporated into InHM as actual evapotranspiration (ET) using a set of sink functions [22]: was a response function of the saturation of the porous medium and the degree of ponding at the land surface, q E max [LT −1 ] was the potential evapotranspiration rate per unit area estimated by the revised Penman-Monteith method, ψ[L] was the pressure head of the subsurface nodes or water depth of the surface nodes, and A b [L 2 ] was the area associated with the surface water equation.

Soil Parameters
Soils were classified as one layer of surface soil (0~20 cm) and two layers of subsurface soil (20~50 cm and below 50 cm). The surface soil layer was furtherly divided into six types according to land covers. Several soil parameters influencing hydrologic response and sediment transport were determined by field measurements or derived from the literature (Table 4). For example, the saturated hydraulic conductivity (Ksat) values for CDS scenario were measured in April 2014, while Ksat values for the other three scenarios were obtained via model calibration. Derived from the previous studies based on soil texture, the damping coefficient, raindrop turbulence factor, and rainfall intensity exponent were, respectively, set to 600 m −1 , 0.25, 1.6 [25,28]. Manning's roughness coefficient and rainsplash coefficient were obtained from model calibration. Surface erodibility coefficient for the first 3 scenarios was calibrated to 0.050. The surface erodibility coefficients for CDS scenario, lower than calibrated value, were derived from the work by Gao et al. [15] to

Soil Parameters
Soils were classified as one layer of surface soil (0~20 cm) and two layers of subsurface soil (20~50 cm and below 50 cm). The surface soil layer was furtherly divided into six types according to land covers. Several soil parameters influencing hydrologic response and sediment transport were determined by field measurements or derived from the literature (Table 4). For example, the saturated hydraulic conductivity (K sat ) values for CDS scenario were measured in April 2014, while K sat values for the other three scenarios were obtained via model calibration. Derived from the previous studies based on soil texture, the damping coefficient, raindrop turbulence factor, and rainfall intensity exponent were, respectively, set to 600 m −1 , 0.25, 1.6 [25,28]. Manning's roughness coefficient and rainsplash coefficient were obtained from model calibration. Surface erodibility coefficient for the first 3 scenarios was calibrated to 0.050. The surface erodibility coefficients for CDS scenario, lower than calibrated value, were derived from the work by Gao et al. [15] to represent the current surface condition. van Genuchten [29] equation was employed to describe the soil water characteristics and the parameters of the equation (Table 5) for loess soils were derived from infiltration experiments conducted in WMG catchment. According to the soil sample data in the catchment, the median diameter of the soil was set to 0.05 mm to represent the uppermost homogeneous loess soil and a single species particle with a particle density of 2650 kg·m −3 was used for all soil layers for the sediment transport simulation. The soil cohesion coefficient was 0.30 [24] for all surface soil except that of the check dam body, which was assigned a larger cohesion coefficient (i.e., 0.60) due to compaction.  [29].; c Parameter related to the pore-size distribution [29].; d Values were all calibrated in WMG BVP.; e Values were all measured in April 2014. Table 6 compares the observed and simulated peak discharges and the time to peak discharges of water and sediment of the ten rainfall-runoff events. EF values for water discharge and sediment discharge were calculated, separately. The EF values in the four calibration events were all higher than 0.70, and the model produced the best simulation results in event 4, which was characterized as low rainfall intensity and rainfall amount with low water and sediment discharges. The average EF values of water discharge and sediment discharge for the six validation events were all higher than 0.55, illustrating that model performances were acceptable both in water discharge simulation and sediment discharge simulation. The observed and simulated hydrographs ( Figure 5) and sedigraphs ( Figure 6) for the 10 events matched reasonably well, also indicating a good representation of the hydrologically-driven sediment transport processes.  4 were used for calibration, and events 5 to 10 were used for validation.; b Q pk referred to the peak discharge, and Q s , pk was the peak sediment discharge.; c T pk referred to the time to Q pk from the start of the event, and T s , pk was the time to Q s , pk from the start of the event.; d EF h referred to the EF value for hydrologic-response simulation, and EF s was the EF value for sediment-transport simulation.   4 were used for calibration, and events 5 to 10 were used for validation.; b Qpk referred to the peak discharge, and Qs, pk was the peak sediment discharge.; c Tpk referred to the time to Qpk from the start of the event, and Ts, pk was the time to Qs, pk from the start of the event.; d EFh referred to the EF value for hydrologic-response simulation, and EFs was the EF value for sediment-transport simulation.    Table 7 provides the InHM simulated GDG water-balance components for each year of the four scenarios. Precipitation of each year was redistributed into three water-balance components (i.e., outflow, storage change, and evapotranspiration). Surface outflow referred to the surface runoff at the outlet of GDG gully, while subsurface outflow included lateral flow leaving the downstream outlet and vertical drainage. Inspection of Table 7 showed that the percentage of each component in the four scenarios was obviously different:

Water Redistribution in Long-Term Simulations of GDG Catchment
 Surface runoff decreased significantly after the construction of GDG1 dam (i.e., averagely reduced by 12.74%) and remained a low percentage in the following 12 simulation years.  Surface storage responded in two different ways in the four scenarios. Surface storage dramatically increased by 8.94% and 18.37% immediately at the end of the first year in SD and EDS scenarios, indicating a large amount of stormwater during the rainy seasons of 1955 and 1959 was retained behind check dams. Comparing to the large growth in 1955 and 1959, small increases of surface storage occurred in the following years of SD and EDS scenarios (i.e., 1956, 1957, 1958, 1960, and 1962). PD and CDS scenarios had similarly slight surface storage increases (i.e., average changes were +0.90% and +1.53%, respectively) for different reasons: (1) surface runoff left the gully without the interception of dams in PD scenario; (2) surface water was retained behind dams and infiltrated into subsurface quickly in CDS scenario. (See the expanding permeable sedimentary layers with relatively high infiltration rate in Figure 7)  Subsurface outflow increased to 7.17% after the construction of GDG1 dam. The mean annual subsurface outflow was 33.5 mm (7.17%) and 54.0 mm (11.49%) for SD and EDS scenarios, increasing by 4.32% after the construction of four more dams. The subsurface outflow remained stable in CDS scenario, fluctuating slightly between wet years and dry years.  Subsurface storage started increasing significantly in 1956, and fluctuating slightly between dry years and wet years. The check dams' role in transforming from a surface reservoir to subsurface reservoir was obvious when comparing the storage-change tendencies of the surface and subsurface.  Table 7 provides the InHM simulated GDG water-balance components for each year of the four scenarios. Precipitation of each year was redistributed into three water-balance components (i.e., outflow, storage change, and evapotranspiration). Surface outflow referred to the surface runoff at the outlet of GDG gully, while subsurface outflow included lateral flow leaving the downstream outlet and vertical drainage. Inspection of Table 7 showed that the percentage of each component in the four scenarios was obviously different:

Water Redistribution in Long-Term Simulations of GDG Catchment
• Surface runoff decreased significantly after the construction of GDG1 dam (i.e., averagely reduced by 12.74%) and remained a low percentage in the following 12 simulation years.  , 1956, 1957, 1958, 1960, and 1962). PD and CDS scenarios had similarly slight surface storage increases (i.e., average changes were +0.90% and +1.53%, respectively) for different reasons: (1) surface runoff left the gully without the interception of dams in PD scenario; (2) surface water was retained behind dams and infiltrated into subsurface quickly in CDS scenario. (See the expanding permeable sedimentary layers with relatively high infiltration rate in Figure 7) • Subsurface outflow increased to 7.17% after the construction of GDG1 dam. The mean annual subsurface outflow was 33.5 mm (7.17%) and 54.0 mm (11.49%) for SD and EDS scenarios, increasing by 4.32% after the construction of four more dams. The subsurface outflow remained stable in CDS scenario, fluctuating slightly between wet years and dry years.

•
Subsurface storage started increasing significantly in 1956, and fluctuating slightly between dry years and wet years. The check dams' role in transforming from a surface reservoir to subsurface reservoir was obvious when comparing the storage-change tendencies of the surface and subsurface.
• Evapotranspiration, the largest term in GDG water balance, showed a slight decline from PD to CDS scenarios. It could be explained by the enhancement of infiltration in the sedimentary areas, which reduced the residence time of rainfall as surface water. In SD and EDS scenarios, dry years showed higher ET proportions than wet years.  Evapotranspiration, the largest term in GDG water balance, showed a slight decline from PD to CDS scenarios. It could be explained by the enhancement of infiltration in the sedimentary areas, which reduced the residence time of rainfall as surface water. In SD and EDS scenarios, dry years showed higher ET proportions than wet years.

Impacts of Check Dams on Groundwater in Long-Term Simulations of GDG Catchment
To further capture the influences of check dams on groundwater and the sedimentary processes, the water table profiles and surface elevations along the gully of each year were compared (Figures 8-11). The impacts of check dams on the GDG1 dam-controlled reach (i.e., A0-A1 reach in Figure 3) in different stages (i.e., SD, EDS, and CDS) are compared in Figures 8-10. Figure 11 compared the water table changes and sedimentary processes of the two-dam reach (i.e., A1-A2 reach in Figure 3) in EDS and CDS scenarios. Perusal of the water table changes lead to the following results: • The water table rose sharply at 200 m upstream from the dam heel of GDG1 and dropped sharply at the dam toe at the end of 1955 (see the green dash line with square symbols in Figure 8), indicating the formation a "subsurface reservoir". • Surface ponding was found in the first year of SD scenario (also see the green line in Figure 8) and disappeared in the following years.

•
The position where water table started rising moved upstream to the location 320 m away from the dam heel of GDG1 at the end of 1958 (see the orange dash line with diamond symbols in Figure 8) and remained there in EDS scenario (Figure 9), indicating the expansion of "subsurface reservoir" referred above.

•
The water table downstream the dam averagely increased by 2.0 m and 0.6 m from the 1st year to the 4th year in SD and EDS scenarios, respectively. However, the water table in the upstream reach remained stable in the SD scenario but averagely rose by 1.0 m in the EDS scenario (see Figures 8 and 9). • Another difference in water table profiles between SD scenario and EDS scenario was the sensitivity to climate. Comparing to the SD scenario, the water table behind GDG1 showed less extreme changes between wet years and dry years in the EDS scenario (see Figures 8 and 9).

•
The water table profiles in the CDS scenario ( Figure 10) were different from those in the former two scenarios. First, the dam's impact that caused the water table rising sharply disappeared in CDS scenario. Second, the water table along the gully had risen by 3~5 m both in the upstream reach and in the downstream reach. Third, the water table at the dam toe was relatively high and might expose as springs, which was actually observed during the field survey in April 2014.

•
The water table in front of GDG2 and GDG4 in EDS scenario showed similar profile to the water table in front of GDG1 in SD scenario (see the green dash lines with square symbols in Figures 8  and 11). However, the intermediate reach of the two dams (i.e., the reach ranging from 880 m to 980 m away from the outlet) experienced a higher water table lifting (3.0 m averagely), indicating a promotion effect of water table rising by two adjacent check dams (Figure 11).
Water 2019, 11, x; doi: www.mdpi.com/journal/water          Table 8 summarizes the eroded sediment mass of each year and the proportion of different zones. The mean erosion moduli, calculated by dividing the eroded mass to the area of GDG gully, was 233.95 t·ha −1 ·a −1 , 226.18 t·ha −1 ·a −1 , 206.77 t·ha −1 ·a −1 and 121.74 t·ha −1 ·a −1 for the four scenarios, respectively. Table 8 shows that the eroded sediment mass decreased apparently in CDS scenario. The main reason for the decrease was that slope erosion was directly alleviated by terraced farming, which replaced nearly 60% of slope farmlands with terraced farmlands. For example, 1.17 × 10 4 t sediment (39%) was eroded from slope farmlands in 1957 (SD scenario with 354.70 mm precipitation). Compared to a similar rainfall condition in 2010 (CDS scenario with 355.00 mm precipitation), only 0.40 × 10 4 t sediment was eroded from the remaining slope farmlands and the terraced farmlands. Another reason was that gully erosion was indirectly alleviated by the existence of check dams, which formed expanding and elevating sedimentary fields in the channel. For example, the total amount of eroded sediment mass from the two gully-zones (Gully_G, Gully_S) decreased from PD scenario to CDS scenario (Table 8).  Table 9 compares the simulated and measured residual deposition heights (∆h) in different stages. Inspection of the surface elevation changes in Figures 8-11 and Table 9 helped to revisit the sedimentary processes of GDG gully:

Impacts of Check Dams on Sediment Transport in Long-Term Simulations of GDG Catchment
• Deposition in front of GDG1 occurred quickly in SD scenario (1955)(1956)(1957)(1958). The 10 m high check dam was nearly fully deposited and facing the risk of dam-break at the end of 1958 ( Figure 8).

•
Being heightened by 13 m before the rainy season of 1959, GDG1 was prevented from being over-deposited during the rainy season of 1959 ( Figure 9). Less sediment silted in front of GDG1 in EDS scenario due to the existence of the four newly-built dams. • GDG2 and GDG4, both of which were 8 m high, experienced quick deposition in EDS scenario (see the green solid line with square symbols in Figure 11). GDG3 and BTG, both of which located on the outlet of tributary gullies, also intercepted large amounts of sediment in EDS scenario (Table 9). • Comparing to PD scenario, a sedimentary layer which has an average thickness of 8 m formed in front of GDG1 dam at the beginning of CDS scenario (see the black solid line in Figure 10). The flat sedimentary area with relatively high soil water content was turned into productive croplands. Long-term alternative erosion and sedimentary formed the undulating terrain in the downstream reach of GDG1. More sediment deposited in the upper reach 600 m away from the outlet rather than at the dam heel for GDG1. • GDG1, which formed productive croplands, had a residual deposition height (∆h) of 7.10 m at the end of CDS scenario. However, the other four check dams all had been almost fully-filled and were vulnerable to floods. • Table 9 also shows that most of the simulated ∆h values were higher than the observed values (except for GDG3 in 1962, GDG2 and GDG4 in 2013), indicating that the sediment volumes deposited in front of dams were underestimated in most simulations.  Table 6 show underestimations of peak sediment discharges, which were most obvious in event 6 with long-duration high rainfall intensities (i.e., 120-150 mm/h). Most of the simulated residual deposition heights ( h) were larger than the observed ones (8 out of 11), indicating an underestimation of deposited sediment behind check dams (Table 9). For example, the simulated ∆h (2.10 m) was nearly two times the observed one (1.10 m) for GDG3 dam in 2013. The most likely reason is that serious collapses occurred on the hillslopes near the two dams. Soils used for constructing check dams were obtained from nearby hillslopes on the two sides of check dams, making the slopes steeper and easier to induce collapse or landslide. According to the inquiries from local farmers, there were several collapses in WMG catchment during the large floods in 27 June 2013 (sediment yield 16,318 tons), one of which occurred on the hillslope near GDG3 dam. Another reason might be the construction of the road on the dam crest. Extra soil from dam crest was poured into the sedimentary area during road construction.

Figures 5 and 6 and
Except for hydraulic erosion, gravity erosion, which usually occurred in the form of landslide or collapse, could dramatically increase sediment yields in a single storm. Since the sediment transport component of InHM was first developed to apply in the hydraulic-erosion-dominating catchments, the sediment transport processes induced by landslides or collapses are not supported yet. Inspiringly, several studies have recently shown the potential of physics-based simulations to forecast landslides and collapses from hydrogeological perspectives (i.e., subsurface fluid-pressures changes in failure-prone locations) [30][31][32]. Future works will be focused on combining the two erosion processes, both of which are related to hydrologic responses, to more accurately predict sediment yields and estimate the sedimentary processes in the Loess Plateau.

Conclusions
A comprehensive physics-based model InHM, after model calibration and validation, was employed in a small gully catchment with a more than 50-year-old check dam system. The simulated residual deposition heights reasonably match with the observed values, indicating the ability of InHM in check dam planning and management. The impacts of check dams on the hydrological response and landforms were investigated and the results are summarized as follows.
(1) Check dams do change the water redistribution in catchment-scale. GDG1 check dam near the gully outlet can effectively reduce surface outflow and increase subsurface outflow. The four other check dams located in the upstream can protect GDG1 from floods risk and promote the redistribution of water. (2) The construction of check dam can significantly change the water table profile along the gully.
A surface reservoir behind the dam will be formed in the first few years and gradually transformed into a subsurface reservoir, resulting in relatively high soil water content in the sedimentary areas. After more than 50 years of operation, the water table along the gully has been averagely elevated by 3-5 m. Adjacent check dams have a promoting effect on water table rising.
(3) Check dams intercept surface water and force sediment in the water to deposit. Gully erosion can be alleviated indirectly due to the formation of expanding sedimentary areas.
The simulations reported herein demonstrate that physics-based simulation can provide a framework for better understanding the impacts (both on hydrological response and landform evolution) of check dams in the Loess Plateau. The study is like a "revisit" to the hydrologic and geomorphic changes that occurred after the construction of these dams, and a prediction of what will happen after long-term operation of the dam system, which could be a useful reference in guiding future check dam construction and management. However, as with previous InHM simulations studying the impacts of human activities on the environment, this study also suffered from the difficulty of validation because of lacking observation data that are not only accurate and credible, but also of the correct kind for distributed simulations [22,24]. Although the residual deposition height (∆h) can be used as a reference to prove the simulated sedimentary processes, more detailed measurements on catchment-scale erosion/deposition are needed to aid future physics-based distributed simulations of hydrologic-response-driven geomorphic evolution processes. Furthermore, gravity erosion should be considered in future InHM simulations and more integrated long-term continuous observations such as the groundwater table distribution along the gully and subsurface outflow rates are needed to further guide the search for a comprehensive understanding of hydrologic responses.  The calculated daily potential evapotranspiration was then incorporated into InHM as actual evapotranspiration (ET) using a set of sink functions [21]: where Q E b [L 3 T −1 ] represents the volumetric evapotranspiration rate, α(ψ)[-] is a response function of the saturation of the porous medium and the degree of ponding at the land surface, q E max [LT −1 ] is the potential evapotranspiration rate per unit area estimated by the revised Penman-Monteith method, ψ[L] is the pressure head of the subsurface nodes or water depth of the surface nodes, and A b [L 2 ] is the area associated with the surface water equation.
The first-order coupling between the surface and subsurface continua, driven by pressure head gradients, occurs via a thin soil layer of thickness, a s in Equation (A3). The control volume finite-element method is employed to discretize the equations in space. Each coupled system of nonlinear equations is solved implicitly using Newton iteration. A more detailed description of the hydrologic-response module of InHM can be found in [19].

Appendix A.2 Sediment-Transport Module
Depth-integrated multiple-species sediment transport, restricted to the surface continuum, is calculated for each sediment species by: C sed max i = 0.05 q s q 3 * g 2 d sed i ψ s (γ sed i − 1) 2 (A11) where α sed i [L 3 T −1 ] is the hydraulic erosion transfer coefficient for species i, C sed max i [L 3 L −3 ] is the concentration at equilibrium transport capacity for species i, q * [LT −1 ] is the local shear velocity, d sed i [L] and γ sed i [-] are the particle diameter and specific gravity, respectively; A[L 2 ] is the area associated with the node in Equation (A12), v sed i [LT −1 ] is the particle settling velocity, ξ[-] is a coefficient related to turbulence in the surface water due to raindrop impact, ϕ[L −1 ] is an erodibility coefficient related to surface properties and texture, and χ i [-] is the particle erodibility factor (ranging from zero to one).