A Novel Method for Watershed Best Management Practices Spatial Optimal Layout under Uncertainty

Abstract

Watershed Best management Practices (BMPs) spatial optimal layout would be affected by uncertainty, and there are still three problems which are worthy of studying in the present studies of watershed BMPs spatial optimal layout under uncertainty: (1) how to integrate multiple uncertainties in optimization model effectively; (2) how to avoid subjective weight in multi-objective uncertainty model; (3) how to develop more elastic schemes for uncertainty impact. To solve the mentioned problems, this study takes Zhegao river watershed, China as an example, interval stochastic fuzzy fractional programming (ISFFP) integrated with SWAT hydrology model is applied for BMPs spatial optimal layout in watershed to reduce non-point source (NPS) pollution. The result shows that the ISFFP method could solve the problems effectively, and the method could be adapted to different types of uncertainty, also the method has seldom been used in uncertainty BMPs spatial optimal layout, and the method is worth of popularization.

1. Introduction

The spatial optimization layout of watershed BMPs could alleviate watershed non-point source (NPS) pollution by setting the categories and number of BMP facilities in the watershed space optimally [1]. It is also proven to be an effective method for NPS pollution control in the watershed. However, uncertainty factor, which widely exists in nature and human society, influences the expected effect of watershed BMP spatial optimization layout [2]. Studying the watershed BMP spatial optimization layout under uncertainty is meaningful to reduce the uncertainty disturbance and ensure the treatment effect of the NPS pollution.

Literature Review

In many cases, the BMP spatial optimal allocation is a multi-objective problem. The problems which are considered by decision makers or planners always include how to maximize the pollution control effect and how to minimize the cost or the establishment. Two major methods are developed for the problem, one is plan-based method [3], and the other is the optimization algorithm-based method [4]. For the plan-based method, the schemes of BMP spatial layout are designed with expert knowledge or by the previous field studies. However, the method could not always achieve the optimal category selection and placement for BMPs. For the optimization algorithm-based method, with optimization programming and mathematical algorithm, the obtained schemes are more reasonable and are closer to the optimal results, and the methods include genetic algorithm (GA), non-dominated sorting genetic algorithm (NSGA-II), strength Pareto evolutionary algorithm 2 (SPEA2) and so on [5,6]. The methods of recent related studies are always hydrology models integrated with optimal algorithms. Hydrology models are used for simulating the distribution of NPS pollution, and optimal algorithms are used for determining optimal schemes on the basis of NPS pollution simulation. Several distributed hydrological model, such as SWAT (Soil and Water Assessment Tool), HSPF (Hydrological Simulation Program-FORTRAN), and AnnAGNPS (Annualized Agricultural Non-Point Source), have been widely used for hydrology and water quality simulation [7,8]. For example, SWAT could inform adaptive water management by facilitating quantitative analysis of different components of the water condition within a watershed. The model could simulate hydrological process and water quality under the impacts of water and land management practices [9,10]. The SWAT model always is integrated with GA, NSGA-II or SPEA2 for BMPs spatial optimization layout.

Uncertainties exist inevitably in nature and human society. From the mathematical viewpoint, uncertainty value could be classified into interval number, fuzzy number and stochastic number [11]. Uncertainty would influence the rationality of the scheme of BMPs spatial optimal allocation by affecting the NPS pollution treatment efficiency and the cost of BMPs facility. For example, the P treatment process of BMPs refers to many factors such as BMPs reaction time, temperature, season and plant species, microbial species, and others [12,13]. Variance in any one of them would cause changes in pollutant removal efficiency. The distribution pattern of most BMPs efficiency values were not evident through visual inspection, because the uncertainties which influence the pollution treatment efficiency of BMPs are always stochastic variances; therefore, it was always assumed to follow a normal distribution. In the related studies, the P treatment efficiency of BMPs which pass the normality test could be regarded as stochastic numbers of normal distribution, and the numbers which fail the test could be regarded as the interval number [14]. For the economic cost, the economic cost of BMPs is subject to market fluctuations, and the price always follows fuzzy number distribution [15].

Most BMP spatial optimal layouts are under certain condition, and very few related studies are under uncertainty. Four problems should be solved when conducting uncertainty BMP spatial optimal layout.

(1)

How to identify and represent uncertainty. Most of the uncertainties in BMP spatial optimal layout include the efficiency uncertainty in BMPs and the uncertainty in cost and budget. The root of efficiency uncertainty is variance factors, which include temperature, rainfall, season, vegetation form, and microbial species, and the distribution pattern of the uncertainty and the uncertainty interval values have been measured by the researchers using experiment methods in related studies [16]. In mathematical models, uncertainty could be in forms of stochastic number, fuzzy number, and interval number. Economic uncertainty stems from market fluctuation and the possible variance of government grant [17]. Economic uncertainty interval could be determined by analyzing or forecasting market price fluctuation. The price uncertainty always follows fuzzy distribution, and budget uncertainty is usually expressed in interval number or in multiple scenarios.

(2)

How to identify the uncertainty influence on BMPs spatial optimal layout. Uncertainty factors could be regarded as the independent variable in BMP spatial optimal layout, and the dependent variable could be seen as system pollution control efficiency and the total costs, when the optimal system follow linear mathematics characteristics, the distribution interval of the dependent variable could be obtained according to the independent variable [18]. However, in reality, the optimal systems are always complex non-linear mathematics model with interaction effects between the parameters, in this case, Monte Carlo method is usually used to judge the impact of uncertainty, through multiple and randomly set data within the interval of the independent variable, multiple corresponding results could be obtained, and then the distribution interval of system results could be considered [19].

(3)

How to integrate uncertainty into BMP spatial optimal layout. In the mathematics model for BMP spatial optimal layout model, uncertainty could be integrated into the optimal model with digital form and in uncertainty scenario form. Genetic algorithm (GA) is widely used for uncertainty in BMPs layout [20]. It can integrate uncertainty factors in mathematical form. However, with GAs, all the related variables need to be set as uncertain values, but in practice, except for part of the variables are uncertain numbers, the others are definite variables.

(4)

How to avoid subjective weight in multi-objective optimization for uncertainty BMP spatial optimal layout. The related studies always use Analytic Hierarchy Process (AHP) method or expert evaluation method to set the weight for each objective programming; however, subjectivity is inevitable in the setting [21].

Thus, the previous studies could identify and represent uncertainty, as well as identifying the influence of uncertainty on the results of BMP layout system effectively. Three problems deserve further study: (1) the problem of integrating multiple uncertainty factors; (2) the problem of avoiding subjective weight in multi-objective programming; (3) the problem of setting elastic scheme, Pareto approximation algorithm is widely used in the related multi-objective studies, although many solutions are produced in calculation processing, and the final chosen solution is the optimum result, the result is still unique. However, in reality, problems that the decision maker would modify in the scheme provisionally or the chosen BMP facility is limited by the geological constraint. Therefore, it is effectual to set an elastic interval for schemes to process the influence of uncertainty, that is, that the obtained schemes are reasonable and could manage the influence of uncertainty when the schemes are in the interval.

Based on it, the purpose of the study is to establish a universal method for watershed BMPs spatial optimal layout under uncertainty, and the three problems will be solved: (1) how to integrate multiple kinds of uncertainty into optimization model; (2) how to avoid subjective weight in multi-objective optimal model; (3) how to develop more elastic schemes to cope with the impact of uncertainty. The developed method could help the decision maker or the planner deal better with complex problems or uncertainty in BMPs spatial optimal layout.

2. Description of Area Studied

In this study, a small watershed named Zhegao river watershed (31°36′8″ N–31°55′40″ N, 117°34′43″ E–117°55′57″ E) is selected as the study area, Zhegao river watershed is located on the north side of Chaohu Basin, and it is also located in Hefei city, Anhui province, China. The total length of Zhegao River is 35 km. The river flows from north to south into Chaohu Lake. In this study, the middle and lower reaches of the Zhegao River (which has an area of 431.2 km²) are selected as the study area. This area is higher in the north and lower in the south. It is located in a plain polder area, which refers to a plain river network or lakeside and other low-lying waterlogging areas. This area is formed through embankments, sluices, and pumping stations. The purpose of constructing a polder area is to resist floods and waterlogging. The study area is characterized by a north subtropical humid monsoon climate, with annual precipitation of 1000–1158 mm, annual evaporation of 1469–1629 mm, and rainfall concentrated mostly in summer [22]. At present, the main part of the study area is agricultural land, and the rest consists of grassland, woodland, and residential land.

The eastern Chaohu Lake is in the state of mild eutrophication yearly. The water pollution of Chaohu Lake are mainly excessive total phosphorus and total nitrogen. A total number of 7 phylas, 42 generas and 93 species phytoplankton plants were identified in Chaohu Lake, among which green algae are the dominant species (50 species). The phytoplankton community structure is unstable and susceptible to the external environment. A total number of 4 phylas, 35 generas and 70 species zooplantkter were identified in Chaohu Lake, with the largest number of rotaworm species, 16 genera and 32 species. Chaohu Lake has abundant zooplankton diversity and good stable community structure [23].

The study area is listed in Figure 1.

The study area faces problems of eutrophication and algae growth, and P is a key factor. P comes from NPS pollution caused by local agricultural activities.

The purpose of this study is to reduce the total P release of the study area through BMP spatial optimal allocation. Taking July, when NPS pollution is worst, as an example, this study analyzes the schemes of BMP spatial allocation layout under different reduction targets of P pollution, specifically, 20%, 40%, and 60%.

The problems need to be solved in the research include: (1) setting objective of optimization layout: maximizing the NPS pollution control efficiency and minimizing the cost; (2) achieving the goal of NPS pollution control; (3) setting place for BMPs facility installation, and Selecting the category and the number of BMPs facility; (4) representing the uncertainty in the model; (5) developing a flexible optimal scheme with floating range. Additionally, the problems are the specific form of the problems of uncertainty BMPs spatial optimal layout in practical case.

3. Materials and Methods

To solve the aforementioned problems, the study develops a theoretical framework to process BMPs spatial optimal layout under uncertainty. The framework includes: (1) Applying hydrological model for simulating NPS pollution emission. (2) BMPs facility selection and uncertainty analysis. (3) Based on the analysis above, applying mathematical optimal model for BMPs spatial optimal layout under uncertainty.

The research framework and the flow diagram are detailed as follows:

The procedure of the study is detailed as follows:

(1)

Applying hydrological model for simulating NPS pollution emission. SWAT model is applied to discriminate the sub watershed in the study area and to simulate the distribution of NPS pollution. This task includes collecting basic data on the study area, establishing the basic database of the research area (e.g., DEM, soil, and land use), using the SWAT model to divide the study area into sub-basins, and simulating the P emissions in each sub-basin.

(2)

BMPs facility selection and uncertainty analysis. These tasks involve the selection of suitable BMPs facilities (e.g., vegetation buffer zones, Ponds system, wetland) and facility parameters (e.g., scale, depth, width and so on) and the analysis of the uncertainty of BMPs (cost uncertainty and P treatment efficacy uncertainty).

(3)

Applying mathematical optimal model for BMPs spatial optimal layout under uncertainty. The integrated interval stochastic fuzzy fractional programming (ISSFP) model is applied for the BMP spatial optimal layout under uncertainty. The results of optimal schemes and NPS pollution reduction effect and total cost could be obtained.

The flow diagram is listed in Figure 2.

ISFFP

Formally, a linear fractional programming is defined as the problem of maximizing/minimizing a ratio of affine function over a polyhedron, and it could be written as [24,25]:

$$\mathrm{maxmize}\frac{{\mathrm{c}}^{\mathrm{T}}\mathrm{x}+\mathsf{\alpha }}{{\mathrm{d}}^{\mathrm{T}}\mathrm{x}+\mathsf{\beta }}$$

(1)

$$\mathrm{Subject} \mathrm{to} \mathrm{Ax}\le \mathrm{b}$$

where $\mathrm{x}$ represents the vector of variables to be determined, $\mathrm{c},\mathrm{d}\in {\mathrm{R}}^{\mathrm{n}}$ are vectors of (know) coefficients, $\mathsf{\alpha },\mathsf{\beta }\in \mathrm{R}$ are constants, $\mathrm{A}\in {\mathrm{R}}^{\mathrm{m}\times \mathrm{n}}$ is a (know) matrix of coefficients, $\mathrm{b}\in {\mathrm{R}}^{\mathrm{m}}$ are vectors of (know) coefficients, $\mathrm{T}$ denotes the transpose.

The general uncertainty fractional model could be expressed as:

$${\mathrm{Maxf}}^{\pm }=\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{C}}_{\mathrm{i}}^{\pm }{\mathrm{x}}_{\mathrm{i}}^{\pm }+{\mathsf{\alpha }}^{\pm }\right)}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{d}}_{\mathrm{i}}^{\pm }{\mathrm{x}}_{\mathrm{i}}^{\pm }+{\mathsf{\beta }}^{\pm }\right)}$$

(2)

$$\begin{gathered} {\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{A}}_{\mathrm{ij}}^{\pm }{\mathrm{x}}_{\mathrm{i}}^{\pm }\le {\mathrm{B}}_{\mathrm{j}}^{\pm } \mathrm{i}=1,2\dots \mathrm{n}; \mathrm{j}=1,2\dots \mathrm{m} \\ {\mathrm{x}}_{\mathrm{i}}^{\pm }\ge 0 \mathrm{i}=1,2\dots \mathrm{n} \end{gathered}$$

where ${\mathrm{A}}_{\mathrm{ij}}^{\pm },{ \mathrm{B}}_{\mathrm{j}}^{\pm },{ \mathrm{C}}_{\mathrm{i}}^{\pm },{ \mathrm{x}}_{\mathrm{i}}^{\pm }$ are the uncertainties, and the uncertainties could be interval numbers or stochastic numbers. In the uncertainty model above, the two types of uncertainties are all in the interval number form with the deterministic and closed boundary, in which the stochastic number could be transformed into the interval number by its probability distribution, and the ‘$-$’ represents the lower bound, and the ‘$+$’ represents the upper bound.

Besides the interval number and stochastic number, another type of uncertainty number is fuzzy number. In the case of the three types of uncertainties appear in the objective function or in the constraints, the uncertainty optimization programming could be expressed as:

$${\mathrm{Maxf}}^{\pm }=\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left\{\left({\mathrm{C}}_{\mathrm{i}}^{\pm }+{\tilde{\mathrm{C}}}_{\mathrm{i}}^{\prime }\right)·{\mathrm{x}}_{\mathrm{i}}^{\pm }+{\mathsf{\alpha }}^{\pm }\right\}}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left\{\left({\mathrm{d}}_{\mathrm{i}}^{\pm }+{\tilde{\mathrm{d}}}_{\mathrm{i}}^{\prime }\right)·{\mathrm{x}}_{\mathrm{i}}^{\pm }+{\mathsf{\beta }}^{\pm }\right\}}$$

(3)

Subject to

$$\begin{gathered} {\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{A}}_{\mathrm{ij}}^{\pm }{\mathrm{x}}_{\mathrm{i}}^{\pm }\le {\mathrm{B}}_{\mathrm{j}}^{\pm } \mathrm{i}=1,2\dots \mathrm{n}; \mathrm{j}=1,2\dots \mathrm{m} \\ {\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\tilde{\mathrm{A}}}_{\mathrm{ig}}^{\prime }{\mathrm{x}}_{\mathrm{i}}^{\pm }\le {\mathrm{B}}_{\mathrm{g}}^{\pm } \mathrm{i}=1,2\dots \mathrm{n}; \mathrm{g}=1,2\dots \mathrm{m} \\ {\mathrm{x}}_{\mathrm{i}}^{\pm }\ge 0 \mathrm{i}=1,2\dots \mathrm{n} \end{gathered}$$

wherein ${\tilde{\mathrm{A}}}_{\mathrm{ig}}^{\prime }$ and ${\tilde{\mathrm{C}}}_{\mathrm{i}}^{\prime }$ are the fuzzy numbers, and ${\tilde{\mathrm{A}}}_{\mathrm{ig}}^{\prime }$ and ${\tilde{\mathrm{C}}}_{\mathrm{i}}^{\prime }$ could be triangle fuzzy numbers or trapezoidal fuzzy numbers. If ${\tilde{\mathrm{A}}}_{\mathrm{ig}}^{\prime }$ and ${\tilde{\mathrm{C}}}_{\mathrm{i}}^{\prime }$ are trapezoidal, the fuzzy numbers could be expressed as ${\tilde{\mathrm{A}}}_{\mathrm{ig}}^{\prime }=\left({\mathrm{A}}_{\mathrm{ig}}^{-},{ \mathrm{A}}_{\mathrm{ig}1},{ \mathrm{A}}_{\mathrm{ig}2},{ \mathrm{A}}_{\mathrm{ig}}^{+}\right)$ and ${\tilde{\mathrm{C}}}_{\mathrm{i}}^{\prime }=\left({\mathrm{C}}_{\mathrm{i}}^{-},{ \mathrm{C}}_{\mathrm{i}1},{ \mathrm{C}}_{\mathrm{i}2},{ \mathrm{C}}_{\mathrm{i}}^{+}\right)$.

The $\mathsf{\alpha }\text{-}\mathrm{cut}$ method could be used to represent the distribution interval of fuzzy number. The distribution interval of fuzzy numbers ${\tilde{\mathrm{A}}}_{\mathrm{ig}}^{\prime }$ and ${\tilde{\mathrm{C}}}_{\mathrm{i}}^{\prime }$ could be represented as $\left[\left(1-\mathsf{\alpha }\right)·{\mathrm{A}}_{\mathrm{ig}}^{-}+\mathsf{\alpha }·{\mathrm{A}}_{\mathrm{ig}1}, \left(1-\mathsf{\alpha }\right)·{\mathrm{A}}_{\mathrm{ig}}^{+}+\mathsf{\alpha }·{\mathrm{A}}_{\mathrm{ig}2}\right]$ and $\left[\left(1-\mathsf{\alpha }\right)·{\mathrm{C}}_{\mathrm{i}}^{-}+\mathsf{\alpha }·{\mathrm{C}}_{\mathrm{i}1}, \left(1-\mathsf{\alpha }\right)·{\mathrm{C}}_{\mathrm{i}}^{+}+\mathsf{\alpha }·{\mathrm{C}}_{\mathrm{i}2}\right]$. Wherein $\mathsf{\alpha }$ is the membership.

4. Study Process

4.1. Discriminating the Sub Watershed and Simulating P Distribution in Each Sub Watershed

4.1.1. Model Introduction

The SWAT model, which was used many times in related studies, is applied in the study [26]. The SWAT model can divide the studied watershed into a sub-watershed and simulate the transformation and migration of P in each sub watershed. The SWAT model could calculate the amount of P emission and P reduction, P emission is mainly from inorganic fertilizers and manure used in agricultural activity, and P reduction is mainly from crop absorption, surface runoff, flow measurement, infiltration, and soil erosion [27]. The SWAT model can fully consider the changes in regional land use, and soil and agricultural tillage in the simulation of P transformation and migration [28].

4.1.2. Database Preparation

The SWAT model consists of spatial and meteorological databases. The spatial database includes a digital elevation model (DEM), soil, land use, and others. The elevation data are the 30 m resolution DEM data provided by the international scientific data service platform. The land-use database was derived from the interpretation data of watershed Landsat TM images in 2018. There are 6 types of land use in the river basin in 2018: these are paddy fields, dry fields, woodland, grassland, water bodies, residential areas, and paddy fields in the main area. The soil database uses the soil type (1:1 million) provided by Nanjing soil as input data for the simulation. The Zhegao river basin consists of yellow brown, yellow brown loam, paddy, coarse bone, limestone, and rinsing paddy soil. The main area mainly consists of paddy soil. The meteorological attribute database consists of the data obtained from a local meteorological observatory, including daily water drop, daily maximum/minimum temperature, daily solar radiation, wind speed, and daily average relative humidity.

4.1.3. Spatial Analysis of NPS Pollution Emission

(1)

Discriminating sub watershed. DEM is used for discriminating river system and sub watershed [29].

(2)

Discriminating hydrologic response unit (HRU). Land use type, soil type, slope and so on are used for discriminating HRU in each sub watershed. The amount of P emission in each HRU could be calculated by SWAT model, and the total amount of P emission in each sub watershed could also be get.

4.1.4. Parameter Validation and Calibration

When no observation data are available, the calibration and validation of the parameters could refer to the parameters of other watersheds at the same latitude and natural conditions [30]. As there is no long-term effective hydrological observation data in the study area, this study transplanted the parameters of the Xinanjiang River Basin under the same latitude and natural conditions [31]. At the same time, the previous studies in the case area are considered as a reference [32], and the relevant parameters needed for the simulation of NPS P pollution are finally obtained (

Table 1. Validated value of SWAT modeling parameters.

Parameter		Validated Value	Parameter	Validated Value	Parameter	Validated Value
CN2.mgt	FRST	59.54	CANMX.hru	12.616	ALPHA_BF.gw	0.1282
	RICE	71.88	ESCO.hru	0.1188	GW_DELAY.gw	36.182
	PAST	61.12	GWQMN.gw	1274.91	SMTMP.bsn	0.0418
	URML	58.06	Usle.mgt	0.1367	CMN.bsn	0.002
	WATR	85.33	Spexp.bsn	0.76	NPERCO.bsn	0.26
	AGRL	73.00	SPCON.bsn	0.0543	PSP.bsn	0.7
PPERCO.bsn		15	CH_K2.rte	38.951	BC2.bsn	2
BC4.bsn		0.01	TIMP.bsn	0.29	BC3.bsn	0.23
AI1.wwq		0.08	PHOSKD.bsn	165	AI2.wwq	0.02
CH_N2.rte		0.1597	BC1.bsn	0.22	RCDCO.bsn	0.05

).

4.1.5. Relevant Data of Sub-Watershed

Through SWAT simulation, the study area is divided into 31 sub-basins. The area of each sub-basin, surface water quantity, and P release amount are shown in

Table 2. The area of each sub-basin, surface water quantity, and P release amount in each basin.

Serial Number	Area (ha)	Volume (m3)	P Emission (kg)
1	2905.30	949,596.18	13,460.24
2	2735.05	812,091.78	11,301.24
3	1.15	472.45	4.73
4	372.37	208,756.04	1399.38
5	2551.61	863,261.32	11,665.97
6	234.10	71,532.67	1402.94
7	4097.82	1,034,086.06	17,456.73
8	1501.46	429,672.59	6184.51
9	64.43	25,940.94	427.14
10	3476.11	746,042.74	14,842.99
11	131.39	49,412.60	1044.83
12	1829.82	481,754.04	7434.54
13	2754.15	404,336.85	6742.16
14	3053.00	730,766.05	12,886.71
15	652.70	198,354.89	3957.96
16	1374.90	448,548.42	8901.12
17	1461.62	509,054.26	7154.65
18	1465.23	324,709.64	6275.58
19	1373.02	443,786.87	7658.69
20	1.72	708.07	8.65
21	986.46	162,490.48	3440.79
22	1004.91	499,378.59	3939.24
23	1735.96	692,615.28	10,700.49
24	1641.29	657,912.61	8895.81
25	500.81	203,230.43	3374.49
26	1390.07	555,748.78	8261.17
27	1678.92	668,359.80	8780.73
28	170.41	69,879.21	1317.25
29	36.97	15,198.90	270.08
30	1199.74	488,726.56	7341.22
31	741.22	299,609.15	4965.44

.

4.2. BMP Selection and Relevant Uncertainty Analysis

4.2.1. Analysis of BMP Character

Based on the analysis of preceding research, three types of BMPs vegetation buffer zone, ponds system and constructed wetlands are selected as the study objects. Because there is heavy rain and flat terrain in the study area, and the three BMPs are suitable for layout for rainwater processing.

(1)

Vegetation buffer zone

The shore vegetation buffer zone is usually composed of trees and other vegetation that climb to the slopes on both sides of the riverbank. The main ways to intercept P in the vegetation filtration zone include plant absorption and soil adsorption [33]. The way for plants to remove P is root absorption. Sediment removes P by absorbing it, and vegetation filtration zones can intercept P, especially granular P, by intercepting sediment [34,35,36].

(2)

Pond system

A pond system uses soil microbial plant system to intercept, deposit, absorb, and transform P through physical, chemical, and biological processes to achieve efficient purification of P pollution [37]. Furthermore, the pond system promotes the growth of green plants and achieves resource utilization and innocuity of P through the biogeochemical cycle of nutrients and water [38].

(3)

Constructed wetlands

Constructed wetlands remove P through the combined action of substrates, aquatic plants, and microorganisms [39,40]. The substrate is the filler, and its main way to remove P is through adsorption, that is, when runoff flows through the constructed wetland, the substrate purifies and removes P from the runoff through certain physical and chemical pathways such as absorption, filtration, ion exchange, and complexation reaction [41]. Aquatic plants can transfer inorganic P to organic components of plants through plant absorption and assimilation, and the P absorbed by plants can be removed by regular harvesting. The micro-organisms can transfer the organ P into P phosphate and also increase the solubility of the organic P through the metabolic activity of P bacteria. In this way, the P in the runoff is removed [42].

4.2.2. Relevant Parameters of BMPs

The volume of NPS treatment of the wetland and pond system depends on its volume and permeability, and the vegetation buffer zone depends on its surface area. In this study, five different scales of BMP facilities were set for the three types of BMPs. The widths of the coastal vegetation zones are 2, 5, 10, 15, and 20 m. With the increase of the width, the area of the vegetation buffer increase correspondingly, and it would affect crop planting if the width is too wide, then the upper limit of the buffer’s width is set as 20 m. The basis of width setting refers to the related study [43].

The depth of the Ponds system is set to 1.6 m. The number of surface areas is set to 5, namely, 0.2%* A_sub.i, 0.4%* A_sub.i, 0.6%* A_sub.i, 0.8%* A_sub.i, and 1.0%* A_sub.i, where A_sub.i represents the area of the i-th subbasin. The depth of the wetland is set to 0.7 m. The number of surface area types is set to 5, namely, 0.2%* A_sub.i, 0.4%* A_sub.i, 0.6%* A_sub.i, 0.8%* A_sub.i, and 1.0%* A_sub.i, where A_sub.i represents the area of the i-th subbasin. The depth of the ponds system and wetlands are reference to the related study [44,45]. (Note: The area of the pond system and wetland does not cover only one-unit BMPs, but the total area of the BMPs is no more than the set value).

Although the capacities of ponds and wetlands are constant values, and the ponds and wetlands involve evaporation and infiltration of water bodies, the formula for calculating the volume of treated water in the ponds and wetlands is as follows: BMP surface area × (depth of facility + permeability of facility + evaporation).

The water evaporation in July is 0.134 m, the permeability of the Pond is 3.25 m, and that of the wetland is 1.38 m [46].

The categories of each BMPs are listed in

Table 3. The categories of each BMPs.

Category	1	2	3	4	5
Vegetation buffer	2 m	5 m	10 m	15 m	20 m
Ponds system	0.2%* Asub.i	0.4%* Asub.i	0.6%* Asub.i	0.8%* Asub.i	1.0%* Asub.i
Wetlands	0.2%* Asub.i	0.4%* Asub.i	0.6%* Asub.i	0.8%* Asub.i	1.0%* Asub.i

.

Additionally, the volume of P treatment of the ponds system and the wetlands are calculated as follows

$${\mathrm{Tr}}_{\mathrm{j}}={\mathrm{ReP}}_{\mathrm{j}}·{\mathrm{EfP}}_{\mathrm{j}}$$

(4)

$${\mathrm{ReP}}_{\mathrm{j}}={\mathrm{TP}}_{\mathrm{i}}·{\mathrm{ReS}}_{\mathrm{j}}/{\mathrm{TS}}_{\mathrm{i}}$$

(5)

$${\mathrm{ReS}}_{\mathrm{j}}={\mathrm{A}}_{\mathrm{j}}·\left({\mathrm{DeP}}_{\mathrm{j}}+\mathrm{Ev}+{\mathrm{In}}_{\mathrm{j}}\right)$$

(6)

where

${\mathrm{Tr}}_{\mathrm{j}}$	The volume of treated P of j-th BMPs
${\mathrm{ReP}}_{\mathrm{j}}$	The volume of retained P of j-th BMPs
${\mathrm{EfP}}_{\mathrm{j}}$	The efficiency of P treatment of j-th BMPs
${\mathrm{TP}}_{\mathrm{i}}$	Total volume of P emission in i-th sub basin
${\mathrm{ReS}}_{\mathrm{j}}$	The retained surface water of j-th BMPs
${\mathrm{TS}}_{\mathrm{i}}$	Total volume of surface water in i-th sub basin
${\mathrm{A}}_{\mathrm{j}}$	The area of j-th BMPs
${\mathrm{Dep}}_{\mathrm{j}}$	The depth of j-th BMPs
$\mathrm{Ev}$	The evaporation rate of the study area in a month
${\mathrm{In}}_{\mathrm{j}}$	The infiltration rate of j-th BMPs in a month

4.2.3. The Uncertainty Analysis

Different type of uncertainty corresponds to the different distribution interval. The distribution interval of stochastic number depends on the confidence probability setting, and the interval of fuzzy number depends on the α-cut setting.

The P treatment efficiency of BMPs could be regarded as the variables with stochastic number and interval number properties. The distribution of P process efficiency of BMP is listed in accordance with relevant studies [47,48]. The uncertain distribution of cost is listed in accordance with relevant studies [49]. Additionally, the confidence probability and α-cut in the study are all set as a single value.

The uncertainty of P treatment efficiency and cost of the three BMPs are listed in

Table 4. The uncertainty of P treatment efficiency and cost of the three BMPs.

	P Treatment Efficiency	Cost (EUR/ha)
Ponds system	80~90%	[14.7, 34.3]
Wetlands	25~90%	[53.9, 80.7]
Vegetation buffer (2 m)	30.09~38.06%	[451.66, 1053.86]
Vegetation buffer (5 m)	41.00~50.92%	[451.66, 1053.86]
Vegetation buffer (10 m)	54.22~65.68%	[451.66, 1053.86]
Vegetation buffer (15 m)	69.00~78.95%	[451.66, 1053.86]
Vegetation buffer (20 m)	77.31~82.96%	[451.66, 1053.86]

.

Additionally, the volume of P which could be processed by BMPs in different grades could be calculated according to the formula and the related parameters. The values are listed in attached tables.

4.3. Construction of Optimization Model

The objective of the model is P treatment maximization and cost minimization under uncertainty. The constraints of the model include total volume of P treatment, allowed area of BMP installation, and number of installed BMPs in each sub-watershed.

4.3.1. Objective Function

P treatment maximization:

$$\mathrm{Maximize}={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{m}}{\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{TP}}_{\mathrm{i}}·{\mathrm{x}}_{\mathrm{ij}}·{\mathrm{TrE}}_{\mathrm{ij}}^{\pm }\hspace{1em}\mathrm{i}=1,2,\dots \mathrm{m}; \mathrm{j}=1,2,\dots \mathrm{n}$$

(7)

${\mathrm{TP}}_{\mathrm{i}}$: The total volume of P emission in the i-th sub basin.

${\mathrm{Tr}}_{\mathrm{ij}}$: The P treatment efficiency of the j-th BMPs facility in the i-th sub basin.

Cost minimization:

$$\mathrm{Maximize}={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{m}}{\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{x}}_{\mathrm{ij}}·{\mathrm{C}}_{\mathrm{j}}^{\pm }\hspace{1em}\mathrm{i}=1,2,\dots \mathrm{m}; \mathrm{j}=1,2,\dots \mathrm{n}$$

(8)

${\mathrm{C}}_{\mathrm{j}}$: the cost of the j-th BMPs facility.

The two above objective functions include maximum function and minimum function, and the two functions could be integrated by fractional programming. The integrated function is as follows.

$$\mathrm{Maxf}\left(\mathrm{x}\right)=\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{m}}{{\displaystyle \sum }}_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{TP}}_{\mathrm{i}}·{\mathrm{x}}_{\mathrm{ij}}·{\mathrm{TrE}}_{\mathrm{ij}}^{\pm }}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{m}}{{\displaystyle \sum }}_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{x}}_{\mathrm{ij}}·{\mathrm{C}}_{\mathrm{j}}^{\pm }}\hspace{1em}\mathrm{i}=1,2,\dots \mathrm{m}; \mathrm{j}=1,2,\dots \mathrm{n}$$

(9)

4.3.2. Constraint Function

The volume of the treated P in the i-th sub basin cannot be less than the specific proportion of the total volume of P emission in the same sub basin.

$${\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{m}}{\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{TP}}_{\mathrm{i}}·{\mathrm{x}}_{\mathrm{ij}}·{\mathrm{TrE}}_{\mathrm{ij}}^{\pm }\ge {\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{TP}}_{\mathrm{i}}·\mathrm{SP}\hspace{1em}\mathrm{i}=1,2,\dots \mathrm{m}; \mathrm{j}=1,2,\dots \mathrm{n}$$

(10)

SP: specific proportion of the total volume of P emission in the i-th sub basin, and the SP in the study are 20%, 40% and 60%, respectively.

The volume of the treated P in the i-th sub basin cannot be more than the total volume of P emission in the same sub basin.

$${\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{m}}{\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{TP}}_{\mathrm{i}}·{\mathrm{x}}_{\mathrm{ij}}·{\mathrm{TrE}}_{\mathrm{ij}}^{\pm }\le {\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{TP}}_{\mathrm{i}}\hspace{1em}\mathrm{i}=1,2,\dots \mathrm{m}; \mathrm{j}=1,2,\dots \mathrm{n}$$

(11)

SP: specific proportion of the total volume of P emission in the i-th sub basin.

The allowed area of BMPs installation in the i-th sub basin cannot be more than 1% of the area of the same sub basin.

$${\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{x}}_{\mathrm{ij}}·{\mathrm{A}}_{\mathrm{ij}}\le 1\%·{\mathrm{A}}_{\mathrm{sub}.\mathrm{i}}\hspace{1em}\mathrm{i}=1,2,\dots \mathrm{m}; \mathrm{j}=1,2,\dots \mathrm{n}$$

(12)

${\mathrm{A}}_{\mathrm{sub}.\mathrm{i}}$: the area of the i-th sub basin.

None negative for the variables

$${\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{x}}_{\mathrm{ij}}\ge 0\hspace{1em}\mathrm{i}=1,2,\dots \mathrm{m}; \mathrm{j}=1,2,\dots \mathrm{n}$$

(13)

01 setting for BMPs installation

$${\mathrm{x}}_{\mathrm{ij}}=0|1\hspace{1em}\mathrm{i}=1,2,\dots \mathrm{m}; \mathrm{j}=1,2,\dots \mathrm{n}$$

(14)

4.4. Uncertainty Scenario Analysis

Three scenarios of P reduction objectives which include 20%, 40%, 60% P reduction target were set in the study, and the upper limit scenario as well as the lower limit scenario are set in each scenario of P reduction target.

The upper limit value of the BMPs facility cost and the lower limit value of P pollution control are considered in the lower limit scenario, while the values are opposite in the upper scenario. The lower limit scenario and the upper limit scenario represent the boundaries of the uncertainty results.

5. Results

(1)

Through integrating ISFFP model and SWAT, the uncertainty BMP spatial optimization layout schemes are obtained. The schemes include the upper scenarios and the lower scenarios for 20%, 40%, and 60% P reduction targets in July. In each sub watershed, the number of the allocated BMPs is only one (the results are shown in

Table 5. The spatial optimization layout of BMPs in each scenario.

No.	20H	20L	40H	40L	60H	60L
1	Pond 1%	Pond 1%	Buffer 20 m	Buffer 20 m	Pond 1%	Buffer 15 m
2	Pond 1%	Pond 1%	Pond 1%	Pond 1%	Pond 1%	Buffer 20 m
3	Buffer 20 m	Buffer 20 m	Pond 1%	Pond 1%	Pond 1%	Pond 1%
4	Buffer 10 m	Buffer 20 m	Pond 1%	Buffer 20 m	Buffer 20 m	Buffer 20 m
5	Wet 0.2%	Pond 1%	Pond 1%	Pond 1%	Pond 1%	Buffer 20 m
6	Pond 1%	Buffer 20 m	Pond 1%	Buffer 20 m	Buffer 20 m	Buffer 15 m
7	Pond 1%	Pond 1%	Pond 1%	Pond 1%	Pond 1%	Pond 1%
8	Pond 1%	Pond 1%	Pond 0.8%	Buffer 20 m	Pond 1%	Buffer 20 m
9	Buffer 20 m	Pond 1%	Pond 0.8%	Pond 1%	Buffer 20 m	Pond 1%
10	Pond 1%	Pond 1%	Pond 1%	Pond 0.8%	Pond 1%	Pond 1%
11	Buffer 15 m	Pond 1%	Pond 1%	Pond 1%	Buffer 20 m	Buffer 20 m
12	Pond 1%	Pond 1%	Buffer 15 m	Pond 1%	Buffer 20 m	Pond 1%
13	Pond 1%	Pond 1%	Pond 1%	Pond 1%	Pond 1%	Pond 1%
14	Pond 1%	Pond 1%	Pond 1%	Pond 1%	Pond 1%	Buffer 20 m
15	Pond 1%	Pond 1%	Pond 1%	Pond 1%	Pond 1%	Pond 1%
16	Pond 1%	Pond 1%	Pond 1%	Pond 1%	Pond 1%	Buffer 20 m
17	Pond 1%	Pond 1%	Buffer 20 m	Pond 0.8%	Buffer 20 m	Buffer 20 m
18	Pond 1%	Pond 1%	Pond 1%	Pond 1%	Pond 1%	Buffer 15 m
19	Pond 1%	Pond 1%	Buffer 20 m	Pond 1%	Buffer 20 m	Buffer 20 m
20	Pond 1%	Pond 1%	Pond 1%	Pond 1%	Pond 1%	Buffer 20 m
21	Pond 1%	Pond 1%	Pond 1%	Pond 1%	Pond 1%	Buffer 20 m
22	Buffer 20 m	Buffer 20 m	Buffer 20 m	Buffer 20 m	Buffer 20 m	Pond 1%
23	Pond 1%	Pond 1%	Buffer 20 m	Buffer 20 m	Buffer 20 m	Buffer 20 m
24	Pond 1%	Pond 1%	Buffer 20 m	Buffer 20 m	Pond 1%	Buffer 20 m
25	Pond 1%	Pond 1%	Pond 1%	Buffer 15 m	Buffer 15 m	Buffer 15 m
26	Pond 1%	Pond 1%	Buffer 20 m	Buffer 20 m	Buffer 20 m	Buffer 20 m
27	Pond 1%	Pond 1%	Buffer 15 m	Buffer 15 m	Buffer 20 m	Buffer 20 m
28	Buffer 20 m	Pond 1%	Buffer 20 m	Buffer 20 m	Buffer 20 m	Buffer 20 m
29	Buffer 20 m	Buffer 20 m	Buffer 20 m	Buffer 20 m	Pond 1%	Buffer 20 m
30	Pond 1%	Buffer 20 m	Buffer 15 m	Buffer 20 m	Buffer 20 m	Buffer 20 m
31	Buffer 20 m	Buffer 20 m	Pond 1%	Buffer 20 m	Pond 1%	Buffer 20 m

Note: The number of BMPs installed in each sub basin is 1, respectively.

). The total cost and total volume of P under each scenario are shown in

Table 6. The total cost and the total volume of the treated P.

	20%		40%		60%
	+	−	+	−	+	−
Total cost (EUR)	3990.58	6585.58	5699.22	10,955.14	6573.14	17,509.48
Total treated P (kg)	38,904.15	40,483.60	89,900.35	80,739.23	80,622.07	121,235.31

. Statistics on the number of all BMP facilities under each scenario are shown in

Table 7. The number of installed BMPs in each scenario.

Type of BMPs		20%		40%		60%
		+	−	+	−	+	−
wetland	0.20%	1	0	0	0	0	0
	0.40%	0	0	0	0	0	0
	0.60%	0	0	0	0	0	0
	0.80%	0	0	0	0	0	0
	1.00%	0	0	0	0	0	0
Pond	0.20%	0	0	0	0	0	0
	0.40%	0	0	0	0	0	0
	0.60%	0	0	0	0	0	0
	0.80%	0	0	2	3	0	0
	1.00%	22	23	17	14	17	8
Vegetation buffer	2 m	0	0	0	0	0	0
	5 m	0	0	0	0	0	0
	10 m	1	0	0	0	0	0
	15 m	1	0	3	2	0	4
	20 m	6	8	9	12	14	19
Total number of BMPs		31	31	31	31	31	31

and Figure 3.

(2)

Table 4 shows that under the condition of the same area, whether in the upper or lower limit, according to the P treatment capacity, green buffer zone > Ponds system > wetland. Taking the sub watershed 1 as an example,

Table 8. The BMP treatment effect of sub watershed 1.

Type of BMPs	Scenario	0.2%	0.4%	0.6%	0.8%	1%
Pond	−	328.14	656.27	984.41	1312.54	1640.68
	+	369.15	738.31	1107.46	1476.61	1845.76
Wetland	−	45.59	91.18	136.76	182.35	227.94
	+	164.12	328.23	492.35	656.47	820.59
Vegetation buffer	−	4038.07	5518.70	7268.53	9287.56	10,364.38
	+	5114.89	6864.72	8883.76	10,633.59	11,172.00

corresponds to the BMP treatment effect of sub watershed 1, and the P control effect of each BMP facility can be observed.

Table 9. The ratio of P treatment capacity of the green buffer zone to the wetland and pond system.

Ratio		0.2%	0.4%	0.6%	0.8%	1%
Vegetation buffer/Wetland	−	88.57	60.53	53.15	50.93	45.47
	+	31.17	20.91	18.04	16.20	13.61
Vegetation buffer/Pond	−	12.31	8.41	7.38	7.08	6.32
	+	13.86	9.30	8.02	7.20	6.05

shows the ratio of P treatment capacity of the green buffer zone to the wetland and ponds system.

Table 8 and Table 9 show that the P treatment capacity of the green buffer zone is much larger than that of the Ponds system and wetland under the same scenario. A 1% scenario is used as an example, the P treatment capacity of the green buffer zone is 6.05 times as much the capacity of the pond in the upper limit scenario and 6.32 times in the lower limit scenario. The ratio of the green buffer zone to wetland is 13.61 times in the upper limit and 45.47 times in the lower limit.

6. Discussion

(1)

ISFFP model is used effectively to solve the three problems of BMPs multi objective spatial optimization layout in introduction section. (1) ISFFP could reflects and integrates the uncertainty factors which are in the forms of stochastic, fuzzy and interval, and the uncertainties mainly existed in P treatment efficiency and economic cost of BMPs. (2) Subjective problem, which existed in the weight setting between multi objective programming, could be avoided by using ISFFP. (3) The schemes, which are in the interval form, are get in the different scenarios, and the intervals represent the range of the reasonable schemes.

(2)

The other specific problems in the case are also solved by ISFFP model. (1) The objectives of P emission treatment maximization and cost minimization are all achieved. (2) The targets of the amount of the P treatment are all reached or exceed in each scenario, and the total treated P are [38,904.15, 40,483.60] kg in 20% scenario, [89,900.35, 80,739.23] kg in 40% scenario, and [80,622.07, 121,235.31] kg in 60% scenario, respectively. (3) The specific amounts and types of BMPs are allocated in each sub watershed.

(3)

ISFFP model is developed for BMPs spatial optimization layout, the results shows that the different schemes for BMPs spatial optimization layout are developed according to the different objectives of water environment treatment, and also all of the objectives are achieved in the study.

(4)

With the increase of the targets of water environment treatment, the more BMPs facilities with higher P treatment efficiency as well as higher costs are applied in BMPs spatial layout schemes, and total costs increase accordingly.

(5)

In the study area, there are different amounts of P emission in each sub watershed, and the types of the installed BMPs in each sub watershed are not the same. The P emission depends on the features of each sub watershed, and the features include area, agrotype, land type and so on. The BMPs with more P treatment efficiency are installed in the sub watershed with the higher P emission.

(6)

In this study, 31 BMPs were assigned under each scenario. However, no upper limit is set for the installed number and the types of BMPs per subbasin in the constraints. The reason is that one-unit BMP is sufficient to cope with the P pollution of the subbasin. Thus, it does not need more than one-unit BMP in a subbasin, and it can be attributed to the strong pollution control ability of the green buffer zone, which can absorb 77.31–82.96% of P. The pollution control target, which is under this limit, could be achieved through 20 m of green buffer zones. This study does not limit the budget. In the lower limit scenario of 60% P treatment scenario, the highest budget is 17,509.48 EUR. If the budget is lower than this amount, then it will limit the setting of the high-cost green buffer zone. This condition would also lead to the simulation results in which more ponds and wetlands are installed, and three BMPs are low effective but also low cost of low-effective but also low-cost ponds and wetlands.

(7)

According to the results, with the increase of the P pollution reduction target, the number of installed BMPs with higher pollution control effect is increase. In the lower limit scenarios of 20%, 40%, and 60% P reduction target, the number of green buffer zones installed are 6, 9, and 14, respectively. In the upper limit scenario, the assigned amounts are 8, 12, and 19. As the total number of BMPs in each scheme is the same, the number of green buffer zones with the highest treatment effect increases, and the number of other types of BMP decreases accordingly.

(8)

As the upper limit value of the BMP facility cost and the lower limit value of P pollution control quantity are considered in the lower limit scenario, then on the premise of completing the pollution control target of each grade, the total cost in the upper limit scenario is less than the total cost in the lower limit scenario.

In the 20%, 40%, and 60% scenarios, the total cost of the upper and lower limit scenarios are [6585.58 (lower scenario), 3990.58 (upper scenario)], [10,955.14 (lower scenario), 5699.22 (upper scenario)], and [17,509.48 (lower scenario), 6573.14 (upper scenario)], respectively.

(9)

As the lower limit of the BMP pollution control efficiency is set in the lower limit scenario, and the upper limit is set in the upper limit scenario, the number of green buffer zones with high cost and high pollution control effect in the upper limit scenario under each control target is also smaller than that in the lower limit scenario. This condition means that solving the same volume requires a smaller number of BMP facilities with high P pollution treatment efficiency and high cost.

(10)

As the cost and P treatment efficacy in this study are considered as the upper and lower limit values, the uncertain values of the cost and treatment efficiency of BMPs are also expressed as interval numbers. Thus, the results are in the form of the upper and lower limit scenarios. The results represent the upper and lower limits of the corresponding schemes under uncertainty impact, which means that the schemes are reasonable when their results are in the interval range.

(11)

The developed method could be extended to other areas. The SWAT model could be applied to simulate NPS emission in different types of land, and the ISFFP model could be used to reflected multiple types of uncertainties, and the model could change according to the types of uncertainties in actual situations.

(12)

The limitation of the study includes two aspects, one is that only two objectives of maximization and minimization are set, and the model would not avoid subjective weight setting completely when there are more objectives such as minimizing land use for BMPs setting, maximizing N emission treatment, and so on. How to apply more appropriate model for the problem would be studied in the future. The other is that the scheme is allocating the appropriate BMPs in each sub watershed; however, in practice, many other problems need to be considered, such as the natural condition for BMPs setting, cultivated land occupied for the setting, and so on, and the related problems also needs considering.

7. Conclusions

(1)

The innovation of the study is that it introduces a new method for solving the uncertainty in BMP spatial optimal layout, and the ISFFP method integrated with SWAT model has rarely been used for this purpose. The advantages of the method include the following: (1) it can reflect multiple uncertainty characters; (2) it could process the weight setting in maximum and minimum programming; and (3) it could achieve flexible schemes with alternative boundaries.

(2)

SWAT model is used in the study to discriminate each sub watershed and simulate P distribution in each watershed. Determining the location of BMP installation based on the results of SWAT model simulation is reasonable.

(3)

The ISFFP model in the study could be converted into different models according to the types of uncertainties.

(4)

According to the spatial layout scheme that corresponds to the upper and lower limit scenarios under different P pollution reduction targets, the upper and lower limit scenarios represent the limit of uncertainty impact on the scheme. The scheme is reasonable when it is in the interval. The total costs in the research results and practical terms are interval numbers.

(5)

The developed method provides the schemes that correspond to the upper and lower limit scenarios. The method can provide more reasonable schemes for decision makers under uncertain conditions.

(6)

Given that identifying the uncertainty distribution mode in practice is always difficult, the advantage of the method used in the study is that it does not need all variables as the uncertainty number nor does it need to set the average value for uncertainty numbers. Therefore, the method developed in this study can be used better in practical condition.