Modeling Sediment Yields and Stream Stability Due to Sediment-Related Disaster in Shihmen Reservoir Watershed in Taiwan

Accurate and reliable estimates of sediment yields from a watershed and identification of unstable stream reaches due to sediment-related disaster are crucial for watershed management, disaster prevention, and hazard mitigation purposes. In this study, we added hydrodynamic and sediment transport modules in a recently developed model to estimate sediment yields and identify the unstable stream reaches in a large-scale watershed (> 100km2). The calibrated and verified models can well reproduce the flow discharge and sediment discharge at the study site, the Shihmen Reservoir Watershed in Taiwan, during several typhoon events. For the scenario applications, the results revealed that the contribution (> 96%) of landslides on sediment supply is much more significant than compared to soil erosion (< 4%). The sediment contribution from the upstream of the hydrological station-Yufeng is approximately 36–55% of the total sediment supply for the rainfall events of 25, 50, 100, and 200 years return period. It also indicates that 22–52% of sediment still remain at foot of the slope and the streams, which become a potential source for sediment hazards in the future. Combining with the bed erosion and deposition depths, flow-induced shear stress from the SRH-2D model, and probability of slope failure within 250 m of stream reaches, the relatively stability of stream reaches can be identified. The results could provide the water resource authorities for reference to take precautionary measures in advance on the stream reaches with high-degree instability.


Introduction
In tropical regions such as Taiwan, high intensity rainfall events from typhoons or plum rains often generate large amounts of sediments eroded from the upstream watershed [1].These sediments are then brought to the downstream reservoir areas, and the decrease of flow velocity causes sediment deposit, which would reduce the reservoir capacity and shorten the lifespan of reservoirs [2,3].Therefore, it is important to accurately estimate the sediment yields suffering intense rainfall events for the purposes of disaster prevention, hazard mitigation, and watershed management.Since sediment yield from a watershed is involved with numerous parameters such as local topography, soil, lithology, meteorological conditions, vegetation cover, land use, watershed morphology, and flow hydraulics [4,5], it is very complicated and highly variable in both time and hydromorphological location of watersheds.
Nowadays, three approaches are commonly used to estimate the sediment yield, which are field measurements, empirical-statistical models, and physically-based soil erosion-deposition models [6].For field measurements, sediment yield is estimated based upon the relationship between the suspended sediment concentration and flow discharge (sediment rating curve) established at hydrological stations.However, the sediment rating curve during the flood period do not normally collapse into a single-value relationship, rather, they often exhibit a hysteretic loop [7].Thus, the simplified relations between discharge and sediment yield cannot reflect temporal variations of sediment yield and generate significant estimate errors during the rising and falling stages of high flow events [8].Empirical-statistical models ignore the physical processes of sediment supply, transport, and delivery and use the data from laboratory or field measurements to obtain the regression equations to firstly estimate the soil erosion [9].Then, the soil erosion multiplies by the Sediment Delivery Ratio (SDR) to acquire the sediment yield.However, the SDR value is with high level of uncertainty and could become the major error source.The commonly-used empirical models are Universal Soil Loss Equation (USLE) [10], Modified Universal Soil Loss Equation (MUSLE) [11], and Revised Universal Soil Loss Equation (RUSLE) [12].Physically-based soil erosion-deposition models adopt physical equations such as EUROSEM [13], HSPF [14], KINEROS [15], LASCAM [16], SHESED [17], and SWAT [18].Nevertheless, the physically-based model needs to quantify numerous geographic and hydrological factors in a watershed.Though these models have provided acceptable estimates on sediment yield due to soil erosion, relatively few models include the landslide effects on total sediment yield.
Landslides play an important role in determining sediment yield from a watershed [19,20].Several studies attempted to establish the landslide component on total sediment yield [21][22][23][24][25][26].For example, Burton and Bathurst [21] developed a physically-based model-SHETRAN to estimate the landslide-related sediment yield.Schuerch et al. [22] monitored the sediment yield in a small catchment and compared the different sediment source supplied to the channels.Hsu et al. [23] integrated the HSPF, TRIGRS [24], and FLO-2D [25] models to estimate the sediment yield from soil erosion, shallow landslide, and debris flow.Chen et al. [26] proposed an integrated model of Grid-based Sediment Production and Transport Model (GSPTM) at a watershed scale to evaluate the dynamic of sediment production and transport in the Shihmen Reservoir watershed in Taiwan.In their model, rainfall-runoff processes, landslide potential, sediment production from landslide and soil erosion, debris flow and mass movement, and sediment transport were included.
This study aimed to estimate sediment yields and assess stream stability due to sedimentrelated disaster in a large-scale watershed (>100 km 2 ).An integrated model originally proposed by Chen et al. [26] was adopted, but two hydrodynamic and sediment transport models were added in the model to better quantify the flow field and calculate suspended sediment concentration, sediment yields and sediment erosion and deposition patterns in streambeds.To identify the unstable reaches in streams, an evaluation system based on the integrated model that can provide the bed erosion and deposition patterns, flow-induced shear stress on the bed and bank of streams, and possibility of slope failure near the streams was developed.The integrated model was then applied to the Shihmen Reservoir watershed to evaluate its performance through calibration and validation and demonstrate its feasibility and potential to estimate the sediment yields.Finally, to identify the response of sediment yields and stream stability with different rainfall events, scenarios with four rainfall return periods ranging from 25 years to 200 years were used in the developed model.

Study Area
Figure 1 shows the map and sub-watersheds of the study area-the Shihmen reservoir watershed.The Shihmen Reservoir watershed is located on the middle reaches of the Dahan River, one of the largest and most important reservoirs in Taiwan.The reservoir was built in 1964 for multiple purposes such as industrial and domestic water supply, irrigation, hydro-power generation, and flood prevention Water 2019, 11, 332 3 of 22 and recreation and had an original total capacity of 3.09 × 10 8 m 3 .The watershed area of the reservoir watershed is 760.2 km 2 with elevation ranging from 158 m (Shihmen Reservoir dam site) to 3524 m and average slope angles of about 30 degrees.The upper reach of the Shihmen Reservoir watershed has deep valleys and vulnerable geology characteristics, and, after the Chi-Chi earthquake occurred in 1999, several landslide areas appeared [27].In 2004, Typhoon Aere hit Taiwan with a total precipitation of 1600 mm in the Shihmen Reservoir watershed, which caused turbidity currents in the reservoir areas, brought approximately 2.8 × 10 7 m 3 sediments into the reservoir, and interrupted the water supply for 17 days [28,29].This typhoon event caused about 10% loss of the original reservoir storage capacity.According to bathymetry survey conducted in 2016, the deposited sediments occupied nearly 34% of the total storage capacity of the Shihmen Reservoir since 1964.To sustain the reservoir's life and ensure water supply, the mid-and long-term sediment desilting projects were proposed and implemented.The mid-term sediment management plan was to modify the existing sluicing outlets [30].The modification project was finished in 2013 and had increased the sediment desilting rate from about 30% to 45% under the typhoon Aere condition [30].For the long-term sediment management plan, sediment bypass tunnels-proposed as an important sediment remove measure,-will be built [31].In addition, Ronghua Dam (grey box in Figure 1), built in 1984, was within the Shihmen Reservoir Watershed and 27 kilometers upstream of the Shihmen Dam in the Dahan River gorge.Its main function is to prevent sediments from moving and silting in the Shimen Reservoir.The original storage capacity of the dam was 1.67 × 10 7 m 3 but has been filled by sediments at present [32].
Water 2018, 10, x FOR PEER REVIEW 3 of 22 watershed area of the reservoir watershed is 760.2 km 2 with elevation ranging from 158 m (Shihmen Reservoir dam site) to 3,524 m and average slope angles of about 30 degrees.The upper reach of the Shihmen Reservoir watershed has deep valleys and vulnerable geology characteristics, and, after the Chi-Chi earthquake occurred in 1999, several landslide areas appeared [27].In 2004, Typhoon Aere hit Taiwan with a total precipitation of 1,600 mm in the Shihmen Reservoir watershed, which caused turbidity currents in the reservoir areas, brought approximately 2.8 × 10  sediments into the reservoir, and interrupted the water supply for 17 days [28,29].This typhoon event caused about 10% loss of the original reservoir storage capacity.According to bathymetry survey conducted in 2016, the deposited sediments occupied nearly 34% of the total storage capacity of the Shihmen Reservoir since 1964.To sustain the reservoir's life and ensure water supply, the mid-and long-term sediment desilting projects were proposed and implemented.The mid-term sediment management plan was to modify the existing sluicing outlets [30].The modification project was finished in 2013 and had increased the sediment desilting rate from about 30% to 45% under the typhoon Aere condition [30].
For the long-term sediment management plan, sediment bypass tunnels-proposed as an important sediment remove measure,-will be built [31].In addition, Ronghua Dam (grey box in Figure 1), built in 1984, was within the Shihmen Reservoir Watershed and 27 kilometers upstream of the Shihmen Dam in the Dahan River gorge.Its main function is to prevent sediments from moving and silting in the Shimen Reservoir.The original storage capacity of the dam was 1.67 × 10  but has been filled by sediments at present [32].

Proposed Model
In this study, three individual models were integrated to estimate sediment yield from the Shihmen Reservoir watershed, which include: (i) Hydrological model, i.e., the rainfall-runoff models, to provide the upstream inflow discharge; (ii) hillside sediment production model to estimate the amount of sediments brought from the hillslope to the rivers; (iii) hydraulic and sediment transport models to obtain the suspended sediment concentration, flow-induced shear stress and bed erosion and deposition in the downstream rivers; and (iv) an evaluation system to assess the stability of stream reaches.Figure 2 illustrates the flowchart for the integrated model to estimate sediment yields and assess the stability of stream reaches in the study areas.The four individual models are briefly described below.In this study, three individual models were integrated to estimate sediment yield from the Shihmen Reservoir watershed, which include: (i) Hydrological model, i.e., the rainfall-runoff models, to provide the upstream inflow discharge; (ii) hillside sediment production model to estimate the amount of sediments brought from the hillslope to the rivers; (iii) hydraulic and sediment transport models to obtain the suspended sediment concentration, flow-induced shear stress and bed erosion and deposition in the downstream rivers; and (iv) an evaluation system to assess the stability of stream reaches.Figure 2 illustrates the flowchart for the integrated model to estimate sediment yields and assess the stability of stream reaches in the study areas.The four individual models are briefly described below.

Hydrological Model
The conceptual Tank model was adopted to simulate the rainfall-runoff processes, thus providing the behavior of rainwater in surface runoff, soil infiltration, and seepage [33][34][35].The Tank model describes the amount of water storage in a watershed with a series of tanks with the outlets on the side and bottom of the tanks [36].In the study, we referred to the three-layer Tank model developed by [37], where the total water depth can be divided into three layers, i.e., the surface ( ), intermediate ( ), and base layers ( ).The sum of the water depths in the three tanks is called the Soil Water Index (SWI), i.e.,  +  +  .The Tank model is suitable to simulate surface runoff, reservoir flood discharging and flood diversion [26].

Hydrological Model
The conceptual Tank model was adopted to simulate the rainfall-runoff processes, thus providing the behavior of rainwater in surface runoff, soil infiltration, and seepage [33][34][35].The Tank model describes the amount of water storage in a watershed with a series of tanks with the outlets on the side and bottom of the tanks [36].In the study, we referred to the three-layer Tank model developed by [37], where the total water depth can be divided into three layers, i.e., the surface (H 1 ), intermediate (H 2 ), and base layers (H 3 ).The sum of the water depths in the three tanks is called the Soil Water Index (SWI), i.e., H 1 + H 2 + H 3 .The Tank model is suitable to simulate surface runoff, reservoir flood discharging and flood diversion [26].

Hillside Sediment Production Model
The hillside sediment production model includes the sediment from landslide areas and soil erosion.The models based upon two different mechanisms for sediment production are described as follows.
(i) Landslide-induced sediment production: In this study, logistic regression (or binary regression analysis), a widely used statistical approach, was adopted to determine the landslide potential [38].The advantage of logistic regression is that the dependent variable y is binary, e.g., landslide occurrence (y = 1) or nonoccurrence (y = 0) in this study, and the independent (or called explanatory) variables are x, the environmental factors possibly to affect the landslide occurrence [36,39].The logistic regression model has been successfully applied to estimate the probability of landslide occurrence [38,40,41].Some detailed verification of applying the logistic regression model to landslide prediction can refer to literature [38,40,41].The logit model from a logistic regression has the following form: where y is the dependent variable, x i is the i-th explanatory variable, a is a constant, b i is the i-th regression coefficient, and e is the error term.The logit of y is the natural logarithm of the odds, which is: where p is the probability of the occurrence of y, and p/(1 − p) is the odds.Combining Equations ( 1) and ( 2), the probability p can be rewritten as: Herein, we adopted the logit model developed by [29] to identify the potential areas for landslide.Twelve related environmental factors (explanatory variables)-eleven were numerical and the other one was categorical-were adopted.The eleven numerical variables are elevation, slope, sin of aspect, cos of aspect, profile curvature, plan curvature, topographic wetness index (TWI), soil water index (SWI), distance to river, distance to ridge, and distance to road, respectively.The categorical variable is geological conditions.The geological data with a scale of 1:50000 published in 2000 was provided by the Central Geological Survey Bureau, Taiwan.There are twelve geological conditions in the Shihmen Reservoir watershed [26].The coefficients for the variables are determined from the landslide inventory during heavy rainfall events.For Shihmen Reservoir watershed, the logistic regression function is expressed as [26] logit(y) = −8 × By using Equation (3), the logit(y) can be converted to the probability p value of landslide occurrence.According to the results from previous literature [e.g., 26,38], the threshold of the p value set to 0.5 can obtain the best area concordance ratio (>70%), which is the overlapped landslide areas between the prediction and the total landslide area on the digital map divided by total landslide area on the digital map.Therefore, when the p value in this study was larger than 0.5, the cell was classified Water 2019, 11, 332 6 of 22 as a landslide cell.On the other hand, when the probability p of landslide occurrence was less than 0.5, the cell was regarded as stable.
Based upon the analysis of landslide potential, the location and area of landslide can be obtained.Secondly, the landslide area needs to be converted to the landslide volume.In this study, the landslide volume-area relationship was given by [42]: where V L is the landslide volume, A L is the landslide area, and a and r are the coefficients related to the geological properties for the selected watershed.In the Shihmen Reservoir Watershed, the landslide volume-area relationship was obtained as shown below [26]: where R 2 is the determination coefficient.
To estimate the sediment amount delivered to streams, the Sediment Delivery Ratio (SDR) between landslide runout distance and the landslide sediment brought to streams is required.In this study, the landslide runout distance L was given by Ikeya's formula [43]: where θ is the landslide slope (degrees) and D (m) is the distance between the centroid of the landslide area and the nearest downslope stream channel.Landslide migration distance L (m) represents the maximum possible moving distance of the sediment produced in a newly added landslide area.Once the SDR value was obtained, we could estimate the sediment amount from landslide into streams.
(ii) Sediment production from soil erosion: In addition to the landslide-derived sediment production, the soil loss from sheet and rill erosion also needs to be taken into account.In this study, MUSLE, an improved soil erosion over USLE and RUSLE, was adopted.In general, MUSLE can be expressed as follows [11]: where V s is sediment yield to the streams in metric tons, V e f f is the effective runoff volume in m 3 (the summation of the runoff obtained from the Tank model in each grid for a given rainfall event), Q p is the peak runoff flowrate in m 3 (the maximum value from the runoff hydrography obtained in the Tank model), K m is the soil erodibility factor, L and S are the slope length and gradient factor, C is the cover management, and P is the erosion control practice factor.In this study, the grid size derived from the digital elevation model (DEM) data provided by Center for Space and Remote Sensing Research, National Central University, Taiwan, was 40 m × 40 m, and the runoff and soil erosion were calculated in each grid and then summed up to the total runoff and soil erosion amount.The K m factor adopted ranged from 0.004 to 0.042 t-hr MJ −1 mm −1 [44,45].The L × S values are given by: where X is the slope length in meters; θ is the angle of the slope in degrees; and m is 0.5 (the slope is 5% or more), 0.4 (on slopes of 3% to 5%), 0.3 (on slopes of 1% to 3%), and 0.2 (on slopes less than 1%) [36].
To determine the slope length X, the DEM data were used to delineate the watershed topography.
In addition, a specified threshold needs to be assigned to distinguish between overland and channels.
Water 2019, 11, 332 As suggested by Liu et al. [46], a threshold of 100 m can be used as a default value to evaluate the slope length in most areas.The C factor was converted from a land-use map in [45].Assuming there is no protection measure, the P value is set to be 1 [44].

Flow and Sediment Transport Models in Streams
In this study, two hydrodynamic and sediment transport models were used.The first one was quasi two-dimensional (2-D) hydraulic and sediment routing models called "Network of Stream Tube Model for Alluvial River Simulation (NETSTARS)", which can provide the suspended sediment concentration and subsequently estimate sediment yields in the study areas.Though the NETSTARS model solves the 1-D flow and sediment governing equations, the NETSTARS applies the concept of stream tubes-imaginary tubes bounded by streamlines-to perform lateral variability in hydraulics and morphodynamics [47].However, the NETSTARS cannot simulate second currents, recirculating flow, or lateral erosion of stream banks.The second one: Sedimentation and River Hydraulics-Two-Dimensional River Flow Modeling (SRH-2D) model developed by Bureau of Reclamation, Department of the Interior, U.S. was applied to provide the spatial-temporal flow field and bed erosion and deposition patterns in detail for the downstream reach of the Shihmen Reservoir Watershed.The model solves the depth-averaged St. Venant equations and uses a flexible mesh that may contain arbitrarily shaped cells [48].In addition, the SRH-2D model adopts very robust and stable numerical schemes with a seamless wetting-drying algorithm.Therefore, the SRH-2D model can simultaneously model all flow regimes (sub-, super-, and trans-critical flows), both steady and unsteady flows, and handle flows over dry surfaces [49].Furthermore, the mobile-bed sediment transport module in the SRH-2D model can predict vertical stream bed changes by tracking multi-size, non-equilibrium sediment transport for suspended, mixed, and bed loads, and for cohesive and non-cohesive sediments [50].

Stability of Stream Reaches
In this study, three indicators associated with (i) erosion and deposition depth on channel beds, (ii) shear stress on channel wet boundary (channel bed and channel bank), and (iii) probability of slope failure were used to evaluate the stability in different stream reaches.The first two indicators can be obtained from the SRH-2D model, and the third one was calculated using the landslide model.The evaluation system including the three indicators for stream stability are described below.

(i) Erosion and deposition depth on channel beds:
According to simulation results from the SRH-2D model, the erosion and deposition depths on channel beds can be obtained.When the erosion or deposition depths range between −1 m (erosion) and 1 m (deposition), the stream is in the safest conditions.As the erosion or deposition depths become larger, the threat to stability of stream and adjacent hydraulic structures such as embankment, dam, groin, etc. is also increasing.Therefore, the indicator S d according to bed erosion or deposition depths that vary from 0 to 1 is given in Table 1 to reflect the threat level on stream stability.(ii) Flow-induced shear stress on channel wet boundary (channel bed and channel bank): The flow-induced shear stress τ w on the channel boundary can be used to evaluate the erosion potential to the streambed and stream bank.Herein, the maximum shear stress τ w,max during the intense rainfall events output from the SRH-2D model was divided into six classes, which were The corresponding value S τ for the flow-induced shear stress τ w,max to the six classes is listed in Table 2.When the slope near the riverbank fails, the streambed stability and adjacent hydraulic structures in streams may be affected or even damaged.Therefore, the probability of slope failure needs to be taken into account when we analyzed the stream stability.According to the landslide inventory collected in Typhoon Aere, 2004, the averaged landslide migration distance L in the Shihmen Reservoir watershed is approximately 250 m [51].It can be inferred that within 250 m on both sides of the stream bank, the slope failure may affect stability of streambed and hydraulic structures in streams.Therefore, the probability of slope failure was estimated using the logistic regression method in the 250 m range around streams.Similar to S d and S τ , the factor S p is dependent on the probability of slope failure in the 250 m range around the streams.Five classes of probability of slope failure p were proposed to determine the value of the factor S p , as listed in Table 3.The stability factor Stab can be obtained by simply taking the average of the above-mentioned three indicators in streams, which is Since there was no sufficient history data about the streambank failure or damage of river revetment, it is unlikely to give a specific threshold of the Stab value to determine if the stream would become unstable at the present stage.The magnitude of this value represents the relative stability of different regions of the stream course.The larger the Stab value, the more chances the stream reaches could become unstable during intense rainfall events.Therefore, the Stab value could help us identify the stream reaches with high chances to become unstable suffering heavy rain.The boundary conditions such upstream flow discharge, downstream water level, and sediment particle size of suspended and bed materials for numerical modelling were provided by the Water Resources Agency of Taiwan.The rainfall data were obtained from the Central Weather Bureau of Taiwan.The bathymetry was interpolated from a river cross-section survey (the interval between each cross-section is 500-600 m) conducted by the Northern Region Water Resource Office, the Water Resources Agency of Taiwan.In 2009, two river cross-section surveys were conducted on January 12 and August 29, respectively.For the case of January 12, 2009, the survey was from St.3 to St.8, while for the case of August 29, 2009, the survey was from St.2 to St.8, i.e., extension from St.2 to St.3.

Results and Discussion
In this section, the hydrodynamic and sediment transport models were firstly calibrated and verified to ensure their capability on providing accurate and reliable flow and sediment transport information.Typhoon Sinlaku and Morakot events were chosen for model calibration and verification.Typhoon Sinkalu stroke Taiwan in September 2008 with heavy rain and strong winds and produced 887 mm rainfall in total, the maximum cumulative rainfall for the past 15 years in the study areas.Typhoon Morakot swept Taiwan from August 7-9, 2009, carried record-breaking rainfall that the maximum five days cumulative rainfall reached 3,004.5 mm in Southern Taiwan, and triggered massive landslides in mountainous area.Though Typhoon Morakot only brought 510 mm rainfall in the study area, it is still valuable for understanding the spatiotemporal variations of flow discharge and sediment movement.The developed model was then applied for different rainfall scenarios to obtain flow and sediment discharge, and sediment erosion and deposition on streambeds.Finally, the stability of stream reaches were evaluated.

Sediment Yielding Prediction Using NETSTARS Model
This study firstly simulated the river discharge and sediment discharge during Typhoon Sinlaku in the Lofu Bridge (St.8 on Figure 1), which is regarded as the reservoir entrance, to calibrate the NETSTARS model.Then, the Typhoon Morakot event was used to verify the calibrated model.The following simulation conditions were used in the calibration and verification processes: (1) Simulation reach was between the Yufeng Station (St.2 on Figure 1) to the Lofu Bridge (St.8 on Figure 1) of the Dahan river.The river cross-section data input to the NETSTARS model was from the survey (from St.2 to St.8) conducted on 29 August 2009; (2) Upstream boundary conditions were inflow river discharge and sediment discharge, which were obtained by the tank model and the sediment production model, respectively; (3) Downstream boundary conditions were water levels measured in the Hsiayun station (St.7 on Figure 1).The Manning's n values were initially based on the mean sediment diameter d 50 and then adjusted to calibrate the developed model until an acceptable agreement between the measured and modelled flow and sediment discharges was reached.The calibrated Manning's n values range from 0.027-0.039.Figure 3 provides the calibration and verification results of sediment discharge in the Lofu Bridge during Typhoons Sinlaku and Morakot.The simulation and measurement data showed a similar trend, and the difference of the peak sediment discharge between the simulation and measurement data was approximately 7-15%.Considering the difficulties in the measurement and prediction of sediment transport, the model performance was reasonable and acceptable for engineering applications [52,53].

Streambed Eroison and Depostion Patterns Using the SRH-2D Model
In addition to the NETSTAR model, the SRH-2D model was used to simulate the erosion and deposition patterns in the downstream reaches of the Shihmen Reservoir Watershed, where the bathymetry data (from St.3 to St.8 on Figure 1) were interpolated from the river cross-section survey.The total length of simulated reach was approximately 28 km with an average slope of 0.013.The bathymetry for the simulated reach and the seven cross-sections used to compare with the measurements are shown in Figure 4.There were totally 14074 unstructured triangular grids in the computational domain, which is better fit to the geometry of the study reach.The boundary conditions about upstream inflow, lateral inflow from Shankuang Creek (St.4 on Figure 1), downstream water level, bed changes during the study period, and sediment particle size of suspended and bed materials were provided by the Water Resources Agency of Taiwan.The Manning's n values calibrated in the NETSTARS model was used in the SRH-2D model.The simulation time was from January 12, 2009 to August 29, 2009, covering the date of Typhoon Morakot passage and corresponding to the date of the two cross-section survey conducted.In the SRH-2D model, four sediment transport capacity equations could be selected.Based upon suggestion from Lai et al. [54], the equation proposed by Parker [55] was adopted in this study.In the equation (see Equation 10 in Lai et al. [54]), the hiding factor and reference non-dimensional shear stress (reference Shield's parameter) need to be determined in advance.The hiding factor accounts for the processes of non-uniform sediment movement, where the larger particles have higher chance exposed to the flow, and the required critical shear stress for sediment incipient motions is reduced [56].On the contrary, the smaller particles more likely to be covered by larger particles would be initiated by larger critical shear stress from the sediment beds.The reference non-dimensional shear stress was used to calculate the initiation of motion of sediment in a fluid flow [50].According to Lai et al. [54], the hiding factor and reference non-dimensional shear stress were set to 0.905 and 0.0385 for streams in Taiwan.

Streambed Eroison and Depostion Patterns Using the SRH-2D Model
In addition to the NETSTAR model, the SRH-2D model was used to simulate the erosion and deposition patterns in the downstream reaches of the Shihmen Reservoir Watershed, where the bathymetry data (from St.3 to St.8 on Figure 1) were interpolated from the river cross-section survey.The total length of simulated reach was approximately 28 km with an average slope of 0.013.The bathymetry for the simulated reach and the seven cross-sections used to compare with the measurements are shown in Figure 4.There were totally 14074 unstructured triangular grids in the computational domain, which is better fit to the geometry of the study reach.The boundary conditions about upstream inflow, lateral inflow from Shankuang Creek (St.4 on Figure 1), downstream water level, bed changes during the study period, and sediment particle size of suspended and bed materials were provided by the Water Resources Agency of Taiwan.The Manning's n values calibrated in the NETSTARS model was used in the SRH-2D model.The simulation time was from January 12, 2009 to August 29, 2009, covering the date of Typhoon Morakot passage and corresponding to the date of the two cross-section survey conducted.In the SRH-2D model, four sediment transport capacity equations could be selected.Based upon suggestion from Lai et al. [54], the equation proposed by Parker [55] was adopted in this study.In the equation (see Equation 10 in Lai et al. [54]), the hiding factor and reference non-dimensional shear stress (reference Shield's parameter) need to be determined in advance.The hiding factor accounts for the processes of non-uniform sediment movement, where the larger particles have higher chance exposed to the flow, and the required critical shear stress for sediment incipient motions is reduced [56].On the contrary, the smaller particles more likely to be covered by larger particles would be initiated by larger critical shear stress from the sediment beds.The reference non-dimensional shear stress was used to calculate the initiation of motion of sediment in a fluid flow [50].According to Lai et al. [54], the hiding factor and reference non-dimensional shear stress were set to 0.905 and 0.0385 for streams in Taiwan.when the overflow occurred at the Ronghua Dam (Figure 5a). Figure 5b shows that the SRH-2D model is also able to simulate and identify the flow field in the main channel and adjacent floodplains (white color).Figure 6 displays the flow field and the bed erosion and deposition patterns of the simulated reaches at 10 am on August 8, 2009, where the position values represent bed erosion and the negative values denote bed deposition.The maximum bed erosion and deposition depths can reach 6 m and 2 m, respectively, thus implying the significant redistribution of bottom sediments during typhoon events.Furthermore, the areas with greater flow velocity usually suffer to bed erosion.Figure 5 provides the flow field and contours of flow velocity at 10 am.On August 8, 2009 (during the Typhoon Morakot event).The results showed that the maximum velocity reached 20 ms −1 when the overflow occurred at the Ronghua Dam (Figure 5a). Figure 5b shows that the SRH-2D model is also able to simulate and identify the flow field in the main channel and adjacent floodplains (white color).Figure 6 displays the flow field and the bed erosion and deposition patterns of the simulated reaches at 10 am on August 8, 2009, where the position values represent bed erosion and the negative values denote bed deposition.The maximum bed erosion and deposition depths can reach 6 m and 2 m, respectively, thus implying the significant redistribution of bottom sediments during typhoon events.Furthermore, the areas with greater flow velocity usually suffer to bed erosion.
The comparisons between the measurement and simulation at seven cross sections are shown in Figure 7.In Section 1 (Figure 7a), the difference of the maximum erosion depth between the simulation and measurements is 1.5 m, and the simulations and measurements both show erosion on the right bank of the stream.For Sections 2 and 3 (Figure 7b,c), although the maximum difference between the simulation and measurements in bed elevation reached 2 m, the trend of sediment erosion and deposition between the measurement and simulation is still similar.In Section 4 (Figure 7d), the simulation and measurement both exhibit the obvious bed erosion at the right bank of the channel.The simulation and measurement are close with each other, indicating that the simulation can well reproduce the sediment transport processes.For Section 5 (Figure 7e), on the left bank of the channel, the simulation and measurement have the same trend on bed changes.However, the simulation displays approximately 10 m erosion on the left bank of the channel.For Section 6 (Figure 7f), the bed erosion and deposition between the simulations and measurements are similar.Due to the channel shrinkage at Section 7 (Figure 7g), some discrepancy between the simulation and measurement was found.The discrepancy between the numerical results and measurements can be attributed to several possible aspects: (1) The bathymetry in the computational domain was interpolated from river cross-section survey; (2) the two-dimensional depth-averaged model could not completely reflect the three-dimensional flow patterns in the study site; and (3) the temporal and spatial variations of the Manning's n values were not taken into account.Overall, the differences between the measurements and numerical results are still acceptable in terms of engineering applications.The SRH-2D model still can reproduce the bed erosion and depositional patterns close to the measurement in a general sense.The comparisons between the measurement and simulation at seven cross sections are shown in Figure 7.In Section 1 (Figure 7a), the difference of the maximum erosion depth between the simulation and measurements is 1.5 m, and the simulations and measurements both show erosion on the right bank of the stream.For Sections 2 and 3 (Figure 7b and c), although the maximum difference  The comparisons between the measurement and simulation at seven cross sections are shown in Figure 7.In Section 1 (Figure 7a), the difference of the maximum erosion depth between the simulation and measurements is 1.5 m, and the simulations and measurements both show erosion on the right bank of the stream.For Sections 2 and 3 (Figure 7b and c), although the maximum difference between the simulation and measurements in bed elevation reached 2 m, the trend of sediment interpolated from river cross-section survey; (2) the two-dimensional depth-averaged model could not completely reflect the three-dimensional flow patterns in the study site; and (3) the temporal and spatial variations of the Manning's n values were not taken into account.Overall, the differences between the measurements and numerical results are still acceptable in terms of engineering applications.The SRH-2D model still can reproduce the bed erosion and depositional patterns close to the measurement in a general sense.Water 2019, 11, 332 14 of 22

Sediment Yield Predictions Under Different Rainfall Return Periods
The developed model was used to estimate sediment yields of the Shihmen Reservoir Watershed under different rainfall scenarios.Herein, the rainfall data were designed as a 24-hour hyetograph with the return periods of 25, 50, 100, and 200 years, respectively.By using the logistic model with the given hyetograph, the probability p distribution of landslide occurrence can be obtained.Figure 8a provides the landslide probability map with the rainfall return period of 200 years.The results reveal that the landslide-prone areas are mostly located on concave slopes or near the stream courses.As mentioned above, the p value of 0.5 (50%) was set as a threshold, and the p value larger than 0.5 can be regarded as a landslide area.Therefore, the reactivated landslide areas can be estimated and provided in Table 4.It is shown that the Shihmen, Shankuang, and Kaoyi watersheds are on much lower levels of reactivated landslide areas than the rest of the watersheds.The possible explanations are as follows.For the Shihmen watershed, the slope-the most important factor for landslide-is gentlest among the seven sub-watersheds.Furthermore, the smaller elevation in the Shihmen watershed can also lead to lighter rainfall and lessen weathering effects on soil.For the Shankuang and Kaoyi watersheds, the slope and elevation are moderate among these watersheds.However, there are less developments of agriculture, business, and industry in the Shankuang and Kaoyi watersheds, thus indicating better conservation of soil and water and less chances for landslide occurrence.On the other hand, for the Hsiayun watershed, although the elevation and slope are also moderate, the region has become a popular site for tourism, mountain training, and camp centers.The newly-built road and associated facility for tourists may have caused the undercuts of slope, thus increasing the frequency of landslide occurrence [57].By applying Equation ( 6), the landslide areas can be converted to the landslide volume.The soil erosion from the hillslope was calculated using the MUSLE in Equation (9). Figure 8b,c provide the distribution of landslide-driven erosion and soil erosion in the Shihmen Reservoir Watershed.Table 5 lists the landslide-driven erosion and soil erosion in the Shihmen Reservoir Watershed.The erosion from landslide accounts for a large proportion (> 96%) of sediment production, while the soil erosion plays an insignificant role on sediment supply.According to the field measurement from Tsai et al. [58] and Cs 137 measurement from Chiu [59], the ratio of soil loss from landslide and watershed soil erosion was approximately 10:1, close to the results in this study.As the rainfall amount increases, landslides will occur more frequently, thus leading to the decreasing contribution of soil erosion to total sediment supply.
from landslide accounts for a large proportion (> 96%) of sediment production, while the soil erosion plays an insignificant role on sediment supply.According to the field measurement from Tsai et al. [58] and Cs 137 measurement from Chiu [59], the ratio of soil loss from landslide and watershed soil erosion was approximately 10:1, close to the results in this study.As the rainfall amount increases, landslides will occur more frequently, thus leading to the decreasing contribution of soil erosion to total sediment supply.According to the landslide and soil erosion data and applying Equation ( 8) to estimate the sediment delivered to streams, the flow and sediment discharges in the Shihmen Reservoir Watershed could be estimated from the NETSTARS model.For the 25-year, 50-year, 100-year, and 200-year return period rainfall events, the sediment discharge from the upstream of the Yufeng station was the greatest.Figure 9 displays the 24-hour 200-year return period design hyetograph, simulated river discharge, and sediment discharge at the upstream inlet (Yufeng station, St.2 on Figure 1) and downstream outlet (Lofu Bridge, St.8 on Figure 1).The sediment contribution from the upstream watershed of the Yufeng station was approximately 55% of the total sediment yield for the 200-year return period rainfall.As the rainfall amount decreased to the 25-year return period level, the ratio of the sediment brought from the Yufeng station to the total sediment supply also declined to 36%.The second greatest sediment contribution was from the upstream of the Shankuang station (St.3 on Figure 1), where the sediment contribution was about 7~9% to the total sediment supply in the Shihmen Reservoir Watershed.In addition, the simulation results also show that 48% (25-year return period rainfall)~78% (200-year return rainfall) of the sediment produced by landslide and hillslope soil erosion were brought to the reservoir areas, while 22-52% of sediment still remained at foot of the slope and the streams.The amount of sediment delivered to the downstream increased with the increasing rainfall and resulting flow discharge.However, there are still certain amounts of sediments left in the upstream watershed.If the next intense rainfall event occurs, the induced flow discharge could reactivate these remaining sediments, which could possibly lead to secondary sediment disasters such as turbidity current in the reservoir areas.
hillslope soil erosion were brought to the reservoir areas, while 22-52% of sediment still remained at foot of the slope and the streams.The amount of sediment delivered to the downstream increased with the increasing rainfall and resulting flow discharge.However, there are still certain amounts of sediments left in the upstream watershed.If the next intense rainfall event occurs, the induced flow discharge could reactivate these remaining sediments, which could possibly lead to secondary sediment disasters such as turbidity current in the reservoir areas.

Stream Stability under Different Rainfall Return Periods
According to the flow and sediment discharge from the three-layer Tank and sediment production models with the combination of the SRH-2D model, the erosion and deposition patterns in the downstream of the Shihmen Reservoir Watershed with different rainfall return periods can be obtained.Figure 10 shows the sediment erosion and deposition amount of the study streams for the 200-year return period rainfall.The results reveal that the amount of sediment deposited in streams is much more than the amount of sediment eroded from streams, indicating that the net amount of sediment is increased in streams.This corresponds to the conclusion drawn from the NETSTARS model that certain amounts of sediments are left in the upstream watershed.In addition, Figure 10 also displays that the depth of sediment deposition in most stream reaches is larger than 3 m.

Stream Stability under Different Rainfall Return Periods
According to the flow and sediment discharge from the three-layer Tank and sediment production models with the combination of the SRH-2D model, the erosion and deposition patterns in the downstream of the Shihmen Reservoir Watershed with different rainfall return periods can be obtained.Figure 10 shows the sediment erosion and deposition amount of the study streams for the 200-year return period rainfall.The results reveal that the amount of sediment deposited in streams is much more than the amount of sediment eroded from streams, indicating that the net amount of sediment is increased in streams.This corresponds to the conclusion drawn from the NETSTARS model that certain amounts of sediments are left in the upstream watershed.In addition, Figure 10 also displays that the depth of sediment deposition in most stream reaches is larger than 3 m. Figure 11 provides the flow-induce shear stress, probability of slope failure, and the stability factor  over the stream reaches.The results show that the maximum flow-induced shear stresses fall in a wide range from less than 500 N/m 2 to larger than 2000 N/m 2 (see Figure 11a).Larger shear stresses usually occur at the narrow streams or concave banks of streams.The probability of slope failure also show great spatial difference (Figure 12b).Comparing Figure 10, 11a and 11b, the distributions of the three indicators used for stream stability estimation are not identical.Only the maximum erosion depth is relevant to the maximum flow-induced shear stress.The three indicators can be regarded almost independent with each other.Therefore, it is reasonable to use the three indicators to evaluate the stream stability.According to the  values, the relative stream stability is shown in Figure 11c, where the higher values imply that the stream has greater probability to become unstable under intense rainfall events.The locations of buildings (purple color in Figure 11c) and stream revetments (grey color in Figure 11c), also provided in Figure 11c, reveal that some people live near the study stream but there are very few revetments constructed along the stream.If the stream reaches near the buildings, they could easily become unstable, and stream protection engineering such as revetments may be necessary.During heavy rainfall periods, the local residents living near the unstable reaches may need to be evacuated.Therefore, Figure 11c can also be a reference to the water resources authorities as to which stream reaches are priorities for improvement and management.Figure 11 provides the flow-induce shear stress, probability of slope failure, and the stability factor Stab over the stream reaches.The results show that the maximum flow-induced shear stresses fall in a wide range from less than 500 N/m 2 to larger than 2000 N/m 2 (see Figure 11a).Larger shear stresses usually occur at the narrow streams or concave banks of streams.The probability of slope failure also show great spatial difference (Figure 11b).Comparing Figure 10, Figure 11a,b, the distributions of the three indicators used for stream stability estimation are not identical.Only the maximum erosion depth is relevant to the maximum flow-induced shear stress.The three indicators can be regarded almost independent with each other.Therefore, it is reasonable to use the three indicators to evaluate the stream stability.According to the Stab values, the relative stream stability is shown in Figure 11c, where the higher values imply that the stream has greater probability to become unstable under intense rainfall events.The locations of buildings (purple color in Figure 11c) and stream revetments (grey color in Figure 11c), also provided in Figure 11c, reveal that some people live near the study stream but there are very few revetments constructed along the stream.If the stream reaches near the buildings, they could easily become unstable, and stream protection engineering such as revetments may be necessary.During heavy rainfall periods, the local residents living near the unstable reaches may need to be evacuated.Therefore, Figure 11c can also be a reference to the water resources authorities as to which stream reaches are priorities for improvement and management.

Conclusions
In this study, hydrodynamic and sediment transport modules were added in an integrated model to estimate sediment yields and identify the possibly unstable stream reaches in the Shihmen Reservoir Watershed, a large-scale watershed (>100 km 2 ) in Taiwan.Through calibration and verification, the hydrodynamic and sediment transport model successfully reproduced river flow patterns, flow discharge, sediment discharge, and streambed changes during typhoon events.Therefore, the developed model is feasible and capable to provide the flow-induced boundary shear stress and accurate estimates of sediments yield caused by two mechanisms such as landslide and soil erosion and flow patterns in a large-scale watershed (>100 km 2 ).Combining with the bed erosion and deposition depths, flow-induced shear stress, and probability of slope failure within 250 m of stream reaches, an evaluation system was proposed to assess the relative stream stability.
The model was then applied to different scenarios such as rainfall events of 25, 50, 100, and 200 years return period in the study areas.The results show that the sediment induced by landslide accounts for more than 96% of sediment supply in the Shihmen Reservoir Watershed.The sediment contribution from the upstream of the Yufeng station is approximately 36-55% of the total sediment yield for different scenarios.They also indicate that 22-52% of sediment still remain at foot of the slope and the streams, which will become a potential source for hazard in the future.Using the evaluation system, the relatively stability of stream reaches can be identified, which could provide the water resource authorities for reference to take precautionary measures on the stream reaches with high-degree instability before the hazards occur.

Conclusions
In this study, hydrodynamic and sediment transport modules were added in an integrated model to estimate sediment yields and identify the possibly unstable stream reaches in the Shihmen Reservoir Watershed, a large-scale watershed (>100 km 2 ) in Taiwan.Through calibration and verification, the hydrodynamic and sediment transport model successfully reproduced river flow patterns, flow discharge, sediment discharge, and streambed changes during typhoon events.Therefore, the developed model is feasible and capable to provide the flow-induced boundary shear stress and accurate estimates of sediments yield caused by two mechanisms such as landslide and soil erosion and flow patterns in a large-scale watershed (>100 km 2 ).Combining with the bed erosion and deposition depths, flow-induced shear stress, and probability of slope failure within 250 m of stream reaches, an evaluation system was proposed to assess the relative stream stability.
The model was then applied to different scenarios such as rainfall events of 25, 50, 100, and 200 years return period in the study areas.The results show that the sediment induced by landslide accounts for more than 96% of sediment supply in the Shihmen Reservoir Watershed.The sediment contribution from the upstream of the Yufeng station is approximately 36-55% of the total sediment yield for different scenarios.They also indicate that 22-52% of sediment still remain at foot of the slope and the streams, which will become a potential source for hazard in the future.Using the evaluation system, the relatively stability of stream reaches can be identified, which could provide the water resource authorities for reference to take precautionary measures on the stream reaches with high-degree instability before the hazards occur.

Figure 1 .
Figure 1.Sub-watersheds of the study area-Shihmen reservoir watershed, Taiwan.X(m) and Y(m) are based upon the TWD97/TM2 coordinate system.The inset figure shows the longitude and latitude of the study areas.

Figure 1 .
Figure 1.Sub-watersheds of the study area-Shihmen reservoir watershed, Taiwan.X(m and Y(m) are based upon the TWD97/TM2 coordinate system.The inset figure shows the longitude and latitude of the study areas.

Figure 2 .
Figure 2. Framework for the proposed integrated model.

Figure 2 .
Figure 2. Framework for the proposed integrated model.

Water 2019 ,Figure 3 .
Figure 3. Calibration and verification results for the sediment discharge.

Figure 3 .
Figure 3. Calibration and verification results for the sediment discharge.

Figure 4 .
Figure 4.The bathymetry and seven cross-sections used to compare with the measurements for the Sedimentation and River Hydraulics-Two-Dimensional River Flow Modeling (SRH-2D) model.The coordinate is based upon the TWD97/TM2 system.

Figure 5
Figure 5 provides the flow field and contours of flow velocity at 10 am.On August 8, 2009 (during the Typhoon Morakot event).The results showed that the maximum velocity reached 20 when the overflow occurred at the Ronghua Dam (Figure5a).Figure5bshows that the SRH-2D model is also able to simulate and identify the flow field in the main channel and adjacent floodplains (white color).Figure6displays the flow field and the bed erosion and deposition patterns of the simulated reaches at 10 am onAugust 8, 2009, where the position values represent bed erosion and the negative values denote bed deposition.The maximum bed erosion and deposition depths can reach 6 m and 2 m, respectively, thus implying the significant redistribution of bottom sediments during typhoon events.Furthermore, the areas with greater flow velocity usually suffer to bed erosion.

Figure 4 .
Figure 4.The bathymetry and seven cross-sections used to compare with the measurements for the Sedimentation and River Hydraulics-Two-Dimensional River Flow Modeling (SRH-2D) model.The coordinate is based upon the TWD97/TM2 system.

Water 2018 , 22 Figure 5 .
Figure 5. Flow field and contours of flow velocity at 10 am on August 8, 2009 (during the Typhoon Morakot event).

Figure 6 .
Figure 6.Flow field and erosion and deposition patterns at 10 am on August 8, 2009 (during the Typhoon Morakot event).

Figure 5 . 22 Figure 5 .
Figure 5. Flow field and contours of flow velocity at 10 am on August 8, 2009 (during the Typhoon Morakot event).

Figure 6 .
Figure 6.Flow field and erosion and deposition patterns at 10 am on August 8, 2009 (during the Typhoon Morakot event).

Figure 6 .
Figure 6.Flow field and erosion and deposition patterns at 10 am on August 8, 2009 (during the Typhoon Morakot event).

Figure 7 .
Figure 7. Comparisons between the measurement and simulation at seven cross sections.Note: the black line ' ' indicates the initial bed location (12 January 2009); the grey line ' ' is the measured bed location (29 August 2009); the dash line ' ' is the simulated bed location (29 August 2009).

22 (g) Section 7 Figure 7 .
Figure 7. Comparisons between the measurement and simulation at seven cross sections.Note: the black line ' ' indicates the initial bed location (12 January 2009); the grey line ' ' is the measured bed location (29 August 2009); the dash line ' ' is the simulated bed location (29 August 2009).

22 (g) Section 7 Figure 7 .
Figure 7. Comparisons between the measurement and simulation at seven cross sections.Note: the black line ' ' indicates the initial bed location (12 January 2009); the grey line ' ' is the measured bed location (29 August 2009); the dash line ' ' is the simulated bed location (29 August 2009).
left bank (m) is the simulated bed location (29 August 2009).

Figure 8 .
Figure 8. Results for the return periods of 200 years in the Shihmen Reservoir Watershed.

Figure 8 .
Figure 8. Results for the return periods of 200 years in the Shihmen Reservoir Watershed.

Figure 9 .
Figure 9. Result for the return periods of 200 years in the Shihmen Reservoir Watershed.The blue solid and dot lines represent the flow discharge at the downstream (Lofu Bridge) and upstream (Yufeng Station), and the brown solid and dot lines denote the sediment discharge at the downstream (Lofu Bridge) and upstream (Yufeng Station).

Figure 9 .
Figure 9. Result for the return periods of 200 years in the Shihmen Reservoir Watershed.The blue solid and dot lines represent the flow discharge at the downstream (Lofu Bridge) and upstream (Yufeng Station), and the brown solid and dot lines denote the sediment discharge at the downstream (Lofu Bridge) and upstream (Yufeng Station).

Figure 10 .
Figure 10.The sediment erosion (a) and deposition (b) amount of the study stream for the 200-year return period rainfall.

Figure 10 .
Figure 10.The sediment erosion (a) and deposition (b) amount of the study stream for the 200-year return period rainfall.

Figure 11 .
Figure 11.Flow-induced shear stress (a), probability of slope failure (b), and stream stability factor (c) of the study stream for the 200-year return period rainfall.

Figure 11 .
Figure 11.Flow-induced shear stress (a), probability of slope failure (b), and stream stability factor (c) of the study stream for the 200-year return period rainfall.

Table 1 .
S d for erosion and deposition depth on streambeds.

Table 2 .
S τ for flow-induced shear stress on streambed boundary

Table 3 .
S p for probability of slope failure near the streambank.

Probability of Slope Failure S p Value Possibility to Damage the River Structures
In this study, the digital elevation model (DEM) data with the grid size of 40 m × 40 m was provided by Center for Space and Remote Sensing Research, National Central University, Taiwan.

Table 4 .
Predicted reactivated landslide areas under different rainfall return periods.

Table 5 .
Predicted sediment production under different rainfall return periods.

Table 5 .
Predicted sediment production under different rainfall return periods.