Identifying critical source areas of nonpoint source pollution in a watershed with SWAT–ECM and AHP methods

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 NH₄-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.

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 km². 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 NH₄-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.

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 NH₄-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.

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 NH₄-N and TP loads from both rural sewage and livestock manure.

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., NH₄-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.

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 R² (>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, NH₄-N, and TP in the DRW.

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 NH₄-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.

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., NH₄-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.

In Figure 5 are the average loads of NH₄-N and TP in each subbasin, according to the SWAT–ECM modeling. The NH₄-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 NH₄-N and TP loads revealed distinct spatial patterns. In the case of NH₄-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 NH₄-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 NH₄-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.

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 NH₄-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 NH₄-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 NH₄-N and TP loads (Khanal et al. 2018).

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 NH₄-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, NH₄-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 NH₄-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 NH₄-N.

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 NH₄-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 NH₄-N. Seven subbasins (113, 118–120, 122–124) were found to overlap as classified CSAs, pointing to the strong correlation between NH₄-N and TP. These overlapping CSAs' subbasins occupied 9.39% of the DRW's area, accounting for 23% of its NH₄-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).

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).

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 NH₄-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 NH₄-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 NH₄-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.

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.

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.

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.

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