Evaluation of the SEdiment Delivery Distributed (SEDD) Model in the Shihmen Reservoir Watershed

: The sediment delivery ratio (SDR) connects the weight of sediments eroded and transported from slopes of a watershed to the weight that eventually enters streams and rivers ending at the watershed outlet. For watershed management agencies, the estimation of annual sediment yield (SY) and the sediment delivery has been a top priority due to the inﬂuence that sedimentation has on the holding capacity of reservoirs and the annual economic cost of sediment-related disasters. This study establishes the SEdiment Delivery Distributed (SEDD) model for the Shihmen Reservoir watershed using watershed-wide SDR w and determines the geospatial distribution of individual SDR i and SY in its sub-watersheds. Furthermore, this research considers the statistical and geospatial distribution of SDR i across the two discretizations of sub-watersheds in the study area. It shows the probability density function (PDF) of the SDR i . The watershed-speciﬁc coe ﬃ cient ( β ) of SDR i is 0.00515 for the Shihmen Reservoir watershed using the recursive method. The SY mean of the entire watershed was determined to be 42.08 t / ha / year. Moreover, maps of the mean SY by 25 and 93 sub-watersheds were proposed for watershed prioritization for future research and remedial works. The outcomes of this study can ameliorate future watershed remediation planning and sediment control by the implementation of geospatial SDR w / SDR i and the inclusion of the sub-watershed prioritization in decision-making. Finally, it is essential to note that the sediment yield modeling can be improved by increased on-site validation and the use of aerial photogrammetry to deliver more updated data to better understand the ﬁeld situations.


Introduction
Land degradation has been confirmed as a threat to agricultural productivity worldwide [1]. The gradual dissipation of quality topsoil transported as sediments along hill slopes and deposited into rivers and reservoirs has impacted multiple agricultural areas across the world. Sediment delivery by soil erosion and landslides has also led to the decreased longevity and failure of reservoirs and other water infrastructure. Moreover, due to the importance of reservoirs and dams to the water security of large populations in the developed world, the impact of sedimentation has been highlighted as a major issue in need of research worldwide over the past half a century and into the future. These processes influence many of the United Nations' (UN) Sustainable Development Goals (SDGs) [1][2][3]. Heavy sustained sedimentation leads to a severe influence on industries and economic growth for agriculture-based communities (SDG 8,9,10), the sustainability of cities and communities (SDG 11), and the availability of clean water downstream (SDG 6). These processes also pollute rivers, streams, and the sea by the transport of chemicals with sediments damaging marine flora and fauna and influencing the sustainability of marine environments (SDG 14).
Southeast Asia presents an even more advanced dilemma as the environmental and geomorphological characteristics across the region displays a strong influence on the rate of soil erosion, and the region's population and economies are already suffering from the impact [1]. After Typhoon Gloria impacted Taiwan, the Shihmen reservoir watershed was one of the many watersheds heavily inundated by sediments from soil erosion and landslides. This event was an example of typhoon events affecting Taiwan due to its tropical monsoon climate. The impact of these events creates frequent and high magnitude sediment disasters and long-term sediment inundation. The effect of extreme meteorological events on Shihmen and other reservoirs and the threat to urban and agricultural water supply quality throughout the nation has encouraged the Government of Taiwan to facilitate soil conservation research and countermeasures within the watershed through the introduction of soil conservation policies and programs. Taiwan's geographical location places it in a precarious position that geomorphologically has a high frequency of highly erodible soil materials and, additionally, its tropical monsoon climate brings frequent extreme rainfall or typhoon events which have led to massive sediment loads, dissociated by mass movements, and soil erosion, to be transported from slopes of watersheds into rivers and streams, affecting water supply in many communities. These events increase the siltation of reservoirs, such as the Shihmen Reservoir in Northern Taiwan, which diminishes the country's capacity to sustainably manage its water supply and poses a threat to the long-term sustainability of the dam and the populations it serves [4,5].
Therefore, sediment modeling and monitoring in Taiwan is of paramount importance to improve watershed conservation and sustainable water management by understanding how sediments are dissociated, translocated, and transported down the slopes of watersheds, entering streams as massive sediment loads. These sediments deposit along the waterway with detrimental impacts to river flow and along slopes, increasing future sediment disaster hazards and clogging reservoirs and dams, affecting water supplies to nearby populated districts and farms.

Sediment Delivery Ratio (SDR)
Researchers have established that the sediment yield (SY) and siltation found in rivers and dams are not equivalent to the gross sediments that have dissociated from the slopes through land degradation processes. The sediment delivery ratio (SDR) was established as the relationship between the SY and the gross erosion (A, in some research denoted as Soil Loss, SL). The SDR acts as a reducer to equate the maximum amount of soil erosion to the sediments delivered to the outlet. SDR is a continuous variable with a range from zero to one. In the initial stages of soil erosion research, SDR was evaluated for an entire watershed and a single outlet, such as a reservoir, denoted as SDR w . However, as computational power has improved and research methods have become more developed, SDR is discretized to sub-watersheds, micro watersheds, and even more minute hydrological units within a watershed. This has been denoted as a geospatial SDR i and is defined as the fraction of the gross soil loss (SL) from the hydrological unit, i, that reaches a flow pathway [6].
Many methods for determining the SDR have been developed over the past six decades. Long-term monitoring of SY in test plots and different regions were the starting point in SDR research and established multiple geophysical and hydrological factors that affect SDR. Wu et al. [7] stated that there are three prominent methods of estimating SDR, which are: (1) the ratio between soil erosion and sediment yield, (2) empirical formula relationships between SDR and its contributing factors, and (3) spatial modeling and numerical modeling using hydraulic and hydrological theories. These include models for the calculation and classification of SDR based on regional and local observed data.
Li [8] determined from the comparison of the USLE results and the outlet sediments at the Shihmen Reservoir that the SDR was 0.49 from the annual siltation amount of 71.2 t/ha/year (between 2011 and 2015), which includes sheet and rill erosion, gully erosion, and landslides. There have been other studies of the SDR in the Shihmen Reservoir watershed but these have mainly focused on the post-landslide erosion SDR or event-based SDR. Tsai et al. [9] defined the SDR of the landslide areas in the watershed from 0.42 to 0.78 for eroded areas and zero for areas where aggradation occurred. Moreover, Tsai et al. [10] employed the model of Bathurst et al. [11] of the sediment delivery of landslides. They determined that the sediment delivery of landslides in the Shihmen Reservoir watershed ranged from 0.75 to 0.89 (1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998), 0.48 to 0.64 (1998)(1999)(2000)(2001)(2002)(2003)(2004), and 0.40 to 0.58 for typhoon Aere (2004). Lin et al. [12] explored the impact of different factors (such as curve number and Manning's) on the SDR in the Shihmen Reservoir watershed.
SWCB [13] previously employed DEMs of difference (DoD) method to determine a correlation between SDR and the entire watershed (SDR w ) in the Shihmen Reservoir watershed and concluded the relationship to be SDR = 41.23W −0.017 Area where W Area is the watershed area, and the SDR w was 0.368. Chen et al. [14] used the Grid-based Sediment Production and Transport Model (GSPTM) to determine SDR w to be 0.299 (rounded to 0.30) for a simulation of the impact of the Typhoon Morakot on the Yufeng and Xiuluan sub-watersheds, and Chiu et al. [15] extended the GSPTM to identify unstable stream reaches. Additionally, Liu [16] determined that the total soil loss (also commonly-termed gross soil erosion or the total amount of soil erosion) was 90.6 t/ha/year by converting measured average erosion depth of erosion pins [17,18] to the amount of soil erosion. This study emphasizes the SDR of recent years and assumes SDR w = 0.49.

Spatial Discretization
This research compares the implementation of the SEDD with grid-cell unit analysis across the watershed and two discretization levels of sub-watersheds (25 and 93 SW). Spatial discretization of watersheds in hydrological models has been a critical issue for many researchers within the remote sensing, physical modeling, hydrology, and mathematical computations fields [19]. Wood et al. [20] elaborated that watersheds are areas made of infinite minute points of related hydrological processes (for example, evaporation, infiltration, and runoff). Watersheds can be further divided into sub-watersheds or sub-basins that are continuous assemblages or spatial objects that drain towards a specific point, contains similar vegetation layers and/or linked slope, or elevation characteristics, and are commonly used in the study of landslides, soil erosion, hydrology, and other natural processes and are specific areas within a watershed. There are many methods of defining watersheds, sub-watersheds, and even more minute forms of discretization, such as hydrological units depending on the characteristics of the areas [19][20][21][22].
Contextually, at the time of the development of many soil erosion models and concepts from the 1970s to the 1990s, computational cost, time, and power were limitations that hindered the discretization of watersheds to more elements [19]. Therefore, many older models were developed with sub-watershed as the smallest hydrological unit. This has continued even to the 2000s with the SWAT (Soil and Water Assessment Tool) model, and research has supported the object area discretization in soil erosion, LULC (Land Use and Land Cover), and sedimentation models and has criticized the use of areal discretization. SWAT is a comprehensive soil erosion model that is popular in the field of hydrology. In this study, we utilized the first stage in the SWAT analysis, the discretization of a study area into sub-basins/sub-watersheds, to develop 25-subwatershed divisions. The two inputs used in this analysis was the surface topography and the land use type [23]. Previously in the Shihmen Reservoir watershed, grid cells and slope units have also been employed in soil erosion research [16].

Study Objectives
Geomorphological changes, climate change, and land degradation processes are highly influential to soil erosion studies. Therefore, it is important to highlight the time scale. The prioritization of watersheds or sub-watersheds for soil erosion, landslide, and conservation works has been a long-standing topic in geomorphology, hydrology, disaster management, and environmental planning research. In the past, the gross soil loss (often derived by the USLE model) and its contributing factors have been the primary focus for prioritization of soil conservation planning. This study explored the sediment yield value to the prioritization of sub-watersheds for remediation and conservation.
This research incorporated and considered the implementation of the SEDD model in the Shihmen Reservoir watershed and the research objectives, methods, and assumptions are as follows: 1.
The main focus of this study was on the advancement of the soil erosion research in the Shihmen Reservoir watershed to include determination of the watershed-specific coefficient, β, using known SDR w values, and developing the SDR i relationship and SEDD model for the watershed, applying the findings to the assessment of the sediment yield (SY) by sheet and rill erosion in the watershed.

2.
In this study, the non-point source sheet and rill erosion calculated by the USLE is the main target land degradation process. Landslides, gully erosion, and channel type erosions are not considered.

3.
The datasets employed are assumed to be from the same period and are the most recent available data for this study area.

4.
This study, additionally, explores the statistical properties of the SDR, SL, and SY of the watershed using two sub-watershed discretizations (25 and 93) of the Shihmen Reservoir watershed and compares the results.

5.
This study determines the watershed prioritization of soil conservation using the SY determined by the SEDD.
This study is a theoretical and statistical analysis of the SDR distribution based upon the SEDD model using real data from the Shihmen Reservoir watershed.

Study Area
The Shihmen Reservoir watershed (elevation map shown in Figure 1a) is a 759.53 km 2 watershed found in Northern Taiwan between latitudes 24.426 • N and 24.861 • N, and longitudes 121.192 • E and 121.479 • E. Local authorities use the reservoir to regulate 23% of the water supply for the nearby communities, industries, and agriculture. The effective storage volume of the reservoir is 207 million m 3 . Still, it has been affected by heavy sedimentation, and some of the additional check dams established to diminish sedimentation have been affected or damaged in past years. The established subdivision of the watershed by local government is into five sub-watersheds, namely Ku-Chu, Yu-Feng, San-Kuang, Pai-Shih, Tai-Kang sub-watersheds. However, in this study, the entire watershed is delineated into 25 sub-watersheds ( Figure 1b) and 93 sub-watersheds (Figure 1c), as explained in Section 3.1. The watershed is traversed by the Tahan River and its irregular pattern of tributaries that spread across the steep, high slopes of the sub-watersheds (angle 20 • to 85 • with an elevation between 220 m and 3527 m). The subtropical monsoon climate and the annual rainfall between 2200 mm and 2800 mm in the watershed encourages a densely vegetated state, which is mostly natural and artificial coniferous and broadleaved trees, and the trend in rainfall pattern has increased over the past decades to create concerns of soil erosion (Figure 1d) [9,31].

Methodology
The SEDD model developed by Ferro and Minacapilli [24] analyzes the spatially distributed SY by determining the parameters of SDR and SL found using any soil loss model. The SDR is found by a relationship between the travel time (tp,i) and a watershed-specific coefficient, β. β cannot be determined analytically but is determined iteratively, by balancing the relationship between a known SY, SL, and SDR. The SY is evaluated for the entire watershed and the two discretizations, 25 and 93 subwatershed models, to compare the impact of aggregation on the model outputs.
The flow chart of the SEDD model implemented in this study is distributed among four stages as shown in Figure

Methodology
The SEDD model developed by Ferro and Minacapilli [24] analyzes the spatially distributed SY by determining the parameters of SDR and SL found using any soil loss model. The SDR is found by a relationship between the travel time (t p,i ) and a watershed-specific coefficient, β. β cannot be determined analytically but is determined iteratively, by balancing the relationship between a known SY, SL, and SDR. The SY is evaluated for the entire watershed and the two discretizations, 25 and 93 sub-watershed models, to compare the impact of aggregation on the model outputs.
The flow chart of the SEDD model implemented in this study is distributed among four stages as shown in Figure 2: Determine the study area and input the DEM.

2.
Determine the slope, flow direction, and river/stream map from a derived depression-less DEM.

3.
Use the flow direction and river/stream map as input to determine the flow distance, which is then used with the slope to determine the travel time.

4.
Find the six parameters of the USLE equation and estimate SL. 5.
Determine the SDR using measured outlet SY and estimated SL.
Back to stage 1.2:

6.
Determine β from the SDR, SL, and the travel time found.

7.
Return to the DEM and discretize the study area into sub-watersheds/hydrological units or any other applicable unit of analysis.
Re-estimate the SDR for the parent watershed and each sub-watershed. 9.
Estimate the sediment yield (SY) and perform any further analysis necessary.
Sustainability 2020, 12, x FOR PEER REVIEW 6 of 20 1. Determine the study area and input the DEM. 2. Determine the slope, flow direction, and river/stream map from a derived depression-less DEM. 3. Use the flow direction and river/stream map as input to determine the flow distance, which is then used with the slope to determine the travel time.
Stage 2: 4. Find the six parameters of the USLE equation and estimate SL. 5. Determine the SDR using measured outlet SY and estimated SL.
Back to stage 1.2: 6. Determine β from the SDR, SL, and the travel time found.
Stage 3: 7. Return to the DEM and discretize the study area into sub-watersheds/hydrological units or any other applicable unit of analysis.
Stage 4: 8. Re-estimate the SDR for the parent watershed and each sub-watershed. 9. Estimate the sediment yield (SY) and perform any further analysis necessary. SDR evaluation for SEDD modeling requires the development of three (3) main parameters: the slope, hydraulic length, and the β coefficient. The "flow distance" module of the ArcGIS software was employed to develop the horizontal flow paths and required the use of the "flow accumulation" and "flow direction" modules. The slope map and other component datasets of the USLE/SEDD SDR evaluation for SEDD modeling requires the development of three (3) main parameters: the slope, hydraulic length, and the β coefficient. The "flow distance" module of the ArcGIS software was employed to develop the horizontal flow paths and required the use of the "flow accumulation" and "flow direction" modules. The slope map and other component datasets of the USLE/SEDD model were developed in the ArcGIS platform.
This study relied upon the R programming language. Specifically, the "raster," "e1071," and "LaplacesDemon" libraries were employed for statistical analysis, conversion of data from raster to point datasets, and other analysis, while the "ggplot2" libraries of R were employed for data visualization [32][33][34][35].
The USLE model in this study was developed in the ArcGIS model builder framework updating the model originally developed by Jhan [36], Yang [37], Li [8], and Liu [38].
The USLE model was developed by Wischmeier and Smith [39,40] and is designed to predict the average annual SL by sheet and rill erosion. It has been continuously improved upon and localized by many researchers over the past decades. The USLE model can be denoted as follows: where A m is the computed soil loss per unit area (t/ha/year), R m is the rainfall and runoff factor (MJ-mm ha −1 h −1 year −1 ), K m is the soil erodibility factor (t-h MJ −1 mm −1 ), L is the slope-length factor (dimensionless), S is the slope-steepness factor (dimensionless), C is the cover and management factor, and P is the support practice factor.

Sediment Distributed Delivery Model (SEDD)
The SEDD utilizes the calculation of the flow length for SDR developed by Ferro and Minacapilli [24], which includes considerations for slope and the roughness/runoff coefficient. SY, the sediment yield by the annual scale, for a basin is discretized into N i morphological units and is measured in tonnes (t): where A mi is the computed soil loss per unit area (t/ha/year) for the morphological unit i, SU i is the area of the morphological unit i (ha), and N i is the number of morphological units over the watershed area (dimensionless). Ferro and Minacapilli [24] concluded that the SDR w can be physically defined using watershed morphological data without the consideration of the soil loss model. SDR i in Equation (2) and the relationship between the SDR w and SDR i was defined as follows (Equations (3) and (4)): where t p,i is the travel time from ith morphological unit to the nearest stream reach (m), β is the roughness and runoff coefficient for the watershed (m −1 ), s p,i is the slope of the hydraulic flow path (m/m), l p,i is the length of the hydraulic flow path (m), λ i,j and S i,j are the length and slope of each morphological unit i localized along the hydraulic path j, and N p is the number of morphological units localized along the hydraulic path j [6,24,[26][27][28][29].
In this study, the parent watershed (the Shihmen Reservoir watershed) was discretized into two discrete sets of children sub-watersheds, specifically 25 and 93 sub-watersheds. The 25 sub-watersheds were delineated using the SWAT model with the input of the DEM and river course of the Shihmen Reservoir watershed. The output was 25 hydrologically connected sub-watersheds, as shown in Figure 1b. The 93 sub-watersheds were delineated using the ArcGIS Watershed module utilizing an input of the flow direction and a stream map derived from flow accumulation. The output was 93 hydrologically connected sub-watersheds, as shown in Figure 1c.
The SEDD model was implemented as 10 m grid data, with over seven million data points and aggregated into sub-watersheds. This study derived the grid-based SEDD model and aggregated the dataset into sub-watersheds to examine the impact this can have on the output SY data. A grid-cell-based analysis is more computationally expensive, and larger watersheds can be quite time consuming or hardware-expensive but may provide alternative results. On the other hand, the hydrological unit analysis is less time consuming to produce, has a lower computational cost, and can garner results that are aggregated by specific areas. However, there can be discrepancies due to oversimplification, some vast.

Sediment Yield Determination
Flow direction (FD) algorithms estimate the flow of a specified material from a source cell (S1) continuously to the next neighboring cells (N1) until a stopping point or limitation point (river or stream in Figure 3) has been reached, such as a stream, river, or dam. The SDR model considers the hydraulic flow of sediments through each morphological/hydrological unit to the trunk river or its tributaries (in this scenario, the Tahan river, and its tributaries) using the single flow direction algorithm, deterministic D8. The distance traveled by particles has been termed the "flow path length" and has varying definitions and methods of calculation. The flow distance module in ArcGIS was used to define the horizontal flow distance to the river using the D8 flow direction. unit analysis is less time consuming to produce, has a lower computational cost, and can garner results that are aggregated by specific areas. However, there can be discrepancies due to oversimplification, some vast.

Sediment Yield Determination
Flow direction (FD) algorithms estimate the flow of a specified material from a source cell (S1) continuously to the next neighboring cells (N1) until a stopping point or limitation point (river or stream in Figure 3) has been reached, such as a stream, river, or dam. The SDR model considers the hydraulic flow of sediments through each morphological/hydrological unit to the trunk river or its tributaries (in this scenario, the Tahan river, and its tributaries) using the single flow direction algorithm, deterministic D8. The distance traveled by particles has been termed the "flow path length" and has varying definitions and methods of calculation. The flow distance module in ArcGIS was used to define the horizontal flow distance to the river using the D8 flow direction.
The relative time taken to reach the sink from the source has been termed the "travel time." Moreover, Jain and Kothyari [27] state that the travel time from a cell i is equivalent to the total travel time through each of the Np cells along its flow path until terminating at a channel. The flow path length and the slope in conjunction with a watershed-wide β coefficient are used to derive, firstly, the travel time and, secondly, the SDR within the SEDD model. The SY was first derived using a watershed-wide grid of the USLE and SDRi, and then compared against the discretization of the 25 sub-watersheds model and the 93 sub-watersheds model using the mean, median, and mode for each sub-watershed.

Results
The results of SEDD modeling in the Shihmen Reservoir watershed and its sub-watersheds are shown in the following sections.

The β Coefficient
Using GIS, we could calculate the flow lengths, slope, and travel time needed in Equations (3) and (4) for the Shihmen Reservoir watershed. At the same time, soil loss could be computed by USLE, and the SDRw is known to be 0.49 previously determined from the annual siltation of the reservoir [8]. Therefore, the only unknown is Equation (4) was the β coefficient. By assuming an initial value of β, we can solve for the β coefficient iteratively to satisfy the Equation (4). We found the β coefficient to be 0.00515. This is the first result of applying the SEDD model to the Shihmen Reservoir watershed using the iterative method, and a number like this has never been obtained before. The determined β coefficient is watershed specific. The relative time taken to reach the sink from the source has been termed the "travel time." Moreover, Jain and Kothyari [27] state that the travel time from a cell i is equivalent to the total travel time through each of the N p cells along its flow path until terminating at a channel. The flow path length and the slope in conjunction with a watershed-wide β coefficient are used to derive, firstly, the travel time and, secondly, the SDR within the SEDD model.
The SY was first derived using a watershed-wide grid of the USLE and SDR i , and then compared against the discretization of the 25 sub-watersheds model and the 93 sub-watersheds model using the mean, median, and mode for each sub-watershed.

Results
The results of SEDD modeling in the Shihmen Reservoir watershed and its sub-watersheds are shown in the following sections.

The β Coefficient
Using GIS, we could calculate the flow lengths, slope, and travel time needed in Equations (3) and (4) for the Shihmen Reservoir watershed. At the same time, soil loss could be computed by USLE, and the SDR w is known to be 0.49 previously determined from the annual siltation of the reservoir [8]. Therefore, the only unknown is Equation (4) was the β coefficient. By assuming an initial value of β, we can solve for the β coefficient iteratively to satisfy the Equation (4). We found the β coefficient to be 0.00515. This is the first result of applying the SEDD model to the Shihmen Reservoir watershed using the iterative method, and a number like this has never been obtained before. The determined β coefficient is watershed specific.

Grid-Based Travel Time, SDR, SL, and SY
The following Figures 4-7 depict the travel time, sediment delivery ratio (SDR), soil loss (SL), and sediment yield (SY) determined using the SEDD model.
The probability density functions (PDF) of the grid-cell-based SDR, SL, and SY of the Shihmen Reservoir watershed have asymmetric, right-skewed distributions (which will be shown later in Section 4.3 with two representative sub-watersheds). This study determined in the Shihmen Reservoir watershed for SDR: the mean value of the distribution was 0.474 (not the same as SDR w ), the median was 0.459, and the mode was 0.317. For the SL of the entire watershed, the mean was estimated as 87.07 t/ha/year, and the median was 71.74 t/ha/year. Lastly, the mean SY for the entire watershed was 42.08 t/ha/year, and the median was 29.80 t/ha/year.

Grid-Based Travel Time, SDR, SL, and SY
The following Figures 4-7 depict the travel time, sediment delivery ratio (SDR), soil loss (SL), and sediment yield (SY) determined using the SEDD model.         The probability density functions (PDF) of the grid-cell-based SDR, SL, and SY of the Shihmen Reservoir watershed have asymmetric, right-skewed distributions (which will be shown later in section 4.3 with two representative sub-watersheds). This study determined in the Shihmen Reservoir watershed for SDR: the mean value of the distribution was 0.474 (not the same as SDRw), the median was 0.459, and the mode was 0.317. For the SL of the entire watershed, the mean was estimated as 87.07 t/ha/year, and the median was 71.74 t/ha/year. Lastly, the mean SY for the entire watershed was 42.08 t/ha/year, and the median was 29.80 t/ha/year.

Comparison by Sub-Watershed Discretization on SDR, SL, and SY
This study considered the distribution of SDR, SL, and SY in the Shihmen Reservoir watershed discretized into 93 sub-watersheds. From Figure 8, it is evident that both the median and the mean of SDR by 93 SW had very similar distributions with a central tendency centered around 0.500.

Comparison by Sub-Watershed Discretization on SDR, SL, and SY
This study considered the distribution of SDR, SL, and SY in the Shihmen Reservoir watershed discretized into 93 sub-watersheds. From Figure 8, it is evident that both the median and the mean of SDR by 93 SW had very similar distributions with a central tendency centered around 0.500.
The range of mean SL by sub-watershed was between 15.60 t/ha/year (SW #8) and 151.21 t/ha/year (SW #45) and the range of the median of SL was 1.15 t/ha/year (SW #5) and 136.36 t/ha/year (SW #69). The evaluation of the median of the SL for 93 SW showed that the median tended to be between 60 and 90 t/ha/year, and over 20 sub-watersheds had a median between 60-70 t/ha/year. In contrast, the mean SL by sub-watersheds were more widely spread and had a number of peaks between 60 and 120 t/ha/year with over 20 sub-watersheds valued at 80 to 100 t/ha/year. The range of the mean SY by sub-watershed in the 93 sub-watershed discretization was between 10.28 t/ha/year (SW#1) and 92.97 t/ha/year (SW#73) while the range of the median was 0.38 t/ha/year (SW#8) and 88.11 t/ha/year (SW #73). For the 25 SW, the mean SY was between 10.52 t/ha/year (SW #3) and 73.92 t/ha/year (SW #8) and the range of the median was between 5.50 t/ha/year (SW #3) and 49.43 t/ha/year (SW #16).
To investigate the mean, median, and mode values by sub-watershed discretization (25 SW and 93 SW) and the correlation of each aspect of sediment yield analysis (SDR, SL, and SY), this study created the boxplots displayed in Figure 9a-e. This study additionally developed the PDF plot of each sub-watershed to investigate similarities and differences between the sub-watersheds-a sample of these plots from the 93 SW are displayed in Figure 10. The variance in the shapes of the probability density functions of SDR, SL, and SY was evident. The range of mean SL by sub-watershed was between 15.60 t/ha/year (SW #8) and 151.21 t/ha/year (SW #45) and the range of the median of SL was 1.15 t/ha/year (SW #5) and 136.36 t/ha/year (SW #69). The evaluation of the median of the SL for 93 SW showed that the median tended to be between 60 and 90 t/ha/year, and over 20 sub-watersheds had a median between 60-70 t/ha/year. In contrast, the mean SL by sub-watersheds were more widely spread and had a number of peaks between 60 and 120 t/ha/year with over 20 sub-watersheds valued at 80 to 100 t/ha/year. The range of the mean SY by sub-watershed in the 93 sub-watershed discretization was between 10.28 t/ha/year (SW#1) and 92.97 t/ha/year (SW#73) while the range of the median was 0.38 t/ha/year (SW#8) and 88.11 t/ha/year (SW #73). For the 25 SW, the mean SY was between 10.52 t/ha/year (SW #3) and 73.92 t/ha/year (SW #8) and the range of the median was between 5.50 t/ha/year (SW #3) and 49.43 t/ha/year (SW #16).
To investigate the mean, median, and mode values by sub-watershed discretization (25 SW and 93 SW) and the correlation of each aspect of sediment yield analysis (SDR, SL, and SY), this study created the boxplots displayed in Figure 9a-e. This study additionally developed the PDF plot of each sub-watershed to investigate similarities and differences between the sub-watersheds-a sample of these plots from the 93 SW are displayed in Figure 10. The variance in the shapes of the probability density functions of SDR, SL, and SY was evident.

Discussion
Past studies have shown the distribution of soil erosion of the Shihmen Reservoir watershed, and this paper estimates the geospatially distributed sediment delivery ratio (SDR i ) under a beta determined by the recursive method of the SEdiment Delivery Distributed (SEDD) model and a known watershed-wide SDR w .

The Importance of SDR
The SDR plays an essential role in soil erosion research as an additional parameter when considering sediment yield, slope conservation/remediation, and sediment control projects for engineers and decision-makers. An SDR of almost 100% denotes an area where the dislodged sediments have a near-perfect chance to reach a nearby river or stream. In comparison, an SDR of 0% within a highly erodible soil may mean that the soil has a high likelihood of depositing before reaching the river or stream. The SDR consideration creates a twofold issue for engineers: the classification of sediment loads entering the stream and the classification of the deposition quantities on the slopes that may lead to future sediment related disasters.
The traditional method of conveying the SDR across a region has historically been studies of small watersheds or empirical equations. The SDR values obtained from these methods were extrapolated to the entire watershed and utilized to calculate the SY from the gross erosion model's output SL. With the advancement of GIS, it is now possible to use morphological units or grid cells to compute the SDR for each location in the watershed. However, unlike SL or SY, it is essential to note that the individual SDR (SDR i ) cannot be averaged directly. The average of the SDR i is not the SDR w because of the definition of SDR. To calculate the SDR w , individual SY has to be computed from the SDR i first. SDR w is then the ratio of the sum of SY divided by the total soil loss. In this study, we calculated the mean of SDR i purely for the sake of statistical analysis.
The determined SDR can be used to determine SY, which is then used as a basis for decision-making, countermeasure designs, and other uses. Therefore, the accuracy of these evaluations is critical to the management of a watershed, such as the Shihmen Reservoir watershed, where significant soil erosion risk is apparent.

The β Coefficient in the Model
The SDR calculation within the SEDD model utilizes one unique parameter requiring expert opinion or field experiments, β. The β coefficient in the model is calculated using three methodologies, the recursive method, the trial-and-error method, and the field experiments method, using Cs 137 [6,[42][43][44]. The β coefficient is a watershed-wide constant, but in some cases, this parameter is evaluated for each spatial discretization unit (for sub-watersheds). Additionally, some researchers have subjectively set a value using known published values of β [45]. The "recursive approach" or "recursive fitting approach" as used in this study defines β using known values of SDR w . When β is determined, SDR i can be determined from Equation (3). Sometimes, SDR w is obtained from empirical equations. For example, Vanoni's method [46] was used to estimate β as a SEDD parameter [6]. Burguet et al. [44] evaluated β over two olive orchard catchments using different annual C factors and R factors, determining the median β from multiple events and improving the assessment of the β parameter for their two study catchments. Other studies have utilized an additional factor, the vegetation or land-cover parameter, α [27,47], but Lopez-Vicente and Navas [29] concluded that the SDR in the SEDD model has limited sensitivity to the land use type. The slope and length of the flow paths have more influence on the model. Additionally, Lai [48] explored the use of the SEDD in the Shihmen Reservoir watershed but determined the β coefficient by setting it to values between 0 and 200 using the Jain and Kothyari's variant of the original model which introduces the coefficient a i (also denoted as k i ) to consider different land uses based on expert opinions from a table introduced by Haan [49]. However, this study utilizes the original model developed by Ferro and Minacapilli [24] because the use of the other model introduces high spatial variance in the a i coefficient [27] and increases uncertainty to the SDR model. This is because the expert opinions have to be utilized to compare Haan's table to land use types not considered originally or are geographically dissimilar. This study determines the β = 0.00515 using the known SDR w value, while Lai [48] determined β = 8.5. For these reasons, this study has been distinguished as a more accurate implementation of the SEDD model.
The SDR w (0.49) used in this study was derived by Li [8], which was calculated from reservoir sedimentation, and the soil loss included sheet and rill erosion, gully erosion, and landslides sediments. This recursive method is discussed thoroughly by Porto and Walling [42] and is employed in regions lacking data to obtain reliable values for β. The recursive method of the SEDD model in our study was developed within the model builder of the ArcGIS software. It extended the works of the model originally developed by Jhan et al. [50] and Chen et al. [51].

SEDD Model Using the Grid-Cell-Based Analysis
The SEDD model determines the sediment yield at the outlet from each hydrological unit or grid cell unit. First, the USLE equation is applied to the study area, and the gross erosion derived is then reduced by the SDR to the equivalent sediment yield at the outlet. The SDR equation within the SEDD model uses physical and hydrological parameters (slope and flow distance, respectively) to determine the sediment delivery to channels and eventually to the reservoir outlet.
The correlation between the SY, SL, and the ratio SDR is evident in the previously shown Figures 5-7. A visual comparison between Figures 4 and 5 can easily distinguish that the SDR values are inversely related to the travel time. As the travel time or distance to a river channel increases, the probability of eroded particles entering the stream decreases and consequently SDR decreases.
In Figure 5, the SDR of the watershed ranged from 0 to 1, and the areas where sediment delivery has a higher likelihood to reach channels or the outlet is focused mainly around the rivers and streams distributed throughout the watershed. The soil loss evaluation by USLE has been previously discussed in Chen et al. [51] and Liu et al. [16]. This study builds upon the model of Liu et al. [16] by increasing the number of rainfall stations in the analysis of R m to 41 stations in and around the watershed [38]. Significantly, the SL is 0-40 t/ha/year at the lower elevation in the northern areas of the watershed that surround the reservoir. In contrast, in the southern reaches of the watershed of higher elevations, the SL increases to 80 t/ha/year or more and often exceeding 100 t/ha/year. The delivery to the channel is limited by the SDR, as shown in Figure 5. Therefore, the results of SY in these areas follow the pattern of the SDR with increases around the river channels. The influence of the SDR on SY is evident from visual inspection of the maps.

Mean, Median, and Mode of 25/93 Sub-Watersheds
This study determined the SDR and SY of the SEDD model. It analyzed the probability density function of the SDR, SL, and SY for a discretization of the watershed into 25 sub-watersheds and 93 sub-watersheds. The distributions of SDR, SL, and SY for the watershed and sub-watersheds were found to be all non-normal distributions.
For each sub-watershed, we calculated its mean, median, and mode of SDR, SL, and SY from its distributions. The results were then aggregated into Figure 9 using boxplots. Additionally shown in Figure 9a-c are the mean and median of the entire watershed (horizontal dash lines). Comparing the 25 sub-watersheds with the 93 sub-watersheds, we can see that the range (between the ends of whiskers) of 25 SW is always smaller than that of the 93 SW, no matter what the distribution type is (SDR, SL, or SY). This means that with finer discretization, there is more variation, and it is more likely to obtain extreme values. Interestingly, the interquartile ranges (IQR) of the means of SDR, SL, and SY of 25 SW and 93 SW are about the same. It can also be noted that although the mean values of SDR and SL of 93 SW have broader ranges than those of the 25 SW, the resulting mean values of SY for both the 25 and 93 SW have about the same range. Figure 9d shows the distribution of the watershed level SDR, SL, and SY. Compared with Figure 9b,c, the ranges of the watershed SL and SY are larger than those of both sub-watershed discretizations. Figure 9e shows the skewness of the SDR, SL, and SY datasets of the 25 SW and 93 SW in comparison with the skewness of the grid-cell-based analysis of the entire watershed. The difference in the skewness of the SDR between all three discretizations (Shihmen as a whole, 25 SW, and 93 SW) is small. Their range (SDR) is much smaller than that of the SL and SY. Additionally, the skewness of the SL and SY are similarly distributed. The lower quartile range is more dominant in the 93 SW, while the 25 SW shows a generally even distribution between the upper and low quartiles. The skewness of the entire watershed SDR was 0.153, for SL 8.359 and for sediment yield 9.134. The upper whiskers of the 93 SW for SL and SY are significantly longer than the lower whiskers. Sixty-two of the 93 sub-watersheds (66.7%) for SDR were positively skewed (right-skewed), meaning the majority of SDR were below the 0.500 while for the 25 sub-watersheds, 17 were positively skewed (68.8%). For the SL, 24 of 25 sub-watersheds (96%) were positively skewed (right-skewed), and 81 of 93 sub-watersheds (87%) were, too. All of the 25 sub-watersheds (100%) and all 93 sub-watersheds were positively skewed for the PDF of SY. The similar positive relationship between the skewness of SL and SY shows that the SL is more influential on the SY than the SDR is. Both example cases in Figure 10  Applying Equation (4) to each of the 25 or 93 sub-watersheds, we also calculated the sub-watershed SDR w, as shown in Figure 11 (not averaging of SDR i ). It can be seen that the discretization of 93 sub-watershed shows a broader range of SDR w values than the 25 sub-watersheds. These SDR w values are essential if a sub-watershed level study is to assess detailed SY features for each sub-watershed. The SDR w can also serve as one of the criteria to select the most critical sub-watersheds to monitor.
Sustainability 2020, 12, x FOR PEER REVIEW 16 of 20 Figure 9e shows the skewness of the SDR, SL, and SY datasets of the 25 SW and 93 SW in comparison with the skewness of the grid-cell-based analysis of the entire watershed. The difference in the skewness of the SDR between all three discretizations (Shihmen as a whole, 25 SW, and 93 SW) is small. Their range (SDR) is much smaller than that of the SL and SY. Additionally, the skewness of the SL and SY are similarly distributed. The lower quartile range is more dominant in the 93 SW, while the 25 SW shows a generally even distribution between the upper and low quartiles. The skewness of the entire watershed SDR was 0.153, for SL 8.359 and for sediment yield 9.134. The upper whiskers of the 93 SW for SL and SY are significantly longer than the lower whiskers. Sixty-two of the 93 sub-watersheds (66.7%) for SDR were positively skewed (right-skewed), meaning the majority of SDR were below the 0.500 while for the 25 sub-watersheds, 17 were positively skewed (68.8%). For the SL, 24 of 25 sub-watersheds (96%) were positively skewed (right-skewed), and 81 of 93 sub-watersheds (87%) were, too. All of the 25 sub-watersheds (100%) and all 93 sub-watersheds were positively skewed for the PDF of SY. The similar positive relationship between the skewness of SL and SY shows that the SL is more influential on the SY than the SDR is. Both example cases in Figure  Applying Equation (4) to each of the 25 or 93 sub-watersheds, we also calculated the sub-watershed SDRw, as shown in Figure 11 (not averaging of SDRi). It can be seen that the discretization of 93 sub-watershed shows a broader range of SDRw values than the 25 sub-watersheds. These SDRw values are essential if a sub-watershed level study is to assess detailed SY features for each sub-watershed. The SDRw can also serve as one of the criteria to select the most critical sub-watersheds to monitor.

Sub-Watershed Prioritization
Soil Conservation and countermeasures against sedimentation of river waterways are a costly endeavor and involve significant budgetary concerns for nations such as Taiwan. Climate change and its effects, such as the increase of global land surface temperatures and changing precipitation patterns, pose increased concerns for water security for large population centers such as Taipei City and Taoyuan. These budgetary concerns reinforce the need for sustainable watershed management and evidence-based decision making for targeted soil conversation and soil erosion mitigation programs. Therefore, sub-watershed prioritization can improve the effectiveness of these interventions while

Sub-Watershed Prioritization
Soil Conservation and countermeasures against sedimentation of river waterways are a costly endeavor and involve significant budgetary concerns for nations such as Taiwan. Climate change and its effects, such as the increase of global land surface temperatures and changing precipitation patterns, pose increased concerns for water security for large population centers such as Taipei City and Taoyuan. These budgetary concerns reinforce the need for sustainable watershed management and evidence-based decision making for targeted soil conversation and soil erosion mitigation programs. Therefore, sub-watershed prioritization can improve the effectiveness of these interventions while controlling budgets. In this study, the three measurements of central tendency (mean, median, and mode) of SDR, SL, and SY were explored. The SEDD model introduces a physically-based methodology for determining the SDR that is useful in understanding what probability of soil will enter the stream or outlet while modeling changing rainfall or land cover levels.
The prioritization by mean SY for 25 and 93 sub-watersheds is shown in Figure 12, each with five levels of prioritization from 0 to 20 t/ha/year (low) to 80 to 100 t/ha/year (high) derived from the ranges previously discussed. It is evident that there are striking similarities between both of these prioritization schemes. Generally speaking, the SY at the outlet is low, but in the eastern and southern extremes of the watershed, there are specific sub-watersheds with medium to high sediment yield values. controlling budgets. In this study, the three measurements of central tendency (mean, median, and mode) of SDR, SL, and SY were explored. The SEDD model introduces a physically-based methodology for determining the SDR that is useful in understanding what probability of soil will enter the stream or outlet while modeling changing rainfall or land cover levels.
The prioritization by mean SY for 25 and 93 sub-watersheds is shown in Figure 12, each with five levels of prioritization from 0 to 20 t/ha/year (low) to 80 to 100 t/ha/year (high) derived from the ranges previously discussed. It is evident that there are striking similarities between both of these prioritization schemes. Generally speaking, the SY at the outlet is low, but in the eastern and southern extremes of the watershed, there are specific sub-watersheds with medium to high sediment yield values.

Conclusions
This study applied the SEdiment Delivery Distributed (SEDD) model to the Shihmen Reservoir watershed in Taiwan. The main merit of the study lies in the first attempt to derive the sediment yield by soil erosion for the entire Shihmen Reservoir watershed and its sub-watersheds. By using the recursive method, this study determined the SEDD β coefficient to be 0.00515 and predicted the spatial distributions (maps) of travel time, SDRi, soil erosion (soil loss), and sediment yield of the Shihmen Reservoir watershed. The resulting average of SY for the Shihmen Reservoir watershed was 42.08 t/ha/year and by sub-watershed 10.52-73.92 t/ha/year for 25 sub-watersheds and 10.28-92.97 t/ha/year for 93 sub-watersheds.
This study also recommended the use of mean SY by sub-watershed in sub-watershed prioritization for future soil conservation decision-making. The model presented takes advantage of all currently available data for the Shihmen reservoir watershed. However, it is essential to note that the sediment yield modeling can be improved by increased on-site validation and the use of aerial photogrammetry to deliver more updated data to better understand the field situations. Considering the implications of climate change, there is also a great need for further research on the sediment yield by soil erosion and other land degradation sources and the impacts of changing precipitation regimes. Using geospatial models of sediment yield as guidance can help the local government to better im-

Conclusions
This study applied the SEdiment Delivery Distributed (SEDD) model to the Shihmen Reservoir watershed in Taiwan. The main merit of the study lies in the first attempt to derive the sediment yield by soil erosion for the entire Shihmen Reservoir watershed and its sub-watersheds. By using the recursive method, this study determined the SEDD β coefficient to be 0.00515 and predicted the spatial distributions (maps) of travel time, SDR i , soil erosion (soil loss), and sediment yield of the Shihmen Reservoir watershed. The resulting average of SY for the Shihmen Reservoir watershed was 42.08 t/ha/year and by sub-watershed 10.52-73.92 t/ha/year for 25 sub-watersheds and 10.28-92.97 t/ha/year for 93 sub-watersheds.
This study also recommended the use of mean SY by sub-watershed in sub-watershed prioritization for future soil conservation decision-making. The model presented takes advantage of all currently available data for the Shihmen reservoir watershed. However, it is essential to note that the sediment yield modeling can be improved by increased on-site validation and the use of aerial photogrammetry to deliver more updated data to better understand the field situations. Considering the implications of climate change, there is also a great need for further research on the sediment yield by soil erosion and other land degradation sources and the impacts of changing precipitation regimes. Using geospatial models of sediment yield as guidance can help the local government to better implement engineering and ecological solutions for soil conservation to achieve sustainable land management (SLM), thereby reducing sediment yield risks and creating a well-balanced solution for all stakeholders.