Soil Erosion Modeling and Comparison Using Slope Units and Grid Cells in Shihmen Reservoir Watershed in Northern Taiwan

Soil erosion is a global problem that will become worse as a result of climate change. While many parts of the world are speculating about the effect of increased rainfall intensity and frequency on soil erosion, Taiwan’s mountainous areas are already facing the power of rainfall erosivity more than six times the global average. To improve the modeling ability of extreme rainfall conditions on highly rugged terrains, we use two analysis units to simulate soil erosion at the Shihmen reservoir watershed in northern Taiwan. The first one is the grid cell method, which divides the study area into 10 m by 10 m grid cells. The second one is the slope unit method, which divides the study area using natural breaks in landform. We compared the modeling results with field measurements of erosion pins. To our surprise, the grid cell method is much more accurate in predicting soil erosion than the slope unit method, although the slope unit method resembles the real terrains much better than the grid cell method. The average erosion pin measurement is 6.5 mm in the Shihmen reservoir watershed, which is equivalent to 90.6 t ha−1 yr−1 of soil erosion.


Introduction
Soil erosion is the detachment of soil particles occurring in a landscape. Compared with wind erosion, water erosion is a much more critical process that produces enormous sediments in Taiwan's watersheds. There are many factors that contribute to the high soil erosion rate of Taiwan. The first and foremost is Taiwan's hilly topography. Situated at the collision boundary of the continental Eurasian plate and the Philippine Sea plate, the mountain range of Taiwan is being lifted at a rate of 0.2-18.5 mm yr −1 [1], which creates some of the most rugged terrain in East Asia. The second problem is Taiwan's thick soil deposits due to heavy weathering as a consequence of the island's tropical to subtropical weather. Finally, the last problem that exacerbates the soil erosion conditions in Taiwan is the heavy rain brought by monsoons and typhoons. Due to the amount and intensity of the rainfall, the rainfall erosivity value of Taiwan is many times higher than the global average obtained by Borrelli et al. [2], making water erosion a very serious concern in Taiwan. In fact, not just in Taiwan, the concern of soil erosion has been evident in many published studies in different parts of the world, including Europe [3,4], the US and Mexico [5,6], mainland China [7,8], India [9,10], Africa [11,12], Australia [13,14], etc. A steady increase in the number of studies on soil erosion in the last few years has also prompted a few regional and global studies that provide a bigger picture of the emerging crisis [2,[15][16][17]. If we factor in the threat of climate change, we will likely face more intense precipitation occurring at a higher frequency in the near future. In order to be prepared for the possible global increase in soil erosion, studying existing high soil erosion areas such as Taiwan's watersheds may provide insights into a preventive measure against a massive increase of global soil loss.
In this study, we selected a representative watershed in Taiwan as the study area. The watershed was monitored by 550 erosion pins installed at 55 locations from 8 September 2008 to 10 October 2011. Two variations (grid cells and slope units) of the soil erosion model based on the universal soil loss equation (USLE) were used to compute the soil loss due to sheet and rill erosion. The results will be compared with the measurement of erosion pins in different sub-watersheds. For the purpose of this research, we assume that the data layers used in the analysis are representative of the study area even though the data were not gathered at the same time. The goals are to examine the difference between a grid cell analysis and an analysis based on natural land boundaries (slope units), and to determine which is more accurate in predicting erosion pin measurements.

Materials and Methods
This study took a modeling and validation approach. That is, we first developed a model to estimate the amount of soil erosion. Then, we validated the model using field measurements of erosion pins. The soil erosion model that we built is based on the digital elevation model (DEM) acquired by the Central Geological Survey (CGS) of Taiwan in 2013 with airborne LiDAR (light detection and ranging). The horizontal resolution of the DEM was 10 m. The soil erosion was computed using the USLE with two different analysis units: grid cells and slope units.

Universal Soil Loss Equation
The universal soil loss equation (USLE) was used to compute the sheet and rill erosion by water. Despite its shortcomings, the USLE and its successor RUSLE (Revised Universal Soil Loss Equation) is still the most common soil erosion model for data scarce catchments. It is reported that at least 90% of soil erosion studies conducted worldwide use RUSLE-based models [18]. According to Wischmeier and Smith, the USLE "is an erosion model designed to compute longtime average soil losses from sheet and rill erosion" [19]. The equation is (in metric units): where Am = computed soil loss per unit area (t ha −1 yr −1 ) Rm = rainfall and runoff factor (MJ-mm ha −1 hr −1 yr −1 ) Km = soil erodibility factor (t-hr MJ −1 mm −1 ) L = slope-length factor S = slope-steepness factor C = cover and management factor P = support practice factor

Study Area
The Shihmen reservoir watershed in northern Taiwan was selected as the study area. The reservoir itself is the third largest in Taiwan, and was constructed with an earth-and-rock fill embankment. As the reservoir supplies municipal, industrial, and irrigation water to northern Taiwan, it is often the focus of public attention. The location of the Shihmen reservoir watershed is shown in Figure 1.
Marked on Figure 1 are the locations of field observations as photographed on 16 March 2018. Figure 2a is an example of soil erosion turned into a shallow landslide near the Hsia-Yun Creek after typhoon Nesat in 2017. The site also exhibits a typical recurrence of the alternation of erosion and landslides, which are common in the area. Figure 2b shows the Ronghua Dam, which is situated 27 km upstream of the Shihmen reservoir, and was constructed in 1983 with a capacity of 12.4 million m 3 . The Ronghua Dam serves as a hydropower facility, but its main function was as a check dam to protect the Shihmen reservoir (Shihmen Dam). As can be seen from Figure 2b, the Ronghua Dam (at the junction of two mountains near the center of the image) is almost filled with sediments. Figure 2c is one of the 55 sites monitored with erosion pins (10 pins per site) near Bao-Li-Ku Creek. The measurements were taken with a caliper. A detailed description of the installation of the erosion pins can be found in Lin et al. [20]. Finally, Figure 2d shows a series of check dams at Su-Le Creek. The check dams are filled with sediments and have completely lost their protection function. As can be seen from these field photos, the sediment problem in the Shihmen reservoir watershed is very alarming, and a good model for estimating the amount of soil erosion in the watershed is crucial. March 2018. Figure 2a is an example of soil erosion turned into a shallow landslide near the Hsia-Yun Creek after typhoon Nesat in 2017. The site also exhibits a typical recurrence of the alternation of erosion and landslides, which are common in the area. Figure 2b shows the Ronghua Dam, which is situated 27 km upstream of the Shihmen reservoir, and was constructed in 1983 with a capacity of 12.4 million m 3 . The Ronghua Dam serves as a hydropower facility, but its main function was as a check dam to protect the Shihmen reservoir (Shihmen Dam). As can be seen from Figure 2b, the Ronghua Dam (at the junction of two mountains near the center of the image) is almost filled with sediments. Figure 2c is one of the 55 sites monitored with erosion pins (10 pins per site) near Bao-Li-Ku Creek. The measurements were taken with a caliper. A detailed description of the installation of the erosion pins can be found in Lin et al. [20]. Finally, Figure 2d shows a series of check dams at Su-Le Creek. The check dams are filled with sediments and have completely lost their protection function. As can be seen from these field photos, the sediment problem in the Shihmen reservoir watershed is very alarming, and a good model for estimating the amount of soil erosion in the watershed is crucial.

Grid Cells and Slope Units
In this study, the rainfall factor R was obtained from Lu et al. [21], and the soil erodibility factor K was obtained from Wann and Hwang [22]. The average R factor was 12,798 MJ-mm ha −1 hr −1 yr −1 , or 6.4 times the global average of 2000 MJ-mm ha −1 hr −1 yr −1 [23]. The K factor ranged from 0.004 to 0.042 t-hr MJ −1 mm −1 [24]. The C factor was converted from a land use survey in 2004 [25] and a conversion table between the land uses and the C factors [26], and the average C value was 0.016. Finally, the P factor was assumed to be 1, as is usually done in the literature when there are no protection measures.
The topographic factors L and S were derived from the CGS DEM. They were computed using Equations (2) and (3) from [19]: where λ = slope length in meters m = 0.5 if percent slope is 5% or more = 0.4 on slopes of 3.5% to 4.5% [19]

Grid Cells and Slope Units
In this study, the rainfall factor R was obtained from Lu et al. [21], and the soil erodibility factor K was obtained from Wann and Hwang [22]. The average R factor was 12,798 MJ-mm ha −1 hr −1 yr −1 , or 6.4 times the global average of 2000 MJ-mm ha −1 hr −1 yr −1 [23]. The K factor ranged from 0.004 to 0.042 t-hr MJ −1 mm −1 [24]. The C factor was converted from a land use survey in 2004 [25] and a conversion table between the land uses and the C factors [26], and the average C value was 0.016. Finally, the P factor was assumed to be 1, as is usually done in the literature when there are no protection measures.
The topographic factors L and S were derived from the CGS DEM. They were computed using Equations (2) and (3) from [19]: where λ = slope length in meters m = 0.5 if percent slope is 5% or more = 0.4 on slopes of 3.5% to 4.5% [19] or 3% to 5% [26] = 0.3 on slopes of 1% to 3% = 0.2 on uniform gradients of less than 1% where θ = the angle of slope The original USLE handbook makes no reference to measuring the slope length in horizontal projection [19,27], but the RUSLE handbook specifically specifies such measurements [28]. Some studies (such as Oliveira et al. [29]) use slope lengths parallel to the surface in their computations.
GIS (Geographic Information System) and soil erosion modeling usually go hand in hand, especially for large areas. As the study area is too vast to be covered in a single plot or unit, the Shihmen reservoir watershed is divided into many smaller units to facilitate computation and analysis. In GIS, the easiest way to accomplish the division of a study area is by using square grid cells, as shown in Figure 3a. Many spatial analysis tools in GIS support the use of gridded datasets. It is a superior method to a station-based approach, but could sometimes provide misleading information [30]. Still, many studies of soil erosion use grid cells in their analysis [9]. On the other hand, the grid cells do not conform to any natural land boundaries and the results may be biased. This is especially problematic for computing the slope-length factor L, which is the most difficult factor to calculate among the six factors of the USLE. To approximate the shape of the hillslopes and compute the slope length λ in Equation (2), different methods have been proposed, such as Hickey [31] and Sextante [32], the unit contributing area [33,34], and a combination of the unit contributing area and LS-TOOL [35]. In addition, Foster and Wischmeier [36] propose dividing a slope into segments. The original USLE handbook makes no reference to measuring the slope length in horizontal projection [27,19], but the RUSLE handbook specifically specifies such measurements [28]. Some studies (such as Oliveira et al. [29]) use slope lengths parallel to the surface in their computations.
GIS (Geographic Information System) and soil erosion modeling usually go hand in hand, especially for large areas. As the study area is too vast to be covered in a single plot or unit, the Shihmen reservoir watershed is divided into many smaller units to facilitate computation and analysis. In GIS, the easiest way to accomplish the division of a study area is by using square grid cells, as shown in Figure 3(a). Many spatial analysis tools in GIS support the use of gridded datasets. It is a superior method to a station-based approach, but could sometimes provide misleading information [30]. Still, many studies of soil erosion use grid cells in their analysis [9]. On the other hand, the grid cells do not conform to any natural land boundaries and the results may be biased. This is especially problematic for computing the slope-length factor L, which is the most difficult factor to calculate among the six factors of the USLE. To approximate the shape of the hillslopes and compute the slope length λ in Equation (2), different methods have been proposed, such as Hickey [31] and Sextante [32], the unit contributing area [33,34], and a combination of the unit contributing area and LS-TOOL [35]. In addition, Foster and Wischmeier [36] propose dividing a slope into segments. In this study, the first approach used was the grid cell method. We first assumed that hillslopes could be divided into uniform 10 m by 10 m squares, and the amount of soil erosion could be evaluated on a cell-by-cell basis. For the slope length factor, two kinds of procedures were employed to calculate the slope lengths-in directions parallel to the slope surfaces [29] and in a horizontal projection. Then for contrast, we relaxed the assumption regarding the computing soil erosion of each cell individually, and used the SAGA (System for Automated Geoscientific Analyses) software to calculate the LS factors.
To contrast the grid cell approach and be more representative of the field situation, slope units (sometimes called terrain units) were also used in this study. In contrast to grid cells, slope units are based on natural terrains of the study area, as shown in Figure 3b. The natural water divides between distinct geographic areas were used as the boundaries of the units. In this method, each slope unit is a close representation of an individual slope, which makes the application of the USLE more suitable. Note that the slope units depend on the breaks of topography. The soil properties and groundwater conditions may be the same or very similar between adjacent units.
Both approaches have their respective advantages and disadvantages. For the grid cell approach, its advantage lies in its simplicity and applicability. However, it is too restrictive on the slope length for the assumption that we made. On the other hand, the slope unit approach places no restriction on the length of the slope, but the approach is difficult to implement and automate. In this study, both the grid cell approach and the slope unit approach were used, and their results In this study, the first approach used was the grid cell method. We first assumed that hillslopes could be divided into uniform 10 m by 10 m squares, and the amount of soil erosion could be evaluated on a cell-by-cell basis. For the slope length factor, two kinds of procedures were employed to calculate the slope lengths-in directions parallel to the slope surfaces [29] and in a horizontal projection. Then for contrast, we relaxed the assumption regarding the computing soil erosion of each cell individually, and used the SAGA (System for Automated Geoscientific Analyses) software to calculate the LS factors.
To contrast the grid cell approach and be more representative of the field situation, slope units (sometimes called terrain units) were also used in this study. In contrast to grid cells, slope units are based on natural terrains of the study area, as shown in Figure 3b. The natural water divides between distinct geographic areas were used as the boundaries of the units. In this method, each slope unit is a close representation of an individual slope, which makes the application of the USLE more suitable. Note that the slope units depend on the breaks of topography. The soil properties and groundwater conditions may be the same or very similar between adjacent units.
Both approaches have their respective advantages and disadvantages. For the grid cell approach, its advantage lies in its simplicity and applicability. However, it is too restrictive on the slope length for the assumption that we made. On the other hand, the slope unit approach places no restriction on the length of the slope, but the approach is difficult to implement and automate. In this study, both the grid cell approach and the slope unit approach were used, and their results were compared. Since the grid cell approach is commonly used in GIS analysis, it does not require further explanation. The following will explain how slope units were delineated in this study.
Although rarely used in soil erosion, the slope unit approach has often been applied to the study of landslides. The slope units were created by inverting the DEM to intersect with the original DEM. The intersections form the required slope units [37]. Unfortunately, the procedure is not fool-proof, and it requires manual inspection and editing. The editing includes the correction of slope angles and slope directions (Figure 4), removal of peculiar areas ( Figure 5), filling in of missing areas (Figure 6), and disregarding shadowed areas and river surfaces (Figure 7). The use of slope units is also called the object-oriented approach [18,38]. were compared. Since the grid cell approach is commonly used in GIS analysis, it does not require further explanation. The following will explain how slope units were delineated in this study. Although rarely used in soil erosion, the slope unit approach has often been applied to the study of landslides. The slope units were created by inverting the DEM to intersect with the original DEM. The intersections form the required slope units [37]. Unfortunately, the procedure is not fool-proof, and it requires manual inspection and editing. The editing includes the correction of slope angles and slope directions (Figure 4), removal of peculiar areas ( Figure 5), filling in of missing areas (Figure 6), and disregarding shadowed areas and river surfaces (Figure 7). The use of slope units is also called the object-oriented approach [18,38].     were compared. Since the grid cell approach is commonly used in GIS analysis, it does not require further explanation. The following will explain how slope units were delineated in this study. Although rarely used in soil erosion, the slope unit approach has often been applied to the study of landslides. The slope units were created by inverting the DEM to intersect with the original DEM. The intersections form the required slope units [37]. Unfortunately, the procedure is not fool-proof, and it requires manual inspection and editing. The editing includes the correction of slope angles and slope directions (Figure 4), removal of peculiar areas ( Figure 5), filling in of missing areas (Figure 6), and disregarding shadowed areas and river surfaces (Figure 7). The use of slope units is also called the object-oriented approach [18,38].    were compared. Since the grid cell approach is commonly used in GIS analysis, it does not require further explanation. The following will explain how slope units were delineated in this study. Although rarely used in soil erosion, the slope unit approach has often been applied to the study of landslides. The slope units were created by inverting the DEM to intersect with the original DEM. The intersections form the required slope units [37]. Unfortunately, the procedure is not fool-proof, and it requires manual inspection and editing. The editing includes the correction of slope angles and slope directions (Figure 4), removal of peculiar areas ( Figure 5), filling in of missing areas (Figure 6), and disregarding shadowed areas and river surfaces (Figure 7). The use of slope units is also called the object-oriented approach [18,38].

Results
This is the first time that a slope unit-based approach is applied to soil erosion modeling in the Shihmen reservoir watershed. The outcomes are interesting compared to the grid cell approach (10 m by 10 m). Here we will discuss our findings in four parts.

Delineation of the Study Area
The result of the delineation of the slope units is shown in Figure 8. All in all there are 19,275 slope units generated in the Shihmen reservoir watershed with an average size of 3.87 ha. It is obvious from Figure 8 that the slope units form a more realistic representation of the terrain of the study area. Therefore, the subsequent modeling of soil erosion based on slope units is likely to be more convincing.

Results
This is the first time that a slope unit-based approach is applied to soil erosion modeling in the Shihmen reservoir watershed. The outcomes are interesting compared to the grid cell approach (10 m by 10 m). Here we will discuss our findings in four parts.

Delineation of the Study Area
The result of the delineation of the slope units is shown in Figure 8. All in all there are 19,275 slope units generated in the Shihmen reservoir watershed with an average size of 3.87 ha. It is obvious from Figure 8 that the slope units form a more realistic representation of the terrain of the study area. Therefore, the subsequent modeling of soil erosion based on slope units is likely to be more convincing.

Results
This is the first time that a slope unit-based approach is applied to soil erosion modeling in the Shihmen reservoir watershed. The outcomes are interesting compared to the grid cell approach (10 m by 10 m). Here we will discuss our findings in four parts.

Delineation of the Study Area
The result of the delineation of the slope units is shown in Figure 8. All in all there are 19,275 slope units generated in the Shihmen reservoir watershed with an average size of 3.87 ha. It is obvious from Figure 8 that the slope units form a more realistic representation of the terrain of the study area. Therefore, the subsequent modeling of soil erosion based on slope units is likely to be more convincing.

Comparison of Topographic Factors L and S
The statistics of the L and S factors are shown in Table 1. Two important observations can be made. First, the ranges and the averages of the S factors were almost the same for both approaches. This is reasonable in that the steepness of a slope should not be altered simply due to the shape or the size of the unit of analysis (at least not significantly). This intuition was supported by our results. Second, as anticipated, the slope length factors (L) were considerably different for both approaches. As shown in Table 1, the average L derived from the slope units was 3.7 times the average L of the grid cells. As USLE is multiplicative, the L factor alone could increase the amount of soil erosion to 3.7 times of the original amount.

Comparison of Soil Erosion
The computed amounts of soil erosion for (1) the grid cell approach; (2) the grid cell approach with a horizontal projection; (3) the grid cell approach using SAGA; and (4) the slope unit approach are shown in Figure 9 using the internationally accepted scale (a, c, e, g) and our own classes (b, d, f, h). As can be seen from Figure 9, the result of the grid cells was pixel-like, whereas the result of the slope units was irregular in shape. All approaches show that the study area has high soil erosion, but the amounts of soil erosion are different. The average amount of soil erosion for the grid cell approach was 97.0 t ha −1 yr −1 ; for the grid cell approach with a horizontal projection it was 82.1 t ha −1 yr −1 ; for the grid cell approach using SAGA it was 280.1 t ha −1 yr −1 ; and for the slope unit approach it was 416.5 t ha −1 yr −1 . In other words, the slope unit method has an amount of soil erosion 4.3 times higher than the grid cell approach. Since most of the remaining five USLE factors are similar in value, this result is in line with the observed ratio of 3.7 between the average L factors of the two approaches.

Comparison of Topographic Factors L and S
The statistics of the L and S factors are shown in Table 1. Two important observations can be made. First, the ranges and the averages of the S factors were almost the same for both approaches. This is reasonable in that the steepness of a slope should not be altered simply due to the shape or the size of the unit of analysis (at least not significantly). This intuition was supported by our results. Second, as anticipated, the slope length factors (L) were considerably different for both approaches. As shown in Table 1, the average L derived from the slope units was 3.7 times the average L of the grid cells. As USLE is multiplicative, the L factor alone could increase the amount of soil erosion to 3.7 times of the original amount.

Comparison of Soil Erosion
The computed amounts of soil erosion for (1) the grid cell approach; (2) the grid cell approach with a horizontal projection; (3) the grid cell approach using SAGA; and (4) the slope unit approach are shown in Figure 9 using the internationally accepted scale (a, c, e, g) and our own classes (b, d, f, h). As can be seen from Figure 9, the result of the grid cells was pixel-like, whereas the result of the slope units was irregular in shape. All approaches show that the study area has high soil erosion, but the amounts of soil erosion are different. The average amount of soil erosion for the grid cell approach was 97.0 t ha −1 yr −1 ; for the grid cell approach with a horizontal projection it was 82.1 t ha −1 yr −1 ; for the grid cell approach using SAGA it was 280.1 t ha −1 yr −1 ; and for the slope unit approach it was 416.5 t ha −1 yr −1 . In other words, the slope unit method has an amount of soil erosion 4.3 times higher than the grid cell approach. Since most of the remaining five USLE factors are similar in value, this result is in line with the observed ratio of 3.7 between the average L factors of the two approaches.   Figure 10 shows the distribution of the 55 sites monitored by erosion pins (10 pins per site), and Table 2 shows the averages of the erosion pin measurements from 8 September 2008 to 10 October 2011 [39]. We divided the entire Shihmen reservoir watershed into five smaller sub-watersheds as shown in Figure 10. The sub-watersheds are Ku-Chu, Yu-Feng, San-Kuang, Pai-Shih, and Tai-Kang sub-watersheds. Then, we computed the average depth of the soil erosion in each sub-watershed using both the grid cell approach and the slope unit approach in Equation (4):

Validation with Erosion Pin Measurements
where At = total amount of soil erosion (t) d = erosion depth (mm) area = area of watershed (ha) γ = unit weight of soil (t m −3 )  Figure 10 shows the distribution of the 55 sites monitored by erosion pins (10 pins per site), and Table 2 shows the averages of the erosion pin measurements from 8 September 2008 to 10 October 2011 [39]. We divided the entire Shihmen reservoir watershed into five smaller sub-watersheds as shown in Figure 10. The sub-watersheds are Ku-Chu, Yu-Feng, San-Kuang, Pai-Shih, and Tai-Kang sub-watersheds. Then, we computed the average depth of the soil erosion in each sub-watershed using both the grid cell approach and the slope unit approach in Equation (4):

Validation with Erosion Pin Measurements
where At = total amount of soil erosion (t) d = erosion depth (mm) area = area of watershed (ha) The results are shown in Table 2. Here we assume γ = 1.4 t m −3 to compute the equivalent erosion depths in each sub-watershed and the entire Shihmen reservoir watershed. As can be seen from Table 2, the erosion depth of the entire Shihmen reservoir watershed computed by the grid cell approach is very close to the average measurement at the erosion pins (6.9 mm vs. 6.5 mm), whereas the erosion depth computed by the slope unit method was many times over the average erosion pin measurement. Individually speaking, the grid cell method was good in some sub-watersheds (the Ku-Chu sub-watershed and Tai-Kang sub-watershed) but not so good in the others (the Yu-Feng sub-watershed, San-Kuang sub-watershed, and Pai-Shih sub-watershed). If we convert the measured average erosion depth to the amount of soil erosion, 6.5 mm is equivalent to 90.6 t ha −1 yr −1 . The results are shown in Table 2. Here we assume γ = 1.4 t m −3 to compute the equivalent erosion depths in each sub-watershed and the entire Shihmen reservoir watershed. As can be seen from Table 2, the erosion depth of the entire Shihmen reservoir watershed computed by the grid cell approach is very close to the average measurement at the erosion pins (6.9 mm vs. 6.5 mm), whereas the erosion depth computed by the slope unit method was many times over the average erosion pin measurement. Individually speaking, the grid cell method was good in some sub-watersheds (the Ku-Chu sub-watershed and Tai-Kang sub-watershed) but not so good in the others (the Yu-Feng sub-watershed, San-Kuang sub-watershed, and Pai-Shih sub-watershed). If we convert the measured average erosion depth to the amount of soil erosion, 6.5 mm is equivalent to 90.6 t ha −1 yr −1 .

Discussion and Conclusions
In this study, we modeled the soil erosion of the Shihmen reservoir watershed in northern Taiwan with two different analysis units, both of which were based on a 10 m resolution DEM acquired through airborne laser scanning. The first was a grid cell approach, and the second was a slope unit approach. As grid cells do not resemble real terrain, it was anticipated that the slope unit approach would substantially improve the estimation of soil erosion. However, contrary to our expectations, the result was the opposite. The grid cell method predicted an average erosion depth of

Discussion and Conclusions
In this study, we modeled the soil erosion of the Shihmen reservoir watershed in northern Taiwan with two different analysis units, both of which were based on a 10 m resolution DEM acquired through airborne laser scanning. The first was a grid cell approach, and the second was a slope unit approach. As grid cells do not resemble real terrain, it was anticipated that the slope unit approach would substantially improve the estimation of soil erosion. However, contrary to our expectations, the result was the opposite. The grid cell method predicted an average erosion depth of 6.9 mm, very close to the actual average erosion depth of 6.5 mm (based on erosion pin measurements at 55 locations). On the other hand, the slope unit method yielded an average L factor 3.7 times the average L factor of the grid cell method, and the slope unit method produced an average amount of soil erosion 4.3 times higher than that of the grid cell method. Note that the grid cell approach using SAGA also did not yield results closer to the measured erosion depths. A cut-off slope angle of 50% or 55% was suggested by McCool et al. [40], Liu et al. [41], and Panagos et al. [34], but was not carried out in this study.
In addition, when we examined the five sub-watersheds of the Shihmen reservoir watershed, we observed good matches between the predictions of the grid cell method and the measurement of the erosion pins for some sub-watersheds but not others. The discrepancy suggests that there are variations between sub-watersheds that cannot be properly accounted for by the grid cell approach. Fortunately, the variations seemed to cancel one another out to produce an average soil erosion (measured) that was very close to the prediction of the grid cell method. On reflection, the grid cell method simply makes more accurate predictions than the slope unit method (despite its perceived sound modeling basis).
The USLE has been criticized for its many limitations and proneness to misuse. Compared with physical-based soil erosion models, it certainly seems that USLE grossly oversimplifies the soil erosion problem and should not be used (with the availability of better models). However, the frequent lack of the necessary data to run a physical-based soil erosion model almost always requires something that is as simple as the USLE. The fact that reliance has to be placed on less data-demanding models such as the USLE suggests that there is still a gap between the ideal and the current situation. In this study, we found that despite the many inherent inaccuracies of the USLE and the simplifications that we made, the USLE still predicted soil erosion at a watershed scale reasonably well. Therefore, the USLE will probably remain one of the very few soil erosion models that can be applied to large regions with scarce data.