Optimizing Non-Point Source Pollution Management: Evaluating Cost-Effective Strategies in a Small Watershed within the Three Gorges Reservoir Area, China

: Non-point source (NPS) pollution poses a significant threat to the water environment, yet controlling it at the watershed scale remains a formidable challenge. Understanding the characteristics and drivers of nitrogen (N) and phosphorus (P) outputs at the watershed scale, along with identifying cost-effective best management practices (BMPs), is crucial for effective pollution control. In this study, we utilized the Wangjiaqiao watershed within the Three Gorges Reservoir Area (TGRA) as a case study to explore the characteristics of N and P load outputs and their dominant drivers by combining the SWAT model and a geographic detector. Based on our analysis of N and P loads within the watershed, we employed the entropy weight method to evaluate the reduction efficiency and cost-effectiveness of 64 BMP scenarios, encompassing seven measures (vegetative filter strips, parallel terraces, 10% fertilizer reduction, 30% fertilizer reduction, residue cover tillage, grass mulching, and returning farmland to forest) and their combinations. Our findings revealed the following:


Introduction
In recent years, the safety of the global river water environment has faced significant threats due to the excessive discharge of nutrients, including nitrogen (N) and phosphorus (P), into water bodies [1][2][3].Studies indicate that 30-50% of surface water sources are affected by non-point source (NPS) pollution originating from upstream watersheds [4,5], triggering a chain reaction of eutrophication, reduced dissolved oxygen, and ecosystem dysfunction [6,7].To counteract further deterioration, countries worldwide have implemented strategies aimed at reducing NPS loads from watersheds.[8].Among these strategies, best management practices (BMPs) are commonly employed interventions to enhance water quality.However, before implementing BMPs, several challenges must be addressed.Watershed complexity influences N and P input, transport, and retention processes, which are affected by various natural elements (e.g., topography, soil properties, climate, and hydrology) and anthropogenic factors (e.g., land use; fertilizer application) [9][10][11].There exists a complex response between these potential drivers and NPS pollution load outputs [12,13].Considering the variations in geographic conditions across regions, understanding the NPS pollution load output characteristics, examining spatial differentials and the drivers of NPS pollution load loss, and conducting simulations and assessments of watershed BMPs for cost-effectiveness are essential prerequisites for effective watershed management.Such endeavors bear significant implications for the subsequent planning and execution of successful watershed pollution control BMPs.
Previously, various modeling tools have been utilized to simulate and predict NPS pollution conditions in watersheds [14,15].Some studies [16,17] have relied on the export coefficient model to estimate NPS loads.However, this empirical model has shortcomings, as it overlooks the unique conditions of the study area (e.g., topography, soil, vegetation, and rainfall), leading to difficulties in ensuring prediction accuracy and an inability to explain pollutant transport and transformation processes [18].By contrast, mechanistic models, rooted in physically dynamic process mechanisms, offer advantages in capturing complex nonlinear interactions during nutrient transport transformations [19].Various dynamic process-based mechanistic models have been developed, including SWAT, HSPF, APEX, AnnAGNPS, etc. [20][21][22].Among them, the SWAT (Soil and Water Assessment Tool) model has gained widespread adoption due to its robust simulation algorithms for the hydrologic cycle, sediment transport, and nutrient transport, as well as its comprehensive database of agricultural management practices [23][24][25].However, the nonlinear interactions of nutrient transport and transformation processes in hydrologic models can be challenging to intuitively understand and may require other methods to account for N and P load drivers [21].Previous studies have often used correlation analysis [26], multiple linear regression [27], and redundancy analysis [28] to analyze factors affecting water quality changes based on sampling data from a few monitoring sites.However, these studies often overlook the spatial variability of pollutant loads, making it difficult to assess the impact of drivers from the spatial perspective of the watershed.By contrast, the geographic detector [29] offers improved accuracy in capturing the spatial distribution of pollutant loads by considering the characteristics of geospatial data.Moreover, Geodetector can qualitatively and quantitatively analyze the effects of single factors and detect the explanatory power of two-factor interactions [30].Thus, combining SWAT modeling and Geodetector can address watershed-scale climatic, geographic, and agroeconomic complexities while detecting the spatial heterogeneity of NPS loads and revealing the driving forces behind them [31].
Implementing BMPs in watersheds is an effective way to mitigate the transport of NPS pollutants from land to surface waters [32].These BMPs are commonly classified into three categories: structural (such as vegetative filter strips and parallel terraces), nonstructural (such as fertilizer reduction and no-till), and landscape management (such as returning farmland to forest) [33].While the efficacy of BMPs in reducing NPS loads in watersheds has garnered significant attention, previous studies have primarily focused on individual BMP effectiveness, lacking systematic research on the effectiveness of combined BMPs.Moreover, the effectiveness of BMPs often varies regionally [34,35], indicating that BMPs effective in one watershed may not yield the same results in another.Traditional field experiments are constrained by a limited monitoring capacity due to time and space constraints, rendering them challenging to apply in watershed decision-making [36].One of the strengths of the SWAT model is its versatility as an effective watershed planning tool, offering a flexible framework for modeling various BMPs [37].However, beyond focusing solely on the environmental benefits, a greater challenge lies in integrating costs and environmental benefits to better evaluate BMP effectiveness and target configuration options in specific watersheds [32].Indeed, most prior studies have primarily focused on environmental benefit analyses, with inadequate exploration of the trade-offs between pollution control and the cost investments of different BMP scenarios [38][39][40].Several studies [41] have explored using genetic algorithms coupled with SWAT models for environmental and economic multi-objective optimization strategies to derive BMPs for each hydrologic response unit (HRU).Despite the perceived sophistication, this approach remains in its infancy, demanding substantial computational resources and time for solutions.Additionally, it often overly emphasizes spatial heterogeneity, leading to overly complex outcomes challenging practical application [32].It is crucial to recognize that complexity in methods does not necessarily equate to greater efficiency in problem-solving.In terms of cost-benefit evaluation indices, the entropy weight method emerges as a suitable multi-attribute decision-making method based on information entropy theory, apt for multivariate comprehensive evaluation [42].To mitigate the subjective influence of arbitrarily determined indicator weights on cost-effectiveness evaluations, the entropy weight method objectively determines weights by fully assessing the actual situation of each indicator's data, thereby reflecting internal differences between the data [43].Wu et al. [42] concluded that combining the SWAT model with the entropy weight method serves as a suitable alternative analytical tool for assessing BMP effectiveness in watershed management.
The Three Gorges Reservoir Area (TGRA) serves as a critical water source for China's water allocation, yet it confronts significant ecological sensitivity, enduring persistent threats from NPS pollution [44].This study centers on the Wangjiaqiao watershed, a representative agricultural watershed within the TGRA.This study aims to (1) analyze the characteristics of N and P load outputs and investigate the primary drivers and interactions influencing N and P outputs within the watershed; (2) evaluate the NPS pollution reduction efficiency of BMPs both individually and in combination; and (3) conduct a comprehensive cost-benefit analysis aimed at identifying optimal implementation strategies tailored to the study area.The outcomes of this study contribute to a deeper understanding of NPS pollution characteristics within the watershed, providing valuable insights to decisionmakers for optimizing BMP strategies while considering cost-benefit trade-offs.

Study Area
The Wangjiaqiao watershed (31 1), with an area of 16.7 km 2 , is located in Zigui County in the eastern part of the Three Gorges Reservoir Area (TGRA) of China [45].Characterized by mountainous terrain, its elevation ranges from 201 to 1060 m.The watershed exhibits a typical subtropical monsoon climate, with an average annual temperature of 17.3 • C.Over the period from 2012 to 2022, the average rainfall stood at 864 mm, with 70% of the precipitation concentrated from May to September.According to the Chinese soil classification system, the major soil types in this area are purple soils, paddy soils, and yellow soils.The primary land use types are forest, orchard, cropland, and shrub forest.Vegetation predominantly consists of secondary broad-leaved mixed species, such as pine, cypress, and citrus.Cultivation within the area primarily focuses on wheat, rice, and corn.The watershed is sparsely populated, devoid of industrial point sources or grazing land distribution.Since the 1980s, a considerable portion of forested land has undergone conversion to citrus cultivation, making citrus cultivation the primary economic activity in the region [27].To enhance citrus production, a substantial amount of fertilizer is applied.primarily focuses on wheat, rice, and corn.The watershed is sparsely populated, devoid of industrial point sources or grazing land distribution.Since the 1980s, a considerable portion of forested land has undergone conversion to citrus cultivation, making citrus cultivation the primary economic activity in the region [27].To enhance citrus production, a substantial amount of fertilizer is applied.

SWAT Model Set Up
The input data (Table 1) required by the SWAT model include the digital elevation model (DEM), land use data (Figure S1, Table S1), soil data (Figure S1, Table S2), meteorological data, runoff and water quality data, and an agricultural management plan (Table S3).The land use status within the study area remained unchanged in 2012-2022; thus, the land use data for the intermediate year 2017 was selected.It was assumed that meteorological variables, including precipitation, air temperature, and solar radiation, exhibited homogeneity owing to the small scale of the watershed [46].In setting up the SWAT model crop database, considering the need for model simplification, wheat and corn rotations were primarily considered for sloping fields and terraces, with single-season rice designated for paddy fields and citrus selected for orchards.To closely replicate the water system, the study area was subdivided into 103 sub-watersheds using a threshold of 8 ha.Following the subdivision of the sub-watersheds, the thresholds for land use, soil, and slope were set at 0%, 0%, and 0%, respectively, to retain the full details required for generating 2109 hydrologic response units (HRUs).

SWAT Model Set Up
The input data (Table 1) required by the SWAT model include the digital elevation model (DEM), land use data (Figure S1, Table S1), soil data (Figure S1, Table S2), meteorological data, runoff and water quality data, and an agricultural management plan (Table S3).The land use status within the study area remained unchanged in 2012-2022; thus, the land use data for the intermediate year 2017 was selected.It was assumed that meteorological variables, including precipitation, air temperature, and solar radiation, exhibited homogeneity owing to the small scale of the watershed [46].In setting up the SWAT model crop database, considering the need for model simplification, wheat and corn rotations were primarily considered for sloping fields and terraces, with single-season rice designated for paddy fields and citrus selected for orchards.To closely replicate the water system, the study area was subdivided into 103 sub-watersheds using a threshold of 8 ha.Following the subdivision of the sub-watersheds, the thresholds for land use, soil, and slope were set at 0%, 0%, and 0%, respectively, to retain the full details required for generating 2109 hydrologic response units (HRUs).

SWAT Model Calibration and Validation
The year 2012 was designated as the warm-up period for the model to ensure its full operation within the hydrologic cycle during simulation.Calibration was executed in three phases: initially, flow calibration; subsequently, TN calibration; and finally, TP calibration.The research calibrated and validated runoff data over a decade, spanning 2013 to 2022, with monthly intervals.A 5-year calibration phase (2013-2017) was designated for tuning sensitive parameters, whereas a subsequent 5-year validation phase (2018-2022) was employed to evaluate the SWAT model's suitability within the study watershed.Given the limited availability of data, nutrient calibration (TN; TP) was based on the entire N and P dataset.The sensitivity analysis of the model parameters was conducted using the SUFI-2 (Sequential Uncertainty Fitting version 2) algorithm through SWAT-CUP.We selected 20 parameters for runoff and nutrient transport processes based on Ding et al. [47] and Lee et al. [35].Table 2 presents the sensitivity ranking and parameter values.To ensure the precision of model simulation, eight iterations were conducted, each comprising 500 runs.During each iteration, we gradually adjusted the parameter ranges until they converged to their optimal values (Table 2).To evaluate the performance of the SWAT model, we calculated the coefficient of determination (R 2 ), the Nash-Sutcliffe efficiency coefficient (NSE), and a Percent Bias (PBIAS) [48,49]: where n is the total number of data records; obs(i) and sim(i) are the measured and simulated values of the time step i; and obs and sim represent the average of the measured and simulated values, respectively.In this study, we followed the threshold values of the indicators recommended by Moriasi et al. [48] for evaluation, as shown in Table 3.

Geographic Detector
A geographic detector typically analyzes the effects of correlated factors by detecting spatial differentiation and identifying the main drivers [50].Referring to existing studies [18,21,51], we identified nine potential drivers (land use, soil, elevation, slope, fertilizer application, surface runoff, USLE_LS factor, USLE_K factor, and distance from the water system) and extracted sample points by constructing a fishing net to spatially obtain data on these factors.The factor detection module quantitatively identifies the contribution of each driver to N and P loads, while the interaction detection module evaluates the combined effect of the two drivers (X1 ∩ X2).A q-value was employed to quantify the magnitude of influence of a single factor or multiple interacting factors, expressed as follows: where h = 1, • • • , m is the stratification of a single factor, X, or the crossed strata of multiple factors, X; N h and N are the number of units in the stratum, h, and the entire area, respectively; and σ 2 h and σ 2 are the variances of the dependent variable in stratum h and the entire area, respectively.The q-value ranges from 0 to 1, indicating the degree of deterministic influence from a single factor or multiple interactions.The significance of the q-value was determined using an F test with a significance level of 0.05.

BMP Implementation Setup and Efficiency Evaluation
Considering that agricultural land is a major source of NPS pollution loads, this study selected commonly used BMPs for simulation [52], including structural measures: vegetative filter strips (VFSs) and parallel terraces (PTs); non-structural measures: fertilizer reduction (FR), residue cover tillage (RCT), and grass mulch (GM); and a landscape management measure: returning farmland to forest (RF).More details can be found in the SWAT input/output file documentation [52].Specific parameter adjustments are described in Table S4, along with the cost of each measure.As FR30 and FR10 are BMPs of the same type, and FR30 is more effective, after single-BMP simulation, the other six BMPs, except FR10, are simulated in scenarios with combination numbers ranging from 2 to 6, resulting in a total of 64 scenarios.Combinations of BMPs are denoted as follows: VFSs and PTs are combined as VFS_PT; VFSs, RCT, and FR30 are combined as VFS_RCT_FR30.The reduction rate, r, is calculated based on the BMP simulation results using the following equation: where LOAD pre-bmp and LOAD post-bmp denote the annual pollution loads at the watershed outlet before and after the implementation of BMPs, respectively.

Entropy Weight Method
The entropy method utilizes an objective scoring approach to reduce the subjective influence on the weights of various indicators [42].Considered indicators include economic costs (annual investment costs for each BMP scenario), TN reduction efficiency, and TP reduction efficiency.By calculating the comprehensive benefit value, z, for different BMP scenarios and ranking them, decision-makers can obtain an optimized watershed BMP implementation plan.The procedural steps for calculation are outlined below: (i) Constructing the BMPs' attribute decision matrix, A: where a ij represents the jth evaluation factor of the ith BMPs, and n and m are the number of BMPs and evaluation factors, respectively.
(ii) Decision matrix standardization: The decision matrix, A, was standardized using the range standardization method to mitigate the impact of physical magnitude on the overall evaluation outcomes.
For the benefit factor, the standardized formula is For the cost factor, the standardized formula is By standardizing the benefit and cost factors, a standardized decision matrix, R, consisting of n BMPs and m evaluation factors was obtained: Land 2024, 13, 742 8 of 21 (iii) Decision matrix normalization, R: (iv) Information entropy, E j , for each evaluation factor: where rij ln rij = 0 when rij = 0. (v) Weight vector, w, of each evaluation factor: (vi) Comprehensive evaluation index, z, for each BMP:

Model Performance in Simulating Flow and TN/TP Loads
Figure 2 displays the calibration and validation results for the flow and nutrient loads.According to monthly simulation results, the flow calibration period exhibited better performance compared with the validation period, with respective R 2 values of 0.85 and 0.70, NSE values of 0.85 and 0.7, and PBIAS values of 6.57% and −10.58%.Regarding monthly TN loads, the calibration results were satisfactory, with an R 2 of 0.94, an NSE of 0.91, and a PBIAS of −27.31%.Similarly, for TP loads, the calibration yielded favorable results, with an R 2 of 0.96, an NSE of 0.92, and a PBIAS of −27.64%.The evaluation criteria were suitably relaxed due to the higher uncertainties inherent in nutrient calibration, as recommended by Moriasi et al. [48].Overall, the model evaluation indices indicate that the SWAT model demonstrates satisfactory performance in simulating flow and nutrient loads in the study area, rendering it suitable for investigating the effects of BMPs.

Distribution Characteristics
Figure 3 shows the multi-year average N and P loading intensities in the watershed.The spatial distribution of TN and TP loads exhibits significant variation across different locations, demonstrating clear spatial heterogeneity.The maximum TN output intensity in the watershed reaches 56.57kg/ha, while the maximum TP loading intensity reaches 1.05 kg/ha.According to the statistical results presented in Table 4, the average annual output intensities of TN and TP in the entire watershed are 13.25 kg/ha and 0.11 kg/ha, respectively.The average annual outputs of TN and TP in the watershed export amount to 20,281.91 kg and 168.38 kg, respectively.Concerning volume and intensity, TN pollution surpasses TP, indicating a more severe pollution level.Moreover, in terms of pollutant forms, organic nitrogen loss predominates in TN, while mineral phosphorus loss dominates in TP.
Land 2024, 13, 742 9 of 21 0.91, and a PBIAS of −27.31%.Similarly, for TP loads, the calibration yielded favorable results, with an R 2 of 0.96, an NSE of 0.92, and a PBIAS of −27.64%.The evaluation criteria were suitably relaxed due to the higher uncertainties inherent in nutrient calibration, as recommended by Moriasi et al. [48].Overall, the model evaluation indices indicate that the SWAT model demonstrates satisfactory performance in simulating flow and nutrient loads in the study area, rendering it suitable for investigating the effects of BMPs.

Distribution Characteristics
Figure 3 shows the multi-year average N and P loading intensities in the watershed.The spatial distribution of TN and TP loads exhibits significant variation across different locations, demonstrating clear spatial heterogeneity.The maximum TN output intensity in the watershed reaches 56.57kg/ha, while the maximum TP loading intensity reaches 1.05 kg/ha.According to the statistical results presented in Table 4, the average annual  4), notable variations in nutrient loading intensities are evident for each type.The sloping field exhibited the highest TN loading intensity, with an average loading intensity of 39.06 kg/ha, which was 2.95 times higher than the watershed's mean value.Conversely, the land use type with the highest TP loading intensity was the paddy field, displaying a loading intensity 3.09 times higher than the watershed's mean value.In terms of the proportion of total load (Figure 4), Orchards demonstrated the highest percentages of TN and TP outputs, both exceeding 50% of the overall watershed loss.Forest TN and TP outputs accounted for 12.53% and 13.7%, respectively, despite covering the largest proportion of the watershed area at 44.66%.The residential area, constituting 3.08% of the watershed area, contributed 6.75% and 2.03% of the TN and TP output, respectively.Agricultural lands, including orchards, sloping fields, terraces, and paddy fields, accounted for 42.3% of the watershed area.However, they contributed 76.62% and 78.31% of the TN and TP output, respectively.These results indicate that agricultural lands significantly contribute to TN and TP loads.When comparing the nutrient loads across different land use types (Table 4), notable variations in nutrient loading intensities are evident for each type.The sloping field exhibited the highest TN loading intensity, with an average loading intensity of 39.06 kg/ha, which was 2.95 times higher than the watershed's mean value.Conversely, the land use type with the highest TP loading intensity was the paddy field, displaying a loading intensity 3.09 times higher than the watershed's mean value.In terms of the proportion of total load (Figure 4), Orchards demonstrated the highest percentages of TN and TP outputs, both exceeding 50% of the overall watershed loss.Forest TN and TP outputs accounted for 12.53% and 13.7%, respectively, despite covering the largest proportion of the watershed area at 44.66%.The residential area, constituting 3.08% of the watershed area, contributed 6.75% and 2.03% of the TN and TP output, respectively.Agricultural lands, including orchards, sloping fields, terraces, and paddy fields, accounted for 42.3% of the watershed area.However, they contributed 76.62% and 78.31% of the TN and TP output, respectively.These results indicate that agricultural lands significantly contribute to TN and TP loads.

Driving Factors
The results of the factor detector are depicted in Figure 5.In determining the spatial distribution of TN, land use emerged as the most dominant driver (0.626), followed by nitrogen fertilizer application (0.502) and surface runoff (0.395).The relative importance of the other six influences ranged from 0 to 0.101.For TP, the top three dominant factors were surface runoff (0.71), land use (0.28), and phosphorus fertilizer application (0.258).
The findings of the interaction detection (Figure 6) reveal that the explanatory power of the interaction between land use and slope, soil type, surface runoff, USLE_LS, and USLE_K for TN loads exceeded 0.74.Notably, this explanatory power was significantly stronger compared with when these factors acted alone.Moreover, the explanatory power of surface runoff ∩ fertilizer application was found to be 0.743.As for TP, the q-values of surface runoff interacting with other factors ranged from 0.722 to 0.886, with the highest

Driving Factors
The results of the factor detector are depicted in Figure 5.In determining the spatial distribution of TN, land use emerged as the most dominant driver (0.626), followed by nitrogen fertilizer application (0.502) and surface runoff (0.395).The relative importance of the other six influences ranged from 0 to 0.101.For TP, the top three dominant factors were surface runoff (0.71), land use (0.28), and phosphorus fertilizer application (0.258).

Efficiencies of BMP Scenarios
Figure 7a,b illustrate the reduction effects of individual BMPs over the simulation period.Table 5 presents the average annual reduction rates.The TN reduction rates for each measure in descending order are PTs (28.12%),VFSs (18.75%),FR30 (16.34%),RCT (7.06%), RF (5.81%),GM (5.79%), and FR10 (5.38%).As for the TP reductions, they are PTs (37.69%),VFSs (32.53%),GM (15.81%),FR30 (15.77%),RCT (13.03%),FR10 (5.35%), and RF (5.21%), in descending order.Among the individual BMPs, parallel terraces proved to be the most effective, followed by VFSs, both of which are structural BMPs.The reduction rates of PTs and VFSs exhibit high variability.Specifically, TN and TP reduction rates in PTs range from 4.17% to 36.79% and 23.61% to 48.58%, respectively.TN and TP reduction rates in VFSs range from 7.79% to 22.63% and 22.27% to 50.39%, respectively.Conversely, the reduction efficiency of FR10 and RF was lower, with an average reduction rate of less than 6% for both TN and TP.All BMPs, except for FR and RCT, were more effective in reducing TP than TN.
Figure 7c,d depict the variation in N and P reduction efficiency for different combinations.For TN loads, the average reduction rates for individual BMPs and combined BMPs with 2, 3, 4, and 5 measures were 12.46%, 24.91%, 34.61%, 42.15%, and 48.47%, respectively.The measures with the largest reduction rates in each group were PTs (28.12%),PT_FR30 (39.21%),VFS_PT_FR30 (48.05%),VFS_PT_RF_FR30 (49.3%), and VFS_PT_RF_RCT_FR30 (53.15%).For TP load, the mean reduction rates for the individual The findings of the interaction detection (Figure 6) reveal that the explanatory power of the interaction between land use and slope, soil type, surface runoff, USLE_LS, and USLE_K for TN loads exceeded 0.74.Notably, this explanatory power was significantly stronger compared with when these factors acted alone.Moreover, the explanatory power of surface runoff ∩ fertilizer application was found to be 0.743.As for TP, the q-values of surface runoff interacting with other factors ranged from 0.722 to 0.886, with the highest being the surface runoff ∩ fertilizer application.In summary, the loss of N and P loads in the watershed is influenced significantly by land use type and fertilizer application resulting from human activities, as well as surface runoff among natural factors.Conversely, the effects of six factors, namely, USLE_K, soil type, slope, USLE_LS, elevation, and distance from the water system, were observed to be relatively weak.
crease.However, there are also scenarios where BMPs with fewer measure combinations demonstrate reduction rates that approximate or even surpass the reduction effects achieved by BMPs with a higher number of combinations.For instance, PT_RCT demonstrated higher rates of TN and TP reduction, reaching 34.85% and 49.24%, respectively, compared with VFS_RCT_FR30, which achieved reductions of 31.32% for TN and 48.08% for TP.The combined implementation of BMPs proves to be less effective than simply summing up the pollutant reduction rates of individual BMPs.For instance, in the case of PT_FR30, the combination resulted in a reduction rate of 39.21%.However, the combined reduction rates for PTs and FR30 were 5.25% higher than that of PT_FR30.As the number of measure combinations increases, range of reduction rate fluctuation gradually decreases and stabilizes.

Cost-Effectiveness of BMP Scenarios
Table 5 presents the cost-effectiveness of different individual BMPs.In the simulation, FR30 is the most effective BMP with a rating of 0.51.Compared with other measures, FR30 offers cost-saving benefits.Additionally, VFSs and PTs also boast high ratings of 0.47 and 0.45, respectively-the benefit values of FR10, GM, RCT, and RF range from 0.24 to 0.31.The data in Figure 8 illustrate the distribution of all simulated scenarios.Among the combinations of BMPs with 2, 3, 4, and 5 measures, the BMPs with the highest combined environmental cost-benefit values are PT_FR30 (0.68), VFS_PT_FR30 (0.8), VFS_PT_RCT_FR30 (0.81), and VFS_PT_RF_RCT_FR30 (0.76).When all six individual BMPs are combined simultaneously, the benefit value is 0.69.In general, the evaluated values of these combined measures are significantly greater than FR30, indicating their superior cost-effectiveness in reducing NPS pollution in the watershed.
The top ten combinations, with the highest combined benefit values, are enclosed in the red dashed circle area of Figure 8, and their corresponding results are detailed in Table  Figure 7c,d depict the variation in N and P reduction efficiency for different combinations.For TN loads, the average reduction rates for individual BMPs and combined BMPs with 2, 3, 4, and 5 measures were 12.46%, 24.91%, 34.61%, 42.15%, and 48.47%, respectively.The measures with the largest reduction rates in each group were PTs (28.12%),PT_FR30 (39.21%),VFS_PT_FR30 (48.05%),VFS_PT_RF_FR30 (49.3%), and VFS_PT_RF_RCT_FR30 (53.15%).For TP load, the mean reduction rates for the individual BMPs and the combined BMPs with 2, 3, 4, and 5 measures were 17.91%, 35.35%, 47.88%, 57.18%, and 64.04%, respectively.The measures with the highest reduction rates in each group were PTs (37.69%),VFS_PT (57.84%),VFS_PT_RCT (65.57%),VFS_PT_RCT_FR30 (68.58%), and VFS_PT_RCT_GM_FR30 (69.59%).Upon combining all six individual BMPs, the TN reduction rate was 53.6%, and the TP reduction rate was 69%.As the number of combinations increases, the average reduction rates of TN and TP gradually increase.However, there are also scenarios where BMPs with fewer measure combinations demonstrate reduction rates that approximate or even surpass the reduction effects achieved by BMPs with a higher number of combinations.For instance, PT_RCT demonstrated higher rates of TN and TP reduction, reaching 34.85% and 49.24%, respectively, compared with VFS_RCT_FR30, which achieved reductions of 31.32% for TN and 48.08% for TP.The combined implementation of BMPs proves to be less effective than simply summing up the pollutant reduction rates of individual BMPs.For instance, in the case of PT_FR30, the combination resulted in a reduction rate of 39.21%.However, the combined reduction rates for PTs and FR30 were 5.25% higher than that of PT_FR30.As the number of measure combinations increases, the range of reduction rate fluctuation gradually decreases and stabilizes.

Cost-Effectiveness of BMP Scenarios
Table 5 presents the cost-effectiveness of different individual BMPs.In the simulation, FR30 is the most effective BMP with a rating of 0.51.Compared with other measures, FR30 offers cost-saving benefits.Additionally, VFSs and PTs also boast high ratings of 0.47 and respectively-the benefit values of FR10, GM, RCT, and RF range from 0.24 to 0.31.The data in Figure 8 illustrate the distribution of all simulated scenarios.Among the combinations of BMPs with 2, 3, 4, and 5 measures, the BMPs with the highest combined environmental cost-benefit values are PT_FR30 (0.68), VFS_PT_FR30 (0.8), VFS_PT_RCT_FR30 (0.81), and VFS_PT_RF_RCT_FR30 (0.76).When all six individual BMPs are combined simultaneously, the benefit value is 0.69.In general, the evaluated values of these combined measures are significantly greater than FR30, indicating their superior cost-effectiveness in reducing NPS pollution in the watershed.
Land 2024, 13, x FOR PEER REVIEW 14 of 22 6. Their comprehensive benefit values range from 0.7 to 0.81, with cost fluctuations spanning −114,300 to 421,600 CNY.The TN reduction rates range from 35.11% to 53.15%, and the TP reduction rates range from 49.51% to 69.59%.Notably, the scenario with the lowest cost ranks ninth, while the scenario with the highest TN reduction rate ranks third, and the scenario with the highest TP reduction rate ranks fifth.The most effective combination of measures is VFS_PT_RCT_FR30, which has a combined benefit value of 0.81 and a cost input of only 187,900 CNY, resulting in a 52.37% reduction in TN and a 68.58% reduction in TP.The combination ranked second is VFS_PT_FR30, which has a combined benefit value of 0.8 and reduces TN and TP by 48.05% and 61.95%, respectively, with a cost input of 25,800 CNY.VFS_PT_RCT_FR30 and VFS_PT_FR30 provide similar benefits.However, the TN and TP reduction rates of VFS_PT_RCT_FR30 are slightly higher than those of VFS_PT_FR30 (approximately 1.08 and 1.11 times, respectively).Despite this, the cost of VFS_PT_RCT_FR30 is 7.3 times higher than that of VFS_PT_FR30.The 9th and 10th ranked combinations showed negative cost values but were less effective in reducing N and P than the other combinations.They could be considered for use when funds are scarce.The top ten ranked combinations have measures appearing in the following order of frequency: FR30, PTs, VFSs, RCT, GM, and RF.This further illustrates the value of FR30, PTs, and VFSs for local applications.The top ten combinations, with the highest combined benefit values, are enclosed in the red dashed circle area of Figure 8, and their corresponding results are detailed in Table 6.Their comprehensive benefit values range from 0.7 to 0.81, with cost fluctuations spanning −114,300 to 421,600 CNY.The TN reduction rates range from 35.11% to 53.15%, and the TP reduction rates range from 49.51% to 69.59%.Notably, the scenario with the lowest cost ranks ninth, while the scenario with the highest TN reduction rate ranks third, and the scenario with the highest TP reduction rate ranks fifth.The most effective combination of measures is VFS_PT_RCT_FR30, which has a combined benefit value of 0.81 and a cost input of only 187,900 CNY, resulting in a 52.37% reduction in TN and a 68.58% reduction in TP.The combination ranked second is VFS_PT_FR30, which has a combined benefit value of 0.8 and reduces TN and TP by 48.05% and 61.95%, respectively, with a cost input of 25,800 CNY.VFS_PT_RCT_FR30 and VFS_PT_FR30 provide similar benefits.However, the TN and TP reduction rates of VFS_PT_RCT_FR30 are slightly higher than those of VFS_PT_FR30 (approximately 1.08 and 1.11 times, respectively).Despite this, the cost of VFS_PT_RCT_FR30 is 7.3 times higher than that of VFS_PT_FR30.The 9th and 10th ranked combinations showed negative cost values but were less effective in reducing N and P than the other combinations.They could be considered for use when funds are scarce.The top ten ranked combinations have measures appearing in the following order of frequency: FR30, PTs, VFSs, RCT, GM, and RF.This further illustrates the value of FR30, PTs, and VFSs for local applications.The accurate comprehension of the output characteristics and primary drivers of N and P loads within watersheds is essential for the scientific development of pollution control strategies.The transport of N and P is influenced by various natural and anthropogenic factors [9], such as topography, precipitation, runoff, soils, land use, and fertilizer application.By employing the SWAT model in combination with Geodetector analysis, this study identified land use, fertilizer application, and surface runoff as the main factors driving spatial variations in N and P loads.The analysis revealed that land use predominantly dictates N distribution, whereas P distribution is primarily influenced by surface runoff due to its lower solubility compared with N. Furthermore, this study unveiled significant interaction effects among various factors.Notably, the combined effects of the runoff factor intersecting with land use and fertilizer application exhibited pronounced enhancement, amplifying their explanatory capacity.The study also revealed the significant enhancement effects of other factor interactions, indicating that, in actual processes, complex interactions between multiple factors collectively determine the spatial patterns of NPS pollution in watersheds [53].
The sources and transport processes of NPS pollution are highly dependent on hydrologic processes, with precipitation and runoff being the driving conditions for N and P transport [18,54,55].Precipitation and surface runoff processes generally occur simultaneously, and the greater the precipitation and surface runoff, the more favorable the transport of sediment and nutrients [45].Meteorological factors, such as precipitation, are assumed to be homogeneous within small watersheds; thus, the spatial exploration of their effects is neglected.However, runoff factors exhibit variability due to disparities in underlying surface characteristics, including slope, slope length, land use, soil properties, and other pertinent factors [21].The underlying surface serves as the interface for nutrient transport, with surface runoff's erosive and scouring forces escalating with increased slope and slope length, consequently amplifying nutrient transport capacity [56].Research has shown significant differences in NPS loading between various soil types [57].Soil types possess unique hydrological, biological, and background content traits, profoundly influencing soil erosion resistance, infiltration capacity, and flow production capacity [58].Land use impacts the generation and transport of N and P loads, which is related to factors affecting vegetation cover and agricultural activities [59][60][61].Wang et al. [19] found that the proportion of agricultural land use significantly affects nutrient loading.Vegetation cover affects the underlying surface characteristics, hydrological processes, and nutrient cycling, further influencing the storage and retention capacity of nutrients in different land use types [62].Woodlands and shrubs can absorb and sequester nutrients to some extent, which helps mitigate the effects of pollution due to higher vegetation cover and less anthropogenic disturbance [7,63].In contrast, agricultural land is plowed year-round, often exhibiting low vegetation cover or even bare ground.Crop growth relies on regular inputs of N and P fertilizers.Fertilizer application directly alters the initial soil nutrient content across various watershed regions.Elevated soil nutrient content increases the potential risk of agricultural NPS pollution formation through surface runoff and leaching loss [64].Tan et al. [18] used the export coefficient model and Geodetector to assert that human activities have a dominant role in N and P output compared with natural factors.This statement implies that changes in regional N and P loads are strongly influenced by anthropogenic factors.
It is worth noting that the weak explanatory power of elevation observed in this study may be attributed to the influence of spatial scale effects on identifying influential factors [65,66].Shen et al. [31] reported elevation as the predominant driver of spatial variance in the nutrient budget within the TGRA.The TGRA is 3832 times larger than the study watershed, with more drastic changes in elevation and a greater concentration of human activities at larger scales.Higher elevations typically exhibit greater vegetation cover and reduced human intervention, while lower elevations, particularly along rivers, are often urbanized or dedicated to agriculture [67].In the studied watershed, despite the substantial range in elevation (201-1060 m), the impact of human activities does not vary significantly across different elevation bands.
In conclusion, exploring the potential drivers of N and P loads offers valuable insights into the characteristics of NPS load output within the watershed.This exploration can also aid in the selection of BMPs and delineation of implementation areas.The main objective is to reduce N and P loss from agricultural land, specifically orchards and cropland, within the watershed.This can be achieved by implementing changes in land use, reducing N and P inputs, controlling topography, and mitigating runoff.Additionally, a current limitation is that the nine drivers selected for this study may not fully account for the variations in N and P load outputs.Including additional potential drivers in future studies could enhance the explanatory power of N and P load outputs.

Analysis of Reduction Efficiency of Watershed BMPs
In different study areas, the effectiveness of identical measures can vary due to differences in climatic, geographical, and agro-economic conditions [33].Leh et al. [34] observed a TN reduction of 1.2% to 4.9% from vegetated filter strips in an area with only 6% agricultural land.In contrast, Teshager et al. [68] modeled filter strips in a watershed with 78% agricultural land area and achieved reductions of approximately 50%.In our watershed, where agriculture covers 42.3% of the land area, filter strips can reduce TN by 7.79% to 22.63%.This highlights the significance of assessing and applying measures locally.
This study highlighted the notable reduction effects of structural measures such as terraces and vegetative filter strips, aligning with previous research by Sun et al. [32] and López-Ballesteros et al. [33].However, even the most effective single measure, terracing, achieved a maximum reduction rate of only 28.12% for TN and 37.69% for TP.This underscores the challenge of attaining high reduction rates with individual measures alone, emphasizing the importance of integrating different measures [69][70][71].Indeed, the study revealed that N and P reductions increased gradually with the combination of measures.It is essential to recognize that, while combined BMPs yielded better reduction effects than single measures, their impact was not purely additive due to interactions between the measures, as noted by Liu et al. [72].When employing combined BMPs, pollutant concentrations diminish from upstream to downstream post-treatment with upstream BMPs.Consequently, downstream reduction efficiencies may be lower compared with their standalone implementation [73].For instance, when reducing fertilizer use, the pollutant concentration in the vegetated filter strip treatment decreases, which, in turn, reduces the effectiveness of the filter strips in reducing pollutants [35].
Our study revealed that the reduction efficiency of TP caused by the implemented BMPs generally surpassed that of TN.This discrepancy can be attributed to the distinctive characteristics of elemental N and P in terms of form and transport transformation within water bodies, as well as the control mechanism of pollutants based on measures such as terracing, vegetative filter strips, and conservation tillage.Unlike P, N tends to be more soluble, with a portion lost through surface runoff, while another portion may leach into deeper soil layers, eventually reaching surface water through lateral flow or groundwater [74].The mechanism of action of BMPs primarily relies on the physical interception or deposition of sediment, whereas treating dissolved matter in flow often requires biological and chemical reaction processes, which are less efficient and take longer to treat [75].Consequently, controlling nitrogen proves to be more challenging than phosphorus due to these factors.

Analysis of Cost-Effectiveness of Watershed BMPs
To optimize BMPs, it is necessary to balance two interrelated and conflicting objectives: reducing pollution rates and minimizing input costs [76].Since the Grain for Green project is an environmental measure strongly advocated by China's ecological restoration project [77], our study simulated reforestation in line with national policy.However, considering the cost and pollution control perspective, the overall benefit of reforestation is not significant.Wu et al. [42] investigated the effect of reforestation on soil erosion control and found that it reached 32.16%.Nonetheless, the implementation cost was considered prohibitively high, resulting in poor overall benefits.Similarly, Sun et al. [32] excluded returning farmland to forest due to its high cost.Furthermore, this study shows that-due to the limited area of cultivable land meeting the requirements, which is only 40.77 ha, and the reduction rates of N and P being less than 6%-it is justifiable to rank the comprehensive benefit value at the lower end of the list.This study's comprehensive assessment of cost-effectiveness suggests that a rational mix of measures is necessary to achieve higher cost-effectiveness.It observes that combinations of BMPs with higher reduction rates often entail greater expenses, reflecting the inherent trade-off between pollution reduction rates and input costs [76].The top-ranked solutions identified in this study offer valuable insights and references for decision-making amid varying economic and environmental constraints [78].Decision-makers can select appropriate cost-effective options by making trade-offs based on specific economic costs and water quality objectives.Among the top ten combinations, vegetative filter strips, parallel terraces, and a 30% reduction in fertilizer were the most frequently recommended.The combination of all three BMPs ranked second.Vegetative filter strips are an economically inexpensive and effective method of preventing soil erosion and enhancing agricultural productivity.They have been recommended for implementation in many studies [34,35,79].Given the mountainous terrain of the study area, the transformation of sloping land holds particular significance [80].Moreover, FR30 emerged as the highest-rated individual BMP, offering both cost savings and substantial reductions in nutrient inputs from the source.Research indicates that agricultural fertilizer application in the TGRA is generally excessive [81,82].Wang et al. [82] emphasized the importance of reducing N fertilizer input to mitigate net anthropogenic N input.Consequently, these findings suggest that combining the three measures above could be a suitable approach for implementation in the study area, as well as other watersheds within the TGRA.

Conclusions
In this study, we established a SWAT model for the Wangjiaqiao watershed to analyze the characteristics of nutrient load outputs.Utilizing a geographic detector, we identified the principal driving factors impacting the spatial distribution of nutrient loads within the watershed.Furthermore, we conducted simulations of various BMP scenarios and explored the feasibility of the watershed cost-benefit assessment framework through the entropy weight method.Our findings revealed spatial heterogeneity in N and P loads, with land use, fertilizer application, and surface runoff as primary drivers, exhibiting interactive enhancement effects.Structural BMPs demonstrated superior N and P reduction efficiency compared with individual measures.Combined BMPs showed greater reduction efficiency, but some combinations yielded less than the sum of individual measures.SWAT combined with the entropy weight method proved effective for selecting BMP schemes.The comprehensive benefits of combined BMPs were generally higher, highlighting the need for diverse strategies to address NPS pollution.The combination of vegetative filter strips, parallel terraces, and a 30% reduction in fertilizer application emerged as a promising approach, offering a balance between reduction benefits and economic feasibility, particularly for small watersheds within TGRA.
In conclusion, this study deepens our understanding of NPS pollution outputs and their underlying driving factors, providing valuable insights decision-makers can use to evaluate and optimize BMP strategies.In light of the limited available data, we exerted the utmost effort to calibrate the runoff and N/P load.Nevertheless, uncertainties stemming from model input data, structure, and parameters persist, indicating the necessity for further research in this domain.The predominant pollution sources in the study watershed are agricultural non-point sources, whereas contributions from residential activities and livestock farming are relatively minor and insignificantly impact these research findings; thus, they were appropriately disregarded.Due to the localized nature of optimization schemes and regional variations in BMP cost-effectiveness, future research should validate the reproducibility of our framework in basins with different scales and characteristics.Moreover, refining the BMP design around NPS pollution drivers, expanding the BMP application database, and developing corresponding simulation methods are crucial steps forward.Furthermore, there is a need to strengthen strategic research on managing stormwater runoff, improving land use structures, and enhancing agricultural production management at the basin scale to address the governance of NPS pollution issues.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/land13060742/s1, Figure S1.Land use map (a) and soil map (b).Table S1.Description of the land use types and codes.Table S2.Soil properties in the study watershed.Table S3.Agricultural practices schedule in the study watershed.Table S4.Individual BMP implementation scenario setting.
Author Contributions: R.C.: methodology, software, data curation, writing-original draft.Y.W.: supervision, conceptualization, project administration.H.L.: conceptualization, methodology, data curation.Z.W.: writing-review and editing, software, data curation.L.M.: data curation, investigation.J.Z.: writing-review and editing, investigation.J.L.: writing-review and editing, investigation.Z.Y.: methodology.Y.Z.: methodology.D.L.: investigation.All authors provided feedback on the analyses and contributed to writing the manuscript.All authors have read and agreed to the published version of the manuscript.

Figure 1 .
Figure 1.Location of (a) the Three Gorges Reservoir Area, (b) Zigui County, and (c) the Wangjiaqiao watershed.

Figure 1 .
Figure 1.Location of (a) the Three Gorges Reservoir Area, (b) Zigui County, and (c) the Wangjiaqiao watershed.

Figure 2 .
Figure 2. Observed and simulated monthly (a) streamflow, (b) TN, and (c) TP loads at the outlet.TN/TP has only a calibration period due to a lack of sufficient data for validation.

Figure 2 .
Figure 2. Observed and simulated monthly (a) streamflow, (b) TN, and (c) TP loads at the outlet.TN/TP has only a calibration period due to a lack of sufficient data for validation.

Figure 3 .
Figure 3. Spatial distribution of multi-year average (a) TN load and (b) TP load at HRU level.

Figure 3 .
Figure 3. Spatial distribution of multi-year average (a) TN load and (b) TP load at HRU level.

Figure 4 .
Figure 4.The proportion of area, TN load, and TP load to the total in different land use types.

Figure 4 .
Figure 4.The proportion of area, TN load, and TP load to the total in different land use types.

Figure 5 .
Figure 5.The q-values of potential drivers of (a) TN and (b) TP.

Figure 6 .
Figure 6.The q-values for the interaction of potential drivers of (a) TN and (b) TP.

Figure 5 .
Figure 5.The q-values of potential drivers of (a) TN and (b) TP.

Land 2024 , 22 Figure 5 .
Figure 5.The q-values of potential drivers of (a) TN and (b) TP.

Figure 6 .
Figure 6.The q-values for the interaction of potential drivers of (a) TN and (b) TP.

Figure 6 .
Figure 6.The q-values for the interaction of potential drivers of (a) TN and (b) TP.

Figure 7 .
Figure 7. (a) TN and (b) TP reduction efficiency of individual BMPs; (c) TN and (d) TP reduction efficiency of BMPs with a different number of combinations.

Figure 7 .
Figure 7. (a) TN and (b) TP reduction efficiency of individual BMPs; (c) TN and (d) TP reduction efficiency of BMPs with a different number of combinations.

Figure 8 .
Figure 8. TN/TP reduction rates, costs, and comprehensive evaluation indices for all BMPs.

Figure 8 .
Figure 8. TN/TP reduction rates, costs, and comprehensive evaluation indices for all BMPs.

Table 1 .
Data used to set up and evaluate the SWAT model in the Wangjiaqiao watershed, China.

Table 2 . The sensitivity and value of parameters included in the calibration procedure. Variable Parameter Description Range Method Rank Value
Note: "r" means the initial value of the parameter multiplied by (1 + calibration value); "v" means the initial value of the parameter is replaced by the calibration value.

Table 3 .
Performance criteria for recommended statistics (R 2 , NSE, and PBIAS) for a monthly time step.

Table 4 .
Summary of pollution loads on different land use types, kg/ha.

Table 4 .
Summary of pollution loads on different land use types, kg/ha.

Table 5 .
The reduction efficiency, cost, and comprehensive benefit value of individual BMPs.

Table 5 .
The reduction efficiency, cost, and comprehensive benefit value of individual BMPs.

Table 6 .
The top 10 BMPs in comprehensive benefits of environmental costs ranking.