Abstract

Identification of critical source areas (CSAs) is pivotal for the management of nonpoint source (NPS) pollution of watersheds. Most studies focus on source (S) factors and ignore the driving (D) factors of such pollution. The Soil and Water Assessment Tool (SWAT) model and the export coefficient method (ECM) were incorporated to quantify the S factors of ammonia nitrogen (NH4–N) and total phosphorus (TP) as NPS pollution. Specifically, S factors coupled with D factors, including precipitation, slope, soil and land use, were regarded as multi–factors. Moreover, the analytical hierarchy process (AHP) method was adopted to determine the respective weights of multi-factors after overlaying the factor maps to identify the CSAs. These CSAs accounted for 23.86% of the total area, and generated 54.94% of NH4-N and 42.59% of the TP loads. In contrast with single and multi-factors, we found that using multi–factors having differing weights was more accurate for identifying CSAs. Our study results indicate this approach is reasonable for CSAs' identification in watersheds, and it can provide insights into different pollution sources and migration, thus providing a sounder basis for future decision–making.

HIGHLIGHTS

  • SWAT-ECM was incorporated to quantify nutrient loads in the watershed.

  • Both source and driving factors were used to identify CSAs.

  • AHP method can calculate the weight of different factors.

  • Considering multi-factors can improve the accuracy of CSAs identification.

Graphical Abstract

Graphical Abstract
Graphical Abstract

INTRODUCTION

Discussions regarding the deterioration of water quality, which is mainly caused by point source (PS) and NPS pollution, have dominated the aquatic environmental research agenda in recent years. NPS pollution is challenging to identify because of its spatial and temporal variability (Wang et al. 2018). In the past decades, much effort has been directed to the control of PS pollution and this has achieved promising results, whereas NPS pollution has gradually become a worldwide concern and daunting challenge to mitigate.

Previous studies have proven that some minority areas contribute disproportionately higher pollution loads in river basins (Giri et al. 2016). These particular areas, called CSAs (Winchell et al. 2015), produce more sediment and nutrient loss due to the combination of specific social activities and agricultural production activities, meteorological and hydrological conditions, as well as soil, land use, and topographic conditions (Lamba et al. 2016). Pollution originating from CSAs is closely related to the water quality of the whole basin. For cost-effective and efficient control of NPS pollution, it is extremely crucial to identify CSAs in watersheds, especially where manpower and financial resources are limited (Dong et al. 2018).

Accurately locating CSAs plays a significant role in the control of NPS pollution, for which sediments and nutrient loads have formed the focus of most previous studies (Khanal et al. 2018). Nitrogen and phosphorus are the decisive elemental nutrients that cause eutrophication and the deterioration of water quality (Jin et al. 2017). Although nutrient loads can greatly influence the determination of CSAs' location (Lou et al. 2016), relying on it alone may incorrectly identify CSAs, as this approach only emphasizes the outcomes and ignores the underlying processes involved. D factors, such as PCP, land use /land cover (LULC), soil, and slope, can further influence the identification of CSAs considerably (Hahn et al. 2014). Hence, explicitly considering both S and D factors is helpful for comprehensively evaluating CSAs in terms of their source and driving processes, an approach more conducive to watershed management by decision makers.

The determination of nutrient load is the critical step in identifying CSAs. The approaches to estimating nutrient load chiefly include applying the ECM (Johnes 1996) and physical–based watershed models. The ECM was derived from the unit load method, and its use prevails in areas where resources and available data are both limited (Yan et al. 2019). The accuracy of its calculation is acceptable and the method has been successfully applied in many studies.

Widely accepted, the physical-based watershed models are utilized to assess the nutrient load of NPS pollution (Li et al. 2018). However, implementing different models will obtain different results (Merwade et al. 2008). The SWAT model (Arnold et al. 1998) in particular is now one of the most commonly accepted models to predict the location of watershed CSAs (Giri et al. 2016). In China, the available monitoring data for water quality and quantity is both small and of a short duration; hence, the spatial distribution of target pollutants is difficult to calibrate and verify (Ongley et al. 2010). Therefore, combining a simple model with a complex model to simulate pollution load is likely more effective.

In this paper, the Dasha River Watershed (DRW) was selected as the research area. The objectives of this study are three-fold: (1) To use a SWAT–ECM incorporated method to estimate the nutrient loads of NH4-N and TP, based on the data availability and water quality monitoring indicators of the research area; (2) To apply the analytical hierarchy process (AHP) method to determine multi-factors' weights and build a CSAs discriminant model identifying CSAs spatially, by overlaying the maps of different factors; (3) To explore the influence of different factors on CSAs' identification and comparison of the different identification methods.

MATERIALS AND METHODS

Study watershed

The proposed work was conducted in the DRW, located in the northern part of Henan Province and the southern part of Shanxi Province (35°03′05″–35°34′45″N, 112°59′05″–113°44′35″E; Figure 1). The Dasha River's total length is 115.5 km and it has a drainage area encompassing 2,688 km2. The watershed has an elevation of 23–1,595 m (above sea level) and an annual PCP of 550–600 mm. The land use types in the DRW are mainly farmland, woodland, residential land, and grassland, representing 48.86%, 31.81%, 10.5% and 6.95% of the area, respectively. At the Xiuwu Hydrological Station, situated in the middle and lower reaches of the Dasha River, the 2010 through 2016 years of monitoring records demonstrates that this river's water quality is heavily polluted. In 2015, for instance, its monthly average concentration of NH4-N and TP reached 3.04 mg/l and 0.41 mg/l respectively. Therefore, the study of this watershed's pollution sources is indispensable to mitigate this problem before it worsens.

Figure 1

Location and hydrographic overview of the Dasha River Watershed (DRW) in China.

Figure 1

Location and hydrographic overview of the Dasha River Watershed (DRW) in China.

Computational framework

The computational framework for the identification of CSAs in this study, as based on S and D factors, is illustrated in Figure 2. The S factors were estimated using SWAT–ECM, and these sources included living sources, livestock breeding, fertilization, and cultivation. The D factors maps are presented in Figure 3, and their corresponding data details are given in Table 1. First, SWAT and ECM were incorporated to estimate the nutrient loads of NH4-N and TP. Second, thematic maps were established based on the selected factors, and the AHP method was adopted to calculate the different weight of each factor. Finally, the S and D factors' respective maps were spatially overlaid in ArcGIS v10.3 to determine the CSAs in the watershed.

Table 1

Basic data of the D factors

TypeSource
PCP Ten meteorological stations nearby and one within the DRW 
Soil Institute of Soil Science, Chinese Academy of Sciences; the Harmonized World Soil Database (HWSD) soil types and nutrient content were obtained from the China Soil Database (http://vdb3.soil.csdb.cn/
LULC FROM–GLC (Finer Resolution Observation and Monitoring of Global Land Cover) data, from Tsinghua University 
Slope Based on DEM data with a resolution of 30 m in ArcGIS; the data set was provided by the Geospatial Data Cloud site, Computer Network Information Center, Chinese Academy of Sciences (http://www.gscloud.cn
TypeSource
PCP Ten meteorological stations nearby and one within the DRW 
Soil Institute of Soil Science, Chinese Academy of Sciences; the Harmonized World Soil Database (HWSD) soil types and nutrient content were obtained from the China Soil Database (http://vdb3.soil.csdb.cn/
LULC FROM–GLC (Finer Resolution Observation and Monitoring of Global Land Cover) data, from Tsinghua University 
Slope Based on DEM data with a resolution of 30 m in ArcGIS; the data set was provided by the Geospatial Data Cloud site, Computer Network Information Center, Chinese Academy of Sciences (http://www.gscloud.cn
Figure 2

Study method framework for identifying CSAs.

Figure 2

Study method framework for identifying CSAs.

Figure 3

(a) Slope map, (b) LULC map, (c) annual PCP map, (d) soil type map. Notes: FLc = Calcaric Fluvisols; FLs = Salic Fluviosls; GLk = Calcic Gleysols; CMc = Calcaric Cambisols; CMe = Eutric Cambisols; LPk = Rendzic Leptosols; LVk = Calcic Luvisols; RGc = Calcaric Regosols; VRe = Eutric Vertisols.

Figure 3

(a) Slope map, (b) LULC map, (c) annual PCP map, (d) soil type map. Notes: FLc = Calcaric Fluvisols; FLs = Salic Fluviosls; GLk = Calcic Gleysols; CMc = Calcaric Cambisols; CMe = Eutric Cambisols; LPk = Rendzic Leptosols; LVk = Calcic Luvisols; RGc = Calcaric Regosols; VRe = Eutric Vertisols.

ECM and SWAT

Due to its high computational efficiency, ECM can be easily embedded in an optimization framework. In this paper, the ECM was used to estimate the pollutant emissions from different agricultural sources (e.g., livestock breeding, rural activities). The equation for it is as follows:
formula
(1)
where, L is the total nutrients load in the basin (kg·a–1); n refers to the type of land use or livestock and population; Ei is the output coefficient of nutrients in the i type land use (kg·hm–2·a–1) or livestock, population output coefficient (kg· ind–1·a–1); Ai is the area (hm2) or the number of livestock and population in the i type of land use; and Ii is the input quantity of nutrients in the i type.

The pollutant export coefficients (Table 2) were obtained from the first national pollution census (ONPC 2007; GNPC 2009). Guidelines for Chinese environmental and economic accounting (Yu et al. 2009) were used to calculate the NH4-N and TP loads from both rural sewage and livestock manure.

Table 2

Pollutant export coefficients for TP and NH4-N of livestock and wastewater that were used in this study

Livestock production source (g head–1/feather–1 d–1)
Rural activities source (g people–1 d–1)
TypeCowBeef cattleSwineSowGoatPoultryDomestic sewageBlack water
TP 62.46 10.52 5.99 11.18 1.23 0.23 0.08 2.06 
NH4–N 14.13 2.63 1.67 1.91 1.32 0.01 0.06 0.12 
Livestock production source (g head–1/feather–1 d–1)
Rural activities source (g people–1 d–1)
TypeCowBeef cattleSwineSowGoatPoultryDomestic sewageBlack water
TP 62.46 10.52 5.99 11.18 1.23 0.23 0.08 2.06 
NH4–N 14.13 2.63 1.67 1.91 1.32 0.01 0.06 0.12 

SWAT was implemented here for the DRW to simulate the migration and attenuation of pollutants. Details of the model setting can be found in literature Neitsch et al. (2002). The detailed river network vector relationships and subbasin divisions are shown in Supplementary Material, Figure S1 and Figure S2. The hydrological information of the Xiuwu Hydrological Station (e.g., streamflow) came from the Henan Provincial Department of Water Resources, while the water quality information (e.g., NH4-N and TP concentrations) was obtained from the Department of Ecology and Environment of Henan Province. In the field surveys, we determined the planting time, crop type, fertilizer application and so on. Agricultural management practices (Table 3) were derived from the results of field surveys and the first pollution source survey. The basic input data required for SWAT can be found in Supplementary Material, Table S1.

Table 3

Schedules of agricultural activities during the wheat–corn rotation in the DRW

CropBasal fertilizationTillagePlantingTopdressing 1Topdressing 2Topdressing 3Harvest
Corn June 8 June 8 June 8 June 15 July 5 – September 28 
750 kg ha–1 (18–12–15) – – 150 kg ha–1 (urea) 225 kg ha–1 (urea) – – 
Wheat October 2 October 2 October 2 October 15 December 5 April 10 June 1 
900 kg ha–1 (15–20–12) – – 60 kg ha–1 (urea) 75 kg ha–1 (urea) 185 kg ha–1 (urea) – 
CropBasal fertilizationTillagePlantingTopdressing 1Topdressing 2Topdressing 3Harvest
Corn June 8 June 8 June 8 June 15 July 5 – September 28 
750 kg ha–1 (18–12–15) – – 150 kg ha–1 (urea) 225 kg ha–1 (urea) – – 
Wheat October 2 October 2 October 2 October 15 December 5 April 10 June 1 
900 kg ha–1 (15–20–12) – – 60 kg ha–1 (urea) 75 kg ha–1 (urea) 185 kg ha–1 (urea) – 

Note: 18–12–15 and 15–20–12 refer to N-P2O5-K2O ratio in the applied fertilizer.

Model calibration and validation

In the present study, observed data of the streamflow and NH4-N and TP concentrations were used to calibrate and validate SWAT's operation on a monthly time step. The simulation period was 01/01/2016–31/12/2019, with the warming-up period from 01/01/2015 to 31/12/2015. Due to the limited data available, the streamflow and NH4-N, and TP were calibrated and validated in different periods. Streamflow calibration and validation was performed from 01/01/2017 to 31/12/2018 and 01/01/2019 to 31/12/2019. The calibration and validation of TP and NH4-N concentrations were conducted for 01/01/2016 to 31/12/2016 and 01/01/2017 to 31/12/2017, respectively. The Nash–Sutcliffe efficiency (NSE) and the coefficient of determination (R2) are commonly chosen to evaluate the performance of hydrological models (Equations (2) and (3)).
formula
(2)
formula
(3)
where, S is the predicted value, is the observed value, is the mean of predicted value, is the mean of observed value, and n is the total number of observations. NSE values can range from 1 to, where a value of 1 indicates a perfect model fit; R2 ranges from 0 to 1, with higher values indicating less error variance.

According to several literature reviews (Santhi et al. 2001; Moriasi et al. 2007; Van Liew et al. 2007), our study's calibration and validation values for the NSE (>0.6) and R2 (>0.5) were within the acceptable range (Supplementary Material, Figure S3). These two statistics suggested that our SWAT model has favorable predictive accuracy for the flow, NH4-N, and TP in the DRW.

Identification of CSAs

The locations of CSAs are determined by different identification approaches; hence, we should apply scientific and reasonable methods according to local site conditions and data availability. In China, most water management practices are based on a single pollution factor. However, the Ecological Protection and Restoration Project of Mountains, Rivers, Forests, Fields, Lakes and Grasses, recently implemented by the Chinese government, emphasizes systematic and comprehensive management practices. In this study, the nutrient yield per unit area of each subbasin was analyzed in descending order for the period of 01/01/2016–31/12/2019; then, results from the literature were adopted as reference (Ghebremichael et al. 2010; Liu et al. 2016) and combined with the management of this area, to define as CSAs those regions whose NH4-N or TP load cumulative contribution was at least 30%. In this process of multi-factors identification of CSAs at the subbasin scale, the interaction among these multi-factors are emphasized. By following the natural interruption method and the Jenks optimization method (Jenks 1967), the CSAs' index, which was calculated by following steps, was divided into three levels. Levels I and II were defined as the Non-CSAs, and level III defined as the CSAs.

Step 1. Data normalization

Using different factors and data dimensions will affect the evaluation results when trying to get the approximate correct solution of a Gaussian process (Dimitriadis et al. 2021). Therefore, to enable the comparison and spatial analysis of different variables, it is necessary to normalize the data set first.
formula
(4)
where, f is the normalized data of a given factor (S or D); a and b are the lower and upper limits of its normalization range, respectively, namely 0 and 1; the min(xi) and max(xi) are the major controllers. These factors are evaluated in an image format, and then the same spatial resolution and geographic coordinate system are assigned.

Step 2. Determine the weights of the multi-factors

The respective weights of differing factors, accordingly, will determine the appraisal and evaluation results. The AHP method is implemented to determine such weights (Wu et al. 2013); it was introduced by Saaty (1980) to solve complex decision–making problems.

The principal steps to determine and assign a weight by the AHP method are as follows: (i) Select the influencing factors, and build the hierarchical structure (Figure 4). By analyzing the primary controlling factors that affect the CSAs, the research object is divided into three levels: a target layer A (i.e., CSAs identification), a criterion layer B (i.e., the S factor layers, D factor layers), and a decision layer C (i.e., NH4-N, TP, PCP, soil, LULC, and slope). (ii) According to the 1–9 scale method (Saaty 1980) (Supplementary Material, Table S2), the two major factors at the criterion level and the six major factors at the decision–making level are scored by experts, and the relative importance of each factor is then evaluated according to the experience of farmers, protection experts, and extension personnel, among others. (iii) Establish a judgment matrix by quantitatively evaluating each factor and calculating the respective weights of primary control factors.

Figure 4

The AHP hierarchical structure model.

Figure 4

The AHP hierarchical structure model.

Step 3. Identify CSAs with multi–factors

After the weight of each factor has been determined, the normalized thematic map is registered into a new figure, by using the spatial image overlay processing function in ArcGIS. Next, the topology is reconstructed to form a new relation–attributes table, and the distribution map of CSAs is then established accordingly. To construct the CSAs' index model by overlaying thematic map, this formula was used:
formula
(5)
where, C is CSAs' index; W is the weight of a major controlling factor; f is the normalized data calculated in Equation (4); j is the given factors; i is the subbasin. When applied, the results demonstrate that C varies between 0 and 1, for which a larger value indicates a stronger critical contribution.

RESULTS AND DISCUSSION

Spatial analysis of nutrient loads

In Figure 5 are the average loads of NH4-N and TP in each subbasin, according to the SWAT–ECM modeling. The NH4-N and TP loads respectively varied from 0.07 to 19.41 kg/ha/yr (Figure 5(a)) and 3.91 to 37.63 kg/ha/yr (Figure 5(b)). The yield maps of NH4-N and TP loads revealed distinct spatial patterns. In the case of NH4-N, it was characterized by great spatial variation and the area with its largest loss was in south of the DWR, namely in subbasins 97, 118–122, 124. The human population there is primarily agricultural, whose living sources (sewage and livestock), along with the pollution load generated by fertilization and cultivation, are the prime NH4-N load sources at over 60% (Babaei et al. 2019). The land use is chiefly agricultural land, with the lowest (50%) in subbasin 97 and the highest (95%) in subbasin118. The areas noted as having less severe NH4-N losses are subbasins 5, 7, 8, 13 and 26, all located in the upper reaches of the Dasha River. The land use there is mostly forest land, consisting of sparsely populated areas with less human activities (Hua et al. 2019). The percentage of agricultural land across different subbasins was closely related to N load (Li et al. 2018), whereas forest area was negatively correlated with it.

Figure 5

Spatial distribution of nutrient loads in the DRW.

Figure 5

Spatial distribution of nutrient loads in the DRW.

The TP load was highest in the south of the basin, which included subbasins 117–120 and 122–124. The TP load contributed by living sources is relatively small, this differing starkly from that of NH4-N. The livestock production and agricultural activities were the principal components generating the TP load. The area where the least loss occurred is the middle of the basin, namely subbasins 36, 52, 86, 87, 99 and 107. This area is mainly urban, with less agricultural activities and livestock breeding and a higher human population density in the cities and towns.

Considering the analysis on the whole, the combination of different factors may lead to different losses of NH4-N and TP in the DWR. Land use types are closely related to the distribution of nutrient loads (Wu et al. 2015). In this respect, those areas featuring high agricultural activity and population density and livestock breeding areas tend to have higher NH4-N and TP loads (Khanal et al. 2018).

Identification of CSAs with the S factors

Mapping the spatial distribution of nutrient loads is beneficial for identifying CSAs in watershed planning. Figure 6 reveals the location and scope of CSAs as determined by the two kinds of nutrient loads investigated for the DRW. For NH4-N, a threshold of 30% was operated to distinguish its CSAs, which were found primarily concentrated in 12 subbasins (97, 108, 111–113, 118–124) in the south of the basin. There, NH4-N output ranged from 19.41 kg/ha/yr to 13.46 kg/ha/yr; mean output was 15.03 kg/ha/yr and covered 12.92% of the basin. In the case of TP, unlike for NH4-N, its CSAs are principally concentrated in the south-to-east area of the basin but also distributed in the upper and middle parts of the basin, with a total of 10 subbasins (113, 114, 116–124) accounting for 15.9% of the basin marked by a relatively continuous overall distribution. Seven subbasins were identified as overlapping CSAs for both TP and NH4-N.

Figure 6

CSAs' identification based on the S factors.

Figure 6

CSAs' identification based on the S factors.

Figure 7 depicts both the nutrient load rate and cumulative percentage of the subbasin and the corresponding cumulative watershed area percentages, in which the load rate is arranged from high to low (according to the unit subbasin nutrient load). The CSAs as determined by the single factor loads of NH4-N and TP are slightly different, however. Under the same contribution rate (i.e., 30%), the basin portion covered by TP was larger (Figure 7), at 2.98% more than the coverage of NH4-N. Seven subbasins (113, 118–120, 122–124) were found to overlap as classified CSAs, pointing to the strong correlation between NH4-N and TP. These overlapping CSAs' subbasins occupied 9.39% of the DRW's area, accounting for 23% of its NH4-N and 20.10% of its TP loads. However, the ranking of overlapping areas in each factor is not entirely the same, because their driving factors may differ (Niraula et al. 2013). The reason behind this disparity between them is the result of many complex interactions, including land use type, slope, management, and hydrology, among others (Pongpetch et al. 2015). These results highlight the necessity of conducting analyses of multi-factors (as opposed to single factors).

Figure 7

Cumulative percentages and rates of NH4-N (a) and TP (b) loads compared for their coverage in the watershed area.

Figure 7

Cumulative percentages and rates of NH4-N (a) and TP (b) loads compared for their coverage in the watershed area.

Identification of CSAs with multi–factors

The weights of each factor calculated by AHP method are presented in Table 4. The random consistency ratio (CR) was <0.10, and the judgment matrix was deemed consistent, in having passed the consistency test, whose calculation result was acceptable (Saaty 1980). Using these weights, the multi-factors CSAs could then be determined by overlaying thematic maps in ArcGIS (Figure 8).

Table 4

Weights of different indexes to the general goal

Influential factorsNH4-N loadTP loadPCPLULCSoilSlope
Weights 0.375 0.375 0.1139 0.0657 0.0352 0.0352 
Influential factorsNH4-N loadTP loadPCPLULCSoilSlope
Weights 0.375 0.375 0.1139 0.0657 0.0352 0.0352 

Note: CR = 0.00384; randomness index = 0.23.

Figure 8

Identification of CSAs based on multi-factors by the AHP method.

Figure 8

Identification of CSAs based on multi-factors by the AHP method.

In all, 24 subbasins were identified as CSAs: 76, 89, 92, 95, 97, 101, 104, 106, 108, 109, 111–124. The calculated CSAs' index in Equation (5) ranged from 0.52–0.92, these found predominantly distributed in the south of the watershed. In terms of nutrient loads, the area of CSAs accounted for 23.86% of the total area of DRW, being responsible for 54.94% of its NH4-N and 43.38% of its TP loads, respectively.

The comparison demonstrated that the number of CSAs subbasins determined by multi-factor method was 11 more than that determined by NH4-N, and 14 more than that determined by TP, which was located in the area of relatively high rainfall. Among the land use types in subbasins of 92, 95, 104, 106 and115, the proportion of land under cultivation was more than 80%, but they were not identified as CSAs according to their nutrient loads. The soil type was Calcaric Cambisols in subbasins of 76, 89 and 109, and their hydrological group was B. Subbasins of 76, 89 and 101 had a much larger slope than that of CSAs identified by the nutrient loads. Obviously, it is not possible to precisely identify CSAs via the S factor alone. The influence of D factor can be analyzed as follows.

The LULC is directly related to NPS pollution and affects the composition of CSAs (Bello et al. 2018). Here, the CSAs were generally distributed in cultivated land (73.22%) and where the largest populations among rural settlements (289 villages) occurred in the study area. The principal nutrient sources are excessive applications of fertilizers, domestic sewage discharge, and agricultural intensification, especially the raising of livestock. Furthermore, the way crops are cultivated greatly reduces the soil and water conservation of sloping farmland, and the poor farming practices result in more land under cultivation than is necessary. Similar to other research findings (Shen et al. 2020), the application of chemical fertilizer on such cultivated land coupled to the improper treatment of animal husbandry manure has significantly increased the loads of NH4-N and TP in the watershed.

The field survey results for the DRW showed that the actual amount of fertilizer applied by local farmers is much higher than that recommended by the government. For the sake of maintaining a high grain yield only, the utilization rate of resources is much reduced, which greatly undermines environmental security (Zhang et al. 2016). According to the Henan government gazette, the ‘Toilet Revolution’ is underway, with the proportion of septic tanks at rural resident discharge points now reaching 76%, leaving the domestic sewage produced by 24% of farmers still directly discharged into nearby waters without any treatment. This has become the predominant contribution to the nutrition load of living sources in CSAs of the DRW.

According to the HWSD soil classification, Calcaric Cambisols and Calcaric Fluvisols are the main soil types in CSAs, whose soil texture is largely loam. The premie factors relevant to soil in CSAs are land use changes and pollutants' composition (Bello et al. 2018). High nutrient application rates increase nutrient availability and fosters nutrients' accumulation in soil. Erosion of soil nutrients via decay and diagenesis is a crucial source of nutrient load (Markewitz et al. 2001). If cultivated too frequently, fertilization increases and bare soil exposure increases, such that the risk of nutrient losses from the soil becomes higher.

Slope determines the runoff volume and runoff time into the stream or river (Ghebremichael et al. 2010). The slope of CSAs is mostly 0–2° and distributed in the middle and south of the watershed, which is conducive to agricultural activities. The northern part of the watershed is characterized by mountainous landforms and nature reserves. The soil organic matter content in the steep slope areas is low and the crop productivity there is generally weak (Shi et al. 2004). Meanwhile, large tracts of forest and natural grassland can slow and limit the transportation of pollutants into rivers (Xin et al. 2017); these land uses can be identified as non–CSAs.

PCP is a major D factor due to its great spatiotemporal uncertainty, and changes to PCP regimes will lead to the distributional change of the CSAs (Tsuzuki 2015). The distribution of heavy-PCP areas is consistent with that of CSAs. Since the rainfall runoff is larger in the west than the east, high rainfall events promote the migration of nitrogen- and phosphorus-related forms of sediment to surface waters, especially when extreme rainfall events occur (Yang et al. 2016), because water largely flows from land to rivers along the surface and near surface horizontal flow paths (Johnes 1996). Our result is consistent with the findings reported by Huang et al. (2015) and Lou et al. (2016), in that PCP is evidently a key factor driving the emergence and spatial pattern of CSAs. The dynamic changes of runoff formation and nutrient loss are closely related to rainfall, soil properties, and land use type.

Recommendations

Identifying CSAs in river watersheds can provide support for farmers, researchers, and government departments to tackle and mitigate NPS pollution. For NPS pollution control in the DRW, the underlying principle is zoning management: the CSAs are the high–priority regions, and other regions are low–priority ones. Among the CSAs, those in the south of the watershed are the main ones, being home to a multitude of rural residential areas; hence, an efficient rural sewage collection and treatment facility maintenance system should be gradually established there. At the same time, the popularization rate of science- and technology-informed agriculture remains generally low. The survey results of some farmers in the watershed evinced that they still employ traditional methods to fertilize and cultivate the land. Instead, we recommend they choose the type, quantity, and frequency of fertilization application according to the test results of two primary soils in CSAs and the types of crops planted, so as to achieve a more efficient, targeted fertilization that improves the overall fertilizer utilization rate by crops. Livestock breeding in CSAs is pronounced and scattered, making it demanding to effectively deal with when trying to mitigate its pollution load. Therefore, the proportion of large-scale livestock breeding operations should be increased, and sound treatment rate of waste should be improved. In other regions, grass planting should be strengthened to improve vegetation coverage to lessen erosion risks and retain nutrients in soil.

CONCLUSIONS

To exactly identify CSAs at the subbasin scale, this study quantitatively calculated S factors and overlaid the D factors with different weights. From the results we drew the following conclusions:

  • 1.

    When identifying CSAs, researchers should take both S and D factors into account. The incorporated SWAT–ECM approach was used to calculate S factors in the studied watershed, whose results were proven to be reliable. Among D factors most relevant to this study area, the four factors PCP, soil, slope, and LULC were selected.

  • 2.

    According to the calculation results of multi-factors assigned differing weights, three regions having differentiated levels were divisible by utilizing the AHP model (constructed in ArcGIS), and those of the highest level were identified as the CSAs. This division provides a comparatively scientific basis for pollution source identification and watershed management.

  • 3.

    Differing weights of multi-factors can be quantified according to different factors in the study area. Importantly, CSAs can be identified more precisely for specific research areas, so that more context-dependent and pragmatic governance suggestions could be put forward.

For future research, we propose a focus upon more precise calculations of certain watershed factors, such as its nutrient load and high-resolution assessment of driving factor maps. The application and accuracy of this method now awaits evaluation in other river basins.

ACKNOWLEDGEMENTS

This work was supported by China National Natural Science Foundation (41430318, 41572222, 41602262, 41702261), National Key R&D Program of China (2016YFC0801800), Beijing Natural Science Foundation (8162036), Fundamental Research Funds for the Central Universities (2010YD02), Innovation Research Team Program of Ministry of Education (IRT1085) and State Key Laboratory of Coal Resources and Safe Mining. The authors would like to thank the editors and the anonymous reviewers for their valuable comments and suggestions.

DATA AVAILABILITY STATEMENT

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

REFERENCES

Arnold
J. G.
,
Srinivasan
R.
,
Muttiah
R. S.
&
Williams
J. R.
1998
Large area hydrologic modeling and assessment part 1: model development
.
J. Am. Water Resour. Assoc.
34
,
73
89
.
Ghebremichael
L. T.
,
Veith
T. L.
&
Watzin
M. C.
2010
Determination of critical source areas for phosphorus loss: Lake Champlain Basin, Vermont
.
Trans. ASABE
53
,
1595
1604
.
GNPC (The Group of Pollution Discharging Coefficient in the National Aquaculture Pollution Census) (In Chinese.)
2009
The Emission Coefficient Manual of Aquaculture Pollution in the First National Pollution Census
.
Hahn
C.
,
Prasuhn
V.
,
Stamm
C.
,
Milledge
D. G.
&
Schulin
R.
2014
A comparison of three simple approaches to identify critical areas for runoff and dissolved reactive phosphorus losses
.
Hydrol. Earth Syst. Sci.
18
,
2975
2991
.
Hua
L.
,
Li
W.
,
Zhai
L.
,
Yen
H.
,
Lei
Q.
,
Liu
H.
,
Ren
T.
,
Xia
Y.
,
Zhang
F.
&
Fan
X.
2019
An innovative approach to identifying agricultural pollution sources and loads by using nutrient export coefficients in watershed modeling
.
J. Hydrol.
571
,
322
331
.
Jenks
G. F.
1967
The data model concept in statistical mapping
.
Int. Yearb. Cartogr.
7
,
186
190
.
Lamba
J.
,
Thompson
A. M.
,
Karthikeyan
K. G.
,
Panuska
J. C.
&
Good
L. W.
2016
Effect of best management practice implementation on sediment and phosphorus load reductions at subwatershed and watershed scale using SWAT model
.
Int. J. Sediment Res.
31
,
386
394
.
Li
W.
,
Zhai
L.
,
Lei
Q.
,
Wollheim
W. M.
,
Liu
J.
,
Liu
H.
,
Hu
W.
,
Ren
T.
,
Wang
H.
&
Liu
S.
2018
Influences of agricultural land use composition and distribution on nitrogen export from a subtropical watershed in China
.
Sci. Total Environ.
642
,
21
32
.
Lou
H.
,
Yang
S.
,
Zhao
C.
,
Shi
L.
,
Wu
L.
,
Wang
Y.
&
Wang
Z.
2016
Detecting and analyzing soil phosphorus loss associated with critical source areas using a remote sensing approach
.
Sci. Total Environ.
573
,
397
408
.
Markewitz
D.
,
Davidson
E. A.
,
Figueiredo
R. D.
,
Victoria
R. L.
&
Krusche
A. V.
2001
Control of cation concentrations in stream waters by surface soil processes in an Amazonian watershed
.
Nature
410
,
802
805
.
Merwade
V.
,
Olivera
F.
,
Arabi
M.
&
Edleman
S.
2008
Uncertainty in flood inundation mapping: current issues and future directions
.
J. Hydrol. Eng.
13
(
7
),
608
620
.
Moriasi
D. N.
,
Arnold
J. G.
,
Van Liew
M. W.
,
Bingner
R. L.
,
Harmel
R. D.
&
Veith
T. L.
2007
Model evaluation guidelines for systematic quantification of accuracy in watershed simulations
.
Trans. ASABE
50
,
885
900
.
Neitsch
S. L.
,
Arnold
J. G.
,
Kiniry
J. R.
,
Srinivasan
R.
&
Williams
J. R.
2002
Soil and Water Assessment Tool Theoretical Documentation Version 2005. Blackland Research Center, Texas Agricultural Experiment Station, Temple, TX
.
Niraula
R.
,
Kalin
L.
,
Srivastava
P.
&
Anderson
C. J.
2013
Identifying critical source areas of nonpoint source pollution with SWAT and GWLF
.
Ecol. Model.
268
,
123
133
.
ONPC (Office of the Leading Group of the First National Pollution Source Survey of the State Council)
2007
Working Manual of the First National Pollution Source Survey
.
China Environmental Science Press
,
Beijing
.
Pongpetch
N.
,
Suwanwaree
P.
,
Yossapol
C.
,
Dasananda
S.
&
Kongjun
T.
2015
Using SWAT to assess the critical areas and nonpoint source pollution reduction best management practices in Lam Takong River Basin, Thailand
.
EnvironmentAsia
8
,
41
52
.
Saaty
T. L.
1980
The Analytic Hierarchy Process
.
McGraw–Hill
,
New York
.
Santhi
C.
,
Arnold
J. G.
,
Williams
J.
,
Hauck
L. M.
&
Dugas
W. A.
2001
Application of a watershed model to evaluate management effects on point and nonpoint source pollution
.
Trans. ASABE
44
,
1559
1570
.
Winchell
M. F.
,
Folle
S.
,
Meals
D.
,
Moore
J.
,
Srinivasan
R.
&
Howe
E. A.
2015
Using SWAT for sub–field identification of phosphorus critical source areas in a saturation excess runoff region
.
Hydrol. Sci. J.–J. Sci. Hydrol.
60
,
844
862
.
Xin
Z.
,
Jintian
C.
,
Yuqi
L.
&
Lei
W.
2017
Geo–cognitive computing method for identifying ‘source–sink’ landscape patterns of river basin non–point source pollution
.
Int. J. Agric. Biol. Eng.
10
(
5
),
55
68
.
Yu
F.
,
Wang
J.
,
Cao
D.
&
Jiang
H.
2009
Guideline for Chinese Environmental and Economic Accounting
.
China Environmental Science Press
,
Beijing
.
Zhang
W.
,
Cao
G.
,
Li
X.
,
Zhang
H.
,
Wang
C.
,
Liu
Q.
,
Chen
X.
,
Cui
Z.
,
Shen
J.
,
Jiang
R.
,
Mi
G.
,
Miao
Y.
,
Zhang
F.
&
Dou
Z.
2016
Closing yield gaps in China by empowering smallholder farmers
.
Nature
537
(
7622
),
671
674
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY-NC-ND 4.0), which permits copying and redistribution for non-commercial purposes with no derivatives, provided the original work is properly cited (http://creativecommons.org/licenses/by-nc-nd/4.0/).

Supplementary data