Regionalization of Agricultural Nonpoint Source Pollution over China with a Combination of Qualitative and Quantitative Method

: Agricultural nonpoint source pollution has been a serious problem in China; however, currently a lack of basic data and quantitative analysis hinders control and reduction of agricultural nonpoint source pollution. Therefore, it is necessary to explore a regionalization method in the study of nationwide agricultural nonpoint source pollution over China. This paper proposes a method of combining both quantitative calculation and qualitative analysis. Based on agricultural nonpoint source pollution mechanism, we ﬁrst proposed the natural environment index, which was calculated from relief degree of land surface, thermal humidity index, water resources quantity and precipitation index, and land cover index. Second, we proposed basic agricultural environment index, which was calculated based on the area of cultivated land use and the quality of integrated soil fertility. Third, we simpliﬁed the spatial distribution of natural environment and basic agricultural environment with the method of choropleth map classiﬁcation, thematic map series, and gravity centers curve. Fourth, we conducted a qualitative analysis for both the natural environment and basic agricultural environment by overlaying the classiﬁcation and existing regionalization maps to reveal the intra-region homogeneity and inter-region heterogeneity with a high reliability. The regionalization method used in this study resulted in a nationwide regional zoning of agricultural nonpoint source pollution over China, and China can be divided into 10 regions, which can be a trustworthy reference for agricultural nonpoint source pollution study and management.


Introduction
In recent years, with livestock pollution playing a more and more prominent role in agricultural nonpoint source (NPS) pollution [1], agricultural nonpoint source (NPS) pollution has become a particularly serious problem for water quality management in China [2]. Although many studies have focused on agricultural NPS pollution estimation [3], evaluation [4][5][6], and management [7][8][9][10][11], studies on the entirety of China have rarely been carried out. This is because when Chinese scholars pay attention to the problem of agricultural NPS pollution, at first, the research area is limited only to the region where agricultural NPS pollution is located [8,9,11]. Under this circumstance, it is necessary to develop an environmental management policy on agricultural NPS pollution on a large scale, including livestock pollution. However, due to the heterogeneous problems and challenges [12], it is not sufficient to apply an uniform environmental policy on a national scale [13]. China's management policy is data, are the calculation of the mean value of every eight days and its annual average value. Data from these two sources were used to calculate the land cover index (LCI).
Data of the paddy fields and dryland area were derived from the land use datasets. Cropping index, based on the Chinese Farming Systems Regionalization, was rasterized by vector data, which were obtained from the Chinese Academy of Agricultural Sciences, the Institute of Agricultural Resources and Regional Planning. Data from the three sources were used to calculate the cultivated land use (CLU) [25].
Soil total nitrogen, soil total phosphorus, soil total potassium, and soil organic matter data were obtained from the National Science and Technology Infrastructure Platform-Data Sharing Infrastructure of Earth System Science, the Chinese 1:4,000,000 soil fertility mass distribution atlas. The data are based on the second national soil survey results of the Chinese 1:4,000,000 and 1:1,000,000 soil maps. All above data were used to calculate the integrated soil fertility (ISF).
Many existing regionalization maps can provide abundant qualitative information. These regionalization maps were used in this research to explore partitioning scheme and identify strategic boundaries. Therefore, we collected data from the Physical Regionalization of China (PR) [26], Overall Regionalization of Agriculture in China (ORA) [27], the Comprehensive Livestock Regionalization of China (CLR) [28], and the Discharge Reduction Districting of Livestock of China (DRDL) [29]. The PR, ORA, and CLR were obtained from the National Sharing Platform of Agriculture Science Data. DRDL is the districting strategy used in the Ministry of Environmental Protection of People's Republic of China.

Agricultural NPS Pollution Mechanism
Agricultural NPS pollution, man-made impact on the natural environment (NE), is caused by the agricultural activities such as planting and breeding. NPS pollution is tightly linked to land use activities that determine the sources and magnitude of pollutant loadings to stream water [30]. Fertilizer application and farmland management behavior have an important influence on NPS pollution in rural areas [31]. Soil NPS pollutant contents under different vegetations had significant temporal dynamics [11]. The prioritizing management strategies based on the contribution of sources across climate and hydrology variability will be critical for controlling nonpoint source pollution [32]. For NPS pollution, the most relevant social indicators are those related to behaviors and factors that influence water quality improvement or protection [33]. Therefore, relevant factors of agricultural NPS include land use patterns [30], farming habits [31], topography [34], soil [11,34], vegetation [11], climate [32], hydrology [11,32], and socio-economic factors [33].
Such factors can be classified into natural environment, basic agricultural environment (BAE), socio-economic environment (SE), agricultural cultivation and breeding situation (ACB), agricultural NPS pollution generation (APG), emission reduction measures (ERM), and agricultural NPS pollution discharge (APD) [35]. The agricultural NPS mechanism is shown in Figure 1. NE and BAE represent regional natural background and its basic conditions of agricultural production, such as regional topography, geomorphology, hydrology, biology, soil, and agricultural habits. NE and BAE reflect the macro-scale differentiation. Thus, NE and BAE factors are selected in NE and BAE represent regional natural background and its basic conditions of agricultural production, such as regional topography, geomorphology, hydrology, biology, soil, and agricultural habits. NE and BAE reflect the macro-scale differentiation. Thus, NE and BAE factors are selected in Sustainability 2020, 12, 405 4 of 16 agricultural NPS pollution regionalization studies. SE and ACB represent human activity intensity in various areas, such as population, economy, energy consumption, crops cultivation, livestock, and poultry breeding. SE and ACB reflect the meso-scale differentiation under the existing macro ones. APG, ERM, and APD represent human disturbance to NE under existing human activity intensity. Livestock and poultry-breeding pollution is the principal part of agricultural NPS pollution. It is also consistent with the above-mentioned mechanism. Therefore, factors in agricultural NPS pollution mechanism should be considered in livestock and poultry breeding pollution regionalization studies.

NE Index System
As an important human activity, agricultural production is closely related to human survival and living environment. For the nationwide study, livestock distribution is closely linked to its product demand distribution, which is closely linked to population distribution. The viewpoint was also stated by Gerber, Chilonda, Franceschini, and Menzi [36]. Feng, Yang, Zhang, and Tang [37] analyzed the natural conditions of human habitat environment and selected topography, climate, hydrology, and vegetation as evaluation factors in the NE suitability study for human settlements in China. In order to simplify the calculation process, weight of variables in the habitat environment index (HEI) formula uses mean values from Feng's research, and the meaning of variables in the formula is the same as that in Feng's paper [37]. The NE index system in this study includes four single index models and one comprehensive model as follows: (1) Relief degree of land surface (RDLS) is a synthetic representation of the altitude and relief amplitude of a region. The RDLS model can be formulated as: where ALT is the regional average elevation calculated at grid size; Max(H) and Min(H) represent the highest and lowest altitude; P(A) and A are the flat area and total area of the region, A equals 25 km 2 in this paper.
(2) Thermal humidity index (THI) characterizes regional climatic condition. The THI model can be formulated as: where t is the average temperature of each month and f is the monthly average relative humidity.
(3) Water resources quantity and precipitation index (WRI) describes the regional water resource conditions. The WRI model can be formulated as: where P is normalized precipitation; W a is normalized water area.
(4) Land cover index (LCI) reflects regional land use status. The LCI model can be formulated as: where NDVI is the normalized difference vegetation index of the region; LT i is the weight of various land use types. (5) Habitat environment index (HEI) is a comprehensive model, characterizing the NE suitability for human settlements in each region of China. The HEI model can be formulated as:

BAE Index System
BAE is an indispensable basis for agricultural production, which reflects the regional utilization of arable land and the regional soil quality. There are two single index models and one comprehensive model in the BAE index system, as follows: (1) The area of cultivated land use (CLU) index is the area of arable land under the conditions of the region-specific farming systems (multiple cropping index), reflecting the regional utilization of arable land. The arable land is selected and also considered as the consumption area of livestock manure. The CLU model can be formulated as: where CLU is cultivated land use index; Sp is the surface area of paddy fields; Sd is the surface area of dry land; MCI is the multiple cropping index.
(2) The quality of integrated soil fertility (ISF) index characterizes the relative quality of soil fertility in various regions, including soil total nitrogen, total phosphorus, total potassium, and organic matter.
where ISF is the quality of integrated soil fertility index; SOM is single factor fertility quality grade of soil organic matter; STN is single factor fertility quality grade of soil total nitrogen; STP is single factor fertility quality grade of soil total phosphorus; STK is single factor fertility quality grade of soil total potassium. The ISF index model was similarly proposed by a previous study [38].
(3) Integrated basic agricultural environment (IBAE) index is a comprehensive model, characterizing the basic agricultural conditions in all regions of China.
where IBAE is the integrated BAE index; NCLU and NISF are normalized CLU and ISF.

Indicators Normalization
To facilitate the comparison among all indicators, each indicator was normalized. In general, NE indicators' normalization method is similar to Feng [37], while the BAE indicators were done with positive normalized methods.

Regionalization Method
The multi-level division of quantitative indicators can provide more information for making qualitative determination, but the number of division levels should be determined in a subjective way. In this study, choropleth map classification (CMC), thematic map series (TMS), and gravity centers curve (GCC) methods were used to determine the levels of NHEI and NIBAE, as similarly done by Fu [1] and Ge and Feng [39]. Both CMC and TMS methods can clearly show the spatial distribution of quantifiable indicators for identifying their spatial patterns.
It should be noted that the value of NHEI and NIBAE can be divided into 16 levels from low to high according to the 16 levels in the choropleth map classification method. For NHEI, level 1-16 can be interpreted as the natural environment ranging from bad to good. For NIBAE, level 1-16 can be interpreted as the basic agricultural conditions ranging from barren to rich. The methods of TMS and GCC can help us to find the spatial distribution classification of NHEI and NIBAE. Exploring spatial distribution and combination relationship between different levels of the two indicators can help us define priorities for agricultural land management practices to mitigate NPS pollution. More attention should be paid to NPS treatment in areas with good natural environment and rich BAE conditions. Figure 2a shows the spatial distribution of the normalized habitat environment index in China at 1 km grid size. The NHEI in East China is larger than that in West China, and the NHEI in South China is larger than that in West China. Figure 2b indicates the spatial patterns of CMC for NHEI, where the NHEI values were equally divided into 16 levels, providing abundant information on the spatial distribution of NHEI. Figure 2c shows the spatial patterns of NHEI at each level determined by TMS,  Figure 2a shows the spatial distribution of the normalized habitat environment index in China at 1 km grid size. The NHEI in East China is larger than that in West China, and the NHEI in South China is larger than that in West China. Figure 2b indicates the spatial patterns of CMC for NHEI, where the NHEI values were equally divided into 16 levels, providing abundant information on the spatial distribution of NHEI. Figure 2c shows the spatial patterns of NHEI at each level determined by TMS, depicting the gradual transition of NHEI from Northwest China (low value) to Southeast China (high value). Based on the 16 levels, the gravity center of each level was identified and connected ( Figure 2d).  Table 1 lists the Euclidean distances between adjacent gravity centers of NHEI. According to the spatial distribution and distances between adjacent gravity centers, demarcation can be selected from 3-4, 8-9, 9-10, 10-11, 11-12. Demarcations are the distance breakpoint between the centers of gravity of different levels of NHEI, which can help us reclassify NHEI. In order to simplify the classification, level 1-3, 4-8, 9-10, 11, 12-16 are the five new levels. The adjusted classification map is given in Figure 2e.   Table 1 lists the Euclidean distances between adjacent gravity centers of NHEI. According to the spatial distribution and distances between adjacent gravity centers, demarcation can be selected from 3-4, 8-9, 9-10, 10-11, 11-12. Demarcations are the distance breakpoint between the centers of gravity of different levels of NHEI, which can help us reclassify NHEI. In order to simplify the classification, level 1-3, 4-8, 9-10, 11, 12-16 are the five new levels. The adjusted classification map is given in Figure 2e.  Figure 3a shows the spatial distribution of the normalized integrated basic agricultural environment in China at 1 km grid size. The high value of NIBAE is mostly distributed in East China, Northeast China, and Central China. The low value of NIBAE is virtually distributed in West China. Figure 3b shows the spatial patterns of CMC for NIBAE. By dividing NIBAE values into 16 equal levels, the values are gradually increased from level 1 to level 16. The CMC map provides sufficient information on the spatial distribution of NIBAE. As shown in Figure 3c, TMS indicates the spatial patterns of NIBAE at each level. It is shown in Figure 3d that GCC of NIBAE, gravity centers of 16 levels, are calculated by mean center tools in ArcMap. GCC of NIBAE is formed by connecting every gravity center in order of the level (1-16). China, Northeast China, and Central China. The low value of NIBAE is virtually distributed in West China. Figure 3b shows the spatial patterns of CMC for NIBAE. By dividing NIBAE values into 16 equal levels, the values are gradually increased from level 1 to level 16. The CMC map provides sufficient information on the spatial distribution of NIBAE. As shown in Figure 3c, TMS indicates the spatial patterns of NIBAE at each level. It is shown in Figure 3d that GCC of NIBAE, gravity centers of 16 levels, are calculated by mean center tools in ArcMap. GCC of NIBAE is formed by connecting every gravity center in order of the level (1-16).  Table 2 lists the Euclidean distances between adjacent gravity centers of NIBAE. According to the spatial distribution and distance between adjacent gravity centers, demarcation can be selected from 1-2, 10-11, 15-16. Demarcations are the distance breakpoints between the centers of gravity of different levels of NIBAE, which can contribute to the reclassification of NIBAE. In order to simplify the classification, four new levels, level 1, 2-10, 11-15, 16 have been formed. The adjusted classification map is shown in Figure 3e.

Regionalization Map of Agricultural Nonpoint Source Pollution
According to the need of national management, Chinese scholars have done a lot of national level regionalization research by using a variety of spatial factors. These regionalizations show the macro differentiation at the national level and provide the national regionalization in the fields of  Table 2 lists the Euclidean distances between adjacent gravity centers of NIBAE. According to the spatial distribution and distance between adjacent gravity centers, demarcation can be selected from 1-2, 10-11, 15-16. Demarcations are the distance breakpoints between the centers of gravity of different levels of NIBAE, which can contribute to the reclassification of NIBAE. In order to simplify the classification, four new levels, level 1, 2-10, 11-15, 16 have been formed. The adjusted classification map is shown in Figure 3e.

Regionalization Map of Agricultural Nonpoint Source Pollution
According to the need of national management, Chinese scholars have done a lot of national level regionalization research by using a variety of spatial factors. These regionalizations show the macro differentiation at the national level and provide the national regionalization in the fields of natural, agricultural, and livestock in terms of various spatial factors. The present regionalizations, as a reference for qualitative analysis, are national management references commonly used by the government. We can see that the Overall Regionalization of Agriculture in China and the Comprehensive Livestock Regionalization of China better match the differentiation of the two indicators; therefore, we combined NHEI, NIBAE, ORA, and CLR for the final agricultural NPS regionalization in two steps as follows.  The first step was to select the Overall Regionalization of Agriculture in China districting scheme as the original reference and then to merge regions according to the variation of the normalized habitat environment index and the normalized integrated basic agricultural environment. In this step, relatively homogeneous areas were selected as references for merging zones. NIBAE matched the boundary of the Overall Regionalization of Agriculture in China, the Physical Regionalization of China and the Comprehensive Livestock Regionalization of China in West China, thus NIBAE was selected as the reference for West China. Correspondingly, NHEI was selected as the reference for South China along with both NHEI and NIBAE for Central China, East China, and North China. The second step was to adjust the regions depending on CLR in areas with larger NIBAE and NHEI variations, mostly in East China and North China.
The Agricultural NPS Pollution regional zoning map is shown in Figure 5.  The Agricultural NPS Pollution regional zoning map is shown in Figure 5.   Table 3 lists the grid number percentage of all levels in each region (row). Each agricultural NPS pollution region has its major types and is different from other regions, reflecting the inter-regional differences. Table 4 lists regional properties of the 10 regions. Figure 5. Map of the agricultural NPS pollution regional zoning in China. Table 3 lists the grid number percentage of all levels in each region (row). Each agricultural NPS pollution region has its major types and is different from other regions, reflecting the inter-regional differences. Table 4 lists regional properties of the 10 regions.  Underdeveloped, a crisscross area of agriculture and animal husbandry. NPS pollution regional zoning overlay maps. The first step was to make a qualitative comparison between regionalization maps. The agricultural NPS pollution regional zoning map, representing the normalized habitat environment index variation of level 2, level 3, and level 4, is better than the Physical Regionalization of China in eastern monsoon region (Figure 4b), especially in VI, VII regions. For the normalized integrated basic agricultural environment, the agricultural NPS pollution regional zoning map, representing the differentiation of each level, also overmatches PR, especially in the boundary at all levels.

Regional Properties of Each Agricultural NPS Pollution Region
For the Discharge Reduction Districting of Livestock, only South Region, Eastern Region, and Northeast Region are compliant with the normalized habitat environment index variations ( Figure  4d). All six regions are not clearly in accordance with the normalized integrated basic agricultural environment differentiation. Thus, the agricultural NPS pollution regional zoning map better represents the NHEI and NIBAE differentiation than DRDL. Figure 6. Agricultural NPS pollution regional zoning overlay maps. Table 5 lists the grid number percentage of all regions in each level (line). Several regions have similar grid number percentage in the same level of the normalized habitat environment index and the normalized integrated basic agricultural environment, reflecting the regional similarity.

Similarity within Regions
1. For NHEI level 1, region I, III, and region II, V, and region IV, VII-X are similar pairs; 2. For NHEI level 2, region II, VI, and region IV, V, and region VII-X are similar pairs; 3. For NHEI level 3, region II, VI, VIII, and region III, VII, and region IX, X are similar pairs; 4. For NHEI level 4, region I, VI, and region II, IV, and region III, V, IX are similar pairs; 5. For NHEI level 5, region I, II, IV, V, VI is a similar pair. 6. For NIBAE level 1, region I, III, and region VI, X, and region VIII, IX are similar pairs; 7. For NIBAE level 2, region II, III, and region IV-VI are similar pairs; 8. For NIBAE level 3, region I, VI, and region IV, IX, and region V, X are similar pairs; 9. For NIBAE level 4, region I, III, VI is a similar pair. Figure 6. Agricultural NPS pollution regional zoning overlay maps.
The first step was to make a qualitative comparison between regionalization maps. The agricultural NPS pollution regional zoning map, representing the normalized habitat environment index variation of level 2, level 3, and level 4, is better than the Physical Regionalization of China in eastern monsoon region (Figure 4b), especially in VI, VII regions. For the normalized integrated basic agricultural environment, the agricultural NPS pollution regional zoning map, representing the differentiation of each level, also overmatches PR, especially in the boundary at all levels.
For the Discharge Reduction Districting of Livestock, only South Region, Eastern Region, and Northeast Region are compliant with the normalized habitat environment index variations ( Figure 4d). All six regions are not clearly in accordance with the normalized integrated basic agricultural environment differentiation. Thus, the agricultural NPS pollution regional zoning map better represents the NHEI and NIBAE differentiation than DRDL. Table 5 lists the grid number percentage of all regions in each level (line). Several regions have similar grid number percentage in the same level of the normalized habitat environment index and the normalized integrated basic agricultural environment, reflecting the regional similarity.

1.
For NHEI level 1, region I, III, and region II, V, and region IV, VII-X are similar pairs; 2.
For NHEI level 2, region II, VI, and region IV, V, and region VII-X are similar pairs; 3.
For NHEI level 3, region II, VI, VIII, and region III, VII, and region IX, X are similar pairs; 4.
For NHEI level 4, region I, VI, and region II, IV, and region III, V, IX are similar pairs; 5.
For NHEI level 5, region I, II, IV, V, VI is a similar pair. 6.
For NIBAE level 1, region I, III, and region VI, X, and region VIII, IX are similar pairs; 7.
For NIBAE level 2, region II, III, and region IV-VI are similar pairs; 8.
For NIBAE level 3, region I, VI, and region IV, IX, and region V, X are similar pairs; 9.
For NIBAE level 4, region I, III, VI is a similar pair. Therefore, qualitative comparison and quantitative analysis show that regional zoning map reflects not only the similarities but also the differences between regions. The key agricultural NPS pollution risk regions are region IV (Northeast region), region VII (Huang-Huai-Hai region), region VIII (Southwest region), region IX (The Middle and Lower Reaches of the Yangtze River region), and region X (South China region).

Conclusions
NE and BAE index system is a good way to describe natural background and its basal condition of agrarian production, which represents one of the fundamental mechanisms in the agricultural NPS pollution. However, dividing the index into several levels needs a subjective methodology, which could be addressed with choropleth map classification, thematic map series, and gravity centers curve methods. By overlaying the existing regionalization boundaries, the spatial distribution of NE and BAE could be further improved for regionalization. By taking advantage of the series of above-mentioned methods and quantitative calculation and qualitative analysis, it is feasible to produce a nationwide agriculture NPS regionalization map in China.
The entirety of China can be divided into 10 regions with high reliability. Each region has its distinctive characteristics. Five of these regions are identified as agricultural NPS pollution risk areas, all of which have high NHEI and NIBAE level. That is to say, the NE index system and BAE index system can help us identify agricultural NPS pollution risk areas by having one of NHEI and NIBAE at a high level. By observing these five regions, it has been found that these regions are distributed in the east coast of China and some areas extending inland. This area is the southeast part of China, the main agricultural area and the main population gathering area. Therefore, the result of NE and BAE index system revealed the common factors of agricultural NPS pollution:

1.
There exists an environment (with higher HEI) that is suitable for human aggregation. That is to say, there are three conditions: Flat terrain (with smaller RDLS), humid climate (with higher THI and WRI), and luxuriant vegetation (with higher LCI).

2.
There exist a good agricultural environment (with higher IBAE) and developed agriculture. To be specific, there are two conditions: More cultivated land (with higher CLU) and more fertile soil (with higher ISF).
This scheme can function as a trustworthy reference for agricultural NPS studies and management. Therefore, how to apply the result to policy decision-making is of great significance. Firstly, managers can use the adjusted classification map for NHEI and NIBAE. The calculation of these indicators involves the multi-year average data, which are close to the mathematical expectation of the mean value of the elements and represent the equilibrium of the regional elements. The strictness of the policy declined with the decrease of NHEI or NIBAE. Secondly, managers can utilize the agricultural NPS pollution regional zoning map, which can warn managers to focus more on the agricultural development in the five risk areas, especially in regions with high NHEI level. These areas are the South China Region (X), the Middle and Lower Reaches of the Yangtze River Region (IX), the Southwest Region (VIII) and the Huang-Huai-Hai Region (VII). It should be noted that the most common feature of these areas is population aggregation.