A Combination Model for Quantifying Non ‐ Point Source Pollution Based on Land Use Type in a Typical Urbanized Area

: Water pollution poses threats to urban environments and subsequently impacts the ecological health and sustainable development of urban areas. Identifying the spatiotemporal variation in non ‐ point sources (NPS) pollution is a prerequisite for improving water quality. This paper aimed to assess the NPS pollution load and then recognized the spatiotemporal characteristics of the pollution sources in a typical urbanized area. A combination model based on land use type was used to simulate the NPS pollution load. The results showed the following: (1) ponds and farmlands had higher pollution production intensities than other land use types, but the intensity and magnitude of pollution emissions were generally greater in urban areas; (2) monthly and annual total nitrogen (TN) and total phosphorus (TP) emissions had the same pattern as rainfall, and TN and TP emissions accounted for 56.2% and 58.0%, respectively, of the total in summer; (3) TN pollution was more serious than TP pollution in the study area, especially in farmlands; (4) urban runoff (UR) and livestock and poultry breeding (LPB) were the main sources of NPS, TN and TP emissions in the study area. If these NPS pollutants cannot be removed from this area, a large amount of freshwater is needed to dilute the current rivers to meet the requirement of the fourth category of China national environmental quality standards for surface water. This problem is serious in the control of polluted rivers in many cities throughout China.


Introduction
China has experienced rapid urbanization since the 1990s [1]. Unfortunately, urban population growth and land use changes have led to the degradation of flowing watercourses and the discharge of excessive wastewater [2]. Such water pollution, which threatens drinking water, agricultural irrigation and the ecological landscape, has become a key issue restricting urban development [3]. Consequently, the Chinese government has regarded urban water pollution with great importance and has invested substantial amounts of energy and money accordingly [4,5]. Nevertheless, some cities still face serious water pollution problems, especially rapidly developing cities where land use patterns are changing fast and sewage management and treatment are lacking [6][7][8]. To further control and improve the water quality situation in such cities, accurately recognizing of the sources of pollutants is of considerable significance.
In general, water pollution is caused mainly by point sources and non-point sources (NPS) pollution. The former, such as industrial and domestic pollution, is relatively easy to collect into the pipe networks and process in the sewage treatment facilities, whereas the latter, which is widely distributed and originates from multiple-source, is more difficult to monitor and control and has tremendous impact on the water quality of lakes and rivers [9][10][11]. Appropriately, many studies have focused on the assessment of NPS pollution loads. The most common and concerning sources of NPS pollution primarily include urban runoff (UR) [12], agricultural runoff (AR) [13], livestock and poultry breeding (LPB) [14] and aquaculture [15]. Ongley et al. [13] estimated the pollution statuses of agricultural and rural areas and their contribution to the total water pollution throughout China. Li et al. [16] simulated the NPS pollution load in the city of Baoding and analysed the effects of NPS pollution on Baiyangdian Lake. Due to the complexity of NPS components, empirical methods are often used to comprehensively evaluate all NPS pollutants, but these methods broadly lack detailed spatial and temporal descriptions [11]. Some studies have attempted to simulate the temporal variation in and the spatial distribution of pollutants, but most research has focused mainly on a certain type of NPS pollution, such as UR or AR [8,17,18]. In addition, because LPB and aquaculture do not correspond exactly to a specific type of land use, their spatial distributions are relatively difficult to describe [15,[19][20][21].
Accurate estimations of NPS pollution are important for the processing of NPS pollution loads [3,13,16]. The most common methods for evaluating NPS pollution loads include empirical evaluation methods and physics-based models [17,22]. The export coefficient model, which is based on a number of statistical materials to calculate pollution load, is the most common empirical method [23]. The methodology behind the export coefficient model is easy to understand, but this approach relies heavily on the accuracy of the empirical data and cannot explain the dynamic processes responsible for generating the pollution loads. Moreover, the use of empirical data rather than realtime data means that the results obtained using this method usually have a large time step (monthly, seasonal or annual). In contrast, physical-based models are employed primarily to establish the relationship of rainfall runoff pollution loads and simulate the NPS pollution loads during the rainfall runoff process [7,24]. These models generally contain good physical mechanisms, and some can provide information regarding the spatial and temporal distribution characteristics of NPS pollution loads [7,8,12]. Some common NPS models are shown in Table 1. Most of these models described above are expert in simulating specific pollution sources or pollutions of certain land use patterns. Urban rainstorm events are often the focus of attention in urban areas; for this purpose, models such as SWMM and STORM are commonly used [18,26,28]. Additionally, some models, such as SWAT and HSPF, are well suited to the evaluation of NPS pollution in a watershed or farmland [17,27,29], but are not applicable to urbanized areas. These NPS evaluation models generally encounter problems with insufficient water quantity or quality modules. Lee et al. [26] compared the water quality simulation effects of SWMM and HSPF in urban areas and found that SWMM is suitable for almost only urban areas, whereas HSPF can be applied to only homogenous land uses. Some studies have tried to improve upon existing models or to couple multiple models to form a new model; the resulting models can acquire better simulation results in diversified land use patterns. Lin et al. [30] assessed the effects of NPS phosphorus on soil by integrating the sediment delivery distributed (SEDD) and pollution load (PLOAD) models. Yang et al. [31] coupled the Xinanjiang model and SWAT to assess the NPS pollution load around Songtao Reservoir on Hainan Island. Nevertheless, land uses in urbanized areas encompass the characteristics of both urban and rural areas, and, thus, a variety of NPS types need to be considered. In addition, most of these models were developed in the United States or Europe; hence, considering the climate and soil distribution differences between these regions and China and other regions, some improvements must be implemented when using these models.
This study investigated a typical urbanized area, the Dafeng Basin. The study area has experienced rapid development with an urban population which has more than doubled in recent decades, but little attention has been paid to the water pollution therein, which seriously affects the quality of peopleʹs lives and hinders sustainable social and economic development within the basin. The most problematic pollutant is domestic sewage, although NPS pollution is also significant; both sources aggravate the existing burden on the water environment. In recent years, the local government has increased its investments in water pollution control and management; however, measured data from monitoring sites shows that despite the many measures that have already been implemented, the levels of TN and TP pollutants still do not meet the requirement of water quality standards (GB 3838-2002 standard in China). Since there is no TN standard for rivers in this standard, we use ammonia nitrogen standard to replace it. In this study, we evaluated the spatiotemporal dynamics of NPS pollution of the Dafeng Basin. On the basis of an evaluation of the previous literature, no model was found to be completely applicable to this type of region. Therefore, we combined multiple lumped models to develop an NPS assessment model that would be well suited to the study area. The main NPS pollutants in the study area were considered, including UR, AR, LPB, aquaculture and others, and their connections with land use were established. The combined model consists of two modules: water quantity and water quality, focusing on TN and TP, and considering the dynamic simulation of pollutant accumulation and erosion of pollutants under different land use conditions. Via the proposed model based on land use type, the spatiotemporal variation in the NPS pollution load were simulated and used to provide a scientific basis for regional water pollution control.

Study Area
The study area is a typical urbanized basin that covers an area of approximately 74.6 km 2 with a permanent population of about 185 thousand, of which has nearly 82.9% were urban residents in 2019. This area includes the urban area of Dafeng District and its surrounding suburbs, which border the Yellow Sea to the east ( Figure 1). The terrain is very flat; the elevation is relatively high, at 10.24 m, in the south and relatively low, at 8.79 m, in the north. The soil types in the study area are all moist soil, and most rivers have been channelized. According to the land use data in 2018 (Figure 2), the proportion of urban area (including urban land and highway) in the study area was 27.5%, while rural areas (including rural land and rural road) accounted for 19.5%, cultivated field (including irrigated land and paddy field) accounted for 30.5%, and other land use types accounted for 22.5%. The climate of this area, which is suitable for the growth of thermophilic crops, belongs to the transition zone between the subtropical and warm temperate zones. According to data acquired in the last 60 years, the average annual rainfall is 1075 mm, and the average annual evaporation was approximately 800 mm. The precipitation from June to September accounts for 63% of the precipitation throughout the whole year. In 2017, the average annual temperature was 15.7 °C, the daily maximum temperature was 38.6 °C and the daily minimum temperature was −6.8 °C. A truecolour image showing the distribution of sub-basins is provided in Figure 3.

Data Used in This Study
All input data employed for the NPS simulations included the following: (1) meteorological data from the period 2000 to 2018, including precipitation and evaporation; (2) vector data on land use types ( Figure 2) and hydrographic networks ( Figure 3); (3) watershed data obtained manually using digital elevation model (DEM) data and high-resolution imagery ( Figure 3); (4) soil type distributions and relevant parameters were from the China Soil Characteristics Dataset (CSCD) [32]; (5) discharge of industrial pollutants, release of endogenous pollutants, and domestic pollutants and loss of NPS pollutants obtained from on-line monitoring data, laboratory test data and field survey data, respectively; and (6) socioeconomic data obtained from local governmental statistical yearbooks. The data and their sources are summarized in Table 2.

Methods for Simulating NPS Pollutants
The combination model is a spatial distribution model ( Figure 4). The main modules include the water quantity module, which describes the processes of runoff yield and flow concentration, and the water quality module, which depicts the processes responsible for the generation of pollutants and their discharge into rivers.

Water Quantity Module
(1) Runoff yield Runoff was simulated according to different types of underlying surfaces. There were three main types: permeable soil, impervious surface and water surface.
For permeable soil, the grid-based Xinanjiang model was used to simulate runoff. The Xinanjiang model can be applied to runoff simulations for soil [33][34][35], but it is not suitable for impervious surfaces and water surfaces [31,35].
Impermeable surface does not consider infiltration, using the complete runoff yield method. The formula used is expressed as follows: where R is the runoff yield (mm); P is the rainfall amount (mm); K is the evapotranspiration conversion coefficient; and EI is the evaporation measured from the evaporating dish (mm). The water tank model, where runoff generation occurs only after the water tank is full, was used for water surfaces. The formula is as follows: where is the water storage at the end of time t. WM represents the saturated water storage, and its formula is as follows: where is the initial water storage and PE is the net rainfall (rainfall minus evaporation) at the end of time t.
Parameter values of the runoff yield simulation in the model are shown in the Table 3. (2) Flow concentration The flow concentration process consists of grid confluences and sub-basin confluences. The confluence of a grid is related to the grid size; when the grid size is small, the confluence time is very short, and the confluence process could be ignored. In this paper, the grid size is 30 m. Therefore, we assumed that the grid confluence is equal to its runoff yield. The sub-basin confluence was divided into two parts: the surface and subsurface confluences. The kinematic wave method was used to calculate the confluence time for the surface runoff [36], while the water tank model was used for the subsurface runoff. For each sub-basin, the cumulative time from each grid to the sub-basin outlet was calculated as follows: where ∆ is the duration when the flow of the i-th grid reaches the outlet of the corresponding subbasin (the value is an integer) (s); L is the ground flow length (m); n is the roughness (the values of different land use types are derived from some references and are shown in Table 4 [37,38]); is the net rain intensity (m/s); and S is the slope.
According to the above formula, the surface runoff in each period through the outlets of the subbasin was calculated as follows: , where is the surface runoff into the river (mm); is unit line of surface flow (m 3 s −1 ); I is the number of all grids in the sub-basin; A is the total area of the basin (km 2 ).

Water Quality Module
(1) Pollutant generation The pollution load generated during the rainfall process includes two states: the sorbed and the dissolved. In this paper, we simplified the exchange potential between sorbed and dissolved phases.
The sorbed pollution load originates predominantly from soil erosion.
where W is the generated sorbed NPS pollution load (tons); is the sediment load in the surface runoff (tons); and is the nitrogen or phosphorus nutrient content in the soil (g/g). is the enrichment factor, which was calculated as follows: where A is the grid area (hm 2 ); and is the surface runoff depth (mm). Soil erosion was calculated using the Modified Universal Soil Loss Equation (MUSLE).
where is the peak flow (m 3 s −1 ); K is the soil erodibility factor; is the vegetation cover and management factor; M is the maintenance measure factor; L is the topographic factor; and is the coarse debris factor. K and were calculated with soil data, whereas Ls was calculated with DEM and slope data. The determinations of and M were referred mainly to other research results and an actual investigation of the study area (Table 4).
In addition, the processes responsible for generating dissolved pollution of AR, UR and LPB are closely related to precipitation. Aquaculture pollution directly flows into the river, and thus has a minimal relationship with precipitation. The NPS pollution driven by precipitation was calculated using the cumulative wash-off model. The original cumulative wash-off model assumes that, during a rainfall period, the pollution load source is continuously attenuated during the rainfall process.
where SS indicates the source intensity of the field rainfall at the i-th moment (tons km −2 ); k is the ground wash coefficient, including the ground wash coefficients of TN and TP (kN and kP) ( Table 4); and is the cumulative runoff depth (mm) at the i-th moment of the rainfall event. As rainfall continues to fall on the field, the pollutant concentration decays exponentially with runoff depth.
Due to the influence of pollutant recharge, pollution sources accumulate. For urban area and LPB, the source of the grid pollutants will change at different times during the rainfall process. For farmland, the source of grid pollutants is replenished during both rainfall process and the dry (field drainage) period. Because the dominant pollution load is discharged during storm events, we ignore the dry period in this paper. The source of the pollutants at the i-th moment is therefore as follows: where 1 i l  is the pollutant replenishment source at time i-1 (tons km −2 h −1 ).
Then, the amount of pollution at the i-th time is as follows: Where Wdi is the dissolved pollution load generated at the i-th time. AR, UR and LPB were all calculated using this method, and the total amount of pollution was calculated as follows: where Wdis is the dissolved pollution load in a grid (tons); j is the type of NPS pollution (UR, AR and LPB); Wdi is the pollution load per unit area of type j at the i-th time (tons km −2 ).
(2) Pollutant emissions into a river The NPS pollution emissions entering a river are closely related to the confluence. For a single grid, the uniform water tank model was used to calculate the surface water and groundwater pollutant concentrations. For the t-th time (t 1 : where VS and VG are the surface water and groundwater storage capacities, respectively (mm); Q and Qs are the total runoff and surface runoff, respectively (mm); q and qs are the total outflow and surface discharge, respectively (m 3 s −1 ); CS is the concentration of TN or TP in the surface water (mg L −1 ); CG is the concentration of TN or TP in the groundwater (mg L −1 ); and Trj and Txf are the production of sorbed nitrogen or phosphorus and dissolved nitrogen or phosphorus, respectively (tons); A is the area of each sub-basin (km 2 ). The initial VS, VG, CS and CG were set in the model. Settling, volatilization and degradation occur during the migration and transformation of pollutants. According to the processing module of the distributed waste load model (DWLM) [39], a processing unit can be generalized into five types: sewage treatment plant, sewer, purification tank, surface water and soil (or underground) water. Drawing upon this idea, this study considered two basic types of surface water and groundwater. The treatment coefficients were adjusted according to the types of land use-that is, the purification rates of the pollutants.

FS
where T is the amount of pollutants into the river (tons); FS (including FNS and FPS of TN and TP respectively) is the treatment coefficient of surface water, and the treatment coefficients of the different land use types differed (Table 4). FG (including FNG and FPG of TN and TP respectively) was the groundwater treatment coefficient, and FNG and FPG were set to 0.85 and 0.95, respectively.

Rationality Analysis of the Simulation Results
(1) Simulation results of water quantity and quality Table 5 shows the results of the runoff simulation in the study area from 2000 to 2018. The average annual precipitation was 1100.8 mm. The soil water storage difference between the final and initial moment was 5.3 mm. ΔR represents the water balance error and the ratio of ΔR to P was less than 5%, which means that the simulation results satisfied the water balance. The total runoff coefficient was 0.51, and the surface runoff coefficient was 0.47. The surface runoff coefficient and total runoff coefficient both followed the same order: urban land > other > farmland. Annual average simulation results of NPS pollution production from the period 2000-2018 are shown in Table 6. From the perspective of soil erosion, the study area was only slightly eroded, and the amount of farmland erosion was the largest. These results are consistent with the spatial distribution data of soil erosion in China, which was provided by Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC) (http://www.resdc.cn). Dissolved nitrogen and phosphorus were the main pollution sources, accounting for 86.8% and 66.9% of the total, respectively. Compared with that of sorbed nitrogen, the proportion of sorbed phosphorus was relatively higher, reaching 33.1%. (2) Analysis of the rationality of the results in 2018 Table 7 shows the results of the various types of pollution loads in 2018. Among them, pollution load was calculated based on multi-source data. Domestic pollution load was obtained mainly from field survey data. The industrial pollution load was derived from on-line monitoring data provided by local government departments. The endogenous release pollution load was referred to laboratory test data; and the NPS load was simulated by the proposed combination model.
In the study area, the annual local water volume and the annual external water supply in 2018 were 4.71×10 7 m 3 and 4.21×10 6 m 3 , respectively. Based on the data and the results of all pollution loads emission into rivers, we estimated that the annual average of TN and TP concentration were 5.27 mg L −1 and 0.70 mg L −1 , respectively, in 2018. These estimated TN and TP concentrations were between the measured values and slightly above the median (Figure 5), and hence the estimated concentrations were considered reasonable.  (

3) Comparison of the simulation results between two methods
We used the export coefficient model, which is the most common method for estimating NPS pollution, to analyse the rationality of the results from the combination model. The annual average simulation results of NPS pollution emissions into rivers from the combination model and the export coefficient model were consistent regarding both the order of magnitude and the load distribution ( Figure 6). The TN and TP differences between the two results were -1.99 tons and −0.44 tons, respectively; their relative errors were less than 10% of the simulation results for the total amount and for each component.

Spatial Distribution of the Simulated NPS Loads
The annual average TN production per grid from the period 2000 to 2018 is shown in Figure 7a. The value ranged from 0 to 7.82 kg year −1 for a 30×30 m grid-that is, 0-8.80 tons km −2 year −1 . Ponds and farmlands (including irrigated lands, paddy field and orchards) had higher TN production intensities than other areas. The TN emissions from the four sources of pollution (UR, LPB, AR and other) for each sub-basin discharged into the river are shown in Figure 7b. The sub-basins with a high unit load being discharged into the river were located mainly in the north-western part of the study area. Among them, UR accounted for the largest proportion in most city-center sub-basins, and LPB had considerable impacts on some sub-basins in most suburban sub-basins. The NPS TP production per grid and emissions into the river each sub-basin are shown in Figures 8a, b. The TP production intensity was prominent in ponds and some farmland grids. The maximum value was 2.57 kg year −1 for a 30×30 m grid-that is, 2.85 tons km −2 year −1 . The total TP emissions and unit load for the north-western sub-basins were both higher than those of the other sub-basins. UR accounted for the largest proportion among these sub-basins. Aquaculture had relatively large impacts on sub-basins 2 and 7.

Temporal Variations in the Simulated NPS Loads
(1) NPS loads during rainfall events Figure 9 shows the NPS TN and TP pollution loads resulting from two rainfall events. The sorbed and dissolved TN and TP were both consistent with the rainfall. The rainfall intensity in May was much greater than that in January. Dissolved nitrogen and phosphorus dominated when the rainfall intensity was low; as the rainfall intensity increased, the proportions of sorbed nitrogen and phosphorus increased, especially that of sorbed phosphorus. (2) Seasonal variation in NPS pollution The NPS TN and TP loads into rivers displayed obvious differences on the monthly scale ( Figure  10). High TN and TP emissions occurred mainly in summer (June, July and August), representing 56.2% and 58.0% of the total TN and TP loads, respectively, especially in July. Although there was more precipitation in August, the TN emissions were lower in August than in June. Overall, the relationship between TN emissions and rainfall was consistent with that between TP emissions and rainfall.  Figure 11 shows the monthly average concentrations of the NPS TN and TP emissions into rivers in each sub-basin without considering point and internal sources. The area-averaged concentrations of TN and TP were relatively low from July to October. The maximum concentrations of TN and TP were 1.24 mg L −1 and 0.19 mg L −1 , respectively (Figure 11(a)). Sub-basins 33 and 37 represented the concentrations of typical urban area and farmland, respectively (Figure 11b, c). The monthly magnitude and its variation in the concentration of farmland were both higher than that of urban area. The typical farmland was characterized by a poorer water quality from February to April and a better water quality from July to October. This phenomenon was related to the pollution caused by livestock and poultry breeding in the farmland sub-basin.

Comparison with Other Studies
Water pollution has received widespread attention, and many studies have explored the spatial and temporal variations in sediment and pollutants. Table 8 shows the simulation results of some of these studies from China. The simulated results of the sediment, TN and TP loads in our study area are comparable to the results reported in these case studies. In particular, the sediment load was above 30 tons km −2 year −1 in agricultural areas, e.g., the Sanjiang Plain, Hainan Island and our research area. However, there were some differences in the TN loads of different regions; nevertheless, the loads were generally greater than 0.3 tons km −2 year −1 in urban areas and agricultural areas. Additionally, the TP loads of different regions varied greatly; in some regions, the TP loads was as low as 0.04 tons km −2 year −1 , while the TP load reached 1.7 tons km −2 year −1 in other regions.

Suggestions for Water Pollution Control
This study calculated various types of pollution loads. Among them, domestic pollution was the main source of pollutants, followed by NPS pollution. However, domestic pollution can be addressed by building pipe networks and corresponding sewage treatment plants. In contrast, NPS pollution is relatively difficult to control.
(1) Pollution load control priority According to the simulation results of NPS pollution emissions, we prioritized the processing of pollution loads. The TN emission intensity decreased in the order of UR > LPB > AR > other, and the TP emission intensity decreased in the order UR > other > LPB > AR (Table 7 and Figure 10). Evidently, UR and LPB were the main sources of TN and TP emissions. Therefore, we should first consider reducing UR and LPB. To reduce the UR pollution load, it is fundamental to improve the urban ecological environment by increasing the urban greening rate and constructing sponge cities [42]. For LPB, the management strategy could be further optimized by returning waste to farmlands as much as possible, based on the pollution load carrying capacity of the land [43]. In addition, farmland strategies require improving the utilization rate of fertilizers [21] to reduce the pollution emissions. For aquaculture, some aquatic plants can be planted to absorb NPS pollutants [19].
The concentrations of NPS TN and TP in rivers were simulated ( Figure 11). The concentrations of TN in sub-basin 37 exceeded the fourth category standard during the dry season, while that of TP below the same standard all year round (Figure 11c). The load reduction pressure of TN was greater than that of TP; hence, for farmlands, the focus should be on reducing TN. In addition, the concentrations of TN and TP emissions were relatively high during the dry season. Thus, the reduction of pollutants during the dry season should be of particular concern.
(2) Measures of water quality improvement Measures to improve water quality mainly include reducing pollution emissions and diluting pollution loads. The water quality compliance requirements in the study area are fourth categorythat is, TN < 1.5 mg L −1 and TP < 0.3 mg L −1 . Based on the discharged pollutant loads in 2018, it is necessary to reduce TN by at least 193.54 tons and TP by at least 20.39 tons to meet these standards. If we adjust the supply water instead of reducing pollutants, at least 1.29×10 8 m 3 of pure water will need to be supplied to meet this standard. If the water used to dilute the pollutants satisfies the third category standard, it will be necessary to supply at least 3.87×10 8 m 3 of water. Among them, the amount of water needed to dilute NPS pollution with pure water or the third category standard water accounted for 16.9%. Certainly, this issue is serious and very common in the control of polluted rivers in many cities throughout China.

Uncertainties and Limitations
The study area is a man-made river basin, the boundaries of which are controlled by water conservancy engineering, and the river flow is controlled artificially. Unfortunately, only observation of the water level and not the flow volume are available. Therefore, a comparison of the simulated runoff with measured data cannot be provided. Here, as an alternative, we performed water balance analysis and the rationality analysis of the runoff coefficient to demonstrate the rationality of the simulation results. For such a rain-sewage confluence area, many sewage outlets lack measured (let alone continuously monitored) data, which makes it very difficult to verify the simulated results directly and therefore introduces some uncertainties. In this paper, an on-site investigation was conducted to mine available information, including the layout of sewage networks and their operation situations, the numbers of the overflow sewage outlets, and the discharges of industrial wastewater. The simulated NPS pollution results, in addition to the surveys of domestic pollution sources, online monitoring of industrial pollution sources, and laboratory tests of sediment contamination, were compared with the measured water quality. The simulated average concentration was consistent with the measured water quality. In addition, we compared the simulated results with those of the export coefficient model; their differences were less than 10% for the simulated total amount and for each component. Although it is impossible to provide further validation, these analyses still represent progress for a developing city with such a complex situation and limited data.
The four sources identified in this paper, namely UR, LPB, AR and other sources, were the four main NPS pollutants in the study area. The source data employed in this paper were obtained mainly from reference materials and local investigations, and some of these data were at monthly or annual scales. Although we considered the pollution load factors and made some adjustments for the pollution sources at different times, some uncertainties in the hourly and daily data remained. Of course, the temporal uncertainty in source data has certain impacts on hourly and daily simulated results but has little effect on monthly and annual-scale results. Therefore, the monthly and annualscale simulation results are reliable.
Many parameters were involved in the model. Among them, the runoff generation module of the Xinanjiang model included 13 parameters (Table 3), which were mainly set based on the relevant literature. Some parameters significantly affected the water quantity simulation. During the parameter debugging process of the Xinanjiang model, we found that W0, K and IMP were relatively sensitive. Table 4 shows several important parameters in the water quality simulation process. The vegetation cover and management factor and the maintenance measure factor were both key parameters for calculating the sediment and sorbed pollution production. Moreover, the ground wash coefficients of TN and TP (kN and kP) had considerable influences on the dissolved pollution load. Finally, the surface water treatment coefficients (FNS and FPS) and groundwater treatment coefficients (FNG and FPG) affected the total amounts of simulated NPS loads into rivers.

Conclusions
In this paper, a combination model was proposed and used to simulate the NPS TN and TP pollution load in the study area. We evaluated the spatiotemporal distribution of NPS pollution load. According to the results, we discovered the following: (1) Ponds and farmlands had higher TN and TP production intensities than the other land use types, and the unit emissions in the northwestern region, which is mainly urban land and contains a great deal of aquaculture and LPB, were relatively higher; (2) The variations in the monthly and interannual TN and TP loadings were consistent with the variations in rainfall. The emissions of TN and TP accounted for 56.2% and 58.0%, respectively, of the total in summer; (3) NPS TN pollution was more serious than NPS TP pollution in the study area, especially in farmlands, and the concentrations of TN and TP emissions were relatively higher during the dry season; (4) UR and LPB were the main sources of NPS TN and TP emissions in the study area. If NPS pollution cannot be removed from the study area, at least 2.19×10 7 m 3 of pure water or 6.56×10 7 m 3 of the third category standard water needs to be supplied to dilute the current rivers to meet the required fourth category standard. Funding: This research was supported by the National Natural Science Foundation of China, "Sponge cities within managed basins: integrated responses of hydrology and environmental quality to constructing a sustainable city" (Project No. 71961137007); and the special fund of basic scientific research operating expenses of Ministry of Water Resources of the People's Republic Of China, "Research on the comprehensive management model and key technologies of regional water problems" (Project No. 519013).

Conflicts of Interest:
The authors declare no conflicts of interest.