Quantifying the Contribution of Agricultural and Urban Non-Point Source Pollutant Loads in Watershed with Urban Agglomeration

: Urban agglomeration is a new characteristic of the Chinese urbanization process, and most of the urban agglomeration is located in the same watershed. Thus, urban non-point source (NPS) pollution, especially the characteristic pollutants in urban areas, aggravates NPS pollution at the watershed scale. Many agricultural studies have been performed at the watershed scale; however, few studies have provided a study framework for estimating the urban NPS pollution in an urban catchment. In this study, an integrated approach for estimating agricultural and urban NPS pollution in an urban agglomeration watershed was proposed by coupling the Soil and Water Assessment Tool (SWAT), the event mean concentration (EMC) method and the Storm Water Management Model (SWMM). The Hun-Taizi River watershed, which contains a typical urban agglomeration and is located in northeastern China, was chosen as the study case. The results indicated that the per unit areas of total nitrogen (TN) and total phosphorus (TP) in the built-up area simulated by the EMC method were 11.9% and 23 times higher than the values simulated by the SWAT. The SWAT greatly underestimated the nutrient yield in the built-up area. This integrated method could provide guidance for water environment management plans considering agricultural and urban NPS pollution in an urban catchment.


Introduction
Point source pollution and non-point source (NPS) pollution have been identified as key triggers of deteriorating water quality. In recent decades, point source pollution has been fairly controlled [1][2][3], while NPS pollution has been recognized as a critical cause of water quality degradation [4]. Agricultural and urban NPS pollution are the two main sources of NPS pollution [5]. Agricultural NPS pollution is a major reason for eutrophication of streams and lakes, contributing up to 50% of the total nitrogen (TN) and total phosphorus (TP) in China [6,7], likely due to excessive pesticide and fertilizer application [8,9]. In 2018, a total of 59% of the Chinese mainland population lived in urban areas, and the proportion is expected to increase to 80% by 2050 [10]. China has the highest population in the world and has experienced unprecedented urbanization over the last several decades; fourteen urban agglomerations have formed in recent years [11]. Rapid urbanization has led to extensive land use changes [12], converting natural land to diverse urban impervious surfaces such as roofs and roads [13,14]. The increasing proportion of impervious surfaces can dramatically alter the urban hydrology cycle by decreasing latitudes of 40 • 27 to 42 • 19 N (Figure 1a,c). The watershed is composed of the Hun River, Taizi River and Daliao River, which flow mainly from northeast to southwest and enter the Bohai Sea in Yingkou city. It has a temperate continental monsoon climate. The average annual temperature from 1989-2018, in the study area was 8.7 • C, and the average annual precipitation was 700 mm. Approximately 70% of the precipitation in the watershed mostly occurs from June to September during monsoons. Forestlands are located mainly in the northeastern and southeastern regions and cover approximately 45% of the study area. Vegetation is dominated by mixed broad-leaved and coniferous forests, while dry farmland and paddy fields in the southwest cover 35% of the region (Figure 1b). The northeastern part of the watershed consists of low hills, while the middle and western parts are mainly alluvial plains. There are four main reservoirs (Dahuofang, Guanyinge, Shenwo and Tang), which are the most important drinking and irrigation water sources for Liaoning Province, with a total storage capacity of 6.9 × 10 9 m 3 . Business zone (BUZ), high density residential zone (HRZ), low density residential zone (LRZ), educational zone (EDZ), industrial zone (INZ), agricultural zone (AGZ), green space zone (GSZ), public service zone (PSZ), preservation zone (PRZ), development zone (DEZ), and water (WTR); S1, S2, S3, S4 and S5 were the rainfall sampling sites in BUZ, EDZ, HRZ, INZ and PRZ, respectively.
The study area covers most of the Central Liaoning Urban Agglomeration, which is concentrated in heavy industrial areas including six cities (Fushun, Shenyang, Benxi, Liaoyang, Anshan, and Yingkou) (Figure 1d). The urbanization rate of the study area was 57.84%, with urban and rural populations of 1.33 × 10 7 and 0.82 × 10 7 , respectively, in 2019. The agricultural crop-sown area was 1.4 × 10 4 km 2 , the total grain output was 8.48 × 10 6 t, and a large gross domestic product (GDP) was produced, which accounted for 47.3% of Liaoning Province in 2019 [40]. In recent years, urban sprawl has accelerated in the Hun-Taizi River watershed due to the policy of "Revitalization of Old Industrial Bases in Northeast China" proposed in 2003. Owing to rapid urban development, impervious surfaces are increasing significantly. Water quality sampling in the Hun-Taizi River watershed A total of 25 water quality sections were selected in the outlet of the main stream of each sub-basin, which was divided by the SWAT (Figure 1b). A total of 225 effective water samples were collected separately at each monitoring section by three sampling events with two or three parallel samplings in July and October 2018 and in April 2019. To avoid the impact of rainfall on water quality, three sampling events were conducted on sunny days. The concentrations of TN and TP were tested within 24 h of collection using alkaline potassium persulfate digestion-UV spectrophotometric and ammonium molybdate spectrophotometric methods, respectively [8,40,41].

2.
Rainfall runoff samplings in different urban functional zones The six cities in the study area are all cities with heavy industry located in the same watershed and climate zone. Five urban functional zones of Shenyang city, including the business zone (BUZ), educational zone (EDZ), high density residential zone (HRZ), industrial zone (INZ) and preservation zone (PRZ), were selected to monitor rainfall runoff (Figure 1d), which is sufficient to be representative of the whole study area. Runoff samples from the underlying surfaces of roofs, roads and grasslands were collected. Roof runoff samples were collected from building drainpipe outlets. Road runoff samples were collected at rain grates. Grassland runoff samples were collected at the outlets of sloping lawns. All rainfall runoff samples were collected using 500 mL polyethylene bottles by researchers with extensive sampling experience. Once runoff was produced, samples were collected. Once rainfall runoff started, samples were collected every 10 min for the first 60 min, every 30 min during a period of 60 to 180 min, and then every 60 min beyond 180 min until the runoff disappeared [19]. A total of 218 samples were collected during three rainfall events (8 July 2018, 14 August 2018 and 2 September 2018) in the five different functional zones in 2018 (Table S1), and 157 samples were collected during five rainfall events in 2012 [19]. The rainfall depth throughout the events was recorded using an automatic gauge with 0.25 mm/tip (RainLog 2.0).
All rainfall runoff samples were tested in the laboratory for total suspended solids (TSS), COD, TN, TP and six typical heavy metals (Pb, Cd, Cr, Cu, Ni, and Zn) within 24 h. The procedure outlined in the Standard Methods of the Examination of Water and Wastewater was adopted to measure the TSS concentrations [38]. The potassium dichromate method was used to test COD [35], and inductively coupled plasma-atomic emission spectroscopy (PerkinElmer, Shelton, Connecticut, USA, ICP-MS NexION 350X) was used to test six heavy metals [42].

1.
Geographic dataset The description of the geographic dataset (topographic data, land use and soil map) used to establish the SWAT is shown in Table S2. The digital elevation model (DEM) with 30 m resolution and a soil map with a 1:50,000 scale were collected. The land use map in 2017 was interpreted with Landsat Operational Land Imager (OLI) images. Seven land use categories were classified with an accuracy of 89.3% compared to 270 field survey points in 2017, including built-up area, water area, grassland, forestland, dry farm land, wetland and paddy field (Figure 1c).
Landsat 8 OLI images in 2017 were used to manually interpret the built-up area of Shenyang, Fushun, Benxi, Liaoyang, Anshan and Yingkou (Figure 1d). The land use map of the built-up area was interpreted based on the Sentinel 2 images in 2017 (10 m × 10 m), including construction land, road, grassland, forestland, farmland, water area and bare land ( Figure 2). The accuracy of interpretation was 88.5% compared with 360 field survey points in 2017. The roof dataset of the six cities was extracted from the QuickBird images (0.61 × 0.61 m), which was collected in 2013 using the software package Brista and updated in 2017 [39]. The land use map of the urban area in 2017 was updated by the roof map. For further analysis, the land use types were reclassified as roof, road (including road and other construction land), pervious surface (grassland, forestland, farmland and bare land) and water area ( Figure 3).

2.
Weather and hydrological dataset A detailed description of the weather and hydrological dataset used to establish the SWAT is shown in Table S2. The daily climate data used in the SWAT from 2001 to 2017 were collected from the China Meteorological Administration. The monthly hydrological data from 2001 to 2014 and water quality data from 2001 to 2005 were collected from the local government. The average annual precipitation of different rain intensities, which was used to estimate the urban NPS pollution, was calculated by the daily rainfall of the six cities during 1989-2018. The rain intensity was classified into light rain (0.1-9.9 mm), moderate rain (10-24.9 mm) and heavy rain (≥25 mm) based on the Grade of Precipitation (GB/T 28592-2012) ( Figure 4).

Agricultural NPS Pollution
The SWAT model, developed by the United States Department of Agriculture (USDA), was chosen as it is widely used to assess agricultural and watershed NPS pollution loads [41,43,44]. Detailed information about the SWAT is available in previous literature [45][46][47][48].

Urban NPS Pollution
The average annual urban NPS pollution loads can be calculated as follows [49]: where L i,j (t) is the average annual NPS pollution loading from different underlying surfaces and rain intensity during 1989-2018; EMC i,j (mg/L) is the mean EMC values of different underlying surface and rain intensity; P j (mm) is the average annual precipitation of different rain intensity during 1989-2018; ϕ i is the runoff coefficient of roof (0.9), road (0.8) and pervious surfaces (0.15) [45]; A i (km 2 ) is the area of roof, road and pervious surfaces; i the is diverse underlying surfaces of roofs, roads and grassland; and j is the different rain intensity of light rain, moderate rain and heavy rain.
The EMC values were used to calculate the average concentration of specific pollutants throughout the 8 rainfall events, including 3 rainfall events collected from different functional zones in 2018 and 5 from the EDZ in 2012. The average EMC values were used to represent the pollutant concentration for different rain intensities and underlying surfaces of the 6 cities. The EMC was computed as follows [33]: where EMC (mg/L) is the event mean concentration; M (mg) is the total mass load of a specific pollutant from different underlying surfaces during a storm event; V (L) is the total runoff volume; C t (mg/L) is the concentration of a specific pollutant at time t; Q t (m 3 /s) is the flow rate of stormwater runoff at time t; T (min) is the total runoff time; and ∆t (min) is the time interval. The SWMM, which was developed by the Environmental Protection Agency (EPA was chosen to simulate the urban NPS pollution loads in the absence of rainfall runoff samples of different rain intensities or underlying surfaces [46]. The SWMM has been extensively used in simulating hydrology and water quality in small urban areas [50]. Detailed descriptions of the SWMM can be found in the related literature [24,47].

Schematic of the Integrated Method Framework
The integrated agricultural and urban NPS study framework was proposed based on the SWAT, the EMC method, and the SWMM ( Figure 5). The data collected in the study area were used to quantify the agricultural NPS pollution loads. Geographic data and monthly monitoring data were adopted to establish, calibrate and validate the SWAT. Rainfall runoff sampling data from diverse underlying surfaces of different functional zones and rain intensities were used to calculate the EMC values of TSS, COD, TN, TP and heavy metals. In the absence of rainfall runoff sampling data, the geographic data and the rainfall runoff sampling data of the built-up area were used to establish, calibrate and validate the SWMM, which was adopted to simulate the NPS pollution loads. EMC values were calculated as the total load divided by the runoff volume. The urban NPS pollutant load was calculated using average EMC values of different underlying surfaces and rain intensities; the area of roof, road and pervious surfaces; the average annual precipitation of different rain patterns; and the runoff coefficient.

Agricultural NPS Pollution in the Hun-Taizi River Watershed
It is essential to properly calibrate and validate the SWAT before implementation in the simulation of agricultural NPS pollution. The coefficient of determination (R 2 ) and the Nash-Sutcliffe efficiency coefficient (E NS ) were chosen as the criteria to evaluate the model's performance [48]. The ranges of R 2 and E NS for streamflow were 0.64-0.99 and 0.58-0.98, respectively, at the nine hydrological stations ( Figure S1 and Table S3). The ranges of R 2 and E NS for TN were 0.63-0.95 and 0.36-0.87, respectively, and for TP, these values were 0.66-0.96 and 0.55-0.9, respectively, at the four water quality stations ( Figure S2 and Table S3). The calibrated and validated SWAT was able to successfully simulate agricultural NPS pollution loads.
Water quality samples from the Hun-Taizi River watershed and the simulated concentrations of TN and TP were used to further verify the SWAT simulation results. The concentrations of TN and TP were calculated by the simulated TN and TP yield divided by the runoff volume at each monitoring section. The R 2 values for TN and TP were 0.46 and 0.66 (p < 0.01), respectively ( Figure 6). There was a significant correlation between the NPS pollutant concentration simulated by the SWAT and the observed data. The average annual nutrient yield of TN was 54,200 t/y and that of TP was 6380 t/y during 2003-2017 in the Hun-Taizi River watershed. The per unit areas of agricultural land for TN and TP were 4.43 t/km 2 and 0.99 t/km 2 , respectively. The spatial distributions of TN and TP were similar, while TN was more widely distributed in the upstream of the basin (Figure 7). Agricultural NPS pollution was higher in the midstream and downstream of the watershed. Compared with other studies, the average annual yields of TN and TP were 7315 and 867 t/y, respectively, in the Xin'anjiang catchment [51], which were much lower than those in this study. This may be because the agricultural land area in this study area was 9984.33 km 2 and that in the Xin'anjiang catchment was 922.8 km 2 . The contribution of agricultural land to TN and TP was 5.47 and 0.58 t/km 2 in the Xin'anjiang catchment [51] and 3.11 and 0.83 t/km 2 in the Three Gorges watershed, respectively [52]. The average annual nutrient yields varied from basin to basin, while the contribution of the agricultural land was similar.

Characteristics of the EMC in Urban NPS Monitoring
The SWMM model parameters were adopted as shown in Tables S4 and S5. The rainfall runoff samples collected in 2012 and 2018 were used to calibrate and validate the SWMM. The values of R 2 and E NS for the SWMM for different underlying surfaces in the 14 August 2018 rainfall event in the PRZ were 0.78-0.86. The ranges of R 2 and E NS for the SWMM to simulate roof and road pollutant loads under heavy rain in different functional zones were 0.67-0.95. The R 2 and E NS for the SWMM to simulate grassland pollutant loads under different rain intensities and different functional zones were 0.58-0.83. The calibrated and validated SWMM was able to successfully simulate urban NPS pollution loads.
Road runoff had higher EMC values of TP, TN, COD and TSS than roof runoff (Figure 8a-c). The concentrations of TP, TN and TSS in road runoff for light rain were 1.03, 3.50 and 373 mg/L, respectively, and the concentrations for heavy rain were 0.7, 2.25 and 422 mg/L, respectively, which exceeded class V of the Environmental Quality Standard for Surface Water (TP 0.4 mg/L and TN 2.0 mg/L) (GB3838-2002) and the second category of the Integrated Wastewater Discharge Standard (TSS 200 mg/L) (GB8978-1996). The EMC of TSS in road runoff for moderate rain was 1210 mg/L, which was approximately 5 times higher than the second category of the Integrated Wastewater Discharge Standard (GB8978-1996). The concentrations of TP in pervious surfaces for light rain, moderate rain and heavy rain were 0.54, 1.83 and 0.86 mg/L, respectively, which exceeded class V of the Environmental Quality Standard for Surface Water (GB3838-2002). This suggests that greenbelts, parks and residential green spaces in urban areas are other major sources of TP. The concentration of TN was also higher in roof runoff, which implies that nitrogen deposition is a source of TN. Analyzing the EMC values for heavy metals under different rain intensities and underlying surfaces showed that road runoff had the highest concentration (Figure 8d-f). The concentrations of Pb and Cr in road runoff for light rain were 168 and 155 µg/L, respectively, and for moderate rain, they were 120 and 108 µg/L, respectively, which exceeded class V of the Environmental Quality Standard for Surface Water (100 µg/L) (GB3838-2002). The concentration of heavy metals for different rain intensities in descending order was light rain, moderate rain and heavy rain, which indicated that with increased rainfall, the runoff volume increased and further decreased the concentration of pollutants. This may be due to the first flush effect.
The total urban NPS pollution loads in the six cities for TP, TN, COD and TSS were 266, 853, 2250 and 192,000 t/y, respectively (  The total average loads of Cd, Ni, Pb, Cr, Cu and Zn in the six cities were 0. 35, 9.93, 25.9, 28.4, 36.1 and 139 t/y, respectively ( Figure 10). Shenyang city had the highest heavy metal loads of Cd, Ni, Pb, Cr, Cu and Zn, which were 0.09, 2.41, 6.08, 7.07, 8.79 and 33 t/y, respectively (Figure 10a). The contributions of roads to Cd, Ni, Pb, Cr, Cu and Zn in the six cities were 0.62, 18.07, 47.91, 50.7, 65.46 and 255.48 kg/km 2 , respectively, and those of roofs to Cd, Ni, Pb, Cr, Cu and Zn were 0.09, 1.43, 1.25, 6.64, 5.63 and 11.62 kg/km 2 . Roads had the highest heavy metal pollution, which accounted for more than 95% of the total amount of pollutants (Figure 10b). The spatial distribution of the annual urban NPS pollutant loads showed that the roads located in the central part of Shenyang city had the highest NPS pollution (Figure 11). This is consistent with the statistical analysis results.

Contribution of Agricultural and Urban NPS Pollution Loads
The annual agricultural TN and TP pollution loads simulated by the SWAT in this watershed were 44,200 and 9910 t/y, respectively (Table 1). According to the EMC and SWMM, urban NPS pollution contributed 853 t/y TN and 266 t/y TP in this watershed. The proportion of agricultural and urban TN yield was 98.1% and 1.9%, respectively, and the proportion of TP yield was 97.4% and 2.6%, respectively. This indicates that agricultural land was the main source of NPS pollution in the Hun-Taizi River watershed. This may be because the area of agricultural land and the built-up area were 9984.33 km 2 and 1141.68 km 2 , respectively, accounting for 35.99% and 4.21% of the study area. Table 1. The average annual and per unit area yield of total nitrogen (TN) and total phosphorus (TP) from different land use types simulated by the Soil and Water Assessment Tool (SWAT) and the event mean concentration (EMC) method. The agricultural TN and TP yields simulated by the SWAT were 4.43 and 0.99 t/km 2 , respectively, and the urban TN and TP pollution calculated by the EMC method was 0.75 and 0.23 t/km 2 , respectively. The urban NPS pollution of TN simulated by the EMC method was 11.9% higher than the value simulated by the SWAT. The contribution of TP from urban areas simulated by the EMC method was 23 times higher than the values simulated by the SWAT. The SWAT greatly underestimated the TP load.

Estimating the Characteristic Pollutants of Urban NPS Pollution in an Urban Agglomeration
Most current NPS pollution research mainly studies nitrogen and phosphorus pollution at the watershed scale [46,53,54], while very few researchers focus on evaluating heavy metal pollution. The quantity of heavy metals is estimated by combining the SWAT with sediment geochemistry [55,56] or a heavy metal module incorporated into the SWAT at the watershed scale [57,58]. We found that the average heavy metal loads (Pb 227, Cu 316, Cr 248, Ni 87, Cd 3.05 and Zn 1220 g/ha) were much higher than previous study cases in the agricultural watershed (Table 2) [56][57][58]. This indicates that the characteristic pollutants in the urban NPS pollution were much higher than those in the agricultural watershed. However, heavy metals were not included in the routine water quality monitoring program of each hydrological section. The SWAT cannot simulate heavy metal pollutant loads at the watershed scale due to the lack of long time-series monitoring data. Therefore, in this study, the characteristic pollutants in urban NPS pollution were not compared with these pollutants at the watershed scale. Rapid urbanization and increased impervious surfaces have changed the hydrological cycle and increased the amount of dust deposited with heavy metals [13,19]. With rainfall, heavy metals accumulated in the urban underlying surfaces are washed into receiving waters, which can seriously affect public health and threaten environmental quality. Controlling urban heavy metal pollution is an important aspect of surface water quality improvement, especially in an urban agglomeration watershed. Above all, urban NPS pollution, especially the characteristic pollutants in urban areas, should be considered when calculating the NPS pollution in an urban agglomeration watershed as this can enrich our knowledge of the NPS pollutant constituents at the watershed scale. Liuyang River -----0-500 [58] Heavy metal pollutant loads in the Muskeg River watershed were the main values in the spatial distribution in the Liuyang River and were the Zn loads in the surface flow.

Integrated Approach for Estimating Agricultural and Urban NPS Pollutant Loads
In this study, an integrated approach was proposed for estimating the agricultural and urban NPS pollution in an urban agglomeration watershed by coupling the SWAT, EMC method and the SWMM. At the watershed scale, most studies have focused on agricultural NPS pollution, especially TN and TP, neglecting urban NPS pollution. With rapid urbanization, urban NPS pollution loads have increased in recent years [36]. The urban NPS pollutant loads were quantitatively estimated by sampling data from diverse underlying surfaces of different functional zones and rain intensities, which included more observed data and better accuracy than similar research that only sampled in the same underlying surface or functional zone [13,51,59]. A total of 225 water quality samples were used to verify the spatial distribution of agricultural NPS pollution of TN and TP. This provided a more efficient and reasonable method to verify the SWAT simulated results. It was more reliable than studies that calibrated and validated the SWAT only for short periods or did not validate the tool due to data limitations [53,54,60,61]. We found that the TN and TP yields in the built-up area simulated by the EMC method were 11.9% and 23 times higher than the values simulated by the SWAT. This suggests that the SWAT underestimated nutrient yields in the urban area. This may be because the SWAT has a limited capability to simulate the location problem of impervious areas within an urban catchment [62]. Our results indicate that the SWAT coupled with the EMC method and SWMM was an effective tool to estimate the agricultural and urban NPS pollutant loads in an urban agglomeration basin. The framework proposed in this study can help researchers in estimating agricultural and urban NPS pollutant loads in urban agglomeration basins and can also help watershed planners develop effective management strategies to mitigate NPS pollution.

Uncertainty Analysis and Implication
The simulation uncertainty of the integrated model was a major concern even though an overall good simulation of nutrient yields was obtained in this study [8]. Generally, integrated modeling involves two kinds of uncertainty: (1) uncertainty due to inaccurate input and calibration data; (2) systematic model uncertainty regardless of each correct model [20,63]. In this study, three years (2003)(2004)(2005) of observed water quality data were used to calibrate and validate the SWAT; water quality stations were not located in watershed outlets as they would have introduced uncertainty to the simulation results [8,57]. The SWMM used to calculate the urban NPS pollutant loads generated by the diverse underlying surfaces might have introduced prediction errors due to lack of drainage system data [50]. Thus, long-term continuous monitoring at the outlet of the watershed, the urban drainage system, and the simulation uncertainty of the integrated models should be considered in future studies.
The urban pollutant loads were estimated by the area of underlying surfaces, EMC values, mean annual rainfall and runoff coefficient. We found that the contribution of urban TN pollution simulated by the SWAT and EMC methods was 0.67 and 0.75 t/km 2 , respectively, while the total annual pollution loads differed greatly. This may be because the land use map used to estimate urban NPS pollution had a higher resolution than the land use map used to simulate agricultural NPS pollution. In this study, the runoff coefficient we used was obtained from the related literature rather than actual values monitored from the different underlying surfaces of the different functional zones [59]. Rainfall runoff samples were collected in Shenyang city, but they were not obtained in the other five cities.
The average values of the EMC and runoff coefficient were used. This neglected the spatial variation in the six cities, which may bring some uncertainties in calculating the urban NPS pollutant loads [29,33]. The NPS pollutant loads of built-up areas in the SWAT were only substituted with urban NPS pollution load in the EMC and SWMM. In the future, complete model coupling should be further developed.

Conclusions
In this study, an integrated approach including the SWAT, measured EMC values and the SWMM was used to quantitively estimate the agricultural and urban NPS pollution loads in the Hun-Taizi River watershed located in the Central Liaoning Urban Agglomeration. The proportions of agricultural and urban NPS pollution loads for TN were 98.1% and 1.9%, respectively, and for TP, they were 97.4% and 2.6%, respectively. Estimating urban heavy metal loads is an important aspect of surface water quality improvement and can enrich NPS pollutant constituents at the watershed scale. The TN and TP yields in the built-up area simulated by the EMC method were 11.9% and 23 times higher than the values simulated by the SWAT, respectively. This implies that the SWAT underestimated nutrient yields in the urban area. The SWAT and EMC method coupled with the SWMM was an effective tool for estimating agricultural and urban NPS pollutant loads that could be applied to other urban catchments for water quality management.