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

**Figure 1.** Study area.


**Table 1.** Top 10 typhoons that influenced the study area during 2004–2015.

(Source: Typhoon Database, Central Weather Bureau, Taiwan [31]).

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

**Figure 2.** Locations of weather stations, sub-basins, and gaging station.

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.

**Figure 3.** Land use, soil, and slope distribution.

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

$$\mathbf{Q\_{surf}} = \frac{\left(\mathbf{R\_{day}} - \mathbf{I\_a}\right)^2}{\left(\mathbf{R\_{day}} - \mathbf{I\_a}\right) + \mathbf{S}} \tag{1}$$

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.

**Figure 4.** Sub-basin distribution.

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

$$\text{S}\_{\text{t}} = \text{S}\_{\text{0}} + \sum\_{i=1}^{t} \text{R}\_{\text{day}} - \text{Q}\_{\text{surf}} - \text{E}\_{\text{a}} - \text{W}\_{\text{secp}} - \text{Q}\_{\text{gw}} \tag{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.

$$\mathbf{S} = 11.8 \times (\mathbf{Q\_{surf}} \times \mathbf{q} \times \mathbf{A})^{0.56} \times \mathbf{K} \times \mathbf{C} \times \mathbf{P} \times \mathbf{LS} \times \mathbf{CFRG} \tag{3}$$

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.

$$\mathbf{c} \mathbf{u} \mathbf{c}^{\text{seq,ch,mx}} = \mathbf{c}\_{\text{sp}} \mathbf{v}\_{\text{ch,pk}} \mathbf{x}^{\text{spec}} \tag{4}$$

where concsed,ch,mx is the maximum sediment concentration (t/m3); csp is linear coefficient; vch,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.

$$\text{Conc}\_{\text{scd,ch,mx}} = \left(\frac{\text{a} \cdot \text{v}\_{\text{ch}} \text{}^{\text{b}} \cdot \text{y}^{\text{c}} \cdot \text{S}^{\text{d}}}{\text{Q}\_{\text{in}}}\right) \cdot \left(\frac{\text{W} + \text{W}\_{\text{btm}}}{2}\right) \tag{5}$$

where concsed,ch,mx is the maximum sediment concentration (tons/m3); vch 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.

$$\mathbf{C}\_{\mathbf{W}} = \mathbf{M} \mathbf{Y}^{\mathbf{W}} \tag{6}$$

where Cw is sediment concentration by weight; Ψ is universal stream power; M and N are coefficients. The universal stream power (Ψ) is calculated as follows.

$$\Psi = \left\{ \left( \mathbf{S}\_{\mathbf{g}} - 1 \right) \mathbf{g} \cdot \mathbf{depth} \cdot \boldsymbol{w}\_{50} \cdot \left[ \log\_{10} \left( \frac{\text{depth}}{\text{D}\_{50}} \right) \right] \right\}^{0.5} \tag{7}$$

where Sg is relative density of solid (2.65); g is acceleration of gravity (9.81 m/s2); depth is flow depth (m); ω<sup>50</sup> 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,

$$\begin{array}{l} \log \mathcal{C}\_{\text{W}} \ = 5.435 - 0.286 \log \frac{\omega \gamma\_{50} \mathcal{D}\_{\text{N}}}{\upsilon} - 0.457 \log \frac{\mathcal{V}\_{\text{v}}}{\omega \gamma\_{50}} \\ \quad + \left(1.799 - 0.409 \log \frac{\omega \gamma\_{50} \mathcal{D}\_{\text{N}}}{\upsilon} - 0.314 \log \frac{\mathcal{V}\_{\text{v}}}{\omega \gamma\_{50}}\right) \log \left(\frac{\mathcal{V}\_{\text{ch}} \mathcal{S}}{\omega \varepsilon\_{50}} - \frac{\mathcal{V}\_{\text{c}} \mathcal{S}}{\omega \varepsilon\_{50}}\right) \end{array} \tag{8}$$

Gravel equation for median size (D50) between 2 mm and 10 mm,

$$\begin{array}{ll} \log \mathcal{C}\_{\text{W}} &= 6.681 - 0.633 \log \frac{\omega\_{\text{50}} \mathcal{D}\_{\text{50}}}{\nu} - 4.816 \log \frac{\mathcal{V}\_{\text{\*}}}{\omega\_{\text{50}}} \\ &+ \left(2.784 - 0.305 \log \frac{\omega\_{\text{50}} \mathcal{D}\_{\text{50}}}{\nu} - 0.282 \log \frac{\mathcal{V}\_{\text{\*}}}{\omega\_{\text{50}}}\right) \log \left(\frac{\mathcal{V}\_{\text{ch}} \mathcal{S}}{\omega\_{\text{50}}} - \frac{\mathcal{V}\_{\text{c}} \mathcal{S}}{\omega\_{\text{50}}}\right) \end{array} \tag{9}$$

where Cw is weight concentration (ppm); ω<sup>50</sup> is fall velocity of the median size sediment (m/s); υ is kinematic viscosity (m2/s); <sup>V</sup><sup>∗</sup> is shear velocity (m/s); Vch is mean channel velocity (m/s); Vcr is critical velocity (m/s); S is energy slope (m/m).
