Identification of Critical Source Areas of Nitrogen Load in the Miyun Reservoir Watershed under Different Hydrological Conditions

The spatiotemporal distribution of critical source areas (CSAs) will change with hydrological conditions. In this study, the CSAs of nitrogen load under different hydrological conditions in the Chaohe River watershed were identified using the cumulative pollution load curve method determined from the nitrogen pollution simulated using the Soil and Water Assessment Tool (SWAT) model. The results showed that: (1) The order of factors impacting nitrogen load intensity is as follows: fertilization intensity, rainfall, runoff, land use type, slope type, and soil type. (2) The primary and secondary CSAs are concentrated in the upper and lower areas of the watershed, where cultivated land (8.36%) and grassland (52.55%) are more abundant. The potential pollution source areas are concentrated in the upper and middle areas of the watershed, where cultivated land (6.99%), grassland (42.37%), and forest land (48.18%) are evenly distributed. The low-risk source areas are concentrated in the middle and left bank of the watershed, where forest land (67.65%) is dominant and the vegetation coverage is highest. The research results have significance for improving the accuracy of the implementation of best management practices, and can provide a reference for the formulation of drinking water protection policies for Beijing.


Introduction
Non-point source (NPS) pollution refers to water pollution caused by dissolved or particulate pollutants from unspecified areas flowing into the receiving water body (such as rivers, lakes, reservoirs, and bays) through runoff processes under the influence of precipitation [1]. The main pollution sources include soil erosion, application of pesticides and fertilizers, rural livestock manure and garbage, irrigation using farmland sewage, urban surface runoff, forest surface runoff, and atmospheric wet and dry deposition [2,3]. A large number of studies have shown that NPS pollution is the main cause of deterioration and eutrophication of surface water bodies such as rivers, lakes, and reservoirs [4][5][6]. Especially in agricultural areas, water pollution caused by excessive nitrogen input has become the most serious water quality problem [7][8][9]. Research on the control and management of NPS pollution has important theoretical and practical significance for protection of the water environment in river basins and ensuring the safety of urban drinking water.
Compared with point source pollution, the temporal and spatial difference in NPS pollution is obvious because of the characteristics of the varying sources, intermittent discharge, and random occurrence. Therefore, it is more difficult, complicated and expensive to monitor, control, and manage information in the statistical yearbook from 2000 to 2010 and the field research by the authors, it was found that there are few industrial point sources in the basin, and the nutrient loads mainly come from NPS pollution such as livestock and poultry breeding, agricultural planting, and village residential life. Furthermore, the climatic type of the Chaohe River basin is a continental monsoon climate with four distinct seasons and obvious wetness, the hydrological conditions differ greatly in different periods. After 2007, due to changes in natural conditions and increased human activities, NPS pollution in the Chaohe River basin has shown an upward trend [52]. This basin is related to drinking water safety in Beijing, so this situation is becoming increasingly worrying [53,54]. There is an urgent need to do research related to protection of the water environment of the basin.

Data Sources
The spatial data includes Digital Elevation Model (DEM) data, land use/cover types, soil types, etc. Attribute data includes weather/climate data (precipitation, temperature, solar radiation, wind The geomorphic setting of the basin is mainly mountainous hills. The terrain is high in the northwest and low in the southeast, with an elevation range of 106-2209 m. In order to ensure the rationality and accuracy of the simulation, the land use and land cover map ( Figure 1) used is for the middle year (2000) of the simulation period. Land use types mainly include forest (54.79%, FRST), pasture (34.69%, PAST), agricultural land (7.67%, AGRL), residential land (1.26%, URBN), water (0.02%, WATR), and institutional land (1.54%, UINS). The classification is based on the land use classification system of the Chinese Academy of Sciences (CAS) Resource and Environment Data Cloud Platform, which divides land use types into the foregoing six first-level types and 25 more detailed second-level types. The year of the soil types map also is 2000, and its classification standard is the "Soil Map of the People's Republic of China at the Scale of 1:1,000,000" compiled and published by the Chinese Soil Census Office in 1995. The soil types in the study area mainly include 26.99% cinnamon soil (CIN), 18.67% fluvo-aquic soil (FLA), 18.14% aeolian sandy soil (AES), 13.17% brown soil (BRO), 9.42% coarse soil (COA), stony soil (STO), black soil (BLA), clay soil (CLA), meadow soil (MEA), chestnut soil (CHE), and paddy soil (PAD). The slope types are dominated by slopes greater than 15 • , accounting for 80.06% of the total area of the basin.
The basin is a typical agricultural area. The main crops are corn, soybeans, peanuts, wheat and chestnuts, and fertilizer application and farmland management are extensive. Combining information in the statistical yearbook from 2000 to 2010 and the field research by the authors, it was found that there are few industrial point sources in the basin, and the nutrient loads mainly come from NPS pollution such as livestock and poultry breeding, agricultural planting, and village residential life. Furthermore, the climatic type of the Chaohe River basin is a continental monsoon climate with four distinct seasons and obvious wetness, the hydrological conditions differ greatly in different periods. After 2007, due to changes in natural conditions and increased human activities, NPS pollution in the Chaohe River basin has shown an upward trend [52]. This basin is related to drinking water safety Sustainability 2020, 12, 964 5 of 22 in Beijing, so this situation is becoming increasingly worrying [53,54]. There is an urgent need to do research related to protection of the water environment of the basin.

Data Sources
The spatial data includes Digital Elevation Model (DEM) data, land use/cover types, soil types, etc. Attribute data includes weather/climate data (precipitation, temperature, solar radiation, wind speed, and relative humidity), hydrology data, water quality data, soil physicochemical properties, and socio-economic data. The data required for the SWAT model are listed in Table 1. The hydrological data and water quality data for this study are from the Xiahui Hydrological Monitoring Station. It is located in Sub-basin 39, which is the outlet of the Chaohe River basin. The frequency of collecting hydrological data is daily, and the frequency of collecting water quality data is monthly. Samples are collected according to surface water collection standards; when collecting water quality data, a mixed sampling method was used, each time about 200 mL of water sample was taken, sealed and taken to the laboratory for testing. There were a variety of water quality indicators tested, including total nitrogen, nitrate nitrogen, ammonia nitrogen, total phosphorus, and so on. But the main nutrients discussed in this study are total nitrogen and nitrate nitrogen. FOSS FIAstar 5000 was used to determine the concentration value, the detection method used was flow injection analysis. It is necessary to use the LOADEST model to establish a regression model of hydrological and water quality data, and calculate the monthly nitrogen load. This load is used with hydrological data as the monitored "data" during model calibration and verification.
The DEM can directly reflect the topographical features of the basin, and its accuracy determines the accuracy of the representation of the watershed in the SWAT modeling [55,56]. Through the hydrological analysis module, SWAT calculates the slope direction of the watershed, identifies the flow direction of the water, and generates the river network. On this basis, sub-basins are characterized and the hydrological parameters of each sub-basin are calculated. In this study, the Chaohe River basin was divided into 39 sub-basins. Sub-basin 8 has the maximum area of 282.46 km 2 , and Sub-basin 22 has the minimum area of 2.43 km 2 . The SWAT model was then built using basic data such as land use, soil type, and meteorology. The minimum thresholds for land use types, soil types, and slope types were set to 5%, 10%, and 10%, respectively, and the watershed was divided into 623 HRUs.

Data Processing
The data in the foregoing database was input into the SWAT model to simulate the Chaohe River basin. Based on parameter sensitivity analysis and the analysis of the spatial and temporal characteristics of nitrogen load, the model was calibrated and verified using the monitored hydrology data and water quality concentration data. Then multivariate correlation analysis was used to calculate the correlation coefficient between various factors and nitrogen load intensity, and the cumulative pollution load curve method [43] was used to determine the zonation of the nitrogen load in the basin under different hydrological conditions. The hydrological conditions here refer to the combination of three types of hydrological years-wet year, normal year, and dry year-and two types of hydrological periods-flood season and non-flood season-forming a total of six hydrological conditions. Finally, combined with the key factors impacting nitrogen load intensity, priority was given to the zoning of the wet year and flood season to finalize the CSAs for nitrogen load in the Chaohe River basin.

Model Parameter Analysis, Calibration, and Verification
The SWAT model [57] is a continuous simulation model developed to predict the effects of watershed management measures on water quality, sediment yield, and pollutants. The SWAT model has a large number of parameters, and the main parameters that have a decisive effect on the simulation results are screened by sensitivity analysis to calibrate the model to improve the efficiency and accuracy of the simulation. This study used the global sensitivity analysis method of Latin Hypercube Sampling (LHS) for parameter uncertainty analysis [58]. The specific process is: sampling the parameter values using the LHS method within the parameter value range, and running the model to obtain the sensitivity of the parameter, which is mainly calculated using the multiple linear regression system of the objective function and parameter values. The sensitivity of the parameter is determined by the absolute value of the statistical value t and the size of the p value. The larger the absolute value of t and the smaller the p value, the more sensitive the parameter is. The advantage of this method is that its analysis process is based on the global optimization over the full parameter space, which avoids the linear assumption of parameters in traditional analysis methods.
This paper selected the monthly average runoff and water quality data of the Xiahui station of the Chaohe River basin from 1985 to 2010 to calibrate and verify the model; 1985-1989 is the spin-up period for the model, 1990-2001 is the calibration period of the model, and 2002-2010 is the verification period of the model. The SUIF-2 optimization algorithm [59] in SWAT-CUP 2012 software was used to analyze and optimize the parameters related to runoff and nitrogen load. Based on the parameter sensitivity analysis and the analysis of the spatial and temporal characteristics of nitrogen load, the model was calibrated and verified using the monitored hydrology and water quality concentration data. After calibration, the simulated values of the runoff and nitrogen load at the outlet of the basin have the best consistency with the measured data.
The Nash-Sutcliffe efficiency (Ens) [60], coefficient of determination (R 2 ) [61] and relative error (RE) are used to evaluate the results of model calibration and verification. Ens is used to measure the degree of fit between the measured data and the simulated value. R 2 is used to measure the consistency of the trend between the measured data and the simulated value. RE is the ratio of absolute error (difference between the measured data and the simulated data) to the measured value. The calculation equations for the three indicators are: where P i is the predicted (simulated) value for the constituent being evaluated, O i is the observed (measured) value, P is the mean of predicted value, O is the mean of observed value, and n is the total number of observations. According to the model efficiency evaluation index of Moriasi et al. [62], when Ens ≥ 0.5, and at the same time, the RE of runoff simulation is within 25% for the flow, the model is considered satisfactory. If the RE of sediment and nutrient simulation is within 55% and 70%, respectively, the simulation results are considered satisfactory.

Statistical Analysis Method
Agricultural NPS pollution results from continuous and dynamic hydrological and water quality processes, including the rainfall runoff process, soil erosion process, surface pollutant leaching and precipitation process, and soil pollutant infiltration process. Therefore, regional NPS pollution is affected by the combined effects of rainfall conditions, surface hydrological conditions, and underlying surface characteristics. The process of pollutant loss is different for different underlying surfaces. The underlying surface conditions play an important role in rainfall infiltration and runoff and nutrient output. The key factors in the underlying surface conditions are land use types, soil types, and slope types. In recent years, with the development of society and the economy, the impact of human disturbance on NPS pollution in river basins has gradually increased. The representative behavior of human disturbance is excessive fertilization during agricultural planting, which has become an item that cannot be ignored in the identification of key factors impacting NPS pollution. Therefore, rainfall, runoff, fertilization amount, land use type, soil type, and slope type were selected as the representative factors influencing NPS pollution, and correlation analysis was done between these factors and the intensity of nitrogen load.
To calculate the relations between the impact factors and nitrogen load intensity for different hydrological years, the monthly mean value of rainfall and runoff in each sub-basin was taken. The fertilization intensity of each sub-watershed was obtained by dividing the fertilization amount by the cultivated land area. The proportion of the area of forest, grassland, and cultivated land was calculated based on the land use distribution map, and the proportion of different land use types in the total area of each sub-basin was obtained. The area ratio of different soil types and slope types was calculated in the same way. The results are shown in Figure 2. the impact of human disturbance on NPS pollution in river basins has gradually increased. The representative behavior of human disturbance is excessive fertilization during agricultural planting, which has become an item that cannot be ignored in the identification of key factors impacting NPS pollution. Therefore, rainfall, runoff, fertilization amount, land use type, soil type, and slope type were selected as the representative factors influencing NPS pollution, and correlation analysis was done between these factors and the intensity of nitrogen load.
To calculate the relations between the impact factors and nitrogen load intensity for different hydrological years, the monthly mean value of rainfall and runoff in each sub-basin was taken. The fertilization intensity of each sub-watershed was obtained by dividing the fertilization amount by the cultivated land area. The proportion of the area of forest, grassland, and cultivated land was calculated based on the land use distribution map, and the proportion of different land use types in the total area of each sub-basin was obtained. The area ratio of different soil types and slope types was calculated in the same way. The results are shown in Figure 2. Correlation analysis (CA) is applicable to data analysis involving a few variables. CA is the analysis of variables that have a causal relation in the population. CA can describe the closeness of the relation between objective things and express it with appropriate statistical indicators. Using multivariate correlation analysis to study various factors related to nitrogen load intensity, the Correlation analysis (CA) is applicable to data analysis involving a few variables. CA is the analysis of variables that have a causal relation in the population. CA can describe the closeness of the relation between objective things and express it with appropriate statistical indicators. Using multivariate correlation analysis to study various factors related to nitrogen load intensity, the degree of influence of various factors on nitrogen load intensity can be comprehensively analyzed, providing a basis for NPS pollution control in the Miyun Reservoir watershed.
Statistical Product and Service Solutions 25 (SPSS 25) data analysis software was used in the foregoing process, and the correlation coefficients between various factors and nitrogen load intensity under different hydrological conditions were calculated. The Pearson correlation coefficient was used to describe the linear relation between the selected variables. The Pearson correlation coefficient r has the following properties: (1) −1 ≤ r ≤ 1, the greater the absolute value of r, the stronger the correlation between the two variables. (2) 0 < r ≤ 1, indicates that there is a positive correlation between the two variables. If r = 1, it indicates that there is a complete positive correlation between the variables.
(3) −1 ≤ r < 0, indicates a negative correlation between the two variables. If r = −1, it indicates that there is a complete negative correlation between the variables. (4) r = 0, indicates no linear correlation between the two variables.

Zonation Method for Nitrogen Load
Accurate zoning of the watershed nitrogen load is an important basis for the configuration of management practices. The traditional zoning method is based on the load results calculated with the model. The multi-year average pollution load of each sub-basin is classified according to the natural breakpoint method, so that the sub-basin with the highest pollution load is identified and used as the target area for the application of management practices. However, this method will mask the spatial and temporal differences in nitrogen load in different hydrological years and hydrological periods, resulting in deviations in the identification results of CSAs from the actual watershed. Therefore, the implementation scheme of the management practices may be low in accuracy, and it is difficult to obtain the best pollution control effect.
Therefore, in order to make the division of risk source areas more accurate and reliable, the sub-basins with high nitrogen load intensity are promoted to CSAs, and the sub-basins with low nitrogen load intensity are the low-risk source areas. The cumulative pollution load curve method [43] was used to comprehensively zone the nitrogen load in the watershed under different hydrological conditions ( Figure 3). The abscissa indicates that the nitrogen load intensity of each sub-basin under each hydrological scheme is ranked from low to high, and the ordinate indicates the cumulative sum of the nitrogen load of each sub-basin as a percentage of the total nitrogen load [43]. This method is not simply dividing the nitrogen load of each sub-basin under different hydrological conditions; instead, the nitrogen load of each sub-basin is accumulated as a percentage of the total load until the accumulated percentage reaches the point to divide the risk source areas. Compared with the traditional natural breakpoint method for risk source areas division, this method uses the idea of sub-basin load accumulation, which can effectively prevent the risk source areas from being divided into several uniform parts, and, thus, relatively improves the scientific nature of the results. Finally, utilizing 20%, 40%, 60%, and 80% as the dividing points, the basin nitrogen load is divided into five risk level zones, which are the first-level critical source area, the second-level critical source area, the potential pollution source area, and the two low-level risk source areas.
conditions; instead, the nitrogen load of each sub-basin is accumulated as a percentage of the total load until the accumulated percentage reaches the point to divide the risk source areas. Compared with the traditional natural breakpoint method for risk source areas division, this method uses the idea of sub-basin load accumulation, which can effectively prevent the risk source areas from being divided into several uniform parts, and, thus, relatively improves the scientific nature of the results. Finally, utilizing 20%, 40%, 60%, and 80% as the dividing points, the basin nitrogen load is divided into five risk level zones, which are the first-level critical source area, the second-level critical source area, the potential pollution source area, and the two low-level risk source areas.

Model Parameter Analysis, Calibration and Verification
The SWAT-CUP software's parameter transformation mainly includes: replace (v), add (a), and multiply (r) [58]. In order to ensure the accuracy of the model simulation, it is necessary to run the SAWT-CUP program multiple times, and iterate the newly recommended parameter range every time the simulation is run. In this way, after several simulations, the final parameters will gradually converge to the optimal value. Parameter sensitivity, calibration and verification results are listed in Table 2. Note: To ensure the reliability of the model and the repeatability of the research, the parameters with the best value less than 1 are kept to 3 decimal places.
Sustainability 2020, 12, 964 10 of 22 The letter "t" is used here to represent the unit of measurement "ton". It can be seen from Figure 4 that SWAT simulation yields good results. The Ens of runoff is 0.85, the R 2 is 0.82, and the RE is 12.74%. The deviation of the total nitrogen simulation is higher than that of the runoff simulation, with Ens being 0.81, R 2 is 0.80 and RE is 14.56% in the complete simulation period. The deviation of the nitrate nitrogen simulation is higher than that of the total nitrogen simulation, with Ens being 0.79, R 2 is 0.77 and RE is 17.47% in the complete simulation period. Considering the guidelines of Moriasi et al. [62], the results of the model evaluation parameters show that the SWAT model achieved very good simulation results for the Chaohe River basin, and the accuracy meets the requirements for accurate NPS pollution simulation.

Temporal Characteristics of Nitrogen Load
In Figure 5, the abscissa represents the study period; the left y-axis represents the total nitrogen and nitrate nitrogen loads each year, the unit is tons/year (t/a); and the right y-axis represents the annual runoff, the unit is cubic meters/year (m 3 /a). It can be seen from Figure 5 that the annual runoff from 1990 to 2010 in the Chaohe River basin is between 5.67 × 10 7 −2.94 × 10 8 m 3 , and multi-year average annual runoff is 1.40 × 10 8 m 3 . During the 21 years of the study period, the annual load of total nitrogen is between 239.71-1033.58 t and the annual load of nitrate nitrogen is between 208.03-901.16 t, and the multi-year average annual loads of total nitrogen and nitrate nitrogen are 539.07 t and 471.90 t, respectively. The nitrate nitrogen load accounts for most of the total nitrogen

Temporal Characteristics of Nitrogen Load
In Figure 5, the abscissa represents the study period; the left y-axis represents the total nitrogen and nitrate nitrogen loads each year, the unit is tons/year (t/a); and the right y-axis represents the annual runoff, the unit is cubic meters/year (m 3 /a). It can be seen from Figure 5 that the annual runoff from 1990 to 2010 in the Chaohe River basin is between 5.67 × 10 7 −2.94 × 10 8 m 3 , and multi-year average annual runoff is 1.40 × 10 8 m 3 . During the 21 years of the study period, the annual load of total nitrogen is between 239.71-1033.58 t and the annual load of nitrate nitrogen is between 208.03-901.16 t, and the multi-year average annual loads of total nitrogen and nitrate nitrogen are 539.07 t and 471.90 t, respectively. The nitrate nitrogen load accounts for most of the total nitrogen load, with an average annual proportion of 87.54%. The three outputs-runoff, and total nitrogen and nitrate nitrogen loads-are significantly lower after 1999 than before 1999. The main reason is that the implementation of the policy of returning farmland to forest and grassland began in 1999 and the "Beijing-Tianjin Sandstorm Source Control Project Soil and Water Conservation Project" began in 2000 [63,64]. The forest coverage in the Chaohe River basin increased from 54.14% in 2000 to 76.96% in 2012. The rapid increase of forest land has a great impact on the water cycle process of the basin, which attenuates the direct response of surface runoff to precipitation while increasing evaporation [48,[65][66][67], thus, greatly reducing the runoff from the basin [68]. Different indicators of wet and dry periods correspond to different methods of wet and dry evaluation, and there are many different methods for selecting wet and dry indicators, such as frequency analysis, fuzzy clustering, Gamma distribution, and mean standard deviation. This paper uses the frequency analysis method to divide the different hydrological years. According to the multi-year monitoring data of various rainfall stations in Miyun Reservoir basin, the annual average rainfall of each station was arranged in descending order, and using Weibull's experience frequency equation, the mathematical expectation of rainfall in different years from each station can be obtained. Referring to Zhang [69] for the division of different hydrological years in the Miyun Reservoir basin, the criteria for the classification of wet and dry (wet year = 25%; normal year = 50%; dry year = 75%), it was determined that the years 1996, 1997, and 2002 were the wet year, normal year, and dry year, respectively. It can be seen from Figure 6 that in the different hydrological years, the total nitrogen and nitrate nitrogen loads in the basin are mainly concentrated in the flood season, that is, from June to October [49,54]. The total nitrogen load in the flood season accounts for 78.51%, 60.38%, and 55.57% of the full year load of the wet year, normal year, and dry year, respectively. The nitrate nitrogen load in the flood season accounts for 79.52%, 67.08% and 53.24% of the full year load of the wet year, normal year, and dry year, respectively. The annual total nitrogen load in the wet year is 1.53 and 3.70 times that of the normal year and dry year, respectively. The annual nitrate nitrogen load in the wet year is 1.46 and 3.75 times that of the normal year and dry year, respectively [70]. It can be seen that the total nitrogen and nitrate nitrogen load in the basin is very different in the different hydrological years and hydrological periods. It can be considered that the wet year is the key year for nitrogen load in the basin, and the nitrogen load is mainly concentrated in the flood season.  Different indicators of wet and dry periods correspond to different methods of wet and dry evaluation, and there are many different methods for selecting wet and dry indicators, such as frequency analysis, fuzzy clustering, Gamma distribution, and mean standard deviation. This paper uses the frequency analysis method to divide the different hydrological years. According to the multi-year monitoring data of various rainfall stations in Miyun Reservoir basin, the annual average rainfall of each station was arranged in descending order, and using Weibull's experience frequency equation, the mathematical expectation of rainfall in different years from each station can be obtained. Referring to Zhang [69] for the division of different hydrological years in the Miyun Reservoir basin, the criteria for the classification of wet and dry (wet year = 25%; normal year = 50%; dry year = 75%), it was determined that the years 1996, 1997, and 2002 were the wet year, normal year, and dry year, respectively. It can be seen from Figure 6 that in the different hydrological years, the total nitrogen and nitrate nitrogen loads in the basin are mainly concentrated in the flood season, that is, from June to October [49,54]. The total nitrogen load in the flood season accounts for 78.51%, 60.38%, and 55.57% of the full year load of the wet year, normal year, and dry year, respectively. The nitrate nitrogen load in the flood season accounts for 79.52%, 67.08% and 53.24% of the full year load of the wet year, normal year, and dry year, respectively. The annual total nitrogen load in the wet year is 1.53 and 3.70 times that of the normal year and dry year, respectively. The annual nitrate nitrogen load in the wet year is 1.46 and 3.75 times that of the normal year and dry year, respectively [70]. It can be seen that the total nitrogen and nitrate nitrogen load in the basin is very different in the different hydrological years and hydrological periods. It can be considered that the wet year is the key year for nitrogen load in the basin, and the nitrogen load is mainly concentrated in the flood season. nitrate nitrogen load in the flood season accounts for 79.52%, 67.08% and 53.24% of the full year load of the wet year, normal year, and dry year, respectively. The annual total nitrogen load in the wet year is 1.53 and 3.70 times that of the normal year and dry year, respectively. The annual nitrate nitrogen load in the wet year is 1.46 and 3.75 times that of the normal year and dry year, respectively [70]. It can be seen that the total nitrogen and nitrate nitrogen load in the basin is very different in the different hydrological years and hydrological periods. It can be considered that the wet year is the key year for nitrogen load in the basin, and the nitrogen load is mainly concentrated in the flood season.

Spatial Characteristics of Nitrogen Load
Based on the rainfall and nitrogen load data of the Chaohe River basin from 1990 to 2010, and using the cumulative pollution load curve method to partition them, the spatial distribution characteristics of nitrogen load in the Chaohe River basin were analyzed. It can be seen from Figure 7 that the average annual rainfall of each sub-basin is between 198.95-625.94 mm, with an overall average of 442.94 mm. The average annual runoff depth of each sub-basin is between 13.34-71.70 mm, with an overall average of 36.02 mm. The distribution of the rainfall and runoff depth is more consistent with the distribution of the total nitrogen load, and the consistency of the distribution of rainfall and runoff depth is relatively weak for the distribution of nitrate nitrogen load.

Preliminary Zonation of Nitrogen Load in Different Hydrological Conditions
The cumulative pollution load curve method [43] was used to classify the risk source areas of nitrogen load in the basin for the flood season of the wet year. Under the other five hydrological conditions, the classification of risk source areas of nitrogen load is based on the grade of the flood season of the wet year. It can be seen from Figure 8 that the maximum difference in nitrogen load intensity during the flood season occurs in the wet year, between 0-3.19 kg/ha, the normal year is between 0-1.87 kg/ha, and the dry year is between 0-0.81 kg/ha. In terms of spatial distribution, the first-level CSAs in the flood season of the wet year mainly are concentrated in the upstream and According to the multi-year average total nitrogen load intensity of the Chaohe River basin, the spatial difference of total nitrogen load is significant (p = 0.012), and the load intensity is between Sustainability 2020, 12, 964 13 of 22 0.42-2.31 kg/ha, with an average of 1.29 kg/ha. The Sub-basins 3, 8, 31, 33, and 35 have the highest nitrogen load intensity, with an average of 2.13 kg/ha, and mainly are concentrated in the upstream and downstream of the basin. That is, Xiaobazi and Kulongshan in the northwest, upper reaches of the basin, and the downstream towns of Fujiadian, Bakeshiying, and Anchungoumen. The reason for the high intensity is that these areas have high slopes and a high proportion of grassland area. As can be seen from Figure 2a,c, the grassland area in these sub-basins accounts for 58.83% of the total area, the slopes of these sub-basins are basically above 15 • , accounting for 90% of the total area, and Sub-basin 35 is even close to 100%. At the same time, the field investigations done by the authors found that there are several mining companies in the area, which result in greater damage to the surface vegetation, especially the vegetation coverage in mountainous areas. Therefore, soil erosion is very likely to occur due to surface runoff when rainfall occurs. The foregoing points undoubtedly increase the risk of nitrogen loss in the sub-basins. Sub-basins 6, 13, 21, 30, 34, and 38 are next in nitrogen load intensity which is between 1.59-1.79 kg/ha, with an average of 1.67 kg/ha. These sub-basins are mainly concentrated in the middle and lower reaches of the basin. The moderate nitrogen loads are mainly the result of the gentle terrain of these sub-basins, their proximity to the river, and their large cultivated area and frequent agricultural activities [71,72]. The total nitrogen load intensity in other areas is less than 1.59 kg/ha. The spatial differences in the nitrate nitrogen load in the Chaohe River basin also is significant (p = 0.023), and the load intensity is between 0.20-1.74 kg/ha, with an average of 1.07 kg/ha. The nitrate nitrogen load in Sub-basins 11, 13, 17, 21, and 22 is the highest, with an average of 1.68 kg/ha. It can be seen that this is different from the areas with the highest total nitrogen load intensity. The areas with the highest nitrate nitrogen load intensity mainly are concentrated in the middle reaches of the basin. The reason is that the nitrate nitrogen load mainly results from the use of fertilizers. Further, the livestock and poultry breeding industry in Sub-basins 11 and 13 is relatively high, and the livestock manure basically is applied to the local farmland. The cultivated land area of Sub-basins 17, 21, and 22 is relatively large and the application amount of chemical fertilizer is high. Sub-basins 16,24,25,30,34, and 37 are next in nitrate load intensity with the load of nitrate nitrogen between 1.46 and 1.52 kg/ha, with an average of 1.49 kg/ha. These sub-basins are mainly concentrated in the middle and lower reaches of the basin, and the nitrate loads mainly result from cultivated land, but the chemical fertilizer application volume in these sub-basins is moderate and there are few livestock and poultry breeding industries. The nitrate nitrogen load intensity in other areas is less than 1.46 kg/ha.

Preliminary Zonation of Nitrogen Load in Different Hydrological Conditions
The cumulative pollution load curve method [43] was used to classify the risk source areas of nitrogen load in the basin for the flood season of the wet year. Under the other five hydrological conditions, the classification of risk source areas of nitrogen load is based on the grade of the flood season of the wet year. It can be seen from Figure 8 that the maximum difference in nitrogen load intensity during the flood season occurs in the wet year, between 0-3.19 kg/ha, the normal year is between 0-1.87 kg/ha, and the dry year is between 0-0.81 kg/ha. In terms of spatial distribution, the first-level CSAs in the flood season of the wet year mainly are concentrated in the upstream and downstream of the basin, which are Sub-basins 3, 8, 31, 33, and 35, and the load intensity is between 2.35-3.19 kg/ha, and the average is 2.81 kg/ha. The secondary CSAs' load intensity is between 1.81-2.34 kg/ha; the potential pollution source areas load intensity is between 1.55-1.80 kg/ha; the sub-basins with a load intensity below 1.54 kg/ha are classified as low-risk source areas. The CSAs in the flood season of the normal year mainly are concentrated in the upper and middle reaches of the basin, which are Sub-basins 4, 6, 7, 10, 12, 13, 14, 16, and 24, with a load intensity ranging from 1.55 to 1.87 kg/ha, with an average of 1.66 kg/ha. The difference in the load intensity of the remaining sub-basins is small, ranging from 0 to 1.54 kg/ha; according to the division standard of the flood season of the wet year, all the sub-basins are classified as low-risk source areas. In the flood season of the dry year, the difference in nitrogen load intensity of each sub-basin is smaller, between 0-0.81 kg/ha. Combined with the division standard of the flood season of the wet year and the actual situation of the basin, since the load intensity is too small, it is not necessary to classify it, and all of the sub-basins can be classified as low-risk source areas. mainly concentrated in the flood season [69]. Therefore, there is little difference in the intensity of nitrogen load during the non-flood season in different hydrological years. In the non-flood season of the wet year, the nitrogen load intensity of Sub-basins 3 and 8 is slightly higher than that of other sub-basins, which are 1.19 and 1.61 kg/ha, respectively. Other sub-basins have low nitrogen load intensity, which are divided into low-risk source areas. In the non-flood season of the normal year and dry year, the maximum nitrogen load intensity is 0.69 and 0.32 kg/ha, respectively, which are already quite low. Similar to the case of flood season of the dry year, the nitrogen load intensity is too small, so there is no need for zoning, and all sub-basins are classified as low-risk source areas. The zonation results of nitrogen load intensity under different hydrological conditions are shown in Figure 8. The major land use and land cover data (Table 3) of each sub-basin in the Chaohe River basin were comprehensively used to calculate the proportion of various land use types in different risk source areas. As can be seen from Table 3, the proportion of grassland, cultivated land and forest land in the first and second CSAs accounts for 52.55%, 8.36%, and 36.27%, respectively, of the total area of the region. The proportions of these three types of land use in potential pollution source areas For many years, the nitrogen load intensity of the basin has been shown as wet year > normal year > dry year, and the nitrogen load intensity in the non-flood season is less than that in the flood season. As the Chaohe River basin belongs to the semi-humid and semi-arid climate, the rainfall is mainly concentrated in the flood season [69]. Therefore, there is little difference in the intensity of nitrogen load during the non-flood season in different hydrological years. In the non-flood season of the wet year, the nitrogen load intensity of Sub-basins 3 and 8 is slightly higher than that of other sub-basins, which are 1.19 and 1.61 kg/ha, respectively. Other sub-basins have low nitrogen load intensity, which are divided into low-risk source areas. In the non-flood season of the normal year and dry year, the maximum nitrogen load intensity is 0.69 and 0.32 kg/ha, respectively, which are already quite low. Similar to the case of flood season of the dry year, the nitrogen load intensity is too small, so there is no need for zoning, and all sub-basins are classified as low-risk source areas. The zonation results of nitrogen load intensity under different hydrological conditions are shown in Figure 8.
The major land use and land cover data (Table 3) of each sub-basin in the Chaohe River basin were comprehensively used to calculate the proportion of various land use types in different risk source areas. As can be seen from Table 3, the proportion of grassland, cultivated land and forest land in the first and second CSAs accounts for 52.55%, 8.36%, and 36.27%, respectively, of the total area of the region. The proportions of these three types of land use in potential pollution source areas are 42.37%, 6.99%, and 48.18%, respectively. Similarly, the proportions of these three types of land use in low-risk source areas are 22.81%, 6.56%, and 67.65%, respectively. It can be seen that the regional nitrogen load intensity has a significant positive correlation with the ratio of grassland and cultivated land area, and a significant negative correlation with the regional forest land area ratio. This finding is consistent with related research [48,67].

Identification of Key Impact Factors
It can be seen from Table 4 that under different hydrological conditions, the most significant factor affecting the nitrogen load intensity in the basin is the fertilization intensity. In addition, the correlation coefficient indicates that the nitrogen load intensity is positively correlated with the rainfall, runoff, percentage of cultivated land, and percentage of grassland. Due to the large amount of chemical fertilizer applied in the agricultural production process, the dissolved nitrate nitrogen is lost through leaching from the soil during a rainfall-runoff event. Additionally, the cultivated land needs long-term tillage, which leads to the migration of a large amount of organic nitrogen adsorbed to sediment reaching the surface water body during rainfall erosion. These results are consistent with other studies [73,74]. There is less dissolved nitrate nitrogen in grassland soil, so the nitrogen load from grassland is dominated by organic nitrogen in the adsorbed state. In the previous section, the correlation between the proportion of land use types in the sub-basin and the nitrogen load intensity was fully analyzed. The two analyses echo each other and fully prove the validity of this conclusion. Moreover, under different hydrological conditions, the nitrogen load intensity is positively correlated with the slope area ratio of 0 • -5 • , 5 • -10 • , and greater than 15 • , and the proportion of aeolian sandy soil area, because cultivated land and grassland are the main target areas for fertilization, and it is basically distributed in the area where the slope is relatively gentle (0 • -10 • ). Further, under the same rainfall intensity, the aeolian sandy soil is more likely to erode than other soil types, and for the higher the slope (>15 • ), it is easier to produce runoff and erosion, which provides an important carrier for nitrogen. Finally, there is a significant negative correlation between the nitrogen load intensity and the forest area ratio; this result is consistent with the results found by Bachman et al. [75]. As the forest land in the Chaohe River basin has been dominated by woody plants for many years, and many dwarf shrubs can cover the surface, these plants can retain more nitrogen and protect nutrients from runoff.
According to the foregoing analysis, under different hydrological conditions, the correlation coefficient of a single impact factor will slightly fluctuate. However, overall, the basic order of all the impact factors remains the same. For example, the correlation coefficients of fertilization intensity in wet year, normal year, and dry year are 0.878, 0.856, and 0.902, respectively. Although there are fluctuations, it is still the factor that has the greatest impact on the whole, followed by rainfall, runoff, land use type, slope type, and soil type. These results are consistent with related research [76,77]. Therefore, it is important to scientifically limit the type of fertilizer, fertilizer amount, fertilization period, and fertilization method, and to strengthen the implementation of water and soil conservation measures in river basins and related management practices for water quality protection for water sources.

Comparison of Zonation Methods of Risk Source Areas
There are various methods for zoning the nitrogen load intensity in the basin, and the natural breakpoint method and cumulative pollution load curve method are the most representative. This study mainly uses the cumulative pollution load curve method. Its four division points can effectively place the area with low load intensity in the low-risk area and the area with high load intensity to the high-risk area, so the scientific nitrogen load zoning results can be accurately obtained.
Liu et al. [42] used the natural breakpoint method (Jenks) for nitrogen load partitioning, and this method is generally suitable for classifying data with an uneven distribution but not biased at both ends. It is required that the data should be different, but the fluctuation range should not be too large, that is, the variance of the data sequence should not be too large or too small. In general, the area of the CSAs only accounts for a small part of the basin, but it contributes most of the load, thus, the difference in nitrogen load intensity between different regions should be obvious. For the Chaohe River basin, the variance of the data sequence is large and distributed to both ends. Thus, the natural breakpoint method is not applicable to the zonation of the nitrogen load intensity of the basin. Even if there are more special areas, that is, the physical and chemical properties of each area in the basin are relatively uniform, and the nitrogen load intensity is relatively close, at this time, the data sequence distribution is concentrated and the variance is small. In this situation, if the natural breakpoint method is used for the zonation of the nitrogen load intensity, the number of sub-basins in each risk source area will be evenly divided. At this time, the difference among the nitrogen load intensities of each sub-basin is not obvious, so the CSAs derived with the natural breakpoint method are not representative. In summary, the natural breakpoint method cannot obtain the comprehensive zoning result of nitrogen load intensity that best meets the actual situation of the Chaohe River basin. The cumulative pollution load method is the comprehensive zoning method of nitrogen load that is suitable for the actual situation of the Chaohe River basin.
There are some differences between the results of the risk source areas zonation of this study and those of Geng et al. [76]. The main reason is that the use of the hydrological conditions is different, and the method of dividing the risk source areas also is different. The wet year, normal year, and dry year selected for the two studies were 1996, 1997, and 2002, but Geng et al. [76] did not separate the flood season and non-flood season when the risk source areas were determined; rather the whole year's nitrogen load was used for comparison. Considering that the annual load does not fully demonstrate the spatial and temporal differences in nitrogen load during different hydrological periods in a year, and Geng et al.'s final nitrogen load zoning did not integrate the results of risk source areas zonation under different hydrological conditions, the results of their study cannot reflect the impact of changes in hydrological conditions on the zonation of CSAs. Furthermore, the method used to divide the risk source areas is the natural breakpoint method, and the shortcomings of this method have been previously discussed in detail. So, the CSAs obtained from Geng et al. [76] are different from the actual situation in the river basin.

Nitrogen Load Final Zonation Results
In the flood season and non-flood season, the total nitrogen load in the wet year is higher than that in the normal year and the dry year. Therefore, when the watershed nitrogen load is subdivided, the partitioning results of the annual nitrogen load of the wet year should be given priority. For example, in the flood season of the normal year, a sub-basin may be classified as a potential source area, but in the flood season of the wet year, it is classified as a primary source area, then the sub-basin ultimately belongs to the first-level CSAs. At the same time, the nitrogen load in the basin is mainly concentrated in the flood season. The nitrogen load in the flood season is several times that of the non-flood season. Therefore, in the final zonation of nitrogen load, the partitioning results of the flood season should also be observed. For example, in the non-flood season of the wet year, a sub-basin may be classified as a potential source area; however, in the flood season of the wet year, it is classified as a primary source area, therefore the sub-basin ultimately belongs to the first-level CSAs. Combined with the administrative division characteristics of the Chaohe River basin, the watershed is divided into three types of regions (four parts), which are the key governance area, mild control area, and ecological conservation area. The final zonation results of nitrogen load in the Chaohe River basin are shown in Figure 9.

Conclusions
(1) The total nitrogen load in the basin is dominated by the nitrate nitrogen load, and there are obvious differences in the spatial distribution of total nitrogen and nitrate nitrogen loads. The areas with high total nitrogen load intensity are mainly concentrated in the upstream and downstream of the basin, while the areas with high nitrate nitrogen load intensity are mainly concentrated in the middle reaches of the basin. Furthermore, the distribution of the rainfall and runoff depth is more consistent with the distribution of total nitrogen load, and the distribution consistency is relatively weak with the distribution of nitrate nitrogen load.
(2) The fertilization intensity is the most important factor affecting nitrogen load. Rainfall, runoff, land use type, slope type, and soil type are secondary factors that affect nitrogen load. In the cultivated land, grassland, and other fertilizer application areas, scientifically limiting fertilizer type, fertilizer amount, fertilization period, and fertilization method have great importance for the prevention and control of nitrogen load.
(3) The nitrogen load in the Chaohe River watershed can be divided into five risk-level zones, of which the primary CSAs (key governance area) and the secondary CSAs (key governance area) mainly are concentrated in the upstream and downstream of the watershed, where cultivated land (8.36%) and grassland (52.55%) account for a relatively high proportion of the total area. In this area, agricultural cultivation and farming activities are frequent, and soil erosion and nutrient loss are prone to occur. The potential pollution source areas (mild control area) mainly are concentrated in the upstream and midstream areas with relatively balanced proportions of cultivated land (6.99%), grassland (42.37%), and forest land (48.18%), as well as medium-intensity nitrogen loads occurring in this area. The low-risk source areas (ecological conservation area) mainly are concentrated in the middle reaches and left bank of the basin. This area is dominated by forest land (67.65%), and the vegetation coverage rate is the highest, so the soil erosion and the nitrogen load intensity are the lowest.
(4) Based on the division of the non-point source pollution risk source areas and the actual situation in the basin, the policy recommendations and mitigation measures to reduce the non-point source pollution load of the basin are given here. Specifically, cultivated land belonging to key governance areas with severe nutrient loss can conditionally be returned to forest land. If such a large-scale project cannot be achieved, then fertilization time, method, and amount can be scientifically arranged to minimize the NPS pollution degree. The mild control area is distributed along the river, and its main recommended control measure is a riparian vegetation filter belt. In areas with lighter NPS pollution, such as the ecological conservation area, closed management measures may be applied to avoid damage to forest resources and avoid the old path of pollution

Conclusions
(1) The total nitrogen load in the basin is dominated by the nitrate nitrogen load, and there are obvious differences in the spatial distribution of total nitrogen and nitrate nitrogen loads. The areas with high total nitrogen load intensity are mainly concentrated in the upstream and downstream of the basin, while the areas with high nitrate nitrogen load intensity are mainly concentrated in the middle reaches of the basin. Furthermore, the distribution of the rainfall and runoff depth is more consistent with the distribution of total nitrogen load, and the distribution consistency is relatively weak with the distribution of nitrate nitrogen load.
(2) The fertilization intensity is the most important factor affecting nitrogen load. Rainfall, runoff, land use type, slope type, and soil type are secondary factors that affect nitrogen load. In the cultivated land, grassland, and other fertilizer application areas, scientifically limiting fertilizer type, fertilizer amount, fertilization period, and fertilization method have great importance for the prevention and control of nitrogen load.
(3) The nitrogen load in the Chaohe River watershed can be divided into five risk-level zones, of which the primary CSAs (key governance area) and the secondary CSAs (key governance area) mainly are concentrated in the upstream and downstream of the watershed, where cultivated land (8.36%) and grassland (52.55%) account for a relatively high proportion of the total area. In this area, agricultural cultivation and farming activities are frequent, and soil erosion and nutrient loss are prone to occur. The potential pollution source areas (mild control area) mainly are concentrated in the upstream and midstream areas with relatively balanced proportions of cultivated land (6.99%), grassland (42.37%), and forest land (48.18%), as well as medium-intensity nitrogen loads occurring in this area. The low-risk source areas (ecological conservation area) mainly are concentrated in the middle reaches and left bank of the basin. This area is dominated by forest land (67.65%), and the vegetation coverage rate is the highest, so the soil erosion and the nitrogen load intensity are the lowest.
(4) Based on the division of the non-point source pollution risk source areas and the actual situation in the basin, the policy recommendations and mitigation measures to reduce the non-point source pollution load of the basin are given here. Specifically, cultivated land belonging to key governance areas with severe nutrient loss can conditionally be returned to forest land. If such a large-scale project cannot be achieved, then fertilization time, method, and amount can be scientifically arranged to minimize the NPS pollution degree. The mild control area is distributed along the river, and its main recommended control measure is a riparian vegetation filter belt. In areas with lighter NPS pollution, such as the ecological conservation area, closed management measures may be applied to avoid damage to forest resources and avoid the old path of pollution first and then treatment.