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

**Honglei Tang 1, Qihua Ran <sup>1</sup> and Jihui Gao 2,\***


Received: 7 May 2019; Accepted: 29 May 2019; Published: 2 June 2019

**Abstract:** Check dams are among of the most widespread and effective engineering structures for conserving water and soil in the Loess Plateau since the 1950s, and have significantly modified the local hydrologic responses and landforms. A representative small catchment was chosen as an example to study the influences 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 significantly alter the water redistribution in the catchment and influence the groundwater table in different 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.

**Keywords:** check dam; hydrologic response; sediment transport; InHM; Loess Plateau

#### **1. 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 <sup>×</sup> 1010 m3 of sediments and formed 3340 km2 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–10] at the basin scale. For example, Xu et al. [6] applied Soil and Water Assessment Tool (SWAT) in a 7725 km2 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 discharges and sediment yields at the outlet. However, other hydrologic features such as groundwater table and subsurface storage changes and the sedimentary processes behind the dams may be also influenced by check dams, and sometimes, much more important (e.g., the primary drinking water source for local residents is the groundwater in many catchments). Hydrologic-response change phenomena such as near-surface Ksat increase [11], hillslope-channel decoupling [12], groundwater recharge increase [13] were reported in many catchments with check dams around the world. Water resources in semi-arid areas are precious and the hydrologic-response changes induced by check dams should be better understood. Huang et al. [14] evaluated the impacts of a 30-year-old check dam on water redistribution and the results showed that infiltration was enhanced in the sedimentary field.

Check dams are often constructed as a system in which individual dam is operated in different purposes, and it is important to assess the impacts of a check dam system on the environment. Considering that hydraulic erosion plays an important role in most parts of the Loess Plateau, simulating hydrologic response and sediment transport simultaneously is essential to further understand the influences of check dams. A large amount of sediment eroded from slopes is deposited along gullies due to the interception of check dams. However, the useful lifetime of many small-size check dams was shortened by rapid sedimentary processes behind dams. Rapid sedimentation is expected when a check dam is mainly used as a productive dam to quickly create farmlands but should be avoided for check dams used for preventing floods. The sedimentary processes in dam-controlled gullies, though seldom reported, are important because local people's lives and property are often threatened by dam-break events caused by overtopping floods. In fact, a large number of check dams were destroyed by floods due to inappropriate position or rapid deposition in the Loess Plateau [6]. Thus, an appropriate forecast of the sedimentary processes behind check dams is necessary for future check dam planning and management.

The objectives of this simulation-based study were to capture, as best as possible, the phenomena of hydrologic-response changes and the sedimentary processes caused by a check dam system and to demonstrate that physics-based simulation can be a useful tool for predicting the sedimentary processes induced by the widely launched soil and water conservation measures in future planning and management. The Integrated Hydrology Model (InHM), a physics-based distributed hydrologic model with sediment-transport capabilities, was employed to "revisit" what had happened in the first few years after the construction of check dams and compare the water table changes and sedimentary processes of the current dam-system with early dam-system. This work focused on a small gully catchment with a developed 5-dam system, located in Loess Plateau. Annual continuous simulations were conducted to capture the changes induced by check dams.

#### **2. Materials and Methods**

#### *2.1. Study Site and Data*

Wangmaogou catchment (WMG) is a loess hilly catchment located near the outlet of Wudinghe Basin, China (Figure 1). The area of the catchment is 5.97 km2, and the elevation ranges from 940 to 1188 m a.s.l. Under the impacts of severe soil and water erosion, landscape featured as interlaced gullies are formed gradually, with a gully density of 4.3 km/km2. The geological structure of WMG is relatively simple: (1) the uppermost layer is 20~30 m deep homogeneous loess soil; (2) the second layer is 50~100 m deep red soil, most of which emerges in the gully head; (3) the underlying bedrock is Triassic sandy shale. WMG catchment has a typical semiarid continental climate. The annual potential evapotranspiration is around 800 mm, while the mean annual precipitation is 513 mm, more than 70% of which is received during the rainy season from June to September [15].

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 18,000 t km−<sup>2</sup> 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.

**Figure 1.** (**a**) Location of WMG catchment (the red star); (**b**) check dam distributions and measurement locations of soil characteristics and hydrologic data. Polygons with different colors represent different gullies and the red dash line depicts the boundary of GDG catchment.

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

landuse change of GDG was relatively simple during the 50 years, with a large area of slope farmlands being replaced by terraced farmlands.

**Figure 2.** Landuse types of WMG catchment: (**a**) 1962; (**b**) 2010. The red dash line presents the boundary of GDG catchment.


**Table 2.** Characteristics of the 5 check dams in GDG catchment.


<sup>1</sup> Completion year, all dams were completed before the rainy season.; <sup>2</sup> h is the residual deposition height, measured in April 2014.; <sup>3</sup> GDG1 dam was heightened by 13 m in 1959.

#### *2.2. 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–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.

#### *2.3. 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–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.

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

$$EF = \left[\sum\_{i=1}^{n} \left(O\_i - \overline{O}\right)^2 - \sum\_{i=1}^{n} \left(S\_i - O\_i\right)^2\right] \Big/ \sum\_{i=1}^{n} \left(O\_i - \overline{O}\right)^2\tag{1}$$

where *Oi* was the observed value (water discharge or sediment discharge), *O* was the average value of the observed values, and *Si* 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 *Ksat* for various soil types) to construct a more realistic boundary-value problem for the following simulations.

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

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.

**Figure 3.** 3D mesh for EDS scenario, the inset is a downstream view from position A2.


**Table 4.** Surface parameters used in InHM for the four scenarios.

<sup>a</sup> The letter behind the underline represents different landuse types: G represents grasslands; F represents slope farmlands; D represents check dams; S represents sedimentary fields caused by dams; T represents terraced farmlands; <sup>b</sup> Manning's roughness coefficient [20], calibrated from the first stage.; <sup>c</sup> Immobile water depth [20], identical for all surface zones.; <sup>d</sup> Surface residual saturation [20], identical for all surface zones.; <sup>e</sup> Average height of non-discretized micro-topography [20], identical for all surface zones.; <sup>f</sup> Rainsplash coefficient [18], calibrated from the first stage.; <sup>g</sup> Rainsplash depth dampening factor [27], identical for all surface zones.; <sup>h</sup> Rain intensity exponent [27], identical for all surface zones.; <sup>i</sup> Rain-induced turbulence coefficient [27], identical for all surface zones.; <sup>j</sup> Surface erodibility coefficient [27], calibrated from the first stage. The values for CDS scenario were derived from Gao et al. [15].; <sup>k</sup> Sedimentary fields are usually used as productive farmlands after averagely 2-year deposition, increasing manning's roughness coefficient to 0.12.

#### *2.4. Model Settings and Parameters*

#### 2.4.1. 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:

$$\psi\_{t=0} = \begin{cases} \ 0.992 \times z\_{surf} - z, & z\_{surf} \le 1060 \\ -0.34 \times z\_{surf} - z + 1411.92, & z\_{surf} > 1060 \end{cases} \tag{2}$$

where ψ*t*=0[*L*] was the initial pressure head of the simulation for a certain node, *zsur 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:

$$\psi\_{t=0} = \begin{cases} \ 0.999 \times z\_{surf} - z, & z\_{surf} \le 1060 \\\ -0.34 \times z\_{surf} - z + 1419.34, & z\_{surf} > 1060 \end{cases} \tag{3}$$

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 (ET0) 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]:

$$Q\_b^E = a(\psi) q\_{\text{max}}^E A\_b \tag{4}$$

where *Q<sup>E</sup> <sup>b</sup>* [*L3T*−*1*] represented the volumetric evapotranspiration rate, <sup>α</sup>(ψ)[–] was a response function of the saturation of the porous medium and the degree of ponding at the land surface, *qE* 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 *Ab*[*L2*] was the area associated with the surface water equation.

**Figure 4.** Measured cumulative rainfall for the 14 simulating years. The inset is the annual potential evapotranspiration (ET0) calculated by the revised Penman-Monteith method.

#### 2.4.2. 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 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−<sup>3</sup> 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.


**Table 5.** Soil parameters used in InHM for the four scenarios.

<sup>a</sup> Residual soil-water content.; <sup>b</sup> Parameter related to the inverse of the air-entry pressure [29].; <sup>c</sup> Parameter related to the pore-size distribution [29].; <sup>d</sup> Values were all calibrated in WMG BVP.; <sup>e</sup> Values were all measured in April 2014.

#### **3. Results and Discussion**

#### *3.1. Results of Calibration and Validation in WMG Catchment*

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.


**Table 6.** Comparisons of observed versus simulated results for 10 rainfall-runoff events in the first simulation stage.

<sup>a</sup> Event number. Events 1 to 4 were used for calibration, and events 5 to 10 were used for validation.; <sup>b</sup> Qpk referred to the peak discharge, and Qs, pk was the peak sediment discharge.; <sup>c</sup> 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.; <sup>d</sup> EFh referred to the EF value for hydrologic-response simulation, and EFs was the EF value for sediment-transport simulation.

**Figure 5.** Observed and simulated hydrographs for the ten events in the first simulation stage. (**1**)–(**4**) were events for calibration, (**6**)–(**10**) were events for validation.

**Figure 6.** Observed and simulated sedigraphs for the ten events in the first simulation stage. (**1**)–(**4**) were events for calibration, (**6**)–(**10**) were events for validation.

#### *3.2. Water Redistribution in Long-Term Simulations of GDG Catchment*

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:


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


**Table 7.** Water redistribution in long-term simulations of GDG catchment.

**Figure 7.** Snapshots of infiltration rate at the time of peak discharge of four events with similar rainfall characteristics extracted from the four scenarios, respectively. (**a**) PD scenario; (**b**) SD scenario; (**c**) EDS scenario; (**d**) CDS scenario.

#### *3.3. 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:


**Figure 8.** Surface elevation changes and water table changes of GDG1 dam-controlled reach during SD scenario. Note that the simulation results in the end of month 11 of each year were used for comparison.

**Figure 9.** Surface elevation changes and water table changes of GDG1 dam-controlled reach during EDS scenario. Note that the simulation results in the end of month 11 of each year were used for comparison.

**Figure 10.** Surface elevation changes and water table changes of GDG1 dam-controlled reach during CDS scenario. Note that the simulation results in the end of month 11 of each year were used for comparison.

**Figure 11.** Surface elevation changes and water table changes of GDG2-GDG4 dam-controlled reach in EDS scenario and CDS scenario. Note that the simulation results in the end of month 11 of each year were used for comparison.

#### *3.4. Impacts of Check Dams on Sediment Transport in Long-Term Simulations of GDG Catchment*

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<sup>−</sup>1, 226.18 t·ha−1·a<sup>−</sup>1, 206.77 t·ha−1·a−<sup>1</sup> and 121.74 t·ha−1·a−<sup>1</sup> 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 <sup>×</sup> 104 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 <sup>×</sup> 10<sup>4</sup> 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 8.** Eroded sediment mass of the four simulation scenarios.

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:



**Table 9.** Observed versus simulated residual deposition height (Δh).

<sup>a</sup> O refers to the observed Δh.; <sup>b</sup> S refers to the simulated Δh.; <sup>c</sup> D refers to Difference = [(Simulated-Observed)/Observed] × 100.

#### *3.5. Model Performance Evaluation: Influence of Gravity Erosion*

Figures 5 and 6 and 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.

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

#### **4. 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.


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

**Author Contributions:** Conceptualization, H.T.; methodology, J.G.; software, Q.R.; validation, H.T.; formal analysis, H.T., J.G.; investigation, H.T.; writing—original draft preparation, H.T.; writing—review and editing, J.G.; supervision, Q.R., J.G.; funding acquisition, Q.R.

**Funding:** This study was funded by the National Key Research and Development Program of China (2016YFC0402404), and the National Natural Science Foundation of China (51679209, 51379184).

**Acknowledgments:** Special thanks were given to Suide Soil and Water Conservation station for providing essential field data.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A The Integrated Hydrologic Model (InHM)**

#### *Appendix A.1 Hydrologic-Response Module*

The 3D subsurface flow in variably saturated porous medium is estimated in InHM by:

$$\nabla \cdot f^a \stackrel{\rightharpoonup}{q} \pm q^b \pm q\_{ps}^c = f^v \frac{\partial \phi S\_w}{\partial t} \tag{A1}$$

$$
\overrightarrow{q} = -k\_{rw} \frac{\rho\_w \mathcal{g}}{\mu\_w} \overrightarrow{k} \nabla(\psi\_p + z) \tag{A2}
$$

where *<sup>f</sup> <sup>a</sup>*[–] is the area fraction related to each continuum, *q* [*LT*−*1*] is the Darcy flux, *qb*[*T*−*1*] is a specified rate source/sink, *qe ps*[*T*−*1*] (equal to −*q<sup>e</sup> sp*) is the rate of water exchange between the porous medium and surface continua, *f <sup>v</sup>*[–] is the volume fraction associated with each continuum, φ[*L3L*−*3*] is porosity, *Sw*[*L3L*−*3*] is the water saturation, *t*[*T*] is time, *krw*[–] is the relative permeability, ρ*w*[*ML*−*3*] is the density of water, *<sup>g</sup>*[*LT*−*2*] is the gravitational acceleration, <sup>μ</sup>*w*[*ML*<sup>−</sup>*1T*−*1*] is the dynamic viscosity of water, *k* [*L2*] is the intrinsic permeability vector, ψ*p*[*L*] is the pressure head, and *z*[*L*] is the elevation head.

The diffusion wave approximation to the depth-integrated shallow wave equations is employed in InHM to estimate the 2D transient surface flow on the land surface (both overland and open channel), with the conservation of water on the land surface expressed by:

$$\nabla \cdot \psi\_s^{\text{multi}} \stackrel{\rightharpoonup}{q}\_s \pm a\_s q^b \pm a\_s q\_{sp}^e = \frac{\partial (S\_{w\_s} H\_t + \psi\_{im})}{\partial t} \tag{A3}$$

$$\overrightarrow{\dot{q}}\_s = -\frac{\left(\psi\_s^{\text{mobile}}\right)^{2/3}}{\overrightarrow{n}\,\Phi^{1/2}}\nabla(\psi\_s + z)\tag{A4}$$

where <sup>ψ</sup>*s*[*L*] is the surface water depth, *<sup>q</sup> <sup>s</sup>*[*LT*−*1*] is the surface water velocity calculated by a two-dimensional form of the Manning's equation, *as*[*L*] is the surface coupling length scale, *qb*[*T*−*1*] is the source/sink rate (i.e., rainfall/evaporation), *q<sup>e</sup> sp*[*T*−*1*] is the rate of water exchange between the surface continua and porous medium, *Ht*[*L*] is the average height of non-discretized surface microtopography, *n*[*TL*−*1*/*3*] is the Manning's surface roughness tensor, and Φ[–] is the friction (or energy) slope; ψ*mobile <sup>s</sup>* [*L*] and ψ*im*[*L*] refer to the water depth exceeding and held in depression storage, respectively.

The calculated daily potential evapotranspiration was then incorporated into InHM as actual evapotranspiration (ET) using a set of sink functions [21]:

$$Q\_b^E = a(\psi) q\_{\text{max}}^E A\_b \tag{A5}$$

where *Q<sup>E</sup> <sup>b</sup>* [*L3T*−*1*] represents the volumetric evapotranspiration rate, <sup>α</sup>(ψ)[–] is a response function of the saturation of the porous medium and the degree of ponding at the land surface, *qE* 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 *Ab*[*L2*] 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, *as* 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:

$$\frac{\partial \mathbb{C}\_{\rm sed}}{\partial t} = -\nabla \cdot \overrightarrow{q}\_{s} \mathbb{C}\_{\rm sed} + \frac{1}{V\_{\rm uv}} (\mathbf{e}\_{\rm sd} + \sum\_{j=1}^{\rm BC} q\_{s\_{j}}^{b} \mathbb{C}\_{\rm sd\_{j}}^{\*}) \tag{A6}$$

$$
\varepsilon\_{\rm sed} = \varepsilon\_{\rm s} + \varepsilon\_{\rm l} \tag{A7}
$$

where *Csed*[*L3L*−*3*] is volumetric sediment concentration, *<sup>q</sup> <sup>s</sup>*[*LT*−*1*] is the depth-averaged surface water velocity, *Vw*[L3] is the volume of water at the node, *esed*[*L3T*−*1*] is the volumetric rate of soil erosion and/or deposition, *qb sj* [*L3T*−*1*] is the rate of water added/removed via the *j*th boundary condition, *C*∗ *sedj* [*L3L*−*3*] is the sediment concentration of the water added/removed via the *j*th boundary condition, *BC* is the total number of boundary conditions. The net erosion rate is the sum of the rainsplash erosion rate *es*[*L3T*−*1*] and the hydraulic erosion rate *eh*[*L3T*−*1*]. The rainsplash erosion rate of each species in Equation (A7) is calculated as:

$$\mathfrak{e}\_{\mathfrak{s}\_i} = \begin{cases} \sigma\_i \mathfrak{c}\_f F(\psi\_s) (\cos(\theta) \cdot r)^b A\_{3D} & , \quad q > 0 \\ 0 & , \quad q < 0 \end{cases} \tag{A8}$$

$$F(\psi\_s) = \exp(-c\_d \psi\_s) \tag{A9}$$

where σ*i*[–] is the source fraction of species *i*, *c <sup>f</sup>*[*(TL*−*1) <sup>b</sup>*−*1*] is the rainsplash coefficient, *b*[–] is the rainfall intensity exponent, θ[–] is the angle of the element from horizontal, *r*[*LT*−*1*] is the rainfall intensity, *A*3*D*[*L2*] is the three-dimensional area associated with the node, and *q*[*LT*−*1*] is the sum of rainfall intensity and infiltration intensity; *F*(ψ*s*)[–] is a damping function, related to surface water depth ψ*s*[*L*] and the surface water damping-effectiveness coefficient *cd*[*L*−*1*], to represent the reduction in splash erosion with increasing surface water depth. The hydraulic erosion rate in Equation (A7) estimated as:

$$
\alpha\_{l\_i} = \alpha\_{\text{seld}\_i} (\mathbb{C}\_{\text{scd}\_{\text{max}\_i}} \sigma\_i - \mathbb{C}\_{\text{scd}\_i}) \tag{A10}
$$

$$\mathcal{C}\_{\text{sol}\_{\text{max}\_i}} = 0.05 \frac{\overrightarrow{q}\_s q\_\*^3}{g^2 d\_{\text{scd}\_i} \psi\_s \left(\mathcal{V}\_{\text{scd}\_i} - 1\right)^2} \tag{A11}$$

α*sedi* = *A* ⎧ ⎪⎪⎨ ⎪⎪⎩ 2*vsedi* <sup>ξ</sup> , *Csed* <sup>&</sup>gt; *Csed*max*<sup>i</sup>* (*erosion*) ϕ *<sup>q</sup> <sup>s</sup>*ψ*s*χ*<sup>i</sup>* , *Csed* <sup>&</sup>lt; *Csed*max*<sup>i</sup>* (*deposition*) (A12)

where α*sedi* [*L3T*−*1*] is the hydraulic erosion transfer coefficient for species *<sup>i</sup>*, *Csed*max*<sup>i</sup>* [*L3L*−*3*] is the concentration at equilibrium transport capacity for species *<sup>i</sup>*, *<sup>q</sup>*∗[*LT*−*1*] is the local shear velocity, *dsedi* [*L*] and γ*sedi* [–] are the particle diameter and specific gravity, respectively; *A*[*L2*] is the area associated with the node in Equation (A12), *vsedi* [*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).

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
