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

: In Taiwan, the steep landscape and highly vulnerable geology make it di ﬃ cult 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 streamﬂow 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 o ﬃ cial SWAT model (version 664). The SWAT-Twn model coupled with the simpliﬁed Bagnold equation could estimate sediment export more accurately and signiﬁcantly reduce the overestimated sediment yield by 65.7%, especially in highly steep areas. Furthermore, ﬁve di ﬀ erent sediment transport methods (simpliﬁed 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 modiﬁed universal soil loss equation (MUSLE) factors and carefully evaluate the sediment transportation equations in SWAT.


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

Study Area and Data
The Chenyulan watershed, located in central Taiwan, has an area of 448 km 2 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 m 3 /s and 1860 m 3 /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)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(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  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.

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

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: where Q surf is surface runoff (mm); R day is daily precipitation (mm); I a 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. 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: 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).
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. 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).
where S t is soil water content (mm); S 0 is initial soil water content (mm); R day is the daily precipitation (mm); Q surf is the daily surface runoff (mm); E a is the daily evapotranspiration (mm); W seep is the daily percolation (mm); Q gw 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.
where S is soil erosion (t); Q surf is surface runoff (mm/ha); q is peak runoff (m 3 /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.

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, 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 (4) where conc sed,ch,mx is the maximum sediment concentration (t/m 3 ); 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.
where conc sed,ch,mx is the maximum sediment concentration (tons/m 3 ); v ch is mean flow velocity (m/s); y is mean flow depth (m); S is energy slope (m/m); Q in is water entering the reach (m 3 ); a, b, c and d are regression coefficients for different bed materials; W is channel width at the water level (m); W btm 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.
where C w is sediment concentration by weight; Ψ is universal stream power; M and N are coefficients. The universal stream power (Ψ) is calculated as follows.
where S g is relative density of solid (2.65); g is acceleration of gravity (9.81 m/s 2 ); depth is flow depth (m); ω 50 is fall velocity of median size particles (m/s); D 50 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 (D 50 ), unit stream power (V ch S), 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 (D 50 (9) where C w is weight concentration (ppm); ω 50 is fall velocity of the median size sediment (m/s); υ is kinematic viscosity (m 2 /s); V * is shear velocity (m/s); V ch is mean channel velocity (m/s); V cr is critical velocity (m/s); S is energy slope (m/m).

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].
where R eq is rainfall erosivity factor; P r 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 (R eq ). R = 11.8(Q·q·A) 0.56 (11) where R is runoff factor; Q is surface runoff (mm/ha); q is peak runoff (m 3 /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, Buildingornon − exposedground, C = 0.01 Barren, C = 1.0 (13)  (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.
where S is the slope factor of HRUs in MUSLE, θ is the slope of HRUs.

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 (R 2 = 0.79) (Equation (19)). ln(V) = 0.687 ln(A) + 2.326 (19) where V is estimated landslide volume (m 3 ); A is landslide area (m 2 ).

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 (R 2 = 0.79) (Equation (19)).
ln( ) = 0.687 ln( ) + 2.326 (19) where V is estimated landslide volume (m 3 ); A is landslide area (m 2 ). 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)).
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.

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

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., R 2 , NSE, PBIAS, and RSR), suggested by [60] with the model performance standards (Table 4). R 2 is the coefficient of determination, presenting the linear correlation between simulated and observed data. The value of R 2 closer to 1 indicates a higher correlation.
where Y sim is the simulated data; Y obs is observed data; Y mean 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].
Percent bias (PBIAS) presents whether the simulated data are overestimated or underestimated. When PBIAS is greater than 0, the simulation is underestimated [62].
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].

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.

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.

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.

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.

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.

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 (Tables 6 and 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,

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 (Tables 6 and 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 R 2 , NSE and smaller RSR, especially the SWAT-Twn performance had better statistical criteria values (R 2 = 0.74, NSE = 0.66, RSR = 0.58). Therefore, we used SWAT-Twn for further sediment calibration and validation.

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)(2005)(2006)(2007)(2008)(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).

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)(2005)(2006)(2007)(2008)(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. Table 8. The sensitivity analysis of sediment-related parameters (p-value < 0.05).
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 R 2 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.
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%.
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 R 2 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.

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

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][17][18][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. 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][17][18][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 subbasins. 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 subbasin area.

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 D 50 (>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].
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.

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 (R 2 and PBIAS) during 2004-2016 ( Figure 11 and Table 12). It should be noted that the simulation results are the same in Figures 9 and 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 R 2 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.

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 (R 2 and PBIAS) during 2004-2016 ( Figure 11 and Table 12). It should be noted that the simulation results are the same in Figures 9 and 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 R 2 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.

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.