10. Dong, Y.Q.; Zhuang, X.H.; Lei, T.W.; Yin, Z.; Ma, Y.Y. A method for measuring erosive flow velocity with simulated rill. *Geoderma* **2014**, *232*–234, 556–562. [CrossRef]

11. Bewket, W.; Sterk, G. Assessment of soil erosion in cultivated fields using a survey methodology for rills in the Chemoga watershed, Ethiopia. *Agric. Ecosyst. Environ.* **2003**, *97*, 81–93. [CrossRef]

12. Aksoy, H.; Unal, N.E.; Cokgor, S.; Gedikli, A.; Yoon, J.; Koca, K.; Inci, S.B.; Eris, E. A rainfall simulator for laboratory-scale assessment of rainfall-runoff-sediment transport processes over a two-dimensional flume. *Catena* **2012**, *98*, 63–72. [CrossRef]

13. Zheng, F.L. A research on method of measuring rill erosion amount. *Bull. Soil Water Conserv.* **1989**, *9*, 41–49.

14. Casalí, J.; Loizu, J.; Campo, M.A.; De Santisteban, L.M.; Álvarez-Mozos, J. Accuracy of methods for field assessment of rill and ephemeral gully erosion. *Catena* **2006**, *67*, 128–138. [CrossRef]

15. Shen, H.; Zheng, F.; Wen, L.; Han, Y.; Hu, W. Impacts of rainfall intensity and slope gradient on rill erosion processes at loessial hillslope. *Soil Tillage Res.* **2016**, *155*, 429–436. [CrossRef]

16. Balaguer-Puig, M.; Marqués-Mateu, Á.; Lerma, J.L.; Ibáñez-Asensio, S. Estimation of small-scale soil erosion in laboratory experiments with Structure from Motion photogrammetry. *Geomorphology* **2017**, *295*, 285–296. [CrossRef]

17. Berger, C.; Schulze, M.; Rieke-Zapp, D.; Schlunegger, F. Rill development and soil erosion: A laboratory study of slope and rainfall intensity. *Earth Surf. Proc. Land.* **2010**, *35*, 1456–1467. [CrossRef]

18. Heng, B.C.P.; Chandler, J.H.; Armstrong, A. Applying close range digital photogrammetry in soil erosion studies. *Photogramm. Rec.* **2010**, *25*, 240–265. [CrossRef]

19. Nouwakpo, S.; Huang, C.-H.; Frankenberger, J.; Bethel, J. A simplified close range photogrammetry method for soil erosion assessment. In *Proceedings of the 2nd Joint Federal Interagency Conference, Las Vegas, NV, USA, 27 June–1 July 2010*.

20. Vinci, A.; Brigante, R.; Todiscoa, F.; Mannocchi, F.; Radicioni, F. Measuring rill erosion by laser scanning. *Catena* **2015**, *124*, 97–108. [CrossRef]

21. Afana, A.; Solé-Benet, A.; Pérez, J.L. Determination of soil erosion using laser scanners. In *Proceedings of the 19th World Congress of Soil Science, Soil Solutions for a Changing World*, Brisbane, Australia, 1–6 August 2010.

22. Longoni, L.; Papini, M.; Brambilla, D.; Barazzetti, L.; Roncoroni, F.; Scaioni, M.; Ivanov, V.I. Monitoring riverbank erosion in mountain catchments using terrestrial laser scanning. *Remote Sens.* **2016**, *8*, 241. [CrossRef]

23. Dabek, P.B.; Patrzalek, C.; Ćmielewski, B.; Żmuda, R. The use of terrestrial laser scanning in monitoring and analyses of erosion phenomena in natural and anthropogenically transformed areas. *Cogent Geosci.* **2018**, *4*, 1437684. [CrossRef]

24. Thoma, D.P.; Gupta, S.C.; Bauer, M.E.; Kirchoff, C.E. Airborne laser scanning for riverbank erosion assessment. *Remote Sens. Environ.* **2005**, *95*, 493–501. [CrossRef]

25. Filin, S.; Baruch, A.; Morik, S.; Avni, Y.; Marco, S. Characterization of land degradation processes using airborne laser scanning. In *Proceedings of the International Archives of the Photogrammetry, Remote Sensing and Spatial Information Science*, Kyoto, Japan, 9–12 August 2010.

26. ASTM. *Standard Test Methods for Specific Gravity of Soil Solids by Water Pycnometer*; D854-14; ASTM: West Conshohocken, PA, USA, 2014.

27. ASTM. *Standard Test Method for Sieve Analysis of Fine and Coarse Aggregate*; ASTM C136/C136M-14; ASTM: West Conshohocken, PA, USA, 2015.

28. ASTM. *Standard Test Methods for Determining the Amount of Material Finer than 75-µm (No. 200) Sieve in Soils by Washing*; D1140-14; ASTM: West Conshohocken, PA, USA, 2014.

29. ASTM. *Standard Test Method for Particle-Size Analysis of Soils*; ASTM D422-63; ASTM: West Conshohocken, PA, USA, 2007.

30. ASTM. *Standard Test Methods for Laboratory Determination of Water (Moisture) Content of Soil and Rock by Mass*; D2216-10; ASTM: West Conshohocken, PA, USA, 2010.

31. ASTM. *Standard Practices for Preserving and Transporting Soil Samples*; D4220/D4220M-14; ASTM: West Conshohocken, PA, USA, 2014.

32. Eijkelkamp Soil & Water. Available online: https://en.eijkelkamp.com/products/field-measurement-equipment/rainfall-simulator.html (accessed on 3 July 2018).

33. Nciizah, A.D.; Wakindiki, I.I.C. Rainfall pattern effects on crusting, infiltration and erodibility in some South African soils with various texture and mineralogy. *Water SA* **2014**, *40*, 57–64. [CrossRef]
34. Salles, C.; Poesen, J.; Sempere-Torres, D. Kinetic energy of rain and its functional relationship with intensity. *J. Hydrol.* **2002**, *257*, 256–270. [CrossRef]

35. Hung, C.-Y.; Capart, H. Rotating laser scan method to measure the transient free-surface topography of small-scale debris flows. *Exp. Fluids* **2013**, *54*, 1544. [CrossRef]

36. Hung, C.-Y. Relation between Debris Flow Rheology and Fan Deposit Morphology Investigated Using Small-Scale Experiments. Master's Thesis, National Taiwan University, Taipei, Taiwan, June 2011.

37. ASTM. *Standard Test Methods for Liquid Limit, Plastic Limit, and Plasticity Index of Soils*; D4318-10; ASTM: West Conshohocken, PA, USA, 2010.

38. Bu, C.-F.; Wu, S.-F.; Yang, K.-B. Effects of physical soil crusts on infiltration and splash erosion in three typical Chinese soils. *Int. J. Sediment Res.* **2014**, *29*, 491–501. [CrossRef]

39. Di Prima, S.; Concialdi, P.; Lassabatere, L.; Angulo-Jaramillo, R.; Pirastru, M.; Cerdà, A.; Keesstra, S. Laboratory testing of Beerkan infiltration experiments for assessing the role of soil sealing on water infiltration. *Catena* **2018**, *167*, 373–384. [CrossRef]

40. Lin, L.-L.; Tung, H.-P.; Wan, H.-S. *Handbook of Soil Physics Experiments*; Department of Soil and Water Conservation, National Chung Hsing University: Taichung City, Taiwan, 2002.

41. Shen, H.O.; Zheng, F.L.; Wen, L.L.; Lu, J.; Han, Y. An experimental study on rill morphology at loess hillslope. *Acta Ecol. Sin.* **2014**, *34*, 5514–5521. [CrossRef]

42. Niu, Y.B.; Gao, Z.L.; Li, Y.H.; Luo, K.; Yuan, X.H.; Du, J.; Zhang, X. Rill morphology development of engineering accumulation and its relationship with runoff and sediment. *Trans. Chin. Soc. Agric. Eng.* **2016**, *32*, 154–161.

43. He, J.-J.; Lu, Y.; Gong, H.-L.; Cai, Q.-G. Experimental study on rill erosion characteristics and its runoff and sediment yield process. *Shuili Xuebao* **2013**, *44*, 398–405.

44. Chen, J.J.; Sun, L.Y.; Cai, C.F.; Liu, J.T.; Cai, Q.G. Rill erosion on different soil slopes and their affecting factors. *Acta Pedol. Sin.* **2013**, *50*, 281–288.

45. Liu, Q.Q.; Li, J.C.; Chen, L.; Xiang, H. Dynamics of overland flow and soil erosion (II): Soil erosion. *Adv. Mech.* **2004**, *34*, 493–506.

46. Yao, C.; Lei, T.; Elliot, W.J.; McCool, D.K.; Zhao, J.; Chen, S. Critical conditions for rill initiation. *Trans. ASABE* **2008**, *51*, 107–114. [CrossRef]

<figure><img src="image_72.png" /></figure>

© 2018 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/).
<figure><img src="image_73.png" /></figure>

Article

# Effects of Climate Change and Human Activities on Soil Erosion in the Xihe River Basin, China

Shanshan Guo, Zhengru Zhu and Leting Lyu *

College of Urban and Environment, Liaoning Normal University, Dalian 116029, China;
shanshanguo126@126.com (S.G.); Zhengruzhu@gmail.com (Z.Z.)

* Correspondence: lvleting@lnnu.edu.cn; Tel.: +86-411-8425-8364

Received: 27 June 2018; Accepted: 13 August 2018; Published: 15 August 2018

<figure><img src="image_74.png" /></figure>

check for
updates

**Abstract:** Climate change and human activities are the major factors affecting runoff and sediment load. We analyzed the inter-annual variation trends of the annual rainfall, air temperature, runoff and sediment load in the Xihe River Basin from 1969–2015. Pettitt's test and the Soil and Water Assessment Tool (SWAT) model were used to detect sudden changes in hydro-meteorological variables and simulate the basin hydrological cycle, respectively. According to the simulation results, we explored spatial distribution of soil erosion in the watershed by utilizing ArcGIS10.0, analyzed the average soil erosion modulus by different types of land use, and quantified the contributions of climate change and human activities to runoff and sediment load in changes. The results showed that: (1) From 1969–2015, both rainfall and air temperature increased, and air temperature increased significantly (p < 0.01) at 0.326 °C/10 a (annual). Runoff and sediment load decreased, and sediment load decreased significantly (p < 0.01) at 1.63 × 10⁵ t/10 a. In 1988, air temperature experienced a sudden increased and sediment load decreased. (2) For runoff, R² and Nash and Sutcliffe efficiency coefficient (Ens) were 0.92 and 0.91 during the calibration period and 0.90 and 0.87 during the validation period, for sediment load, R² and Ens were 0.60 and 0.55 during the calibration period and 0.70 and 0.69 during the validation period, meeting the model's applicability requirements. (3) Soil erosion was worse in the upper basin than other regions, and highest in cultivated land. Climate change exacerbates runoff and sediment load with overall contribution to the total change of –26.54% and –8.8%, respectively. Human activities decreased runoff and sediment load with overall contribution to the total change of 126.54% and 108.8% respectively. The variation of runoff and sediment load in the Xihe River Basin is largely caused by human activities.

**Keywords:** climate change; human activities; soil erosion; SWAT model; Xihe River Basin

## 1. Introduction

Soil erosion can be differentiated into natural erosion, which is a geographical phenomenon affecting the surface of the Earth, and human activities-accelerated erosion. The speed of soil erosion arising from human activities significantly exceeds the rate of soil formation. Soil erosion leads to the loss of soil nutrients, decline in soil fertility [1], deterioration of the ecological environment, and changes to the ecological landscape [2]. Soil erosion can also result in the accumulation of sediment in river channels, rising of river beds, sediment accumulation in reservoirs, and increased risk of flooding [3]. The dangers of soil erosion are compelling, and it is among the most significant global environmental concerns [4,5]. Scholars have conducted extensive research in various regions around the world with the aim of distinguishing the factors that influence soil erosion and developing new control methods [6–14]. The Soil and Water Assessment Tool (SWAT) [15,16] is a large-scale watershed hydrological model developed by the Agricultural Research Center of the US Department of Agriculture. The model contains three sub-models, namely the hydrological process, soil erosion,
and pollution load and can simulate runoff and sediment load, nonpoint source pollution, and pesticide
diffusion [17]. The tool is widely applied today in studying the water and sediment simulation and
soil erosion of the basin [18–23]. Using the SWAT model, Wu et al. [24] quantified the impacts of
climate change and human activities on runoff and sediment in the Yanhe River Basin and found that
climate change and human activities are the main factors affecting runoff and sediment, respectively.
Alighalehbabakhani et al. [25] used the SWAT model to evaluate sediment accumulation rates within
the Lake Rockwell and Ballville reservoirs and concluded that sediment accumulation efficiency
is closely related to the characteristics of the inflow, reservoir inflow time, and reservoir size. Yen
et al. [26] built the Arroyo Colorado watershed SWAT model and studied the effects of the four
different sediment transport functions to predict sediment. They demonstrated that without complete
observation data, four sediment transport equations might produce similar sediment yield and the
default model might not necessarily simulate the best results. At present quantitatively distinguish the
effects of climate change and human activities on runoff and soil erosion has attracted much attention
from scholars and identifying the impacts of each factor has important applications value in controlling
soil erosion [27].

China is one of the countries with the most serious soil erosion in the world. The soil erosion
types are various, and the erosion process is complex. It mainly occurs in the Loess Plateau of the
middle and upper reaches of the Yellow River, the hilly areas in the middle and upper reaches of the
Yangtze River, and the Northeast Plain area. According to the second national remote sensing survey
of soil erosion in 2001, the area of water and wind erosion in China is 3.569 million km², accounting for
37% of the land area. These phenomena seriously affect China's ecological environment construction,
hindering social and economic development.

The continuously increasing population of the basin and the rapid development of companies
such as Benxi Steel Group and Bei Steel Group along the Xihe River have increased domestic as well
as industrial and agricultural water use, resulting in a reduction of runoff in the basin. In addition,
the accelerated urbanization rate has increased the urban land area, the ground is hardened and the
soil is not easily eroded. From 2006–2008, the Nanfen District of the Xihe River Basin implemented
numerous control measures including automation of dams, flood dikes, and dredging, therefore
sediment discharge from the upper stream into the channel was intercepted and cleared, resulting in
a significant reduction in sediment load [28]. The Nanfen tailings pool has also achieved ecological
benefits. The forest surface has increased, effectively decreased runoff and soil erosion in the tailings
reservoir area. In 2012, to restore the natural environment of the Xihe River Basin, the Nanfen District
implemented the “Blue Water Project” for Xihe River and its tributaries, including river closure
and protection. A land area of 1299 ha of land was preserved, and 264 ha of forest were planted,
which effectively reduced soil erosion and sediment transported downstream.

The overall objective of this study is to investigate changes in runoff and sediment load from
climatic change and human activities in Xihe River Basin of Northeast China. The primary objectives
are: (1) to identify abrupt change points in annual hydro-meteorological series from 1969–2015; (2) to
apply a SWAT model in the Xihe River Basin; (3) to evaluate the impacts of climate change and
human activities on runoff and soil erosion. The results can provide valuable reference information for
controlling soil erosion in this basin.

**2. Materials and Methods**

2.1. *Research Area Overview*

The Xihe River Basin is in the eastern region of Liaoning Province (China), centered in Benxi City and is a tributary of the Taizihe River, which is located between 40°47′–41°16′ N and 123°32′–123°59′ E. The diagram of the research area is shown in Figure 1. Its source lies in Fengcheng County, Dandong City, Liaoning Province and the basin has a total area of 1047 km². It is the main source of water for industrial, agricultural, and residential consumption in Benxi City. The area is dominated by
mountains ranging from 75 to 1157 m a.s.l (above sea level). The characteristic of the area is high lands in the east and south, and low lands in the center and northwest. The east belongs to the Changbai Mountains. The region is in the semi-humid and semi-arid monsoon climate zones. The months with the highest and lowest temperature are July and January respectively, with an annual average temperature and annual rainfall of 8.4 °C and 772.18 mm, respectively. There are five types of land uses in the basin: woodland, cultivated land, towns, water areas, and unused land. The soil types are Calcaric Regosols, Haplic Phaeozem, Haplic Luvisols, Gleyic Luvisol, and Eutric Cambisols.

<figure><img src="image_75.png" /><figcaption>Figure 1. Research Area Overview.</figcaption></figure>

## 2.2. Data Sources and Research Methodology

### 2.2.1. Data Sources

SWAT requires topographical, land use, soil, and meteorological inputs. Thematic layers and climatic data were developed using the sources specified in Table 1. The land use map and soil map are shown in Figure 2.

<table><caption>Table 1. Input data of Soil and Water Assessment Tool (SWAT) model in the Xihe River Basin.</caption><thead><tr><th>Data Type</th><th>Data Input</th><th>Accuracy</th><th>Data Source</th></tr></thead><tbody><tr><td rowspan="3">Spatial data</td><td>Digital Elevation Model (DEM)</td><td>30 m</td><td>Geospatial Data Cloud.<br><a href="http://www.gscloud.cn/">http://www.gscloud.cn/</a><br>16 August 2017</td></tr><tr><td>Land use map</td><td>30 m</td><td>Geospatial Data Cloud Remote Sensing Map</td></tr><tr><td>Soil type map</td><td>1000 m</td><td>Harmonized World Soil Database, HWSD</td></tr><tr><td>Meteorological data</td><td>Precipitation, Temperature, Wind, Solar radiation, Humidity</td><td>Daily<br>1969–2015</td><td>China Meteorological Data Center. <a href="http://data.cma.cn/">http://data.cma.cn/</a><br>2 October 2017</td></tr><tr><td>Hydrological data</td><td>Runoff, Sediment</td><td>Daily</td><td>Hydrological Yearbook</td></tr><tr><td>Soil data</td><td>Soil moisture density, Effective water content</td><td>1969–2015</td><td>Harmonized World Soil Database, HWSD</td></tr></tbody></table>
<figure><img src="image_76.png" /><figcaption>Figure 2. Land use map (a) and Soil map (b).</figcaption></figure>

### 2.2.2. Analysis of the Changed Characteristics of Hydro-Meteorological Elements

In this study, linear regression analysis and the Pettitt sudden change point detection method [29] were used to analyze the inter-annual variation trends and sudden changes in rainfall, air temperature, runoff, and sediment load in the Xihe River Basin from 1969–2015. The Pettitt sudden change detection method is a non-parametric test that can quantify the statistically significant level of change points while varying the hydrological and meteorological factors, which has been widely used by scholars [30–32]. This method uses the Mann–Whitney statistical function $U_{t,N}$, and considers the samples $x_1, x_2, \dots, x_t$ and $x_{t+1}, \dots, x_N$ to be independently and identically distributed. $N$ represents the sample size. $U_{t,N}$ is calculated as:

$$U_{t,N} = U_{t-1,N} + \sum_{i=1}^{N} \text{sgn}(x_t - x_i) (t = 2, \dots, N) \quad (1)$$

If time $t$ is satisfied by $k_t = \max_{1 \le t < N} |U_{t,N}|$ then $t$ represents a sudden change point and the formula for the significant level of $p$ of the change point is as follows:

$$p = 2 \exp\left(\frac{-6k_t^2}{N^3 + N^2}\right) \quad (2)$$

If $p \le 0.05$ then the detected variations are considered statistically significant.

According to the sudden change analysis, the time series of runoff and sediment load from 1969–2015 were divided into two periods. The period prior to the sudden change was considered as ‘the baseline period’ and used as the calibration period. The period after the sudden change was considered as ‘the impacted period’ and used as the validation period [33–35].

### 2.2.3. SWAT Model Set Up

The SWAT model is a physically based distributed model that can evaluate the impact of land management practices on water resources, sediments and agricultural chemical yields in large and complex basins. The major components of SWAT include climate, hydrology, sediment movement, crop growth, nutrient cycling and agricultural management. SWAT model divide the basin into numerous
sub-basins based on topography, and then further subdivided into a series of Hydrologic Response Units (HRU) with unique soil, land cover and slope characteristics [16].

In this study, we chose the Soil Conservation Service curve number (SCS-CN) procedure that SWAT model provides [36] to estimate surface runoff. SWAT calculates the peak runoff rate using a modified rational method. The potential evapotranspiration is estimated by using the Penman-Monteith method [37]. It is calculated by the following formula:

$$ \lambda E = \frac{\Delta \cdot (H_{net} - G) + \rho_{air} \cdot C_p \cdot (e_z^o - e_z)/r_a}{\Delta + \gamma \cdot (1 + r_c/r_a)} \quad (3) $$

where $\lambda E$ is latent heat flux density (MJ/(m²·d)); $E$ is evaporation rate (mm/d); $\Delta$ is saturation vapor pressure-the slope of the temperature curve; $H_{net}$ is net radiation (MJ/(m²·d)); $G$ is ground heat flux density (MJ/(m²·d)); $\rho_{air}$ is air density (kg/m³); $C_p$ is specific heat at fixed air pressure (MJ/kg·°C); $e_z^o$ is saturated vapor pressure at z height (kPa); $e_z$ is the vapor pressure at z height (kPa); $\gamma$ is psychrometric constant (kPa/°C); $r_c$ is the impedance of the vegetation canopy (s/m); $r_a$ is the diffusion impedance of the air layer (s/m).

The SWAT model uses the Modified Universal Soil Loss Equation (MUSLE) to simulate sediment yield [38]. Sediment transport is a function of deposition and degradation, which are determined through comparing the sediment concentration and maximum sediment concentration [39]. Further details of hydrological and sediment transport processes can be found in the SWAT Theoretical Documentation [39].

We used the ArcSWAT 2009 interface to establish and parameterize the model. First, we chose a threshold drainage area of 22.13 km², which is recommended by SWAT model, to delineate the basin. And the basin was divided into 27 sub-basins. Furthermore 960 HRUs were generated according to the thresholds of land use (0%), soil type (0%), and slope class (0%), respectively. The HWSD database provides partial physical attribute data for the soil, the effective field water content (SOL_AWC) and other attributes were calculated using the soil water characteristic software “SPAW”. The minimum infiltration rate was calculated according to the empirical formula introduced in literature [40], and the soil erosion factor USLE_K was calculated according to the basic principle of SWAT [39]. Finally, we set up the SWAT model using daily meteorological data from 1969–2015, and calibrated the parameters using runoff and sediment load data of the baseline period, therefore the natural rainfall-runoff process was simulated.

### 2.2.4. Model Sensitivity Analysis, Calibration and Validation

There are many parameters in SWAT model that affect the results of runoff and sediment load simulations. In order to select the main parameters to improve the usability of the model, the sensitivity analysis of the parameters should be carried out first. We use the sequential uncertainty fitting algorithm (SUFI-2) algorithm [41] in SWAT-CUP 2012 to carry out sensitivity analysis. The ‘baseline period’ was used as the calibration period, and the ‘impacted period’ was used as the validation period. The model performance was assessed by the coefficient of determination ($R^2$) and Nash and Sutcliffe efficiency coefficient (Ens) [42]. $R^2$ is the square of the correlation coefficient, and the closer it is to 1, the higher the degree of agreement between the simulated value and the observed value is. Ens is a statistical measure that determines the relative magnitude of the residual variance compared to the measured data variance, and the closer it is to 1, the closer the simulated value is to the observed value. When $R^2 \ge 0.6$ and $Ens \ge 0.5$, the simulation results may be deemed satisfactory [43].

### 2.2.5. Quantitative Analysis of the Effects of Climate Change and Human Activities on Runoff and Sediment load

Runoff and sediment load are simultaneously affected by climate change and human activities, the total changes in runoff and sediment load is the difference between the baseline period and the impacted period of the observed value, and was calculated as follows:

$$ \Delta Q_t = Q_v - Q_b \quad (4) $$
where ΔQ_t represents total change, while Q_v and Q_h represent the observed annual average values of the baseline period and the impacted period, respectively. Change in total runoff and sediment load can be separated into climate change and human activities, as follows:

$$ \Delta Q_t = \Delta Q_c + \Delta Q_h \qquad (5) $$

where ΔQ_c represents hydrological changes attributable to climate change, and ΔQ_h represents hydrological changes attributable to human activities, calculated as follows:

$$ \Delta Q_c = Q_{vs} - Q_{bs} \qquad (6) $$

where Qvs represents the SWAT model simulated baseline period average annual values, and Qbs represents the SWAT model simulated impacted period average annual values. The difference between Qvs and Qbs represents the amount of runoff and sediment load attributable to climate change. According to Equation (5), the difference between ΔQ_t and ΔQ_c is the amount of change attributable to human activities. The percent contributions of change in runoff and sediment load attributable to climate change and human activities can be calculated as follows:

$$ P_c = \frac{\Delta Q_c}{\Delta Q_t} \times 100\% \qquad (7) $$

$$ P_h = \frac{\Delta Q_h}{\Delta Q_t} \times 100\% \qquad (8) $$

# 3. Results

## 3.1. Temporal Trends and Sudden Changes Analysis of Hydrological and Meteorological Variables

### 3.1.1. Inter-Annual Variation Trends in Hydrological and Meteorological Variables

The variation trends of annual precipitation, temperature, runoff, and sediment load in the Xihe River Basin from 1969–2015 can be found in Figure 3a,d. In the basin, the annual precipitation showed an insignificant increasing trend and the annual runoff showed an insignificant decreasing trend. The annual average precipitation and runoff was 772.18 mm and 3.16 × 10⁸ m³, respectively. The average annual temperature was 8.4 °C and annual temperature showed a significantly increasing trend (p < 0.01) at a rate of 0.326 °C/10 a. The annual average sediment load was 5.02 × 10⁵ t and annual sediment load showed a significantly decreasing trend (p < 0.01) at a rate of 1.63 × 10⁵ t/10 a.

<figure><img src="image_77.png" /><figcaption>Figure 3. Cont.</figcaption></figure>
<figure><img src="image_78.png" /><figcaption>**Figure 3.** Inter-annual variation trend of the annual (a) precipitation, (b) temperature, (c) runoff and (d) sediment load from 1969–2015 in the Xihe River Basin.</figcaption></figure>

### 3.1.2. Hydrological and Meteorological Variables Sudden Changes Analysis

The results of the Pettitt sudden changes analysis of annual precipitation, temperature, runoff, and sediment load in the Xihe River Basin from 1969–2015 can be found in Figure 4a–b. The results showed no sudden changes in precipitation and runoff. The results of the temperature tests $|U_{t,N}|$ showed that the largest year was 1988 ($p < 0.01$). Prior to the change point, the annual average temperature was 7.88 °C, and afterwards it was 8.78 °C. Following the change point, the annual average temperature increased by 11.56%. Sediment load in 1988 also experienced a sudden change ($p < 0.05$). Prior to the change point, the annual average sediment load was $7.89 \times 10^5$ t and after the change point was $3.75 \times 10^5$ t. After the change point, the annual average sediment load decreased by 52.47%. To study the effects of climate change and human activities on runoff and sediment load, the study period was divided into two periods, the baseline period 1971–1987 and the impacted period 1988–2015.

<figure><img src="image_79.png" /><figcaption>**Figure 4.** Results of sudden changes in annual (a) precipitation, (b) temperature, (c) runoff and (d) sediment load from 1969–2015 in the Xihe River Basin.</figcaption></figure>
## 3.2. SWAT Model Simulation Results

The results of parameters sensitivity analysis for runoff and sediment load are showed in Table 2. The most sensitive parameters for runoff and sediment load are the base flow alpha factor (ALPHA_BF) and sediment movement linear compensation coefficient (SPCON), respectively. What needs illustration is that the final range of USLE_P is 0.56–1; this corresponds to the farming practices in the Xihe River Basin. In this watershed, the cultivated land is mainly distributed along the river. Contour tillage with conservative farming methods is most common.

<table><caption>Table 2. Parameters sensitivity analysis results in the Xihe River Basin.</caption><thead><tr><th>Parameter Name and Document Type</th><th>Parameter Description</th><th>Process</th><th>Rank</th><th>Beg. Range</th><th>End Range</th></tr></thead><tbody><tr><td>v_ALPHA_BF.gw</td><td>Base flow α coefficient</td><td></td><td>1</td><td>0.0~1.0</td><td>-0.13~0.62</td></tr><tr><td>r_CN2.mgt</td><td>Runoff curve parameter</td><td></td><td>6</td><td>-0.2~0.2</td><td>-0.06~0.23</td></tr><tr><td>v_SLSUBBSN.hru</td><td>Average grade length</td><td></td><td>9</td><td>10~150</td><td>14.03~22.41</td></tr><tr><td>v_ESCO.hru</td><td>Soil evaporation coefficient</td><td>Runoff</td><td>4</td><td>0.0~1.0</td><td>0.06~0.69</td></tr><tr><td>r_SOL_K(.sol)</td><td>Saturated water conductivity coefficient</td><td></td><td>10</td><td>0.0~25</td><td>0.88~2.66</td></tr><tr><td>r_SOL_Z(.sol)</td><td>Soil depth</td><td></td><td>8</td><td>-0.5~0.5</td><td>0.10~0.31</td></tr><tr><td>v_CH_K2.rte</td><td>River effective water conductivity coefficient</td><td></td><td>2</td><td>-0.01~150</td><td>47.97~144.52</td></tr><tr><td>v_RCHRG_DP.gw</td><td>Water table permeation</td><td></td><td>3</td><td>0.0~1.0</td><td>0.26~0.79</td></tr><tr><td>v_EPCO.hru</td><td>Material transpiration coefficient</td><td></td><td>5</td><td>0.0~1.0</td><td>0.37~1.12</td></tr><tr><td>v_BIOMIX.mgt</td><td>Biological mix efficiency coefficient</td><td></td><td>7</td><td>0.0~1.0</td><td>0.03~0.68</td></tr><tr><td>v_SPCON.bsn</td><td>Sediment movement linear compensation coefficient</td><td></td><td>1</td><td>0.0001~0.01</td><td>0.0057~0.0097</td></tr><tr><td>v_CH_COV2.rte</td><td>River coverage coefficient</td><td>Sediment load</td><td>2</td><td>0.0~1.0</td><td>0.63~1.00</td></tr><tr><td>v_SPEXP.bsn</td><td>Sediment movement index coefficient</td><td></td><td>4</td><td>0.0~1.5</td><td>1.20~1.42</td></tr><tr><td>v_CH_ERODMO().rte</td><td>River erosion coefficient</td><td></td><td>3</td><td>0.0~1.0</td><td>0.75~1.00</td></tr><tr><td>v_USLE_P.mgt</td><td>Water and land conservation measure factors</td><td></td><td>5</td><td>0.0~1.0</td><td>0.56~1.00</td></tr></tbody></table>

Notes: Beg. Range shows beginning range. .gw shows groundwater hydrology; .rte shows river hydrology; .hru shows HRU conventional data entry; .mgt shows HRU management; .sol shows soil data entry; .bsn shows conventional basin parameters; v shows parameter changes by the specified value; r shows parameter changes by the original value, 1+ specified value.

According to the results of sudden change analysis, 1969–1987 and 1986–2015 were considered as the calibration and validation periods. The preheating period for both the calibration and validation periods was 2 years, and the runoff was calibrated first, followed by sediment load. The $R^2$ and Ens of the runoff simulation in the calibration period were 0.92 and 0.91, respectively, while the corresponding values of the validation period were 0.90 and 0.87 with high simulation accuracy. The $R^2$ and Ens for sediment load simulation in the calibration period were 0.60 and 0.55, respectively, the corresponding values of the validation period were 0.70 and 0.69, which simulation effect are
satisfactory. The comparison of the simulated and observed values of the runoff and sediment load is shown in Figure 5a,b.

<figure><img src="image_80.png" /><figcaption>**Figure 5.** (a) runoff and (b) sediment load calibration and validation periods simulation results at Xihe River Basin during the period 1971–2015.</figcaption></figure>

### 3.3. Xihe River Basin Soil Erosion Conditions

Figure 6 shows the spatial distribution of annual average soil erosion modulus in the Xihe River Basin. The study area lies in the northeast black soil region. According to the Classification Standard for Soil Erosion SL 190-2007 [44], the allowable soil loss in the region is 200 t/(km²·a). Soil erosion modulus under 200 t/(km²·a) is classified as micro erosion, and soil erosion modulus between 200 and 2500 t/(km²·a) is classified as mild erosion. As can be seen in Figure 6, only sub-basin regions 1, 2, 5, 7, 9, and 15 are within the allowable range of soil erosion and have an area of 279.063 km², occupying 25.23% of the study area. Other areas of sub-basin experience more erosion, representing an area of 826.937 km², occupying 74.77% of the study area. Furthermore, the most severe soil erosion is in sub-basin region 19 in the Xihe River Basin, which has an average soil erosion modulus of 595.4 t/(km²·a). Overall, the upper part of the basin has a higher erosion rate than the lower part.
<figure><img src="image_81.png" /><figcaption>**Figure 6.** Soil erosion modulus spatial distribution in the Xihe River Basin.</figcaption></figure>

Table 3 shows the annual average soil erosion volume and modulus by types of land use. Cultivated land exhibited the highest levels of erosion, 3.6561 × 10⁵ t, followed by forest land, and lastly grassland. The highest levels of soil erosion modulus are also exhibited by cultivated land, 1856.56 t/(km²·a), followed by grassland, then urban land, unused land and lastly forest land.

<table><caption>**Table 3.** Xihe River Basin soil erosion by types land use.</caption><thead><tr><th>Land Use Type</th><th>Area/km²</th><th>Soil Erosion Volume/10⁴ t</th><th>Average Soil Erosion Modulus (t·km<sup>-2</sup>·a<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Cultivated land</td><td>194.5</td><td>36.11</td><td>1856.56</td></tr><tr><td>Forest land</td><td>865.06</td><td>1.11</td><td>12.83</td></tr><tr><td>Grassland</td><td>2.78</td><td>0.0099</td><td>35.61</td></tr><tr><td>Unused land</td><td>7.3</td><td>0.0158</td><td>21.64</td></tr><tr><td>Urban land</td><td>22.04</td><td>0.0633</td><td>28.72</td></tr><tr><td>Water</td><td>14.32</td><td>-</td><td>-</td></tr></tbody></table>

### 3.4. Effects of Climate Change and Human Activities on Runoff and Sediment Load

Table 4 shows the contributions of climate change and human activities on runoff and sediment load. Compared to the baseline period, annual average runoff decreased by 29.95 m³/s in the impact period. Climate change contributed an increase of 7.95 m³/s accounting for -26.54% of the total change, and human activities contributed a decrease of 37.9 m³/s accounting for 126.54% of the total change. Annual average sediment load decreased by 4.1 × 10⁵ t in the impact period, to which climate change and human activities contributed an increase of 0.36 × 10⁵ t (-8.8% of the total) and a decrease of 4.45 × 10⁵ t (108.8% of the total), respectively. The overall decrease in runoff and sediment load can be attributed to human factors, which exceeded the effect of climate change and may be considered a main factor impacting the hydrological cycle of the Xihe River Basin.
<table><caption>**Table 4.** Relative contributions of climate change and human activities to runoff and sediment load.</caption><thead><tr><th></th><th>Period</th><th>Q<sub>v</sub></th><th>Q<sub>s</sub></th><th>ΔQ<sub>t</sub></th><th>ΔQ<sub>c</sub></th><th>ΔQ<sub>h</sub></th><th>P<sub>c</sub>(%)</th><th>P<sub>h</sub>(%)</th></tr></thead><tbody><tr><td rowspan="2">Runoff (m<sup>3</sup>/s)</td><td>Calibration</td><td>136.94</td><td>119.62</td><td rowspan="2">-29.95</td><td rowspan="2">7.95</td><td rowspan="2">-37.9</td><td rowspan="2">-26.54</td><td rowspan="2">126.54</td></tr><tr><td>Validation</td><td>106.99</td><td>127.57</td></tr><tr><td rowspan="2">Sediment load (10<sup>5</sup> t)</td><td>Calibration</td><td>7.18</td><td>3.52</td><td rowspan="2">-4.1</td><td rowspan="2">0.36</td><td rowspan="2">-4.45</td><td rowspan="2">-8.8</td><td rowspan="2">108.8</td></tr><tr><td>Validation</td><td>3.08</td><td>3.88</td></tr></tbody></table>

Notes: Q_v shows the calibration and validation period observed values, Q_s shows simulated values.

## 4. Discussion

Accuracy of runoff simulation of the SWAT model is better than sediment load simulation. The same phenomenon has also been found in other studies [45,46]. In the SWAT model, runoff and peak rate are the main input factors in estimating the sediment discharge equation, but when SWAT model estimates runoff, precipitation intensity, frequency and duration are neglected [47], which would lead to an underestimation of the runoff, finally, the sediment load simulation results in large differences.

The average soil erosion modulus of the Xihe River Basin is 334.06 t/(km²·a), which is similar to the counties and cities in Chaoyang city in the Liaoning Province, the results of which were studied by Yi Kai et al. [48] using the RUSLE (Revised Universal Soil Loss Equation) model. However, our result is lower than that study. The difference is not only because of the difference between climatic conditions and human activities in these two regions, but is also related to the different methods applied. In SWAT, soil erosion is calculated by the USLE (Universal Soil Loss Equation) model. Firstly, the USLE model only considers the changes in rainfall erosivity factor, and does not consider other changes in related to soil erosion, such as the vegetation coverage, the RUSLE model is more complete. Secondly, the RUSLE model obtained parameters through observing or taking similar empirical values, such as soil and water conservation measures factor and vegetation cover factor, without calibrating and validating. The SWAT model selected parameters by sensitivity analysis and calibrated these parameters with observed runoff and sediment load. Then soil erosion was calculated using the best simulated value.

The method of quantifying the contribution rate of climate change and human activities to runoff and sediment load in changes is simple, the disadvantage of it is that it cannot reflect the changes of runoff and sediment load caused by the changes of soil type, topography and geomorphology. During the study period, the contribution rates of human activities to runoff and sediment load reduction were 126.54% and 108.8%, respectively, which were far greater than the impact of climate change. This result is similar to the Jinghe River Basin, an arid inland basin in northwest China [49]. Climate change usually impacts runoff variation periodically and lastingly, whereas human activities impacts runoff variation suddenly and directionally [50]. The results could be useful for understanding the variation and driving factors for runoff and sediment load reduction in the basin.

## 5. Conclusions

This study analyzed the inter-annual validation trends and sudden changes in hydrological and meteorological variables in the Xihe River Basin from 1969–2015. The SWAT model was used to simulate the hydrological cycle and the spatial distribution characteristics of soil erosion. The study quantified the contributions of climate change and human activities to runoff and sediment load in changes. Which are summarized as follows:

(1) During the period 1969–2015, annual rainfall and temperature exhibited an increasing trend. Whereas annual runoff and sediment exhibited a decreasing trend. A significant trend was observed in the annual temperature and sediment load. In 1988, the annual temperature exhibited a sudden increase, but sediment load showed a sudden decline.

(2) The evaluating indicators showed that the R² and Ens of the SWAT model rose to an acceptable level after calibration, which indicated that the hydrological cycle of the Xihe River Basin could
be accurately simulated using the SWAT model. The simulation effect of runoff was better than that of the sediment load.

(3) Soil erosion in the upper part of the Xihe River Basin was more serious than that of the lower part, and the soil erosion modulus for cultivated land was highest. The contribution of human activities to runoff and sediment load changes was greater than that of climate change. Therefore, human activities is the primary factor that affects the hydrological cycle in the Xihe River Basin. In addition, it should be noted that the upstream should be paid more attention. Finally, the specific human activities that affect runoff and sediment load, as well as their modes of action, require further study.

**Author Contributions:** Conceptualization, S.G. and L.L.; Methodology, S.G. and L.L.; Software, S.G. and L.L.; Validation, L.L. and Z.Z.; Formal Analysis, S.G.; Investigation, S.G.; Resources, L.L. and Z.Z.; Data Curation, S.G.; Writing—Original Draft Preparation, S.G.; Writing—Review & Editing, L.L.; Visualization, S.G.; Supervision, L.L. and Z.Z.; Project Administration, L.L.; Funding Acquisition, L.L.

**Funding:** This research was funded by National Natural Science Foundation of China grant number 41701208.

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

## References

1. Pimentel, D. Soil Erosion: A Food and Environmental Threat. *Environ. Dev. Sustain.* **2006**, *8*, 119–137.
   [CrossRef]

2. Munro, R.N.; Deckers, J.; Haile, M.; Grove, A.T.; Poesen, J.; Nyssen, J. Soil landscapes, land cover change and erosion features of the Central Plateau region of Tigrai, Ethiopia: Photo-monitoring with an interval of 30 years. *Catena* **2008**, *75*, 55–64. [CrossRef]

3. Hu, S.; Cao, M.M.; Liu, Q.; Zhang, T.Q.; Qiu, H.J.; Liu, W.; Song, J.X. Comparative study on the soil conservation function of InVEST model under different perspectives. *Geogr. Res.* **2014**, *33*, 2393–2406.
   [CrossRef]

4. Qiao, Y.; Qiao, Y. Fast soil erosion investigation and dynamic analysis in the loess plateau of China by using information composite technique. *Adv. Space Res.* **2002**, *29*, 85–88.

5. Dotterweich, M. The history of human-induced soil erosion: Geomorphic legacies, early descriptions and research, and the development of soil conservation—A global synopsis. *Geomorphology* **2013**, *201*, 1–34.
   [CrossRef]

6. Hamel, P.; Chaplin-Kramer, R.; Sim, S.; Mueller, C. A new approach to modeling the sediment retention service (InVEST 3.0): Case study of the Cape Fear catchment, North Carolina, USA. *Sci. Total Environ.* **2015**, *524–525*, 166–177. [CrossRef] [PubMed]

7. Khanal, S.; Parajuli, P.B. Evaluating the Impacts of Forest Clear Cutting on Water and Sediment Yields Using SWAT in Mississippi. *J. Water Resour. Prot.* **2013**, *5*, 474–483. [CrossRef]

8. Li, T.; Liu, K.; Ma, L.Y.; Bao, Y.B.; Wu, L. Evaluation on Soil Erosion Effects Driven by Land Use Changes over Danjiang River Basin of Qinling Mountain. *J. Nat. Resour.* **2012**, *23*, 2249–2256. [CrossRef]

9. Nicu, I.C. Natural risk assessment and mitigation of cultural heritage sites in North-eastern Romania (Valea Oii river basin). *Area* **2018**. [CrossRef]

10. Chen, S.X.; Yang, X.H.; Xiao, L.L.; Cai, H.Y. Study of Soil Erosion in the Southern Hillside Area of China Based on RUSLE Model. *Resour. Sci.* **2014**, *36*, 1288–1297.

11. Golosov, V.; Koiter, A.; Ivanov, M.; Maltsev, K.; Gusarov, A.; Sharifullin, A.; Radchenko, I. Assessment of soil erosion rate trends in two agricultural regions of European Russia for the last 60 years. *J. Soils Sediments* **2018**. [CrossRef]

12. Panagos, P.; Borrelli, P.; Poesen, J.; Ballabio, C.; Lugato, E.; Meusburger, K.; Montanarella, L.; Alewell, C. The new assessment of soil loss by water erosion in Europe. *Environ. Sci. Policy* **2015**, *54*, 438–447. [CrossRef]

13. Labrière, N.; Locatelli, B.; Laumonier, Y.; Freycon, V.; Bernoux, M. Soil erosion in the humid tropics: A systematic quantitative review. *Agric. Ecosyst. Environ.* **2015**, *203*, 127–139. [CrossRef]

14. Keesstra, S.; Pereira, P.; Novara, A.; Brevik, E.C.; Molina, C.A.; Alcántara, L.P.; Jordán, A.; Cerdà, A. Effects of soil management techniques on soil water erosion in apricot orchards. *Sci. Total Environ.* **2016**, *551–552*, 357–366. [CrossRef] [PubMed]
15. Arnold, J.G.; Allen, P.M.; Bernhardt, G. A comprehensive surface-groundwater flow model. *J. Hydrol.* **1993**, *142*, 47–69. [CrossRef]

16. Arnold, J.G.; Srinivasan, R.; Muttiah, R.S.; Williams, J.R. Large area hydrologic modeling and assessment part I: Model development. *J. Am. Water Resour. Assoc.* **1998**, *34*, 73–89. [CrossRef]

17. Hajigholizadeh, M.; Melesse, A.; Fuentes, H. Erosion and Sediment Transport Modelling in Shallow Waters: A Review on Approaches, Models and Applications. *Int. J. Environ. Res. Public Health* **2018**, *15*, 518. [CrossRef] [PubMed]

18. Cornelissen, T.; Diekkrüger, B.; Giertz, S. A comparison of hydrological models for assessing the impact of land use and climate change on discharge in a tropical catchment. *J. Hydrol.* **2013**, *498*, 221–236. [CrossRef]

19. Abdelwahab, O.M.M.; Ricci, G.F.; De Girolamo, A.M.; Gentile, F. Modelling soil erosion in a Mediterranean watershed: Comparison between SWAT and AnnAGNPS models. *Environ. Res.* **2018**, *166*, 363–376. [CrossRef] [PubMed]

20. Aghakhani Afshar, A.; Hassanzadeh, Y. Determination of Monthly Hydrological Erosion Severity and Runoff in Torogh Dam Watershed Basin Using SWAT and WEPP Models. *Iran. J. Sci. Technol. Trans. Civ. Eng.* **2017**, *41*, 221–228. [CrossRef]

21. Hasan, Z.A.; Hamidon, N.; Yusof, M.S.; Ghani, A.A. Flow and sediment yield simulations for Bukit Merah Reservoir catchment, Malaysia: A case study. *Water Sci. Technol.* **2012**, *66*, 2170–2176. [CrossRef] [PubMed]

22. Xiao, J.C.; Luo, D.G.; Wang, Z.Z. Soil Erosion Simulation in Fuhe Basin Based on SWAT Model. *Res. Soil Water Conserv.* **2013**, *20*, 18–24.

23. Zuo, D.P.; Xu, Z.X.; Yao, W.Y.; Jin, S.Y.; Xiao, P.Q.; Ran, D.C. Assessing the effects of changes in land use and climate on runoff and sediment yields from a watershed in the Loess Plateau of China. *Sci. Total Environ.* **2016**, *544*, 238–250. [CrossRef] [PubMed]

24. Wu, J.W.; Miao, C.Y.; Yang, T.T.; Duan, Q.Y.; Zhang, X.M. Modeling streamflow and sediment responses to climate change and human activities in the Yanhe River, China. *Hydrol. Res.* **2018**, *49*, 150–162. [CrossRef]

25. Alighalehbabakhani, F.; Miller, C.J.; Selegean, J.P.; Barkach, J.; Sadatiyan Abkenar, S.M.; Dahl, T.; Baskaran, M. Estimates of sediment trapping rates for two reservoirs in the Lake Erie watershed: Past and present scenarios. *J. Hydrol.* **2017**, *544*, 147–155. [CrossRef]

26. Yen, H.; Lu, S.L.; Feng, Q.Y.; Wang, R.Y.; Gao, J.G.; Brady, D.; Sharifi, A.; Ahn, J.; Chen, S.T.; Jeong, J.; et al. Assessment of Optional Sediment Transport Functions via the Complex Watershed Simulation Model SWAT. *Water* **2017**, *9*, 76. [CrossRef]

27. Dong, L.H.; Xiong, L.H.; Yu, K.X.; Li, S. Research advances in effects of climate change and human activities on hydrology. *Adv. Water Sci.* **2012**, *23*, 278–285.

28. Yang, C.L.; Ma, X.P.; Hou, W.; Li, F.Y.; Liu, Q.; Li, Y.; Cheng, Z.H.; Kong, W.J. Soil Physical and Chemical Properties of Riparian Zone Along Xi River. *Sci. Technol. Rev.* **2012**, *30*, 61–66. [CrossRef]

29. Pettitt, A.N. A non-parametric approach to the change point problem. *J. R. Stat. Soc.* **1979**, *28*, 126–135. [CrossRef]

30. Zuo, D.P.; Xu, Z.X.; Wu, W.; Zhao, J.; Zhao, F.F. Identification of Streamflow Response to Climate Change and Human Activities in the Wei River Basin, China. *Water Resour. Manag.* **2014**, *28*, 833–851. [CrossRef]

31. Liu, M.F.; Gao, Y.C.; Gan, G.J. Long-Term Trends in Annual Runoff and the Impact of Meteorological Factors in the Baiyangdian Watershed. *Resour. Sci.* **2011**, *33*, 1438–1445.

32. Wang, J.; Gao, Y.; Wang, S. Assessing the response of runoff to climate change and human activities for a typical basin in the Northern Taihang Mountain, China. *J. Earth Syst. Sci.* **2018**, *127*, 1–15. [CrossRef]

33. Li, Y.Y.; Chang, J.X.; Wang, Y.M.; Jin, W.T.; Guo, A.J. Spatiotemporal Impacts of Climate, Land Cover Change and Direct Human Activities on Runoff Variations in the Wei River Basin, China. *Water* **2016**, *8*, 220. [CrossRef]

34. Li, P.F.; Mu, X.M.; Holden, J.; Wu, Y.P.; Irvine, B.; Wang, F.; Gao, P.; Zhao, G.J.; Sun, W.Y. Comparison of soil erosion models used to study the Chinese Loess Plateau. *Earth Sci. Rev.* **2017**, *170*, 17–30. [CrossRef]

35. Yin, J.; He, F.; Xiong, Y.J.; Qiu, G.Y. Effects of land use/land cover and climate changes on surface runoff in a semi-humid and semi-arid transition zone in northwest China. *Hydrol. Earth Syst. Sci.* **2017**, *21*, 183–196. [CrossRef]

36. USDA Soil Conservation Service (SCS). *National Engineering Handbook Section 4 Hydrology*; USDA: Washington, DC, USA, 1972.
37. Monteith, J.L. Evaporation and the environment. In *The State and Movement of Water in Living Organism*; Cambridge University Press: Swansea, UK, 1965.

38. William, J.R. Chapter 25: The EPIC model. In *Computer Models of Watershed Hydrology*; Singh, V.P., Ed.; Water Resources Publications, LLC: Highlands Ranch, CO, USA, 1995; pp. 909–1000.

39. Neitsch, S.L.; Arnold, J.G.; Kiniry, J.R.; Williams, J.R. *Soil and Water Assessment Tool Theoretical Documentation Version 2009*; Technical Report No. 406; Texas Water Resources Institute: College Station, TX, USA, 2011.

40. Che, Z.H. The empirical formulas and curves of soil permeability system number are discussed. *Water Resour. Hydropower Northeast China* **1995**, *19*, 17–19. [CrossRef]

41. Abbaspour, K.C.; Yang, J.; Maximov, I.; Siber, R.; Bogner, K.; Mieleitner, J.; Zobrist, J.; Srinivasan, R. Modelling hydrology and water quality in the pre-alpine/alpine Thur watershed using SWAT. *J. Hydrol.* **2007**, *333*, 413–430. [CrossRef]

42. Nash, J.E.; Sutcliffe, J.V. River flow forecasting through conceptual models part I-a discussion of principles. *J. Hydrol.* **1970**, *10*, 282–290. [CrossRef]

43. Moriasi, D.N.; Arnold, J.G.; Van Liew, M.W.; Bingner, R.L.; Harmel, R.D.; Veith, T.L. Model Evaluation Guidelines for Systematic Quantification of Accuracy in Watershed Simulations. *Am. Soc. Agric. Biol. Eng.* **2007**, *50*, 885–900.

44. SL 190-2007. *Classification Standard for Soil Erosion*; WaterPower Press: Beijing, China, 1977.

45. Zeng, Y.; Wei, L. Runoff and Sediment Simulation in Purple Hilly Area Based on SWAT Model. *J. Geo-Inf. Sci.* **2013**, *15*, 401–407. [CrossRef]

46. Pang, J.P.; Liu, C.M.; Xu, Z.X. Streamflow and Soil Erosion Simulation Based on SWAT Model. *Res. Soil Water Conserv.* **2007**, *14*, 88–93.

47. Assouline, S.; Ben-Hur, M. Effects of rainfall intensity and slope gradient on the dynamics of interrill erosion during soil surface sealing. *Catena* **2006**, *66*, 211–220. [CrossRef]

48. Yi, K.; Wang, S.Y.; Wang, X.; Yao, H.L. The Characteristics of Spatial-temporal Differentiation of Soil Erosion Based on RUSLE Model: A Case Study of Chaoyang City. *Sci. Geogr. Sin.* **2015**, *35*, 365–372.

49. Dong, W.; Cui, B.; Liu, Z.; Zhang, K. Relative effects of human activities and climate change on the river runoff in an arid basin in northwest China. *Hydrol. Process.* **2013**, *28*, 4854–4864. [CrossRef]

50. Miao, C.; Yang, Y.; Liu, B.; Gao, Y.; Li, S. Streamflow changes and its influencing factors in the mainstream of the Songhua River basin, Northeast China over the past 50 years. *Environ. Earth Sci.* **2011**, *63*, 489–499. [CrossRef]

<figure><img src="image_82.png" /></figure>

© 2018 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/).
Article

# The Role of Attenuation and Land Management in Small Catchments to Remove Sediment and Phosphorus: A Modelling Study of Mitigation Options and Impacts

Russell Adams<sup>1,*</sup>, Paul Quinn<sup>1</sup>, Nick Barber<sup>2</sup> and Sim Reaney<sup>2●</sup>

Received: 22 August 2018; Accepted: 10 September 2018; Published: 12 September 2018

<figure><img src="image_83.png" /></figure>

check for
updates

**Abstract:** It is well known that soil, hillslopes, and watercourses in small catchments possess a degree of natural attenuation that affects both the shape of the outlet hydrograph and the transport of nutrients and sediments. The widespread adoption of Natural Based Solutions (NBS) practices in the headwaters of these catchments is expected to add additional attenuation primarily through increasing the amount of new storage available to accommodate flood flows. The actual type of NBS features used to add storage could include swales, ditches, and small ponds (acting as sediment traps). Here, recent data collected from monitored features (from the Demonstration Test Catchments project in the Newby Beck catchment (Eden) in northwest England) were used to provide first estimates of the percentages of the suspended sediment (SS) and total phosphorus (TP) loads that could be trapped by additional features. The Catchment Runoff Attenuation Flux Tool (CRAFT) was then used to model this catchment (Newby Beck) to investigate whether adding additional attenuation, along with the ability to trap and retain SS (and attached P), will have any effect on the flood peak and associated peak concentrations of SS and TP. The modelling tested the hypothesis that increasing the amount of new storage (thus adding attenuation capacity) in the catchment will have a beneficial effect. The model results implied that a small decrease of the order of 5–10% in the peak concentrations of SS and TP was observable after adding 2000 m³ to 8000 m³ of additional storage to the catchment.

**Keywords:** runoff; suspended sediment; phosphorus; water quality modelling; mitigation measures; flooding

## 1. Introduction

It is becoming widely accepted that Nature Based Solutions (NBS—Nature Based Green Instructure Solutions) [1] and “Natural Flood Management” (NFM) (defined by the United Kingdom Environment Agency as part of Working with Natural Processes) [2] can have a positive impact in terms of reducing flooding, most observably by lowering the peak discharge of the flood hydrograph to enable this [3–7]. The construction of different types of “soft engineered” measures (or features) in headwater catchments has become an established part of this strategy [2,8]. Previous studies of the performance of features have concentrated primarily on the attenuation capabilities of these features in terms of reducing flooding, e.g., Belford Burn [9] and Pickering Beck [4,10] (in the U.K.). Moreover, the improvement in water quality (quantifiable by a reduction in concentrations and/or loads of nutrients and sediments) brought about by the construction of features in rural catchments has been studied in the U.K. [11–13], Irish Republic [14], and in New Zealand [15], but this issue has generally received less attention than the mitigating benefits of NFM in terms of reducing storm events.

<footnote n="1">¹ School of Engineering, Newcastle University, Newcastle Upon Tyne NE1 7RU, UK; p.f.quinn@ncl.ac.uk</footnote>

<footnote n="2">² Department of Geography, Durham University, Durham DH1 3LE, UK; nick_barbie@hotmail.com (N.B.); sim.reaney@durham.ac.uk (S.R.)</footnote>

<footnote>* Correspondence: russell.adams@ncl.ac.uk; Tel.: +44-191-208-5773</footnote>
Avery [16] coined the term “rural sustainable drainage systems” (RSuDS) to reflect the trend to construct wetlands and other types of features in rural catchments, and since the focus of this study is on the use of mitigation features to improve water quality rather than NFM per se, the term RSuDS henceforth is not used. The term “runoff attenuation feature” (RAF) has also been used in the literature [9,11]; essentially, an RSuD is a type of RAF that adds attenuation to a ditch or channel in which it is constructed [8]. However, the design purpose of RSuDs by default leads to sediment trapping and P stripping and associated “buffering” of other chemicals and microbial pollutants. Therefore, their ability to store and hence attenuate larger flood flows is lower. Often RAFs are designed to target flood flows primarily and thus the ideal solution is to both trap sediment and attenuate larger flow. This study explores the key role of adding attenuation to a catchment to target both water quantity and quality issues.

Environmentalists are interested in spatial patterns because they are essential in the scaling-up from localised measurements to larger spatial scales in order to provide assessments of mitigation impacts on pollution at the catchment, regional, or national scale for policy purposes [17]. In terms of addressing the impacts of mitigation features, few studies have attempted to assess the effect of mitigation features at the catchment scale [18]. These impacts have often been monitored by local water sampling and the measurement of runoff entering and leaving the features. For example, water quality, sediments, and nutrients have commonly been measured by automatic water quality samplers in order to collect data on concentrations and (if flow was measured) loads [12]. However, an important research question still remains as to whether impacts measured at the experimental scale are observable downstream, where monitoring points are often located (e.g., U.K. Environment Agency weekly sampling sites). Longer term monitoring programs (e.g., the Irish Agricultural Catchments (IAC) programme [19] and the Demonstration Test Catchments (DTC) programme [20,21]) are required if the larger scale impacts of these features are to be detected, if this is indeed possible with existing monitoring networks, and to address climatic issues (e.g., floods and droughts). Thus, the evidence that shows the rewards of adding attenuation capacity to a catchment are needed by end users, which is one goal of the DTC programme [21].

It is known from measuring water levels to observe runoff events before and after their construction that mitigation features have the potential to add attenuation to ditches and/or headwater streams [4,5]. These features can either: (i) divert water from channels via draw-off structures to separate storage areas or disconnected channels (classified as “off-line”), or (ii) temporarily detain runoff using “in-line” interventions located within low-order streams and ditches [11]. In-line features involve direct intervention in the channel or ditch itself such as the creation of artificial barriers, such as large woody debris (LWD) and engineered log jams (ELJs), and can also be applied in combination with off-line features [22]. Off-line features include riparian buffer strips, swales (vegetated channels), and ponds [5,13,23].

In terms of sediments and nutrients (principally total phosphorus (TP) and suspended sediment (SS)), the evidence from the case studies [11,24] is that these mitigation features can trap significant quantities of particulates with attached, insoluble forms of P. However, Barber [11] stated that few studies [13,23,25] had addressed either the effectiveness of mitigation features at the catchment scale, and suggested that, in order to address the requirement for urgent action with respect to meeting water quality targets in U.K. agricultural catchments, further research was required.

The natural attenuation of SS and bound nutrients, including forms of nitrogen (N) and P in riparian channel systems, is less-widely studied; however, one U.S. study [26] did highlight an important ecosystem function where the channels retained N and P exported from row crop fields during baseflow conditions, thus preventing higher exports into the estuary downstream in their catchment located in South Carolina, USA.

The primary aim of the study is to address the impacts of land management by altering hydrological flow paths and the overall catchment attenuation capacity on flow rates and nutrient losses. The modelling study described below demonstrates whether the impact of adding mitigation
features, i.e., RAFs, at the headwater scale can be observable further downstream at a larger
measurement scale. This modelling allows an estimate of how much attenuation can be achieved and
the corresponding loss of productive land that may be required. A secondary aim is to explore whether
the chosen model can simulate improvements to land management in a catchment that are designed to
reduce losses of sediment and P in surface runoff. A further important research question poses: "Are
there any significant differences in the performance of different types of features?" This can be limited
by the available data on the performance of mitigation features [11,24].

A modelling case study to pursue the above aims used data collected from a catchment-based field programme in northern England [11,24] to parameterise the Catchment Runoff Attenuation Flux Tool (CRAFT) [27,28]. Scenarios were considered where: (i) RAFs were simulated by the model (by adding attenuation storage and trapping sediment and associated particulate forms of P) in order to see if the first aim can be achieved by investigating the impact on both the runoff hydrographs and the time series of P and SS concentrations modelled at the catchment outlet and also at the outlet of a mitigated sub-catchment, and (ii) land management options are applied instead to the same sub-catchment in order to reduce surface runoff from fields in the catchment and associated losses of sediment and particulate forms of P in surface runoff.

**2. Methods**

The methodology underpinning this study can be summarized by the following steps:

1. Develop the CRAFT model (at Newcastle University, Newcastle Upon Tyne, UK) to simulate nutrient and sediment fluxes at the catchment scale and add the capability to attenuate these fluxes in the surface runoff [27,28].

2. Calibrate the CRAFT to the existing runoff, nutrient, and sediment data collected in October 2011–September 2012 and establish a baseline scenario [28].

3. Evaluate the performance of a series of demonstration RAFs in the Mitigation sub-catchment of the NBC (mass of sediment and P trapped) [24] and use the information to inform the future scenarios.

4. Run simulations of the NBC using CRAFT with additional attenuation to represent three scenarios of land management.

5. Interpret the results of these scenarios in terms of (i) local and catchment scale impacts on runoff, TP, and SS; and (ii) compare the modelled changes in TP and SS yield with the measured reductions from step 3.

2.1. *Description of the CRAFT*

The Catchment Runoff Attenuation Flux Tool (CRAFT) [27,28] was selected for the case study.
In this study, use was made of the model's attenuation store, which attenuates the surface runoff
generated from rainfall excess. The store is connected, therefore, to only one of the three runoff
pathways in CRAFT [27]. The store mimics the physical processes of attenuation using the minimum
information requirement [29] (MIR approach) in this case by utilizing a linear storage–discharge
relationship [28]. The attenuation storage depth was calculated by the model at each timestep so the
total volume of storage in the catchment was obtained by multiplying this by catchment area. It was
necessary to specify the maximum storage depth in the model. If the depth in the store exceeded this
value, then the excess runoff was added to the outflow from the store. The rate of drainage of the
attenuation store was controlled by the K<sub>LAG</sub> parameter; this is equivalent to the reciprocal of lag time
1/k in the storage discharge relationship for a linear routing model (S = kQ) [30].

In terms of definitions of various P species, the CRAFT can output the following forms of P:
particulate (insoluble, unreactive) P (PUP), particulate (insoluble, reactive) P (PRP), soluble reactive P
(SRP) in groundwater, and fast subsurface flow. Therefore, in terms of the modelled concentrations
and fluxes the P forms are combined so that: TRP = PRP + SRP and TP = TRP + PUP.
## 2.2. *Description of Case Study*

The catchment modelled in the case study (Newby Beck) is one of the three instrumented during the EdenDTC project, located in the River Eden catchment in northwest England [20,31–33]. It contains both control and Mitigation sub-catchments of similar sizes (approximately 2 km²) out of a total area of 12.5 km². Capitals will be used to denote the monitored Mitigation sub-catchment to distinguish it from the modelled one. The main source of sediment in the wider Eden catchment was found to be bank erosion (exacerbated by livestock accessing watercourses and causing poaching), and also caused by streams undercutting the banks leading to instability [34]. Reference [34] found that good vegetation cover ensured that erosion by overland flow entrainment was likely to be a minor source of suspended sediment in the channels; however, field observations made in the EdenDTC project [32] have identified field sources of suspended sediment that become active during runoff events (e.g., tracks and ditches that become transport pathways).

In terms of the potential for future construction of RAFs (primarily as RSuDS), the areas that need to be treated within a sub-catchment can by determined by either: (i) using the current guidelines in terms of what proportion of the catchment area should be treated (this area varies typically from at least 1% (as advocated in New Zealand [15]) for riparian wetland coverage in dairy catchments), and up to 5% (advocated by Quinn [5]); or (ii) calculating the volume of attenuation storage that is required in the catchment to achieve targeted water quality improvements. In Europe, the Water Framework Directive [35] has been responsible for setting water quality targets for rivers and lakes and such targets can be used to examine what improvements are required.

Figure 1 shows a map of the Eden catchment showing the Newby Beck catchment (NBC) and the location of the monitoring sites for rainfall and flow in the catchment (these have been monitored as part of the EdenDTC project). Contours at 20 m intervals depict the topography. The dashed line shows the outline of the Mitigation sub-catchment.

<figure><img src="image_84.png" /><figcaption>**Figure 1.** Map of the Eden catchment showing the NBC sub-catchment studied as part of the EdenDTC project, and the location of the monitoring sites for rainfall and flow in the sub-catchments. Contours at 20 m intervals depict the topography.</figcaption></figure>
The monitoring data collected at the NBC outlet (i.e., the entire catchment) at Morland have been summarized previously [20,28,31,32]; the dataset comprises Q, turbidity plus TRP, and TP concentrations measured every 15–30 min but for modelling purposes these values have been converted to hourly values. Observed totals (runoff, rainfall depth and nutrient and sediment loads) are shown here in Table 1 for the time period of interest in this study (April–September 2012). Observed TP concentrations and calculated TP yields were based on samples taken by the bankside equipment and included a small component of dissolved unreactive or organic P (DUP), which was not modelled and assumed to be negligible. Predicted TP export (yields) data are also available from an export coefficient-based modelling study [36], and for the geoclimatic region including the NBC this predicted baseline TP export (year 2000) to be 1.39 kg·ha⁻¹·year⁻¹ P, of which 52.5% originated from diffuse sources.

At both catchments, turbidity was measured at 15-min intervals using an YSI 6600 V2 multi-parameter sonde (YSI Incorporated, Yellow Springs, OH, USA). A strong relationship (from regression analysis) between turbidity and SS (the latter was measured from grab samples collected by an autosampler at the NBC outlet during storm events) was identified [31,32], enabling turbidity to act as a proxy for suspended sediment. These data were used for the calculation of an “observed” yield from the catchment for modelling purposes.

<table><caption>**Table 1.** Observations at the NBC April–September 2012 from monitoring data.</caption><thead><tr><td>Observation</td><td>Value</td></tr></thead><tbody><tr><td>Catchment Area (km²)</td><td>12.5</td></tr><tr><td>Rainfall (mm)</td><td>686</td></tr><tr><td>Runoff (mm)</td><td>303</td></tr><tr><td>TP Yield (kg·ha⁻¹)</td><td>0.73</td></tr><tr><td>SS Yield (t·km⁻²)</td><td>18.1</td></tr><tr><td>TP Load (kg)</td><td>908</td></tr><tr><td>SS Load (t)</td><td>229</td></tr><tr><td>TP mean Concentration (mg P·L⁻¹)</td><td>0.077</td></tr><tr><td>SS mean Concentration (mg·L⁻¹)</td><td>4.3</td></tr></tbody></table>

Nine mitigation features have been constructed in the Mitigation sub-catchment that target surface runoff from two farms located in the headwater area, with a combination of swales, small ponds, and ditches designed to intercept runoff from farm tracks and fields and divert this into the features [24]. The mass of sediment and nutrients trapped by the features has been calculated from the accumulated sediment and then analyzing the removed sediment for nutrient content. These data have been collected on an annual basis since late 2014 and allowed the loads to be estimated from five of the nine features.

In the NBC surface, runoff pathways (including ditches and drains) represent the major runoff pathway for exporting sediment attached P via sediment transport [31,32]. The best fit for a transfer function model for TP load from the NBC was a single store model with a sole quick flow pathway [33]. Therefore, this flow pathway was targeted by the mitigation features in order to develop a strategy that adds attenuation to the outlet hydrograph and load time series of P. Note that no mitigation features were in place in 2012, so the observed data represent baseline conditions prior to any intervention being made.

## 2.3. Mitigation Modelling Approach

The CRAFT has already been calibrated and validated on the NBC [28]. In this earlier study, several simulations of runoff and P were carried out in order to test different hypotheses of conceptual models for the entire NBC, these differed primarily in whether any attenuation of the surface runoff flow pathway was included in the model. The scenario chosen for use in this study as a baseline was the “lagged” [28] in which the modelled hydrographs have added attenuation representing the natural storage in the catchment during runoff events.
The modelling strategy followed the approach of using simulations where the volume of attenuation storage is selected a priori. Implementing the proposed measures would likely require terrain analysis to identify runoff pathways and Critical Source Areas (CSAs) in order to select suitable sites for constructing off-line features [29,37–39].

A series of curves, such as those shown in Figure 2 below, can be plotted that relate a representative set of outputs from the modelling (e.g., Qp i.e., peak flow, or the load, or concentration, of TP or SS) to the degree of attenuation or storage in the catchment, where the left axis represents baseline conditions with no additional features present to add storage capacity. Therefore, the origin represents a minimum degree of attenuation and a small amount of storage from the existing floodplain and riparian areas.

<figure><img src="image_85.png" /><figcaption>**Figure 2.** Sketch of relationship between load or concentration and degree of attenuation and storage (both plotted on the x-axis) in a hypothetical catchment indicating how adding storage can achieve a desired target.</figcaption></figure>

Adding features can then be shown to add storage capacity and thus decrease the value of the output(s) until an optimum target is reached. In terms of a conceptual model of adding attenuation storage to a catchment, it may be expedient to consider the UK Water Framework Directive Targets (WFD-UKTAG) [40] for achieving reductions in TRP concentrations in order to improve the ecological status of the catchments (site-specific values were reported by Ockenden [32]). The two curves shown in Figure 2 could represent different types of RAF with a higher or lower optimum design storage capacity that is required to meet the required optima (indicated by the intersections of the curves and the dashed vertical red lines at points A and B). A similar set of curves relating particulate P (PP) load reduction at the catchment outlet to the managed proportion of the catchment for two large catchments in Austria and Hungary were developed [41] using the PhosFate model. The modelled load reductions were achieved by adopting best management practice (BMP) interventions over part of the catchments. In the U.K., a national scale modelling platform based on the export coefficient method simulated several scenarios of nutrient load reductions [36]. The scenarios relevant to the NBC were based on:
(i) on-farm mitigation measures, and (ii) farming practices modified to comply with WFD targets. In (i), the TP export from sheep was reduced by 50% and P fertilizer loads applied to grass and arable crops were reduced by 50%, and in (ii), TP loads from cattle farming were reduced by 25% in addition to the reductions in (i). These scenarios predicted that TP exports would reduce by 22.1% and 32.9% respectively (based on data from the year 2000).

The following model simulations investigated whether the target of reducing P and SS concentrations could be met by adding storage capacity. It could be argued that this added storage capacity is actually offsetting the loss of natural attenuation storage caused by deforestation and agricultural intensification. What the final target should be is debatable but the premise that more storage capacity gives a better status is required as part of a long-term plan to reach WFD-UKTAG [40] status, for example. Even as a basic estimate a target of 10 mm of new storage capacity over 1 km² would require 10,000 m³. This would require a storage pond of 1 m depth with a surface area of 100 m by 100 m. Implementing NBS [1,2] would suggest that this could be spread throughout the catchment in an RAF network, and hopefully soil improvement and buffer zones would add to the storage capacity as well. It is proposed that ditch management is a primary basis for the first two scenarios, hence for 1 m of storage depth, a minimum of 1 km of ditch with an effective width of 10 m would be needed. To gain this storage capacity, a ditch would need to be widened and barriers to flow constructed; Figure 3 is an actual example of how a traditional narrow "V" shaped ditch can be substantially modified (Netherton Burn catchment [11]).

<figure><img src="image_86.png" /><figcaption>**Figure 3.** An example of a modified ditch with flow barriers acting as both a sediment trap and flood flow storage in Netherton Burn (photo from Barber 2013, credit to N. Barber).</figcaption></figure>

## CRAFT Mitigation Modelling

The three mitigation scenarios were modelled using the CRAFT and discussed below. Figure 4 shows sketches indicating both the flow pathways under each scenario (and the baseline) as well as the design of the RAFs themselves (in Scenarios 2 and 3). For simplicity, the area of the mitigated sub-catchment in the modelling was set to 10% of the NBC area (1.25 km²), which was slightly smaller than the actual Mitigation sub-catchment (so this will be referred to as the "mitigated sub-catchment").
<figure><img src="image_87.png" /><figcaption>**Figure 4.** Schematic describing Runoff Pathways and design sketches of RAFs in Scenarios 1–3, and Runoff Pathways in the Baseline simulation. Note that Scenario 3 has a widened ditch eating into the buffer zone.</figcaption></figure>

Scenario 1 (Figure 4 top pane): Managing the topsoil in agricultural areas (primarily grassland) to improve the soil health and reduce SS and P loads transported in surface runoff (indicated by the thick blue lines). A review of crop and soil management techniques [42] designed to reduce the quantities of P and SS lost to watercourses found that firstly in upland areas, one of the main mechanisms responsible for P and SS losses is erosion by overland flow detaching soil particles with attached P during rainfall events. Second, in arable fields, the following options can reduce P losses: reduced tillage strategies, application of soil conditioners, crop rotation, crop management, and planting either catch or cover crops. In grasslands, reducing stocking densities and preventing highly trafficked areas from developing (e.g., around feed lots or water troughs) in order to avoid poaching of the soil are the main options available to reduce P losses from the surface runoff. This scenario assumed that the entire mitigated sub-catchment was treated in this way by increasing the maximum infiltration capacity of the soil in the model and evaluated the results in terms of changes to P fluxes. It was assumed that some “pollution swapping” could take place between particulate and soluble forms of P, hence the total P loads may not reduce as much as expected before the “experiment” was carried out. This was essentially a 100% intervention to the land use management across the mitigated sub-catchment. The green lines indicate the major flow pathways under this scenario, which was dominated by infiltration and subsurface flow.

Scenario 2 (Figure 4 middle pane): Treating the mitigated sub-catchment’s watercourses with riparian buffer strips by adding 2000 m³ of additional storage to the sub-catchment and trapping 40% of the SS and attached P transported via surface runoff (overland flow). As a first estimate, it was assumed that up to 10% of the sub-catchment area (i.e., 12.5 ha) was set aside for these measures. This was a 10% land use change intervention option idealising the green/ecological corridor with
minimum management other than sediment removal. As the soil was still degraded, overland flow was still dominant during events.

Scenario 3 (Figure 4 lower pane): Adding engineered RAFs to address food flows and to trap and remove 80% of the SS and attached P transported via the surface runoff (overland flow). These features were assumed to be a combination of offline storages and inline ditches that added 8000 m³ of additional storage to the sub-catchment. As a first estimate, it was assumed that 5% of the mitigated sub-catchment area (i.e., 6.25 ha) was set aside and modified for these measures [5]. This was a reduced area, but it required a widening of the ditch to give the increased volume of storage needed to cope with larger storm events. Essentially, the NFM component was being optimised here, and thus a higher level of maintenance may be needed to keep the total infrastructure operating at its optimum.

Scenarios 2 and 3 were modelled as follows: considering that there are two parameters that can be varied in the CRAFT to represent the attenuation and trapping, these are added by the features in the mitigation sub-catchment only, namely (i) $K_{LAG}$ and (ii) a removal or trapping “efficiency” $K_{RE}$, that apply to the modelled pathways of (a) TP (i.e., PRP and PUP) and (b) SS transported by the surface runoff pathway. Their values were pre-selected for Scenarios 2 and 3 (see Table 2) based on expert judgement and evidence from the field experiments conducted in the mitigation sub-catchment.

<table><caption>**Table 2.** Details of the scenarios modelled in the NBC: parameter Values and storages (S). Baseline values are shown for comparison.</caption><thead><tr><th>Scenario</th><th>K<sub>LAG</sub> (h<sup>-1</sup>)</th><th>S (Total) (mm)</th><th>S (Total) (m<sup>3</sup>)</th><th>S (Added) (m<sup>3</sup>)</th><th>K<sub>RE</sub> (-)</th><th>Other Parameters</th></tr></thead><tbody><tr><td>Baseline</td><td>0.75</td><td>4.86</td><td>6075</td><td>0</td><td>0</td><td>No changes</td></tr><tr><td>1</td><td>0.75</td><td>4.86</td><td>6075</td><td>0</td><td>0</td><td>Increased Infiltration Capacity</td></tr><tr><td>2</td><td>0.83</td><td>6.46</td><td>8075</td><td>2000</td><td>0.4</td><td>No changes</td></tr><tr><td>3</td><td>0.93</td><td>11.3</td><td>14,075</td><td>8000</td><td>0.8</td><td>No changes</td></tr></tbody></table>

The parameter values used for the baseline were:

$K_{LAG} = 0.75 h^{-1}$ (representing natural attenuation, from the “lagged” simulation [28]).

$K_{RE} = 0$ (for the baseline with no trapping).

The unmitigated portion of the catchment was modelled with CRAFT using the calibrated baseline parameter set [28]. There was no additional attenuation of flow, sediments, and nutrients, or removal simulated in this portion of the catchment, or in the fast subsurface and deep groundwater flow pathways in any of the scenarios anywhere.

In the CRAFT, nutrients (in this case PUP and PRP fluxes) and SS were modelled slightly differently from runoff in that their removal was permitted from the attenuation store, which represents the ability of the modelled features to remove (i.e., trap) particulates. Therefore, the mass balance Equations (1) to (3) for the stores become (for the components of PRP, PUP and SS transported by the surface runoff (SR) pathway, which is indicated by “SR” in parentheses):

$$PRPL(SR)_{out} = PRPL(SR)_{in} (1 - K_{LAG})(1 - K_{RE}) \quad (1)$$

$$PUPL(SR)_{out} = PUPL(SR)_{in} (1 - K_{LAG})(1 - K_{RE}) \quad (2)$$

$$SSL(SR)_{out} = SSL(SR)_{in} (1 - K_{LAG})(1 - K_{RE}) \quad (3)$$

where $PRPL_{out}$, $PUPL_{out}$, and $SSL_{out}$ were the loads per time step of PRP, PUP, and SS from the mitigated sub-catchments, respectively. CRAFT outputs were the specific discharge (i.e., runoff depth) and specific yield (i.e., load/unit area). Therefore, the model outputs were scaled up by the areas of the mitigated ($A_{mit}$) and unmitigated areas of the catchments ($A_{umit}$) to obtain the total flow and loads
from the entire catchment. SRP loads from the fast subsurface and slow groundwater flow components were added to the loads from the surface runoff pathways in a mass balance to compute a total TP load. Concentrations of TP and SS at the catchment’s or sub-catchment’s outlet were calculated by dividing the total loads by the total flow. The flow from the catchment $Q_{catch}$ was the sum of the flow components from both the mitigated and unmitigated areas, where the suffixes “mit” and “umit” denote these respectively

$$Q_{catch} = Q_{mit} + Q_{umit} \qquad (4)$$

Nutrient (TP) and SS concentrations were calculated in the same way, where TPL, TRPL, and SSL were the loads with suffixes (as above) indicating which part of the catchment (suffix “catch”) the load originated from

$$TPL_{catch} = TPL_{mit} + TPL_{umit} \qquad (5)$$

where

$$TPL = PUPL(SR) + TRPL \qquad (6)$$

$$SSL_{catch} = SSL_{mit} + SSL_{umit} \qquad (7)$$

A simplified, additive mixing model was used to calculate the flows and loads at the outlet, with no additional attenuation added to represent the in-stream reaches (clearly this was a simplification of reality where in-stream routing could further attenuate the outlet hydrograph [43]). This assumption holds for small catchments where the main stem channel length is less than 10 km and the travel time (lag) between hillslope and outlet is of the order of a few hours. The functional unit defined by a CRAFT sub-catchment could equate to a representative elementary area (REA) [44].

Lastly, the maximum volumes of added storage ($V_{add}$) required in the mitigated sub-catchment were supplied by the user for each scenario. The area ($A_{mit}$) of the mitigated sub-catchment is set to 10% of the total catchment area ($A_t$). The model calculated the depth in the attenuation store at each timestep ($D_{add}$), as it works with depths rather than volumes. The attenuation store was empty at the start of the simulation. The maximum storage volume required in the attenuation store $V_{att}$ could therefore be calculated using Equations (8)–(10).

$$A_{mit} = 0.1A_t \qquad (8)$$

$$V_{att} = V_{nat} + V_{add} \qquad (9)$$

where, in general terms, at each timestep t:

$$V(t) = A_{mit} D(t) \qquad (10)$$

Thus, the required model parameter, $D_{max}$ was calculated using Equation (11):

$$D_{max} = V_{att} / A_{mit} \qquad (11)$$

where $D_{max}$ was the maximum depth of water in the attenuation store (m), $V_{nat}$ was the modelled volume of natural storage in the catchment, which comes from the results of calibrated baseline model results (by extracting $D_{max}$ from these).

Since the outflow $Q_{mit}$ from the attenuation store was a function of $K_{LAG}$, the required value of this parameter could be back-calculated from Equation (12). In practice this was obtained by increasing the value of the parameter until the desired value of $D_{max}$ was achieved in the attenuation store.

$$Q_{mit} = V_{att}(1 - K_{LAG}) \qquad (12)$$

The results from these simulations provided the flows and loads of TP and SS from the mitigated sub-catchment (i.e., $Q_{mit}$, TPL$_{mit}$ and SSL$_{mit}$ in Equations (4)–(7)).
### 3. Results

The CRAFT had already been calibrated to baseline conditions during October 2011–September 2012 in the NBC [28]. In this study the parameter values from the “lagged” scenario, where natural attenuation were included, were used for the baseline simulation. The modelled runoff is shown in Figure 5 for comparison against the observed runoff, for the period January 1 to September 30, 2012 at the NBC outlet. Table 1 summarizes the observations recorded during the entire period between April and September 2012 including the three selected events. The SS yield of 18 t km⁻² over a six-month period appears to fit into the middle range of estimates for the larger Eden catchment made at 14 monitoring points, which was 4–73 t km⁻² year [34]. The lowest monitoring point measured yield from a 1373 km² sub-catchment of the River Eden [34].

<figure><img src="image_88.png" /><figcaption>**Figure 5.** Time series plots of observed (“Obs”) and modelled (“Mod”) (baseline) runoff plus rainfall at the NBC outlet from January–September 2012.</figcaption></figure>

#### 3.1. Results from the 2012 Events (Baseline Simulation)

In terms of investigating the temporal scaling of the results, three large events during a wet 5-month period in 2012 were analysed, and these are denoted by the numbered boxes in Figure 5. According to the observations, the annual (2011–2012) exports of TP and SS from the NBC were 1762 kg and 477.5 tonnes respectively [28,32], and from the three events during 2012, the modelled TP load was predicted to be 725 kg and the modelled SS load was 322 tonnes from the baseline simulation, which indicates that these three events alone contributed a significant proportion of the annual P and sediment exports.

#### 3.2. Results from Scenarios 1–3

Table 2 shows the model parameter values used in each of the scenarios along with the values used in the baseline simulation (for comparison). It also shows the total and added storage in the attenuation component of the model in Scenarios 1–3.

The left-hand pane of Figure 6 shows the modelled and observed runoff at the NBC outlet for Event 3. The modelled attenuation storage per unit time and area (expressed as a depth) is also shown by the solid black line, and both modelled runoff and storage are shown for the baseline case only. The right-hand pane of Figure 6 shows the modelled runoff from the mitigated sub-catchment for the (unmitigated) baseline (solid red line) and Scenarios 2 and 3 (dotted and dashed red lines). The black lines show the attenuation storage during the event for the three simulations (same line styles as runoff). Its value is shown on the right-hand axis and Scenario 3 had a far greater effect on the shape of the hydrograph than Scenario 2 due to the much greater additional added storage (8000 m³ vs. 2000 m³, which is the equivalent of storing an additional 6.4 mm vs. 1.6 mm of runoff during events in the attenuation store). The effect was to both flatten and delay the hydrograph peak due to added attenuation. Both runoff and storage in Scenario 1 were identical to the baseline values as no attenuation storage was added to the model in this scenario, so these are not shown for clarity.
<figure><img src="image_89.png" /><figcaption>**Figure 6.** Time series plots of runoff (Q) and storage (S) during Event 3: Left hand (LH) pane shows modelled and observed Q (red) and modelled S (black) for baseline only at the NBC outlet. Right hand (RH) pane shows modelled Q (red) and S (black) at the mitigated sub-catchment outlet with different line types representing the baseline and Scenarios 2 (S2) and 3 (S3).</figcaption></figure>

Figure 7 shows the results from Scenarios 1–3 alongside the baseline (blue) results for the SS (Figure 7a) and TP (Figure 7b) concentrations in event 3. These were extracted at the NBC outlet (left hand panes) and mitigated sub-catchment outlet (right hand panes) and are shown as different coloured lines for each scenario. The effect on TP concentrations during event 3 can be seen in Figure 7b. The two scenarios that added attenuation storage and trapped particulate P (PUP) (Scenario 2 (green) and Scenario 3 (red)) reduced the peak TP concentrations, but Scenario 1 actually increased the TP concentrations during the falling limb event due to an increase in SRP in the fast subsurface flow pathway; however, the maximum TP concentration was reduced by eliminating the surface runoff pathway as a major source of PUP through improved soil management. Hence, the overall outcome for the event was better even though some instantaneous values worsened. Scenarios 1 (purple) and 3 (red) had the greatest effect on TP concentrations during event 3 in terms of reducing them; in the case of Scenarios 2 and 3, this was due to reducing the amount of PUP exported from the mitigated sub-catchment through trapping sediment with attached P.

<figure><img src="image_90.png" /><figcaption>**Figure 7.** (a) LH pane shows the modelled (“Mod”) and observed (“Obs”) SS concentrations at the NBC outlet. RH pane shows the SS concentrations at the mitigated sub-catchment outlet (Scenarios 1–3 shown by colored lines, with baseline in blue). (b) LH pane shows modelled (black) and observed (blue) TP concentrations at the NBC outlet. RH pane shows the TP concentrations at the mitigated sub-catchment outlet (Scenarios 1–3 shown by colored lines, with baseline in blue). Only Event 3 is shown.</figcaption></figure>
In the left-hand panes of Figure 7a,b, the black lines represent observed concentrations and the blue lines represent the modelled concentrations (at the NBC outlet) from the baseline simulation. The model performed reasonably well in capturing the peak concentrations, although there are some errors in timing of the peaks (compared with the blue line). The model underpredicted the TP concentrations slightly.

In terms of P yields, Figure 8 shows the modelled yields of TP and SS from the baseline simulation alongside the yields from Scenarios 1–3 from the mitigated sub-catchment. The bars are split into three coloured segments, each representing the yield transported by each of the three flow pathways (SR = surface runoff, “Fast S/S” = fast subsurface, “Slow G/W” = slow groundwater). The sum of the three segments was thus the total yield from all three pathways added together. The results covered a six-month period in 2012. Therefore, these yields included periods of low flows in addition to events; however, based on field observations and the fluxes transported by the modelled flow pathways, surface runoff during events contributed to the vast majority of the total SS and P losses from this catchment.

<figure><img src="image_91.png" /><figcaption>**Figure 8.** Plots showing modelled yields of TP and SS exported from the mitigated sub-catchment under all three scenarios in comparison with the modelled yields from the baseline.</figcaption></figure>

In Scenario 1, the total yield of SS and TP reduced quite drastically due to the reduction in transport via the fast surface flow pathway. There was some pollution swapping into SRP (transported by the fast subsurface flow pathway) as a result of these interventions increasing the amount of P that could infiltrate through the upper soil layer.

The results from Scenarios 2 and 3, where there has been simulation of the removal of SS and P from the surface runoff pathway by mitigation RAFs, show that there was a considerable reduction in SS TP yields according to the model. This reduction was achieved by trapping SS and attached particulate P, hence also reducing the components of PUP and PRP transported by the surface runoff pathway by 40–80%. These results were for the mitigated sub-catchment (1.25 km²) only and the results scaled up to the NBC outlet are discussed below.

### 3.3. Comparison of Recent Sampling Campaign and Model Results

The most recent available results from the mitigation features were as follows [24]. The mass of TP and SS trapped by the five features ranged from 1.5 to 2.5 kg and 0.5 to 2.8 tonnes, respectively (total mass 9.8 kg of P and 6.5 tonnes of sediment). These totals represent the masses of TP and SS collected over a 9-month period ending in summer 2015. Expressed as a yield of TP (loads per unit
area over the Mitigation sub-catchment), the reductions were 0.06 kg·ha⁻¹ of TP and 0.04 t·ha⁻¹ of SS. At present, only 24% of the sub-catchment area is treated by these features (3% of the NBC).

Model simulations using the parameter combination that was described above in the NBC over a 6-month period in 2012 that included the three events in April–September estimated that the fluxes of P trapped by mitigation features would be 0.23 kg·ha⁻¹ of PUP and 0.06 kg·ha⁻¹ of PRP in Scenario 2, and 0.46 kg·ha⁻¹ of PUP and 0.11 kg·ha⁻¹ of PRP in Scenario 3. The ratio between these yields corresponds closely to the ratios between the values of the parameter K<sub>RE</sub> in Scenarios 2 and 3 (40% and 80%, respectively).

**4. Discussion**

It is important to remember that European agricultural catchments have been heavily modified by centuries of intensive farming practices such that the degree of attenuation and storage have been reduced from prehistoric times when the catchments were in pristine (forested) condition prior to anthropogenic modification.

One pertinent research question relates to the detectability of mitigation features at the catchment scale, which in this case was ten times larger than the sub-catchment where the features were located. The results in terms of reducing the mean and maximum concentrations of SS and TP during Event 3 were analysed and were shown in the upper pane of Figure 9 for the three scenarios. Note that the effects of adding the features on the modelled Qp at the NBC outlet were not evaluated for reasons given above.

<figure><img src="image_92.png" /><figcaption>**Figure 9.** Plots showing the percentage reductions in TP and SS concentrations modelled at the NBC outlet achieved under: LH pane shows the three scenarios (legend indicates the scenario number) relative to the baseline in Event 3 only, and the RH pane shows the improvements under Scenario 3 relative to the baseline only, with Events 1–3 shown for comparison (improvements in Events 1–3 shown with different fill patterns).</figcaption></figure>

Scenario 1 evaluated the soil management improvements covering 10% of the total catchment to reduce surface runoff and particulate P loads transported by this pathway. The event mean and maximum TP concentrations were reduced by 4.7% and 6.4%, respectively, due to reductions in the PUP concentrations transported by the surface runoff pathway. The reductions in the event maximum TP concentrations were 4.3–5.9% during Event 3 at the NBC outlet when both attenuation storage and the removal of P were simulated in Scenarios 2 and 3. The reductions in Events 1 and 2 were similar under the three scenarios to Event 3 and are shown in the lower pane of Figure 9 for Scenario 3 only. These reductions are likely to be conservative estimates since the effect of the in-stream routing and dispersion along the channel reach that connected the mitigated sub-catchment outlet to the main channel outlet were not incorporated into the model structure.

In Scenario 1, there is considerable scope however for reducing the maximum TP concentration and TP yields by this method, as shown by the reduction in the maximum TP concentration in Event 1 of over 6% (Figure 9). In Scenario 2, the value of K<sub>RE</sub> of 40% was quite conservative compared to some
of the recent findings from U.K. field studies [11,12,24,45], but the load reduction was still predicted to be up to 4% at the catchment outlet with a smaller predicted decrease in the mean and maximum event TP concentration (less than those decreases achieved by Scenarios 1 (TP) and 3 (TP)). The model results showed that adding additional attenuation and removing P from 10% of the catchment reduced event loads of TP by up to 8% at the outlet under Scenario 3 (where $K_{RE}$ was 80%), which was the best performing of the three scenarios.

In terms of locating mitigation features of size ca. 2000–8000 m³ in one of the headwater sub-catchments of the NBC, the volumes required in Scenarios 2 and 3 do not seem excessively large if one assumes that a maximum water depth of 1 m is permissible under the relevant legislation (thus requiring 0.2 and 0.8 ha of the catchments, respectively, to be set aside for storage). In the U.K., stringent government guidelines [21] apply to all features greater than 10,000 m³ capacity, which may make constructing larger features expensive and probably uneconomic in terms of a cost–benefit ratio. For comparison, the five recent pilot features constructed in the NBC Mitigation sub-catchment had a total storage capacity of circa 400 m³ and treated only 38.4 ha (less than 0.031% of the total NBC catchment by area). It is important to consider that RAFs should be constructed in other parts of the NBC in order to achieve the desired reduction in P concentrations and loads.

It should be stressed that further research is required to evaluate the performance of these mitigation features when constructed at a larger scale than previously, i.e., up to 4–8 times greater than at the scale (17–34 ha) of Belford Burn [11,12,45] and up to 3 times greater than the 38 ha of the NBC Mitigation sub-catchment [24]. It is vital that trapped material is recovered by removal of sediment in well-designed traps or by removing vegetation.

Field studies are required to ascertain whether Scenario 1 could deliver the results shown here in reality; however, subsurface field drains are known to be a source of high SRP loads from agricultural catchments. Under Scenario 2, the strategic location of ditches and small channels in the landscape may play a crucial role in future management option for surface runoff driven agricultural pollution. The space in and around ditches afforded by buffer zones and fenced off channels could allow the construction of long riparian mitigation zones without significant impact on farming. The ability to address both NFM and pollution targets and create new ecological habitats may justify the investment and maintenance of RAFs [6]. The implication from all three scenarios is that an extensive network of RAFs is needed in farmed landscapes, and this could require a considerable shift in standard farming and environmental payments for schemes to be taken up and maintenance to be ensured. In reality, a mixture of all three scenarios being delivered would be attractive and that will very much depend on the local conditions and the farming community.

It is possible that unforeseen constraints, such as planning regulations, infrastructure restrictions and land ownership issues, may further restrict the widespread adoption of mitigation measures [4,22]. Therefore, the exact location of these flow pathways and identification of suitable sites for constructing features is beyond the scope of this study.

Further reductions in loads could be achieved in any case by targeting a larger percentage of the catchment for the construction of mitigation features, which in rural, pastoral farming-based areas like Cumbria (where the EdenDTC project is based) may be achievable, especially on a seasonal basis under the right conditions and government policies to compensate farmers. However, this policy might not be so attractive if the catchment contains arable land of higher value to the farmers than the rough and improved pasture in the upper Eden [22].

## 5. Conclusions

An attenuation component has been added to the CRAFT model to represent the storage and attenuation of surface runoff during events and also the trapping of particulate forms of P and suspended sediment, which provides a methodology for modelling NBS.

Currently, the removal efficiency of the trapping process has to be supplied by the user based on expert knowledge. The effect of adding attenuation and trapping SS plus particulate forms of P has
been modelled in a small (1.25 km²) sub-catchment. Due to the variability in the removal efficiencies measured in the field from previous studies of NBS, two modelling scenarios with removal efficiencies of 40% and 80% have been evaluated here with the latter clearly providing the best results in terms of improvements in reducing sediment and nutrient (P) fluxes and concentrations.

According to the model results, at the outlet of the entire Newby Beck catchment (ca. 10 × larger), a small reduction in the peak SS and TP concentrations would be observable during an event where surface runoff is the dominant pathway. A third scenario evaluated the impacts of reducing surface runoff to represent improved land management. This scenario reduced the amount of SS and P transported by the surface runoff pathway but increased the amount of SRP transported by the fast subsurface flow pathway. The reduction in both SS fluxes and maximum event concentration was more pronounced since the majority of the modelled flux was predicted to occur via the surface runoff pathway.

This modelling study represents a first step towards developing and testing a fully integrated model that can simulate both natural and added attenuation storage. Mitigation features are represented in the model as an aggregated storage effect rather than through utilizing a physically-based model of each feature individually. Research into the implications of adding new attenuation features to catchments with lighter soils and groundwater dominated runoff as well as adoption and maintenance issues is also needed.

One final conclusion arising from this study is the need to provide guidance to end users. The hypothetical attenuation versus mitigation curve can inform a policymaker on how to set annual targets over time that dictates the total amount of attenuation added to a catchment over a longer period. Hence, these guidelines would indicate how long it would take to reach the environmental target(s) based on an annual financial budget.

**Author Contributions:** Conceptualization, R.A. and P.Q.; Data curation, N.B.; Formal analysis, R.A.; Investigation, R.A.; Methodology, R.A., P.Q., and N.B.; Software, R.A.; Writing—original draft, R.A.; Writing—review and editing, R.A., P.Q., N.B., and S.R.

**Funding:** The Eden Demonstration Test Catchment (Eden DTC) research platform that funded this research and supplied the monitoring data used in this paper is funded by the Department for Environment, Food and Rural Affairs (DEFRA) (project WQ0210) and is further supported by the Welsh Assembly Government and the Environment Agency.

**Acknowledgments:** The Eden DTC (www.edendtc.org.uk) team includes research staff from Lancaster University, University of Durham, Newcastle University (PROACTIVE team), Eden Rivers Trust, Centre for Ecology and Hydrology, British Geological Survey, and Newton Rigg College, and their assistance in undertaking this research is gratefully acknowledged.

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

## References

1. NWRM—Natural Water Retention Measures. European Union. Available online: www.nwrm.eu (accessed on 9 February 2018).

2. Environment Agency Working with Natural Processes to Reduce Flood Risk. Available online: www.gov.uk/government/publications/working-with-natural-processes-to-reduce-flood-riskPublished31/10/2017 (accessed on 9 February 2018).

3. O'Connell, P.E.; Ewen, J.; O'Donnell, G.; Quinn, P. Is there a link between agricultural land-use management and flooding? *Hydrol. Earth Syst. Sci.* **2007**, *11*, 96–107. [CrossRef]

4. Nisbet, T.R.; Marrington, S.; Thomas, H.; Broadmeadow, S.B.; Valatin, G. *Slowing the Flow at Pickering. Final Report for the Department of Environment, Food and Rural Affairs*; Project RMP5455; DEFRA: London, UK, 2011.

5. Quinn, P.; O'Donnell, G.; Nicholson, A.; Wilkinson, M.; Owen, G.; Jonczyk, J.; Barber, N.; Hardwick, M.; Davies, G. *Potential Use of Runoff Attenuation Features in Small Rural Catchments for Flood Mitigation*; Newcastle University, Environment Agency, Royal Haskoning DHV: Newcastle Upon Tyne, UK, 2013.

6. Charlesworth, S.M.; Booth, C.A. *Sustainable Surface Water Management: A Handbook for SUDS*; John Wiley & Sons: Chichester, UK, 2016.
7. SEPA. *Natural Flood Management Handbook*, 2nd ed.; Scottish Environment Protection Agency: Stirling, UK, 2015; Available online: https://www.sepa.org.uk/media/163560/sepa-natural-flood-management-handbook1.pdf (accessed on 9 February 2018).

8. Duffy, A.; Moir, S.; Berwick, N.; Shabashow, J.; D'Arcy, B.; Wade, R. *Rural Sustainable Drainage Systems: A Practical Design and Build. Guide for Scotland's Farmers and Landowner*; CREW: Scotland, UK, 2015; Available online: http://crew.ac.uk/publications (accessed on 6 June 2018).

9. Wilkinson, M.E.; Quinn, P.F.; Welton, P. Runoff management during the September 2008 floods in the Belford catchment, Northumberland. *J. Flood Risk Manag.* **2010**, *3*, 285–295. [CrossRef]

10. Nisbet, T.R.; Thomas, H. *Restoring Floodplain Woodland for Flood Alleviation. Final Report for the Department of Environment, Food and Rural Affairs*; Project SLD2316; DEFRA: London, UK, 2008.

11. Barber, N. Sediment, Nutrient and Runoff Management and Mitigation in Rural Catchments. Ph.D. Thesis, Newcastle University, Newcastle Upon Tyne, UK, 2014.

12. Barber, N.J.; Quinn, P.F. Mitigating diffuse water pollution from agriculture using soft engineered runoff attenuation features. *Area* **2012**, *44*, 454–462. [CrossRef]

13. Deasy, C.; Quinton, J.N.; Silgram, M.; Bailey, A.P.; Jackson, B.; Stevens, C.J. Contributing understanding of mitigation options for phosphorus and sediment to a review of the efficacy of contemporary agricultural stewardship measures. *Agric. Syst.* **2010**, *103*, 105–109. [CrossRef]

14. Doody, D.G.; Archbold, M.; Foy, R.H.; Flynn, R. Approaches to the implementation of the Water Framework Directive: Targeting mitigation measures at critical source areas of diffuse phosphorus in Irish catchments. *J. Environ. Manag.* **2012**, *93*, 225–234. [CrossRef] [PubMed]

15. Wilcock, R.J.; Müller, K.; van Assema, G.B.; Bellingham, M.A.; Ovenden, R. Attenuation of nitrogen, phosphorus and E. coli inputs from pasture runoff to surface waters by a farm wetland: The importance of wetland shape and residence time. *Water Air Soil Pollut.* **2012**, *223*, 499–509. [CrossRef]

16. Avery, L.M. *Rural Sustainable Drainage Systems (RSuDS)*; Environment Agency: Bristol, UK, 2012.

17. Blöschl, G. Scaling in hydrology. *Hydrol. Process.* **2001**, *15*, 709–711. [CrossRef]

18. Biggs, J.; Stoate, C.; Williams, P.; Brown, C.; Casey, A.; Davies, S.; Diego, I.G.; Hawczak, A.; Kizuka, T.; McGoff, E.; et al. *Water Friendly Farming Autumn 2016 Update*; Freshwater Habitats Trust: Oxford, UK; Game & Wildlife Conservation Trust: Fordingbridge, UK, 2016; Available online: https://freshwaterhabitats.org.uk/wp-content/uploads/2016/11/Water-Friendly-Farming-update-2016.pdf (accessed on 18 April 2018).

19. Fealy, R.M.; Buckley, C.; Mechan, S.; Melland, A.; Mellander, P.E.; Shortle, G.; Wall, D.; Jordan, P. The Irish Agricultural Catchments Programme: Catchment selection using spatial multi-criteria decision analysis. *Soil Use Manag.* **2010**, *26*, 225–236. [CrossRef]

20. Owen, G.J.; Perks, M.T.; Benskin, C.M.H.; Wilkinson, M.E.; Jonczyk, J.; Quinn, P.F. Monitoring agricultural diffuse pollution through a dense monitoring network in the River Eden Demonstration Test Catchment, Cumbria, UK. *Area* **2012**, *44*, 443–453. [CrossRef]

21. McGonigle, D.F.; Burke, S.P.; Collins, A.L.; Gartner, R.; Haft, M.R.; Harris, R.C.; Haygarth, P.M.; Hedges, M.C.; Hiscock, K.M.; Lovett, A.A. Developing demonstration test catchments as a platform for transdisciplinary land management research in England and Wales. *Environ. Sci. Process. Impacts* **2014**, *16*, 1618–1628. [CrossRef] [PubMed]

22. Metcalfe, P.; Beven, K.; Hankin, B.; Lamb, R. A modelling framework for evaluation of the hydrological impacts of nature-based approaches to flood risk management, with application to in-channel interventions across a 29 km² scale catchment in the United Kingdom. *Hydrol. Process.* **2017**, *31*, 1734–1748. [CrossRef]

23. Bechmann, M.; Deelstra, J.; Stálnacke, P.; Eggestad, H.O.; Øygarden, L.; Pengerud, A. Monitoring catchment scale agricultural pollution in Norway: Policy instruments, implementation of mitigation methods and trends in nutrient and sediment losses. *Environ. Sci. Policy* **2008**, *11*, 102–114. [CrossRef]

24. Barber, N.J.; Reaney, S.N.; Barker, P.A.; Benskin, C.; Burke, S.; Cleasby, W.; Haygarth, P.; Jonczyk, J.C.; Owen, G.J.; Snell, M.A.; et al. The treatment train approach to reducing nonpoint source pollution from agriculture. In Proceedings of the AGU Fall Meeting, San Francisco, CA, USA, 12–16 December 2016.

25. Kay, P.; Edwards, A.C.; Foulger, M. A review of the efficacy of contemporary agricultural stewardship measures for ameliorating water pollution problems of key concern to the UK water industry. *Agric. Syst.* **2009**, *99*, 67–75. [CrossRef]
26. Ensign, S.H.; McMillan, S.K.; Thompson, S.P.; Piehler, M.F. Nitrogen and phosphorus attenuation within the stream network of a coastal, agricultural watershed. *J. Environ. Qual.* **2006**, *35*, 1237–1247. [CrossRef] [PubMed]

27. Adams, R.; Quinn, P.F.; Bowes, M.J. The Catchment Runoff Attenuation Flux Tool, A minimum information requirement nutrient pollution model. *Hydrol. Earth Syst.* **2015**, *19*, 1641–1657. [CrossRef]

28. Adams, R.; Quinn, P.F.; Perks, M.; Barber, N.J.; Jonczyk, J.; Owen, G.J. Simulating high frequency water quality monitoring data using a catchment runoff attenuation flux tool (CRAFT). *Sci. Total Environ.* **2016**, *572*, 1622–1635. [CrossRef] [PubMed]

29. Heathwaite, A.L.; Quinn, P.F.; Hewett, C.J.M. Modelling and managing critical source areas of diffuse pollution from agricultural land using flow connectivity simulation. *J. Hydrol.* **2005**, *304*, 446–461. [CrossRef]

30. Singh, V.P. Is Hydrology kinematic? *Hydrol. Process.* **2002**, *16*, 667–716. [CrossRef]

31. Perks, M.T.; Owen, G.J.; Benskin, C.M.H.; Jonczyk, J.; Deasy, C.; Burke, S.; Haygarth, P.M. Dominant mechanisms for the delivery of fine sediment and phosphorus to fluvial networks draining grassland dominated headwater catchments. *Sci. Total Environ.* **2015**, *523*, 178–190. [CrossRef] [PubMed]

32. Ockenden, M.C.; Deasy, C.E.; Benskin, C.M.; Beven, K.J.; Burke, S.; Collins, A.L.; Evans, R.; Falloon, P.D.; Forber, K.J.; Hiscock, K.M.; et al. Potential effects of changing climate on catchment processes and nutrient transfers: Evidence from high temporal concentration-flow dynamics in headwater catchments. *Sci. Total Environ.* **2016**, *548*, 325–339. [CrossRef] [PubMed]

33. Ockenden, M.C.; Tych, W.; Beven, K.; Collins, A.; Evans, R.; Falloon, P.; Forber, K.; Hiscock, K.; Hollaway, M.; Kahana, R.; et al. Prediction of storm transfers and annual loads with data-based mechanistic models using high-frequency data. *Hydrol. Earth Syst.* **2017**, *18*, 6425–6444. [CrossRef]

34. Mills, C.F.; Bathurst, J.C. Spatial variability of suspended sediment yield in a gravel-bed river across four orders of magnitude of catchment area. *Catena* **2015**, *133*, 14–24. [CrossRef]

35. European Commission. *Directive 2000/60/EC: Establishing a Framework for Community Action in the Field of Water Policy.* (*The Water Framework Directive*); European Commission: Brussels, Belgium, 2000.

36. Greene, S.; Johnes, P.J.; Bloomfield, J.P.; Reaney, S.M.; Lawley, R.; Elkhatib, Y.; Freer, J.; Odoni, N.; Macleod, C.J.; Percy, B. A geospatial framework to support integrated biogeochemical modelling in the United Kingdom. *Environ. Model. Softw.* **2015**, *68*, 219–232. [CrossRef]

37. Hewett, C.J.; Quinn, P.F.; Heathwaite, A.L.; Doyle, A.; Burke, S.; Whitehead, P.G.; Lerner, D.N. A multi-scale framework for strategic management of diffuse pollution. *Environ. Model. Softw.* **2009**, *24*, 74–85. [CrossRef]

38. Lane, S.N.; Reaney, S.M.; Heathwaite, A.L. Representation of landscape hydrological connectivity using a topographically driven surface flow index. *Water Resour. Res.* **2009**, *45*. [CrossRef]

39. Reaney, S.M.; Lane, S.N.; Heathwaite, A.L.; Dugdale, L.J. Risk-based modelling of diffuse land use impacts from rural landscapes upon salmonid fry abundance. *Ecol. Model.* **2011**, *222*, 1016–1029. [CrossRef]

40. WFD-UKTAG. *UKTAG River Assessment Method—Phosphorus: River Phosphorus Standards*; Water Framework Directive-United Kingdom Technical Advisory Group (WFD-UKTAG): Stirling, UK, 2014.

41. Kovacs, A.; Honti, M.; Zessner, M.; Eder, A.; Clement, A.; Blöschl, G. Identification of phosphorus emission hotspots in agricultural catchments. *Sci. Total Environ.* **2012**, *433*, 74–88. [CrossRef] [PubMed]

42. Schoumans, O.F.; Chardon, W.J.; Bechmann, M.E.; Gascuel-Odoux, C.; Hofman, G.; Kronvang, B.; Dorioz, J.M. Mitigation options to reduce phosphorus losses from the agricultural sector and improve surface water quality: A review. *Sci. Total Environ.* **2014**, *468*, 1255–1266. [CrossRef] [PubMed]

43. Beven, K. On the generalized kinematic routing method. *Water Resour. Res.* **1979**, *15*, 1238–1242. [CrossRef]

44. Wood, E.F.; Sivapalan, M.; Beven, K.; Band, L. Effects of spatial variability and scale with implications to hydrologic modeling. *J. Hydrol.* **1988**, *102*, 29–47. [CrossRef]

45. Wilkinson, M.E.; Quinn, P.F.; Barber, N.J.; Jonczyk, J. A framework for managing runoff and pollution in the rural landscape using a Catchment Systems Engineering approach. *Sci. Total Environ.* **2014**, *468*, 1245–1254. [CrossRef] [PubMed]

<figure><img src="image_93.png" /></figure>

© 2018 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/).
<figure><img src="image_94.png" /></figure>

Article

# Shear Stress-Based Analysis of Sediment Incipient Deposition in Rigid Boundary Open Channels

Necati Erdem Unal*

Department of Civil Engineering, Hydraulics Division, Istanbul Technical University, Maslak, 34469 Istanbul, Turkey; neu@itu.edu.tr; Tel.: +90-212-285-3727

Received: 24 September 2018; Accepted: 4 October 2018; Published: 9 October 2018

<figure><img src="image_95.png" /></figure>

check for
updates

**Abstract:** Urban drainage and sewer systems, and channels in general, are treated by the deposition of sediment that comes from water collecting systems, such as roads, parking lots, land, cultivation areas, and so forth, which are all under gradual or sudden change. The carrying capacity of urban area channels is reduced heavily by sediment transport that might even totally block the channel. In order to solve the sedimentation problem, it is therefore important that the channel is designed by considering self-cleansing criteria. Incipient deposition is proposed as a conservative method for channel design and is the subject of this study. With this aim, an experimental study carried out in trapezoidal, rectangular, circular, U-shape, and V-bottom channels is presented. Four different sizes of sand were used as sediment in the experiments performed in a tilting flume under nine different longitudinal channel bed slopes. A shear stress approach is considered, with the Shields and Yalin methods used in the analysis. Using the experimental data, functionals are developed for both methods. It is seen that the bed shear stress changes with the shape of the channel cross-section. Incipient deposition in rectangular and V-bottom channels starts under the lowest and the highest shear stress, respectively, due mainly to the shape of the channel cross-section that affects the distribution of shear stress on the channel bed.

**Keywords:** incipient deposition; sediment transport; self-cleansing; sewer systems; shear stress; urban drainage system

## 1. Introduction

The sediment transport issue has always been an important scientific and practical problem [1,2], and kept its importance with the changes in hydrology [3–5] emerging with the change in the sediment load of urban watersheds and alluvial streams. As the outlets of urban watersheds, the sewer and urban drainage systems are heavily affected by any change in their watersheds. Therefore, research on sediment transport is continuously needed for sustainable practice in drainage systems.

Sediment deposition is avoided as it causes numerous unwanted problems in urban drainage and sewer systems. It reduces the hydraulic capacity of the channel by decreasing the flow cross-sectional area or blocking the channel. The performance and efficiency of drainage systems is heavily affected by deposition. Additional funds should be invested to keep the system working. Furthermore, sedimentation creates environmental problems. Such problems in drainage and sewer systems could be prevented or minimized by the use of self-cleansing criteria, with which sediment particles deposited at the channel bed start to move [6], or sediment particles suspended within the flow are transported without being deposited [7–9]. Transportation of sediment particles without deposition is preferred in drainage system design to keep the channel bed clean [10–14]. In this regard, incipient deposition is a concept linked to the channel design, and identified as the sediment transport mode in which sediment particles are clustered visibly in certain areas at the channel bed [15]. In other words, sediment particles are transported as bed load or accumulated at the channel bed, but without making a permanent
deposited bed layer [15–17]. Flow velocity is sufficiently low for incipient deposition of sediment. At the incipient deposition, sediment particles in suspension within flow start moving downward to reach the channel bed.

Incipient deposition has been studied in several fixed bed channels by Loveless [15] who assumed that incipient motion and incipient deposition were similar concepts but with a slight difference [18]. The experimental data of Loveless [15] fit the sediment transport models of May [19] and Ackers [20]. Safari et al. [16] studied the incipient deposition concept using the experimental data of Loveless [15]. As a conclusion, velocity at incipient deposition was found to be higher than velocity at incipient motion for non-cohesive sediments. Incipient deposition is observed when flow velocity decreases gradually to a level that allows sediment particles to deposit. In the opposite case, when flow velocity increases gradually, sediment particles with no motion start moving when flow velocity reaches a level high enough for incipient motion. The former is called the incipient deposition velocity, which is higher than the latter, the incipient motion velocity. This is a hysteretic curve with a higher threshold velocity for incipient deposition and lower threshold velocity for the incipient motion. In this paradigm, there is a shear stress (or velocity) threshold below which erosion does not occur, and a lower threshold above which deposition does not occur; erosion and deposition occur simultaneously between the two thresholds [21]. Therefore, results denied the common assumption that incipient deposition and incipient motion are the same.

In order to explain the exact difference between the incipient deposition and incipient motion, Aksoy and Safari [22] performed a preliminary study on incipient motion and incipient deposition in a trapezoidal cross-section channel, and found that the flow has higher shear stress at incipient deposition than at the incipient motion. For the sake of achieving more conclusive results, a new set of experiments with a wider range of sediment size in different channel cross-sections seemed important. Not only because of its importance, but also for the sake of getting experimental data with a wider range in terms of sediment size, channel cross-section, channel slope, discharge, and so forth, Unal et al. [23] constructed a laboratory experimental setup to study the incipient motion and incipient deposition in trapezoidal, rectangular, circular, U-shape, and V-bottom channels, and performed experiments for the self-cleansing design of fixed bed systems.

In this study, experimental data from an indoor laboratory flume was analyzed to understand the incipient deposition of sediment particles within flow. The shear stress approach was considered for the analysis in which the Shields [24] and Yalin [25] methods were used.

## 2. Mechanism of Particle Motion and Methodology

Sediment particles in flow move under the influence of two types of hydrodynamic forces; the first of which has a positive impact through the drag force and the lift force, while the second discourages motion through the buoyed weight of sediment and the resisting force against motion. The drag force should be equal to the resistance force in the sediment threshold condition.

The shear stress- and velocity-based approaches were commonly used in the analysis of sediment threshold and incipient deposition in this study. The velocity-based approach has been applied on the experimental data existing in the literature [26]. The shear stress approach was used in this study for which Shields [24] and Yalin [25] methods are considered.

### 2.1. Shields Method

The shear stress approach used in the incipient deposition of sediment is based on the shear velocity ($u_*$), defined as

$$u_* = \sqrt{\frac{\tau_{id}}{\rho}} \qquad (1)$$
in which $\tau_{id}$ is the bed shear stress under the incipient deposition condition and $\rho$ is the specific mass of water. The dimensionless shear stress is calculated by

$$ \tau_{id}^* = \frac{\tau_{id}}{\rho g d(s-1)}. \qquad (2) $$

in which $g$ is the acceleration due to gravity, $d$ is the median size of sediment particles, and $s$ is the sediment relative mass density. The dimensionless shear stress ($\tau_{id}^*$) is indicated as

$$ \tau_{id}^* = f(Re^*) \qquad (3) $$

where $Re^*$ is the particle Reynolds number ($Re^*$) defined by

$$ Re^* = \frac{u_* d}{\nu} \qquad (4) $$

in which $\nu$ is the kinematic viscosity of water.

Experimental data has been used by many researchers [9,15–17,22,23,26] to determine the functional between $\tau_{id}^*$ and $Re^*$ for practical problems of incipient motion [27,28]. In this study, the same methodology was adopted for the incipient deposition.

Average flow velocity ($V_{id}$) calculated from the incipient deposition experimental data is used to calculate the incipient deposition shear stress by

$$ \tau_{id} = \frac{\lambda \rho V_{id}^2}{8} \qquad (5) $$

in which $\lambda$ is the channel friction factor to be calculated by the Colebrook–White equation [29] as

$$ \lambda = \frac{1}{4 \left[ \log \left( \frac{k_b}{14.8R} + \frac{0.22\nu}{R\sqrt{gRS}} \right) \right]^2} \qquad (6) $$

in which $k_b$ is the roughness height of the bed taken to be the same as the median size of sediment particles ($d$), $R$ is the hydraulic radius of the channel, and $S$ is the slope of the channel bed.

## 2.2. Yalin Method

Yalin [25] suggested a combination of the dimensionless parameters initially proposed by Shields [24] as

$$ \tau_{id}^* = f(D_{gr}) \qquad (7) $$

in which $D_{gr}$ is the dimensionless grain size parameter defined by

$$ D_{gr} = \left[ \frac{(s-1)g d^3}{\nu^2} \right]^{1/3} \qquad (8) $$

The shear velocity is eliminated, and only fluid and sediment characteristics are retained in the formulation.

## 3. Experiments

An experimental setup was configured as in Figure 1 [22,23]. An iron-made support structure was constructed. Twelve meter-long transparent acrylic glass (plexiglass) channels were mounted on the support structure. Five different cross-sections were considered for the channels; they were trapezoidal, rectangular, circular, U-shape, and V-bottom. The surface width of the rectangular, U-shape, and V-bottom channels was 300 mm, while the trapezoidal channel had the same width at
the bottom, and outer angles of 60° at the 30 cm-long side walls. The U-shape and V-bottom channels had a cross-fall of 50 mm longitudinally along the centerline of the bottom. The inner diameter of the circular channel was 290 mm.

<figure><img src="image_96.png" /><figcaption>**Figure 1.** Experimental setup and cross-sections of the channels.</figcaption></figure>

In the experiments, the bed slope was changed between 0.00147 and 0.01106. Four non-cohesive sands were poured into the channel from the sediment feeder placed 3 m upstream of the observation section of the channel. The granulometric curve and the characteristics of the sands are shown in Figure 2 and Table 1, from which it can be seen that they have uniform size distribution. The discharge was measured by an ultrasonic flowmeter (BSUF-TTCL, Bass Instruments, Istanbul, Turkey) with an accuracy better than 1.0% of read. Sediment motion was observed in each experiment in the 4 m-long observation section of the channel, 4 m from the inlet and 4 m from the outlet of the channel (Figure 1). Uniform flow conditions were satisfied in the channel before observations and measurements were done.

<figure><img src="image_97.png" /><figcaption>**Figure 2.** Granulometric curve of sediment.</figcaption></figure>
<table><caption>**Table 1.** Sediment characteristics. *d*: median diameter, *s*: relative mass density, *σ*<sub>g</sub>: geometric standard deviation of sediment particles.</caption><tr><td>d (mm)</td><td>0.15</td><td>0.58</td><td>1.08</td><td>1.52</td></tr><tr><td>s</td><td>2.60</td><td>2.63</td><td>2.56</td><td>2.60</td></tr><tr><td>σ<sub>g</sub></td><td>1.3</td><td>1.6</td><td>1.3</td><td>1.4</td></tr></table>

Experiments started with a flow velocity high enough to achieve the non-deposition condition. The average velocity was adjusted by increasing or decreasing flow discharge into the channel. In the non-deposition condition, sediment particles are prevented from being deposited; that is, sediment particles within flow are in motion. Flow velocity was gradually decreased until incipient deposition was achieved; that is, flow switches from non-deposition to incipient deposition. It was assumed that incipient deposition was satisfied when sediment particles were clustered visibly in certain areas at the channel bottom [15]. In this case, flow velocity is sufficiently low for incipient deposition of sediment. Incipient deposition was observed in the same form in trapezoidal and rectangular channels due to their flat bed. However, in the channels with U-shape and V-bottom, sediment particles were deposited in the center line of the channel with the same cross-fall. The form of the deposition depends on the channel bed. Sediment particles were accumulated on each other along the narrow centerline in the V-bottom channel, while accumulation in the circular and U-shape channels was not that narrow, as due to the wider bed along the centerline width, sediment particles spread over the bed width to make a deposited sediment layer instead.

# 4. Results

## 4.1. Shields Method

The incipient deposition experimental data of the channels are plotted on the Shields diagram (Figure 3), with the upper and lower limits as proposed by Paphitis [30]. Using the experimental data, $\tau_{id}^*$ and $Re^*$ were calculated using Equations (2) and (4), respectively, and functional relationships were developed by curve fitting to the measured data as

$$ \tau_{id}^* = 0.74(Re^*)^{-0.86} \qquad 3.13 < Re^* < 47.61 \qquad r^2 = 0.913 \qquad (9) $$

$$ \tau_{id}^* = 0.32(Re^*)^{-0.97} \qquad 2.36 < Re^* < 29.89 \qquad r^2 = 0.953 \qquad (10) $$

$$ \tau_{id}^* = 0.57(Re^*)^{-0.83} \qquad 3.02 < Re^* < 41.11 \qquad r^2 = 0.972 \qquad (11) $$

$$ \tau_{id}^* = 0.55(Re^*)^{-0.72} \qquad 3.19 < Re^* < 48.11 \qquad r^2 = 0.960 \qquad (12) $$

$$ \tau_{id}^* = 0.79(Re^*)^{-0.78} \qquad 3.46 < Re^* < 52.62 \qquad r^2 = 0.907 \qquad (13) $$

in the range $2.36 < Re^* < 52.62$. for the trapezoidal, rectangular, circular, U-shape, and V-bottom channels, respectively. In Equations (9)–(13), $r^2$, the determination coefficient, shows the goodness of fit of the curves. The data of the five channels are close to each other and partially overlap in some of the cases. It is seen from Figure 3 that, for the 0.15 mm- and 0.58 mm-particle size sands, incipient deposition shear stress remained above the upper limit of the Shields curve for the non-rectangular channels; however, in the rectangular channel, it is on the upper limit curve for the finest sand (the 0.15 mm-particle size sand) and on the average curve for the 0.58 mm-particle size sand. It should be kept in mind that the Shields curve has been developed for the incipient motion of sediment in loose boundary channels. Therefore, Figure 3 indicates, for sand finer than 0.58 mm, that the incipient deposition shear stress in rigid boundary channels is higher than the incipient motion shear stress in loose boundary channels. For sand with a 1.08 mm diameter, the incipient deposition shear stress remains between the upper limit and the average curve in the non-rectangular channels. It is on the lower limit curve in the rectangular channel case. It is also shown in Figure 3 that the incipient
deposition shear stress for the coarsest sand with a 1.52 mm diameter remains between the average and lower limit curves in all the channels other than the rectangular cross-section. It is below the lower limit curve for the rectangular channel. Generally, the incipient deposition shear stress of coarse sediment (1.08 mm and 1.52 mm) is lower than incipient motion shear stress in loose boundary channels. The coarser the sediment, the lower the shear stress under which the sediment particles initiate deposition within flow.

<figure><img src="image_98.png" /><figcaption>Figure 3. Relationship between the incipient deposition shear stress and Reynolds number based on the Shields method (circles show values calculated from the measurements, lines are the fitted equations).</figcaption></figure>

It can also be seen that incipient deposition starts under lower shear stress in the rectangular channel (Figure 3). The incipient deposition shear stress of the trapezoidal, circular, U-shape, and V-bottom channels are close to each other. However, sediment particles in the V-bottom and U-shape channels initiate deposition under higher shear stress. The rectangular channel on the other hand obviously has different performance than the other channels; it allows sediment to move within flow until a lower shear stress is approached.

## 4.2. Yalin Method

Incipient deposition experimental data are plotted on the Yalin diagram in Figure 4. For the five channels, the dimensionless incipient deposition shear stress ($\tau_{id}^*$) and grain size ($D_{gr}$) were calculated by Equations (2) and (8), respectively. Utilizing the incipient deposition experimental data,

$$ \tau_{id}^* = 0.84D_{gr}^{-0.92} \qquad r^2 = 0.913 \qquad (14) $$

$$ \tau_{id}^* = 0.47D_{gr}^{-0.99} \qquad r^2 = 0.953 \qquad (15) $$

$$ \tau_{id}^* = 0.68D_{gr}^{-0.89} \qquad r^2 = 0.972 \qquad (16) $$

$$ \tau_{id}^* = 0.65D_{gr}^{-0.80} \qquad r^2 = 0.960 \qquad (17) $$

$$ \tau_{id}^* = 0.87D_{gr}^{-0.86} \qquad r^2 = 0.907 \qquad (18) $$

are proposed in the range of 3.76 < $D_{gr}$ < 38.06 for the trapezoidal, rectangular, circular, U-shape, and V-bottom channels, respectively. It can be seen that the incipient deposition shear stress remains above the upper limit of the Yalin curve for the trapezoidal, U-shape, and V-bottom channels. It is on the upper limit for the fine sand ($d$ = 0.15 mm) in the rectangular channel. For the medium sand ($d$ = 0.58 mm) it is on and below the average curve in the circular and rectangular channels, respectively. For the coarser sand ($d$ = 1.08 mm), the incipient deposition shear stress remains between the upper limit and the average curve in the trapezoidal, U-shape, and V-bottom channels. It is below the lower
limit for the rectangular channel, and on the average curve for the circular channel. For the coarsest sand ($d = 1.52$ mm), the incipient deposition shear stress is between the average curve and the lower limit in the channels with no rectangular cross-section. It is below the lower limit of the curve for the rectangular channel.

<figure><img src="image_99.png" /><figcaption>**Figure 4.** Relationship between the incipient deposition shear stress and grain size based on the Yalin method (circles show values calculated from the measurements, lines are the fitted equations).</figcaption></figure>

The incipient deposition shear stress calculated from the experimental data (Figure 4) shows that sediment particles within flow in a rectangular channel start deposition under lower shear stress compared to the other channels. In other words, sediment particles within flow in the rectangular channel among the tested channels could be kept moving within flow under the lowest shear stress. The non-rectangular channels all have similar performance. Among the non-rectangular channels, the incipient deposition shear stress is lowest in the circular channel and the highest in the V-bottom channel. As a general result, it is clear that the coarser the sediment, the lower the shear stress under which sediment particles initiate deposition within flow. It can be said that the shape of the channel cross-section significantly affects the incipient deposition shear stress. Consequently, the incipient deposition shear stress is lower in the rectangular channel compared to the other channels.

**5. Discussion**

The results of the shear stress approach (Shields and Yalin methods) show that the channel cross-sectional shape significantly affects the incipient deposition shear stress. This is due to the change in the wall-normal component of the gravitational forces from which the friction force stems [31]. It is important to stress that the concept of incipient deposition is quite different to that of incipient motion. Not only the driving forces for the motion, but also the resisting forces against the motion are radically different in these two concepts. For the incipient motion, the friction coefficient is the static friction coefficient, which is accepted to be equal to the tangent of the friction angle between the bed and the grain (which becomes the angle of repose for a bed composed of identical particles). On the other hand, for incipient deposition, the friction coefficient is the dynamic friction coefficient, which is the tangent of the dynamic friction angle between the bed and the grain. Fredsoe and Deigaard [31] pointed out this important difference, and carried out an analysis for the correction of the transport rate of sediment particles moving on a transverse slope.

It should additionally be emphasized that, since the bed conditions in this study are smooth (i.e., the bed is the so-called "starving bed" where sediment particles travel on an otherwise hydraulically-smooth bed), one could expect that the changes with respect to the bed shape become even more pronounced compared to a regular bed composed of identical particles. Finally, it is possible to state that the critical shear stress curves for the incipient deposition for each of the given channel
cross-sections (trapezoidal, rectangular, circular, U-shape, and V-bottom) are expected to be different from each other, unless the channel is very wide, or unless a separate sophisticated correction for the critical shear stress for the deposition is carried out for each and every case. Such a correction, which was out of the scope of this study, should be expected to involve an averaging of the critical shear stress, possibly by integration across the channel cross-section.

The Shields and Yalin methods provide similar results. The behavior of the channels in terms of the magnitude of the incipient deposition shear stress does not change with the selected method. In both methods, the sediment particles initiate deposition under the lowest shear stress in the rectangular channel, and under the highest shear stress in the V-bottom channel. This is an expected result because the incipient deposition shear stress is calculated from the channel, sediment, and flow characteristics.

Another point worth discussing is the trend in the dimensionless incipient deposition shear stress against the particle Reynolds number in the Shields curve, and the grain size in the Yalin curve. The range of the experiments is $2.36 < Re^* < 52.62$ for the Shields curve and $3.76 < D_{gr} < 38.06$ for the Yalin curve. Within these ranges, both curves gradually decrease with increasing $Re^*$ and $D_{gr}$. Experimental observation fits the general character of the Shields and Yalin curves within the range of the experiments.

It should also be mentioned that the incipient deposition shear stress decreases with increasing channel bed slope. In other words, sediment particles move within flow without being deposited until the shear stress becomes low enough to initiate deposition. Deposition starts under lower shear stress when the channel has a steeper slope. This is an observation for all sediment sizes and channel cross-sections.

One more point to discuss is the similarity with the velocity-based analysis of the incipient deposition performed by Aksoy et al. [26] under which the Novak-Nalluri [32] and Yang [33] methods were used. The rectangular cross-section channel was marked with the lowest velocity at the incipient deposition. Therefore, based on either the shear stress or the velocity approach, rectangular channels seem to be preferable.

Finally, as a general discussion, when the shear stress approach is compared with the velocity-based approach [26], there is a dilemma in the selection of the appropriate approach—shear stress or velocity. In the shear stress approach, two dimensionless parameters are used, namely the shear stress and the particle Reynolds number, both of which are dependent on shear stress; its critical value cannot be determined directly, but can only be calculated through a trial-and-error technique implicitly. In the velocity approach, critical velocity is calculated by an explicit solution. This gives an advantage to the velocity approach in terms of computation. However, as Vanoni [34] demonstrated, in using the velocity approach, flow depth or hydraulic radius, and in general channel cross-section, must be specified.

## 6. Conclusions

Incipient deposition is a different concept as opposed to the common assumption that the incipient deposition and incipient motion are the same. Therefore, the incipient deposition was considered solely in this study. An experimental analysis was performed for fixed bed channels with trapezoidal, rectangular, circular, U-shape, and V-bottom cross-sections. Experimental data was analyzed using the shear stress approach, under which the Shields and Yalin methods were considered to calculate the incipient deposition shear stress. Both methods showed that the incipient deposition starts under lower shear stress in the rectangular channel compared to the non-rectangular channels. This indicates that rectangular channels have higher efficiency of sediment transport as sediment particles deposit at lower velocities. This gives an advantage to rectangular channels in the design of urban drainage and sewer systems or irrigation canals. The trapezoidal, circular, U-shape, and V-bottom channels have similar performance, among which the circular channel has the lowest incipient deposition shear stress and the V-bottom channel has the highest. This makes the circular channel the second most preferable channel after the rectangular channel, and the V-bottom channel the least. Analysis brings
the conclusion that the cross-section of the channel significantly affects the shear stress at incipient deposition. The general trends in the Shields and Yalin curves were traced with the experimental data within the data range tested. One observation is that the incipient deposition shear stress decreases with the increasing channel bed slope; that is, sediment deposition in the channel is delayed when the channel has a steeper slope. The outputs of this study are expected to be considered together with the available literature and be employed for practical use in rigid boundary channel design. Further experiments are encouraged to extend the range of sediment size such that the validity range of the developed equations is increased.

**Funding:** The Scientific and Technological Research Council of Turkey (TUBITAK) Project number 114M283, The Scientific Research Projects (BAP) Unit of Istanbul Technical University Project number 37973.

**Acknowledgments:** The author is thankful to the Reviewers, Guest Editors and Academic Editor for their important comments with which the paper has been well improved.

**Conflicts of Interest:** The author declares no conflict of interest.

## References

1. Graf, W.H.; Acaroglu, E.R. Sediment transport in conveyance systems (Part 1)/A physical model for sediment transport in conveyance systems. *Hydrol. Sci. J.* **1968**, *13*, 20–39. [CrossRef]

2. Bogardi, J. *Sediment Transport in Alluvial Streams*; Akademiai Kiado: Budapest, Hungary, 1974.

3. Montanari, A.; Young, G.; Savenije, H.H.G.; Hughes, D.; Wagener, T.; Ren, L.L.; Koutsoyiannis, D.; Cudennec, C.; Toth, E.; Grimaldi, S.; et al. "Panta Rhei—Everything flows": Change in hydrology and society—The IAHS scientific decade 2013–2022. *Hydrol. Sci. J.* **2013**, *58*, 1256–1275. [CrossRef]

4. Ceola, S.; Montanari, A.; Krueger, T.; Dyer, F.; Kreibich, H.; Westerberg, I.; Carr, G.; Cudennec, C.; Elshorbagy, A.; Savenije, H.; et al. Adaptation of water resources systems to changing society and environment: A statement by the International Association of Hydrological Sciences. *Hydrol. Sci. J.* **2016**, *61*, 2803–2817. [CrossRef]

5. McMillan, H.; Montanari, A.; Cudennec, C.; Savenije, H.; Kreibich, H.; Krüger, T.; Liu, J.; Meija, A.; van Loon, A.; Aksoy, H.; et al. Panta Rhei 2013–2015: Global perspectives on hydrology, society and change. *Hydrol. Sci. J.* **2016**, *61*, 1174–1191. [CrossRef]

6. Bong, C.H.J.; Lau, T.L.; Ghani, A.A.; Chan, N.W. Sediment deposit thickness and its effect on critical velocity for incipient motion. *Water Sci. Technol.* **2016**, *74*, 1876–1884. [CrossRef] [PubMed]

7. May, R.W.; Ackers, J.C.; Butler, D.; John, S. Development of design methodology for self-cleansing sewers. *Water Sci. Technol.* **1996**, *33*, 195–205. [CrossRef]

8. Butler, D.; May, R.; Ackers, J. Self-cleansing sewer design based on sediment transport principles. *J. Hydraul. Eng.* **2003**, *129*, 276–282. [CrossRef]

9. Safari, M.J.S.; Aksoy, H.; Unal, N.E.; Mohammadi, M. Non-deposition self-cleansing design criteria for drainage systems. *J. Hydro-Environ. Res.* **2017**, *14*, 76–84. [CrossRef]

10. Mayerle, R.; Nalluri, C.; Novak, P. Sediment transport in rigid bed conveyances. *J. Hydraul. Res.* **1991**, *29*, 475–495. [CrossRef]

11. May, R.W.P. *Sediment Transport in Pipes and Sewers with Deposited Beds*; Report SR 320; HR Wallingford: Wallingford, UK, 1993.

12. Ota, J.J.; Nalluri, C. Urban storm sewer design: Approach in consideration of sediments. *J. Hydraul. Eng.* **2003**, *129*, 291–297. [CrossRef]

13. Vongvisessomjai, N.; Tingsanchali, T.; Babel, M.S. Non-deposition design criteria for sewers with part-full flow. *Urban Water J.* **2010**, *7*, 61–77. [CrossRef]

14. Ota, J.J.; Perrusquia, G.S. Particle velocity and sediment transport at the limit of deposition in sewers. *Water Sci. Technol.* **2013**, *67*, 959–967. [CrossRef] [PubMed]

15. Loveless, J.H. Sediment Transport in Rigid Boundary Channels with Particular Reference to the Condition of Incipient Deposition. Ph.D. Thesis, University of London, London, UK, 1992.

16. Safari, M.J.S.; Aksoy, H.; Mohammadi, M. Incipient deposition of sediment in rigid boundary open channels. *Environ. Fluid Mech.* **2015**, *15*, 1053–1068. [CrossRef]
17. Safari, M.J.S.; Aksoy, H.; Mohammadi, M. Artificial neural network and regression models for flow velocity at sediment incipient deposition. *J. Hydrol.* **2016**, *541*, 1420–1429. [CrossRef]

18. Task Committee. Sediment transportation mechanics: Initiation of motion. *J. Hydraul. Div. ASCE* **1966**, *92(HY2)*, 291–314.

19. May, R.W.P. *Sediment Transport in Sewers*; Report No. IT 222; Hydraulic Research Station: Wallingford, UK, 1982.

20. Ackers, P. Sediment transport in sewers and the design implications. In *Proceedings of the International Conference on Planning, Construction, Maintenance, and Operation of Sewerage Systems, Reading, UK*, 12–14 September 1984; pp. 215–230.

21. Sanford, L.P.; Halka, J.P. Assessing the paradigm of mutually exclusive erosion and deposition of mud, with examples from upper Chesapeake Bay. *Mar. Geol.* **1993**, *114*, 37–57. [CrossRef]

22. Aksoy, H.; Safari, M.J.S. *Incipient Deposition of Sediment Particles in Rigid Boundary Channels*; Technical Report 113M062; Scientific and Technological Research Council of Turkey (TÜBİTAK): Istanbul, Turkey, 2014. (In Turkish)

23. Unal, N.E.; Aksoy, H.; Safari, M.J.S. *Self-Cleansing Drainage System Design by Incipient Motion and Incipient Deposition-Based Models*; Technical Report 114M283; Scientific and Technological Research Council of Turkey (TÜBİTAK): Istanbul, Turkey, 2016. (In Turkish)

24. Shields, A. *Application of Similarity Principles and Turbulence Research to Bed-Load Movement*; Translation from German; California Institute of Technology: Pasadena, CA, USA, 1936.

25. Yalin, M.S. *Mechanics of Sediment Transport*; Pergamon Press: Oxford, NY, USA, 1972.

26. Aksoy, H.; Safari, M.J.S.; Unal, N.E.; Mohammadi, M. Velocity-based analysis of incipient deposition in rigid boundary channels. *Water Sci. Technol.* **2017**, *76*, 2535–2543. [CrossRef] [PubMed]

27. Ippen, A.T.; Verma, R.P. *The Motion of Discrete Particles along the Bed of a Turbulent Stream: International Hydraulics Convention*; St. Anthony Falls Laboratory: Minneapolis, MN, USA, 1953; pp. 7–20.

28. Novak, P.; Nalluri, C. Sediment transport in smooth fixed bed channels. *J. Hydraul. Div. ASCE* **1975**, *101*, 1139–1154.

29. Butler, D.; May, R.W.P.; Ackers, J.C. Sediment transport in sewers Part 1: Background. *Proc. Inst. Civ. Eng.-Water Marit. Energy* **1996**, *118*, 103–112. [CrossRef]

30. Paphitis, D. Sediment movement under unidirectional flows: An assessment of empirical threshold curves. *Coas. Eng.* **2001**, *43*, 227–245. [CrossRef]

31. Fredsoe, J.; Deigaard, R. *Mechanics of Coastal Sediment Transport*; World Scientific: Singapore, 1992.

32. Novak, P.; Nalluri, C. Incipient motion of sediment particles over fixed beds. *J. Hydraul. Res.* **1984**, *22*, 181–197. [CrossRef]

33. Yang, C.T. Incipient motion and sediment transport. *J. Hydraul. Div. ASCE* **1973**, *99*, 1679–1704.

34. Vanoni, V.A. *Sedimentation Engineering*; American Society of Civil Engineers (ASCE): Reston, VA, USA, 2006.

<figure><img src="image_100.png" /></figure>

© 2018 by the author. 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/).
Article

# Riverbed Migrations in Western Taiwan under Climate Change

Yi-Chiung Chao<sup>1,●</sup>, Chi-Wen Chen<sup>1,2,*●</sup>, Hsin-Chi Li<sup>1,●</sup> and Yung-Ming Chen<sup>1</sup>

Received: 2 November 2018; Accepted: 8 November 2018; Published: 12 November 2018

<figure><img src="image_101.png" /></figure>

check for
updates

**Abstract:** In recent years, extreme weather phenomena have occurred worldwide, resulting in many catastrophic disasters. Under the impact of climate change, the frequency of extreme rainfall events in Taiwan will increase, according to a report on climate change in Taiwan. This study analyzed riverbed migrations, such as degradation and aggradation, caused by extreme rainfall events under climate change for the Choshui River, Taiwan. We used the CCHE1D model to simulate changes in flow discharge and riverbed caused by typhoon events for the base period (1979–2003) and the end of the 21st century (2075–2099) according to the climate change scenario of representative concentration pathways 8.5 (RCP8.5) and dynamical downscaling of rainfall data in Taiwan. According to the results on flow discharge, at the end of the 21st century, the average peak flow during extreme rainfall events will increase by 20% relative to the base period, but the time required to reach the peak will be 8 h shorter than that in the base period. In terms of the results of degradation and aggradation of the riverbed, at the end of the 21st century, the amount of aggradation will increase by 33% over that of the base period. In the future, upstream sediment will be blocked by the Chichi weir, increasing the severity of scouring downstream. In addition, due to the increased peak flow discharge in the future, the scouring of the pier may be more serious than it is currently. More detailed 2D or 3D hydrological models are necessary in future works, which could adequately address the erosive phenomena created by bridge piers. Our results indicate that not only will flood disasters occur within a shorter time duration, but the catchment will also face more severe degradation and aggradation in the future.

**Keywords:** aggradation; CCHE1D; climate change; degradation; dynamical downscaling; flow discharge; migration; riverbed; sediment

## 1. Introduction

Global climate change has led to changes in the spatial and temporal distribution of precipitation, evaporation, and runoff, which in turn contribute to environmental issues such as frequent droughts, increased flooding, waterlogging disasters, and aggravated soil erosion [1]. At the present time, one of the main problems in soil erosion research around the world is to assess the impact of climate change on sediment cycling [2]. It is expected that the increase in global precipitation will affect the extent, frequency, and magnitude of soil erosion and sediment redistribution [3,4]. This could lead to more severe hydrological cycles [5] and increased rainfall erosivity [6]. Nevertheless, contrasting impacts have recently been found in different parts of the world. Zhao et al. [7] observed that owing to anthropogenic causes, the sediment export from the Loess Plateau of China has increased significantly,

<footnote n="1">1. National Science and Technology Center for Disaster Reduction, No. 200, Sec. 3, Beixin Road, Xindian District, New Taipei City 23143, Taiwan; ycchao@ncdr.nat.gov.tw (Y.-C.C.); hsinchi@ncdr.nat.gov.tw (H.-C.L.); ymchen@ncdr.nat.gov.tw (Y.-M.C.)</footnote>

<footnote n="2">2. Center for Spatial Information Science, The University of Tokyo, 5-1-5 Kashiwanoha Kashiwa, Chiba 277-8568, Japan</footnote>

<footnote>* Correspondence: kevin4919@gmail.com; Tel.: +81-4-7136-4306</footnote>
but with the recent climate change, stream flow and sediment load have decreased enormously. Foster
et al. [8] analyzed changes in sediment transport for the Karoo uplands, South Africa, and found that
the sediment yield has generally increased over the past decades, which is related to factors such as
the increase in the frequency of high magnitude rainfall events.

Numerical models have been employed in engineering studies to predict the flow discharge and riverbed migrations along rivers. Several 1-D models have been proposed, and each model has its advantages and limitations [9]. The use of these models requires the proper simulation of hydraulic conditions [10]. The development of 1-D models is mainly focused on numerical aspects of the advection–dispersion equation [11–13] and exchange with dead zones [14–16]. Some 1-D models, such as the SIMCAT [17], QUAL2KW [18], and Multiphysics software COMSOL [19], are limited to steady flow conditions, while other models, such as the OTIS [20], CCHE1D [21], MIKE 11 [22], SD model [23], HEC-RAS [24], and ADISTS [25], can simulate unsteady flows and solute transport. In general, the validation testing focuses on theoretical situations or simplified river geometries in limited spatial and temporal scales [26].

Climate change is an important issue that has attracted growing attention in recent years. The report of the International Panel on Climate Change (IPCC) noted that rises in temperature and sea level are ongoing and that appropriate adaptations to reduce disaster risk are necessary [1]. In the current changing environment, both the probability of the occurrence of a strong typhoon and the rainfall intensity during a typhoon event will increase [27–29]. Many studies have explored the potential effects of different climate change scenarios on runoff and sediment yield in catchments [30–34]. In addition, there are also several studies focused on the impact of a changing climate on fluvial hydro-morphodynamics, as well as on the uncertainties related to future scenarios [35–41]. According to the Taiwan climate change report, the frequency of extreme rainfall events has also increased in the past 20 years in Taiwan [42]. Under the impact of climate change, extreme rainfall events will become more frequent, and riverbed migrations, such as degradation and aggradation, will become more serious, possibly affecting existing constructions and land use, and even human life and property. Therefore, the objectives of this study were to calibrate and validate a general one-dimensional channel network model for the Choshui River in western Taiwan and to assess changes in flow discharge and riverbed caused by typhoon events for the base period (1979–2003) and the end of the 21st century (2075–2099) according to the climate change scenario and dynamical downscaling of rainfall data in Taiwan.

## 2. Study Area

The Choshui River, located in western Taiwan, originates from the saddle between the main peak and east peak of Mount Hehuan. It flows through Nantou County, Chiayi County, Yunlin County, and Changhua County before reaching the Taiwan Strait. Having a mainstream about 187 km long, it is the longest river in Taiwan, and it has an average slope of 0.53%. The main tributaries of Choshui River are the Danda River, Jyunda River, Chenyoulan River, Tonpuze River, and Chingshui River. The Choshui River catchment covers an area of about 3157 km², which is the second largest catchment area in Taiwan. The topography slopes downward from east to west and the area of elevation >1000 m accounts for 64% of the catchment area. The slope decreases from east to west, and the area of slope >29% accounts for 62.5%, with an average slope of 54.5% (Figure 1).

The Choshui River catchment is characterized by friable Miocene and Pliocene sedimentary, metasedimentary, and low-grade metamorphic rocks and Quaternary alluvium deposits [43]. The soil cover in this area is mainly composed of colluvial soil, alluvial soil, and lithosol, with some yellow soil and laterite [44]. The Choshui River catchment can be divided into upstream and downstream areas at the Chunyun Bridge. The upstream area is well-covered with vegetation; the total area of natural forest and artificial forest is about 470 km², accounting for 92% of all land use. In the downstream area, agriculture occupies 1530 km², accounting for 85.5% of the land use, and the construction area is about 56 km², accounting for 3.14%.
The annual average rainfall in the Choshui River catchment is 2200 mm, but the spatial distribution of rainfall varies greatly. The annual rainfall in the upstream area can average 2500 mm, in contrast to 1100 mm in the downstream area and only 830 mm in the coastal area. The flow discharge of the Choshui River varies enormously between wet and dry seasons. About 80% of the total flow discharge is discharged during the wet season (May–October) [45]. According to [46], from the convergence of the Choshui River and the Chingshui River to the estuary, the flow discharges of the five-, 10-, 50-, and 100-year recurrence intervals are 0.65 × 10⁴, 1.11 × 10⁴, 1.44 × 10⁴, 2.27 × 10⁴, and 2.66 × 10⁴ m³/s, respectively.

<figure><img src="image_102.png" /><figcaption>Figure 1. Choshui River catchment of Taiwan.</figcaption></figure>

# 3. Data and Methods

## 3.1. CCHE1D Model

The National Center for Computational Hydroscience and Engineering of the University of Mississippi has developed a general one-dimensional channel network model that can simulate steady and unsteady flows and sedimentation processes in dendritic channel networks, including fractional sediment transport, riverbed aggradation and degradation, riverbed materials, bank erosion, and the resulting channel morphologic changes [21], which is called the CCHE1D model. The CCHE1D model solves the Saint–Venant equation and applies the implicit Preissmann four-point finite-difference and discrete-governing equation. A linear iterative method was developed in the model for the discrete-governing equation, which is solved by the double-sweep method. The one-dimensional hydraulic governing equations are as follows [47,48]:

Continuity equation:

$$ \frac{\partial A}{\partial t} + \frac{\partial Q}{\partial x} = q \qquad (1) $$
Momentum equation:

$$ \frac{\partial}{\partial t} \left( \frac{Q}{A} \right) + \frac{\partial}{\partial x} \left( \frac{\beta Q^2}{2A^2} \right) + g \frac{\partial h}{\partial x} + g (S_f - S_0) = 0 \quad (2) $$

where $x$ and $t$ are the spatial (m) and temporal (s) axes, $A$ is the flow area (m²), $Q$ is the flow discharge (m³/s), $q$ is the side discharge per unit channel length (m²/s), $h$ is the water surface elevation (m), $S_f$ and $S_0$ are the friction slope and riverbed slope (m/m), $\beta$ is the correction factor due to the non-uniformity of velocity distribution in the cross section, and $g$ is the gravitational acceleration (m/s²).

In addition, when the inertial force is small, the diffusion wave model can be used for the hydraulic routing. In general, when the Froude number is less than 0.5, the first two items on the left side of Equation (2) can be ignored. Then the momentum equation can be rewritten into the diffusion wave model, and the equation is as follows:

$$ \frac{\partial h}{\partial x} + S_f - S_0 = 0 \quad (3) $$

Equation (3) can still be used when the Froude number is larger than 0.5 if the simulation case is close to a uniform flow, because the first two items on the left side of Equation (2) will be sufficiently small to be ignored.

In terms of the sediment transport model, the CCHE1D model uses the non-equilibrium sediment transport concept to calculate the non-uniform sediment transport in the river. The equation for the non-equilibrium and non-uniform sediment transport control of total load is expressed as follows:

$$ \frac{\partial(AC_{tk})}{\partial t} + \frac{\partial Q_{tk}}{\partial x} + \frac{1}{L}(Q_{tk} - Q_{t+k}) = q_{lk} \quad (4) $$

where $C_{tk}$ is the section-averaged sediment concentration of size class $k$ (kg/m³), $Q_{tk}$ is the actual sediment transport rate (kg/s), $Q_{t+k}$ is the sediment transport capacity or the so-called equilibrium transport rate (kg/s), $L$ is the non-equilibrium adaptation length of sediment transport (m), and $q_{lk}$ is the side inflow or outflow sediment discharge from bank boundaries or tributary streams per unit channel length (kg/m/s).

The CCHE1D model offers the following four sediment discharge formulas for calculation: (1) SEDTRA module [49]; (2) Wu et al. formula [50]; (3) Modified Ackers and White's 1973 formula [51]; and (4) Modified Engelund and Hansen's 1967 formula [47]. In this study, we used the Wu et al. formula for our analyses because it is the most suitable formula for Taiwan, according to previous studies [52,53].

## 3.2. Parameter Setting for the Model

This study used the measured profile data of the Choshui River (2004, 2008, and 2012), provided by the 4th River Management Office of Taiwan, to build the model and set parameters. We used the measured profile data of 2004 as the initial topography for flood event simulations after 2004. Then the measured profile data of 2008 were used for calibrating the model, and those of 2012 were used for validation. We simulated the mainstream of the Choshui River for about 76 km from the Shuanglong Bridge (5 km upstream of the Yufeng Bridge) to the estuary. The tributary of the Choshui River, the Chingshui River, for about 13 km from the Longmen Bridge to the convergence with the mainstream was also simulated, and the results were used as the input data for the Choshui River.

This study used flow discharge during the wet seasons from 2004 to 2007 as the model calibration and continued to use flow discharge during the wet seasons from 2008 to 2011 as the validation. We used the change of time series in flow discharge during the period of the wet season for simulations. Flow discharge observed at the Yufeng Bridge station was used as the upstream boundary condition of the Choshui River, while flow discharge observed at the Tongtou station (2.5 km upstream of the
Longmen Bridge) was used as the upstream boundary condition of the Chingshui River. However, the Tongtou station was abandoned on 31 August 2009, so for points later than that, we used the flow discharge observed at the Longmen Bridge station (Figure 1). For the downstream, the open boundary method was applied for simulations.

Data of the grain size distribution of the riverbed of the Choshui River were measured in 2007, while those of the Chingshui River were measured in 2008 [46,54]. Manning's roughness coefficients for a riverbed of 0.029–0.044 were used, as recommended in government reports [46,54]. We applied the rating curves of flow discharge and sediment discharge to calculate the input of sediment, and also considered wash load in the calculations (Figure 2). Correlation coefficients ($r^2$) of rating curves for the Choshui River and the Chingshui River are 0.58 and 0.47 respectively. Although the $r^2$ is not very high, this study is mainly focused on riverbed migrations caused by heavy rainfall, which show a better correlation between flow discharge and sediment concentration when flow discharge is high (Figure 2). For the other general parameters for sediment transport simulation, this study set the suspended load adaptation coefficient ($\alpha$) as 1.0, bed load adaptation length ($L$) as 1000 m, and mixing layer as 2.0 m. According to [46], nine groups of representative particle sizes were selected, ranging from 0.001 to 1000 m, and a porosity of 0.4 was used. Yeh et al. [50] have made a sensitivity analysis on parameters of CCHE1D model in the Choshui River and found that the sensitivity of most parameters is not high. Only the thickness of mixing layer is sensitive to the degradation and aggradation of riverbed. This study followed the model parameters suggested by them for the parameter setting to get the most accurate results.

<figure><img src="image_103.png" /><figcaption>**Figure 2.** Rating curves of flow discharge and sediment concentration for (a) the Yufeng Bridge station (1994/2–2013/5), and (b) the Tongtou station (1994/6–2009/8) and the Longmen Bridge station (2009/10–2013/5).</figcaption></figure>

### 3.3. Assessing Degradation and Aggradation Conditions of Rivers under Climate Change

The resolution of the general circulation model (GCM) used for climate estimation and simulation is usually 100–200 km. Consequently, it is difficult to identify the structure of a typhoon with a radius of less than several hundred kilometers. Owing to the low spatial resolution of the model, it is impossible to reasonably simulate the extreme weather system with a physical model, and it is impossible to carry out extreme weather estimation under climate change [55,56]. Therefore, this study used a climate projection simulated from the high-resolution (approximately 20 km horizontal resolution) atmospheric general circulation model (AGCM) developed by the Meteorological Research Institute (MRI) under
the Japan Meteorological Agency (JMA) [57]. The low boundary drive by ensemble mean sea surface temperature (SST) of the climate change scenario of representative concentration pathways 8.5 (RCP8.5) in Coupled Model Intercomparison Project Phase 5 (CMIP5) [58] is called MRI-AGCM. This model has been validated for the simulation of typhoon generation [59,60]. However, the 20-km horizontal resolution still could not necessarily reflect the precipitation intensity caused by topographical effects in Taiwan. Therefore, the MRI-AGCM data were dynamically downscaled using the Weather Research and Forecasting Model (WRF), proposed by the US National Center for Atmospheric Research (NCAR), to obtain more realistic rainfall data for Taiwan with spatial and temporal resolutions of 5 km per hour.

Based on the downscaled rainfall data, there were 87 typhoons during the base period (1979–2003) and will be 43 typhoons in the last quarter of the 21st century. We selected eight typhoons in the base period and four typhoons at the end of the century, which accounted for the top 10% of the total rainfall for the Choshui River catchment among all typhoons in the two periods, to discuss the impact of climate change on riverbed migrations.

# 4. Results and Discussion

## 4.1. *Model Calibration and Validation*

This study used flow discharge during the wet seasons from 2004 to 2007 as the model calibration and continued to use flow discharge during the wet seasons from 2008 to 2011 as the validation. The total times for both periods of calibration and validation were 18,790 and 16,590 h, respectively. Figure 3 shows the comparison result between calibration and validation of the measured and simulated flow discharges at the Chunyun Bridge. We can find that the simulated peak flow discharge during typhoons is underestimated both in calibration and validation. Even so, trends of simulated and measured flow discharge are similar. The maximum percent errors of peak discharge (EQp) are -1.98% in the calibration and -9.50% in the validation. Overall, the coefficient of efficiency (CE) of the CCHE1D model in the calibration and validation are 0.85 and 0.83, respectively, indicating that the CCHE1D model can accurately simulate the rise and decline of peak discharge; therefore, we can use this model for simulating sediment transport under climate change scenarios.

<figure><img src="image_104.png" /><figcaption>**Figure 3.** Calibration and validation of the measured and simulated flow discharge at the Chunyun Bridge.</figcaption></figure>
In terms of degradation and aggradation in the river, this study used the measured profile data of 2004 as the initial topography for flood event simulations after 2004. Then the measured profile data of 2008 were used for calibrating the model and the measured profile data of 2012 were used for validation. Figure 4a shows the calibration result of the riverbed changes during 2004–2008 using the lowest point of each section along the river. The CCHE1D model can reasonably reproduce the process of degradation and aggradation in the river. Figure 4b shows the validation results of the riverbed changes for 2008–2012. It can be found that if this set of parameters is used, although the aggradation in downstream was underestimated, the riverbed changes in degradation and aggradation for most parts could be well validated. This represents that the CCHE1D model can truly grasp characteristics of degradation and aggradation in the Choshui River, so this model can also be used to simulate the degradation and aggradation in the river under climate change scenarios.

<figure><img src="image_105.png" /><figcaption>**Figure 4.** (a) Calibration result of the riverbed changes during 2004–2008 and (b) the validation result of the riverbed changes during 2008–2012 using the lowest point of each section along the river.</figcaption></figure>

## 4.2. Peak Flow Discharge during Typhoons under Climate Change

Based on the downscaled rainfall data, there were 87 typhoons during the base period (1979–2003) and will be 43 typhoons in the last quarter of the 21st century (2075–2099). Table 1 shows the simulation results of peak flow discharge of each of the top 20 typhoons during the two periods, accounting for the top 20 of the total rainfall amounts for the Choshui River catchment among all typhoons, at the Chunyun Bridge. According to the government report [44], the peak flow discharge of the 100-year recurrence interval ($Q_{100}$) in the Choshui River before the Chingshui River merges with it is 25,394 m³/s, which was used as the basis for the design of levees. We can find from the simulation results that the peak flow discharge during two typhoon events in the base period is greater than $Q_{100}$, while the peak flow discharge during three typhoon events at the end of the 21st century is greater than $Q_{100}$. This result indicates that, relative to the base period, the Choshui River is more likely to have a flood overflow at the end of the 21st century.

In addition, the peak flow discharge of the two-year recurrence interval ($Q_2$) in the Choshui River before the Chingshui River merge point is 5829 m³/s [44]. The peak flow discharge caused by two typhoon events in the base period is less than $Q_2$, and the peak flow discharge of six typhoon events at the end of the 21st century is less than $Q_2$. This result suggests that the flood disasters caused by
extreme rainfall events will be less frequent at the end of the 21st century than in the base period. We could assess from the above results that the sediment transport in the Choshui River was more frequent in the base period than it will be at the end of the 21st century; however, at the end of the 21st century, the severity of sediment migration will be more obvious than that in the base period.

<table><caption>**Table 1.** Simulation results of peak flow discharge of each of the top 20 typhoons during the two periods at the Chunyun Bridge.</caption><thead><tr><td>Event Ranking</td><td>Base Period (m³/s)</td><td>End of the 21st Century (m³/s)</td></tr></thead><tbody><tr><td>1</td><td>16,325</td><td>14,000</td></tr><tr><td>2</td><td>27,005 **</td><td>26,873 **</td></tr><tr><td>3</td><td>20,897</td><td>29,513 **</td></tr><tr><td>4</td><td>10,015</td><td>13,278</td></tr><tr><td>5</td><td>25,577 **</td><td>27,600 **</td></tr><tr><td>6</td><td>7422</td><td>15,706</td></tr><tr><td>7</td><td>23,515</td><td>23,062</td></tr><tr><td>8</td><td>8788</td><td>11,561</td></tr><tr><td>9</td><td>23,195</td><td>14,170</td></tr><tr><td>10</td><td>11,830</td><td>22,305</td></tr><tr><td>11</td><td>11,620</td><td>17,313</td></tr><tr><td>12</td><td>14,842</td><td>11,857</td></tr><tr><td>13</td><td>5071 *</td><td>9660</td></tr><tr><td>14</td><td>8931</td><td>8662</td></tr><tr><td>15</td><td>6519</td><td>5045 *</td></tr><tr><td>16</td><td>9048</td><td>2794 *</td></tr><tr><td>17</td><td>7128</td><td>2175 *</td></tr><tr><td>18</td><td>7977</td><td>3716 *</td></tr><tr><td>19</td><td>9437</td><td>3730 *</td></tr><tr><td>20</td><td>5341 *</td><td>1068 *</td></tr></tbody></table>

Here we selected the top 10% typhoons during the base period (1979–2003; $n = 8$) and during the end of the 21st century (2075–2099; $n = 4$), respectively, according to the total rainfall for the Choshui River catchment, for the following discussion on peak flow discharge along the whole reach of the Choshui River. The average peak flow discharges at the Yufeng Bridge in the base period and at the end of the 21st century are 17,443 and 20,916 m³/s, respectively. On average, it takes about 48 and 40 h to reach the peak values. We can notice that the average peak flow caused by extreme rainfall events at the end of the 21st century will increase by 20% from the base period, but the time required to reach the peak at the end of the 21st century will be 8 h shorter than that in the base period. It means that we will face greater flood events and need more urgent response time in the future. Figure 5a shows the simulation result of the most severe flood event at the end of the 21st century (peak flow discharge = 29,513 m³/s at the Chunyun Bridge). During the most severe flood event at the end of the 21st century, the locations where the overflow may occur are the estuary and the confluence of the Choshui River and Chingshui River (Figure 5b,c). Generally, in the case of normal energy transfer, energy transfer from upstream to downstream will be affected by topography and other factors, and energy will be lost, resulting in a lower water level downstream than upstream. The above water level is the result of quantitative simulations; however, the real situation is a dynamic flow. Therefore, the real downstream water level may be lower than the simulation result of the quantitative flow owing to energy losses, but it is still necessary to pay attention to flood control safety.

<footnote>* represents that the peak flow discharge is less than that of the two-year recurrence interval ($Q_2$ = 5829 m³/s).</footnote>

<footnote>** represents that the peak flow discharge is greater than that of the 100-year recurrence interval ($Q_{100}$ = 25,394 m³/s).</footnote>
<figure><img src="image_106.png" /><figcaption>Figure 5. Simulation result of the most severe flood event at the end of the 21st century.</figcaption></figure>

4.3. *Changes in Degradation and Aggradation of the Riverbed under Climate Change*

Changes in degradation and aggradation of the riverbed of the Choshui River under climate
change scenarios are shown in Figure 6. The Choshui River can be divided into three simulation
zones according to the artificial weir and the intersection of the main stream and tributary: (I) the
estuary to the Chunyun Bridge; (II) the Chunyun Bridge to the Chichi weir; (III) the Chichi weir to the
Shuanglong Bridge. It can be found from the results that zone I presents a slight aggradation from
the base period to the end of the 21st century; however, the problem of pier scouring still requires
attention. Zone II is located at lower reaches of the Chichi weir, so the riverbed suffers from strong
changes in degradation and aggradation. At the end of the 21st century, the degree of degradation
will be more serious than that in the base period. Because of the effect of the Chichi weir, zone III
has obvious sediment accumulation in the upper part of the weir, resulting in the aggradation of the
riverbed in the area. However, the location of the Yufeng Bridge shows clear degradation. Overall,
the simulation results show that the entire area of simulated reaches is aggraded and that the amount
of aggradation increases by about 33% at the end of the 21st century as compared with the base period.

Li et al. [61] pointed out that the landslide volume in the Choshui River catchment will increase
by 14% at the end of the 21st century, relative to that of the base period. This should also be reflected in
the results of degradation and aggradation of the riverbed. Landslides in the Choshui River catchment
will increase in the future, and landslide materials will accumulate in rivers. On the other hand,
the increase of peak flow discharge in the future will severely flush the surrounding structures and
downstream areas, which may cause crises such as damaged bridges and dam dumping.
<figure><img src="image_107.png" /><figcaption>**Figure 6.** Changes in erosion and aggradation of the riverbed in the Choshui River under climate change scenarios.</figcaption></figure>

Although degradation and aggradation of the riverbed are natural phenomena of a river, under the effect of climate change, the rainfall pattern will change to short duration-high intensity rainfall in the future [27–29]. This change will result in more serious degradation and aggradation of the riverbed in a catchment that occurs in the present situation. It could be imagined that not only will the response time during flood events be shortened, but also more serious riverbed migrations will occur in the future.

**5. Conclusions**

This study used the flow discharge during the wet seasons from 2004 to 2012 for calibration and validation of the CCHE1D model, a general one-dimensional channel network model that can simulate steady and unsteady flows and sedimentation processes in dendritic channel networks, along the Choshui River in western Taiwan. The CCHE1D model was used for simulating changes in flow discharge and riverbed caused by typhoon events for the base period (1979–2003) and the end of the 21st century (2075–2099). According to the climate change scenario of RCP8.5 and dynamically downscaled rainfall data in Taiwan, from the base period to the end of the 21st century, average peak flow caused by extreme rainfall events will increase by about 20%, but the time required to reach the peak at the end of the 21st century will be 8 h shorter than that in the base period. For the same time period, the amount of aggradation along the Choshui River will increase by about 33% from the base period to the end of the 21st century. Upper and lower parts of the Chichi weir will show large changes in aggradation and degradation, respectively, and the pier may suffer from serious scouring owing to the increase in peak flow discharge in the future. Our results could provide useful information to related authorities for the decision making on plans of public construction and land use in a changing climate. In future works, it is essential to use 2D or 3D hydrological models to adequately address problems such as the erosive phenomena created by bridge piers.

**Author Contributions:** Y.-C.C. collected the data, performed the statistical and simulation analysis, and drafted the manuscript. C.-W.C. conceived of the study, participated in its design and coordination, and drafted the manuscript. H.-C.L. and Y.-M.C. participated in the design of the study and provided methods. All authors read and approved the final manuscript.

**Funding:** This research was supported by the Taiwan Climate Change Projection and Information Platform (TCCIP; MOST106-2621-M-865-001).

**Acknowledgments:** We would like to thank the theme C of the Program for Risk Information on Climate Change of Japan (SOUSEI-C) for providing the MRI-AGCM data.

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

1. IPCC. *Climate Change 2007: The Physical Science Basis*. Summary for Policymakers; Contribution of working group I to the fourth assessment report; The Intergovernmental Panel on Climate Change: Geneva, Switzerland, 2007.

2. Mullan, D.; Favis-Mortlock, D.; Fealy, R. Addressing key limitations associated with modelling soil erosion under the impacts of future climate change. *Agric. Meteorol.* **2012**, *156*, 18–30. [CrossRef]

3. Peizhen, Z.; Molnar, P.; Downs, W.R. Increased sedimentation rates and grain sizes 2–4 Myr ago due to the influence of climate change on erosion rates. *Nature* **2001**, *410*, 891–897. [CrossRef] [PubMed]

4. Pruski, F.F.; Nearing, M.A. Climate-induced changes in erosion during the 21st century for eight US locations. *Water Resour. Res.* **2002**, *38*, [CrossRef]

5. Nearing, M.A.; Jetten, V.; Baffaut, C.; Cerdan, O.; Couturier, A.; Hernandez, M.; Le Bissonnais, Y.; Nichols, M.H.; Nunes, J.P.; Renschler, C.S.; et al. Modeling response of soil erosion and runoff to changes in precipitation and cover. *Catena* **2005**, *61*, 131–154. [CrossRef]

6. Nearing, M.A.; Pruski, F.F.; O'Neill, M.R. Expected climate change impacts on soil erosion rates: A review. *J. Soil Water Conserv.* **USA** **2004**, *59*, 43–50.

7. Zhao, G.; Mu, X.; Wen, Z.; Wang, F.; Gao, P. Soil erosion, conservation and eco-environment changes in the Loess Plateau of China. *Land Degrad. Dev.* **2013**, *24*, 499–510. [CrossRef]

8. Foster, I.D.L.; Rowntree, K.M.; Boardman, J.; Mighall, T.M. Changing sediment yield and sediment dynamics in the Karoo uplands, South Africa; post-European impacts. *Land Degrad. Dev.* **2012**, *23*, 508–522. [CrossRef]

9. Cox, B.A. A review of currently available in-stream water-quality models and their applicability for simulating dissolved oxygen in lowland rivers. *Sci. Total. Environ.* **2003**, *314*, 335–377. [CrossRef]

10. Leibundgut, C.; Maloszewski, P.; Kuülls, C. *Tracers in Hydrology*; Wiley-Blackwell: Hoboken, NJ, USA, 2009.

11. Russell, T.F.; Celia, M.A. An overview of research on Eulerian–Lagrangian Localised Adjoint Methods (ELLAM). *Adv. Water Resour.* **2002**, *25*, 1215–1231. [CrossRef]

12. Rubio, A.D.; Zalts, A.; El Hasi, C.D. Numerical solution of the advection–reaction–diffusion equation at different scales. *Environ. Model. Softw.* **2008**, *23*, 90–95. [CrossRef]

13. Shen, C.; Phanikumar, M.S. An efficient space-fractional dispersion approximation for stream solute transport modeling. *Adv. Water Resour.* **2009**, *32*, 1482–1494. [CrossRef]

14. Bencala, K.E.; Walters, R.A. Simulation of solute transport in a mountain pool-and-riffle stream: A transient storage model. *Water Resour. Res.* **1983**, *19*, 718–724. [CrossRef]

15. Wörman, A.; Packman, A.I.; Johansson, H.; Jonsson, K. Effect of flow-induced exchange in hyporheic zones on longitudinal transport of solutes in streams and rivers. *Water Resour. Res.* **2002**, *38*, 61. [CrossRef]

16. Anderson, E.J.; Phanikumar, M.S. Surface storage dynamics in large rivers: Comparing three-dimensional particle transport, one-dimensional fractional derivative, and multirate transient storage models. *Water Resour. Res.* **2011**, *47*, [CrossRef]

17. UK Environment Agency. *SIMCAT 7.6: A Guide and Reference for Users*; UK Environment Agency: Bristol, UK, 2001.

18. Pelletier, G.J.; Chapra, C.S.; Tao, H. QUAL2Kw, A framework for modeling water quality in streams and rivers using a genetic algorithm for calibration. *Environ. Model. Softw.* **2006**, *21*, 419–425. [CrossRef]

19. Ani, E.C.; Wallis, S.; Kraslawski, A.; Agachi, P.S. Development, calibration and evaluation of two mathematical models for pollutant transport in a small river. *Environ. Model. Softw.* **2009**, *24*, 1139–1152. [CrossRef]

20. Runkel, R.L. *One-Dimensional Transport with Inflow and Storage (OTIS): A Solute Transport Model for Streams and Rivers*; Water-Resources Investigations Report; United States Geological Survey: Reston, VA, USA, 1998; 80p.

21. Vieira, D.A.N. Integrated Modeling of Watershed and Channel Processes. Ph.D. Thesis, University of Mississippi, Oxford, MS, USA, 2004.

22. Danish Hydraulic Institute (DHI). *MIKE 11-A Modelling System for Rivers and Channels*; Danish Hydraulic Institute: Horsholm, Denmark, 2007.

23. Deng, Z.Q.; Jung, H.S. Scaling dispersion model for pollutant transport in rivers. *Environ. Model. Softw.* **2009**, *24*, 627–631. [CrossRef]

24. US Army Corps of Engineers (USACE). *HEC-RAS River Analysis System*; User's Manual—Version 4.1; USACE: Washington, DC, USA, 2010.
25. Launay, M.; Le Coz, J.; Camenen, B.; Walter, C.; Angot, H.; Dramaisa, G.; Faure, J.F.; Coquery, M. Calibrating pollutant dispersion in 1-D hydraulic models of river networks. *J. Hydro-Environ. Res.* **2015**, *9*, 120–132.
[CrossRef]

26. Zerihun, D.; Furman, A.; Warrick, A.W.; Sanchez, C.A. Coupled surface–subsurface solute transport model for irrigation borders and basins. *J. Irrig. Drain. Eng.* **2005**, *131*, 396–406. [CrossRef]

27. Emanuel, K. Increasing destructiveness of tropical cyclones over the past 30 years. *Nature* **2005**, *436*, 686–688.
[CrossRef] [PubMed]

28. Webster, P.J.; Holland, G.J.; Curry, J.A.; Chang, H.R. Changes in tropical cyclone number, duration, and intensity in a warming environment. *Science* **2005**, *309*, 1844–1846. [CrossRef] [PubMed]

29. Tu, J.Y.; Chou, C.; Chu, P.S. The abrupt shift of typhoon activity in the vicinity of Taiwan and its association with western North Pacific–East Asian climate change. *J. Clim.* **2009**, *22*, 3617–3628. [CrossRef]

30. Bussi, G.; Francés, F.; Horel, E.; López-Tarazón, J.A.; Batalla, R.J. Modelling the impact of climate change on sediment yield in a highly erodible Mediterranean catchment. *J. Soils Sediments* **2014**, *14*, 1921–1937.
[CrossRef]

31. Lee, S.Y.; Hamlet, A.F.; Grossman, E.E. Impacts of climate change on regulated streamflow, hydrologic extremes, hydropower production, and sediment discharge in the Skagit river basin. *Northwest Sci.* **2016**, *90*, 23–43.
[CrossRef]

32. Laura, M.; Tartari, G.; Salerno, F.; Valsecchi, L.; Bravi, C.; Lorenzi, E.; Genoni, P.; Guzzella, L. Climate change impacts on sediment quality of Subalpine reservoirs: Implications on management. *Water* **2017**, *9*, 680.
[CrossRef]

33. Ren, Z.; Feng, Z.; Li, P.; Wang, D.; Cheng, S.; Gong, J. Response of runoff and sediment yield from climate change in the Yanhe watershed, China. *J. Coast. Res.* **2017**, *80*, 30–35. [CrossRef]

34. Dahl, T.A.; Kendall, A.D.; Hyndman, D.W. Impacts of projected climate change on sediment yield and dredging costs. *Hydrol. Process.* **2018**, *32*, 1223–1234. [CrossRef]

35. Ashmore, P.; Church, M. *The Impact of Climate Change on Rivers and River Processes in Canada*; Geological Survey of Canada: Ottawa, ON, Canada, 2001; Volume 555.

36. Bronstert, A. Floods and climate change: Interactions and impacts. *Risk Anal. Int. J.* **2003**, *23*, 545–557.
[CrossRef]

37. Pfister, L.; Kwadijk, J.; Musy, A.; Bronstert, A.; Hoffmann, L. Climate change, land use change and runoff prediction in the Rhine–Meuse basins. *River Res. Appl.* **2004**, *20*, 229–241. [CrossRef]

38. Döll, P.; Zhang, J. Impact of climate change on freshwater ecosystems: A global-scale analysis of ecologically relevant river flow alterations. *Hydrol. Earth Syst. Sci.* **2010**, *14*, 783–799. [CrossRef]

39. Guerrero, M.; Nones, M.; Saurral, R.; Montroull, N.; Szupiany, R.N. Parana River morphodynamics in the context of climate change. *Int. J. River Basin Manag.* **2013**, *11*, 423–437. [CrossRef]

40. Piniewski, M.; Teresa, H.K.; Romanowicz, R.; Osuch, M.; Szcześniak, M.; Kardel, I.; Okruszko, T.; Mezghani, A.; Kundzewicz, Z.W. What can we learn from the projections of changes of flow patterns? Results from Polish case studies. *Acta Geophys.* **2017**, *65*, 809–827. [CrossRef]

41. Garijo, C.; Mediero, L. Influence of climate change on flood magnitude and seasonality in the Arga River catchment in Spain. *Acta Geophys.* **2018**, *66*, 769–790. [CrossRef]

42. Hsu, H.H.; Chou, C.; Wu, Y.C.; Lu, M.M.; Chen, C.T.; Chen, Y.M. *Climate Change in Taiwan: Scientific Report 2011* (Summary); National Science Council: Taipei, Taiwan, 2011.

43. Ho, C.S. *An Introduction to the Geology of Taiwan: Explanatory Text of the Geologic Map of Taiwan*; Central Geological Survey, Ministry of Economic Affairs: Taipei, Taiwan, 1988.

44. Water Resources Planning Institute (WRPI). *Review Report on Jhuo-Shuei River Basin Integrated Regulation Plan*; Water Resources Agency, Ministry of Economic Affairs: Taichung, Taiwan, 2011.

45. Water Resources Planning Institute (WRPI). *A Study of Flood Control and Sediment Management Due to Climate Change of Jhuoshuei River*; Water Resources Agency, Ministry of Economic Affairs: Taichung, Taiwan, 2013.

46. Water Resources Planning Institute (WRPI). *Review Report on River Planning of Jhuo-Shuei River*; Water Resources Agency, Ministry of Economic Affairs: Taichung, Taiwan, 2007.

47. Wu, W.; Vieira, D.A. *One-Dimensional Channel Network Model CCHE1D Version 3.0—Technical Manual*; The University of Mississippi, National Center for Computational Hydroscience and Engineering: Oxford, MS, USA, 2002.
48. Wu, W.; Vieira, D.A.; Wang, S.S.Y. One-Dimensional Numerical Model for Nonuniform Sediment Transport under Unsteady Flows in Channel Networks. *J. Hydraul. Eng.* **2004**, *130*, 914–923. [CrossRef]

49. Garbrecht, J.; Kuhnle, R.; Alonso, C. A Sediment Transport Capacity Formulation for Application to Large Channel Networks. *J. Soil Water Conserv.* **1995**, *50*, 527–529.

50. Wu, W.; Wang, S.S.Y.; Jia, Y. Nonuniform sediment transport in alluvial rivers. *J. Hydraul. Res.* **2000**, *38*, 427–434. [CrossRef]

51. Proffitt, G.T.; Sutherland, A.J. Transport of non-uniform sediments. *J. Hydraul. Res.* **1983**, *21*, 33–43. [CrossRef]

52. Yeh, K.C.; Wang, S.Y.; Wu, W.M.; Lien, H.C.; Liao, C.T. *Implement and Application Study of NCCHE's River Migration Models*; Water Resources Planning Institute, WRA: Taichung, Taiwan, 2007.

53. Liao, C.T.; Yeh, K.C.; Huang, M.W. Development and application of 2-D mobile-bed model with bedrock river evolution mechanism. *J. Hydro-Environ. Res.* **2014**, *8*, 210–222. [CrossRef]

54. Water Resources Planning Institute (WRPI). *Report on River Planning of Cing-Shuei River, Tributary of Jhuo-Shuei River*, Water Resources Agency, Ministry of Economic Affairs: Taichung, Taiwan, 2008.

55. Roble, R.G.; Ridley, E.C. A thermosphere-ionosphere-mesosphere-electrodynamics general circulation model (TIME-GCM): Equinox solar cycle minimum simulations (30–500 km). *Geophys. Res. Lett.* **1994**, *21*, 417–420. [CrossRef]

56. Johns, T.C.; Carnell, R.E.; Crossley, J.F.; Gregory, J.M.; Mitchell, J.F.; Senior, C.A.; Tett, S.F.B.; Wood, R.A. The second Hadley Centre coupled ocean-atmosphere GCM: Model description, spinup and validation. *Clim. Dynam.* **1997**, *13*, 103–134. [CrossRef]

57. Mizuta, R.; Yoshimura, H.; Murakami, H.; Matsueda, M.; Endo, H.; Ose, T.; Kamiguchi, K.; Hosaka, M.; Sugi, M.; Yukimoto, S.; et al. Climate simulations using MRI-AGCM3.2 with 20-km grid. *J. Meteor. Soc. Jpn.* **2012**, *90A*, 233–258. [CrossRef]

58. Mizuta, R.; Arakawa, O.; Ose, T.; Kusunoki, S.; Endo, H.; Kitoh, A. Classification of CMIP5 future climate responses by the tropical sea surface temperature changes. *SOLA* **2014**, *10*, 167–171. [CrossRef]

59. Murakami, H.; Wang, Y.; Yoshimura, H.; Mizuta, R.; Sugi, M.; Shindo, E.; Adachi, Y.; Yukimoto, S.; Hosaka, M.; Kusunoki, S.; et al. Future changes in tropical cyclone activity projected by the new high-resolution MRI-AGCM. *J. Clim.* **2012**, *25*, 3237–3260. [CrossRef]

60. Murakami, H.; Hsu, P.C.; Arakawa, O.; Li, T. Influence of model biases on projected future changes in tropical cyclone frequency of occurrence. *J. Clim.* **2014**, *27*, 2159–2181. [CrossRef]

61. Li, H.C.; Shih, H.J.; Wu, T.; Chao, Y.C.; Chen, W.B.; Cheng, C.T.; Chen, D.R.; Jang, J.H.; Chen, Y.M. *The Loss Assessment for Cho-Shui River Basin in the End of Century under Climate Change*; NCDR 104-T09; National Science and Technology Center for Disaster Reduction: Taipei, Taiwan, 2015.

<figure><img src="image_108.png" /></figure>

© 2018 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/).
<figure><img src="image_109.png" /></figure>

Article

# A Shear Reynolds Number-Based Classification Method of the Nonuniform Bed Load Transport

Gergely T. Török<sup>1,2,*</sup>①, János Józsa<sup>1,2</sup> and Sándor Baranya<sup>2</sup>

1 MTA-BME Water Management Research Group, Hungarian Academy of Science—Budapest University of Technology and Economics, Műegyetem rakpart 3, H-1111 Budapest, Hungary; jozsa.janos@epito.bme.hu

2 Department of Hydraulic and Water Resources Engineering, Faculty of Civil Engineering, Budapest University of Technology and Economics, Műegyetem rakpart 3, H-1111 Budapest, Hungary; baranya.sandor@epito.bme.hu

* Correspondence: torok.gergely@epito.bme.hu; Tel.: +36-1-463-2248

Received: 30 November 2018; Accepted: 24 December 2018; Published: 3 January 2019

<figure><img src="image_110.png" /></figure>

check for
updates

**Abstract:** The aim of this study is to introduce a novel method which can separate sand- or gravel-dominated bed load transport in rivers with mixed-size bed material. When dealing with large rivers with complex hydrodynamics and morphodynamics, the bed load transport modes can indicate strong variation even locally, which requires a suitable approach to estimate the locally unique behavior of the sediment transport. However, the literature offers only few studies regarding this issue, and they are concerned with uniform bed load. In order to partly fill this gap, we suggest here a decision criteria which utilizes the shear Reynolds number. The method was verified with data from field and laboratory measurements, both performed at nonuniform bed material compositions. The comparative assessment of the results show that the shear Reynolds number-based method operates more reliably than the Shields-Parker diagram and it is expected to predict the sand or gravel transport domination with a <5% uncertainty. The results contribute to the improvement of numerical sediment transport modeling as well as to the field implementation of bed load transport measurements.

**Keywords:** bed load transport; shear Reynolds number; mixed-size bed material; complex morphodynamics

## 1. Introduction

In terms of typical geomorphological features, rivers can be divided into three main section types: upper, middle and lower course [1]. One of the decisive differences among these river features is the erosion capacity of the flow, which constantly decreases along the river. In fact, it can be stated that the erosion prevails along the upper course rivers, yielding coarser gravel bed material and significant bed load transport. As to the middle course, all the three characteristic sediment transport processes, such as the local erosion of the bed, the transport of coarse and fine particles together with river bed aggradation take place. Finally, the lower course type rivers can be generally characterized with the deposition of fine sediments transported from upstream. Due to the spatially and temporally varying bed load transport capacity of the flow and the changing river planforms (straight, meandering, or braided river channel patterns [2]), the morphological features can also show significant variability.

Many classification methods can be found in the literature to describe the morphological properties of river sections (e.g., [3–7]). The main goal of these tools is to typify and predict the dominant morphological properties and processes, such as defining the predominant sediment transport mode (bed load, or suspended load) [8], the channel pattern type [1,9], the specific bed material grain size [7], the erosion capacity [10], or the bed armor measure [6]. According to these methods, the morphological
properties and processes can be estimated as the function of a few hydro-morphological variables, e.g., the water depth, the channel width, the mean flow velocity, the mean grain size of the bed material, or the longitudinal bed slope.

For instance, there are classification methods, suggesting that the dominant grain size of the bed load (sand or gravel) can be obtained based on the *d*₅₀ grain size of the bed material [8]. In most cases, such a classification is straightforward due to the fact that a dominant bed material fraction can be chosen, considering, e.g., sand-bed (*d*₅₀ < 2 mm) or gravel-bed (*d*₅₀ > 2 mm) streams [8]. However, there are several situations, when such a clear distinction cannot be made as both sand- and gravel-dominated zones appear even in shorter river reaches. Moreover, the locally dominant bed load fraction can alter according to the current flow conditions. So, only the motion of finer particles of a nonuniform bed material can occur during low-water regime. In turn, the higher flood waves increase the sediment transport capacity of the flow, resulting in probably gravel-dominated bed load. That is, because the bed load is mainly fed by the bed itself [11], very different grain sizes can characterize the bed load transport within river sections with mixed-size bed material. These sections are typical to middle-course rivers, at the transition sections in terms of bed material.

## Case Study

Such a river reach is the upper Hungarian Danube reach between rkm 1798 and rkm 1795 (Figure 1). At this section, the middle course resulted in a meandering pattern type, with nonuniform bed material (*d*₅₀ ranges between 0.32 and 70.5 mm) [9,12]. Furthermore, conventional river regulation measures (e.g., groin fields) were installed, which further enhanced the diversity of the morphological properties and processes [13–17]. The following parameters characterize the river: the main channel width at mean water regime ranges between 150 m and 350 m; the average water surface gradient is around 0.0002–0.00025 [17], whereas the characteristic flow discharges are Qₘ = 2000 m³/s (mean flow), Qₛf = 4300–4500 m³/s (range of bankfull discharge), Q₁₀₀ = 10,400 m³/s (100 year flood event) [18] and Q₂₀₀ = ~10,800 m³/s (200 year flood event) [19]. Here, the dominant fraction in the bed material, and hence in the bed load, shows a strongly varying spatial distribution. That is, bed armoring and gravel bed load takes place in the main channel, while clear sand transport can be observed in the near-bank zones and in the groin fields. Finally, the accumulation and erosion of the gravel–sand mixture is detectable at some places, resulting in the formation of side bars (Figure 1).

<figure><img src="image_111.png" /><figcaption>Figure 1. The sketch of the investigated section of the Danube River.</figcaption></figure>

The herein-presented study focuses on this issue: how the dominant fraction of the bed load transport can be determined, when both the river hydrodynamics and the morphological features show strong heterogeneity within the studied reach. When local or reach-scale morphodynamic analysis is