Evaluation of Distribution Properties of Non-Point Source Pollution in a Subtropical Monsoon Watershed by a Hydrological Model with a Modified Runo ff Module

Non-point source pollution (NPS) is difficult to manage for watersheds, due to the scattering of pollution sources and uncertainty over the time it takes to accumulate. Since local agriculture and poultry rearing prevail, NPS occupies a large proportion of local pollution. In this paper, we modified the runoff module of the Soil and Water Assessment Tool (SWAT) to study the distribution properties and the effect of control of NPS in the Binjiang watershed in Southern China. The model was run from 2005 to 2014. The runoff simulation’s accuracy had apparently improved compared to the original model, as the Nash coefficient (Ens) had improved from 0.72 to 0.89, and the determination coefficient (R2) had improved from 0.75 to 0.91, an improvement in accuracy of 23.61% (Ens) and 21.33% (R2). Thus, the modified model (SWAT-m) is more adaptable for the simulation of extensive non-urban watersheds in subtropical monsoon climates. The validated model was used in further analysis. The results show that 82–90% of total nitrogen and 83–89% of total phosphorus were concentrated in the period from April to September annually. Through the introduction of the pollution generation potential parameter Φi, we obtained the descending order of all sub-basins in terms of their pollution generation potential. The critical source areas (CSAs) were found to be the northeast sub-basins in lower terrain that is used mainly for agricultural applications, accounting for 53% of the total watershed. The accumulation time is April to July, occupying 69% of annual generation. The simulation of management measures showed that NPS control has a good effect, with a 15 m filter strip in CSAs. Ammonia nitrogen and total phosphorus can be reduced apparently by 32% and 43%, respectively. The results may provide support for the management of NPS in watersheds under similar conditions.


Introduction
The Binjiang watershed, located in Qingyuan City, Guangdong Province, China, is an important source of local water supply for agriculture.It covers the area of latitude 23 • 44 30 N-24 • 16 04 N, longitude 112 • 32 40 E-113 • 4 50 E, as illustrated in Figure 1.Since local agriculture and poultry rearing prevail and there is abundant precipitation from the subtropical monsoon climate, non-point source pollution (NPS) occupies a large proportion of local water pollution.NPS is a form of pollutant that flows into underground or surface water from dispersed, wide-ranging, and traceable sources [1].

The Modification of the Hydrological Model
The distributed hydrological model Soil and Water Assessment Tool (SWAT) has a modular structure and consists of hydrology, sediment, and chemical subroutines applicable to watershed scales.The hybrid model is spatially based on hydrological response units that include both physical and conceptual methods [17].The hydrologic cycle as simulated by SWAT is based on the water balance equation [18].The in-stream kinetics used in SWAT for nutrient/pollutant routing are adapted from a widely used stream water quality model, QUAL2E [19].In this study, nitrogen and phosphorus are the main nutrients/pollutants.
Surface runoff is flow that occurs along a sloping surface.In the runoff module of SWAT, the surface runoff volume is computed using a modification of the Soil Conservation Service (SCS) curve number method [20].The equation for the SCS curve number method is as follows: R day -I a +S (1) where Qsurf is the accumulated runoff or rainfall excess (mm), Rday is the rainfall depth for the day (mm), Ia are the initial abstractions, which include surface storage, interception, and infiltration prior to runoff (mm), and S is the retention parameter (mm).The retention parameter varies spatially due to changes in soils, land use, management, and slope and temporally due to changes in soil water content.The retention parameter S is defined as: S=25.4 1000 CN -10 (2) The NPS model is an extended hydrological model that includes additional modules for assessing the release, migration, and leakage of pollutants from land sources, soil erosion, and sediment transport [2].Rapid advances in remote sensing and geographical information systems (GISs) during the past three decades have provided an enormous number of opportunities for hydrology research [3].The GIS-based distributed hydrological model Soil and Water Assessment Tool (SWAT) was developed by the USDA (United States Department of Agriculture), and is currently being utilized in several large basin projects, such as the HUMUS (Hydrologic Unit Model of the United States) project [4].This model is also used as a decision tool [5] to predict NPS pollution load [6] and for the assessment [7] of NPS pollution reduction effects of various control measures [8].Through the combined use [9] of land use changes [10] and climate changes [11], it is also used to predict the impact of future climate/land use scenarios [12] on hydrology [13] and pollution in basins [14].However, researchers have usually focused on the total outflow load or concentration per hectare of NPS rather than the NPS's properties, such as the generation potential.Nevertheless, the runoff module of the original SWAT has some limitations.Several researchers have found that the simplification of the surface runoff equation of SWAT does not comply with real runoff generation under certain conditions [15,16].Therefore, we modified the surface-runoff module of SWAT, and the modified SWAT-m model was applied to do this research.As for the NPS outflow load per unit area, it is difficult to tell whether a certain amount, say 1 t/ha annually, is high or low, as this depends on local pollution or different precipitation conditions.The only standard of a certain outflow load per hectare could be a violation of justice when deciding on the NPS critical source areas (CSAs).Thus, we introduced a new dimensionless parameter Φ i to measure the NPS generation potential of each sub-basin impartially.Through this, Water 2019, 11, 993 3 of 21 we can compare the NPS contribution of different sub-basins and pollutant components under various circumstances and, therefore, accurately position NPS CSAs.
Through the combination of remote sensing and GIS, we established this model to study the hydrology and NPS in a subtropical monsoon watershed.Further study of NPS, especially the distribution properties of NPS outflow and generation, the accurate positioning of NPS CSAs, and the impact of different filter strips on NPS pollutants, will contribute to NPS control and management in watersheds under similar conditions.

The Modification of the Hydrological Model
The distributed hydrological model Soil and Water Assessment Tool (SWAT) has a modular structure and consists of hydrology, sediment, and chemical subroutines applicable to watershed scales.The hybrid model is spatially based on hydrological response units that include both physical and conceptual methods [17].The hydrologic cycle as simulated by SWAT is based on the water balance equation [18].The in-stream kinetics used in SWAT for nutrient/pollutant routing are adapted from a widely used stream water quality model, QUAL2E [19].In this study, nitrogen and phosphorus are the main nutrients/pollutants.
Surface runoff is flow that occurs along a sloping surface.In the runoff module of SWAT, the surface runoff volume is computed using a modification of the Soil Conservation Service (SCS) curve number method [20].The equation for the SCS curve number method is as follows: where Q surf is the accumulated runoff or rainfall excess (mm), R day is the rainfall depth for the day (mm), I a are the initial abstractions, which include surface storage, interception, and infiltration prior to runoff (mm), and S is the retention parameter (mm).The retention parameter varies spatially due to changes in soils, land use, management, and slope and temporally due to changes in soil water content.The retention parameter S is defined as: where CN is the curve number for the day.The initial abstractions, I a , are commonly approximated as 0.2S in SWAT, which is the limitation of this module.In the original equation, the relation between I a and S is described by another parameter λ.Here, λ is fixed at 0.2, where we would modify the module.
Many researchers have proved that the fixed initial abstraction 0.2 does not match well with their study area.Researchers have taken measures to find out that the initial abstraction λ ranges from 0 to 0.3 [15,16,20].Woodward [21] also proved by his research that, in the experience of almost 90% of researchers, λ should be less than 0.2.In order to avoid deviation in the simulation caused by a single initial abstraction, this paper adopts the method of a segmentation value to improve the model, that is, a different initial abstraction λ is adopted according to the curve number (CN) value ranges of different regions.Based on data about local precipitation and hydrology events, we have obtained the initial abstraction values under different conditions through the event analysis and back calculation method.For instance, for an area covered by a natural forest with sandy-loam soil, the CN value is about 50-65, and the corresponding initial abstraction is relatively high: λ is about 0.12.For a concentrated residential area with clay soil, the CN values are relatively high, approximately 85-95, and the initial abstraction is relatively low: about 0.05.For cross-sectional conditions, such as clay soil farmland, the CN value is between 65 and 85, and the corresponding initial abstraction is around 0.08.According to the above classification method, we used Equation ( 1) to obtain the flow equation as follows: Through the modification of the surf.f,main.f, and simulate.ffunctions from the source code of SWAT in Visual Studio, we recompiled the modified code into the swatuser.exeprogram for further simulation and research.The modified version of SWAT is denoted SWAT-m.We compare the modified model with the original model in the next section.

Establishment of Databases
Spatial and attribute databases were established to build up the model.The databases include topography, land cover, soil classification, and weather data.
Firstly, topography data were extracted from a digital elevation model (DEM).The data used in this study is from the International Scientific Data Service Platform [22].We performed slope classification and hydrological analysis based on the DEM.
Secondly, we used as a source of land cover condition data the National Earth System Science Data Sharing Infrastructure [23].The classification of land use conditions in the study area is shown in Table 1.Thirdly, the soil classification and soil properties were extracted from the 1:500,000 scale FAO-UNESCO (Food and Agriculture Organization of the United Nations) Digital Soil Map of the new comprehensive Harmonized World Soil Database, as shown in Figure 2. Most of the soil properties that were obtained refer to "Guangdong soil" [24] and Soil Science, Chinese Academy of Science [25].Additionally, daily weather data from six stations in the watershed from 2005 to 2014 were obtained from the China Meteorological Science Data Sharing Service Network [26], which is shown in Figure 3.
Thirdly, the soil classification and soil properties were extracted from the 1:500,000 scale FAO-UNESCO (Food and Agriculture Organization of the United Nations) Digital Soil Map of the new comprehensive Harmonized World Soil Database, as shown in Figure 2. Most of the soil properties that were obtained refer to "Guangdong soil" [24] and Soil Science, Chinese Academy of Science [25].Additionally, daily weather data from six stations in the watershed from 2005 to 2014 were obtained from the China Meteorological Science Data Sharing Service Network [26], which is shown in Figure 3.The remaining data, such as point source, agricultural, and water supply volume data, were extracted from local chronicles and environmental survey reports from the local Environmental Protection Agency.According to the input data requirements, all geographical coordinates of spatial data need to be unified.In this study, the projection coordinate system used was Gauss Kruger Projection from the Beijing 54 Coordinate System.The remaining data, such as point source, agricultural, and water supply volume data, were extracted from local chronicles and environmental survey reports from the local Environmental Protection Agency.According to the input data requirements, all geographical coordinates of spatial data need to be unified.In this study, the projection coordinate system used was Gauss Kruger Projection from the Beijing 54 Coordinate System.

Delineation of Sub-Basins and the Hydrologic Response Unit
To simulate runoff in a watershed with a river network, we should first delineate its sub-basins.The model obtains the flow direction of each grid by the maximum gradient method, with the D8 algorithm, on the DEM, and then counts the number of flows through the grid to estimate the cumulative flow [27].Through the process of water delineation, the study area is divided into 31 sub-basins, which are shown in Figure 3.
Based on the division of the study area into sub-basins, the model applies an algebra calculation to the land use layer, the soil classification, and the slope classification to generate the smallest hydrological response unit (HRU).The slope classification refers to cultivated land gradients: less than 15 degrees; 15-25 degrees; and greater than 25 degrees.According to the local conditions, the multiple HRU threshold of land use, slope gradient, and soil classification was set to 20%, 20%, and 10%, respectively.After the raster calculation, 749 HRUs were generated in the study region.Runoff was predicted separately for each HRU and routed to obtain the total runoff for the watershed.Now, we are ready to run the model.

Calibration and Validation
In this study, the model was calibrated by the SUFI-2 algorithm method through the model SWAT-CUP.The measured data from 2005 to 2009 were used to calibrate the parameters, and the data from 2010 to 2014 were used to validate the model.To evaluate the simulation's accuracy, we applied three indicators: determination coefficient (R 2 ), relative error (Re), and Nash coefficient (Ens) [28,29].

The NPS Generation Potential Parameter
Because the specific value of the outflow load per hectare of NPS is impacted by various kinds of pollutants, the overall pollution levels, and the precipitation conditions in different regions, the only standard of a certain load per hectare could be a violation of justice in an NPS comparison between different districts and various pollutants and when deciding on NPS CSAs.For instance, with the standard of 10kg/ha for NPS CSAs, for some upstream watersheds, maybe only one sub-basin is identified as a CSA, yet for some downstream populated watersheds, maybe all sub-basins are over 10kg/ha and thus constitute CSAs under this criteria, which is improper for CSAs and NPS management.In order to analyze and compare the generation potential of NPS among different regions and facilitate the accurate positioning of NPS CSAs, we introduce a new dimensionless parameter to measure the generation potential of NPS in each sub-basin.The equation is as follows: where Φ i is the generation potential parameter of sub-basin No. i of a certain kind of NPS pollution; i stands for the number of the sub-basin, n is the total number of sub-basins in the study area; C i stands for the average pollutant yield per unit area in the time step, kg/ha monthly or annually; A i stands for the area of sub-basin No. i; and A stands for the total area of all sub-basins.Φ i = 1 is the average level of NPS generation, which means that all sub-basins have the same generation potential of a certain pollutant.The lower the contribution is, the closer the value of Φ i is to zero; conversely, the higher the Water 2019, 11, 993 7 of 21 generation potential is, the higher the value of Φ i is; the maximum is less than n.For different kinds of pollutants, the value of Φ i can be added to show the multiple contributions of a certain sub-basin.

The Theory of Filter Strips for NPS Control
After obtaining the distribution properties of NPS CSAs, this section focuses on pollution control in those CSAs, mainly sub-basins 1-11 upstream in the watershed.Among the control measures in the Best Management Practices (BMPs), such as Terracing, Contouring, Filter Strips, and Artificial Wetlands [7], we chose filter strips, as the NPS control effect is ideal and the construction and operation cost is relatively low, especially in rural areas [8].
A filter strip is a plant-intensive strip that intercepts and filters runoff from an uphill or upstream area.It can decrease the flow rate and increase the infiltration of the area to reduce runoff and pollutants.The filter strip algorithm that SWAT uses is derived from White and Arnold's research methods for vegetation filter strips, which can effectively reduce sediment, pollutants, bacteria, and pesticide components, but do not affect runoff [18,30].A filter strip's trapping of sediments and nutrients in the SWAT management module is calculated as: where Trap ef and Trap ef,sub are the surface and subsurface flow constituent loadings that the filter strip traps, respectively, and Width filstrip is the width of the filter strip.The proper parameter of the filter strip's width can be found in FILTERW.mgt in the management module of SWAT.Based on the hydro-geological conditions, we present a better design of a riparian filter strip through a scenario simulation of a filter strip with a width of 3-30 m, and provide a reference for NPS control in similar clay-loam areas in a subtropical monsoon climate.

Parameter Sensitivity Analysis
In this study, the model was calibrated using the SUFI-2 algorithm method through the model SWAT-CUP.The first 12 sensitive parameters of runoff were obtained through a sensitivity analysis and are shown in Table 2. p-values were used to determine the significance of the sensitivity (values close to zero have more significance).t-Stat provides a measure of sensitivity (larger absolute values are more sensitive) [31].
As for the pollutants, the sensitivity analysis showed that five parameters (nitrate percolation coefficient, benthic source rate for NH 4 -N, rate coefficient for organic N settling in the reach at 20 • C, rate constant for biological oxidation of NH 4 to NO 2 , and rate constant for hydrolysis of organic N to NH 4 ) have a great influence on ammonia nitrogen.The phosphorus absorption coefficient, soil phosphorus separation coefficient, phosphorus permeability coefficient, residual mineralization rate, and rate constant for organic phosphorus transferred to dissolved phosphorus at 20 • C have a great influence on total phosphorus.These parameters were applied in the calibration process.

The Comparison of SWAT and SWAT-m
In order to prove that the improved model has better adaptability, the original SWAT was applied to simulate runoff in the same area as a blank control.The input data were all the same as the inputs into the modified model (SWAT-m).The best parameters of the calibration period (2005)(2006)(2007)(2008)(2009) were applied to the validation period (2010-2014) for verification.The comparison of the runoff simulation accuracy of SWAT and SWAT-m based on three indicators (determination coefficient (R 2 ), relative error (Re), and Nash coefficient (Ens)), is depicted in Table 3.It can be seen from Table 3 that the R 2 of SWAT is 0.82 and 0.75 for the calibration period and the validation period, respectively, which is greater than the standard 0.60 [29].The Nash coefficient (Ens) of the original SWAT is 0.76 and 0.72 for the calibration period and the validation period, respectively, which does not reach the standard of 0.75 [28].These results indicate that the simulation results are not that ideal.Apparently, the simulation accuracy is not as good as SWAT-m, especially for the validation period.From the comparison of the Nash coefficients, the simulation accuracy of SWAT-m has improved by 23.61%.The improvement in simulation accuracy is more obvious in the comparison of relative error, which is as high as 27% in the validation period, and the simulation accuracy of SWAT-m has improved by 48.15%.Therefore, if the original model is used in further NPS simulation, the relative error generated by the runoff simulation will further increase; the principle is similar to the butterfly effect.In summary, the improved model (SWAT-m) has higher runoff simulation accuracy; thus, we will carry out NPS simulation using the modified model.

Results of Simulations with the Modified Model
The model was calibrated using observed values from 2005 to 2009 and with an optimal set of parameters.Results of Re < 10%, R 2 > 93%, and Ens > 91% were obtained.Then, we input the best parameters and reran the model over the validation period (2010-2014).The results of the simulation well-agreed with the observed values, with Re < 15%, R 2 > 89%, and Ens > 90%.The comparison of the simulated and measured values of runoff is shown in Figure 4a.As for pollutants, after calibration and validation of the runoff simulation, SWAT-m was applied to a pollutants simulation in which the calibration and validation time was the same as in the runoff simulation.The comparison of the simulated and measured values of the monthly average ammonia nitrogen and total phosphorus load is shown in Figure 4b,c.The best calibration and validation results for ammonia nitrogen were Re = 38%, R 2 = 0.85, and Ens = 0.71 (calibration) and Re = 39%, R 2 = 0.77, and Ens = 0.69 (validation).For total phosphorus, the best calibration and validation results were Re = 30%, R 2 = 0.78, and Ens = 0.79 (calibration) and Re = 37%, R 2 = 0.75, and Ens = 0.72 (validation).Generally, if the absolute value of Re is less than 15%, the absolute value of R 2 is greater than 0.6, and the absolute value of Ens is greater than 0.75, the simulation result is sufficiently accurate [28].The standard for nutrient assessment is not as high.The absolute value of relative error should be less than 40%, the absolute value of the Nash coefficient (Ens) should be greater than 0.5 [31], and the standard for R 2 is the same as that for runoff [32].We can tell that the simulations for both nitrogen and phosphorus are able to comply with the standard of accuracy, indicating that the modified model is effective for further research in NPS simulation and analysis.

Temporal Distribution
Based on the simulation, we summarized the annual runoff and pollutants discharge.The temporal distribution of ammonia nitrogen (NH 4 -N), total nitrogen (TN), and total phosphorus (TP) in different hydrological years is shown in Figure 5. From the figure, we can see that the pollution load from April to September (rainy season) accounted for a large proportion within a year, with the highest average in May and June.The contribution of NH 4 -N, TN, and TP load during the flood season to the annual load is 81-87%, 82-90%, and 83-89%, with an average of 84%, 85%, and 86%, respectively.The change between different hydrological years was obvious, especially the variation in NPS temporal distribution over 12 months.For instance, the coefficient of variation in the high-flow year, 0.78-0.85exactly, is lower than that in the low-flow year, 0.95-0.99exactly.Correspondingly, the difference between the maximum and minimum value in the low-flow year is the largest.with abundant rainfall; the correlation coefficient (R) of runoff and rainfall over the year is 0.87.Since runoff generated from rainfall is a driver of NPS, the temporal distribution of NPS is affected by the central precipitation time.The average annual rainfall in the watershed is mainly concentrated in the period from April to September, accounting for about 80% of annual precipitation.Therefore, the temporal distribution of NPS also presents the same regularity.It approximately belongs to a normal distribution, and accumulates from April to September.This temporal distribution property of NPS is consistent with results from studies of similar watersheds [33], indicating the close relationship between the temporal distribution of NPS and local precipitation.The reason for this could be attributed to the relationship between precipitation and runoff, which is illustrated in Figure 6.It can be seen from Figure 6 that the runoff distribution shows a significant correlation with rainfall among the time series, and runoff also accumulates in months with abundant rainfall; the correlation coefficient (R) of runoff and rainfall over the year is 0.87.Since runoff generated from rainfall is a driver of NPS, the temporal distribution of NPS is affected by the central precipitation time.The average annual rainfall in the watershed is mainly concentrated in the period from April to September, accounting for about 80% of annual precipitation.Therefore, the temporal distribution of NPS also presents the same regularity.It approximately belongs to a normal distribution, and accumulates from April to September.This temporal distribution property of NPS is consistent with results from studies of similar watersheds [33], indicating the close relationship between the temporal distribution of NPS and local precipitation.As for the reason, since the stream is almost the lowest in topography, it is the destination at which surface runoff converges.Therefore, the NPS from upstream and remote sub-basins tends to converge in sub-basins alongside the stream, which has relatively low-altitude terrain.This is why the total load of NPS accumulates in the sub-basins alongside the main stream, especially in the southern area.

Analysis of the Generation Potential of NPS
The NPS load in sub-basin No. 31, for instance, does not generate itself, but accumulates from upstream sub-basins.Here, the question arises as to where the CSAs are; that is, which sub-basins generated a high level of NPS.

NPS Generation with Topography
Here, we analyze the generation distribution of NPS.It is focused on a sub-basin's own pollution yield, and not the total outflow load that could be derived from other sub-basins.Researchers have usually focused on land cover with the spatial distribution of NPS, and studied the impact of land use change on NPS distribution [10][11][12].Considering the gentle rolling/hilling topography and unchanged land cover of the watershed in the study area, we focused on NPS generation with topography.The main components of TN and TP generation under the topography conditions are depicted in Figure 8. spatial distribution of NPS load is not evenly distributed, and mainly accumulates in the central and southern areas alongside the main stream.
As for the reason, since the stream is almost the lowest in topography, it is the destination at which surface runoff converges.Therefore, the NPS from upstream and remote sub-basins tends to converge in sub-basins alongside the stream, which has relatively low-altitude terrain.This is why the total load of NPS accumulates in the sub-basins alongside the main stream, especially in the southern area.

Analysis of the Generation Potential of NPS
The NPS load in sub-basin No. 31, for instance, does not generate itself, but accumulates from upstream sub-basins.Here, the question arises as to where the CSAs are; that is, which sub-basins generated a high level of NPS.

NPS Generation with Topography
Here, we analyze the generation distribution of NPS.It is focused on a sub-basin's own pollution yield, and not the total outflow load that could be derived from other sub-basins.Researchers have usually focused on land cover with the spatial distribution of NPS, and studied the impact of land use change on NPS distribution [10][11][12].Considering the gentle rolling/hilling topography and unchanged land cover of the watershed in the study area, we focused on NPS generation with topography.The main components of TN and TP generation under the topography conditions are depicted in Figure 8.It can be seen that the spatial distribution of TN and TP generation shares some similarities.For example, the topography value and NPS generation have a negative correlation.In other words, the area where the terrain is low has a relatively high level of NPS yield.Another similarity is that, among the main components of TN and TP, the organic component, ORGN (Organic Nitrogen) and ORGP (Organic Phosphorus), is the dominant component.It is shown in Figure 8 that, regardless of whether the total height of the bar is high or not, the organic component occupies most of the total height: around 60-80%.

Accurate Positioning of NPS CSAs
In order to analyze the generation potential of NPS more thoroughly, the generation potential (Φi) of the 31 sub-basins was calculated and is shown in Table 4.In the table, Φi > 1 indicates that the generation potential of this sub-basin is higher than average; considering the multiple influences of nitrogen and phosphorus, the value of Φi > 2 in the "TN + TP" column means the same.We arranged the 31 sub-basins based on a descending contribution to NPS generation; high generators are It can be seen that the spatial distribution of TN and TP generation shares some similarities.For example, the topography value and NPS generation have a negative correlation.In other words, the area where the terrain is low has a relatively high level of NPS yield.Another similarity is that, among the main components of TN and TP, the organic component, ORGN (Organic Nitrogen) and ORGP (Organic Phosphorus), is the dominant component.It is shown in Figure 8 that, regardless of whether the total height of the bar is high or not, the organic component occupies most of the total height: around 60-80%.

Accurate Positioning of NPS CSAs
In order to analyze the generation potential of NPS more thoroughly, the generation potential (Φ i ) of the 31 sub-basins was calculated and is shown in Table 4.In the table, Φ i > 1 indicates that the generation potential of this sub-basin is higher than average; considering the multiple influences of nitrogen and phosphorus, the value of Φ i > 2 in the "TN + TP" column means the same.We arranged the 31 sub-basins based on a descending contribution to NPS generation; high generators are highlighted in gray.Table 4 shows that the sub-basins with a high level of generation potential are Nos.2-7, 9, 11, 22, and 23, for which the value of Φ i is over 3.0, which means that the NPS generation potential is over 1.5 times the average level.To answer the question of where the CSAs are, a large proportion is from the northeast sub-basins with lower terrain that is mainly used for agricultural applications, accounting for 53% of the total NPS generation of the watershed.The value of Φ 31 in the "TN + TP" column is less than 2.0; therefore, the generation potential at sub-basin No. 31 is less than the average level.
Additionally, we classified the sub-basins into three ranks (relatively low, medium, relatively high) based on Φ i , the generation potential of "TN + TP".In the first rank, Φ i is less than 1.6.This rank generates little pollution; examples are sub-basin Nos. 15, 17, and 26.The terrain is relatively high, and the land cover is mainly thick evergreen forest.Although precipitation is relatively high, the densely distributed vegetation alleviates pollution loss by rainfall.In the second rank, Φ i is between 1.6 and 3.0.This rank generates a medium level of NPS; examples are sub-basin Nos. 18, 19, and 20.The terrain is relatively low.The land cover contains forest, grassland, and small amount of agriculture.There is a medium amount of precipitation.In the third rank, Φ i is greater than 3.0.This rank generates a relatively high level of pollution; examples are sub-basin Nos. 2, 6, and 9.The terrain is low and the precipitation is relatively high.The land cover is mainly farmland and residential.Without densely distributed forest, the last rank generates a large amount of NPS pollution by human activities and rainfall scouring.This is consistent with the results of studies on the impact of land cover change on NPS pollution [12] that show that an increase in farmland or residential areas and a reduction in green vegetation would increase the NPS pollution to a large extent [10][11][12][13][14].These are the CSAs that demand special attention.The generation potential of "TN + TP" and land cover is depicted in Figure 9.

The Main NPS Pollutants in Sub-Basins
Compared with the load per hectare, which is applied widely in NPS studies [4][5][6], the introduction of the NPS generation potential parameter Φi has several advantages.For a certain NPS pollutant or a combination of several pollutants, the generation potential parameter Φi can tell us which are the high generators (CSAs) and the severity of their pollution.As for load per hectare, for instance, researchers have set 5 kg/ha TN as the standard of NPS CSAs when studying the head stream watershed of the Yellow River in China [34].However, this standard could not be applied to

The Main NPS Pollutants in Sub-Basins
Compared with the load per hectare, which is applied widely in NPS studies [4][5][6], the introduction of the NPS generation potential parameter Φ i has several advantages.For a certain NPS pollutant or a Water 2019, 11, 993 15 of 21 combination of several pollutants, the generation potential parameter Φ i can tell us which are the high generators (CSAs) and the severity of their pollution.As for load per hectare, for instance, researchers have set 5 kg/ha TN as the standard of NPS CSAs when studying the head stream watershed of the Yellow River in China [34].However, this standard could not be applied to other watersheds or other pollutants under different pollution conditions.Nevertheless, the introduction of the dimensionless parameter Φ i can help us to make progress.The value of Φ i > 3 in "TN + TP" or Φ i > 1.5 in other pollutants designates the high NPS generators or CSAs, regardless of the background or condition of the watershed.Moreover, for a certain district, the comparison of Φ i of different pollutants can provide information on the main local pollutants.For instance, in Table 4, for sub-basin No. 18, Φ ORGN = 1.339,Φ NH4 = 2.642, and Φ NSURQ = 0.977, which means that the organic nitrogen and NO 3 in the surface runoff are average, while the ammonia nitrogen is over 2 times the average level.Thus, in sub-basin No. 18, the main pollutant is ammonia nitrogen, which demands that focus be placed on relevant control measures.In summary, the sub-basins for which the main pollutant is ammonia nitrogen are Nos. 1, 13, 14, 18, 20, 22-24, 28, and 31, which are located in the central and southern districts alongside the main stream, and account for 26.67% of the whole watershed.The sub-basins for which the main pollutants are sediment and organic nitrogen are Nos.2-7, 9, and 10, which are mainly located in the northern districts for farming or agricultural usage, and account for 16.43% of the total area of the watershed.The sub-basins for which the main pollutant is biological oxygen demand (BOD) are Nos.12, 13, 16, 17, 19, 21, and 25-27, which are mainly located in low-NPS-generation areas with natural vegetation, and account for 35.47% of the total area of the Binjiang watershed.

Temporal Distribution of NPS CSAs
We also analyzed the generation potential of a combination of TN and TP in each sub-basin during a 12-month period, as depicted in Figure 10.As for high generators, the temporal distribution of NPS tends to accumulate from April to July, accounting for 69% of the annual yield, which will assist the temporal NPS control in certain areas.

Temporal Distribution of NPS CSAs
We also analyzed the generation potential of a combination of TN and TP in each sub-basin during a 12-month period, as depicted in Figure 10.As for high generators, the temporal distribution of NPS tends to accumulate from April to July, accounting for 69% of the annual yield, which will assist the temporal NPS control in certain areas.As shown in the figure, the x-axis represents the 31 sub-basins, the y-axis represents the 12 months in a year, and the z-axis represents the generation potential of NPS.It can be seen that the distribution of NPS generation potential in the study area is not balanced, whether on the time scale or the spatial scale.Firstly, the x-axis of the geographical space is from right to left.It can be seen that Water 2019, 11,993 As shown in the figure, the x-axis represents the 31 sub-basins, the y-axis represents the 12 months in a year, and the z-axis represents the generation potential of NPS.It can be seen that the distribution of NPS generation potential in the study area is not balanced, whether on the time scale or the spatial scale.Firstly, the x-axis of the geographical space is from right to left.It can be seen that the upstream area (sub-basin Nos.1-11) has a relatively high potential for NPS generation, followed by the middle districts (sub-basin Nos.22 and 23).As for the central margin and southern part of the watershed, only sub-basin No. 28 has a slightly higher NPS generation potential; the overall NPS generation potential throughout the 12-month period is relatively low.Moreover, from January to December, it can be seen that the highest NPS generation occurs in May and June, and is concentrated mainly from April to July.In addition, for some CSAs (sub-basin Nos.2-6 and 9), the NPS generation potential in February and March is also relatively high.Researchers studying the NPS in farmland in southern China obtained similar results, except for the rainy season.February and March also have high NPS load [33], which could be attributed to the tillage and fertilization of farmland in this period.Therefore, during this time, we should focus on these particular areas, especially those used for farming.

Scenario
In this study, a riparian filter strip was set for the hydrological response unit around the main stream and tributary of the Binjiang river, mainly in the NPS CSAs (sub-basin Nos.[1][2][3][4][5][6][7][8][9][10][11].The river bank's average slope ratio is 1:8.The vegetation is shrub herb in natural grassland.The filter strip's width was designed to be 3-30 m.The riparian buffer strips were added to the management module of SWAT, and their reduction effects on different NPS pollutants were simulated.

Simulation and Best Selection
Based on the above scenario, the improved SWAT model was used to simulate the generation of NPS pollutants with different filter strips.The evaluated pollutants were SED (sediment or suspended particulate), BOD (biochemical oxygen demand), NH 4 (ammonia nitrogen), TN (total nitrogen), and TP (total phosphorus).With the original condition as a blank scenario, we simulated the reduction effects of the riparian filter strips on various pollutants with different structural designs, depicted in Figure 11.
The NPS generation and reduction rate are listed in Table 5.
Water 2019, 11, x FOR PEER REVIEW 18 of 21 to a reduction in nitrogen and BOD loads between 2.41% and 5.75% [37].The reason for this is that inorganic nitrogen migrates with surface runoff in a dissolved state.Organic and mineral phosphorus migrate and adhere to sediment in surface runoff.According to the literature, more than 70% of phosphorus and about 20% of nitrogen (mainly organic nitrogen) accumulate in this way in catchment area [38].Since filter strips have a better sedimentation effect on SED, the TP adsorbed to particulate matter is also removed.Ammonia nitrogen, nitrate nitrogen, and BOD are mainly removed by absorption of plant roots and the degradation of microorganisms.Because of the slope of the river bank, the hydraulic retention time of pollution flowing through a riparian filter strip is insufficient for absorption and degradation; therefore, the reduction effect of TN and BOD is lower.Similarly, Stagge [39] demonstrated the performance of grass swales with experimental methods to measure pollutant concentrations in inflow and outflow, and found that grass swales reduced sediment and TP by 40-85% and 20-45%, respectively, while the reduction in TN and BOD was less than 15%.From the perspective of structure design, there is a certain nonlinear relationship between the design scenarios and the efficiency of the reduction in the pollutants.As shown in the figure, when the width of the buffer strip is less than 15 meters, the reduction effect of NPS pollutants gradually increases; that is, there is a significant positive correlation.When the width of the filter strip is greater than 15 meters, the control effect on NPS pollutants tends towards stability.Taking NH4 and TP as examples, the relationship between average reduction rate (y) and riparian buffer strip width (x) is similar to a logarithmic equation.Based on the data from Table 5, SPSS software was used to solve  As shown above, from the perspective of pollutants, the filter strip has the best control effect on SED; the average annual reduction rate reaches 41-50%.The reduction rate of the average annual output of TP reaches 37-45%.Researchers have studied the effectiveness of grass strips on water quality, and the results showed that the effects of grass strips on the removal of sediments and nutrients were significant: over 20% in sediments and TP [35,36].However, the reduction effect on ammonia nitrogen and TN is relatively low.The average annual reduction rate of TN, for instance, is only 10-22%.Similarly, the implementation of a 25% grass strip in Central Indiana, United States, led to a reduction in nitrogen and BOD loads between 2.41% and 5.75% [37].The reason for this is that inorganic nitrogen migrates with surface runoff in a dissolved state.Organic and mineral phosphorus migrate and adhere to sediment in surface runoff.According to the literature, more than 70% of phosphorus and about 20% of nitrogen (mainly organic nitrogen) accumulate in this way in catchment area [38].Since filter strips have a better sedimentation effect on SED, the TP adsorbed to particulate matter is also removed.Ammonia nitrogen, nitrate nitrogen, and BOD are mainly removed by absorption of plant roots and the degradation of microorganisms.Because of the slope of the river bank, the hydraulic retention time of pollution flowing through a riparian filter strip is insufficient for absorption and degradation; therefore, the reduction effect of TN and BOD is lower.Similarly, Stagge [39] demonstrated the performance of grass swales with experimental methods to measure pollutant concentrations in inflow and outflow, and found that grass swales reduced sediment and TP by 40-85% and 20-45%, respectively, while the reduction in TN and BOD was less than 15%.
From the perspective of structure design, there is a certain nonlinear relationship between the design scenarios and the efficiency of the reduction in the pollutants.As shown in the figure, when the width of the buffer strip is less than 15 meters, the reduction effect of NPS pollutants gradually increases; that is, there is a significant positive correlation.When the width of the filter strip is greater than 15 meters, the control effect on NPS pollutants tends towards stability.Taking NH 4 and TP as examples, the relationship between average reduction rate (y) and riparian buffer strip width (x) is similar to a logarithmic equation.Based on the data from Table 5, SPSS software was used to solve the fitting regression equation y = a•ln(x) + b.The solved nonlinear fitting regression equation is as follows: y NH4 = 0.0694ln(x) + 0.1157 (8) y TP = 0.0827ln(x) + 0.1876.
The fitting degree of this equation is fairly high.The relative error (Re) of NH 4 is only 6.81%, and the determination coefficient (R 2 ) of NH 4 is 0.9804.The relative error (Re) of TP is only 4.16%, and the R 2 of TP is as high as 0.9882.For areas with similar natural conditions (a subtropical monsoon climate, clay-loam soil, and non-urban areas), the fitting equation can be used to estimate the NPS reduction effect of riparian filter strips in the absence of numerical data for simulation.Researchers have studied the width of a filter strip and its pollution reduction effect in large flat watersheds.The results showned that the 15 m filter strip is slightly better than the 10 m filter strip, which is apparently better than the 5 m filter strip.In general, a filter strip that is 10-15 m in width can intercept 3-45% of TN, and 30-80% of TP [36][37][38][39][40]. Similarly, considering the reduction effect and cost, a 15 m filter strip alongside the main stream is suggested, as it can produce a reduction effect on sediment (SED), total nitrogen (TN), and total phosphorus (TP) of up to 48.12%, 19.04%, and 42.33%, respectively.

Conclusions
Because of the limitations to the surface runoff module in SWAT, in order to avoid deviation in simulations caused by a single initial abstraction value, this paper adopted the method of a segmentation value of initial abstraction to improve the model, named SWAT-m.By data collection, we established databases for topography, soil, land use, and weather.Then, we performed runoff and NPS pollutant simulations.The results enabled us to draw the following conclusions: Through the comparison of runoff simulations using SWAT and SWAT-m, we obtained that the modified model has better simulation accuracy.In particular, during the validation period, the evaluation indicators are: Ens = 0.72, R 2 = 0.75, and Re = 27% (SWAT); Ens = 0.89, R 2 = 0.91, and Re = 14% (SWAT-m).The modified model improves the simulation accuracy by 23.61% (Ens) and 21.33% (R 2 ), indicating its better adaptability to simulate extensive watershed pollution under various natural conditions.
As for NPS pollutants, simulation by SWAT-m also complies with the standard for accuracy.Thus, it was applied in further research.For a watershed in a subtropical monsoon climate, the temporal distribution of NPS approximately belongs to a normal distribution, and accumulates from April to September.The contribution of TN and TP in the rainy season to the annual load is 82-90% and 83-89%, with an average of 85% and 86%, respectively.Additionally, the average annual NPS load is not evenly distributed, and mainly accumulates in the central and southern areas alongside the main stream.
As for the spatial distribution of NPS, the topography value and NPS generation have a negative correlation.Through the generation potential parameter Φ i , we obtained the descending order of all sub-basins in terms of their contribution to total NPS generation.The CSAs were found to be mostly northeast sub-basins with lower terrain that is used mainly for agricultural applications, accounting for about 53% of the total generation.The temporal distribution of high NPS generators accumulated from April to July, reaching 69% of the annual load.Also, through a comparison of Φ i between pollutants, we found that the main pollutants of CSAs are organic nitrogen (ORGN), ammonia nitrogen (NH4), and particulate matter (SED); the main pollutant for the central and southern areas alongside the main stream is ammonia nitrogen.
By simulating the impact of a riparian filter strip on NPS pollution in CSAs, it was shown that filter strips efficiently reduce particulate matter (SED) and total phosphorus (TP) by up to 48.12% and 42.33%, respectively.The reduction effect on pollutants and the designed width of the filter strip indicated an approximate logarithmic relationship, which could be applied to similar watersheds without numerical data for simulation.When the buffer width reaches 15 meters, the reduction rate tends towards stability; thus, the recommended design is a 15-m wide filter strip, with a slope ratio of 1:8, that is covered with local natural herbaceous vegetation alongside the river.
Because of limitations in time, researchers, and measured data, we only validated two kinds of pollutants (ammonia nitrogen and total phosphorus).Further research may simulate the other main components of NPS, such as BOD (biochemical oxygen demand).
As a result of the research in this study, we can accurately position the high generators of NPS pollution and its accumulation time.Further research could focus on control measures other than filter strips to decrease NPS.The results of this study are of practical significance to the management of NPS pollution in watersheds under similar conditions.

Figure 1 .
Figure 1.The location of the research area.The upper half is the location of China, the lower half shows the contours of Guangdong Province, and yellow area in the middle depicts the research area.

Figure 1 .
Figure 1.The location of the research area.The upper half is the location of China, the lower half shows the contours of Guangdong Province, and yellow area in the middle depicts the research area.

Figure 2 .
Figure 2. The soil classification of the study area.Streams are depicted by blue lines.

Figure 2 . 21 Figure 3 .
Figure 2. The soil classification of the study area.Streams are depicted by blue lines.Water 2019, 11, x FOR PEER REVIEW 6 of 21

Figure 3 .
Figure 3.The division of the study area into sub-basins and the hydrologic stations in the study area.The topography is described by the gradient color of grass green.Streams are shown by blue lines.Hydrologic stations are denoted by yellow points.Weather stations are depicted by orange pentagons.Sub-basins are numbered in dark blue.
Water 2019, 11, x FOR PEER REVIEW 10 of 21

Figure 4 .
Figure 4.The calibration and validation of (a) runoff, (b) ammonia nitrogen, and (c) total phosphorus.Black points indicate observed values.Blue lines indicate simulated values.

3. 2 .
Temporal and Spatial Distribution of NPS 3.2.1.Temporal Distribution Based on the simulation, we summarized the annual runoff and pollutants discharge.The temporal distribution of ammonia nitrogen (NH4-N), total nitrogen (TN), and total phosphorus (TP) in different hydrological years is shown in Figure 5. From the figure, we can see that the pollution

Figure 4 .
Figure 4.The calibration and validation of (a) runoff, (b) ammonia nitrogen, and (c) total phosphorus.Black points indicate observed values.Blue lines indicate simulated values.

Figure 5 .
Figure 5.The percentage of monthly load of (a) ammonia nitrogen, (b) total nitrogen, and (c) total phosphorus in different hydrological years.

Figure 5 .
Figure 5.The percentage of monthly load of (a) ammonia nitrogen, (b) total nitrogen, and (c) total phosphorus in different hydrological years.

Figure 5 .
Figure 5.The percentage of monthly load of (a) ammonia nitrogen, (b) total nitrogen, and (c) total phosphorus in different hydrological years.

Figure 6 .
Figure 6.The correlation between precipitation and runoff during the simulation period.Figure 6.The correlation between precipitation and runoff during the simulation period.

Figure 6 .
Figure 6.The correlation between precipitation and runoff during the simulation period.Figure 6.The correlation between precipitation and runoff during the simulation period.

3. 2 . 2 .
Spatial Distribution of NPSThe average annual load of NH 4 -N and TP in the sub-basins is shown in Figure7.Apparently, the area alongside the main stream has a high level of NPS load; for example, sub-basin Nos. 18, 20, and 31 have an average annual outflow load of around 160 tons of NH 4 -N and 250 tons of TP.The northern part is at a medium level; the average annual outflow load is about 108 tons of NH 4 -N and 165 tons of TP.The eastern area generates the least outflow load; for example, sub-basin Nos. 10 and 16 have an average annual load that is less than 30 tons of NH 4 -N and 50 tons of TP.To conclude, the spatial distribution of NPS load is not evenly distributed, and mainly accumulates in the central and southern areas alongside the main stream.

Figure 7 .
Figure 7.The spatial distribution of the average annual load of ammonia nitrogen (a) and total phosphorus (b).Sub-basins are numbered in dark blue with gray boundary lines.

Figure 7 .
Figure 7.The spatial distribution of the average annual load of ammonia nitrogen (a) and total phosphorus (b).Sub-basins are numbered in dark blue with gray boundary lines.

Figure 8 .
Figure 8.The spatial distribution of Total Nitrogen (TN) (a) and Total Phosphorus (TP) (b) generation.The bars indicate the average yield of TN or TP per unit area.In the left panel, blue bars denote Organic Nitrogen (ORGN), yellow denotes Ammonia Nitrogen (NH4), and rose red denotes NO3 in surface runoff (NSURQ).In the right panel, green denotes Organic Phosphorus (ORGP), purple denotes Soluble Phosphorus (SOLP), and red denotes Mineral Phosphorus (SEDP).Sub-basins are numbered in red color with gray boundary lines.

Figure 8 .
Figure 8.The spatial distribution of Total Nitrogen (TN) (a) and Total Phosphorus (TP) (b) generation.The bars indicate the average yield of TN or TP per unit area.In the left panel, blue bars denote Organic Nitrogen (ORGN), yellow denotes Ammonia Nitrogen (NH 4 ), and rose red denotes NO 3 in surface runoff (NSURQ).In the right panel, green denotes Organic Phosphorus (ORGP), purple denotes Soluble Phosphorus (SOLP), and red denotes Mineral Phosphorus (SEDP).Sub-basins are numbered in red color with gray boundary lines.
Water 2019, 11, x FOR PEER REVIEW 15 of 21

Figure 9 .
Figure 9.The spatial distribution of non-point source pollution (NPS) generation potential.The bars stand for the NPS generation potential parameter Φi of each sub-basin.Blue denotes TN, and purple denotes TP.

Figure 9 .
Figure 9.The spatial distribution of non-point source pollution (NPS) generation potential.The bars stand for the NPS generation potential parameter Φ i of each sub-basin.Blue denotes TN, and purple denotes TP.
Water 2019, 11, x FOR PEER REVIEW 16 of 21 is biological oxygen demand (BOD) are Nos.12, 13, 16, 17, 19, 21, and 25-27, which are mainly located in low-NPS-generation areas with natural vegetation, and account for 35.47% of the total area of the Binjiang watershed.

Figure 10 .
Figure 10.The generation potential of "TN + TP" in the sub-basins during a 12-month period.

Figure 10 .
Figure 10.The generation potential of "TN + TP" in the sub-basins during a 12-month period.

Figure 11 .
Figure 11.The design of the riparian filter strip and its impact on pollution reduction.

Figure 11 .
Figure 11.The design of the riparian filter strip and its impact on pollution reduction.

Table 1 .
The classification of land use conditions of the study area.

Table 2 .
The sensitive parameters for the runoff simulation.

Table 3 .
The comparison of SWAT and SWAT-m on the runoff simulation accuracy.

Table 4 .
The value of the generation potential parameter Φ i of all sub-basins for different pollutants.

Table 5 .
The NPS outflow load and reduction in different filter strip scenarios.