Evaluation of Best Management Practices for non-point source pollution based on the SWAT model in the Hanjiang River Basin, China

As a distributed hydrological model of the river basin, the SWAT model has gradually become an indispensable and powerful tool in water resources and water environmental protection. It is usually used to evaluate the long-term effects of climate change patterns and land management patterns in complex basins on runoff, sediment and agricultural nutrients (Ding 2018; Himanshu et al. 2019). For example, Sharma & Tiwari (2019) estimated that the monthly average sediment yield, total nitrogen (TN) load and total phosphorus (TP) load of Mayitong Reservoir were 1.53 million tons, 1,834.2 kg and 191.1 kg, respectively. With the application of SWAT model more and more widely, the coupling model is constantly appearing. LV et al. (2020) believed that MLP/MWOA algorithm can significantly improve the accuracy of runoff prediction, and that the SWAT-MLP/MWOA model is the best model for runoff prediction. Liu (2016) coupled Long-Term Hydrologic Impact Analysis (L-THIA) and SWAT model to simulate the non-point source (NPS) pollution load in the Xiangtan River Basin, which better reflected the transport process of nitrogen and phosphorus.

The idea of Best Management Practices (BMPs) was first proposed by the United States Department of Agriculture (USDA) in the mid-1970s to reduce or prevent water pollution. At present, the commonly used management practices mainly focus on engineering management practices, such as the use of riparian vegetation ecosystem and retention ponds to reduce or eliminate the impact and harm on agricultural production (Meng et al. 2013). The study by Liu et al. (2019) in northeastern Indiana showed that a vegetation buffer zone has the most obvious effect on reducing DRP/TP at annual scales and in spring. Sun et al. (2020) found that in the Loess Plateau region, the contribution rate of landscape engineering management practices to runoff reduction and sediment reduction was the largest in the past 20 years. The contribution rate to runoff reduction increased from 54.4% of P₁ (1997–2006) to 94.7% of P₂ (2007–2016), but the contribution rate to sediment reduction had no obvious change (61.6% of P₁ to 62.8% of P₂). Qiu et al. (2020) believed that the increase in the frequency of rainfall events may reduce the efficiency of BMPs in management of Miyun reservoir, and the watershed management should be adjusted according to climate change in the future. Zhang et al. (2013) proposed the idea of ‘risk assessment-planning zoning-separate management’ to evaluate the NPS pollution in the watershed step by step, which provided a scientific plan for the management of NPS pollution sources in the basin.

At present, in the Shaanxi section of the Hanjiang River Basin, there are studies on water environment carrying capacity, analysis of characteristics of vegetation cover change, spatial and temporal variation and source analysis of NPS pollution, CSA identification of NPS pollution, water pollution planning and prevention, ecological compensation, etc. Qi et al. (2020) focused on the dynamic precipitation pattern and the potential runoff production capacity of the region on annual and seasonal time scales in the upper Hanjiang River Basin from 1961 to 2016. Shi et al. (2011) investigated the current situation and pollution sources of water pollution in the Hanjiang River Basin, and found that the main pollutants were chemical oxygen demand (COD), TP and ammonia nitrogen (NH₃⁻N), the main pollution sources were agricultural surface runoff, urban domestic sewage, livestock and poultry breeding. Tang et al. (2018) analyzed the source of agricultural NPS pollution in the Hanjiang River Basin, which showed that the high load areas of equivalent standard pollution were concentrated in the middle of the river basin and based on the rapid clustering method, it was determined that there are six pollution types.

As the largest tributary of the Yangtze River, the Hanjiang River's upstream agricultural NPS pollution, mining, industrial waste water, urban domestic sewage, air pollution, surface runoff, and soil erosion may cause sudden changes in the water quality. Therefore, this paper takes the Hanjiang River Basin controlled by Ankang hydrological station as the study area, analyzes the temporal and spatial distribution characteristics of NPS pollution based on SWAT model, identifies the CSA, and then conducts reduction effects of management practices, puts forward the BMPs for the districts and counties, and provides a certain basis for the treatment of NPS pollution in the Hanjiang River Basin.

Hanjiang River is the biggest tributary in the middle reach of Yangtze River. It originates from the southern foot of the Qinling Mountains in Ningqiang County, Shaanxi Province, northwest of China. The Hanjiang River Basin covers an area of 159,000 km², with a total length of 1,577 km, among which Shaanxi Province is 652 km, accounting for 41.3%. The study area of this paper is the Hanjiang River Basin with the Ankang Hydrological Station as the control section with a drainage area of 35,347 km², as shown in Figure 1.

Influenced by the climate of Qinling Mountains and Daba Mountains, the Hanjiang River Basin belongs to the tropical monsoon climate zone. The area has abundant rainfall, with an annual rainfall of 750–1,100 mm. In winter, the thermal properties of the land and sea are different, the wind blows from the land to the sea, and the precipitation is less. In summer, due to the deflection of the air pressure belt, the wind blows from the ocean, and the precipitation is more. Therefore, the rainfall in the year is mainly concentrated in June to September. In addition, under the influence of vertical regionalism, the annual average temperature is 12–15 °C, the sunshine duration is 1,400–1,900 h, and the frost-free period is 210–270 d. The main soil types are: brown soil, coarse bone soil, yellow brown soil, paddy soil, brown soil, etc. Yellowish brown soil is mainly distributed in low mountains and hilly areas, brown soil is mainly distributed in the high mountains of Qinling Mountains, and paddy soil is distributed in local areas and along the Hanjiang River.

The SWAT model database (Table 1) mainly includes two categories: spatial data and attribute data. Spatial data include digital elevation model (DEM), land use map and soil type map. The DEM is downloaded from the Geospatial Data Cloud (GS Cloud) website with 30 m × 30 m resolution (Figure 1). The land use map from 2017 was from the Global Ecological Environment Remote Sensing Monitoring (GEERSM) platform with 30 m × 30 m resolution (Figure 2). The soil data were from the Harmonized World Soil Database (HWSD), and the grid size was 1 km × 1 km (Figure 3). Attribute data include meteorological data, soil physical and chemical properties, water quality data and agricultural management data. The description and sources of data related to this study are shown in Table 1.

During the construction of the SWAT model, the land use index table was established according to the land use database provided by the model to reclassify the land use (Table 2). At the same time, the index table of soil type was established and the raster data of soil type was reclassified. Since management practices are laid out within a certain slope range in subsequent studies, the slope was divided into five grades (≤10°, 10°–15°, 15°–25°, >25°). Spatial datasets (DEM, land use map, soil type map) required for the SWAT model setup were projected into the same co-ordinate system (The WGS_1984_UTM_Zone_48N) under the ArcGIS platform.

In this study, the study area was divided into 31 sub-basins, based on the user-defined threshold area of 71,436.19 hm² along with outlet points (Ankang Hydrological Station), as shown in Figure 4.

The SWAT model further divided the sub-basins into hydrological response units (HRU), which have the same land use, soil type, and slope. This paper assumes that all land use, soil types and slopes in the sub-basins are reserved and the 31 sub-basins are divided into 2,947 HRUs. After the hydrological and meteorological data (2009–2018) and agricultural management data have been imported, the model can be simulated. In order to ensure the accuracy and stability of the simulation, the warm-up period was set as 2 years (2009 and 2010).

The parameter sensitivity analysis was carried on using the Sequential Uncertainty fitting (SUFI-2) algorithm in SWAT-CUP (SWAT-Calibration and Uncertainty Programs), and the analysis was carried out on a monthly time step. The results of sensitivity analysis were measured by two indexes, t-stat and p-value. After multiple iterations within the set parameter range, the greater the absolute value of t-stat and the closer the p-value is to 0, the stronger is the sensitivity of the parameter. In this study, four indexes of runoff, sediment, TN and TP were analyzed for parameter sensitivity respectively, and the models were iterated 500 times. The sensitivity parameters of each index are shown in Table 3. The results show that the higher sensitivity parameters for the runoff are ALPHA_BF, CN2, SLSUBBSN, GWQMN, CH_N2 and so on. The higher sensitivity parameters for sediment are OV_N, CN2, USLE_P and so on. The higher sensitivity parameters for TN are ESCO, CN2, ERORGN and so on. The higher sensitivity parameters for TP are CN2, ERORGP, ALPHA_BF and so on.

Note: n is the total number of data series, O_i is the observed data, S_i is the simulated data of the SWAT model, Ō is the average of observed data, and S̄ is the average value of simulated data of the SWAT model.

In this study, the requirement of runoff is R² > 0.8, E_NS > 0.7 and Re is within the range of 25%. For sediment and nutrients, R² > 0.6, E_NS > 0.6 and Re is within the range of 40%. Since the water quality data (monthly pollutant concentration) collected for the calibration of the SWAT model includes both point source (PS) pollution and NPS pollution, it was necessary to separate the NPS pollution. Firstly, the measured runoff was divided into surface runoff and base flow by digital filtering technology. Secondly, the monthly pollutant load was obtained by multiplying the average monthly surface runoff by pollutant concentration. Then, the pollutant concentration in the dry season of the year is taken as the PS concentration and multiplied by the base flow to get the PS pollution load of monthly pollutant. Finally, the monthly NPS pollution load is obtained by subtracting the monthly PS pollution load from the monthly total pollutant load.

According to the agricultural management mode, the spatial characteristics of NPS pollution and the topographic characteristics of the Hanjiang River Basin, three types of management practices, namely engineering management practices, non-engineering management practices and landscape management practices were set up. The non-engineering management practices include fertilizer reduction, no tillage and stubble coverage. Engineering management practices include vegetation buffer zone, grassed waterway, contour hedgerow and terraces fields. Landscape management is a practice of returning farmland to forest land. In this paper, the initial status of the model is set to scenario 0, and the settings and parameters of the eight scenarios are shown in Table 4. Among them, the parameter settings of USLE_P and SLSUBBSN are shown in Tables 5 and 6. (Arnold et al. 2012).

Note: R is reduction rate; L_BAS is TN or TP load under the initial scenario; L_BMPs is TN or TP load after introducing BMPs.

The economic benefit evaluation usually considers the cost of BMPs. The multi-attribute decision-making method based on information entropy is used to evaluate the comprehensive benefits of BMPs. The calculation steps are as follows (Sun et al. 2013).

Step 1: Construction of the BMPs attribute decision matrix:

The attribute factors of BMPs are constructed into decision matrix A= (a_ij)_n×m, where a_ij represents j attribute factors in the i-th BMPs; n is the number of BMPs; m is the number of attribute factors.

Step 2: Standardized processing of decision matrix:

The larger the calculation result of the comprehensive attribute value Z, the better the comprehensive control effect of the BMPs.

For runoff and TP indexes, the calibration period is from 2011 to 2014 and the verification period is from 2015 to 2018. However, due to the limited observed data collected, the sediment was calibrated from 2011 to 2013, and the verification period was 2014 and 2016. The TN was calibrated from 2013 to 2014, and the verification period was 2015 and 2018. The calibration and verification results of the SWAT model are shown in Table 7 and Figure 5. The results show that, in the calibration period, the Re of monthly runoff between the simulated values and the measured values was within 25%, the R² and E_NS of monthly runoff was 0.96 and 0.89. The Re of monthly sediment between the simulated values and the measured values was within 10%, the R² and E_NS of monthly sediment was 0.73 and 0.72. The Re of monthly TN between the simulated values and the measured values was within 40%, the R² and E_NS of monthly TN was 0.85 and 0.68. The Re of monthly TP between the simulated values and the measured values was within 30%, the R² and E_NS of monthly TP was 0.90 and 0.79. In the verification period, the Re of monthly runoff between the simulated values and the measured values was within 25%, the R² and E_NS of monthly runoff was 0.82 and 0.82. The Re of monthly sediment between the simulated values and the measured values was within 20%, the R² and E_NS of monthly sediment was 0.77 and 0.70. The Re of monthly TN between the simulated values and the measured values was within 30%, the R² and E_NS of monthly TN was 0.75 and 0.71. The Re of monthly TP between the simulated values and the measured values was within 20%, the R² and E_NS of monthly TP was 0.72 and 0.52. Therefore, the runoff, sediment, TN and TP basically met the accuracy requirements of the model.

In conclusion, the SWAT model can be used to simulate the NPS pollution in the study area. The relevant sensitivity parameters were finally determined, as shown in Table 8.

Rainfall and runoff is an important part of the hydrological cycle and is the driving force and carrier of NPS pollution output. Figure 6 shows the process of rainfall and runoff from 2011 to 2018. Rainfall is mainly concentrated in May to September of each year, and runoff is significantly affected by rainfall, showing an obvious positive correlation. The rainfall amount is large in the wet season, while small in the dry season. The proportion of rainfall accounts for 72.84% of the whole year, and the proportion of runoff is 66.30%. In the study area, the regulation and storage effect of the reservoir has changed the natural runoff process, which results in abnormal discharge.

Figures 7 and 8 shows the simulation results of TN and TP loads at the annual and monthly scales in Ankang Hydrological Station. It can be seen that the TN and TP loads show certain fluctuations under the influence of rainfall, and the change law is basically the same, but the overall TN load is greater than the TP load. The output loads of TN and TP vary greatly from year to year, and the maximum values both appear in 2011, which are 13,022.24 t and 738.29 t. The minimum values both appear in 2016, which are 7,260.35 t and 435.41 t, respectively. The pollution load during the year was mainly concentrated in June to October, and TN and TP loads account for 84.23% and 91.54% of the whole year, indicating that rainfall–runoff is closely related to the generation of NPS pollution load.

The production of NPS pollution has pronounced spatial characteristics. It is closely related with various features of the study area, such as rainfall, land use, soil type, topography, and so on. As a distribution model, the SWAT model combined with ArcGIS can analyze the spatial distribution of the NPS pollution in the study area. This paper analyzes the spatial distribution characteristics of rainfall, runoff, sediment, TN and TP averaged over the years from 2011 to 2018, as shown in Figure 9.

It can be seen from Figure 9(a) that the southeast region has the highest rainfall, such as the 19, 27, 28, 29 and 30 sub-basins, with an average rainfall of about 1,000 mm. The areas with the least rainfall are mainly concentrated in two, four, eight and nine sub-basins, with an average rainfall is about 800 mm. Comparing Figure 9(a) with Figure 9(b), it was found that the rainfall and runoff in the basin showed a certain positive correlation. The rainfall in the Nos. 15, 19, and 27 sub-basins was heavy, and the slope is above 20°, so the corresponding runoff is also relatively large. It can be seen from Figure 9(c) that the distribution of sediment erosion intensity in the basin is extremely uneven. The sub-basin with the highest sediment erosion intensity is as high as 8.140 t/ha, and the lowest sub-basin is only 0.210 t/ha. There is a strong correlation between sediment erosion intensity and runoff, that is, it is mainly concentrated in the area with large rainfall and runoff. In short, the spatial distribution of sediment erosion intensity is related to rainfall, runoff, terrain slope, soil type and other factors.

As is shown in Figure 9(d) and 9(e), we can find that the loss of TN and TP load is mainly concentrated in the southern part of the basin. The distribution of loss intensity of TN load is significantly correlated with runoff, which is mainly distributed in sub-basins Nos. 15, 16, 19, 23 and 27. Nitrogen is dominated by ammonium nitrogen and nitrate nitrogen, both of which are easily soluble in water. Therefore, nitrogen is mainly exported with runoff. The distribution of loss intensity of TP load is consistent with that of sediment erosion intensity, and the most serious area is located in the No. 27 sub-basin. Phosphorus is chemically inactive, not easily soluble in water and is generally adsorbed on the surface of the soil. Therefore, under the action of rainfall–runoff, phosphorus is transported into the river with sediment.

According to the statistical results, the loss of nitrogen and phosphorus in sub-basins 6, 15, 16, 19, 21, 26, 27 and 28 is the most serious, with a loss intensity of TN load of 3.025–4.693 kg/ha and a loss intensity of TP load of 0.113–0.317 kg/ha. The TN and TP loads of the eight sub-basins account for 48.62% and 53.22%, which is nearly half of the output load of the whole basin. Therefore, the eight sub-basins were identified as the CSA of NPS pollution in the study area.

• (1)

  Fertilizer reduction

Figure 10 shows the changes in TN and TP loads under non-engineering management practices, where Scenario 1 is fertilizer reduction. Comparing Scenario 0 and Scenario 1, it was found that the reduction effect of 20% fertilizer reduction on NPS pollution load is not obvious, the reduction rate of TN load was 1.15–2.33%, and the average reduction rate was 1.65%. The reduction rate of TP load was 7.93–11.05%, and the average reduction rate was 9.46%. Sun (2020) simulated chemical fertilizer reductions of 10, 30, and 50% in the Chenghai Basin Bain. The corresponding nitrogen load reduction rates were 0.88, 2.66 and 4.63%, and phosphorus load reduction rates were 0.1, 0.32 and 0.53%, respectively. While Liu et al. (2014) believed that a 20% reduction in chemical fertilizers can reduce 50% of the nitrogen load and 40% of the phosphorus load, and the reduction effect was ideal in Xiangxi River Basin. Therefore, the reduction effect of fertilizer on NPS pollution load has certain regional differences, which is mainly affected by rainfall conditions and crop management methods. The study area is mainly dominated by agricultural production. Long-term excessive fertilization results in the accumulation of farmland nutrients, which makes it impossible to significantly reduce the nitrogen and phosphorus content in the soil over the short term, and has limited effect on the reduction of the NPS pollution load in the basin.

• (2)

  Farming management practices

Scenario 2 and Scenario 3 are no tillage and stubble coverage respectively. It can be seen from Figure 10 that the reduction effect of stubble coverage was significantly better than that of no tillage. The reduction rate of no tillage on TN load was 2.19%–9.44%, while the reduction rate on TP load was negative. No tillage increases the output of dissolved phosphorus, so it can increase the output of TP load to a certain extent. The reduction rate of TN load by stubble coverage was 24.93%–53.52%, and the reduction rate of TP load was 60.68%–78.84%.

• (1)

  Vegetation buffer zone

Figure 11 shows the changes of TN and TP loads under non-engineering management practices, where Scenario 4 is the vegetation buffer zone. The vegetation buffer zone has a better reduction effect on CSA. The reduction rate of TN load by the vegetation buffer zone of 5 m width was 16.70%–31.23%, and the reduction rate of TP load is 27.70%–38.05%. The mechanism of the vegetation buffer zone to filter NPS pollution is mainly physical interception, plant utilization, microbial transformation and soil adsorption. Therefore, the reduction rate of the vegetation buffer zone on the phosphorus load is better than the nitrogen load, which is consistent with the results of this paper, that is, phosphorus mainly exists in the form of particles, and is easier to be intercepted by the vegetation buffer zone (Fu et al. 2019). Therefore, this study comprehensively considers various factors, selects the vegetation buffer zone with a width of 5 m and achieves a more ideal reduction effect on NPS pollution.

• (2)

  Grassed waterway

Among the engineering management practices in the study area, the grassed waterway (ie Scenario 5) has the best control effect on NPS pollution. The reduction rate of TN load by grassed waterway in CSA was 46.71%–80.15%, and the average reduction rate was 65.62%. The reduction rate of TP load is 43.97%–83.46%, and the average reduction rate was 66.77%, indicating that the grassed waterway has the same reduction effect on nitrogen and phosphorus loads.

• (3)

  Contour hedgerow

Contour hedgerow is planted along contour lines to slow down velocity and reduce soil erosion. It can be seen that the reduction rate of TN load by contour hedgerow (Scenario 6) is 20.46%–40.59%, and the average reduction rate was 31.17%. The reduction rate of TP load was 36.16%–55.02%, and the average reduction rate was 44.74%.

• (4)

  Terraces field (<25°)

The terraces field is an effective soil and water conservation practice in controlling soil and water loss on the sloping farmland. After simulating the terraces field (Scenario 7), the reduction rate of TN load in the CSA reached 5.39%–43.70%, and the reduction rate of TP load was 3.83%–31.37%. From the perspective of the scope of reduction, the reduction efficiency of terraces field has obvious spatial differences. Therefore, for areas with large slope farmland, terraces field can effectively control runoff, thereby reducing the output of NPS pollution caused by soil erosion.

In summary, the engineering management practices have a good reduction effect on the nitrogen and phosphorus loads, and the reduction efficiency from high to low is in the order of grassed waterway > contour hedgerow > vegetation buffer zone > terraced field (<25°).

Returning farmland to forest land directly reduces the area of sloping farmland and effectively reduces the risk of soil erosion in steep slope areas, thereby reducing the output of NPS pollution to a certain extent. It can be seen from Figure 12 that the reduction rate of TN load in areas above 25° (Scenario 8) was 22.59%–65.96%, and the average reduction rate was 45.16%. The reduction rate of TP load was 21.87%–84.58%, the average reduction rate was 59.53%. Compared with engineering management practices, its reduction effect is second only to the grassed waterway.

Individual BMPs has a significant difference in the reduction effect of nitrogen and phosphorus loads. Considering the requirements of joint prevention and control of NPS pollution, the management practices with better reduction rate were selected from engineering management practices, non-engineering management practices and landscape management practices, and three combined BMPs were set, as shown in Table 9.

As shown in Figure 13, in the CSA, the three combined BMPs have a considerable reduction effect on nitrogen and phosphorus loads. The reduction rate of TN load was about 60%, and the reduction rate of TP load was more than 80%. Among them, the reduction effect of ‘stubble coverage + grassed waterway + returning farmland to forest land (Scenario 10)’ is the best, with an average reduction rate of 67.34% for TN load and 83.40% for TP load. The second is the reduction effect of ‘stubble coverage + vegetation buffer zone + returning farmland to forest land (Scenario 9)’, with the reduction rates being 63.51% and 83.36%, respectively. Comparing the reduction effect of individual BMPs, it can be found that the control effect of combined BMPs on nitrogen and phosphorus loads is not the sum of the reduction efficiency of individual BMPs, and the combined reduction effect is less than the sum of the individual reduction efficiency. The reason may be that vegetation buffer zone, grassed waterway, and contour hedgerow all intercept runoff and slow down the velocity, thereby controlling the transport of NPS pollution in the runoff. Therefore, the selection of NPS pollution management practices in the basin should be adapted to local conditions, grouped reasonably, and classified to reduce NPS pollution load.

Combined with the relevant research and the subsidy standard of the Yangtze River Basin, the annual input cost of BMPs was determined. The technology of reducing the amount of fertilizer application involves the costs of construction, scientific research, etc., and it is difficult to quantitatively allocate to the unit area. Therefore, this paper does not consider the BMPs of fertilizer reduction (Scenario 1). Combined with the simulation results in Section 3.3, the comprehensive attribute value Z is obtained, as shown in Table 10.

It can be seen from Table 10 that the reduction of NPS pollution load by different BMPs differs greatly. In the evaluation of environmental benefits, the grassed waterway has the most obvious reduction effect on TN, with a reduction of 2.648 kg/ha. The reduction effect of stubble coverage on TP is the most obvious, with a reduction of 0.142 kg/ha. In the evaluation of cost–effectiveness, the input cost of the terraces field is the largest and the grassed waterway is the smallest. For the comprehensive evaluation results, among the single BMPs, the Z value of grassed waterway in the engineering practices is higher, reaching 0.93, followed by the stubble coverage in the non-engineering practices, with the Z value of 0.73. The last is the BMPs of returning farmland to forest land with the Z value of 0.70. Although the combined BMPs have high reduction rates on NPS pollution, their high cost input resulted in a relatively low Z value of 0.81, 0.84, and 0.65.

In this paper, a SWAT model of Ankang section of the Hanjiang River Basin is established based on the collected spatial database and attribute database. The results showed that R², E_NS and R_e meet the requirement of model accuracy, and thus can be used to study the characteristics of NPS pollution. The simulation results of the SWAT model showed that TN and TP loads in flood season are significantly higher than that in the non-flood season, and have a positive correlation with the rainfall. The distribution of loss intensity of TN load was strongly correlated with runoff, and the distribution of loss intensity of TP load is significantly correlated with that of sediment erosion intensity. According to the loss intensity distribution and loss ratio of TN and TP load, the southern region (sub-basins Nos. 6, 15, 16, 19, 21, 26, 27 and 28) was determined to be the CSA of NPS pollution. Considering the agricultural management mode, spatial characteristics and terrain characteristics of NPS pollution loss in Hanjiang River Basin, 11 best management measures are set up in the CSA, and their environmental and comprehensive benefits are evaluated. The results of environmental benefit evaluation showed that the reduction rates of stubble mulch, grassed waterway and returning farmland to forest land were relatively high among the eight individual BMPs,. The reduction rates of three combined BMPs are ‘Stubble coverage + Grassed waterway + Returning farmland to forest land (>25°)’ > ‘Stubble coverage + Vegetation buffer zone + Returning farmland to forest land (>25°)’ > ‘Stubble coverage + Contour hedgerow + Returning farmland to forest land (>25°)’. The comprehensive benefit evaluation results showed that the comprehensive attribute value Z of stubble coverage is the highest, which is 0.93. Followed by the combined BMPs of ‘Stubble coverage + Grassed waterway + Returning farmland to forest land (>25°)’ with the attribute value Z at 0.84.

Overall, stubble coverage is the most effective way to control NPS pollution in the river basin. It has the highest environmental benefits, lower cost input and the highest comprehensive benefits, so it can be widely used. At the same time, in the areas near the river, grassed waterway can be selected, and in the area with small forest area, the measure of returning farmland to forest can be selected. However, in the future, we should comprehensively consider social and economic factors such as regional population, and urban distribution, and strive to maximize the comprehensive benefits of the BMPs.

The study was financially supported by the key research and development project of Shaanxi Province (2019ZDLSF06-01) and the National Natural Science Foundation of China (51879215).

All relevant data are included in the paper or its Supplementary Information.