**Effects of Different Spatial Configuration Units for the Spatial Optimization of Watershed Best Management Practice Scenarios**

#### **Liang-Jun Zhu 1,2, Cheng-Zhi Qin 1,2,3,\*, A-Xing Zhu 1,2,3,4,5, Junzhi Liu 3,4 and Hui Wu <sup>6</sup>**


Received: 22 December 2018; Accepted: 31 January 2019; Published: 2 February 2019

**Abstract:** Different spatial configurations (or scenarios) of multiple best management practices (BMPs) at the watershed scale may have significantly different environmental effectiveness, economic efficiency, and practicality for integrated watershed management. Several types of spatial configuration units, which have resulted from the spatial discretization of a watershed at different levels and used to allocate BMPs spatially to form an individual BMP scenario, have been proposed for BMP scenarios optimization, such as the hydrologic response unit (HRU) etc. However, a comparison among the main types of spatial configuration units for BMP scenarios optimization based on the same one watershed model for an area is still lacking. This paper investigated and compared the effects of four main types of spatial configuration units for BMP scenarios optimization, i.e., HRUs, spatially explicit HRUs, hydrologically connected fields, and slope position units (i.e., landform positions at hillslope scale). The BMP scenarios optimization was conducted based on a fully distributed watershed modeling framework named the Spatially Explicit Integrated Modeling System (SEIMS) and an intelligent optimization algorithm (i.e., NSGA-II, short for Non-dominated Sorting Genetic Algorithm II). Different kinds of expert knowledge were considered during the BMP scenarios optimization, including without any knowledge used, using knowledge on suitable landuse types/slope positions of individual BMPs, knowledge of upstream–downstream relationships, and knowledge on the spatial relationships between BMPs and spatial positions along the hillslope. The results showed that the more expert knowledge considered, the better the comprehensive cost-effectiveness and practicality of the optimized BMP scenarios, and the better the optimizing efficiency. Thus, the spatial configuration units that support the representation of expert knowledge on the spatial relationships between BMPs and spatial positions (i.e., hydrologically connected fields and slope position units) are considered to be the most effective spatial configuration units for BMP scenarios optimization, especially when slope position units are adopted together with knowledge on the spatial relationships between BMPs and slope positions along a hillslope.

**Keywords:** spatial configuration units; best management practices (BMPs); spatial optimization; hydrologic response units (HRUs); hydrologically connected fields; slope positions; watershed process simulation

#### **1. Introduction**

Different best management practices (or beneficial management practices, BMPs for short) scenarios (i.e., spatial configurations of multiple BMPs) at the watershed scale may have significantly different environmental effectiveness, economic efficiency, and practicality [1–6]. They are valuable for decision-making of integrated watershed management to assess the environmental effectiveness and economic efficiency of watershed BMP scenarios and then propose optimal ones. Currently, a popular approach to achieving this target is based on watershed modeling [7–9] coupled with intelligent optimization algorithms [2,4,5,10–13], so-called BMP scenarios optimization. To conduct the BMP scenarios optimization, each individual BMP scenario is created by automatically selecting and allocating BMPs on spatial configuration units (also called BMP configuration units hereafter), which have resulted from the spatial discretization of a watershed at one among different levels (such as subbasins, and hydrologic response units; Figure 1). When a specific type of BMP configuration unit is chosen for the BMP scenarios optimization of a watershed, normally each individual BMP configuration unit in the watershed is allowed to be configured with only one type of BMP (as the situation in this study). Then the effects of the scenario on watershed behavior are simulated by watershed models [4,13]. The simulation result is the basis of the automatic spatial optimization of BMP scenarios. Therefore, the determination of the BMP configuration units becomes the one key issue for the BMP scenarios optimization.

**Figure 1.** Schematic diagram of spatial discretization of a watershed at different levels (such as sub-basins, hydrologic response units (HRUs), farms, hydrologically connected fields, slope position units, and grid cells) and the corresponding best management practice (BMP) scenario examples after allocating multiple types of BMPs based on different kinds of expert knowledge.

Currently, spatial configuration units used in BMP scenarios optimization mainly include five types with different levels, i.e., subbasins [14–17], hydrologic response units (HRUs) [12,18,19], farms [10,11], hydrologically connected fields [5], and slope position units (i.e., landform positions at hillslope scale, such as ridge, backslope, etc.) [2].

Subbasins are normally regarded as relatively closed and independent spatial configuration units which can be further delineated into different levels of spatial configuration units such as landform positions, HRUs, and gridded cells [20]. Thus, subbasins are coarse-grained units for configuring spatially explicit BMPs, and it is suitable to serve as BMP configuration units for situations in which the potential locations of BMPs are predefined within each subbasin [17]. Meanwhile, using subbasins as BMP configuration units can decrease the search space of spatial optimization more than that of the adoption of other more detailed spatial configuration units, and thus can achieve a better optimizing efficiency [15]. However, using subbasins as BMP configuration units is not available for configuring multiple spatially explicit BMPs within one subbasin while BMPs often have different effects on different spatial positions of a hillslope [2,21,22].

HRUs are delineated as hydrologic homogeneous areas according to landuse, soil, and topography (commonly presented through slope percentage) within one subbasin [23]. HRUs are extensively used for the sake of convenience, especially when the Soil and Water Assessment Tool (SWAT) [23] is applied to watershed modeling. However, an HRU is spatially discrete—it may occupy the hillslope from ridge to valley bottom—and is even spatially ambiguous when the three thresholds introduced in SWAT for delineating HRUs (i.e., the minimum percentage of the landuse area over the subbasin area, the soil area over the landuse area, and the slope class area over the soil area) are set as non-zero values [24]. Although setting these thresholds to be non-zero can reduce the count of HRUs and hence improve the computational efficiency of simulations, the inappropriate representation of the study area may introduce some levels of ambiguity in simulations. Therefore, the selected optimal BMP scenarios may not be easy to implement spatially explicitly. Similar to subbasins, HRUs are also inappropriate to serve as BMP configuration units for spatially explicit BMPs. Furthermore, the impact of BMPs of an upslope HRU on a downslope HRU cannot be represented [24], although this impact is important for incorporating the expert knowledge of spatial relationships between BMPs and spatial positions [2,5].

To improve the practicality of BMP scenarios, some studies adopted farms [10,11], or the so-called spatially explicit HRUs that are directly defined by farm boundaries [25] or landuse/landcover field maps [26], as BMP configuration units. These BMP configuration units are spatially one-to-one matched with farm or landuse/landcover fields and thus can be collectively referred to as spatially explicit HRUs [26]. Compared with spatially discrete HRUs, spatially explicit HRUs can make corresponding BMP scenarios more practical to stakeholders (such as farmers, landowners, or land managers). Meanwhile, due to the ignorance of topographic variance inside each farm or landuse/landcover field, such BMP configuration units normally have a smaller count than spatially discrete HRUs for a study area, which means higher optimizing efficiency [1,26]. However, there are still no explicitly defined upstream–downstream relationships between spatially explicit HRUs, which means that the spatial relationships between BMPs and spatial positions cannot be represented effectively [2,5].

Wu et al. [5] proposed hydrologically connected fields with upstream–downstream relationships to be BMP configuration units, which can be delineated by considering spatial topology based on flow directions and a landuse map. With hydrologically connected fields, a set of expert rules of BMP interactions based on the upstream–downstream relationships can be developed in the form of "if-then" rules, i.e., if a field has been configured with one BMP, its adjacent upstream fields should not be configured with BMP, otherwise, BMP will be randomly selected and configured on the adjacent upstream fields according to their landuse type [5]. To the best of our knowledge, the study of Wu et al. [5] is the first work on incorporating the spatial relationships between BMPs and spatial positions into BMP scenarios optimization, so that its cost-effectiveness and optimizing efficiency could be improved. However, hydrologically connected fields may be delineated across multiple landform positions or subbasins, since the hydrological relationship, considered by Wu et al. [5], is built at the watershed scale rather than the hillslope or subbasin scale. This means that hydrologically connected fields have weak spatial relationships to homogeneous functional units at the hillslope scale from the perspectives of physical geography such as geomorphic, soil, and hydrologic conditions [2]. Therefore, hydrologically connected fields also face the shortcoming that the spatially explicit BMPs cannot be represented effectively. A possible way for spatially explicit HRUs and hydrologically connected fields to overcome the above-mentioned shortcoming is to delineate them so as to be small enough patchworks of gridded cells within homogeneous functional units, or even individual gridded cells [27–29]. However, this approach will render the optimization based on such spatial configuration units computationally intensive or even unsolvable [27]. This makes such an approach only suitable for the spatial optimization of one single BMP within a little watershed [28,29], thus having too narrow applicability for normal watershed management which considers multiple BMPs.

More recently, Qin et al. [2] proposed slope position units [30,31] as BMP configuration units, so as to overcome the above-mentioned shortcomings of hydrologically connected fields and other BMP configuration units. Slope position units (e.g., ridge, backslope, and valley), as spatially contiguous and topographically connected units along the hillslope, are inherently related to physical watershed processes [20,30,32], and thus affect the effectiveness of BMPs [2]. Besides, the count of slope position units is comparatively limited, so as to ensure optimizing efficiency. Therefore, the spatial relationships between BMPs and slope positions along the hillslope could be explicitly and effectively considered during the spatial optimization of BMP scenarios. According to the preliminary study [2], this spatial-relationship-considered way of using slope position units as BMP configuration units is effective, efficient, and practical for BMP scenarios optimization, compared to a standard random optimization way of selecting and allocating BMPs randomly on configuration units.

Although many studies have assessed the cost-effectiveness of each individual type of above-mentioned BMP configuration units for the spatial optimization of BMP scenarios, as far as we know, a comparison among main types of BMP configuration units based on the same one watershed model for an area is still lacking. To discuss the effects of different BMP configuration units for watershed BMP scenarios optimization with regard to cost-effectiveness, optimizing efficiency, and practicality, this paper compares four types of BMP configuration units (i.e., HRUs, spatially explicit HRUs, hydrologically connected fields, and slope position units) in the spatial optimization of configuring multiple spatially explicit BMPs for mitigating soil erosion based on a spatially distributed watershed model. Subbasins are not included in this study, because subbasins are too coarse-grained to represent multiple spatially explicit BMPs.

#### **2. Materials and Methods**

#### *2.1. Methodology*

To compare the effects of these four types of BMP configuration units for watershed BMP scenarios optimization, a widely used spatial optimization framework of BMP scenarios based on watershed modeling coupled with intelligent optimization algorithms [2,10,12] was adopted in this study. As shown in Figure 2, the spatial optimization framework of watershed BMP scenarios mainly consists of four components: (1) BMP configuration units for allocating BMPs within the watershed; (2) a BMP knowledge base together with BMP configuration units as inputs to generate and evaluate BMP scenarios; (3) models for evaluating each watershed BMP scenario, including a distributed watershed model that can simulate spatial interactions between spatially explicitly distributed BMPs and effectively assess the environmental effectiveness of each BMP scenario [2,5,33], and a BMP scenario cost model for estimating the economic efficiency of each BMP scenario; (4) a multi-objective optimization component based on an intelligent optimization algorithm such as NSGA-II (Non-dominated Sorting Genetic Algorithm II) [34], which includes initializing BMP scenarios based on BMP configuration units and a BMP knowledge base, and generating new BMP scenarios or proposing optimal ones based on the evaluation results of all current BMP scenarios. These components are elaborated in four subsections (Sections 2.3–2.6) followed by the subsection of the study area and dataset (Section 2.2).

**Figure 2.** Framework comparing the effects of different BMP configuration units for the spatial optimization of watershed BMP scenarios in this study (extended from [2]).

Besides knowledge on the environmental effectiveness and cost-benefit of individual BMPs, BMP configuration knowledge (such as the suitable landuse types/slope positions for individual BMPs, and the spatial relationships between BMPs and spatial positions) is necessary for BMP scenarios optimization (Figure 2). According to the characteristics of BMP configuration units, different kinds of BMP configuration knowledge can be saved in the BMP knowledge base and applied (hereafter referred to as BMP configuration strategy) (Figure 2). For example, the simple rule of suitable landuse types/slope positions is applicable for all types of BMP configuration units, while the expert rules based on upstream–downstream relationships between BMPs can only be applied to the BMP configuration units which have upstream–downstream relationships [5], such as hydrologically connected fields and slope position units. Furthermore, the expert rules based on the spatial relationships between BMPs and slope positions along the hillslope are only applicable for slope position units [2]. Section 2.4 contains more detailed information about the BMP knowledge base and BMP configuration strategies adopted in this study.

To assess the effects of four types of BMP configuration units for watershed BMP scenarios optimization, all feasible combinations of BMP configuration units and BMP configuration strategies were investigated by the same watershed BMP scenarios optimization framework (Figure 2). The results of different BMP configuration units applied with available BMP configuration strategies are discussed from the perspectives of cost-effectiveness, optimizing efficiency, and practicality (Figure 2). The detailed experimental design is described in Section 2.7.

#### *2.2. Study Area and Dataset*

The Youwuzhen watershed (~5.39 km2), located in the typical severely eroded red-soil hilly region in southeastern China, was selected as the study area (Figure 3). Low hills with steep slopes (up to 52.9◦ and with an average slope of 16.8◦) and broad alluvial valleys are the primary geomorphology forms [2]. Forest (59.8%), paddy field (20.6%), and orchard (12.8%) are the primary landuse types in the study area (Figure 4). Additionally, forests in the study area are dominated by secondary or human-made forests with low coverage due to the destruction of vegetation caused by soil erosion and economic development in the past [35] (Figure 4). Soil types in the study area are red soil (78.4%) and paddy soil (21.6%) which can be classified as Ultisols and Inceptisols in US Soil Taxonomy, respectively [36]. Red soil is mainly distributed in the hilly region, while paddy soil is mainly distributed in the broad alluvial valleys with a similar spatial pattern of paddy rice landuse (Figure 4). The climate belongs to the mid-subtropical monsoon moist climate with an annual average temperature of 18.3 ◦C. The annual average precipitation is 1697.0 mm [35].

**Figure 3.** Map of the Youwuzhen watershed in Fujian Province, China (adapted from [2]).

**Figure 4.** Map of landuse in the study area.

The basic spatial data collected for watershed modeling of the Youwuzhen watershed include a gridded Digital Elevation Model, a soil type map, and a landuse type map, which were all unified to be of 10 m resolution [2]. Soil properties were derived from field sampling data [35]. Landuse/landcover related parameters were referenced from the SWAT database [37] and relevant literature [38]. The climate data containing daily meteorological data and precipitation data from 2012 to 2015 were derived from the National Meteorological Information Center of China Meteorological Administration and local monitoring stations, respectively. The periodic site-monitoring streamflow and sediment discharge data of the watershed outlet from 2012 to 2015 were provided by the Soil and Water Conservation Bureau of Changting County, Fujian province, China.

With an accumulated threshold of 0.185 km2 for the study area [35], the Youwuzhen watershed was delineated into 17 subbasins (Figure 3). The streamflow and sediment discharge data were screened by a rule that requires complete records of rainstorms with more than three consecutive days for watershed modeling due to the limited data quality [2]. Finally, the year 2012 was selected as a warm-up period for watershed modeling, the years 2014 and 2015 for calibration, and the year 2013 for validation.

#### *2.3. Delineation of BMP Configuration Units*

#### 2.3.1. HRUs

The QSWAT, an open source user interface for the SWAT model [39], was used to delineate the typical spatially discrete HRUs by overlaying the landuse map, soil map, and the classification map of the slope percentage. The slope percentage was classified into five classes with nearly equal areas according to the quantile classification method, i.e., 0%–9%, 9%–23%, 23–36%, 36%–48%, and larger than 48%. Three threshold values introduced in SWAT to delineate HRUs were all set to 0%, so as to obtain a full spatial coverage of HRUs (Figure 5a). Finally, a total of 355 HRUs were generated in the study area (Figure 5a).

**Figure 5.** Delineations of BMP configuration units of the Youwuzhen watershed: (**a**) HRUs; (**b**) spatially explicit HRUs; (**c**) hydrologically connected fields; (**d**) slope position units.

#### 2.3.2. Spatially Explicit HRUs

The method proposed by Teshager et al. [26] was adopted to delineate the spatially explicit HRUs. First, the landuse map was split by river and road networks—not forest type—in the map. Then, the split landuse map was intersected by the subbasin boundary. After eliminating small polygons, landuse polygons were re-labeled by assigning different codes for the polygons with the same landuse type within a subbasin. For example, three polygons with a landuse of orchard (ORCD) within a subbasin were re-assigned to ORCD1, ORCD2, and ORCD3, respectively. Similarly, the processed landuse map was intersected by the soil map and the landuse-soil polygons with the same soil type within a subbasin were re-labeled according to soil types. A single slope percentage class was implicitly

used in this method to avoid HRU fragmentation and ensure that each HRU matched individual fields in a subbasin. Finally, 166 spatially explicit HRUs were generated (Figure 5b).

#### 2.3.3. Hydrologically Connected Fields

The basic idea of the delineation algorithm of the hydrologically connected field proposed by Wu et al. [5] is to build a gridded cell tree structure based on flow directions and then merge gridded cells with the same landuse type in this tree structure into one field. A threshold value of the minimum size of a field was set to eliminate tiny fields by aggregating them into their downslope large ones. The smaller the threshold value, the more the fields will be delineated. The threshold value was insensitive to the cost-effectiveness of the optimization of BMP scenarios according to the sensitivity analysis conducted by Wu et al. [5]. Therefore, under the consideration of optimizing efficiency and comparability with other BMP configuration units, 103 hydrologically connected fields were delineated with a threshold value of 70 cells (Figure 5c).

#### 2.3.4. Slope Position Units

The same as in the case study of Qin et al. [2], a simple system of three types of slope positions (i.e., ridge, backslope, and valley) was adopted to delineate slope position units. The automated program developed by Zhu et al. [31] was used to derive fuzzy memberships of gridded cells to each slope position and then a crisp classification map of slope positions was generated by the maximization principle [30]. The slope position map was then intersected by the hillslope boundary to ensure each hillslope had a full sequence of slope position units from the top to the bottom of the hillslope. Finally, 105 slope position units within 35 hillslopes were delineated.

#### *2.4. BMP Knowledge Base and BMP Configuration Strategies*

As shown in Table 1, four BMPs that had been widely implemented in Changting County for ecological restoration and soil and water conservation were selected in this study, i.e., Closing measures (CM), Arbor-bush-herb mixed plantation (ABHMP), Low-quality forest improvement (LQFI), and Orchard improvement (OI) [2,40].

The BMP knowledge base used in this study mainly includes two categories of knowledge, i.e., the environmental effectiveness and cost-benefit knowledge used for evaluating BMP scenarios, and the BMP configuration knowledge used with BMP configuration strategies. The four BMPs considered in this study can improve soil properties through long- or/and short-term processes and thus achieve environmental effectiveness (i.e., improving water conservation and mitigating soil erosion) [39]. Therefore, the environmental effectiveness of the four BMPs considered in this study can be mainly represented by the improvements of soil properties and the change of USLE (Universal Soil Loss Equation) factors in the area with these BMPs in the watershed model [2,5,8]. The Fujian Soil and Water Conservation Monitoring Station monitored the dynamic changes in the soil properties such as organic matter, bulk density, and total porosity improved by various BMPs by taking samples annually from 2000 to 2008 [41]. This study assumed that the long-term environmental effectiveness of BMPs can reach relative stability after several years of maintenance from their first establishment [42]. Therefore, the relative changes of the 8-year sampled and derived soil properties (Table 2) were used to represent the environmental effectiveness of BMPs in the watershed modeling. Besides, the relative changes of the conservation practice factor of USLE (i.e., USLE\_P) in Table 2 were adopted from the calibrated SWAT model for this area [35]. The BMP cost-benefit data including initial implementation cost, annual maintenance cost, and annual benefit, were estimated from reports of local government projects [2] (Table 2).


**Table 1.** Brief descriptions of the four BMPs considered in this study (adapted from [2,40]; photos from [35]).

BMP configuration knowledge is normally in the form of rules. They include knowledge on individual BMPs (such as its suitable landuse types, suitable slope positions, and overall environmental effectiveness grade; Table 3) according to the characteristics of individual BMPs [35,41] (Table 1), and knowledge on spatial relationships between BMPs (such as the upstream–downstream relationships between BMP configuration units, see below) [2,5]. The overall environmental effectiveness grade of a BMP ranges from 1 to 5 which represents the improvement degree of mitigating soil erosion for an area with the BMP (the higher the better; Table 3). The overall environmental effectiveness grade can be used to formalize the expert knowledge on the spatial configuration of BMPs along the hillslope scale [2].


**Table 2.** Environmental effectiveness and cost-benefit knowledge of the four BMPs (CM = Closing measures, ABHMP = Arbor–bush–herb mixed plantation, LQFI = Low-quality forest improvement, OI = Orchard improvement) (adapted from [2]).

<sup>1</sup> Environmental effectiveness of BMPs includes soil property parameters (i.e., OM = organic matter, BD = bulk density, PORO = total porosity, and SOL\_K = soil hydraulic conductivity) and universal soil loss equation (USLE) factors (i.e., USLE\_K = soil erodibility factor and USLE\_P = conservation practice factor). Values in this column represent relative changes (i.e., multiplying) to the original property values and thus have no units.

**Table 3.** Suitable landuse types/slope positions and the overall environmental effectiveness grade of the four BMPs (CM = Closing measures, ABHMP = Arbor–bush–herb mixed plantation, LQFI = Low-quality forest improvement, OI = Orchard improvement) (adapted from [2]).


According to knowledge (or rules) used for BMP configuration, four BMP configuration strategies were adopted during assessment of the effects of four types of BMP configuration units for BMP scenarios optimization (Figure 2).


#### *2.5. Watershed Model and BMP Scenario Cost Model*

Note that BMP configuration units are not necessary to be consistent with the basic simulation unit (e.g., gridded cell) of watershed models. To simulate the spatial interactions between spatially explicitly distributed BMPs effectively, a fully distributed watershed model was constructed based on an open-source, modular, and parallelized watershed modeling framework, Spatially Explicit Integrated Modeling System (SEIMS, https://github.com/lreis2415/SEIMS; [33]) in this study. With the flexible modular structure and the parallel-computing middleware [33,43], SEIMS allows users to add their own algorithms in a nearly serial programming manner and to customize parallelized watershed models according to the characteristics of the study area and the application requirements. SEIMS also supports model-level parallel computation for applications which need numerous model runs, such as BMP scenarios optimization in this study.

To evaluate the long-term effects of BMP scenarios on mitigating soil erosion, a SEIMS-based watershed model with gridded cells as the basic simulation unit was constructed to simulate hydrology, soil erosion, plant growth, and nutrient cycling processes at a daily time-step [2]. The hydrology processes considered in this watershed model include interception, surface depression storage, surface runoff, potential evapotranspiration, percolation, interflow, groundwater flow, and channel flow. Soil erosion on hillslopes was estimated by the Modified Universal Soil Loss Equation (MUSLE) [44], and sediment routing in channels was simulated by a simplified Bagnold stream power equation adapted in the SWAT model. The plant growth process and nutrient (i.e., nitrogen and phosphorous) cycling process were also adapted from SWAT. More details about the corresponding simulation algorithms can be found in Qin et al. [2].

Note that the BMP module of SEIMS applies BMP knowledge to update the input parameters on each of the gridded cells configured with BMP suitable for current landuse on the cell before the simulation, while the non-effective configurations are ignored. The non-effective configurations may occur in two circumstances: (1) the BMP configuration unit includes multiple landuse types (e.g., the generalized hydrologically connected fields and slope position units) and some of them with small areas are unsuitable for the currently configured BMP; (2) the BMP configuration unit applied with the random configuration strategy, e.g., the BMP configuration unit with paddy rice landuse configured with Low-quality forest improvement (LQFI) or Orchard improvement (OI) practices.

After conducting the parameter sensitivity analysis using the Morris screening method [45], the SEIMS-based watershed model in the study area was calibrated by an auto-calibration procedure based on the NSGA-II algorithm (a multi-objective optimization algorithm extended from the Genetic Algorithm; [34]) provided in SEIMS [33]. The Morris screening method is a so-called one-step-at-a-time (OAT) global sensitivity analysis method, which was used to qualitatively identify important parameters for the simulation of streamflow and sediment export at the outlet in this study. Then a small number of sensitivity parameters (such as the baseflow exponent and the baseflow recession coefficient for groundwater and soil water capacity) were selected for auto-calibration in this study [33]. The NSGA-II was also applied to BMP scenarios optimization in this study, and is described in detail in Section 2.6. One optimal calibration solution was selected as the baseline scenario (Table 4). The calibration and validation results of streamflow (m<sup>3</sup> s<sup>−</sup>1) and sediment export (kg) at the outlet of Youwuzhen watershed were evaluated by widely-used model performance indicators such as NSE (Nash–Sutcliffe Efficiency), PBIAS (Percent BIAS), and RSR (Root mean Square error–standard deviation Ratio) [46]. According to the general performance ratings for simulations at a monthly time step [46], the model performance is satisfactory when NSE ≥0.50, RSR ≤0.70, and the absolute value of PBIAS ≤25% (55% for sediment). Thus, the streamflow performance is nearly satisfactory, while the sediment is relatively poor because of a few unreasonable peak values (Table 4). Nevertheless, the general trends of hydrographs in the study area can be captured by the calibrated SEIMS-based model according to visual judgement. Therefore, it is acceptable to apply the calibrated SEIMS-based watershed model to BMP scenarios optimization. The annual average sediment yields from 2013 to

2015 of the entire watershed, under each BMP scenario for spatial optimization, were calculated by this model.

**Table 4.** The calibration (2014–2015) and validation (2013) results of the baseline scenario by the Spatially Explicit Integrated Modeling System (SEIMS)-based watershed model. (NSE = Nash–Sutcliffe Efficiency, RSR = Root mean Square error–standard deviation Ratio, PBIAS = Percent BIAS).


Besides the SEIMS-based watershed model for evaluating environmental effectiveness, a simple BMP scenario cost model (Equation (1)) was adopted to calculate the net cost of each BMP scenario according to the cost-benefit knowledge in the BMP knowledge base [2].

$$f\_{\text{net}-\text{cost}}(X) = \sum\_{i=1}^{n} A(\mathbf{x}\_{i}) \times \left\{ \left[ \mathbb{C}(\mathbf{x}\_{i}) + yr \times \left( M(\mathbf{x}\_{i}) - B(\mathbf{x}\_{i}) \right) \right] \right\} \tag{1}$$

where *f* net-cost(*X*) is the net cost of a BMP scenario (represented as *X*); *n* is the count of BMP configuration units; *A*(*x*i) is the area covered by the BMPs implemented in the *i* th configuration unit; *yr* is the years when the effectiveness of BMPs reach stability, which is 8 in this study (see Section 2.4); *C*(*x*i), *M*(*x*i), and *B*(*x*i) are unit costs (CNY 10,000/km2) for initial implementation, annual maintenance, and annual benefit (Table 2), respectively.

#### *2.6. Multi-Objective Optimization by an Intelligent Optimization Algorithm*

The objectives in this study are to minimize the net cost of the BMP scenario (Equation (2)) and maximize the reduction rate of soil erosion. The reduction rate of soil erosion of each BMP scenario is the relative change compared to the baseline scenario (Equation (3)).

$$\text{optimal solutions} = \min\left(-f\_{\text{reduction}-\text{rate}}(X), f\_{\text{net}-\text{cost}}(X)\right) \tag{2}$$

$$f\_{\text{reduction}-\text{rate}}(X) = (v(0) - v(X)) / v(0) \tag{3}$$

where *f*reduction-rate(*X*) is the reduction rate of soil erosion under the BMP scenario *X* compared to that under the baseline scenario; *v*(0) and *v*(*X*) are the total amount of soil erosion (kg) under the baseline scenario and the *X* scenario, respectively.

The NSGA-II [33], which has been successfully applied to many similar studies [2,16–18], was selected as the intelligent optimization algorithm in this study. As shown in Figure 2, the NSGA-II algorithm first initializes the initial BMP scenarios (called "population" in NSGA-II) based on one type of BMP configuration unit and one available BMP configuration strategy as described in Sections 2.3 and 2.4. A BMP scenario (called "individual") is represented as an array with a length equal to the number of BMP configuration units (called "chromosome") with the corresponding reduction rate of soil erosion and net cost values (called "fitness") to be evaluated. Each value of the chromosome (called "gene") stands for one selected BMP type or without BMP on a BMP configuration unit. Then, each individual of the initial population, as the initial "generation" for the following optimization, is evaluated by objective functions, i.e., the reduction rate of soil erosion is assessed by the calibrated watershed model and the net cost estimated by the BMP scenario cost model. The evaluated individuals are sorted, based on the non-domination of fitness, and selected by a specified number as an elite set for each generation which is known as near-optimal Pareto solutions [33]. Then a circular process of regenerating and evaluating BMP scenarios (i.e., new generation) proceeds until a user-assigned maximum generation number is reached. Individuals of current generation (called "offspring") combined with the near-optimal Pareto solutions of former generation (called "parent") are proceeded by non-dominated sorting, selection, crossover, and mutation operations to regenerate new BMP scenarios [33].

When a BMP configuration strategy with knowledge (i.e., SUIT strategy, UPDOWN strategy, or HILLSLP strategy) was adopted, it was incorporated into not only the initialization [5] but also the regeneration of BMP scenarios, i.e., crossover and mutation operations [2]. The crossover operation of the UPDOWN strategy proposed by Wu et al. [5] was extended as follows. The randomly selected BMP configuration unit (i.e., the position of the crossover gene) is first checked to see whether it can ensure that the two generated "children" individuals still conform to the "if–then" rules after the exchange of the subtrees in which the selected unit is the most downstream unit. If the current selected unit fails to meet the condition, the downstream units of the current one will be checked in order, until a qualified unit is reached. The finally reached unit, except for the most downstream unit of the entire study area, will be used as the crossover gene. Under the HILLSLP strategy, the BMPs are configured along each hillslope from bottom to top during the initialization, and the crossover operation is to exchange the randomly selected hillslopes (without breaking the configuration among all slope position units along the same hillslope). Such generated "children" individuals conform to the HILLSLP strategy. As for the mutation operation of the UPDOWN strategy and HILLSLP strategy, the potential BMP types for a BMP configuration unit (i.e., a mutant gene) are first determined according to the rule set of adopted knowledge and the BMP types of its upstream and downstream units (i.e., gene values). Then, a different BMP from the current one will be randomly selected for the mutant gene. In such a way, every BMP scenario generated during spatial optimization is reasonable in terms of the corresponding knowledge, and may result in higher optimizing efficiency.

#### *2.7. Design of Comparison Experiment*

In this study, four types of BMP configuration units (i.e., HRUs, spatially explicit HRUs (EXPLICITHRU), hydrologically connected fields (CONNFIELD), and slope position units (SLPPOS)) and four BMP configuration strategies (i.e., RAND strategy, SUIT strategy, UPDOWN strategy, and HILLSLP strategy) were combined according to their availability (see Section 2.4) for the spatial optimization of multiple spatially explicit BMPs. This means that, in total, 11 experiments of feasible combinations were conducted, i.e., HRU+RAND which means the combination of using HRU as BMP configuration units applied with the random configuration strategy (ex analogia), HRU+SUIT, EXPLICITHRU+RAND, EXPLICITHRU+SUIT, CONNFIELD+RAND, CONNFIELD+SUIT, CONNFIELD+UPDOWN, SLPPOS+RAND, SLPPOS+SUIT, SLPPOS+UPDOWN, and SLPPOS+HILLSLP. The Python script of the BMP scenarios optimization based on slope position units developed by Qin et al. [2] was extended to support multiple types of BMP configuration units and BMP configuration strategies considered in this study. The script built on the model-level parallel framework of SEIMS [32] distributes computing tasks dynamically across a Linux cluster to improve computation efficiency.

The parameter settings of the NSGA-II algorithm were kept the same for all experiments. The initial population size was 480 with a selection rate of 0.8 and a maximum generation number of 100 [5]. The crossover probability and the mutation probability were 0.8 and 0.1, respectively.

To explore the effects of different BMP configuration units for the spatial optimization of spatially explicit BMPs, the results of different BMP configuration units applied with the same BMP configuration knowledge are first compared. Then, the combinations of BMP configuration units and the corresponding optimal configuration strategies are compared. The optimal configuration strategy for a type of BMP configuration unit was considered to be the most reasonable one based on a BMP knowledge base and not necessarily the one with the best non-dominated Pareto solutions from the mathematical perspective. Thus, SLPPOS+HILLSLP, CONNFIELD+UPDOWN, EXPLICITHRU+SUIT, and HRU+SUIT were selected for this comparison according to the available level with the most BMP configuration knowledge for each BMP configuration unit (Section 2.4).

The results of different BMP configuration units will be discussed from the perspectives of cost-effectiveness, optimizing efficiency, and practicality. The near-optimal Pareto solutions plotted as a scatter plot can give a simple and direct interpretation of the convergence and diversity of different spatial optimization results. When the near-optimal Pareto solutions from different combinations are compared, the one with more non-dominated solutions indicates a better cost-effectiveness [33], which means more BMP scenarios with better multi-objective optimization provided as optimized solutions for decision-making. Besides, the hypervolume index [47], which measures the volume (or area for two-dimensions) of objective space covered by a set of near-optimal Pareto solutions, provides a quantitative comparison of the cost-effectiveness considering both convergence and diversity [48]. A higher hypervolume index indicates a better quality of solution.

The changes of the hypervolume index with generations can provide a qualitative estimation of the optimizing efficiency. For an ideal optimization, the hypervolume index will increase rapidly at the beginning of the optimization, then increase slowly, and eventually remain stable. Therefore, the faster the hypervolume index reaches stability, the better the optimizing efficiency. For convenience and consistency of comparison, a criterion is adopted to judge whether the hypervolume index has reached stability in this study, i.e., if the increment rate of the hypervolume index compared to the former generation is lower than 0.1% for three consecutive generations and there are also no three consecutive increment rates greater than 0.1% in the following generations. For comparability in this study, the worst reference point for calculating the hypervolume index of all experiments was set to (300, 0), which represents the net cost being CNY 3 million and the reduction rate of soil erosion being zero.

Note that both the near-optimal Pareto solutions and hypervolume index represent evaluations from a mathematical perspective, which have less practical meaning than the practicality of the spatial distribution of BMP scenarios for decision-making in integrated watershed management [2]. Therefore, the practicality of the spatial distributions of selected near-optimal Pareto solutions from different BMP configuration units with the corresponding optimal BMP configuration strategy will be qualitatively discussed based on the visual interpretation and local experiences of the study area.

#### **3. Results**

#### *3.1. Comparison among Different BMP Configuration Units with the RAND Strategy*

Figure 6 shows the near-optimal Pareto solutions of the 100th generation (Figure 6a) and the hypervolume index with generations (Figure 6b) from four types of BMP configuration units applied with the random configuration strategy (i.e., SLPPOS+RAND, CONNFIELD+RAND, EXPLICITHRU+RAND, and HRU+RAND). Effective but different near-optimal Pareto solutions are generated by all combinations in almost the same solution space, i.e., soil erosion reduction rates range from 0.10 to 0.50 with the net cost range from CNY 0.03 to 1.50 million (Figure 6a). The near-optimal Pareto solutions of SLPPOS+RAND and EXPLICITHRU+RAND are nearly overlapped, while HRU+RAND produced the most non-dominated solutions at the nearly entire solution space and CONNFIELD+RAND was dominated by the other three combinations. HRU+RAND achieved the best overall performance considering convergence and diversity, compared to the other three combinations with the RAND strategy (Figure 6a). All four of these combinations obtained approximately the same values of the hypervolume index with generations, which showed a similar change trend, i.e., the hypervolume index increased rapidly for about the first 20 generations (e.g., with increment rates of the hypervolume index greater than 1%) and then increased slowly until stability was reached at the 46th, 41th, 56th, and 68th generations for SLPPOS+RAND, CONNFIELD+RAND, EXPLICITHRU+RAND, and HRU+RAND, respectively (Figure 6b), among which CONNFIELD+RAND showed the best optimizing efficiency.

**Figure 6.** Comparisons among four types of BMP configuration units applied with the random configuration (RAND) strategy: (**a**) near-optimal Pareto solutions of the 100th generation and, (**b**) the hypervolume index with generations.

As the BMP configuration strategy remained the same, the characteristics of BMP configuration units are the main causes of the differences among the optimization results. Compared to the other three types of configuration units, HRUs in the study area obtained the largest count of units and the most detailed spatial delineation. This means that under the RAND strategy, a small number of HRUs configured with BMPs (i.e., at a low net cost) at some critical erosion zone may result in a relatively higher soil erosion reduction rate (i.e., better cost-effectiveness). The generalized hydrologically connected fields with upstream–downstream relationships at the watershed scale usually have large areas and may occupy most of the upslope positions of subbasins, where the critical erosion zone is often located. Therefore, CONNFIELD+RAND performed less effectively at low net costs than the other three combinations and obtained relatively stable high soil erosion reduction rates with relative greater net costs (Figure 6a). Since, in the study area, both landuse and soil types have a similar spatial pattern to the topography (Section 2.2), and spatially explicit HRUs were overlaid by landuse, soil types, and subbasin boundaries, EXPLICITHRU can represent the spatial distribution of topographic characteristics to some degree (Figure 5b). Considering that the slope position units were totally delineated according to topographic characteristics, the similarity between these two spatial configuration units maybe the reason for the quite similar results between EXPLICITHRU+RAND and SLPPOS+RAND. With consistent crossover and mutation operations during optimization, the hypervolume index trends with generations are very similar among four types of BMP configuration units applied with the RAND strategy (Figure 5b). It might be inferred that the optimizing efficiency under the RAND strategy is negatively correlated with the count of BMP configuration units.

#### *3.2. Comparison among Different BMP Configuration Units with the SUIT Strategy*

Figure 7 shows a comparison among four types of BMP configuration units applied with the SUIT strategy (i.e., SLPPOS+SUIT, CONNFIELD+SUIT, EXPLICITHRU+SUIT, and HRU+SUIT). The solution spaces of the spatial optimization of EXPLICITHRU+SUIT and HRU+SUIT were inclined to reach higher net costs, while SLPPOS+SUIT and CONNFIELD+SUIT comparatively concentrated on the solution spaces with low net costs (Figure 7a). HRU+SUIT and SLPPOS+SUIT obtained the most non-dominated solutions when the corresponding net costs were greater or less than about CNY 1.10 million, respectively. For HRU+SUIT and EXPLICITHRU+SUIT, the soil erosion reduction rate reached relative stability at about 0.52 when the corresponding net costs were greater than CNY 1.60 million (Figure 7a). The solutions of SLPPOS+SUIT obtained very similar hypervolume index values to those from CONNFIELD+SUIT in the first stage of optimization (about the first 35 generations) and then a higher hypervolume index during the following generations (Figure 7b). The EXPLICITHRU+SUIT and HRU+SUIT produced worse diversity than SLPPOS+SUIT and CONNFIELD+SUIT (Figure 7a). This phenomenon can also be observed from the lower hypervolume index from EXPLICITHRU+SUIT and HRU+SUIT (Figure 7b). The SLPPOS+SUIT, CONNFIELD+SUIT, EXPLICITHRU+SUIT, and HRU+SUIT reached stability at the 46th, 36th, 45th, and 60th generations, respectively (Figure 7b), among which the CONNFIELD+SUIT showed the best optimizing efficiency.

**Figure 7.** Comparisons among four types of BMP configuration units applied with the suitable landuse types/slope positions of individual BMPs (SUIT) strategy: (**a**) near-optimal Pareto solutions of the 100th generation, and (**b**) the hypervolume index with generations.

Since the HRUs and spatially explicit HRUs delineated in this study have a single landuse type for each spatial configuration unit, the non-effective configurations caused by BMP configuration units with multiple landuse types (e.g., hydrologically connected fields and slope position units) (Section 2.5) do not exist. Therefore, under the same parameter-settings of the optimization algorithm, EXPLICITHRU+SUIT and HRU+SUIT inherently generate more effective BMP configurations in locations within spatial configuration units than SLPPOS+SUIT and CONNFIELD+SUIT and thus result in a solution space with relatively high net costs.

With the best near-optimal Pareto solutions, the highest hypervolume index, and the satisfied optimizing efficiency, SLPPOS+SUIT obtained the best overall performance compared to the other three combinations with the SUIT strategy. This may be attributed to the knowledge used by the SUIT strategy for slope position units that considers not only suitable landuse types but also suitable slope positions of individual BMPs.

#### *3.3. Comparison between Feasible BMP Configuration Units with the UPDOWN Strategy*

Figure 8 presents a comparison between two types of BMP configuration units applied with the UPDOWN strategy (i.e., SLPPOS+UPDOWN and CONNFIELD+UPDOWN). CONNFIELD+UPDOWN and SLPPOS+UPDOWN obtained most of their non-dominated solutions when the corresponding net costs were greater or less than about CNY 0.25 million, respectively (Figure 8a). SLPPOS+UPDOWN had a narrower solution space than CONNFIELD+UPDOWN and hence had worse diversity and a lower hypervolume index. CONNFIELD+UPDOWN showed a more steady hypervolume index trend with generations and better optimizing efficiency since it reached stability at the 39th generation which is lower than the 55th generation for SLPPOS+UPDOWN (Figure 8b). SLPPOS+UPDOWN reached slightly better hypervolume index stability after about the 70th generation, i.e., a slightly better convergence than CONNFIELD+UPDOWN.

**Figure 8.** Comparisons between two types of BMP configuration units applied with the UPDOWN strategy: (**a**) near-optimal Pareto solutions of the 100th generation, and (**b**) the hypervolume index with generations.

Although the same strategy was applied, the upstream–downstream relationships of slope position units were built at the hillslope scale instead of the watershed scale of hydrologically connected fields. This means that the number of slope position units within one hillslope (i.e., three in this study) is the maximum number of genes allowed to be exchanged during crossover operation in NSGA-II (Section 2.6). This puts the SLPPOS+UPDOWN result in a slightly better convergence and worse diversity than CONNFIELD+UPDOWN. The non-dominated solutions with low net costs of SLPPOS+UPDOWN indicated the better cost-effectiveness of slope position units under a tight budget, while its worse cost-effectiveness with higher net costs than CONNFIELD+UPDOWN may imply that the UPDOWN strategy initially developed for hydrologically connected fields [5] is not the most effective strategy for slope position units.

#### *3.4. Comparison among Different BMP Configuration Units with the Corresponding Optimal Configuration Strategies*

As shown in Figure 9a, HRU+SUIT and SLPPOS+HILLSLP generated the most effective non-dominated solutions when the corresponding net costs were greater and less than about CNY 1.0 million, respectively. According to the hypervolume index (Figure 9b), SLPPOS+HILLSLP reached the stable hypervolume index after the 37th generation, as well as the largest hypervolume index value, compared to the other three combinations. Thus, SLPPOS+HILLSLP showed the best optimizing efficiency. In summary, SLPPOS+HILLSLP showed the best overall performance (Figure 9), followed by EXPLICITHRU+SUIT and CONNFIELD+UPDOWN. Although HRU+SUIT had non-dominant solutions mostly at high net costs, it was still considered to have the worst overall performance because it produced the worst diversity, lowest hypervolume index, and slowest optimizing efficiency (i.e., reaching stability at the 60th generation; Section 3.2).

**Figure 9.** Comparisons among four types of BMP configuration units applied with the corresponding optimal configuration strategies: (**a**) near-optimal Pareto solutions of the 100th generation, and (**b**) the hypervolume index changed with generations.

To compare the practicality of the spatial distribution of optimized BMP scenarios, four BMP scenarios (Figure 10) were selected in the overlapped solution space of the four combinations (i.e., at a soil erosion reduction rate of around 0.40 and net cost of CNY 0.70–1.20 million). The most fragmented spatial delineation of HRUs caused the worst practicality from the perspective of actual watershed management (Figure 10a). EXPLICITHRU+SUIT (Figure 10b) generated similar BMPs regions to HRU+SUIT (Figure 10a), but those from EXPLICITHRU+SUIT had a more concentrated and practical distribution of BMPs than those from HRU+SUIT. The BMP scenarios from CONNFIELD+UPDOWN required the highest net cost to achieve the same soil erosion reduction rate among these selected BMP scenarios, since CONNFIELD+UPDOWN selected the most expensive BMP (i.e., Orchard Improvement) for allocation (Figure 10c). With the adoption of the UPDOWN strategy and a smaller configuration area, the BMP scenario from CONNFIELD+UPDOWN (Figure 10c) achieved better practicality than that from EXPLICITHRU+SUIT (Figure 10b). With the SUIT and UPDOWN strategies, the lack of precise spatial relationships between BMPs and spatial positions at the hillslope scale may cause several inappropriate configurations (Figure 10a–c), e.g., the LQFI (Low-quality forest improvement) configured on ridge areas (Table 1), thus reducing the practicality of the BMP scenario from EXPLICITHRU+SUIT and CONNFIELD+UPDOWN. Depending on the HILLSLP strategy used to adopt BMP knowledge derived from the local experience of integrated watershed management, SLPPOS+HILLSLP concisely and precisely configured CM (Closing measures) and ABHMP (Arbor–bush–herb mixed plantation) on slope position units at the hillslope scale, e.g., CM-ABHMP on the ridge and backslope sequence, and also configured the BMP with the best overall effectiveness grade (i.e., ABHMP according to Table 3) on ridge and backslope. Thus, SLPPOS+HILLSLP obtained the best practicality with the lowest net cost among these selected BMP scenarios, followed by CONNFIELD+UPDOWN and EXPLICITHRU+SUIT.

**Figure 10.** Comparison of the spatial distribution of selected BMP scenarios of the 100th generation from four types of BMP configuration units applied with the corresponding optimal strategies: (**a**) one HRU+SUIT scenario with a soil erosion reduction rate of 0.41 and a net cost of CNY 0.98 million; (**b**) one EXPLICIT+SUIT scenario with a soil erosion reduction rate of 0.40 and a net cost of CNY 0.97 million; (**c**) one CONNFIELD+UPDOWN scenario with a soil erosion reduction rate of 0.40 and a net cost of CNY 1.08 million; and (**d**) one SLPPOS+HILLSLP scenario with a soil erosion reduction rate of 0.40 and a net cost of CNY 0.78 million.

#### **4. Discussion**

From a mathematical viewpoint, all feasible combinations of BMP configuration units and BMP configuration strategies generated effective near-optimal Pareto solutions (Figures 6–9). Different delineations of BMP configuration units have different characteristics such as the number of units, the spatial distribution characteristics, and the spatial relationships with homogeneous functional units from the perspective of physical geography. These differences affect the characteristics of the generated BMP scenarios and the search spaces of spatial optimization, thus resulting in the different cost-effectiveness of near-optimal Pareto solutions and optimizing efficiency when applied with the same BMP configuration strategy (Sections 3.1–3.3).

BMP scenarios generated based on different kinds of BMP configuration knowledge had significant differences in the cost-effectiveness of the near-optimal Pareto solutions, optimizing efficiency, and the spatial distribution of the BMP scenarios. Using the strategies that adopt BMP configuration knowledge (i.e., the SUIT, UPDOWN, and HILLSLP strategies), those non-effective configurations generated by the random configuration strategy (Section 2.5) can be avoided for all BMP configuration units. Thus, with the BMP configuration constraint, better convergence, and worse diversity of near-optimal Pareto solutions were obtained, as well as overall lower values of hypervolume index and better optimizing efficiency (e.g., comparing Figure 7, Figure 8, and Figure 9 with Figure 6). Such phenomena could also be observed when the adopted knowledge extended from simple knowledge to knowledge on the spatial relationships between BMPs and spatial positions, i.e., from the SUIT strategy to the UPDOWN strategy or the HILLSLP strategy.

Besides the narrowed solution space of spatial optimization and improved optimizing efficiency, the adoption of domain (or expert) knowledge considering spatial relationships between BMPs and spatial positions, e.g., the UPDOWN and HILLSLP strategies in this study, can also facilitate the optimized BMP solutions with geographical and practical meanings (Section 3.4). The formal representation of such expert knowledge is often introduced according to the characteristics of specific BMP configuration units, which means that one BMP configuration strategy developed associated with one BMP configuration unit may not be effectively applied to another BMP configuration unit, e.g., the ineffective situation in which the UPDOWN strategy was applied to slope position units (Figures 8 and 9).

The random strategy and simple knowledge-based strategy (e.g., the SUIT strategy) may obtain more effective near-optimal BMPs solutions at higher net costs than the strategies with knowledge on the spatial relationships between BMPs and spatial positions (e.g., the UPDOWN and HILLSLP strategies) (Figure 9a). However, the latter type of strategies can effectively represent the local experience of integrated watershed management with better practicality [2] and hence are more valuable for real applications, especially when the expert knowledge on the spatial relationships between BMPs and slope positions along a hillslope can be considered (Figure 10d).

#### **5. Conclusions**

This paper presents a comparison among the four main types of BMP configuration units (i.e., HRUs, spatially explicit HRUs, hydrologically connected fields, and slope position units) for watershed BMP scenarios optimization based on the same one distributed watershed model. Four BMP configuration strategies, adopting different kinds of BMP knowledge of the study area, were considered, i.e., the random configuration strategy, the strategy with knowledge on suitable landuse types/slope positions of individual BMPs, the strategy based on expert knowledge of upstream–downstream relationships [5], and the strategy with expert knowledge on the spatial relationships between BMPs and slope positions along the hillslope [2]. A total of 11 experiments of feasible combinations were conducted with the same optimization algorithm (i.e., NSGA-II) and then compared.

The comparison showed that different BMP configuration units applied with different configuration strategies had significant differences in near-optimal Pareto solutions, optimizing efficiency, and spatial distribution of BMP scenarios. Generally, the more the expert (or domain) knowledge was considered, the better the comprehensive cost-effectiveness and practicality of the optimized BMP scenarios, and the better the optimizing efficiency. Therefore, BMP configuration units that support the adoption of expert knowledge on the spatial relationships between BMPs and spatial locations [2,5] (i.e., hydrologically connected fields, and slope position units) are considered to be the most valuable spatial configuration units for watershed BMP scenarios optimization and integrated watershed management. Overall, using the slope position units as BMP configuration units with the HILLSLP strategy, the best comprehensive results of BMP scenarios optimization were obtained.

This study provided a useful reference for the spatial optimization of watershed BMP scenarios with multiple spatially explicit BMPs. For those spatially explicit BMPs, more BMP configuration knowledge derived from local management experiences should be summarized and adopted together with the feasible BMP configuration units for such knowledge (e.g., slope position units) during watershed BMP scenarios optimization.

**Author Contributions:** Conceptualization, L.-J.Z., C.-Z.Q., and A.-X.Z.; Methodology, L.-J.Z. and C.-Z.Q.; Software, L.-J.Z., J.L., and H.W.; Investigation, L.-J.Z.; Resources, J.L., C.-Z.Q., and A.-X.Z.; Writing—Original Draft Preparation, L.-J.Z. and C.-Z.Q.; Funding Acquisition, C.-Z.Q., and A.-X.Z.

**Funding:** The work reported was financially supported by the National Natural Science Foundation of China (No. 41871362 and 41431177), the National Key Technology R&D Program (No. 2013BAC08B03-4), the Chinese Academy of Sciences (No. XDA23100503), the Innovation Project of LREIS (No. O88RA20CYA), National Basic Research Program of China (Project No.: 2015CB954102), PAPD, and the Outstanding Innovation Team in Colleges and Universities in Jiangsu Province. Supports to A-Xing Zhu through the Vilas Associate Award, the Hammel Faculty Fellow Award, and the Manasse Chair Professorship from the University of Wisconsin-Madison are greatly appreciated.

**Acknowledgments:** The authors thank Zhi-Biao Chen for kindly permitting us to use his photos of the four types of BMPs considered in this study.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*

## **Evaluating the Effectiveness of Spatially Reconfiguring Erosion Hot Spots to Reduce Stream Sediment Load in an Upland Agricultural Catchment of South Korea**

#### **Kwanghun Choi 1,2,\*, Ganga Ram Maharjan 3,4 and Björn Reineking 1,5**


Received: 30 March 2019 ; Accepted: 30 April 2019; Published: 7 May 2019

**Abstract:** Upland agricultural expansion and intensification cause soil erosion, which has a negative impact on the environment and socioeconomic factors by degrading the quality of both nutrient-rich surface soil and water. The Haean catchment is a well-known upland agricultural area in South Korea, which generates a large amount of sediment from its cropland. The transportation of nutrient-rich sediment to the stream adversely affects the water quality of the Han River watershed, which supports over twenty million people. In this paper, we suggest a spatially explicit mitigation method to reduce the amount of sediment yield to the stream of the catchment by converting soil erosion hot spots into forest. To evaluate the effectiveness of this reconfiguration, we estimated the sediment redistribution rate and assessed the soil erosion risk in the Haean catchment using the daily based Morgan–Morgan–Finney (DMMF) model. We found that dry crop fields located in the steep hill-slope suffer from severe soil erosion, and the rice paddy, orchard, and urban area, which are located in a comparatively lower and flatter area, suffer less from erosion. Although located in the steep hill-slope, the forest exhibits high sediment trapping capabilities in this model. When the erosion-prone crop lands were managed by sequentially reconfiguring their land use and land cover (LULC) to the forest from the area with the most severe erosion to the area with the least severe erosion, the result showed a strong reduction in sediment yield flowing to the stream. A change of 3% of the catchment's crop lands of the catchment into forest reduced the sediment yield entering into the stream by approximately 10% and a change of 10% of crop lands potentially resulted in a sediment yield reduction by approximately 50%. According to these results, identifying erosion hot spots and managing them by reconfiguring their LULC is effective in reducing terrestrial sediment yield entering into the stream.

**Keywords:** DMMF; landscape configuration; landscape ecology; hydrology

#### **1. Introduction**

Agriculture expansion and intensification often lead to severe soil erosion in the course of altering naturally dominated surface configurations [1–3]. The problem is prominent in upland agriculture areas under monsoonal climate because of the disturbed erosion-prone hill-slopes receiving intermittent concentrated heavy rainfall [4,5]. A large amount of surface runoff from heavy rainfall washes out nutrient-rich surface soil from deforested upland agriculture areas and degrades the soil

quality of the agricultural area [6]. Eroded nutrient-rich soil particles cause not only soil quality degradation of the agricultural area but also on- and off-site water deterioration when these particles enter the stream of a catchment [7–9].

The Han River watershed in South Korea experiences extreme downpours that cause severe soil erosion and subsequent water deterioration every summer monsoon season [3,10,11]. These problems are worsening, as upland agricultural areas expand and the intensity of monsoonal rainfall increase due to ongoing climate change [12,13]. The Han River is the primary freshwater source for the Seoul Metropolitan area where over 25 million inhabitants (ca. 50% of the South Korean population) reside. Therefore, soil erosion control in this region is highly relevant to provide clean and usable freshwater resources to the residents [14,15]. With increasing demand for food crops, intensive upland agriculture is expanding in the mountainous upstream regions of the Han River watershed where few agricultural activities had been performed previously [2]. The Haean catchment is one of the largest contributors to sediment in the watershed, where abrupt land use and land cover (LULC) changes have taken place on forested hill-slope areas [11,16,17]. The LULC changes on the erosion-prone hill-slopes of this catchment generate a massive amount of sediment flowing into the river system and eventually deteriorate the water quality of the Han River [3]. Various studies have been conducted in this catchment to understand the sediment redistribution patterns and determine optimal measures to mitigate this problem. Field-level studies have focused on the effect of surface configurations of the dry croplands and their field margins on sediment yields. Arnhold et al. [11] and Ruidisch et al. [16] investigated the effect of plastic mulch applied to dry croplands on surface runoff and sediment yield. Ali and Reineking [5] showed the effectiveness of natural field margin (i.e., vegetated filter strip next to the dry cropland) for preventing off-site sediment yield. They reported that the natural field margin captured sediments more efficiently under the increased rainfall and slope conditions than intensively managed field margins with less dense vegetation cover. Arnhold et al. [17] found that organic farming yielded less sediment than conventional farming because organic farming tends to protect the soil surface by preserving more vegetations that are not cultivated crops.

At the catchment level, the soil and water analysis tool (SWAT) [18] has been widely used to test the effectiveness of various best management practices (BMPs) to reduce the sediment yield under complex terrain and landscape configurations [3,19]. Maharjan et al. [3] showed the effectiveness of catchment-wide cover crop cultivation in the dry croplands to reduce suspended sediment yields entering the stream. Jang et al. [19] projected vegetation filter strip, rice straw mulching, and fertilizer control scenarios to dry croplands of the catchment and found that the application of vegetation filter strips and rice straw mulching was efficient in reducing sediment yields from the catchment. The BMPs suggested in the aforementioned studies are often premised on the compliance of each stakeholder, which is not easily accomplished [20–22]. Different from the BMP approaches relying on stakeholders participation, several studies are paying attention to the importance of the landscape and its spatial configuration, which has a significant impact on ecosystem services and functions, including soil erosion and water quality control [23–25]. Furthermore, these studies showed that ecosystem services and functions often responded non-linearly to the spatial relocation of the agricultural landscape, implying the effectiveness of spatial configuration on enhancing ecosystem services [23,24,26]. Therefore, identifying soil erosion hot spots and assessing the sediment reduction rate by altering the surface configuration of hot spots promise to help establishing cost-effective soil erosion control methods in the catchment.

To consider the spatial context of soil erosion, a spatially explicit and distributed soil erosion model that can simulate the sediment budget of each element, considering the sediment inputs from the upslope areas is needed. Among the various soil erosion models, the daily based Morgan–Morgan–Finney (DMMF) model [15] is one of the most appropriate tools because the model can project soil erosion and deposition explicitly, considering the spatial connectivity, which facilitates the assessment of the impact of the spatial context of landscape on sediment redistribution patterns. Furthermore, the DMMF is suitable for projecting under a monsoon climate, accompanying

concentrated rainfall during a short period [15]. Vegetative filter strips (VFSs) are known as an effective tool for reducing sediment yield from the field or catchment because of their cost-effective surface protecting and sediment trapping capabilities [5,19,25,27–29]. We adopt the forest, which is a type of VFS, as an alternative LULC for soil erosion hot spots to reduce the total sediment yield into the stream of the catchment. In this study, we assessed the importance of the spatial conversion of erosion hot spots into forest on soil erosion control using the spatially explicit daily based Morgan–Morgan–Finney (DMMF) soil erosion model. The detailed objectives are to:


#### **2. Materials and Methods**

#### *2.1. Study Area*

The study was conducted in the Haean catchment (Figure 1). The Haean catchment is a bowl-shaped small mountainous erosion basin (64.4 km2) located in the northeastern part of South Korea (38.277° N, 128.135° E). As an erosion basin, the central area is low and flat, and it becomes higher and steeper toward the boundary. The lowest altitude of the catchment is 339 m, and the highest one is 1321 m [2,3,11,30]. Geologically, the catchment consists primarily of two bedrocks. One is gneiss at the higher elevation near the catchment boundary, and the other is highly weathered granite at the flat central area [2,30]. Differential erosion between the two bedrocks formed the unique bowl-shaped catchment [2]. The major soil type of the catchment is cambisol from weathered granite. The dominant soil texture of the catchment is loamy sand (59.4%) followed by sandy loam (27.5%), and sand (10.5%), which has a high infiltration capacity [3,30].

The climate of the catchment is characterized by cold and dry winter, affected by the continental Siberian high, and hot and humid summer affected by the subtropical North Pacific high [30–32]. The average annual precipitation from 2009 to 2011 is 1599 mm, and almost 70% of the rainfall is concentrated in the three months from June to August [3,11,19,30]. Due to climate change, the period of rain spell, as well as the frequency and intensity of heavy rainfall, has increased in this region [33,34].

The dominant land cover type of the catchment is forest. Forest mainly covers the summit and upper hill-slope areas around the boundary of the catchment, occupying 58% of the entire catchment area. Dry croplands (22%), including bean, cabbage, potato, radish, and ginseng, dominate the lower hill-slope areas adjacent to the forest edge. Rice paddies (8%) and residential areas (3%) (e.g., roads and artificial structures) occupy the flat central area of the catchment. Semi-natural vegetation field (8%), shrublands (1%), and bare surface (5%), including fallow and barren field, cover the remaining areas [35].

The dry croplands have been expanded into the forest that is located in the hill-slope area. Due to the upland agriculture expansion after deforestation, the catchment yields a massive amount of sediment into the stream during the summer monsoon season. The sediment is transported to the Soyang reservoir. This reservoir is the largest reservoir in South Korea as well as the crucial freshwater source for citizens living in the Seoul metropolitan area [3,11,30]. Weather stations and hydrological measurement facilities are installed in the catchment to monitor the climate and stream conditions, and erosion control dams and the reservoir have been constructed to reduce the sediment yield from the catchment [30,36].

**Figure 1.** General description of the study area. Locations of the Soyang lake watershed and Haean catchment in South Korea are described in the lower left figure. In the upper right figure, the topography and stream networks of the study area, with the monitoring sites (red triangles) and weather stations (yellow circles) used for the DMMF model are presented.

#### *2.2. Model Description*

We used the DMMF model [15] to assess the soil erosion risk and simulate the impact of the spatial reconfiguration of erosion hot spots into forest on sediment yield within the Haean catchment The DMMF model was modified from the widely used Morgan–Morgan–Finney (MMF) soil erosion model [37], which has a simple structure while maintaining physical foundations [15,38–41].

The DMMF model has three significant modifications relative to the MMF model: the adoption of a daily time step, the consideration of the effect of impervious ground cover on soil erosion, and the revision of the equations and sequence of the subprocesses for a better physical representation of physical processes, such as surface runoff and sediment redistribution [15,42]. These modifications enable the model to be more suitable for estimating surface runoff and soil erosion on a complex surface terrain under an intensive seasonal rainfall regime than the previous version.

The DMMF model can estimate the amount of surface and subsurface water input from the upslope area and output to the downslope area after hydrological processes for each element (e.g., each grid cell in a raster map). The model also estimates the sediment budget of each element by calculating the amount of sediments flowing into and out of the element. The hydrological processes of the model are determined by rainfall, evapotranspiration, surface/subsurface water inflows, and initial soil water content (Figure 2). After calculating the water budget for the element, the model calculates sediment budgets, considering the amount of sediment input from the upslope areas, rainfall intensity, topography, soil characteristics, surface configurations, and vegetation structures (Figure 3). The detailed input parameters are presented in Table 1 and detailed structure and equations are described in the Appendix A.

**Figure 2.** Schematic hydrological phase of the DMMF model (modified from Figure 3 of Choi et al. [15]).

**Figure 3.** Schematic sediment phase of the DMMF model (modified from Figure 4 of Choi et al. [15]).

In contrast with the SWAT model, which has been frequently applied to this catchment, the DMMF model can estimate the erosion and deposition of an element, considering the interconnectivity with adjacent elements. Therefore, the model can be used to estimate the impact of the spatial reconfiguration of erosion hot spots into forest on sediment yields more explicitly for each element and the entire catchment.


**Table 1.** Input parameters of the daily based Morgan–Morgan–Finney (DMMF) model (modified from Table 1 of Choi et al. [15])

#### *2.3. Model Parameterization*

As shown in Table 1, the DMMF model requires the topography, climate, soil, and LULC datasets to project surface runoff and sediment redistribution patterns of the catchment.

Topography data (i.e., the slope angle (*S*) and grid size of a raster map (*res*)) were derived from the digital elevation model (DEM) with 30 m resolution. The parameter *res* is used to calculate the width (*w*) and length (*l*) of an element that are equivalent to *res* and *res*/ cos(*S*), respectively [15].

Climate data were obtained from two sources. The daily rainfall (*R*) and mean rainfall intensity of a day (*RI*) were obtained from weather stations installed in the catchment, and the evapotranspiration (*ET*) was obtained from remote sensing data provided by the Moderate Resolution Imaging Spectroradiometer (MODIS) [43]. We estimated *R* and *RI* from each weather station and spatially interpolated them using inverse distance weighted (IDW) method, which showed the optimal result on this catchment among four methods such as inverse distance weighted, spline, nearest neighbor, and kriging, according to Shope et al. [30]. For the *ET*, we resampled the 8-day average MODIS/Terra Evapotranspiration data to fit to the DEM of this catchment.

The soil data set covers the texture, depth, hydraulic properties, and detachabilities. The soil texture (i.e., the proportion of clay (*Pc*), silt (*Pz*), and sand (*Ps*) in the surface soil), soil depth (*SD*), and soil hydraulic properties (i.e., saturated soil water content (*θsat*), soil water content at field capacity (*θ f c*), and saturated lateral hydraulic conductivity (*K*) of the entire soil profile) were derived from a 2009 catchment-wide field survey from the TERRECO project (see Table 2 and Figure 4) [30].


**Table 2.** Typical soil characteristics of each represented soil class of the Haean catchment from a 2009 catchment-wide field survey from TERRECO project.

\* *θsat*, *θ f c*, and *K* were estimated with the model ROSETTA Lite v.1.1 [44]. The numbers in parentheses indicate the range of values of soil layers that constitute each represented soil class.

**Figure 4.** Represented soil class from a 2009 catchment-wide field survey from the TERRECO project.

Reference values for soil detachability from Morgan and Duzant [40] were used as the initial values of soil detachability by rainfall (i.e., for clay (*DKc*), silt (*DKz*), and sand (*DKs*)) and by runoff (i.e., for clay (*DRc*), silt (*DRz*), and sand (*DRs*)). We assumed that the initial soil water content of the entire soil profile (*θinit*) is equal to the soil water content at field capacity (*θ f c*) by starting the simulation at three days after the first heavy rainfall of the year, because the excess soil water was usually drained away two or three days after the soil was fully saturated by rainfall.

The LULC types characterize the physical structures of surface and vegetation, which regulate the quantity of surface runoff and runoff velocity. Surface structures incorporate a portion of the impervious cover area (*IMP*), such as plastic mulching and paved facilities, flow depth of surface runoff (*da*), and Manning's roughness coefficient of the soil surface (*n*). Vegetation structures contain the permanent interception of rainfall (*PI*), pervious ground cover (*GC*), canopy cover (*CC*), average vegetation height (*PH*), average diameter of individual plant elements at the surface (*D*), and number of individual plant elements per unit area (*NV*). LULC parameters were derived based on the LULC map of the Haean catchment in the year 2010 from Seo et al. [35] (see Figure 5).

**Figure 5.** LULC classes and their spatial configurations for the Haean catchment in the year 2010 [35].

We classified the original LULCs into 14 categories (i.e., forest, rice paddy, semi-natural, bare soil, ginseng, potato, bean, radish, cabbage, other dry crops, shrub, orchard, urban, and water bodies). Forest, rice paddy, semi-natural, bare soil, ginseng, potato, bean radish, and cabbage are major LULCs that covered more than 1% of the catchment area. Minor LULCs were aggregated into groups of other dry crops, shrub, orchard, urban, and water bodies according to their physical characteristics. We used field measurement data of *CC*, *PH*, *NV*, *IMP*, *da*, and *n* for major dry crops such as bean, cabbage, potato, and radish, whose data were obtained from the field campaign of the TERRECO project, which was also used in Arnhold et al. [17]. The daily forest *CC* was estimated using the average values of 8-day normalized difference vegetation index (NDVI) for forest in the catchment from MODIS [45,46]. The average NDVI values were converted to canopy cover (*CC*), using the equation suggested by Gutman and Ignatov [47]. LULC parameters for rice and ginseng, and the average diameter of individual plant elements (*D*) for major dry crops were obtained from agricultural technology portal provided by Rural Development Administration of South Korea (RDA) [48]. The average LULC parameters of major dry crops were used for the LULC parameters of other dry crops, while the guide values from Morgan and Duzant [40] were adopted for other LULC parameters. Detailed initial parameter settings are presented in Table 3.


**Table 3.** The initial parameter settings for each LULC class.

(*a*) Typical leaf-out and leaf-fall dates of each LULC were presented as day of the year (DOY). For annual crops, the dates represented the typical planting and harvest date of each crop [30]; (*b*) The reference values from Morgan and Duzant [40] were used for the area proportion of the permanent interception of rainfall (*PI*) for each LULC type; (*c*) *IMP* for dry fields are different between cultivation and non-cultivation periods. Values in parentheses represent *IMP* for non-cultivation periods; (*d*) *GC* for dry fields is different before and after harvest. After harvest, crop residues and weeds remained as the ground cover of dry fields, according to dry crop data, from the field campaign of the TERRECO project in 2009. *GC* for rice paddy in cultivation season was set to one reflecting water-filled condition that protected the surface from erosion; (*e*) Because *CC* values varied with time, we made a list of maximum *CC* (*CCmax*). Semi-natural, shrub, and ginseng utilize fixed reference values from Morgan and Duzant [40]; (*f*) We used fixed reference *PH* values from Morgan and Duzant [40] for LULCs of other than dry crops. Maximum *PH* values for dry crops were listed from the field measurement data varying with time; (*g*) We used fixed reference *D* values from Morgan and Duzant [40] for LULCs of other than dry crops. *D* values for dry crops utilized typical crop characteristics from Rural Development Administration of South Korea [48]; (*h*) We used reference *NV* values from Morgan and Duzant [40] for LULCs of other than dry crops and ginseng. *NV* values which were estimated from the field measurement data and Rural Development Administration of South Korea [48] were used for dry crops and ginseng, respectively; (*i*) We assumed shallow rill condition for forest, semi-natural and shrub, and assumed unchannelled flow condition for bare soil, rice paddy, and orchard using values presented in Morgan and Duzant [40]. *da* values for other LULCs derived from furrow heights of the fields, using field measurement data for dry crops and data from the Rural Development Administration of South Korea [48] for ginseng; (*j*) According to the guide values for Manning's *n* from Morgan [49], the values of *n* for natural land covers (i.e., forest, semi-natural, and shrub), crop fields, ginseng, and smooth surfaces (bare soil and urban) are 0.2, 0.1, 0.2, and 0.01, referring to natural range land, average tillage conditions, wheat mulching, and smooth bare soil or asphalt conditions, respectively; \* The permeable black awning screen is generally installed 1.3 m above the ginseng field [48], and it acts as a plant canopy. Therefore, the cover ratio of the screen in the field and height of the screen is utilized for canopy cover (*CC*) and plant height (*PH*) values for ginseng.

#### *2.4. Model Calibration and Validation*

The DMMF model was calibrated and validated for stream discharge and suspended sediment to test its performance in the Haean catchment. The testing was performed utilizing data from the year 2010 when the LULC map, as well as the field-measured stream discharge and suspended sediment data, were well established [30,35]. We confined the testing period from the 67th day of the year (DOY), which is three days after first heavy rainfall of the year, to reduce the uncertainty of initial soil water content by equating it with the soil water content at field capacity. We equalized the two parameters based on the field measurement guidelines for soil water content at field capacity, which recommend soil sampling two or three days after rainfall that is heavy enough to saturate the soil. The three sub-catchments of S1, S2, and S3 (see Figure 1) were selected for model calibration and validation. The data from the S1 and S2 were utilized for two-step calibration, and those from the S3 were used for model validation. Two-step calibration was performed on the forest-related parameters utilizing the data from the S1 site, and the other parameters were calibrated utilizing the data from S2. This calibration method enables us to prevent the significance of forest-related parameters of dominant LULC type in the entire catchment, from overtaking the importance of other parameters, resulting in those parameters being ignored. The DMMF model can estimate the outputs of the surface and subsurface runoff, and the sediment from the elements. However, the measured data are stream discharge and suspended sediments at the outlet of each sub-catchment. Because the model does not consider in-stream processes and the impact of groundwater on the base flow of the stream, it is not

appropriate to directly compare the result from the model with the measured data. To match different comparative objects, we compared the total daily discharge of each site to total daily surface runoff and subsurface interflow flowing into the stream from the model, while adding a constant corresponding to base flow from groundwater. To match the sediment yield from the terrestrial part with the suspended sediments measured at the outlet of each sub-catchment, we should consider the in-stream sediment processes and impact of erosion control facilities. Reflecting sediment deposition on the stream bed load, we assumed that only a part of the terrestrial sediment yield entering the stream was sampled at each measuring point for each sub-catchment. Therefore, we compared the suspended sediments measured from the outlet of each measuring point to the sediment flowing into the stream from the model, multiplied by a constant, reflecting the in-stream sediment process. Our assumptions can be described as below,

$$Q\_m = Q\_s + IF\_s + a\_r \tag{1}$$

$$SL\_m = \beta \times SL\_s.\tag{2}$$

Here, *Qm* represents the measured daily total discharge, and *Qs*, *IFs*, and *α* represent the daily surface runoff, daily subsurface interflow simulated from the DMMF model, and a constant reflecting the base flow from groundwater (unit: m3/s). *SLm* represents the total daily suspended sediments measured at the outlet of each sub-catchment, and *SLs* and *β* represent the terrestrial sediment yield entering the stream from the model simulation and constant representing the in-stream sediment deposition rate, respectively.

#### 2.4.1. Sensitivity Analysis

To select important parameters to be calibrated among unmeasured or highly uncertain parameters, we performed site-specific sensitivity analyses, using the Sobol' method [50–52]. The Sobol' method is a variance based sensitivity analysis that is widely used in environmental and hydrological modeling, such as SWAT and TOPMODEL [53,54]. This method can estimate the total effect of each parameter on the model output, considering the combined effects among parameters. Therefore, the Sobol' method is more suitable for analyzing the sensitivity of non-linear and non-additive models containing many parameters, as opposed to the local or one-at-a-time (OAT) methods [53,55]. The relative sensitivity of parameters is expressed as the Sobol' total index (*SI*)—the ratio of the amount of total variance caused by a parameter to the amount of variance induced from all parameters (i.e., the unconditional variance of the model) [52]. If we have p-dimensional parameter set, the first-order sensitivity of the *i*-th parameter can be described as,

$$S\_i = \frac{V\_{X\_i}(E\_{X\_{-i}}(Y|X\_i))}{V(Y)},\tag{3}$$

where *VX*−*<sup>i</sup>* (*EXi* (*Y*|*X*−*i*)) is the variance of the model solely by *i*-th parameter (*Xi*). Then the total sensitivity of the *i*-th parameter (*SIi*) can be calculated as below,

$$SI\_i = 1 - \frac{V\_{X\_{-i}}(E\_{X\_i}(Y|X\_{-i}))}{V(Y)},\tag{4}$$

where *VX*−*<sup>i</sup>* (*EXi* (*Y*|*X*−*i*)) *<sup>V</sup>*(*Y*) indicates that the sum of first-order sensitivities of all parameters except *<sup>i</sup>*-th parameter. Parameters with large *SI* indicate a relatively high impact on the model output, while those with small *SI* indicate a relatively low impact on the model output.

Because the soil hydraulic parameters (i.e., *θsat*, *θ f c*, and *K*), soil detachabilities (i.e., *DKc*, *DKz*, *DKs*, *DRc*, *DRz*, and *DRs*) and LULC parameters (i.e., *PI*, *IMP*, *GC*, *CC*, *PH*, *D*, *NV*, *d*, and *n*) were not measured or had high uncertainties, their importance was tested on model outputs. Before performing

sensitivity analysis, we set the range of the parameters to be tested. The ranges of soil hydraulic parameters (i.e., *θsat*, *θ f c*, and *K*) were set based on the range of estimated values for each represented soil class (see Table 2). The upper bound of *θ f c* was set as the minimum *θsat*, and the upper bound of *K* was set to 18 times of the maximum *K* to reflect high uncertainty of the parameter [56]. The ranges of the un-measured LULC parameters were set based on the initial parameter settings for each LULC type (see Table 3). We adjusted the parameters using a range of ±100% for the initial parameter settings for each LULC type. If the upper or lower limits of the proportional parameters is out of the range between zero to one, we set the lower limits to zero and the upper limits to one. In this study, SIs for the input parameters were estimated using the "sobolmartinez" function of the "sensitivity" package [57] on R version 3.5.1 [58], a well-established open-source program for statistical computing, providing many analysis packages. We used the default bootstrapping option of the function, employing a sample size of 103.

#### 2.4.2. Calibration

To find the optimal combination of the parameter set, which allows model outputs to explain the measured stream discharge and suspended sediments from each site, we performed two-step calibration. For each step, we adjusted the important parameters with *SI* greater than 0.05 (i.e., contributing 5% of the total variance), and we adjusted the constants for the in-stream processes (*α* and *β*) additionally for sub-catchment S2, where data were measured in the stream outlet. We searched for the optimal combination of the parameter set, using the differential evolution (DE) optimization method [59,60]. The DE algorithm is a heuristic optimization method with an evolution strategy for finding the global optimum value. Requiring few prerequisites for its execution, the algorithm is applicable to non-differential, nonlinear, and multimodal models. As a result, the DE algorithm has been applied to a variety of fields including hydrological model calibration [15,59–63]. We applied the DE algorithm for model calibration using the "DEoptim" package [61,64] on R version 3.5.1 [58]. We used the Nash-Sutcliffe efficiency coefficient (NSE) [65] between model outputs and field-measured data as an objective function for the DE algorithm. To treat NSE values from stream discharge and suspended sediments fairly, we evaluated the NSE values for each measurement and used the average NSE value as the final objective function:

$$F\_{obj} = 1 - \frac{\text{NSE}(Q\_m) + \text{NSE}(SL\_m)}{2},\tag{5}$$

where *Fobj* is the objective function to evaluate the model performance. We ran the function for 10<sup>3</sup> iterations, and ran for three different initial states to try to find the global minimum as an optimum value.

#### 2.4.3. Validation

Using adjusted parameters from calibration steps, model performance was tested for the S3 site, which is located near the catchment outlet. Considering site-specific base flow from groundwater and in-stream sediment processes for the S3 site, we adjusted the constants for the in-stream processes (*α* and *β*). We utilized the NSE, the percent bias (PBIAS), and the coefficient of determination (R2) as statistical criteria for model performance evaluation [66,67]. The function "gof" from the "hydroGOF" package [68] in R version 3.5.1 [58] was used to evaluate statistical criteria.

#### *2.5. Identifying Annual Sediment Redistribution Patterns and Assessing Soil Erosion Risk*

Projecting validated parameters on the DMMF model, we simulated and calculated the annual sediment redistribution patterns of the catchment. Based on the simulated result, we assessed the net soil erosion rate (*SLnet*: t/ha/year) for each element of the catchment. *SLnet* is the net soil erosion for each element, which is the amount of sediment input to each element from upslope elements (*SLin*) subtracted from the amount of sediment output from the element (*SLout*). Soil erosion risk was assessed by using *SLnet* of each element. We classified *SLnet* into five categories, namely tolerable, low, moderate, high, and severe, as shown in Table 4 according to the soil erosion risk categories defined by OECD [69,70] which is one of the internationally used criteria.

**Table 4.** Soil erosion risk categories defined by OECD [69,70].


Based on the net soil erosion rate of the entire catchment, we assessed the soil erosion characteristics for each LULC class. For the assessment, we calculated the mean *SLnet* for each LULC class.

#### *2.6. Evaluation of the Impact of Spatial Reconfiguration of Erosion Hot Spots into Forest*

We assessed the impact of the spatial reconfiguration of erosion hot spots into forest, based on the annual sediment redistribution patterns of the catchment. Erosion hot spots represent elements in which much annual net soil erosion (*SLnet*) occurs. To compare the impact of spatial reconfiguration, we calculated the annual sediment yields being generated from the terrestrial area and entering to the water bodies of the entire catchment (*SYbase*) as a base line condition. *SYbase* is the total amount of sediment yields entering the water bodies of the entire catchment, which is equal to the total amount of *SLin* flowing into water bodies. To increase the robustness of our analysis, we only used the values between the 2.5th percentile and the 97.5th percentile for all the elements in the catchment to exclude the impact of extreme values that can occur from model outputs. The lower extreme values were set to the value of the 2.5th percentile and the upper extreme values were set to the value of the 97.5th percentile. The impact of the spatial reconfiguration of erosion hot spots into forest was evaluated by calculating the total annual sediment yields entering the stream (*SYtot*), using the DMMF model as bare soil and croplands (i.e., bean, cabbage, ginseng, orchard, potato, radish and rice field) being sequentially changed into the forest. We selected forest, the original LULC type before anthropogenic land cover changes, as the alternative LULC to mitigate erosion-prone areas. Similar to the methods Chaplin-Kramer et al. [23] and Chaplin-Kramer et al. [24] which compute ecosystem services by marginally changing forest into agricultural areas, we computed *SYtot* by gradually converting 1% of the bare soil and croplands in the catchment into forest until all bare soil and croplands elements are converted into forest. Based on this result, we presented the total sediment yields (*SYtot*), reduction rate of the sediment yields entering the stream compared to base line condition (*SYbase*), and sediment yield reduction efficiency per conversion area (t/m2).

#### **3. Results**

#### *3.1. Model Performance*

According to the calibration and validation results, the DMMF model showed competitive performance, predicting stream discharge, but showed poorer performance in evaluating the amount of suspended sediments at the outlet of each sub-catchment. We performed two-step calibration by comparing the model outputs to the measured data collected from sub-catchment S1 and S2. The LULC and soil types of sub-catchment S1 are classified as forest and forest soil, according to Tables 2 and 3. The calculated Sobol' index for important parameters, both for stream discharge (*SIQ*) and suspended sediments to the stream (*SISL*), are presented in Table 5.

**Table 5.** List of important parameters from forested site with Sobol' index greater than 0.05 for stream discharge (*SIQ*) and suspended sediment to the stream (*SISL*), and their optimized values from the DE algorithm.


According to the Sobol' index, the amount of stream discharge was highly influenced by the permanent interception of rainfall (*PI*) and lateral soil hydraulic conductivity (*K*), which regulate the amount of rainfall and flow rate of subsurface interflow of the sub-catchment, respectively. Vegetation and surface cover structures (*GC*, *PI*, and *da*), detachability of clay particles (*DRc*), soil water content at field capacity (*θ f c*), and lateral soil hydraulic conductivity (*K*) exhibited a relatively large impact on suspended sediments generated from the sub-catchment. This result indicates that the suspended sediments generated from the sub-catchment are determined by the amount of surface runoff and the erosivity of surface, because *PI*, *K*, and *θ f c* determine the amount of surface runoff by regulating the amount of rainfall and partitioning the rate of surface and subsurface water. Parameters *GC*, *da*, and *DRc* determine the erosivity by surface runoff.

We determined an optimized parameter set by adjusting selected important parameters from sensitivity analysis using the DE algorithm (see Table 5). With the optimized parameter set, the stream discharge and suspended sediment from the model outputs were compared with those from field measurements (see Figure 6).

After calibrating the forest-related parameters, we calibrated the other parameters, based on the measurement data collected from sub-catchment S2. We calculated the relative importance of parameters for both the stream discharge (*SIQ*) and suspended sediments to the stream (*SISL*), using the Sobol' index, and presented them in Table 6.


**Table 6.** List of important parameters (*SI* > 0.05) for stream discharge (*SIQ*) and suspended sediment (*SISL*), and their optimized values from DE algorithm.

**Figure 6.** Calibration (S1 and S2) and validation (S3) result for stream discharge, and suspended sediment measured at the outlet of each sub-catchment.

According to sensitivity analysis, model outputs were highly sensitive to soil hydraulic characteristics of moderate to steep dry field soil and land cover structures of the semi-natural field. In details, the stream discharge of the sub-catchment was highly sensitive to the permanent interception of rainfall (*PI*) of the semi-natural, rice paddy, and other dry crops; the lateral hydraulic conductivity (*K*) of the moderate to steep dry field and flat dry field soils; and the soil water content at field capacity (*θ f c*) of the moderate to steep dry field. This result indicates that stream discharge is highly influenced by the amount of rainfall reaching the ground (*PI*s) and the flow rate of subsurface interflow (*K*s and *θ <sup>f</sup> c*) of this region. The sediment yield to the stream is sensitive to the soil detachability by runoff (*DRc* and *DRz*) of the moderate to steep dry field soil, soil water content at field capacity (*θ f c*) of the moderate to steep dry field soil, flow depth (*da*) of the semi-natural field and bean field, and ground cover ratio (*GC*) of the semi-natural field. This result emphasizes the role of the moderate to steep dry field soil, which is the second largest soil type, following forest soil, and demonstrates the crucial role of the semi-natural field on determining suspended sediment output from the model.

The performance statistics for the calibration and its time series plots of observed versus simulated stream discharge and suspended sediment were presented in Figure 6. For the calibration steps, the NSE values for stream discharge were 0.92 and 0.88 for sub-catchment S1 and S2, respectively. The R2 values for stream discharge were 0.93 and 0.88, respectively, and the PBIAS values for stream discharge were −18.6 and 0.1, respectively. The NSE values for suspended sediment were 0.99 and 0.43 for sub-catchments S1 and S2, respectively. The R2 values for suspended sediment were 0.99 and 0.44, and the PBIAS values for suspended sediment were −6.8 and −22.1 for the sub-catchments, respectively. The site-specific constants reflecting the baseflow from groundwater (*α*) and in-stream sediment deposition rate (*β*) for sub-catchment S2 are 1.75 × <sup>10</sup>−<sup>2</sup> <sup>m</sup>3/s and 4.57 × <sup>10</sup><sup>−</sup>2. In validation steps, the NSE values for stream discharge and suspended sediment were 0.75 and 0.18, respectively, with the site-specific *<sup>α</sup>* and *<sup>β</sup>* being 1.711 m3/s and 6.76 × <sup>10</sup><sup>−</sup>2, respectively. The R2 for discharge and sediment were 0.83 and 0.39, respectively, and the PBIAS for discharge and sediment were 0 and −40.5, respectively. According to the model performance evaluation criteria suggested by Moriasi et al. [67], the DMMF model showed good performance for discharge in both calibration and validation steps. Though there is no clear model performance evaluation criteria suggested for daily time scale sediment result for watershed model due to limited reported data [67], When we apply the performance evaluation criteria for monthly time scale sediment result for watershed scale model, the model might be considered to have a slightly poor performance for sediment during the calibration and validation steps, as the NSE and R<sup>2</sup> values were less than 0.45 and 0.40, respectively.

#### *3.2. Sediment Redistribution Pattern of the Catchment*

Simulating the model with optimized parameters, we calculated the annual net soil erosion rate (*SLnet*) for each element and classified them into five classes–tolerate, low, moderate, high, and severe—as in Figure 7.

According to Figure 7, elements with severe soil erosion (>33 t/ha/year) were concentrated on the dry crop field with moderate to steep slope conditions on the interface with the forest. The estimate of the mean annual net soil erosion rate by each LULC type (Table 7) shows that bare soil and dry crop field suffered from severe soil erosion. On the other hand, forest, rice paddy, orchard, and urban areas showed good sediment capturing capabilities.

**Figure 7.** (**a**) Annual net soil erosion (t/ha/year) of the entire Haean catchment and (**b**) soil erosion class according to the soil erosion risk categories from OECD [69,70].


**Table 7.** Mean annual net soil erosion rate (t/ha/year) and mean slope of each LULC type.

*3.3. Impacts of Conversion of Erosion Hot Spots into Forest on Total Sediment Yield Entering the Stream*

The LULC conversion of erosion hot spots into forest showed a dramatic impact in the reduction of sediment yields entering the stream, as shown in Figure 8.

**Figure 8.** Total annual sediment yields entering the stream (**upper panel**) and sediment yield reduction efficiency per unit conversion area (**lower panel**) through changing bare soil and crop fields into forest sequentially from the area with the highest to the area with the lowest amount of net soil erosion (*SLnet*), sediment inflow to the element (*SLin*), and sediment output from the element (*SLout*).

When each bare soil and crop field element in the catchment was converted into the forest sequentially from the area with the highest soil erosion rate to the area with the lowest soil erosion rate, the amount of total annual sediment yield of the catchment to the stream sharply decreased having a shape similar to an inverted sigmoid function. Changing the 3% of erosion hot spots that have suffered the most from severe soil erosion caused a reduction in sediment yield entering the stream of ca. 10% from the baseline condition (*SYbase*), and a change in 10% of most severe hot spots is expected to reduce sediment yields by ca. 50%. Among the elements *SLnet*, *SLin*, and *SLout*, the altered areas revealed that outputs from the element (*SLout*) proved to be the most effective in reducing the total sediment yield into the stream. A simulation of the sediment yields entering the stream showed that the reducing rate in sediment yield for *SLnet* was less effective than those for *SLout* and *SLin*. Due to total annual sediment yields sigmoidally decreases as bare soil and crop fields begin changed into forest, sediment yield reduction efficiency per unit conversion area increased until ca. 10% of total crop land area converted to forest and then gradually decreased. A simulation of the sediment yield reduction efficiency showed that the element (*SLout*) was most efficient for all conversion intervals.

#### **4. Discussion**

Our findings emphasize the importance of landscape configuration on regulating ecosystem services by showing the effectiveness of spatial reconfiguration of soil erosion hot spots into forest on reducing the amount of sediment yield entering the stream. We simulated the annual sediment redistribution pattern in the Haean catchment, utilizing the daily based Morgan–Morgan–Finney (DMMF) soil erosion model. According to the result, the soil erosion rate varied greatly depending on the topography and LULC type, and the area located on the steep hill-slope, which is adjacent to the forest severely suffered from soil erosion. When reconfiguring the landscape patterns of croplands by sequentially altering erosion hot spots from the most severe to the least severe areas into forest, we found dramatic effects in the reduction of sediment yields entering the stream in this catchment. The reduction rate may reach ca. 50% when the 10% most severe erosion hot spots were altered, and we can expect a reduction rate of over 80% when the ca. 20% most severe erosion hot spots are altered. In the following, we first discuss model performance and limitation, and then potential management implications.

#### *4.1. Model Performance*

The assessment of soil erosion risk and measurement of the effectiveness of the spatial reconfiguration of erosion hot spots in reducing sediment yields entering the stream were based on the calibrated and validated simulations of the DMMF soil erosion model. According to the model performance criteria from Moriasi et al. [67], the DMMF model showed satisfactory performance for predicting stream discharge during the calibration and validation processes, with mean NSE values of 0.90 and 0.75, mean R2 of 0.91 and 0.83, and maximum PBIAS of −18.6 and 0 during calibration and validation steps, respectively. The model showed comparatively poor performance for predicting suspended sediment at the outlet of each sub-catchment, except the small forested site (S1) where the stream does not exist. The mean NSE values were 0.66 and 0.18, mean R<sup>2</sup> were 0.67 and 0.39, and maximum PBIAS were −22.1 and −40.5, respectively. When we compared the model performance statistics of the DMMF model to those from previous studies using soil and water analysis tool (SWAT), the model showed competitive performance in predicting stream discharge but poorer performance in terms of predicting suspended sediments in the stream [3,19]. Maharjan et al. [3] reported that mean NSE values for stream discharge were 0.82 during calibration and 0.45 during validation. In addition, they showed that mean NSE values for suspended sediment were 0.78 and 0.60 during calibration and validation, respectively. Jang et al. [19] also reported mean NSE values for stream discharge of 0.78 and 0.66 during calibration and validation, respectively. They reported mean R<sup>2</sup> for suspended sediment were 0.80 and 0.76 during calibration and validation, respectively. In terms of soil erosion rate for each crop field, the DMMF model estimated that the average annual soil loss of major dry crops

ranged between 79.3 t/ha/year and 763.8 t/ha/year for bean, radish, potato, and cabbage, and the average annual soil loss from whole dry crop fields was 379.7 t/ha/year. Arnhold et al. [17] reported that 30–54 t/ha/year of soil loss occurred in the dry crop fields, including bean, radish, potato, and cabbage, from the plot-level field measurement. Furthermore, Maharjan et al. [3] estimated that 35.5–53.0 t/ha/year of soil loss occurred in the dry crop fields from the SWAT model. When we compared the results from the DMMF model with those from other studies, the amount of soil loss from this study is far greater. The reasons that the DMMF model showed poor performance for predicting suspended sediment in the stream can be analyzed from two perspectives. The first reason involves the discrepancy of data types between the DMMF model and observed data. The observed data were stream discharge and suspended sediment at the outlet of each sub-catchment. On the other hand, the DMMF model can estimate the total sediment yields entering the stream that belongs to each sub-catchment. The DMMF model is efficient for estimating sheet and rill erosion, but it has limitations in estimating in-stream sediment processes such as stream bed deposition, channel erosion, and sediment transport in the stream. Considering the limitations of the model, we use site-specific coefficients, which assume that suspended sediments measured at the outlet are proportional to the sediment yields inflowing into the stream. However, incorporating the quantity and the velocity of stream water discharge, sediment flux, and physical characteristics of channel structures such as gradient, width, depth, and length, into the in-stream sediment process, is complicated [71,72]. The reasons above may lead to a high sediment deposition rate in the stream (i.e., low measured sediment ratio (*β*)), which in turn, causes a high soil erosion rate in the terrestrial area. Because the study sites are affected by monsoon climate, such that its rainfall pattern is not uniform but rather with a lot of extremes, a large amount of sediment is deposited during low rainfall events, and the deposited sediments are washed out by a huge amount of fast stream discharge accompanying heavy rainfall. Temporal lags between the rainfall event and stream discharge are negligible for the Haean catchment, but for suspended sediments, the lags are significant and highly depend on the stream length because of the difference in travel velocities between water and soil particles [73–76]. Therefore, the model performance for predicting stream discharge may be better than that for predicting suspended sediments. The stream widens and deepens as it descends to the lower area, according to Lee [2], and the length of the stream also increases as the size of sub-catchment grows. The uncertainty caused by in-stream processes increases as the size of the sub-catchment grows, which reduces the model performance in predicting suspended sediments in this study. SWAT and USLE-based models are usually calibrated and validated at the fixed spatial area with a different temporal period. Therefore, in-stream sediment processes can be included in the parameters, which may lead to better model performance. However, the DMMF model is a spatially distributed semi-processed model and used the same temporal period with a different spatial area for calibration and validation in this research, so that the in-stream processes cannot be included in the model.

Secondly, many sediment reduction facilities, such as dams for freshwater, debris barrier and culvert systems around crop fields, and road infrastructures, which can affect sediment transport processes, have been installed in the Haean catchment [30,36]. The dam and debris barriers create reservoirs that impede the stream flow and filter out sediments in the facilities. This disrupts the correct evaluation of the model performance for this catchment. Shope et al. [30] showed complex stream networks, including the culvert systems around crop fields and the road infrastructure. The culvert systems extend the travel time of suspended sediments and reduce the runoff and transport velocities of sediments by altering the flow direction abruptly. Increased travel time and decreased transport velocity tend to increase the deposition rate of sediments compared to the condition without the culvert system. The deposited sediments in the culvert flow into the stream by runoff, with sufficient power to wash out. The culvert system is also responsible for the temporal lag between the rainfall event and the presence of suspended sediments in the catchment. Sediment reduction facilities trap a huge amount of sediments, which make the measured sediment ratio (*β*) in this study have very low values. Because of the small *β*, the stream bed deposition rate became too large, and consequently, the overall

erosion rate from terrestrial area increased. To cope with this problem, the in-stream processes will need to be considered more precisely through model improvements.

#### *4.2. Assessment of Soil Erosion Risk and the Effectiveness of Spatial Reconfiguration of Erosion Hot Spots on Reducing Sediment Yield Entering the Stream*

We estimated the annual net soil erosion rate of the entire catchment and assessed the soil erosion risk class according to the OECD criteria. According to this study, soil erosion is concentrated on the hill-slope of the catchment, and the problem is more significant for the bare soil and dry crop fields, such as bean, radish, and potato, in this area. In addition, forest in the valley showed a considerable amount of soil loss, also suffering from erosion due to the concentrated surface runoff and steep slope. Compared with other studies, the soil erosion risk pattern and the average annual soil loss from the DMMF model is qualitatively consistent with the soil erosion risk map from Lee et al. [77], with average climate conditions for the 2010s using the USLE-based SATEEC [78] model. According to this study, urban area, orchard, and rice field showed better performance for sediment capturing capabilities than forest. However, the urban area and rice field are located in the lower and flatter area than forest, so that the sediment inputs from the upslope area tend to be deposited in this area. Furthermore, because the urban area is usually paved with impervious covers, such as concrete and asphalt, and the rice field is filled with water, which acts as a pervious cover that prevents surface erosion, these areas have little soil loss but receive huge input from the upslope area. Though the forest is in a region where the slope is very steep, the average amount of soil loss is smaller compared with other land types, and it also shows excellent sediment capturing capability, in general. Like the other studies, we can conclude that the main cause of severe erosion in the catchment is cropland extension after deforestation at the hill-slope area of the catchment [2,3,11,17,19,77]. We also assessed the effect of spatial reconfiguration of LULCs on reducing sediment yields entering the stream. In this study, the spatial reconfiguration of erosion hot spots into forest showed excellent reduction efficiency in sediment yields entering the stream. We identified that the sediment yields entering the stream were reduced sharply, as crop lands were sequentially changed into forest from the area with the most severe soil loss to the area with the least soil loss. An sigmoidal sediment reduction rate from altering LULCs to forest indicates that forest is not only effective in preventing surface erosion but also effective in capturing sediment input from the upslope area. In addition, the result suggested that altering LULCs based on the amount of sediment output from the element is the most effective way of reducing sediment yields entering the stream. This result is consistent with previous studies that emphasize the effectiveness of vegetative filter strips located at sediment sources such as crop fields [3,5,19,27–29]. The result can also be generalized to consider the effect of riparian vegetation buffer strip on reducing sediment yields entering the stream, located at the interface between crop fields or natural sediment sources and the stream channel [4,79,80]. This study also demonstrated that the sediment yield reduction efficiency initially increased as the first few bare soil area and crop lands with the most severe soil loss were converted into forest. The sediment yield reduction efficiency were maximized when ca. 10% of the area converted, and then the efficiency decreased gradually. These patterns can be explained by two aspects of the forest's sediment yield reduction capability; protecting surface from soil erosion, and capturing sediment inputs. The areas with the most severe soil loss are located at the steep hillslope where surface runoff is concentrated. These areas have a large transport capacity of the runoff, beyond the sediment capturing capability of forest because transport capacity is greater than the available sediment for transport [15]. In these areas, conversion of crop lands into forest can reduce soil loss from the surface but cannot capture sediment inputs from upslope which is larger than surface soil loss. As slope becomes milder and the amount of surface runoff decreases due to gradual conversion of crop lands into forest, transport capacity gradually decreases. Decreased transport capacity caused by decreased slope gradient and surface runoff lets forest capture more sediments, maintaining the surface protecting capability from soil loss. Therefore, the sediment yield reduction capabilities of forest become small and the sediment yield reduction efficiency by changing crop lands into forest

decreases gradually. According to these results, one can reduce sediment yields entering stream efficiently by identifying an optimal percentage of crop land conversion into forest which brings out the best efficiency of sediment yield reduction per unit conversion area.

#### **5. Conclusions**

In this study, we identified the soil erosion risk of Haean catchment spatially explicitly by projecting sediment redistribution patterns using the DMMF model. In addition, we measured the sediment yield reduction efficiency entering the stream by sequentially altering erosion hot spots into forest from that which has the highest soil loss to that which has the lowest soil loss. The DMMF model showed competitive performance estimating stream discharge but exhibited lower performance estimating suspended sediments at each sub-catchment outlet. When we applied the DMMF model to the Haean catchment, the bare soil surface and dry crop fields located on the steep hill-slope of the catchment suffered mostly from severe soil erosion. On the other hand, forest, rice paddy, orchard, and urban areas suffer less from soil erosion. By changing the erosion hot spots from cropland to forest, the overall amount of sediments exporting to the stream of the catchment was effectively reduced. The sediment yield reduction efficiency was maximized when ca. 10% of crop lands were converted to forest. This study implies that one can achieve the goal of reducing sediment yields entering the stream by identifying the location of erosion hot spots and managing the area intensively. Although previous studies showed good mitigation effects of BMPs that require compliance of stakeholders, this may not be easy and takes much time for stakeholders to follow the BMPs, because the degree of acceptance of the policy depends on the situation and tendency of each stakeholder [19]. On the other hand, the spatial reconfiguration approach proposed in this study can reduce the number of stakeholders relevant to soil erosion mitigation measures. However, this approach reduces crop yields because crop lands are converted to non-crop lands to reduce sediment yields from the catchment. In addition, the sediment yield reduction efficiency decreases after a certain point of spatial reconfiguration. Therefore, the two approaches—BMP measures such as cultivating cover crops, mulching surface with straw, and managing field margin naturally, and conversion of crop lands with the more severe soil loss—are complementary measures to reduce sediment yields into the stream.

**Author Contributions:** Conceptualization, K.C. and B.R.; Data curation, K.C. and G.R.M.; Formal analysis, K.C.; Investigation, K.C.; Methodology, K.C.; Resources, G.R.M.; Software, K.C.; Supervision, B.R.; Validation, K.C.; Visualization, K.C.; Writing original draft, K.C.; Writing review & editing, G.R.M. and B.R.

**Funding:** This study is part of the International Research Training Group "Complex Terrain and Ecological Heterogeneity" (TERRECO) funded by the German Research Foundation (DFG) with the grant number [GRK 1565/1]. This study was also supported by the National Research Foundation of Korea (NRF) with the grant number [NRF-2017R1A2B4010460].

**Acknowledgments:** This paper is dedicated to the memory of our wonderful colleague, Sebastian Arnhold, who passed away last year. This publication was funded by the German Research Foundation (DFG), the National Research Foundation of Korea (NRF). This publication was also funded by the German Research Foundation (DFG) and the University of Bayreuth in the funding programme Open Access Publishing. We would like to thank Editage (www.editage.co.kr) for English language editing.

**Conflicts of Interest:** The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

#### **Appendix A. Detailed Structure of the Daily Based Morgan–Morgan–Finney (DMMF) Soil Erosion Model**

Morgan–Morgan–Finney (MMF) model [37] is a conceptual soil erosion model, which estimates the annual soil erosion rate from an area by comparing the amount soil particles detached from the surface (*SS*) and transport capacity of surface runoff (*TC*) [37,38,40]. The first version of MMF model [37] estimated soil erosion rate of an area by comparing the amount of soil particles detached by raindrop impact (*F*) and transport capacity of surface runoff (*TC*). The second version of model, the revised Morgan–Morgan–Finney (RMMF) model [38] started to consider the amount of soil particles generated by surface runoff (*H*). In the third version, the modified Morgan–Morgan–Finney (MMMF) model [40], the interconnectivity of surface runoff, various sub-processes such as the subsurface interflow and gravitational deposition processes, and parameters such as the physical structure of vegetation and surface ground conditions were introduced to calculate transport capacity of surface runoff (*TC*) and the amount of soil particles available for transport (*G*) more physically rigorously [41]. The daily based Morgan–Morgan–Finney (DMMF) soil erosion model [15] is also estimates daily soil loss from an element by comparing transport capacity of surface runoff (*TC*) and the available sediment for transport (*G*). The DMMF model is mainly comprised of hydrological and sediment phases. The hydrological phase determines the amount of surface runoff and subsurface interflow, and the sediment phase determines the amount of sediment budgets of the element.

#### *Appendix A.1. Hydrological Phase*

The effective rainfall (*Reff* ; mm) which is the volume of rainfall reaching the unit surface area of an element is the main driver of hydrological phase. Following the corrected version of the effective rainfall (*Reff*) from Choi et al. [42], *Reff* is calculated as,

*Reff* = *R* × (1 − *PI*) × cos(*S*), (A1)

where *PI* is the proportion of the permanent interception area and *S* is the slope of an element. Similar to MMF model, surface runoff can be generated when the total input of water to the element exceeds the surface water infiltration capacity (*SWc*; mm), which is the soil moisture storage capacity considering the proportion of the impervious area (*IMP*). *SWc* is defined as,

$$\mathcal{SW}\_{\mathbb{C}} = (1 - Imp) \times \left( \mathcal{SW}\_{\text{sat}} - \mathcal{SW}\_{\text{init}} - \frac{\Sigma IF\_{\text{in}}}{A} \right), \tag{A2}$$

where *SWsat* (mm) is the volume of water per unit area when soil is fully saturated, and *SWinit* (mm) is the volume of initial water per unit area that is already existed in the soil. Σ*IFin* (L) is the volume of subsurface water inputs from upslope and *A* (m2) is the area of an element. The amount of the surface runoff (*Q*; mm) is calculated as,

$$Q = R\_{eff} + \frac{\Sigma Q\_{in}}{A} - S\mathcal{W}\_{c\ \prime} \tag{A3}$$

where *Qin* (L is the volume of surface runoff inflow from upslope areas. The amount of water in the soil also flows out from the element as a subsurface interflow (*IFout*; L) when the voludme of soil water budget per unit area (*SW*; mm) of the element exceeds the volume of soil water at field capacity per unit area (*SWf c*; mm). The soil water budget (*SW*) is estimated as,

$$SW = (SW\_{init} + \frac{\Sigma IF\_{in}}{A}) + (R\_{eff} + \frac{\Sigma Q\_{in}}{A} - Q) - ET,\tag{A4}$$

where *ET* (mm) is the volume of water evapotranspirates per unit area from the element. Then the volume of subsurface water flowing out from the element (*IFout*) can be described as,

$$IF\_{\rm out} = K \times \sin(S) \times (SW - SW\_{\rm fc}) \times w \,\,\,\,\,\tag{A5}$$

where *K* (m/day) is the saturated soil lateral hydraulic conductivity and *w* (m) is the width of the element. A part of soil water remains with remaining water content (*θr*; vol/vol) which can be described as,

$$\theta\_r = \frac{\left(SW - IF\_{\text{out}}/A\right)}{1000 \times SD},\tag{A6}$$

where *SD* is the soil depth of the element, and 1000 is the constant to convert meters to millimeters. The *θ<sup>r</sup>* can be changed into *θinit* for the next day.

#### *Appendix A.2. Sediment Phase*

Sediment phase determines the total mass of soil particles which is taken out of the element through three steps: delivery of detached soil particles into the surface runoff, gravitational deposition, and estimation of hhe sediment loss from the element (*SL*) by comparing transport capacity of the runoff (*TC*; kg/m2) and sediment available for tranport (*G*; kg/m2). In the model, soil particles are detached from the surface by raindrop impact and surface runoff. The mass of soil particles detached by raindrops per unit area (*F*; kg/m2) is described as,

$$F = 0.001 \times DK \times P \times (1 - EPA) \times KE \, , \tag{A7}$$

where *DK* (g/J) is the detachability of soil particles by raindrop impact, *P* (%) is the proportion of each soil particle size class (i.e., clay, silt, and sand), *KE* (J/m2) is the kinetic energy of the effective rainfall considering direct throughfall and leaf drainage from the plant, and 0.001 is the unit conversion factor from g to kg. Also, *EPA* is the erosion protected area:

$$EPA = IMP + (1 - IMP) \times GC \,\, , \tag{A8}$$

where *GC* is the proportion of ground cover and *IMP* is the proportion of the impervious area (*IMP*) of the element. The mass of detached soil particles by the surface runoff (*H*; kg/m2) is described as,

$$H = 0.001 \times DR \times P \times Q^{1.5} \times (1 - EPA) \times (\sin(S))^{0.3},\tag{A9}$$

where *DR* (g/mm) is the detachability of soil particles by runoff per unit volume of surface runoff and *Q* is the volume of runoff per unit area, *S* is the slope of the element, and 0.001 is the unit conversion factor from g to kg. Sediment inputs from upslope elements (Σ*SLin*) also flows into surface runoff. The mass of delivered sediments to the surface runoff per unit area (*SS*; kg/m2) is,

$$SS = F + H + \frac{\Sigma SL\_{in}}{A} \,. \tag{A10}$$

A part of sediments delivered to the surface runoff (*SS*) in the runoff settle down to the ground by gravity. The gravitational deposition rate of the suspended sediments (*SS*) in runoff (*DEP*) is,

$$DEP = 0.441 \times N\_{f'} \,\tag{A11}$$

where *Nf* is the particle fall number which is the probabilistic ratio of falling particles [81], The *Nf* can be estimated as,

$$N\_f = \frac{l}{v} \times \frac{v\_s}{d} \, \, \, \, \, \tag{A12}$$

where *v* (m/s) is the velocity of the surface runoff, *vs* is the settling velocity of each particle size class, and *d* (m) is the depth of the surface runoff.

The remaining suspended sediments become available for transport per unit volume of surface runoff per unit area (*G*; kg/m2) and be estimated as,

$$G = SS \times \left(1 - DEP\right). \tag{A13}$$

The part of the availabe sediments for transport (*G*) can flow out from the element according to the transport capacity of the runoff (*TC*; kg/m2) of an element which is determined by the volume of runoff per unit area of an element (*Q*), the slope angle (*S*) and the surface conditions [40]. Due to the physical condition of surface affect runoff velocity, the tranport capacity of runoff can be described using the ratio between actual runoff velocity (*v*) and the reference velocity of the element (*vr*; m/s) [42].

$$TC = 0.001 \times \left(\frac{v}{v\_r}\right) \times Q^2 \times \sin(S) \,. \tag{A14}$$

The reference velocity (*vr*) is,

$$v\_r = \frac{1}{n\_r} \times d\_r^{2/3} \times \sqrt{\tan(S)}\,\,\,\,\,\tag{A15}$$

with 0.015 for Manning's coefficient (*nr*) and 0.005 for runoff depth (*dr*) representing for a standard surface condition. The transport capacity of the runoff (*TC*) and the available sediment for transport (*G*) determines the amount of sediment loss from the element (*SL*) [40,82]. When *TC* is greater than *G*, the surface runoff washes out all the sediments available for transport, otherwise, the amount of sediment (*SL*) which is equal to *TC* can be transported from the element.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
