Non-point source pollution and long-term effects of best management measures simulated in the Qifeng River Basin in the karst area of Southwest China

The deterioration of water quality caused by non-point source (NPS) pollution has become a hot research topic in the world NPS pollution was mainly due to agricultural activities (Carpenter et al. 1998; Li et al. 2015). Excessive nitrogen and phosphorus caused by NPS pollution entering a water body will cause eutrophication and then pose a large threat to the ecosystems (Pease et al. 2010; Schindler et al. 2016). So far, NPS pollution has posed serious challenges to the effective and efficient design and implementation of global management policies (Shortle & Horan 2017). In China, karst landform area accounts for more than one-third of China's land area (Song et al. 2017). The largest and most complete karst landform is located in the southwestern part of China. In China, agricultural NPS pollution has been a major problem since the 1990s (Ongley et al. 2010; Wang et al. 2018a, 2018b). The Qifeng River Basin is a typical river in the karst area of Southwest China, which is located in the agricultural production center of Lingui City. Due to the special geological structure of karst areas and the uncertainty of the applicability of the model, NPS pollution in karst areas has attracted significant attention. Thus, a comprehensive and careful investigation of NPS pollution in karst areas is urgently needed in China.

NPS pollution simulation models mainly include empirical models and mechanism models. Process-based mechanism model provides a useful alternative for pollutant load apportionment, which are capable of simulating the movement of water, as well as the transport and transformation of water pollutants from different pollution sources within a river basin (He et al. 2019). Many process-based mechanisms models, including the Soil and Water Assessment Tool (Merriman et al. 2018) (SWAT) model, Hydrological Simulation Program Fortran (Srinivas et al. 2020) (HSPF) model, Hydrological Predictions for the Environment (HYPE) model (Goran et al. 2010), Integrated Catchments – Nitrogen dynamics (INCA-N) model (Whitehead et al. 2013), and AnnAGNPS model (Wang 2014), have been applied to estimation for the nutrient cycles in diverse watersheds. Although SWAT model has been used for simulating runoff and water quality of NPS pollution in karst areas, large errors may exist according to Malago et al. (2016) and Amatya et al. (2011), therefore, cumbersome calibration processes are required on the base of the characteristics of karst areas to obtain a good simulation on NPS pollution. The outstanding performance of AnnAGNPS model on simulating runoff, sediment, and nutrient load in different watersheds have been acknowledged since it integrates the powerful Geographic Information System (GIS) of spatial analysis to achieve the spatial distribution of NPS pollution in one basin. AnnAGNPS has been used satisfactorily in numbers of cases to simulate runoff, sediment, and nutrient load in non-karst area watersheds (Li et al. 2015; Villamizar & Brown 2016; Karki et al. 2017; Romano et al. 2018; Jiang et al. 2019; Momm et al. 2019). Despite the acknowledged advantages of AnnAGNPS model for simulation of NPS pollution in non-karst watershed, it should be appropriately adjusted according to karst characteristics with a good simulation of NPS pollution. Therefore, it is necessary to discuss the applicability of the AnnAGNPS model in karst areas and to evaluate the robustness of the model in karst areas.

Best management measures (BMPs), by changing or affecting processes such as hydrology, soil erosion, ecology, and nutrient cycling in the watershed, have been one of the most effective means of NPS pollution control (Tian et al. 2010; Dakhlalla & Parajuli 2016). At present, the watershed model (such as SWAT, and AnnAGNPS) can help to select and provide the BMPs decision process that most effectively reduces the pollutant load for the management makers (Mtibaa et al. 2018; Qi et al. 2018). However, the previous studies mainly focused on the single measure with a single scenario simulation and pay less attention to a single measure with multi-scenarios simulation. In addition, short-term simulation cannot fully reflect the reduction effect of management measures due to the hysteresis quality of NPS pollution. The long-term effects of the entire river basin management measures should also attract attention. It is envisioned that a long-term simulation method can cover the various dynamic conditions under which the BMPs performs. Therefore, long-term simulation has better accuracy than event-based evaluation of BMPs. Moreover, using AnnAGNPS to quantitatively investigate the long-term reduction efficiency of BMPs can bring clearer guidance to the treatment of NPS pollution in karst areas' watershed.

The objectives of this study were: (1) to test the applicability of AnnAGNPS model in the karst agricultural basin of the Qifeng River watershed; (2) to analyze the temporal and spatial distribution characteristics of NPS pollution in the watershed; and (3) to assess the effectiveness of nutrient loads reduction quantitatively by a long-term simulation approach to give best management practices for NPS control.

The Qifeng River watershed located in the southwestern Guilin City of Guangxi Province in China (110°13′11.29′′–110°22′23.59′′E, 24°55′36″–25°6′47″N), and the administrative area belongs to Lingui County of Guilin City and Yanshan District of Guilin City (Supplementary information (SI), Figure S1). The study area is a typical karst area in southwestern China. The total area of the Qifeng River watershed is 22,609.71 ha. Rainfall ranges from 1,313.3 mm to 2,452.7 mm, with the average annual rainfall being 1,835.8 mm. The average annual temperature is 19.5 °C. Guilin can be characterized by three seasons (according to rainfall, namely, dry, wet, and normal seasons); these seasons generally occur from January to April, May to August and September to December, respectively.

In the 1980s, the US Department of Agriculture's Agricultural Research Service (USDA-ARS) and the Natural Resources Defense Agency (NRCS) jointly developed the AnnAGNPS model, a continuous simulation of distributed parameter models with daily time-step pollutant loads (AnnAGNPS Technical Processes). It was the continuation of the single event model AGNPS (Young et al. 1989). The model was developed to predict surface runoff, sediment erosion and the transfer of nutrients such as nitrogen and phosphorus in agricultural watersheds. The AnnAGNPS version 5.50 was used in this study (NE.Ref).

The Digital Elevation Model (DEM) used in this study was taken from the Geospatial Data Cloud (NE.Ref) of the Computer Network Information Center of the Chinese Academy of Sciences (SI, Figure S1(d)). Due to the medium size of the study area, a DEM data source with a spatial resolution of 30 m was selected. DEM is utilized to automatically generate the water system trend of the watershed. The characteristic of cells and reaches is determined by using Terrain Optimization Parameters Annualized Agricultural Non-Point Source (TOPAGNPS) and Agricultural Flow (AGFLOW) components of the model. AnnAGNPS cells and reaches in the watershed were delineated based on the critical source area (CSA) and a minimum source channel length (MSCL). The Qifeng River watershed study area was finally divided into 1225 AnnAGNPS cells with 504 reaches by using a CSA of 20 ha and a MSCL of 200 m (SI, Figure S2(a)).

A digital soil map was attained from the Lingui County Land and Resources Bureau. The original map is a vector diagram with a scale 1:400,000. The dominant soil type for each sub-watershed cell was determined through overlaid soil map onto the map of the delineated watershed using Agricultural Geographic Information System (ArcGIS). Paddy soil is the major soil type in the watershed (SI, Figure S2(b)), and it consists of four soil layers with a top layer depth of 20 cm and a total soil depth of 100 cm. The data come from the second national soil census (NE.Ref), and the soil hydrological parameters are calculated using the SPAW software. Detailed properties of soil in the Qifeng River watershed are listed in Table S1.

A digital land use map was attained from the Lingui County Land and Resources Bureau. AnnAGNPS/AGNPS GIS Tools was used to determine the dominant land use in each cell. The land use data used in each AnnAGNPS cell is shown in Figure S2(c) (SI). The cultivated land is divided into paddy fields and dry land, accounting for 19.79 and 25.59% of the study area, respectively. The watershed is mainly agricultural land with fields planted with rice and corn. The garden area accounts for 14.72%, forest land accounts for 25.57%, and the residential land accounts for 11.64%.

Management data is comprised of management field, management schedule data, management operation data, crop data, non-crop data, fertilizer data and pesticide data. Management field data encompass the management schedule implements for a management field, the land use type and the first year of rotation. Management schedule data allow to schedule the events for cropland and non-cropland. Management operation data include the detailed information of farming operation such as the time of disk, fertilize, spray and harvest. The management measures are listed in Tables S2, S3, and S4. The management data were all obtained from the questionnaire result from households and statistical yearbook data from government in this study.

The AnnAGNPS model predicts surface runoff based on the precipitation and evaporation. Meteorological data required in the model consists of six daily parameters, namely, maximum temperature, minimum temperature, precipitation, dew point temperature, sky cover or solar radiation and wind speed. The weather data was attained from China Meteorological Data Sharing Network. The global storm type was type II based on the rainfall distribution at Lingui station.

In this study, runoff, nitrogen, and phosphorus load at the watershed outlet were calibrated to optimize the model parameters input so that the simulated results were more in line with the measured results. The measured data at outlet of watershed were taken from Guilin Hydrology and Water Resource Bureau (NE.Ref). The AnnAGNPS model performance was evaluated quantitatively by using the coefficient of determination (R²), and Nash-Sutcliffe Efficiency (NSE), and Relative error (Re).

According to the previous research (Liang et al. 2020), referring to the results of the best sub-basin division and parameter sensitivity analysis, four management measures (viz., fertilizer reduction (FR), soil and water conservation P factor (P), stubble tillage (ST), and less tillage and no tillage (LT)) were selected and implemented in the AnnAGNPS model. The specific measures of scenario simulation are given in Table S5.

Various studies have shown that the curve number (CN) is the most sensitive parameter for runoff (Chahor et al. 2014; Villamizar & Brown 2016). This study adjusted the CN value for all land use to make the predicted runoff close to the measured runoff. It was found that CN value is significantly higher than the value of CN under non-karst watershed.

As a result, during the calibration period, the simulation annual runoff was 213.77 m³/s, which was similar to the observed value (258.97 m³/s) at the outlet of watershed. For the validation period, the simulated and measured runoff value at the outlet were 240.31 m³/s and 249.6 m³/s, respectively. The Re between simulated and measured runoff were −17.46 and −3.72%, respectively.

Figure 1(a) shows that monthly runoff simulated by the model was consistent with measured value at the outlet during calibration and validation. Though the predicted monthly runoff were under- and over-predicted in some months, the model still showed a good agreement on the whole compared to measured value with an overall R² value of 0.91, NSE value of 0.83, and Re value of −14.94% during the calibration. The statistical analysis during validation period showed very good results with R² = 0.98, NSE = 0.95, and Re = −10.86%, respectively (Liang et al. 2020).

The AnnAGNPS model simulated runoff with a R² value of 0.71, NSE value of 0.89, and Re value of −19.25% during the calibration period, respectively (Figure 1(b)). The statistical analysis during validation period showed good results with R² = 0.72, NSE = 0.96, and Re = −13.69%, respectively (Liang et al. 2020). The model underestimates the simulation of surface water on a daily scale.

The observed monthly TN and TP from 2017 to 2018 were used to calibrate and validate the model (Figure 2(c) and 2(d)). In the calibration period, the model showed a good result of R² (0.88), Re value of −31%, and an acceptable value of NSE (0.56), respectively. After calibration, the R² (0.93) value of the model was increased, but the values of NSE (0.53) and Re (−32%) were slightly lower than the calibration period. The simulated value was below the observed data, indicating a underprediction of TN loads with a quite high R², but a quite low NSE. The model underestimated TN loads both during calibration and validation periods.

The R² value of TP between simulated and observed data at monthly scale was 0.93 both on calibration and validation periods (Figure 2(d)). However, NSE did not demonstrate a very good value of 0.62 and 0.57 on calibration and validation periods, respectively. During calibration and validation periods, the Re value was 29%.

The temporal distribution of simulated runoff, TN, and TP compared with measured precipitation with a coefficient larger than 0.9 is shown in Figure 3. It manifested that the model can well describe the surface runoff generated by rainfall and well described the trend of rising (rainy season) or falling (dry season). TN and TP showed a high correlation (R² > 0.9) with rainfall, which further verified that rainfall was the driving factor and carrying medium of NPS pollution in rivers. The temporal distribution of TN and TP in April–August was significantly higher than other months and accounted for more than 80% of the annual year. The rainy season of the Qifeng River Basin from April to August increased in the surface runoff, so the pollutants increased with the surface runoff.

As can be seen from Figure 4, the spatial distribution of simulated runoff, TN and TP was roughly the same, but the values of pollutants were different, which indicated that TP and TN had the same driving factor and pollution source. The annual average runoff outputs of all cells were between 0.01 and 38,979.18 m³yr⁻¹ with a spatial distribution divided into five levels (Figure 4(a)). The annual average TN and TP outputs of the whole watershed ranged from 1.63 to 51,930.63 kg yr⁻¹, and 0.54 to 18,317.18 kg yr⁻¹, respectively. There were 12 cells with the largest pollution output, accounting for 0.98% of the total, while 696 cells were observed with the lightest pollution, accounting for 56.82%. Overlay analysis of the relationship between land use type and pollutant output was conducted using ArcGIS. The combination of simulated spatial distribution and overlay analysis of NPS pollution showed that the largest pollutant output was observed in orchard land of each sub-basin.

In order to determine the reduction rate of different BMPs to reduce NPS pollution, ten years simulation for four single measures with multi-scenarios were conducted on the entire river basin. The reduction rate of NPS pollution loads for each BMPs scenario was compared to the benchmark scenario.

The trends of annual average reduction rates of sediment, TN, and TP are shown in Figures 5 and 6. Application of the FR had no effect on reducing runoff and sediment but had a strong positive correlation with TN and TP. The average annual reduction rates in scenario 3 of TN and TP were 19.51 and 19.97%, respectively. The long-term reduction rate of TN fluctuated with time, while the reduction rate of TP had steadily increased. Nitrate nitrogen has strong penetrability and will be engulfed into the watershed by stable baseflow and soil runoff. The reduction efficiency of TN is immediately attributed to the weak accumulation of TN in the soil. In contrast, TP is easy to deposit and adhere to the soil due to its high accumulation. In addition, excessive application of phosphate fertilizer each year will cause the accumulation of TP in the soil year by year. Therefore, the treatment of TP takes a long time to show good results.

The 10-year average reduction rates of sediment, TN, and TP were gradually increasing as shown in Figure S3 (SI). The implementation of P factor can effectively reduce the loads of sediment, TN, and TP with the average rates of 48.44%, 29.37%, and 28.83%, respectively. The aim of the P factor was mainly to reduce the nitrogen and phosphorus carried by sediment by reducing sediment erosion. The long-term reduction rates of TN and TP in the simulation period are illustrated in Figure S4 (SI). Scenarios 1 and 2 had a much smaller effect than scenario 3 which showed P factor played a more important reduction efficiency in woodland. Scenario 3 implemented P factors on dry land + paddy field + woodland, and the area was greater than that of dry land and paddy field in scenarios 1 and 2. The trends of TN and TP reduction rates were the same in that the reduction rates of scenarios 1 and 2 were relatively stable, while the reduction rate of scenario 3 was gradually decreasing.

As can be observed in Figure S5 (SI), the trends of 10-year average reduction rates in sediment, TN, and TP by the ST scenario simulation were the same. The maximum 10-year average reduction rates of sediment, TN, and TP reached 16.64%, 17.845%, and 17.80%, respectively. The long-term effect on sediment, TN, and TP displayed that the reduction rate fluctuated with time, and there was no specific trend of increasing or decreasing with time (SI, Figure S6).

The results of LT scenario simulation in sediment, TN, and TP are shown in Figure S7 and Figure S8 (SI). The results indicated that LT could only slightly reduce the output of sediment, TN and TP. Reducing the disturbance area can reduce the oxidation of soil organic matter and reduce the damage of soil aggregates. Different perturbed areas during cultivation will produce different degrees of soil compaction and change soil porosity. In this study, the disturbance area of the soil was reduced, but the reduction of sediment, TN, and TP was not great.

The reduction rates under four single measures are compared in Table 1. It can be seen from the table that the reduction rates of TN and TP were as follows: P factor >FR > ST > LT, while the ranking of the reduction rate of sediment was P factor >ST > LT > FR.

The AnnAGNPS model performed well in runoff estimation in a karst area watershed, given that the measured and simulated value had a great agreement (Figure 1(a) and 1(b)). Owing to the characteristics of karst areas, surface runoff is easy to leak into the ground and the retention amount of runoff is less than that of non-karst areas. By increasing the value of CN, the initial retention is reduced, which was in line with the characteristics of karst areas. It can be considered that the difference between karst region and non-karst region lies in the value of CN. According to Chiew et al. (1993) for the simulation of runoff, it was considered that the model's R² and NSE are greater than 0.5, indicating that the model was acceptable. These relatively high R² and NSE values between simulated and observed monthly and daily runoff indicated that the AnnAGNPS model satisfactorily predicted runoff at the Qifeng River watershed. The results of this study were similar with many existing results in non-karst areas and karst areas using another model (SWAT) (Zhang et al. 2019), in which the simulation of runoff showed good results. Villamizar & Brown (2016) conducted that the AnnAGNPS was applied to simulate surface runoff in the river Cauca with R² = 0.73 and NSE = 0.70. Chahor et al. (2014) implemented the AnnAGNPS model to predict runoff obtaining R² = 0.75 and NSE = 0.75 at a monthly scale in a small Mediterranean agricultural watershed in Spain. A good performance of AnnAGNPS model also were conducted in many watersheds in non-karsts areas, such as Hanalei River Basin (Polyakov et al. 2007), Heigou River watershed (Tian et al. 2010), Taihu Lake basin (Li et al. 2015), and an agricultural watershed in East-Central Mississippi (Karki et al. 2017).

However, the AnnAGNPS model underestimated nitrogen with an acceptable result (NSE > 0.5, R² > 0.9). On the one hand, the under-predicted value was attributed to the limitation of model and the uncertainty of parameters. On the other hand, nutrients were lost to the groundwater aquifer through sewage pits and pipes in karst areas. Yan et al. (2018) used artificial rainfall to simulate the underground leakage in the karst area of Guizhou Province. They found that the leakage rate was slow, but the effect could not be ignored. Song et al. (2017) found that the N-loss load was higher for subsurface runoff than surface runoff in the karst area. Many studies have shown that groundwater nitrate pollution is severe in a karst area. The nitrogen has been transferred to the underground aquifer by subsurface runoff.

Compared with existing studies on nutrient simulation through AnnAGNPS model, this study showed similar or better nutrient simulation results. Yuan et al. (2003) reported that the model overpredicted nitrogen during the winter season and underpredicted nitrogen loadings during crop season. Li et al. (2015) applied AnnAGNPS to simulate phosphorus with a value of R² (0.60 and 0.83) but a lower NSE (0.61 and 0.5) during calibration and validation periods in a typical small watershed in China. Karki et al. (2017) conducted that R² of 0.74 and NSE of 0.54 for monthly phosphorus estimation was obtained in East-Central Mississippi. Nevertheless, the result of the simulation of nitrogen was not satisfactory. A poor performance of the AnnAGNPS model in simulating nutrient loadings was also reported by Suttles et al. 2003; Shamshad et al. 2008; Pease et al. 2010. These researchers indicated that it was difficult to predict the short-term nutrient loadings due to the simplification of nutrients processes and uncertainties in the input parameters. The simulated results represent a portion of the nutrient loading and do not fully estimate the nutrients from all sources. Based on mass conservation, nutrient loading will be affected with any missing input or output parameters in a watershed. Most model evaluations assume absolute quality of the observed data, however, measured data are often accompanied by errors because of different sources of uncertainty in a sample collection, handling and analysis (Baginska et al. 2003). Shamshad et al. (2008) pointed out that an R² of 1 for nutrient loading is largely impossible, because the AnnAGNPS model assumed that all the simulated components reach the outlet at the beginning of the simulation on the second day. Due to the complexity of a karst areas, more simulation experiments are needed to make the AnnAGNPS model better applicable in karst areas in future.

In conclusion, the AnnAGNPS model can be successfully used to estimate runoff and nutrient from agricultural watersheds in the karst area.

The differences in land use type, management measures, rainfall and climate in the watershed determine the temporal and spatial variability of NPS pollution (Emirhüseyinoğlu & Ryan 2020; Marin et al. 2020; Ricci et al. 2020). The temporal distribution of TN and TP has the largest output during the rainy season, accounting for more than 80% of the year. This phenomenon can be explained by the fact that heavy rainfall during the wet season caused strong soil erosion (Gao et al. 2019), and TN and TP loss significantly increased. This was similar to the existing study (Song et al. 2017) conclusions in which N-loss and rainfall are highest during the rainy season (April to September). Yue et al. (2015) stated that approximately 85% of nitrate transports occur during the wet season by isotope analysis of ¹⁵N and ¹⁸O in NO³⁻, which is in good accordance with our study. Existing research showed that TP transport significantly increased with elevated rainfall intensity (Gao et al. 2010). On the spatial scale, the high output of TN and TP were mainly distributed in orchard land due to the increase in fertilizer application and the steep terrain which was prone to soil erosion on the slope. By contrast, the orchard land in a typical small watershed in the Three Gorges Reservoir area was not the highest output area which meant that the same land use with different fertilization management techniques affected the reduction results of NPS pollution (Gao et al. 2019). It was also shown that TN and TP were more derived from agricultural production activities and were affected by human behavior. The driving factors of TN and TP were mainly rainfall runoff and soil erosion, which were affected by hydrological conditions and topographic factors. Therefore, the model has the ability to simulate TN and TP loadings and proved the spatial difference of NPS pollution in karst areas.

Numerous studies have demonstrated that the BMPs have positive effects on the reduction of NPS pollution (Ni & Parajuli 2018; Wang et al. 2018a, 2018b, 2020; Ricci et al. 2020). For the efficiency of FR, the reduction rate of TN and TP were positive and efficient while it was zero for sediment. This result was supported by other studies (Tian et al. 2010; Sun et al. 2013; Jiang 2016; Wang et al. 2016). Overall, the FR measure had a positive effect on reducing the pollutant of TN and TP, but the effect on reducing sediment was zero.

The eco-friendly measures of P factor only can reduce the particulate N and P, and it will increase dissolved nitrogen and phosphorus. These results were similar to other studies (Qiu et al. 2019; Lee et al. 2020). When the P factor was first implemented, the reduction rate of particulate nitrogen and phosphorus was greater than the increased rate of dissolved N and P, but with the passage of time, the reduction rate slightly decreased. In previous studies, it was found that this soil and water conservation tillage would produce paradoxes, which produced a tendency of trade-offs (Geng et al. 2015; Geng & Sharpley 2019).

The reduction efficiency of ST varied with time which indicated that the ST measure was also accompanied by other influencing factors, such as temperature and rainfall. Increasing the surface residual coverage rate was to change the surface crop residues to protect the soil from the loss of soil erosion caused by operations such as rainfall and crop irrigation, thereby reducing the load of TN and TP lost to the water body and protecting the farmland. However, Merriman et al. (2019) proposed that crop cover resulted in an increase in dissolved reactive phosphorus and TP, which were different from our study. This may be explained by the fact that nitrogen and phosphorus loss mechanisms in karst areas were distinct with non-karst areas.

LT was implemented in the basin which produced a small reduction effect on sediment and nutrients. The results were in good consistency with the researches by Huang et al. (2010) and Qiu et al. (2019). Protective tillage methods such as no-till and contour planting have reduced the output of sediment but increased the output of nutrients. A study similar to our results also showed that reduction rate of TN and TP was also very low, only 0.2 and 0.8% (Ding et al. 2020). This can be explained by reasons that smooth surface is potential to cause surface runoff to increase nutrient loss. In addition, limited soil disturbing activities results in more nutrient accumulating in soil layer of cultivation that further increase the amount of nutrient washed out by surface runoff (Qiu et al. 2019). However, Abdelwahab et al. (2014) found that the measure of no tillage can effectively decrease the load of soil erosion in a Mediterranean agricultural watershed in southern Italy. Ricci et al. (2020) also found that no tillage can reduce the sediment which was different from our study. This also demonstrated that less tillage and no tillage were not the BMPs in every watershed, and they need to be adapted to local conditions.

Selection of BMPs in the Qifeng River watershed should be based on the optimized measure reducing sediment and nutrient simultaneously. Although the reduction rate of the P factor was the largest, it decreased over time in the long run. However, the measure of FR cannot reduce the output of sediment but has a positive reduction effect on nutrients. LT has the smallest reduction effect on NPS pollution. ST was the BMPs for the Qifeng River Basin that can reduce NPS pollution more effectively. The BMPs simulation for long-term effect by AnnAGNPS model demonstrated that the reductions of NPS pollution depends on many factors, such as land use, climate, management of cultivation, and so on.

In this study, the AnnAGNPS model was implemented on the watershed to evaluate its suitability. The effectiveness of NPS pollution for four measured multi-scenarios in karst conditions were also estimated. AnnAGNPS model demonstrated ability in simulating NPS pollution after calibration in the Qifeng River watershed. The temporal distribution of TN and TP showed that TN and TP have the largest output during the rainy season. Rainfall runoff was the major influencing factor of NPS pollution output loads. Spatially, TN and TP were mainly distributed in the upper reaches of the sub-basin mainly occurred on orchards and woodlands. TN and TP were affected by hydrological conditions and topographical conditions, and human activities were also the main influencing factors. ST was selected as the BMP for the Qifeng River Basin. Long-term simulation results showed that different management measures have a different trend of reduction rate. Using the BMPs according to local conditions can control NPS pollution effectively.

Conceptualization, LY.L.; software, LY.L.; formal analysis, LY.L.; investigation, GS.P. and LY.L.; resources, HH.Z.; data curation, LY.L.; writing – original draft preparation, LY.L.; writing – review and editing, QL.T.; visualization, Z.L.; supervision, Z.L.; project administration, HH.Z.; funding acquisition, HH.Z.

This research was funded by National Natural Science Foundation of China, grant number 51578171.

The authors declare no conflict of interest.

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