2.1. RUSLE

The USLE was designed to estimate sheet and rill erosion from hillslope area, but not address soil deposition and channel or gully erosion within a watershed (^{Renard et al., 1991}). The RULSE is an index method containing factors that represent how climate, soil, topography, and land use affect rill and interrill soil erosion caused by raindrop impact and surface runoff. The RUSLE equation is as follows:

where Loss is the soil loss(ton/ha·year), R is the rainfall erosivity factor, K is a soil erodibility factor, L is the slope length factor, S is the slope steepness factor, C is a cover management factor, and P is a supporting practices factor.

The R-factor represents the long-term average erosivity of the climate, calculated by the total rainfall energy (E) and the maximum 30 min rainfall intensity (I30) for rainfall events. The R-factor can be calculated as follows :

where E is the total storm kinetic energy (MJ ha^-1), I30 is the maximum 30-min intensity (mm h^-1), N is number of years, ei is the rainfall energy per unit depth of rainfall (MJ ha^-1 mm^-1 h^-1), I is the rainfall intensity (mm h^-1), and ΔV_i is the duration of the increment over which I is constant in hours (h). Guak calculated the R-factors in 2003 using equations (3) ~ (5) for eight rainfall gauge stations monitored every 1 minute within the Lake Imha watersheds (Fig. 1 and Table 2) (^{Guak, 2007}). A grid-based R-factor map was generated using the Spline interpolation method of ArcView 3.0 (^{ESRI, 2002}), as shown in Fig. 2.

Fig. 1. Study area.

Fig. 2. The 30×30 m grid based RUSLE factor maps for the Lake Imha watershed.

The K-factor was determined using the Erickson triangular nomograph method, considering the percentage of sand, silt, and clay in the soil (^{Erickson, 1997}). A 1:2500 scaled soil map including information on soil texture was obtained from the Korean National Academy of Agricultural Science. K-values were allocated to the soil map and a grid-based K-factor map was generated, as shown in Fig. 2.

A grid-based LS-factor map was generated using the SATEEC GIS system (^{Park et al., 2010}) with the Digital Elevation Model (DEM) obtained from the Environmental Geographic System (EGIS) (^{ME, 2014}). The SATEEC GIS system uses Moore and Burch’s method, using the following equation (^{Moore and Burch, 1986}):

where A is the flow accumulation cell size, and θ is the slope angle in degrees.

^{Park (2003)} allocated C values to land use as classified by the Korean Ministry of Environment by referencing documents, for which the C values are shown in Table 3. The land use classification map was obtained from EGIS, and the grid-based C-factor map is shown in Fig. 2.

The P-factor considers the support practice of protecting cultivated areas from soil erosion. Williams and Smith proposed P values considering a combination of land slope and the supporting cropland practices including contouring, contour strip cropping, and terracing (^{Wischmeier and Smith, 1978}). In Korea, rice paddy fields and other agricultural areas can be considered as terrace systems and contour tillage, respectively (^{Guak, 2007}; ^{Park, 2003}). The P values, from Williams and Smith (^{Wischmeier and Smith, 1978}), are shown in Table 4. The p-factor map was generated using DEM and land use map and is shown in Fig. 2.

2.2. Research approaches

The SATEEC GIS system was used to generate the annual soil loss map for 2003. The significant water quality problem, especially high turbidity concentration by soil loss, in the Imha multi-purpose dam occurred in 2003. An overview of the application of the SATEEC GIS system is shown in Fig. 3 (^{Lim et al., 2005}). Total suspended sediment load at the outlet of subwatershed simulated by HSPF was considered as the suspended sediment loads because HSPF considers the deposition or scour in streams. The research was divided by part 1 which was performed by ^{Jeon et al. (2016)} and part 2. This paper concerns part 2. A diagram of the modeling approach is shown in Fig. 5.

Fig. 3. Diagram of SATEEC GIS system application.

Fig. 5. Flow diagram of research approaches.

In Part 1, HSPF was calibrated at the uncalibrated subwatersheds (Fig. 1) by matching with sediment loads at the outlet of the Lake Imha watershed, which were calculated by soil loss and SDR. The process of validation was not performed due to the data limitation. The soil loss was clipped for the subwatersheds and the suspended sediment loads from subwatersheds 4, 9, 13, 15, 36, and 46 was calculated by HSPF (Fig. 1).

In Part 2, SDRs of the six calibrated subwatersheds for 2003 were calculated using the ratio of soil erosion by RUSLE and suspended sediment loads with the calibrated HSPF. Various topographic parameters were used in development of the SDR equations with reference to Table 1. The parameters included area (A, km²), watershed relief (R_f, m), curve number (CN), channel slope (SLP_ch), drainage slope (SLP_dr), the ratio of watershed relief to watershed length (R_f/L_b, m/km), and the ratio of watershed relief to main channel length (R_f/L_w, m/km). ^{Jeon et al. (2016)} optimized the curve number for the combination of land use and hydrologic soil group for the Lake Imha watershed. The calibrated grid based on the CN map by Jeon et al. was used in this study (Fig. 4). The watershed relief was calculated by taking the difference between the highest and lowest elevations using DEM. Channel and drainage slope were calculated using the BASINS Delineation Tool. The SDR equations developed in this study were evaluated by comparing the magnitude of the change in sediment loads within the subwatershed. Statistical analyses were performed with the IBM SPSS Statistics program (ver.21). Sediment loads at the outlet of the Lake Imha watershed were calculated by multiplying soil erosion calculated with RUSLE and the SDR calculated using the SDR equation.

Fig. 4. SCS-CN map of the Lake Imha watershed.

3.1. Development of SDR equation

The application of RUSLE revealed that about 2,070,580 ton/year of soil eroded from Lake Imha watershed during 2003. The LS-factor and soil erosion maps obtained using the SATEEC GIS system are shown in Figs. 6 and 7. The SATEEC GIS system is a useful tool for the application of RUSLE and analysis of the spatial distribution of soil erosion. The SDRs calculated by the ratio of soil loss from RUSLE and sediment load from HSPF for the calibrated subwatersheds are presented in Table 5. The SDRs ranged from 0.025 to 0.172, and decreased with increasing drainage area. The subwatershed parameters and correlation analysis with the SDRs are shown in Tables 6 and 7, respectively. SDRs were strongly correlated with the R_f/L_w, R_f/L_w, CN, and SLP_ch, showing correlation coefficients (R) of 0.95, 0.93, 0.80, and 0.76, respectively. The SDR equations for the Lake Imha watershed were developed, as shown in equations (7) ~ (11), referring to Table 1 and the results of the correlation analysis.

Fig. 6. LS and soil loss map within the Lake Imha watershed using the SATEEC GIS system.

Fig. 7. Soil loss maps for the subwatersheds covered by monitoring gauge stations.

Some SDRs calculated by equation (8) were about zero when the channel slope was gradual, as shown in Fig. 8, indicating that equation (8) was not appropriate for use as the SDR equation for the Lake Imha watershed. The differences of sediment inflow and outflow at outlet 33 calculated with equations (7) and (11) were significant, showing differences of -48% and 28%, respectively (Table 8). Considering the relatively small area of subwatershed 33 (14.8 km²), those differences seem unrealistic. At the point of outlet 33 and 47, equations (9) and (10) generated similar SDRs and were reasonable compared with (7) and (11) so the two equations could be recommended for SDR equation in the Lake Imha watershed. Fig. 9.

Fig. 8. Comparison of the SDRs for six calibrated subwatersheds calculated by HSPF and RUSLE through the SDR equations.

Fig. 9. Sediment mass balance between subwatersheds 33 and 47.

Many researchers have proposed area-based power functional SDR equations that decrease with increasing drainage area (Table 1). However, this type of equation can sometimes have significant error at watersheds that have major tributaries close to the outlet, as is the case for the Lake Imha watershed. The area-based power functional SDR equation in this study would be unrealistically decreased after entering major tributaries that have relatively large drainage areas.

3.2. HSPF calibration for uncalibrated subwatersheds

HSPF was calibrated for uncalibrated subwatersheds using the sediment load calculated by multiplying the soil loss and SDR using equation (10), named the ‘SDR equation’ in Table 9, instead of the observed sediment load. Table 6 shows the sediment loads of Lake Imha by calibrated HSPF (HSPF_cal), by using the calibrated HSPF parameters of the spatial nearest neighbor (HSPF_emp), and by using the default HSPF parameters (HSPF_def). The relative errors were calculated by the difference between the SDR equation and the three kinds of HSPF simulation results. Although the uncalibrated subwatershed area was 37% of the total Lake Imha watershed, HSPF_def and HSPF_emp caused significant errors, showing 112% and 48% of relative errors, respectively. A more reasonable method for the determination of HSPF parameters related to sediment calibration was to use the sediment load determined by RUSLE with the SDR equation. ^{U. S. EPA (2006)} guided the sediment calibration of HSPF, and proposed calibration of HSPF coupled with RUSLE and SDR. However, the SDR equation is very sensitive to site specifications, so the user should carefully select an accurate SDR equation for use.

3.3. Characteristics of suspended sediment inflow to Lake Imha

The simulated yearly suspended sediment inflow to Lake Imha during 1996 ~ 2010, calculated by HSPF, and the statistical analysis are shown in Fig. 10 and Table 10, respectively. The average yearly suspended sediment inflow to Lake Imha was 26,506 ton/year. A significantly higher inflow of suspended sediment to Lake Imha occurred during 2002, 71,900 ton/year. The maximum yearly suspended sediment load was 11 times higher than the minimum loads.

Fig. 10. Yearly suspended sediment inflow to Lake Imha during 2000 ~ 2010.

Statistical analysis of the monthly suspended sediment loads entering Lake Imha calculated by HSPF for 1996 ~ 2000 is shown in Table 11 and Fig. 11. The highest monthly suspended sediment was loaded during July, showing 12,627 ton/month, which was up to 40% of the total yearly load. The second and third highest loads were 7,369 ton/month (28%) during August and 2,996 ton/month (11%) during June. Most of the suspended sediment entering Lake Imha was loaded during June ~ August, amounting to 79% of the total yearly sediment load. This is a common characteristic of runoff and nonpoint source pollution loaded in the Asian summer monsoon climate, for which most rainfall events occur during June ~ August (^{Kettering et al., 2012}; ^{Kim et al., 2014}).

Fig. 11. Monthly suspended sediment inflow to Lake Imha during 1996 ~ 2010.