Next Article in Journal
Can We Calibrate a Daily Time-Step Hydrological Model Using Monthly Time-Step Discharge Data?
Next Article in Special Issue
Evolution Pattern of Tailings Flow from Dam Failure and the Buffering Effect of Debris Blocking Dams
Previous Article in Journal
Application of Poly-γ-Glutamic Acid Flocculant to Flocculation–Sedimentation Treatment of Ultrafine Cement Suspension
Previous Article in Special Issue
Effect of Sediment Load Boundary Conditions in Predicting Sediment Delta of Tarbela Reservoir in Pakistan
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Assessment of Sediment Transport Functions with the Modified SWAT-Twn Model for a Taiwanese Small Mountainous Watershed

Department of Civil and Disaster Prevention Engineering, National United University, Miaoli City 36063, Taiwan
*
Author to whom correspondence should be addressed.
Water 2019, 11(9), 1749; https://doi.org/10.3390/w11091749
Submission received: 4 July 2019 / Revised: 14 August 2019 / Accepted: 19 August 2019 / Published: 22 August 2019
(This article belongs to the Special Issue Modeling of Soil Erosion and Sediment Transport)

Abstract

:
In Taiwan, the steep landscape and highly vulnerable geology make it difficult to predict soil erosion and sediment transportation via variable transport conditions. In this study, we integrated the Taiwan universal soil loss equation (TUSLE) and slope stability conditions in the soil and water assessment tool (SWAT) as the SWAT-Twn model to improve sediment simulation and assess the sediment transport functions in the Chenyulan watershed, a small mountainous catchment. The results showed that the simulation of streamflow was satisfactory for calibration and validation. Before model calibration and validation for sediment, SWAT-Twn with default sediment transport method performed better in sediment simulation than the official SWAT model (version 664). The SWAT-Twn model coupled with the simplified Bagnold equation could estimate sediment export more accurately and significantly reduce the overestimated sediment yield by 65.7%, especially in highly steep areas. Furthermore, five different sediment transport methods (simplified Bagnold equation with/without routing by particle size, Kodoatie equation, Molinas and Wu equation, and Yang sand and gravel equation) were evaluated. It is suggested that modelers who conduct sediment studies in the mountainous watersheds with extreme rainfall conditions should adjust the modified universal soil loss equation (MUSLE) factors and carefully evaluate the sediment transportation equations in SWAT.

1. Introduction

The amount and occurrence time of sediment export from the watersheds are mainly caused by different types of soil erosion, such as channel bed and bank erosion, overland surface erosion, and glacier erosion [1]. As soil erosion causes shear stress by rainfall, surface runoff brings most sediment yields from the overland eroded soil and then is usually transported to the downstream [2]. Moreover, extreme natural disturbances (i.e., typhoons, earthquakes, landslides, and floods) also play as the trigger for abnormal sediment yields in some regions [3], and further change the characteristics of sediment yields and sediment transports in watersheds [4]. Many studies have indicated that climate change has caused higher rainfall intensity and annual precipitation, leading to increases in sediment yields and soil erosion rates [5].
In order to comprehensively quantify the impacts of rainfall, soil erodibility, land use cover, topography, and support practice to sediment yields, the universal soil loss equation (USLE) [6] and the modified universal soil loss equation (MUSLE) [7] have been developed. These two empirical equations have been used to estimate the soil loss in watersheds worldwide [8]. Some other physically-based erosion models were developed in the past decades, such as those used in ANSWERS (Areal Nonpoint Source Watershed Environmental Response Simulation) [9], EPIC (Environmental Policy Integrated Climate) model [10], WEPP (Water Erosion Prediction Project) [11], and EUROSEM (European Soil Erosion Model) [12]. The main difference between MUSLE and USLE is that USLE uses rainfall as an indicator of erosive energy, while MULSE uses the amount of runoff to simulate erosion and sediment yield. Sadeghi et al. [8] reviewed 49 papers of MUSLE application worldwide, and presented that most of MUSLE studies were to estimate the sediment yield on a storm basis (73.91%), monthly basis (2.17%), and annual basis (17.39%), while other applications of MUSLE (6.53%) were studied in soil erosion in storm-wise scale, pollutant estimation, and return periods of annual sediment yield. Therefore, MUSLE can better simulate the sediment yields during single storms with a low level of estimation error [13]. Since USLE and MUSLE were developed based on experiment watersheds in the U.S., [8] indicated that MUSLE would have huge errors without the calibration but can present reliable results for sediment yield on a storm basis after calibration, especially when it is applied under appropriate conditions (i.e., rangeland watersheds, similar climatic conditions as the USA) similar to where the original model was developed.
The process of erosion is usually described in three stages, which are detachment, transport and deposition; and four main types, which are sheet, rill, gully and in-stream erosion [14]. Sheet erosion and rill erosion are often classified as overland flow erosion caused by raindrop or overland flow [2,15]. Both types of erosions are usually considered together in erosion modelling. Gully erosion is the removal of soil along the channels of concentrated flow and is controlled by the thresholds related to slope and catchment area [2]. In-stream erosion involves the direct removal of sediment from stream beds and stream banks, especially during high flow periods huge amount of sediment transported through the stream network originates from the stream channel [14]. Erosion/sediment transport models can be categorized into three types: empirical model, conceptual model, and physics-based model [14]. The empirical models, the simplest ones, are based on the statistical observations of experiment areas in response to the characteristics of sediment transport [16], while conceptual models usually represent the catchment as a series of internal storages with parameters values determined through calibration against observed data [17]. Physics-based models integrate some fundamental physical equations, such as the equations of conservation of mass and momentum for flow, and those equations for sediment [14]. Compared to physics-based models and more complex conceptual models, which usually lack of sufficient spatially distributed input data, empirical and conceptual models are suggested to be combined for presenting the event responsiveness and sensitivity to climate variability [14].
The soil and water assessment Tool (SWAT) is a semi-distributed hydrological watershed model, which can simulate water balance, plant growth, and transport of sediment, nutrients, heavy metals and pesticides. The SWAT model has been used worldwide for simulating the impact of climate change and land use change on streamflow, sediment and nutrients exports, and best management practices (BMPs) on watershed responses in different countries [18,19,20]. Generally, the SWAT model performed well for simulating streamflow, sediment and nutrient exports at various spatio-temporal scales. MUSLE plays an important role of the simulation of continuous sediment loads [21]. Sediment transport is often called total sediment load, which is the sum of bed load and suspended load. The sediment transport is generally modeled through the sediment flow caused by overland flow and the channel erosion [22]. Some studies indicated unsatisfactory sediment simulation of the SWAT model. Addis et al. [23] used the SWAT model to simulate discharge and sediment transport at a small mountainous catchment in Ethiopian plateau, and showed that SWAT performed well for discharge but overestimated daily sediment transport with NSE = 0.07 and −1.76 for calibration and validation, respectively, mainly due to insufficient observed data and the heavy rainfall during simulation period. Bressiani et al. [24] indicated that the SWAT model was suitable for sediment simulation in most areas in Brazil and suggested the change in transmission of sediment should be reflected in the current SWAT version.
As the sediment transport and sediment concentration can affect the nutrient and turbidity in water, the estimation of sediment transport or sediment concentration need to be more reasonable. Arnold et al. [25] advised that SWAT users need to calibrate the discharge and sediment transport or concentration sequentially before simulating the nutrients. However, extreme rainfalls in Taiwan result in serious debris flows and landslides, making it more difficult to simulate the sediment transport in the river. As the runoff and sediment hazards in Taiwan are increasing, Lee et al. [26] indicated that the runoff and sediment yield would increase with the storm events and become more frequently occurred, especially in the small mountainous catchments. Chiu et al. [27] simulated the sediment yield in the Shihmen reservoir in Taiwan, and revealed that natural distributions (i.e., typhoon, storm, and earthquake) have become an important potential source of sediment yield. Therefore, Chang et al. [28] suggested that landslide should be considered in the model to simulate more accurate and reasonable sediment exports. In order to more accurately simulate the sediment transport and sediment yields in a small mountainous watershed in Taiwan, we aim to: (1) integrate Taiwan universal soil loss equation (TUSLE) and the landslide area-volume estimation equation into the SWAT model (version 664) as SWAT-Twn model; (2) examine five sediment transport methods in the SWAT model; (3) provide the calibration and validation experience of sediment transport simulation for the SWAT users.

2. Materials and Methods

2.1. Study Area and Data

The Chenyulan watershed, located in central Taiwan, has an area of 448 km2 with the elevation ranging from 292 m to 3893 m (Figure 1). The Chenyulan watershed has steep terrain that is 49.74% of the total area with slopes steeper than 60%. Typhoons averagely hit Taiwan for 3–4 times in a year, leading large amounts of flow and sediment in the watersheds. Especially in 1996 and 2009, Typhoon Herb and Typhoon Morakot have brought more than 2000 mm of accumulated rainfall in two days, resulting in 3370 m3/s and 1860 m3/s of peak discharge and daily average streamflow, respectively. Moreover, extremely high sediment concentration of 98,499 mg/L was observed at Nemoupu station in the Chenyulan watershed when Typhoon Morakot occurred in 2009. During the study period (2004–2015), several typhoons that influenced the watershed have been recorded (Table 1). Moreover, in 1999, the Chi-Chi earthquake of 7.3 Richter magnitude scale occurred in central Taiwan (Figure 1) and has seriously influenced the watershed landscape to become more fragile and sensitive [29,30].
The environmental data required for the SWAT model include: the digital elevation model (DEM), weather, observed streamflow and sediment concentration, land use and soil data. Moreover, the model simulation can be more precisely to represent the characteristics of the watershed by adding detailed information, such as point source pollutants, reservoirs, agriculture management scenarios [32,33,34,35,36,37]. The DEM data was available from Taiwan Geospatial One Stop (TGOS) [38]. We collected the data (i.e., precipitation, temperature, solar radiation, relative humidity, wind speed) of six weather stations in the watershed during 2003—2015 from the Taiwan Data Bank of Atmospheric and Hydrologic Research (DBAR) [39] (Figure 2). Some precipitation data missing during typhoons were replaced by the typhoon database of Taiwan Central Weather Bureau [31]. The observed daily streamflow and sediment concentration data during 2004–2015 were collected from the annual Hydrological Year Book of Taiwan Republic of China [40]. Most of the gages in the Chenyulan watershed seriously lack in measured data, and the Nemoupu gage is the only station that has complete data from 2004 to 2015 (Figure 2). The measured data during 2004–2009 and 2010–2015 were used for model calibration and validation, respectively. Since the sediment data were not measured continuously for every day, the sediment rating curve was first established to estimate the daily sediment concentration and loads for the model calibration and validation.
The land use data were collected form Second Land Use Survey in 2006 from [41] (Figure 3). Forest and agricultural lands are the two major land uses in the watershed, occupying 74.46% and 14.05% of the total area, respectively. The landslide area was set as barren (land use code in SWAT: BARR), accounting for 2.89% of total area. Soil data was collected from [42]. Approximately 82% of the total area has not been surveyed, thus it was assumed to be darkish colluvial soil in this study [30]. The pale colluvial clay accounts for 12% of the study area. The slope was calculated from the DEM data through the GIS spatial analysis tool. The slope was divided into five classes (0–9%, 9–30%, 30–45%, 45–60%, >60%), occupying 3.29%, 12.25%, 15.1%, 19.62%, and 49.74% of the total area, respectively.

2.2. SWAT Model

2.2.1. Model Description

The soil and water assessment tool (SWAT) model, a semi-distributed watershed-based hydrological model, was developed by the Agricultural Research Service, United States Department of Agriculture in 1994 [43]. The model can simulate the hydrological processes and sediment/nutrient transports at various spatio-temporal scales. In the SWAT model, several sub-basins are first delineated by setting the stream outlets (Figure 2), and further a unique combination of land use, soil and slope forms the basic modeling unit, called the hydrologic response unit (HRU). In this study, a total of 1173 HRUs of 23 sub-basins in the Chenyulan watershed were delineated (Figure 4).
The SWAT model simulates the surface runoff by using the curve number method, developed by the U.S. Soil Conservation Service (SCS). The SCS curve number equation is as follows:
Q surf = ( R day I a ) 2 ( R day I a ) + S
where Qsurf is surface runoff (mm); Rday is daily precipitation (mm); Ia is the initial loss of surface water storage, interception, and percolation (mm); S is retention parameter (mm), which changes with soil type, land use and slope.
In SWAT model, infiltration is calculated by the Green and Ampt infiltration method. For evapotranspiration, three methods of evaporation estimation can be selected, including FAO Penman–Monteith, Hargreaves and Priestly-Taylor method. The hydrologic process is calculated by Equation (2).
S t = S 0 + i = 1 t R day Q surf E a W seep Q gw
where St is soil water content (mm); S0 is initial soil water content (mm); Rday is the daily precipitation (mm); Qsurf is the daily surface runoff (mm); Ea is the daily evapotranspiration (mm); Wseep is the daily percolation (mm); Qgw is the daily baseflow (mm). In SWAT, the soil water content is calculated for different layers defined by the users, and 50% of the evaporative demand is extracted from the top 10 mm of soil [44].
SWAT model uses the modified universal soil loss equation (MUSLE) to estimate soil loss at the HRU scale.
S = 11.8 × ( Q surf × q × A ) 0.56 × K × C × P × LS × CFRG
where S is soil erosion (t); Qsurf is surface runoff (mm/ha); q is peak runoff (m3/s); A is the area of HRU (ha); K is the soil erodibility factor; C is the USLE land use/cover management factor; P is the USLE support practice factor; LS is the topographic factor; CFRG is the coarse fragment factor, which is calculated as the function of percent rock in the first soil layer. CFRG value equals 1 when there is no rock in the first soil layer. A higher percent of rock will result in smaller CFRG value and less soil loss.

2.2.2. Sediment Transport Methods

There are five sediment transport methods in SWAT model (version 664), which are the simplified Bagnold method with/without routing by particle size, Kodoatie method, Molinas and Wu method, and the Yang sand and gravel method. These methods can be applied for rivers of various river bed materials.
• Simplified Bagnold equation (EQN0 and EQN1):
The Bagnold method is the default method in SWAT model. The maximum sediment concentration is calculated as follows.
conc sed , ch , mx = c sp · v ch , pk spexp
where conc sed , ch , mx is the maximum sediment concentration (t/m3); c sp is linear coefficient; v ch , pk is peak flow velocity (m/s); spexp is exponential coefficient.
In EQN0, as the default and only one sediment transport method in SWAT 2005 version, the bed load is limited by the channel cover and erodibility factors, and the sediment carried by channel is always near the calculated maximum transport capacity [44]. Moreover, it does not keep track of sediment pools in various particle sizes. Thus, EQN1, additional stream power equation in SWAT 2016 version, has been incorporated with physics-based approach for channel erosion.
• Kodoatie equation (EQN2):
The Kodoatie equation is an optimized sediment transport equation, based on the non-linear sediment equation [45] and the field observation data.
conc sed , ch , mx = ( a · v ch b · y c · S d Q in ) · ( W + W btm 2 )
where conc sed , ch , mx is the maximum sediment concentration (tons/m3); v ch is mean flow velocity (m/s); y is mean flow depth (m); S is energy slope (m/m); Qin is water entering the reach (m3); a, b, c and d are regression coefficients for different bed materials; W is channel width at the water level (m); Wbtm is bottom width of the channel (m).
• Molinas and Wu equation (EQN3):
Molinas and Wu [46] developed the sediment transport equation based on the universal stream power for rivers of large sand bed. The sediment weight concentration is calculated as follows.
C w = M Ψ N
where Cw is sediment concentration by weight; Ψ is universal stream power; M and N are coefficients. The universal stream power (Ψ) is calculated as follows.
Ψ = { ( S g 1 ) · g · depth · ω 50 · [ log 10 ( depth D 50 ) ] } 0.5
where Sg is relative density of solid (2.65); g is acceleration of gravity (9.81 m/s2); depth is flow depth (m); ω 50 is fall velocity of median size particles (m/s); D50 is median sediment size.
• Yang sand and gravel equation (EQN4):
Developed by Yang [47], the sediment weight concentration was calculated based on the sediment size (D50), unit stream power (VchS), shear velocity ( V * ), fall velocity ( ω 50 ), and kinematic viscosity (υ). The equations are separated for sand and gravel bed material shown as follows.
Sand equation for median size (D50) less than 2 mm,
log C w = 5.435 0.286 log ω 50 D 50 υ 0.457 log V * ω 50 + ( 1.799 0.409 log ω 50 D 50 υ 0.314 log V * ω 50 ) log ( V ch S ω 50 V cr S ω 50 )
Gravel equation for median size (D50) between 2 mm and 10 mm,
log C w = 6.681 0.633 log ω 50 D 50 υ 4.816 log V * ω 50 + ( 2.784 0.305 log ω 50 D 50 υ 0.282 log V * ω 50 ) log ( V ch S ω 50 V cr S ω 50 )
where Cw is weight concentration (ppm); ω 50 is fall velocity of the median size sediment (m/s); υ is kinematic viscosity (m2/s); V * is shear velocity (m/s); Vch is mean channel velocity (m/s); Vcr is critical velocity (m/s); S is energy slope (m/m).

2.3. Taiwan Universal Soil Loss Equation (TUSLE)

In order to present the characteristics of Taiwan soil, we applied the Taiwan universal soil loss equation (TUSLE) [48], revised from the universal soil loss equation (USLE) in the SWAT model. The factor adjustments are shown as follows.
• Soil erodibility factor (K)
The adjusted K factors were based on the soil survey in Taiwan conducted by [49,50]. The K values in the study area range from 0.13 to 0.40 (Table 2), and higher value indicates the soil layer is easier to erode.
• Rainfall erosivity factor (R)
The R factor in USLE is calculated by 30-min maximum rainfall intensity and rainfall kinetic energy [6]. Since the R factor of USLE is complicated to calculate, Chen et al. [48] estimated the R factor for TUSLE by the regression equation developed by the U.S. Department of Agriculture, Agriculture Research Service (Equation (10)) [48].
R eq = 823.8 + 5.213 P r
where Req is rainfall erosivity factor; Pr is annual precipitation (mm).
In MUSLE, the rainfall erosivity factor is replaced by the runoff factor as the runoff factor could better represent the surface runoff and overland sediment transport characteristics [7]. Thus, in this study, we calculated the runoff factor (R) in MUSLE, instead of the rainfall erosivity factor (Req).
R = 11.8 ( Q · q · A ) 0.56
where R is runoff factor; Q is surface runoff (mm/ha); q is peak runoff (m3/s); A is the area of catchment (ha).
• Cover and management factor (C)
The C factor was estimated by non-linear equation with the normalized difference vegetation index (NDVI) to avoid overestimating the C values of area with low soil erosion rate [51] (Table 3).
NDVI 0 , C = ( 1 NDVI 2 ) 1 + NDVI
NDVI < 0 , { Building or non exposed ground ,   C = 0.01 Barren ,     C = 1.0
• Topographic factor (LS)
Wischmeier and Smith [52] established the LS equation (Equation (14)) as the product of L factor (Equation (15)) and S factor. In the SWAT model, the exponential factor (m) in the L factor equation is defined as Equation (16) [44], while TUSLE adopted the classification suggested by [52] that the exponential factor (m) varies with the slope, where m = 0.5, 0.4, 0.3 and 0.2 is used for the average slope greater than 5%, between 3–5%, between 1–3%, and less than 1%, respectively. Thus, TUSLE can increase the underestimated L factor at flatter slope and reduce the overestimated L factor at a steeper slope (Figure 5a). McCool et al. [53] indicated that Wischmeier and Smith’s [52] topographic factor equation could only be suitable for the slope from 0.1% to 18%, and developed S factor equation (Equations (17) and (18)) to more reasonably predict soil loss at steep topography. The comparison of S factor by [52,53] (Figure 5b) showed that the equation would overestimate S factor in areas of steeper slope. Therefore, we combined the L factor equation with m values from TUSLE and the S factor equation by [53] to calculate the LS factor.
LS = ( X 22.13 ) m ( 0.0654 + 4.56 sin θ + 65.4 sin 2 θ )
L = ( X 22.13 ) m
m = 0.6 × ( 1 e 35.835 * θ )
where X = slope length (m), m = exponential factor.
S = 10.8 sin θ + 0.03 , θ < 9 %
S = ( sin θ 0.0896 ) 0.6 , θ 9 %
where S is the slope factor of HRUs in MUSLE, θ is the slope of HRUs.

2.4. Calculation of Landslide Volume

The Chenyulan watershed has suffered by several severe natural disasters, such as Typhoon Herb (in 1996), 921 earthquake (in 1999), and Typhoon Morakot (in 2009), resulting in significant landslides, debris flows, change in landslide characteristics of central Taiwan [4]. Many studies have conducted the survey of landslide characteristic changed after the 921 earthquake [26,27,54,55,56,57]. Due to the lack of landslide calculation in the SWAT model, we integrated the landslide estimation equation developed by [28] into the SWAT model. The landslide volume was calculated by using the correlation between the landslide area and volume (R2 = 0.79) (Equation (19)).
ln ( V ) = 0.687 ln ( A ) + 2.326
where V is estimated landslide volume (m3); A is landslide area (m2).
It was reported that the percentage of landslide and debris flow in total sediment load would be greater than 60% when the daily cumulative precipitation is higher than 350 mm [58]. Therefore, the landslide volume estimation is triggered only when the daily precipitation exceeds 350 mm. The landslide area is read into the transformed landslide volume estimate equation as follows (Equation (20)).
V = e 2.326 · A 0.687
In this study, the landslide area was identified by the land use survey map. The landslide volume is calculated by Equation (19), and then added into the sediment yields of HRUs.

2.5. Model Calibration and Validation

We applied the SWAT-CUP (SWAT Calibration Uncertainty Program) [59] for model sensitivity analysis, calibration and validation. In SWAT-CUP, parameter uncertainty can be analyzed by using Sequential Uncertainty Fitting version 2 (SUFI2), Generalized Likelihood Uncertainty Estimation (GLUE), Particle Swarm Optimization (PSO), Parameter Solution (ParaSol), or Marko Chain Monte Carlo (MCMC). Each uncertainty analysis method can run multiple simulations and find the ranges of best parameters for the study project. Parameters at any particular sub-basin, land use, soil, slope, and even HRU can be individually calibrated to reflect the unique spatial characteristics.
Among these uncertainty analysis methods, SUFI2, ParaSol and GLUE are easier to calibrate the parameters [59]. It is suggested that 3000, 7500, 10,000, 100,000, and 45,000 times of simulation are needed for SUFI2, ParaSol, GLUE, PSO, and MCMC, respectively, in order to get satisfactory simulation results [59]. SUFI2 method was selected in this study as it requires the least number of simulations. In sensitivity analysis, p-value is used to distinguish whether parameters are sensitive or not. Parameters that have p-value smaller than or equal to 0.05 are considered as sensitive parameters for further calibration [59]. After sensitivity analysis, the selection of calibrated ranges and fitted values of the parameters are identified based on the statistical criteria (i.e., R2, NSE, PBIAS, and RSR), suggested by [60] with the model performance standards (Table 4).
R2 is the coefficient of determination, presenting the linear correlation between simulated and observed data. The value of R2 closer to 1 indicates a higher correlation.
R 2 = i = 1 n ( Y i s i m Y m e a n ) 2 i = 1 n ( Y i o b s Y m e a n ) 2
where Ysim is the simulated data; Yobs is observed data; Ymean is the average of observation.
Nash–Sutcliffe efficiency (NSE) presents the residuals of measured data [61], and the value ranges from −∞ to 1. NSE value that equals 1, indicates the simulation is same as the observation, while NSE > 0.5 is acceptable for SWAT model performance [60].
NSE = 1 [ i = 1 n ( Y i o b s Y i s i m ) 2 i n ( Y i o b s Y m e a n ) 2 ]
Percent bias (PBIAS) presents whether the simulated data are overestimated or underestimated. When PBIAS is greater than 0, the simulation is underestimated [62].
PBIAS ( % ) = i = 1 n ( Y i o b s Y i s i m ) i = 1 n Y i o b s × 100
RMSE-observation standard deviation ratio (RSR) is the ratio between root mean square error and standard deviation. The smaller the RSR is, the better simulation performance is [63].
RSR = i = 1 n ( Y i o b s Y s i m ) 2 i = 1 n ( Y i o b s Y m e a n ) 2

2.6. Sediment Load Estimation

In order to calibrate the daily sediment load, we estimated the observed daily sediment load by using the sediment rating curve, which describes the relationship between sediment concentration and water discharge. Since the Chenyulan watershed is a small mountainous watershed in Asia, we adopted the sediment rating curve method suggested by [64], who conducted studies in small mountainous watersheds in Japan.
We collected discontinuous data of sediment concentration and corresponding instant streamflow at the Nemoupu gauging station from 2004 to 2015 with a total of 338 data points. Both streamflow and sediment concentration data were first converted into logarithm and fitted with the linear regression model. A good relationship (r = 0.75) between sediment concentration and streamflow was found (Figure 6). Thus, the sediment rating curve was used for estimating the daily mean sediment concentration with the observed daily streamflow data, and the estimated daily sediment load could be further calculated.

3. Results

In this study, we compared three SWAT models, which are the official SWAT 2016 (version 664), SWAT-TUSLE and SWAT-Twn. The SWAT-TUSLE was modified with TUSLE which calculates C factor based on NDVI and L factors based on the slope. The SWAT-Twn was the integration of SWAT-TUSLE and landslide volume equation. Since we did not modify the streamflow-related equations in SWAT, the streamflow simulations are the same for all these three models. We first calibrated the official SWAT 2016 model for daily streamflow and compared the performance of sediment load simulations from SWAT 2016, SWAT-TUSLE, and SWAT-Twn.

3.1. Model Calibration and Validation for Streamflow

The observed streamflow from 2004 to 2015 was used for model calibration and validation. The streamflow-related parameters for calibration were referred to the previous study [30]. The parameters include curve number (CN2), plant uptake compensation factor (EPCO), surface runoff lag time (SURLAG), baseflow alpha factor (ALPHA_BF), effective hydraulic conductivity in main channel alluvium (CH_K2), and Manning’s “n” value for the main channel (CH_N2). In addition, we included Manning’s “n” value for the tributary channel (CH_N1) and effective hydraulic conductivity in tributary channel alluvium (CH_K1) as they are also sensitive in this study. In order to differentiate the characteristics of the parameters in various land use, slope, and soil, some parameters (i.e., CN2, ALPHA_BF, CH_K1, CH_K2, CH_N1, CH_N2, CH_K1) were individually calibrated for specific land use, slop and soil. Table 5 shows the calibrated ranges and fitted parameter values for daily streamflow simulation. The model did satisfactory and good performance for the calibration and validation, respectively (Figure 7), indicating the fitted streamflow parameters in model could well reflect the runoff characteristics of the Chenyulan watershed.

3.2. Comparison of SWAT 2016 and Modified SWAT Models

After calibrating the daily streamflow, we compared the uncalibrated simulated sediment yields from SWAT 2016, SWAT-TUSLE, and SWAT-Twn to quantify the impacts of using TUSLE and landslide volume equation on sediment yields at HRU and watershed levels (Table 6 and Table 7). It should be noted that we only used the sediment yield data during the streamflow calibration period because the fitted parameter values during validation period are different than those during calibration period. However, the range of parameter values are the same for both calibration and validation periods, and the simulation results can reflect the model uncertainty. Thus, we used the simulated sediment yield data during the calibration period (2004–2009) with calibrated fitted streamflow-related parameters and default sediment-related parameters to demonstrate the difference driven by different models.
The major differences between SWAT 2016 and the two other modified SWAT (SWAT-TUSLE and SWAT-Twn) are the LS factor, which has more influence in steep slope areas (slope > 9%), and the C factor, which was calculated by NDVI resulting various C factor values for different land uses. There were 77.12% and 80.29% of urban and agricultural lands located in areas with steeper slope (slope > 9%). Therefore, with unchanged C factor for urban, sediment yield from urban had decreased by approximate 60% due to the modified LS factors in TUSLE (Table 6). However, sediment yields from agricultural lands did not change significantly by modified SWAT models. It was because the C factor calculated by NDVI is doubled than the SWAT default value (Table 3), compensating the decrease in sediment yields by modified LS factor in TUSLE. Besides urban and agricultural lands, significant changes in sediment yields from forest, grassland and landslide were found. Although LS factors could influence the sediment yields, the C factor of forest and grassland which were changed from 0.001–0.003 to 0.2, played an important role in increases in sediment yields.
In the SWAT-Twn model, the landslide volume estimation is activated when the daily precipitation reaches over 350 mm. It should be noted that landslide volume estimation is only applied to the landslide area, not other land uses. It is obvious that sediment yields from landslide area increased significantly when the landslide volume estimation was activated in SWAT-Twn. Moreover, since forest is the main land use occupying 74.46% of the watershed and the NDVI-calculated C factor of forest is greater than SWAT default C factor, the annual sediment yields from the watershed were increased by 59.9% and 65.7% by SWAT-TUSLE and SWAT-Twn, respectively (Table 7). The increase of 5.8% of sediment yield by SWAT-Twn was due to landslide volume estimation at the landslide areas. It shows that landslide volume estimation should be considered as the major contribution to sediment yields.
Before calibrating the sediment, these models overestimated the daily sediment load in terms of great positive PBIAS values (Figure 8). However, SWAT-TUSLE and SWAT-Twn performed better than SWAT 2016 in terms of greater R2, NSE and smaller RSR, especially the SWAT-Twn performance had better statistical criteria values (R2 = 0.74, NSE = 0.66, RSR = 0.58). Therefore, we used SWAT-Twn for further sediment calibration and validation.

3.3. Model Calibration and Validation for Sediment

The SWAT-Twn model was calibrated and validated for sediment loads with five different sediment transport methods (i.e., EQN0–4). The calibration (2004–2009) and validation (2010–2015) periods for sediment loads were the same as those for streamflow. First, the sensitivity analyses for sediment parameters for different sediment transport methods were examined (Table 8). A total of eight sediment-related parameters were identified as sensitive parameters. Four of these parameters (i.e., SPCON, SPEXP, PRF_BSN, ADJ_PKR) are estimated on basin level (*.bsn), meaning that the parameter values are fixed for the entire watershed; while the rest of parameters (i.e., CH_COV1, CH_BNK_D50, CH_BED_D50) are estimated on reach level (*.rch), which could be varied by spatial and slope conditions.
The linear parameter (SPCON), exponent parameter (SPEXP) and peak rate adjustment factor (PRF_BSN) are only activated for the simplified Bagnold equation (EQN0 and EQN1). The peak rate adjustment factor (ADJ_PKR) was found to be sensitive to EQN0, EQN3 and EQN4. In order to identify the difference in the reach-level parameters, we separated the watershed by slope of 60% as almost half (49.58%) of the Chenyulan watershed is at a slope greater than 60%.
The channel bank vegetation coefficient for shear stress (CH_COV1) at the sub-basins with mean slope greater than 60% was sensitive for EQN2, EQN3 and EQN4, indicating that the vegetation at steeper slope areas would have great influence on sediment load compared to that at flatter slope areas. The median particle diameter of channel bank (CH_BNK_D50) at steeper slope areas was only sensitive for EQN4, while the median particle diameter of channel bed (CH_BED_D50) were sensitive for EQN1, EQN2, and EQN3. It shows that the median particle diameter of channel bank or bed should be measured for increasing the accuracy of sediment simulation with EQN1-4. Although both EQN0 and EQN1 are simplified Bagnold stream power equations, the EQN0 (default in SWAT 2005 version) does not keep track of sediment pools in various particle sizes, while the EQN1 (additional stream power equation in SWAT 2016) has been incorporated with physics-based approach for channel erosion. Moreover, the simulation of channel erosion with EQN0 is not partitioned between stream bank and stream bed. Thus, both CH_BNK_D50 and CH_BED_D50 are not sensitive for EQN0, while the EQN1 was sensitive with bed erosion (CH_BED_D50).
After identifying the sensitive parameters for those five sediment transport methods, the SWAT-Twn was calibrated and validated separately with each sediment transport method (Table 9 and Figure 9). Generally, the simulation results by EQN0 and EQN1 were better than those by other sediment transport methods, in terms of R2 and NSE greater than 0.5. It indicates that Bagnold equation is more suitable for the Chenyulan watershed. Moreover, the SWAT-Twn with EQN2, EQN3 or EQN4 was found to underestimate for peak flows and overestimate for low flows (Figure 9). Thus, the overestimated low flows at the most flow period resulted in great negative PBIAS values.

4. Discussion

4.1. Improvement of SWAT-Twn Model for Watershed Sediment Production

In the nearly-flat upland watersheds where surface runoff seldom occurs, sediment transported by surface runoff is usually significantly over predicted in a model [65]. Since the Chenyulan watershed is mountainous with hilly slope, modeling the sediment yield and transport responses with considering unique hydrological responses is important. Compared to SWAT 2016, the major modifications in SWAT-TUSLE are C factor and LS factor. The main land use in the Chenyulan watershed is forest, for which C factor is 0.2 by NDVI calculation and much greater than the SWAT default C factor (0.001). Although the decrease in LS factor would cause less sediment yields, the increase in sediment yield by the adjusted C factor was much greater than the decrease in sediment yield by the adjusted LS factor. Thus, it was found that some sub-basins (i.e., sub-basins no. 5, 8, 9, 12, 13, 16–19) where forest is dominated (>70%) generated more sediment yields (Figure 10 and Table 10). Moreover, some agricultural-dominated (>30%) sub-basins (i.e., sub-basins no. 1, 3, 7, 10, 11 15) are mostly located in central and northern parts of the watershed, where slope is less steep. Thus, the adjusted LS factor has a relatively small impact on sediment yields from agricultural-dominated sub-basins. Great change in sediment yield was found in the south part of the watershed, where both forest and grassland are dominated, as the adjusted C factors of pasture and urban are greater than the SWAT default ones. Since the landslide volume equation is the additional improvement in SWAT-Twn from SWAT-TUSLE, the sediment yields by SWAT-Twn were increased by 2%, 40% and 3% in the sub-basins no. 1, 7, and 17, respectively, where the landslide area is greater than 5% of the sub-basin area.

4.2. Calibration with Different Sedmient Transport Methods

Based on the sensitivity analysis (Table 8), the calibrated values are listed in Table 11. For EQN0 and EQN1, the calibrated ranges and fitted values of SPCON, SPEXP and PRF_BSN are similar. Moreover, ADJ_PKR shows similar fitted value (0.56–0.63) and calibrated ranges (0.5–1.5) for EQN0, EQN3, and EQN4, indicating the characteristics of peak flow could be well represented in different sediment transport methods. For the simplified Bagnold method (EQN0 and EQN1), channel erodibility is controlled by the channel erodibility factor (CH_COV1) ranging from 0 to 1. The default value (0) indicates non-erosive channel, while the value of 1 indicates no resistance to erosion. However, CH_COV1, which is conceptually similar to the soil erodibility factor used in the USLE equation, was not sensitive and thus we used the default value for the simulation with EQN0 and EQN1. For other physics-based methods (EQN2-4), the channel erodibility is calculated by shear stress, and the CH_COV1 is defined as channel bank vegetation coefficient for estimating critical shear stress [66]. CH_COV1 is sensitive for EQN2, EQN3, and EQN4, however, the fitted values are quite different. The fitted CH_COV1 values indicate that the channel vegetation is between sparse trees (CH_COV1 = 5.40) and dense trees (CH_COV1 = 19.20), sparse trees, and grassy (CH_COV1 = 1.97) for EQN2, EQN3, and EQN4, respectively. The smaller the CH_COV1 value is, the greater the channel erodibility coefficient is [44]. Interestingly, the reflection of CH_COV1 on the channel erodibility is consistent with the suitability of using the Kodatie model (EQN2) for the stream bed material size ranging from silt to gravel, Molinas and Wu model (EQN3) for primarily sand size particles, and Yang sand and gravel model (EQN4) for primarily sand and gravel size particles.
Moreover, due to the lack of the information on channel materials, the calibrated median particle size diameter of channel bank sediment (CH_BNK_D50) or bed sediment (CH_BED_D50) showed various results among different sediment transport methods. For EQN1 and EQN4, the median size of bank and bed sediments were identified as mostly much greater than large aggregate (500 µm), thus the particle size distribution for D50 (>2000 µm) assumed by SWAT is 15% of sand, 15% of silt, 5% of clay and 65% of gravel [44]. For EQN2, the calibrated CH_BED_D50 range is between 250 µm and 2000 µm, indicating the channels are characterized as medium to very coarse sand-bed materials [40]. Similarly, for EQN3, the fitted CH_BED_D50 value is between sand (200 µm) and large aggregate (500 µm) [44]. Therefore, for EQN2 and EQN3, the particle size distribution for D50 (50 to 2000 µm) assumed by SWAT is 65% of sand, 15% of silt, 15% of clay and 5% of gravel.
In SWAT, the particle size distribution is usually used for estimating the bank and bed erosion, and the percentage of sediment (sand, silt, clay, gravel, small aggregate and large aggregate) that gets deposited, is calculated by the fall velocity of median size particle. The amount of sediment transported out of the reach is calculated as follows: Amount of total load entering the reach at the beginning of the time period minus the amount of total load deposited in the reach plus the amount of sediment due to bank and bed erosion in the reach. The total load in the reach considered includes bed and suspended load coming from the foregoing reaches, as well as suspended load originating from the soil erosion of the surrounding sub-basin. Thus, further improvements (i.e., the assumption of particle size distribution for bank and bed erosion, the calculation of fall velocity for wide range of particle size) need to be done to calculate more reliable and accurate sediment loads in SWAT.

4.3. Selection of a Suitable Sediment Transport Method

Due to the fact that sediment loads are extremely high during heavy rainfall events, the comparison between measured and simulated data in linear-scaled plot was difficult to identify their difference in sediment loads of low values. Therefore, we applied the logarithmic scale for the measured and simulated sediment data with different sediment transport methods and compared their statistical criteria (R2 and PBIAS) during 2004–2016 (Figure 11 and Table 12). It should be noted that the simulation results are the same in Figure 9 and Figure 11. The simplified Bagnold method (EQN0 and EQN1) performed better for the low sediment loads with smaller PBIAS values. Although the model has been calibrated for sediment by SWAT-CUP, sediment loads are still overestimated for all sediment transport methods. It is because the SWAT-CUP model tends to adjust the parameters to meet the high observed sediment loads, but results in great PBIAS values. Thus, the SWAT model simulated high sediment loads better than the low sediment loads. Interestingly, when the measured and simulated sediment data are plotted in logarithmic scale, the R2 values increase and PBIAS values decrease (Table 12). The greater improvement in the statistical criteria values shows the higher degree of overestimation when applying the Molinas and Wu model (EQN3) and Yang sand and gravel model (EQN4). By comparing the sediment load data in linear and logarithmic scales, we can better identify the suitable sediment transport method and suggest that further calibration and validation are needed for log-transformed simulated sediment data in the SWAT-CUP model.

5. Conclusions

The Chenyulan watershed has suffered from serious landslide and debris flow induced by heavy rainfall and typhoons. In this study, we integrated the TUSLE and landslide volume estimation into the SWAT model as SWAT-Twn. By evaluating the simulated sediment yields from different land uses, the importance of topographic (LS) factor and NDVI-calculated weighted C factor were identified and landslide volume estimation should be taken into concern. The examination of five different sediment transport methods revealed some important issues. First, the level of sensitivity of sediment-related parameters is different for those sediment transport methods, and parameters (i.e., CH_COV1, CH_BNK_D50, CH_BED_D50) that are estimated on each level, are suggested to be calibrated by spatial and slope conditions. Second, it is more accurate to investigate the channel vegetation (CH_COV1) and measure the particle sizes of channel bank and bed sediment (CH_BNK_D50 and CH_BED_D50). The calibrated parameter values by SWAT-CUP for different sediment transport methods may be misleading. Third, the particle size distribution assumed by SWAT is suggested to be an option that can be edited by users. Furthermore, the calculation of fall velocity is suggested to not be only limited for median particle size as it would be biased for channels of wide range of particle sizes. Last but not the least, like the streamflow simulation in SWAT and SWAT-CUP, an option for the user to compare and plot the sediment simulation in logarithmic scale would provide more insights into sediment calibration. In sum, the SWAT-Twn model performed better than SWAT 2016 and SWAT-TUSLE, as TUSLE calculated less sediment at steep area, resulting reasonable sediment export simulation at low flow condition and landslide volume estimation reflected the real situation. Additional improvements in SWAT and SWAT-CUP need to be made to better predict the sediment yields and loads at mountainous watersheds.

Author Contributions

Conceptualization, C.-M.L. and L.-C.C.; Formal analysis, C.-M.L.; Funding acquisition, L.-C.C.; Methodology, C.-M.L.; Resources, L.-C.C.; Software, C.-M.L.; Supervision, L.-C.C.; Validation, L.-C.C.; Visualization, C.-M.L.; Writing—original draft, C.-M.L. and L.-C.C.; Writing—review & editing, L.-C.C.

Funding

This research was funded by the Ministry of Science and Technology, Taiwan (MOST 107-2625-M-239-001).

Acknowledgments

We thank Wen-Cheng Liu from National United University, and Yung-Chieh Wang and Chia-Jeng Chen from National Chung Hsing University for constructive comments that improved the manuscript.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

References

  1. Garcia-Ruiz, J.M.; Begueria, S.; Nadal-Romero, E.; Gonzalez-Hidalgo, J.C.; Kana-Renault, N.; Sanjuan, Y. A meta-analysis of soil erosion rates across the world. Geomorphology 2015, 239, 160–173. [Google Scholar] [CrossRef] [Green Version]
  2. Freebairn, D.M.; Loch, R.J.; Silburn, D.M. Chapter 9, Soil erosion conservation for vertisols. In Developments in Soil Science; Elsevier: Amsterdam, The Nederlands, 1996; Volume 24, pp. 303–362. [Google Scholar]
  3. Garcia, M.H. Sediment Engineering; American Society of Civil Engineers: Reston, VA, USA, 2008. [Google Scholar]
  4. Dadson, S.J.; Hovius, N.; Chen, H.; Dade, W.B.; Hsieh, M.-L.; Willett, S.D.; Hu, J.-C.; Horng, M.-J.; Chen, M.-C.; Stark, C.P.; et al. Links between erosion, runoff variability and seismicity in Taiwan orogen. Nature 2003, 426, 648–651. [Google Scholar] [CrossRef] [PubMed]
  5. Mullan, D.; Matthews, T.; Vandaele, K.; Barr, I.D.; Swindles, G.T.; Meneely, J.; Boardman, J.; Murphy, C. Climate impacts on soil erosion and muddy flooding at 1.5 versus 2C warming. Land Degrad. Dev. 2019, 30, 94–108. [Google Scholar] [CrossRef]
  6. Wischmeier, W.H.; Smith, D.D. Rainfall energy and its relationship to soil loss. Trans. Am. Geophys. Union 1958, 39, 285–291. [Google Scholar] [CrossRef]
  7. Williams, J.R. Sediment yield prediction with Universal Equation using runoff equation. In Present and Prospective Technology for Predicting Sediment Yields and Sources; USA Department of Agriculture: New Orleans, LA, USA, 1975; pp. 244–252. [Google Scholar]
  8. Sadeghi, S.H.R.; Gholami, L.; Khaledi Darvishan, A.; Saeidi, P. A review of the application of the MUSLE model worldwide. Hydrol. Sci. J. 2014, 59, 365–375. [Google Scholar] [CrossRef]
  9. Beasley, D.D.; Huggins, L.F.; Monke, E.J. ANSWERS: A model for watershed planning. Trans. ASAE 1980, 23, 938–944. [Google Scholar] [CrossRef]
  10. Williams, J.R.; Jones, C.A.; Dyke, P.T. A modeling approach to determining the relationship between erosion and soil productivity. Trans. ASAE 1984, 27, 129–144. [Google Scholar] [CrossRef]
  11. Nearing, M.A.; Foster, G.R.; Lane, L.J.; Finkner, S.C. A process-based soil erosion model for USDA-Water erosion prediction project technology. Trans. ASAE 1989, 32, 1587–1593. [Google Scholar] [CrossRef]
  12. Morgan, R.P.C.; Qionton, J.N.; Simuth, R.E.; Govers, G.; Poesen, J.W.A.; Auerswald, K.; Chisci, G.; Torri, D.; Styczen, M.E. The European Soil Erosion Model (EUROSEM): A dynamic approach for predicting sediment transport form fields and small catchments. Earth Surf. Process. Landf. 1998, 23, 527–544. [Google Scholar] [CrossRef]
  13. Sadeghi, S.H.R.; Mizuyama, T. Applicability of the Modified Universal Soil Loss Equation for prediction of sediment yield in Khanmirza watershed, Iran. Hydrol. Sci. J. 2007, 52, 1068–1075. [Google Scholar] [CrossRef]
  14. Merrtitt, W.S.; Letcher, R.A.; Jakeman, A.J. A review of erosion and sediment transport models. Environ. Model. Softw. 2003, 18, 761–799. [Google Scholar] [CrossRef]
  15. Hairsine, P.B.; Rose, C.W. Modeling water erosion due to overland flow using physical principles: 1. Sheet flow. Water Resour. Res. 1992, 28, 237–243. [Google Scholar] [CrossRef]
  16. Wheater, H.S.; Jakeman, A.J.; Beven, K.J. Progress and directions in rainfall-runoff modeling. In Modelling Change in Environmental Systems; Jakeman, A.J., Beck, M.B., McAleer, M.J., Eds.; John Wiley and Sons: Hoboken, NJ, USA, 1993; pp. 101–132. [Google Scholar]
  17. Abbott, M.B.; Bathurst, J.C.; Cunge, J.A.; O’Connell, P.E.; Rasmussen, J. An introduction to the European Hydrological System-Système Hydrologique Europèen SHE. 1. History and philosophy of a physically-based, distributed modelling system. J. Hydrol. 1986, 87, 45–59. [Google Scholar] [CrossRef]
  18. Storm, D.E.; White, M.J.; Stoodley, S. Modeling Non-Point Source Component for the Fort Cobb TMDL Final Report; Oklahoma Department of Environmental Quality: Oklahoma, OK, USA, 2003.
  19. Pohlert, T.; Huisman, J.A.; Breuer, L.; Frede, H.G. Modeling of river Dill, Germany. Adv. Geosci. 2005, 5, 7–12. [Google Scholar] [CrossRef]
  20. Yevenes, M.A.; Mannaerts, C.M. Seasonal and land use impacts on the nitrate budget and export of a mesoscale catchment in Southern Portugal. Agric. Water Manag. 2011, 102, 54–65. [Google Scholar] [CrossRef]
  21. Arnold, J.G.; Srinivasan, R.; Muttian, R.S.; Williams, J.R. Large area hydrologic modeling and assessment, part I: Model development. J. Am. Water Resour. Assoc. 1998, 34, 73–89. [Google Scholar] [CrossRef]
  22. Prosser, I.P.; Rustomji, P. Sediment transport capacity relations for overland flow. Prog. Phys. Geogr. 2000, 24, 179–193. [Google Scholar] [CrossRef]
  23. Addis, H.K.; Strohmeier, S.; Ziadat, F.; Melaku, N.D.; Klik, A. Modeling streamflow and sediment using SWAT in Ethiopian Highlands. Int. J. Agric. Biol. Eng. 2016, 9, 51–66. [Google Scholar]
  24. Bressiani, D.A.; Gassman, P.W.; Fernandes, J.G.; Garbossa, L.H.P.; Srinivasan, R.; Bonumá, N.B.; Mendiondo, E.M. Review of Soil and Water Assessment Tool (SWAT) applications in Brazil: challenges and prospects. Int. J. Agric. Biol. Eng. 2015, 8, 9–35. [Google Scholar]
  25. Arnold, J.G.; Moriasi, D.N.; Gassman, P.W.; Abbaspour, K.C.; White, M.J.; Srinivasan, R.; Santhi, C.; Harmel, R.D.; van Griensven, A.; Van Liew, M.W.; et al. SWAT: model use, calibration, and validation. Am. Soc. Agric. Biol. Eng. 2012, 55, 1491–1508. [Google Scholar]
  26. Lee, T.-Y.; Huang, J.-C.; Lee, J.-Y.; Jien, S.-H.; Zehetner, F.; Kao, S.-S. Magnified sediment export of small mountainous rivers in Taiwan: chain reaction from increased rainfall intensity under global warming. PLoS ONE 2015, 10, e0138283. [Google Scholar] [CrossRef]
  27. Chiu, Y.-J.; Lee, H.-Y.; Wang, T.-L.; Yu, J.; Lin, Y.-T.; Yuan, Y. Modeling sediment yields and stream stability due to sediment-related disaster in Shihmen Reservoir watershed in Taiwan. Water 2019, 11, 332. [Google Scholar] [CrossRef]
  28. Chang, C.-T.; Harrison, J.F.; Huang, Y.-C. Modeling typhoon-induced alterations on river sediment transport and turbidity based on dynamic landslide inventories: Gaoping River Basin, Taiwan. Water 2015, 7, 6910–6930. [Google Scholar] [CrossRef]
  29. Chiang, L.-C.; Lin, Y.-P.; Huang, T.; Schmeller, D.S.; Verburg, P.H.; Liu, Y.-L.; Ding, T.S. Simulation of ecosystem service responses to multiple disturbances from an earthquake and several typhoons. Landsc. Urban Plan. 2014, 122, 41–55. [Google Scholar] [CrossRef]
  30. Chiang, L.-C.; Chuang, Y.-T.; Han, C.-C. Integrating landscape metrics and hydrologic modeling to assess the impact of natural disturbances on ecohydrological processes in the Chenyulan watershed, Taiwan. Int. J. Environ. Res. Public Health 2019, 122, 266. [Google Scholar] [CrossRef]
  31. Typhoon Database-Taiwan Central Weather Bureau. Available online: https://rdc28.cwb.gov.tw/TDB (accessed on 21 June 2019).
  32. Abbaspour, K.C.; Rouhokahnejad, E.; Vaghefi, S.; Srinivasan, R.; Yang, H.; Klove, B. A continental-scale hydrology and water quality model for Europe: calibration and uncertainty of a high-resolution large-scale SWAT model. J. Hydrol. 2015, 524, 733–752. [Google Scholar] [CrossRef]
  33. Cibin, R.; Trybula, E.; Chaubey, I.; Brouder, S.M.; Volenec, J.J. Watershed-scale impacts of bioenergy crops on hydrology and water quality using improved SWAT model. Glob. Chang. Biol. Bioenergy 2016, 8, 837–848. [Google Scholar] [CrossRef]
  34. Liu, R.; Xu, F.; Zhang, P.; Yu, W.; Men, C. Identifying non-point source critical source areas based on multi-factors at a basin scale with SWAT. J. Hydrol. 2016, 533, 379–388. [Google Scholar] [CrossRef]
  35. Nguyen, H.H.; Reckmagel, F.; Meyer, W.; Frizenschaf, J.; Shrestha, M.K. Modelling the impacts of altered management practices, land use and climate changes on the water quality of the Millbrook catchment-reservoir system in South Australia. J. Environ. Manag. 2017, 202, 1–11. [Google Scholar] [CrossRef]
  36. Jang, S.S.; Ahn, S.R.; Kim, S.J. Evaluation of executable best management practices in Haean highland agricultural catchment of South Korea using SWAT. Agric. Water Manag. 2017, 180, 224–234. [Google Scholar] [CrossRef]
  37. Meng, X.; Wang, H.; Shi, C.; Wu, Y.; Ji, X. Establishment and evaluation of the China Meteorological Assimilation Driving Datasets for the SWAT model (CMADS). Water 2018, 10, 1555. [Google Scholar] [CrossRef]
  38. Taiwan Geospatial One Stop (TGOS). Available online: https://www.tgos.tw (accessed on 21 June 2019).
  39. Taiwan Data Bank of Atmospheric and Hydrologic Research (DBAR). Available online: https://www.dbar.pccu.edu.tw (accessed on 21 June 2019).
  40. Annual Hydrological Year Book of Taiwan, Republic of China. Available online: https://gweb.wra.gov.tw/wrhygis (accessed on 21 June 2019).
  41. Taiwan MAP Service-National Land Surveying and Mapping Center. Available online: https://maps.nlsc.gov.tw (accessed on 21 June 2019).
  42. Taiwan Soil Resource and Cultivate Land Cover Map Searching System. Available online: https://farmcloud.tari.gov.tw/SOA (accessed on 21 June 2019).
  43. Srinivasan, R.; Arnold, J.G. Integration of a basin-scale water quality model with GIS. Water Resour. Bull. 1994, 30, 453–462. [Google Scholar] [CrossRef]
  44. Neitsch, S.L.; Arnold, J.G.; Kiniry, J.R.; Williams, J.R. Soil and Water Assessment Tool Theoretical Documentation; Texas Water Resources Institute Technical Report No. 406; Texas Water Resources Institute: Temple, TX, USA, 2011. [Google Scholar]
  45. Posada-Garcia, L. Transport of sands in Deeps Rivers. Ph.D. Thesis, Colorado State University, Fort Collins, CO, USA, 1995. [Google Scholar]
  46. Molinas, A.; Wu, B. Transport of sediment in large sand-bed rivers. J. Hydraul. Res. 2001, 39, 135–146. [Google Scholar] [CrossRef]
  47. Yang, C.T. Sediment Transport Theory and Practice; The McGaw-Hill Companies, Inc.: New York, NY, USA, 1996. [Google Scholar]
  48. Chen, S.-C.; Wu, C.-Y.; Wu, Y.-L.; Wang, S.-H. Taiwan Universal Soil Loss Equation (TUSLE) based on revised factors and GIS layers—an example from the Shihmem Reservoir Watershed. J. Chin. Soil Water Conserv. 2009, 40, 185–197. (In Chinese) [Google Scholar]
  49. Wan, S.-S.; Huang, J.-I. Soil erosion of slope in Taiwan. J. Chin. Soil Water Conserv. 1989, 20, 127–144. (In Chinese) [Google Scholar]
  50. Hsieh, J.-S.; Wang, M.-G. Major Soil Maps in Taiwan; National Chung Hsing University Soil Survey and Testing Center: Taiwan, 1991. (In Chinese) [Google Scholar]
  51. Lin, W.-T. Automated Watershed Delineation for Spatial Information Extraction and Slopeland Yield Evaluation. Ph.D. Thesis, National Chung Hsing University, Taiwan, 2002. (In Chinese). [Google Scholar]
  52. McCool, D.K.; Brown, L.C.; Foster, G.R.; Mutchler, C.K.; Meyer, L.D. Revised slope steepness factor for the Universal Soil Loss Equation. Trans. ASAE 1987, 30, 1387–1396. [Google Scholar] [CrossRef]
  53. Wischmeier, W.H.; Smith, D.D. Predicting rainfall erosion losses—A guide to conservation planning. In Agriculture Handbook No. 537; USA Department of Agriculture: New Orleans, LA, USA, 1978. [Google Scholar]
  54. Hovius, N.; Stark, C.P.; Allen, P.A. Sediment flux from a mountain belt derived by landslides mapping. Geology 1997, 25, 231–234. [Google Scholar] [CrossRef]
  55. Guzzetti, F.; Ardizzone, F.; Cardinali, M.; Galli, M.; Reichenbach, P.; Rossi, M. Distribution of landslides in the Upper Tiber River basin, central Italy. Geomorphology 2008, 96, 105–122. [Google Scholar] [CrossRef]
  56. Dong, J.J.; Lee, C.T.; Tung, Y.H.; Liu, C.N.; Lin, K.P.; Lee, Y.J. The role of sediment budget in understanding debris flow susceptibility. Earth Surf. Process Landf. 2009, 34, 1612–1624. [Google Scholar] [CrossRef]
  57. Chang, C.-W.; Lin, P.-S.; Tsai, C.-L. Estimation of sediment volume of debris flow caused by extreme rainfall in Taiwan. Eng. Geol. 2011, 123, 83–90. [Google Scholar] [CrossRef]
  58. Liu, G.-F.; Lai, J.-S.; Chen, Y.-C.; Chiu, Y.-J.; Shen, C.-W.; Liang, W.-S.; Chiang, M.-H. A Study of Flood Control and Sediment Management Due to Climate Change of Jhoushuei River; Water Resource Agency: Taiwan, 2015. (In Chinese)
  59. Abbaspour, K.C. SWAT-CUP: SWAT Calibration and Uncertainty Programs—A User Manual; Swiss; EAWAG: Dübendorf, Switzerland, 2015. [Google Scholar]
  60. Moriasi, D.N.; Arnold, J.G.; Van Liew, M.W.; Bingner, R.L.; Haemel, R.D.; Veith, T.L. Model evaluation guidelines for systematic qualification of accuracy in watershed simulation. Trans. ASABE 2007, 50, 282–290. [Google Scholar] [CrossRef]
  61. Nash, J.E.; Sutcliffe, J.V. River flow forecasting through conceptual models 1. A discussion of principles. J. Hydrol. 1970, 10, 282–290. [Google Scholar] [CrossRef]
  62. Gupta, H.V.; Sorooshian, S.; Yapo, P.O. Status of automatic calibration for hydrologic models: Comparison with miltilevel expert calibration. J. Hydrol. Eng. 1999, 4, 135–143. [Google Scholar] [CrossRef]
  63. Singh, J.H.; Knaoo, H.V.; Arnold, J.G.; Demissie, M. Hydrologic modeling of the Iroquois River watershed using HSPF and SWAT. J. Am. Water Resour. Assoc. 2004, 41, 361–375. [Google Scholar]
  64. Sadeghi, S.H.R.; Mizuyama, T.; Miyata, S.; Gomi, T.; Kosugi, K.; Fukushima, T.; Mizugaki, S.; Onda, Y. Development, evaluation and interpretation of sediment rating curves for a Japanese small mountainous reforested watershed. Geoderma 2008, 144, 198–211. [Google Scholar] [CrossRef]
  65. Mitchell, J.K.; Banasik, K.; Hirschi, M.C.; Cooke, R.A.C.; Kalita, P. There is not always surface runoff and sediment transport. In Proceedings of the International Symposium on Soil Erosion Research for the 21st Century, Honolulu, HI, USA, 3–5 January 2001. [Google Scholar]
  66. Julian, J.P.; Torres, R. Hydraulic erosion of cohesive river banks. Geomorphology 2006, 76, 193–206. [Google Scholar] [CrossRef]
Figure 1. Study area.
Figure 1. Study area.
Water 11 01749 g001
Figure 2. Locations of weather stations, sub-basins, and gaging station.
Figure 2. Locations of weather stations, sub-basins, and gaging station.
Water 11 01749 g002
Figure 3. Land use, soil, and slope distribution.
Figure 3. Land use, soil, and slope distribution.
Water 11 01749 g003
Figure 4. Sub-basin distribution.
Figure 4. Sub-basin distribution.
Water 11 01749 g004
Figure 5. Comparison of L factor and S factor.
Figure 5. Comparison of L factor and S factor.
Water 11 01749 g005
Figure 6. Sediment rating curve.
Figure 6. Sediment rating curve.
Water 11 01749 g006
Figure 7. Streamflow calibration and validation results.
Figure 7. Streamflow calibration and validation results.
Water 11 01749 g007
Figure 8. Simulated daily sediment yield by three models with calibrated streamflow.
Figure 8. Simulated daily sediment yield by three models with calibrated streamflow.
Water 11 01749 g008
Figure 9. Sediment calibration and validation by different sediment transport methods.
Figure 9. Sediment calibration and validation by different sediment transport methods.
Water 11 01749 g009
Figure 10. Sediment yields at sub-basins simulated by three models.
Figure 10. Sediment yields at sub-basins simulated by three models.
Water 11 01749 g010
Figure 11. Comparison of simulated and measured sediment loads in logarithmic scale.
Figure 11. Comparison of simulated and measured sediment loads in logarithmic scale.
Water 11 01749 g011
Table 1. Top 10 typhoons that influenced the study area during 2004–2015.
Table 1. Top 10 typhoons that influenced the study area during 2004–2015.
RankTyphoonOccurred PeriodAccumulated Precipitation (mm)Weather Station
1Morakot5–10 August 20091953C1M440
2Matmo21–23 July 20141136C1M440
3Sinlaku11–16 September 20081072.5C0H9A0
4Mindulle28 June–3 July 2004774.5C0H9A0
5Krosa4–10 October 2007771C0H9A0
6Haitang16–20 July 2005760.5C0H9A0
7Saola30 July–3 August 2012693.5C1I310
8Sepat16–19 August 2007660.5C1M440
9Soulil11–13 July 2013624.5C0H9A0
10Jangmi26–29 September 2008593C0H9A0
(Source: Typhoon Database, Central Weather Bureau, Taiwan [31]).
Table 2. K factor used in the Chenyulan watershed.
Table 2. K factor used in the Chenyulan watershed.
Soil TypeAlluvial SoilLithosolDarkish Colluvial SoilRed SoilPale Colluvial SoilTaiwan ClayYellow Soil
K factor0.40.30.360.130.190.20.29
Table 3. Normalized difference vegetation index (NDVI)-calculated weighted C value for different land uses.
Table 3. Normalized difference vegetation index (NDVI)-calculated weighted C value for different land uses.
NDVI-Calculated CNumber of Grids
GrassForestAgricultureLandslide
0.01391461,33400
0.1304967,8198463707
0.28417153,46826,0072362
0.3412059,36514,1411590
0.4408456,68917,4222793
10092058026
Weighted C0.20.20.40.7
SWAT Default C0.0030.0010.20.2
Table 4. Statistical criteria for model performance [60].
Table 4. Statistical criteria for model performance [60].
Model PerformanceNSEPBIAS (%)RSR
FlowSedimentNutrient
Very Good0.75 < NSE ≤ 1.00|PBIAS| < 10|PBIAS| < 15|PBIAS| < 250.00 ≤ RSR ≤ 0.50
Good0.65 < NSE ≤ 0.7510 ≤ |PBIAS| < 1515 ≤ |PBIAS| < 3025 ≤ |PBIAS| < 400.50 < RSR ≤ 0.60
Satisfactory0.5 < NSE ≤ 0.6515 ≤ |PBIAS| < 2530 ≤ |PBIAS| < 5540 ≤ |PBIAS| < 700.60 < RSR ≤ 0.70
UnsatisfactoryNSE ≤ 0.5|PBIAS| ≥ 25|PBIAS| ≥ 55|PBIAS| ≥ 70RSR > 0.70
Table 5. Calibration results for streamflow.
Table 5. Calibration results for streamflow.
ParameterUnitClass 1DefaultCalibrated Value
FittedMin.Max.
CN2-FRST6037 (−38.3%)35 (−41.9%)57 (−4.7%)
RNGE6942 (−38.9%)39 (−44.1%)61 (−12.1%)
AGRL7747 (−38.9%)43 (−44.1%)68 (−12.1%)
EPCO-All10.340.100.44
SURLAG-All420.773.2423.11
ALPHA_BF1/dayDownstream0.0480.310.180.53
Sub. mean slope > 60%0.380.060.43
CH_K2mm/hDownstream0543.86427.89811.25
Sub. mean slope > 60%0546.43386.38744.02
Head stream0515.61293.20579.25
CH_N2-Sub. mean slope > 60%0.0140.130.100.25
Head stream0.370.240.46
CH_K1-Downstream058.8853.76161.54
CH_N1-Sub. mean slope > 60%0.01411.022.5211.08
1 FRST, RNGE, AGRL denote forest, grassland, agricultural land, respectively; three groups of sub-basin by slope are downstream (sub-basin no. 1–3, 5–7, 9–11, 13–15), sub-basin mean slope that greater than 60% (sub-basin no. 4, 8, 12, 16, 18, 19) and head stream (sub-basin no. 17, 20–23).
Table 6. Sediment yields at hydrologic response unit (HRU) level.
Table 6. Sediment yields at hydrologic response unit (HRU) level.
Land UseSediment Yield (t/ha)Difference (%)
SWAT 2016SWAT-TUSLESWAT-TwnSWAT-TUSLESWAT-Twn
Urban315.8130.3130.3−58.7−58.7
Agricultural land1201.01196.01196.0−0.4−0.4
Forest212.6724.4724.4240.7240.7
Grassland271.1933.1933.1244.2244.2
Landslide5461.61949.16872.3−64.325.8
Table 7. Sediment yields at watershed level.
Table 7. Sediment yields at watershed level.
YearSediment Yield (t/ha)Difference (%)
SWAT 2016SWAT-TUSLESWAT-TwnSWAT-TUSLESWAT-Twn
20048381335141159.268.3
20057561295132071.374.6
200610281574163553.159.0
20076711133116468.773.3
200811111561162340.546.2
20097461336138279.285.3
Average8581372142259.965.7
Table 8. The sensitivity analysis of sediment-related parameters (p-value < 0.05).
Table 8. The sensitivity analysis of sediment-related parameters (p-value < 0.05).
ParameterFile NameEQN0EQN1EQN2EQN3EQN4
SPCON*.bsnVV
SPEXP*.bsnVV
PRF_BSN*.bsnVV
ADJ_PKR*.bsnV VV
CH_COV1*.rte V 3V 3V 3
CH_BNK_D50 1*.rte V 3
CH_BED_D50 1*.rte V 2V 3V 3
1 unit: μm; 2 the parameter is sensitive for sub-basins with mean slope < 60%; 3 the parameter is sensitive for sub-basins with mean slope > 60%.
Table 9. Statistical results for sediment calibration and validation.
Table 9. Statistical results for sediment calibration and validation.
Statistical ParameterEQN0EQN1EQN2EQN3EQN4
Cal. 1Val. 2Cal.Val.Cal.Val.Cal.Val.Cal.Val.
R20.740.550.660.660.240.570.280.580.340.57
NSE0.730.730.650.650.220.220.250.250.310.31
PBIAS (%)−89.8−81.4−116.8−74.5−164.1−216.4−226.5−168.7−240.8−314.1
RSR0.520.520.590.590.880.880.860.860.830.83
1 Calibration; 2 Validation.
Table 10. Land use distribution (%) in sub-basins.
Table 10. Land use distribution (%) in sub-basins.
Sub-BasinArea (km2)WaterGrasslandUrbanAgricultural LandLandslideForest
134.6710.5403.1740.435.2940.05
216.037.1501.9829.511.159.2
30.7523.49014.6135.85026.52
430.656.364.750.2717.294.5766.88
524.630.1500.1220.791.377.02
630.045.070.061.0332.23.1158.91
71.8430.920044.8713.2311.46
843.312.0410.620.031.142.9683.54
917.111.897.760.5411.32.1675.71
1014.318.790.282.5132.454.0652.38
110.9718.3405.2543.051.6432.2
1210.306.6300.881.7691.21
1314.731.277.5016.723.870.71
144.765.8503.9823.880.8665.9
159.75.5701.7831.112.3459.68
1627.212.824.560.8611.920.7979.52
17161.0412.3804.226.5375.31
1824.093.054.161.39.682.0280.26
1913.12.494.2106.221.9785.59
2033.650.017.2300.463.5189
2142.491.73.290.143.973.0287.51
2217.7903.5500.161.0595.16
2320.6806.27000.7693.01
Table 11. Calibrated sediment-related parameters for different sediment transport methods.
Table 11. Calibrated sediment-related parameters for different sediment transport methods.
ParameterEQN0EQN1EQN2EQN3EQN4
Range (Fitted)Range (Fitted)Range (Fitted)Range (Fitted)Range (Fitted)
SPCON0.01–0.10 (0.012)0.01–0.10 (0.042)---
SPEXP1.02–1.20 (1.11)1.00–1.40 (1.25)---
PRF_BSN1.00–1.20 (1.02)0.80–1.20 (0.95)---
ADJ_PKR0.50–1.20 (0.61)--0.50–1.20(0.63)0.50–1.50 (0.56)
CH_COV1--7.00-19.00 (10.51)5.00–15.00 (0.13)1.00–9.00 (1.70)
CH_BNK_D50 1----6500–9000 (7956.25)
CH_BED_D50 1-500–6000 (2411)100–1800 (605.75)350–1000 (377.63)-
1 unit: µm.
Table 12. Statistical results for simulated sediment loads in linear and logarithmic scales.
Table 12. Statistical results for simulated sediment loads in linear and logarithmic scales.
MethodPlot ScaleR2PBIAS (%)
EQN0Linear0.71−86.8
Logarithmic0.64−34.3
EQN1Linear0.66−101.6
Logarithmic0.68−36.4
EQN2Linear0.28−182.8
Logarithmic0.68−48.6
EQN3Linear0.31−205.7
Logarithmic0.64−51.3
EQN4Linear0.37−267.1
Logarithmic0.70−55.4

Share and Cite

MDPI and ACS Style

Lu, C.-M.; Chiang, L.-C. Assessment of Sediment Transport Functions with the Modified SWAT-Twn Model for a Taiwanese Small Mountainous Watershed. Water 2019, 11, 1749. https://doi.org/10.3390/w11091749

AMA Style

Lu C-M, Chiang L-C. Assessment of Sediment Transport Functions with the Modified SWAT-Twn Model for a Taiwanese Small Mountainous Watershed. Water. 2019; 11(9):1749. https://doi.org/10.3390/w11091749

Chicago/Turabian Style

Lu, Chih-Mei, and Li-Chi Chiang. 2019. "Assessment of Sediment Transport Functions with the Modified SWAT-Twn Model for a Taiwanese Small Mountainous Watershed" Water 11, no. 9: 1749. https://doi.org/10.3390/w11091749

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop