Spatial and Seasonal Dynamics of Water Environmental Capacity in Mountainous Rivers of the Southeastern Coast, China

The south-east littoral is one of the most populous and developed regions in China suffering from serious water pollution problems, and the Xian-Jiang Basin in the mid of this region is among the most polluted watersheds. Critical information is needed but lacking for improved pollution control and water quality assessment, among which water environmental capacity (WEC) is the most important variable but is difficult to calculate. In this study, a one-dimensional water quality model combined with a matrix calculation algorithm was first developed and calibrated with in-situ observations in the Xian-Jiang basin. Then, the model was applied to analyze the spatial and temporal patterns of WEC of the entire basin. The results indicated that, in 2015, the total pollutant discharges into the river reached 6719.68 t/yr, 488.12 t/yr, and 128.57 t/yr for COD, NH3-N and TP, respectively. The spatial pattern suggested a strong correlation between these water contaminants and industrial enterprises, residential areas, and land-use types in the basin. Furthermore, it was noticed that there was a significant seasonal pattern in WEC that the dry season pollution is much greater than that in the plum season, while that in the typhoon season appears to be the weakest among all seasons. The WEC differed significantly among the 24 sub-basins during the dry season but varied to a smaller extent in other seasons, suggesting differential complex spatial-temporal dependency of the WEC.


Introduction
Coastal regions, where about 70% of the population and industrial capitals worldwide, are full of well-developed metropolitan economic corridors, leading to serious water shortages and environment deteriorations, due to rapid urbanization and industrialization [1,2]. In recent years, total water pollution control as a policyinstrument has been adopted by environmental protection agencies in China to meet the national water quality standards [3] and to ensure drinking water safety. As the premise and foundation of water pollution control, knowledge of water environmental capacity (WEC) plays a critical role [4][5][6]. Like total maximum daily load (TMDL), the concept of WEC refers to the total maximum load of pollutants that a waterbody can accommodate without significantly affecting its environmental function and water quality standards set forth by environmental agencies [7][8][9]. Coastal climatological and orographic characteristics where mountainous land and plain crisscross (influenced by a marine monsoon climate), lead to an obvious dissimilarity in rain-runoff in space and time [10,11] that is likely to result in different WEC. However, existing WEC modeling is insufficient

Study Area
The study area (Xian-Jiang Basin) is located in Ningbo, Zhejiang Province in the southeast of China, with a drainage area of 306.70 km 2 of 15 main rivers ( Figure 1). The area has a subtropical monsoon climate that encompasses a dry season from December to January and a plum rain season from May to June. The time period from July to September is the typhoon season, when rains are mixed with typhoon storms [21]. The basin's terrain transitions from mountain topography in the west to the gentle plain in the coastal east. The basin provides an important source of drinking water for Ningbo City.
In spite of its importance as Ningbo's primary drinking water source, the basin is among the most industrialized cities in China, experiencing serious freshwater pollution problems over the past three decades that threatens water environment and drinking water security [22,23]. The primary water pollutants found in the water are chemical oxygen demand (COD), NH 3 -N, and total phosphorus (TP). The poor water quality was attributed to the industrial plants, residential sewage, agricultural run-off, large-scale livestock culture, and wastewater treatment plant discharges (Figure 1) within the basin.

Hydrological Characteristics of Xian-Jiang Basin
The average monthly rainfall over the last 15 years (2000-2014) recorded at the three rain gauge stations in the Hong-Jia-Ta catchment with the area of 151 km 2 near the basin shows a bimodal pattern ( Figure 2). The monthly average precipitation over the last 15 years ranged from 50 mm in the winter season to 250 mm in the summer season. The heaviest rainfall (more than 248 mm) occurred in June and August; July and September also received high rainfall, while the minimum rainfall occurred in December with only 63.53 mm. Rainfall intensity in the dry season, plum rain and typhoon seasons were 129.42 mm, 346.76 mm, and 604.99 mm respectively, accounting for 8.42%, 22.56%, and 39.36% of the annual total precipitation. Long-term hydrological observations at Hong-Jia-Ta hydrological station in Ningbo City were used as the monthly run-off data (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) in the Hong-Jia-Ta catchment, as it is the only station

Hydrological Characteristics of Xian-Jiang Basin
The average monthly rainfall over the last 15 years (2000-2014) recorded at the three rain gauge stations in the Hong-Jia-Ta catchment with the area of 151 km 2 near the basin shows a bimodal pattern ( Figure 2). The monthly average precipitation over the last 15 years ranged from 50 mm in the winter season to 250 mm in the summer season. The heaviest rainfall (more than 248 mm) occurred in June and August; July and September also received high rainfall, while the minimum rainfall occurred in December with only 63.53 mm. Rainfall intensity in the dry season, plum rain and typhoon seasons were 129.42 mm, 346.76 mm, and 604.99 mm respectively, accounting for 8.42%, 22.56%, and 39.36% of the annual total precipitation.

Hydrological Characteristics of Xian-Jiang Basin
The average monthly rainfall over the last 15 years (2000-2014) recorded at the three rain gauge stations in the Hong-Jia-Ta catchment with the area of 151 km 2 near the basin shows a bimodal pattern ( Figure 2). The monthly average precipitation over the last 15 years ranged from 50 mm in the winter season to 250 mm in the summer season. The heaviest rainfall (more than 248 mm) occurred in June and August; July and September also received high rainfall, while the minimum rainfall occurred in December with only 63.53 mm. Rainfall intensity in the dry season, plum rain and typhoon seasons were 129.42 mm, 346.76 mm, and 604.99 mm respectively, accounting for 8.42%, 22.56%, and 39.36% of the annual total precipitation. Long-term hydrological observations at Hong-Jia-Ta hydrological station in Ningbo City were used as the monthly run-off data (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) in the Hong-Jia-Ta catchment, as it is the only station Long-term hydrological observations at Hong-Jia-Ta hydrological station in Ningbo City were used as the monthly run-off data (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) in the Hong-Jia-Ta catchment, as it is the only station closest to the basin. Both rainfall and runoff appear to have a similar pattern (Figure 3a), and there is a high correlation that can be fit with an exponential function (Figure 3b) [24]. The explanation power or correlation coefficient is 0.87, which suggests that runoff can be approximated by its rainfall patterns. The streamflow in the Xian-Jiang basin and its sub-basins can, therefore, be estimated using precipitation data (2015) from Feng-Hua station within the basin through the hydrologic analogy method (Equation (1)) [25].
where Y i is the ith sub-basin's monthly runoff (m 3 /s); k i is the area of ith sub-catchment (km 2 ); K is the area of Hong-Jia-Ta catchment; x is the monthly precipitation of the basin in 2015 (mm). closest to the basin. Both rainfall and runoff appear to have a similar pattern (Figure 3a), and there is a high correlation that can be fit with an exponential function ( Figure 3b) [24]. The explanation power or correlation coefficient is 0.87, which suggests that runoff can be approximated by its rainfall patterns. The streamflow in the Xian-Jiang basin and its sub-basins can, therefore, be estimated using precipitation data (2015) from Feng-Hua station within the basin through the hydrologic analogy method (Equation (1)) [25].
where is the th sub-basin's monthly runoff (m 3 /s); is the area of th sub-catchment (km 2 ); is the area of Hong-Jia-Ta catchment; is the monthly precipitation of the basin in 2015 (mm).

Calculation and Allocation of Pollution Loads
The pollution sources in the basin were identified mostly from agricultural non-point sources, large-scale breeding farms, domestic sewage, industrial enterprises, and a wastewater treatment plant. The load from the first two sources was computed by using Johnes model (Equation (2)) [26]. The rate loss of various agricultural non-point sources loads (paddy field, dry land, economic forest, and garden) was determined based on previous studies citing analogy basins [27][28][29]. For the domestic sewage, there are primarily three mechanisms to discharge pollutants into rivers, according to a field survey and a review of previous studies: (1) sewage pipes divert effluence into rivers and streams after a treatment by sewage disposal facilities in villages; (2) the effluent is discharged directly into rivers and streams after a simple treatment by septic tanks; and (3) all others are all discharged into a sewage network and eventually enter sewage treatment plants before being discharged into streams and rivers.
The pollution loads for each of the three were calculated in this study. The pollution loads of 36 key industrial enterprises and wastewater treatment plants in the basin were obtained from the environmental statistical data of Ningbo city in 2015. As for pollution loads allocation into rivers, the pollution load from plants (including the wastewater treatment plant) discharged directly to the reaches where their outlets lay.
The load of agricultural sources was apportioned equally over reaches based on the spatial linkage between land-use/land-cover (LULC) ( Figure 4) and river systems.
The domestic source loads were allocated to the corresponding reaches through spatial correlations between residential districts and reaches.
The loads from large-scaled farming were allocated into the appropriate reach which were geographically closest to the farms based on longitude and latitude coordinates.

Calculation and Allocation of Pollution Loads
The pollution sources in the basin were identified mostly from agricultural non-point sources, large-scale breeding farms, domestic sewage, industrial enterprises, and a wastewater treatment plant. The load from the first two sources was computed by using Johnes model (Equation (2)) [26]. The rate loss of various agricultural non-point sources loads (paddy field, dry land, economic forest, and garden) was determined based on previous studies citing analogy basins [27][28][29]. For the domestic sewage, there are primarily three mechanisms to discharge pollutants into rivers, according to a field survey and a review of previous studies: (1) sewage pipes divert effluence into rivers and streams after a treatment by sewage disposal facilities in villages; (2) the effluent is discharged directly into rivers and streams after a simple treatment by septic tanks; and (3) all others are all discharged into a sewage network and eventually enter sewage treatment plants before being discharged into streams and rivers.
The pollution loads for each of the three were calculated in this study. The pollution loads of 36 key industrial enterprises and wastewater treatment plants in the basin were obtained from the environmental statistical data of Ningbo city in 2015. As for pollution loads allocation into rivers, the pollution load from plants (including the wastewater treatment plant) discharged directly to the reaches where their outlets lay.
The load of agricultural sources was apportioned equally over reaches based on the spatial linkage between land-use/land-cover (LULC) ( Figure 4) and river systems.
The domestic source loads were allocated to the corresponding reaches through spatial correlations between residential districts and reaches.
The loads from large-scaled farming were allocated into the appropriate reach which were geographically closest to the farms based on longitude and latitude coordinates.
Therefore, the total contaminants discharged to the river and stream systems are calculated as the summation of all components.
where L j is the discharge of contaminant j in the basin (kg); E ij is the export coefficient of the ith land-use type (kg/hm 2 ) or the number of the i th poultry for contaminant j; A i is the area of the ith land-use type (hm 2 ) or the discharge coefficient of the i th poultry (kg/yr). Therefore, the total contaminants discharged to the river and stream systems are calculated as the summation of all components.
where is the discharge of contaminant in the basin (kg); is the export coefficient of the th land-use type (kg/hm 2 ) or the number of the th poultry for contaminant ; is the area of the th land-use type (hm 2 ) or the discharge coefficient of the th poultry (kg/yr).

Determination of Stream Segmentation and Hydraulic Characteristics
For advantageous analysis of WEC, the watershed area was divided into 25 sub-basins (S01-S25). The hydrologic network was extracted from the digital elevation model (DEM) data, using the function of hydrologic analysis in ArcGIS 10.2 (version 10.2, ESRI Inc., Redlands, CA, USA), where there was one main stream (Xian River) and fourteen tributaries ( Figure 5). Each reach was then subdivided into uniform computational elements for the model. These were 100 m in size, almost ten times more than the average calculation precision, and the reaches were divided into 1350 segmentations in total. The properties of the segmentations (length, width, sides slope, channel gradient, and elevation) were measured or estimated using ArcMap 10.2 (version 10.2, ESRI Inc., Redlands, CA, USA). The serial numbers of sub-basins and computational elements are displayed in Figure 5.

Determination of Stream Segmentation and Hydraulic Characteristics
For advantageous analysis of WEC, the watershed area was divided into 25 sub-basins (S01-S25). The hydrologic network was extracted from the digital elevation model (DEM) data, using the function of hydrologic analysis in ArcGIS 10.2 (version 10.2, ESRI Inc., Redlands, CA, USA), where there was one main stream (Xian River) and fourteen tributaries ( Figure 5). Each reach was then subdivided into uniform computational elements for the model. These were 100 m in size, almost ten times more than the average calculation precision, and the reaches were divided into 1350 segmentations in total. The properties of the segmentations (length, width, sides slope, channel gradient, and elevation) were

Hydrodynamic Process
Every-segmentation in the reaches can be conceptualized as a trapezoidal channel with steady flow. For each channel, the nature of the relationship between flow and water depth was described by the Manning equation, expressed as Equation (3) [30]. To calculate the depth, flow rate, bottom width, riverbed slopes and side slope factors of the elements were required. Bottom width, riverbed slopes, and side slope factors were extracted from the DEM and vector data of the watershed. The flow rate in each element was calculated using the data from the above-mentioned hydrologic analysis. The flow of the 875th element located at the outlet of the reservoir was set based on the monthly average values of the discharge of the reservoir. According to previous studies on the Manning coefficient of roughness, the suggested n is generally in the range of 0.015-0.15, depending on the various types of channels [31,32]. In this paper, empirical n was set as 0.04 based on the properties of the river channels in the basin.
where, H = depth of the river (m); = flow of the river (m 3 /s); = Manning roughness coefficient; = bottom width of the river (m); = riverbed slope ( o ); and , represent both side slope factors, respectively. = 1, 2, ⋯ , , is the number of iterations. When the estimation error

Hydrodynamic Process
Every-segmentation in the reaches can be conceptualized as a trapezoidal channel with steady flow. For each channel, the nature of the relationship between flow and water depth was described by the Manning equation, expressed as Equation (3) [30]. To calculate the depth, flow rate, bottom width, riverbed slopes and side slope factors of the elements were required. Bottom width, riverbed slopes, and side slope factors were extracted from the DEM and vector data of the watershed. The flow rate in each element was calculated using the data from the above-mentioned hydrologic analysis. The flow of the 875th element located at the outlet of the reservoir was set based on the monthly average values of the discharge of the reservoir. According to previous studies on the Manning coefficient of roughness, the suggested n is generally in the range of 0.015-0.15, depending on the various types of channels [31,32]. In this paper, empirical n was set as 0.04 based on the properties of the river channels in the basin. where, H = depth of the river (m); Q = flow of the river (m 3 /s); n = Manning roughness coefficient; B 0 = bottom width of the river (m); S = riverbed slope ( o ); and S S1 , S S2 represent both side slope factors, respectively. k = 1, 2, · · · , i, i is the number of iterations. When the estimation error ε α × 100%) < 0.001%, the iteration process stops. Having obtained the water depth of each element, it was then substituted into Equation (4) to calculate the cross-sectional area of the channel (m 2 ): Finally, the cross-sectional area A C (m 2 ) was applied to Equation (5) to obtain the flow velocity U (m/s) of each element:

Pollutant Degradation Process
In this paper, a one-dimensional water quality model was chosen to simulate pollutant degradation in the rivers. Based on water balance relationships, the fundamental equation for water quality simulation (Equation (6)) can be written as described in Reference [33]: The rapidly flowing rivers are dominant in mountain areas due to elevation differences that transverse velocity gradient which is much greater than vertical velocity gradient in rivers. Effects of longitudinal dispersion of contaminant are negligible compared to advection (the longitudinal dispersion coefficient E = 0) [34][35][36], then the model can be expressed as: where c is the concentration of the water quality (mg/L); u is the flow velocity (m/s) that can be calculated using the Manning equation; t is the time of river water flowing through from the headwater to somewhere (s); x is the distance that river water flows through in time t (m), x(t) = ut, and k 1 is the pollutant degradation coefficient (d −1 ). For initial conditions of x(t) = 0 and c = c 0 , the accurate solution of Equation (7) can be obtained as follows: When the initial concentration (x = 0, c = c 0 ) is known (the headwater referred to the water quality standard of water source protection area in China; the reservoir outlet referred to the actual monitoring concentrations), the concentration of a certain pollutant at x distance can be calculated using Equation (8).

Matrix Cacluation Algorithms
The process of little and nonpoint source input and degradation in the river is illustrated in Figure 6. In this figure, i represents the ith calculated element, i =1,2, · · · , n; Q 10 , and C 10 represent the initial flow rate and concentrations from the headwater; Q 1i , and C 1i represent the incoming flow rate and concentrations from the i − 1th calculated element; Q i , and C i represent the rainfall-runoff and pollutant concentrations draining into the ith calculated element; Q 2i , and C 2i represent the outflow and concentrations of effluent discharging from the ith calculated element into the downstream; k 1i represents the pollutant degradation coefficient between i and i + 1; and t i represents the time consuming of flowing from i to i + 1. the initial flow rate and concentrations from the headwater; , and represent the incoming flow rate and concentrations from the − 1 th calculated element; , and represent the rainfall-runoff and pollutant concentrations draining into the th calculated element; , and represent the outflow and concentrations of effluent discharging from the th calculated element into the downstream; represents the pollutant degradation coefficient between and + 1; and represents the time consuming of flowing from to + 1. Figure 6. Schematic of conceptual model of one-dimensional steady-state river. Figure 6. Schematic of conceptual model of one-dimensional steady-state river.
Based on solving the Streeter-Phelps equation, the changes in contamination concentration due to attenuation from i − 1 to the i element can be described using Equation (9): According to the principle of flow continuity, the balance of flow and concentrations for each node can be drawn using the following equations [33]: Then, Equations (9) and (10) can be derived using the following matrix equation: where Having obtained the needed computational parameters, the little and nonpoint pollutant concentrations with C as the input, and the corresponding pollutant concentrations C 2 as the output for each calculated element could be obtained. The computed pollutant concentrations of each element were applied to Equations (15) and (16) (below) to assess the water environmental capacity of each element.

Water Environmental Capacity Calculation
To calculate the WEC of each element with limited data, a simple and proper computation model is needed. In this study, a section control method (that has been one of the most frequently used algorithms) was adopted that can be expressed as Equation (11) [37], where q l (average discharge flow from the lth outfall) was far less than Q (reach designed flow of element i). The simplified model can be expressed using Equation (12).
where W i , C Si , and C 0i represent the water environmental capacity, the target concentration of the water quality, and the actual concentration of the pollutant of element i, respectively; k i is the degradation coefficient of the pollutant in element i, and t i is the time consuming of flowing through the element i. Eventually, the contamination concentration and WEC of each calculated element i was simulated through the one-dimensional pollutant-water response model mentioned above. The calculation procedure of all models was programmed using MATLAB R2012b (version 8.0, The MathWorks, Natick, MA, USA).

Model Calibration and Verification
Water quality data (COD, NH 3 -N, and TP) of the monitoring section in the basin collected by the Zhejiang Province environmental monitoring station were used to calibrate and validate the water quality model at the beginning of the study. To minimize errors between the simulation results and observed data, the constants of COD, NH 3 -N, and TP degradation (k 1 ) were appropriately adjusted by trial and error, and were 0.1 day −1 , 0.08 day −1 , and 0.08 day −1 , respectively [38,39]. The measured data inserted into the model for calibration and verification were monthly values (the mean of three water samples collecting in a month) from January to December 2015. Data from single months were used for calibration, and data from double months were used to verify the model. A comparison of the calibration results and the field measurements, in general, agreed well with the monitoring data, with the exception of a few values. For example, the simulated results of July and August deviated from the observed values most likely due to the delay of gathering water samples in a torrential downpour period. The root mean square error (RMSE) of COD, NH 3 -N, and TP between the simulated and measured values in the rest of the months of 2015 (excluding July and Augest) were 6.28, 0.39, and 0.08, respectively, and met the requirement of precision. Therefore, the models were acceptable to simulate the water quality and WEC under the condition of limited data. It should be noted that it does not fully validate a watershed model [18].

Contaminant Loadings Analyses
The pollutant quantity inlets into the river directly relates to the degree of contamination in the Xian-Jiang Basin. According to the calculation results, the total emissions into the river of COD, NH 3 -N, and TP in 2015 were 6719.68 t/a, 488.12 t/a, and 128.57 t/a for the Xian-Jiang Basin, respectively. Furthermore, the total contamination loads into the river in 2015 based on the type of pollution source were calculated (Table 1). COD discharge was due to agricultural nonpoint sources, formalization cultivation, domestic sewage, industrial enterprises above the designated size and wastewater treatment plant were 1188.78 t/a, 204.25 t/a, 4519.20 t/a, 65.68 t/a, and 741.77 t/a respectively, and accounted for 17.69%, 3.04%, 67.25%, 0.98%, and 11.04% of the total COD discharge in the watershed, respectively. Domestic sewage and agricultural nonpoint sources have been identified as the most significant sources of COD pollution, while industrial enterprises above the designated size and formalization cultivation were the least significant sources. The discharge and proportion of domestic sewage, agricultural nonpoint sources, wastewater treatment plants, formalization cultivation, and industrial enterprises above the designated size for NH 3 -N were 388.96 t/a (79.69%), >46.03 t/a (9.43%), >26.33 t/a (5.39%), >22.06 t/a (4.52%), and >4.74 t/a (0.97%). The domestic sewage (76.01 t/a) and agricultural nonpoint sources (24.86 t/a) were the two biggest contributors for TP in the watershed, which accounted for 59.02% and 19.34% of the aggregate, respectively. The discharge and proportion of wastewater treatment plant and formalization cultivation were 15.04 t/a (11.70%) and 12.35 t/a (9.61%), respectively. Table 1. The contamination loads and proportion of pollution sources into the basin.

The Spatio-Temporal Analyses of Water Quality
With the Manning formula-SP model simulation, the contaminant concentrations of the river system in the basin were obtained. Figures 7-9 present the spatial distribution of water quality standards of COD, NH 3 -N, and TP in three different seasons, respectively. As shown, the water quality of the pollutants exhibited a seasonal variability, which in the dry season was demonstrably worse than in other seasons. The water quality during the plum rain season and typhoon season was in a comparable condition, and both were superior to the water quality objective in the majority of the reaches. This is because during the dry season, there is much less precipitation than in other seasons, and the main contaminative source is domestic sewage in the basin leading to a small water volume containing masses of pollutants. For the plum rain and typhoon seasons, the volumes were both large enough to digest the pollution loads. Although the plum rain season has more precipitation than the typhoon season, at the same time, runoff and drainage carry more nonpoint source pollutants into the river. The seriously polluted reaches that were Grade III or even worse have similar spatial distribution between plum and typhoon seasons. Especially, for COD, the number of the reaches in plum season were slightly more than in typhoon season in the upstream watershed of the reservoir. Regarding the Xi-Xi River in the middle area of the lower reaches, the water quality in plum season was worse than in typhoon season in the same reaches.
Selecting the most inferior dry season as an example, water quality in the upper reaches of the reservoir was generally better than in the lower reaches, where there is a water source protection area and a few factories and settlements. Among the three kinds of pollutants, COD had the best status upstream of the reservoir, the concentration of which basically met the Grade II standards, with the exception of the two tributaries near the reservoir. However, for NH 3 -N and TP, few tributaries in the upper reaches were Grade III or even worse, where a large number of paddy fields were distributed on both sides (compared with the LULC in 2015) carried plenty of nitrogen and phosphorous into the river. Furthermore, the water quality of TP was more deteriorated than NH 3 -N in the same reaches, partly due to the excessive use of phosphatic fertilizer in China, and partly because phosphorous is usually considered as a limiting nutrient for aquatic bioactivity and is largely removed through biodegradation or sedimentation [19]. The massive flow rate in this area prevented TP degradation and settlement. Therefore, the rice area along the river in the protection zone must be replaced by eco-friendly types. The concentrations of COD, NH 3 -N, and TP increased acutely in the middle area of the lower reaches where industry and resident cluster areas are located. A focus on the water quality in this area showed that the concentration of COD in the main river (Xian River) was much better than that of the tributaries due to water yield. This was especially the case for the Xi-Xi River where almost a whole section of the river was worse than Grade V, which was mainly attributed to industrial wastewater discharge from many machinery manufacturing plants around the river containing heavy concentrations of COD, besides the municipal loads. The concentration of NH 3 -N in this area wholly violated the water quality standard of Grade III including the main river. The top three heaviest polluted reaches of this area included outlet reaches where the river water flows into the Xian River and is the local center of the residential area, and the middle of the Bei-Xi River, as well as the upper reach of the Xi-Xi River where intensive machinery manufacturing plants are located. The situation of TP was similar to NH 3 -N, which corresponded to the same location as the most seriously polluted area. To summarize, a compelling spatial correlation existed between the density of pollutants and the distribution of plants and settlements.

The Temporal Variability of WEC in the Basin
The WEC of 24 sub-basins (excluding the S15 where the reservoir is located) as a group of subjects has been calculated for different hydrological seasons (dry season, plum rain season, and typhoon season), respectively. If the calculations of WEC were positive, it means the environmental capacity remained in surplus and could still accommodate the pollution load into the river with the premise of not violating the water quality standard. The positive values were replaced by zero indicating that no load needed to be cut; negative values indicated that the pollution load exceeded the environmental capacity and were replaced by its absolute values standing for the quantity that needed to be reduced. The independent-samples nonparametric statistical analyses with Mann-Whitney U tests were performed to assess group differences of the reduced WEC among the three hydrological seasons (Figure 10

The Temporal Variability of WEC in the Basin
The WEC of 24 sub-basins (excluding the S15 where the reservoir is located) as a group of subjects has been calculated for different hydrological seasons (dry season, plum rain season, and typhoon season), respectively. If the calculations of WEC were positive, it means the environmental capacity remained in surplus and could still accommodate the pollution load into the river with the premise of not violating the water quality standard. The positive values were replaced by zero indicating that no load needed to be cut; negative values indicated that the pollution load exceeded the environmental capacity and were replaced by its absolute values standing for the quantity that needed to be reduced. The independent-samples nonparametric statistical analyses with Mann-Whitney U tests were performed to assess group differences of the reduced WEC among the three hydrological seasons (Figure 10 As seen, the average reduced WEC of COD in the dry season was nearly eight and 15 times higher than in the plum rain and typhoon seasons, respectively (Mann-Whitney U-test, p < 0.001, As seen, the average reduced WEC of COD in the dry season was nearly eight and 15 times higher than in the plum rain and typhoon seasons, respectively (Mann-Whitney U-test, p < 0.001, one-sided, n = 24), indicating that there was a significant seasonal difference on COD gross reduction in the basin. The COD reduction had no obvious difference between the plum rain and typhoon seasons (Mann-Whitney U-test, p = 0.673, one-sided, n = 24). The dry season group also revealed a striking difference in reduced WEC of NH 3 -N with the plum rain and typhoon season groups (Mann-Whitney U-test, p < 0.001, one-sided, n = 24). Moreover, there was still no significant difference in a comparison of the plum rain season with the typhoon season (Mann-Whitney U-test, p = 0.572, one-sided, n = 24). As for TP, the distribution of the reduced WEC of the plum rain and typhoon season groups were far below that of the dry season group (Mann-Whitney U-test, p < 0.001, one-sided, n = 24). The average reduction of the dry season group (4.44 t/a) was 4.11 and 5.77 times that of the plum rain and typhoon season groups (plum season = 1.08 t/a; typhoon season = 0.77 t/a). No significant change was detected for the reduced WEC of TP between the plum rain and typhoon season groups (Mann-Whitney U-test, p = 0.419, one-sided, n = 24). Thus, the differences of reduced WEC in the basin for the three contaminants between the dry season and the other seasons were highly significant, suggesting that the reduced WEC exhibited a temporal variability that the pollution loads in the dry season were much more serious than other seasons in the study area.

The spatial Variability of WEC in Different Sub-Basins
The most polluted dry season was selected as the study period. The quantity of reduced WEC per unit area of contaminants in the 24 sub-basins was calculated (Table 2) and the descriptive statistics were described as per Table 3. The spatial pattern of hot spots and safety areas within the basin was shown in Figure 11 (Hot spots represent the areas with the top three WEC reductions; Cutting area represent the area needed to be reduced; Safety area represent the area with no need to cut). The distribution pattern of hot spots for COD and TP were exactly the same (S19, S20 and S21) and were the industrial and populated clusters, while they were mainly distributed in populated areas (S24 and S25) for NH 3 -N. These were the key regions of environmental management. Nearly all safety areas were located in upstream watershed of the reservoir where the water source conservation area is located. The average reduction of 24 sub-basins for the COD, NH 3 -N, and TP pollutants was 24.07 t/a, 1.12 t/a, and 0.73 t/a, respectively. Among them, the largest reduction for COD was in S19 with up to 132.95 t/km 2 ·a −1 ; however, there were five sub-basins with a minimum of 0.00 t/km 2 ·a −1 , including S01, S05, S10, S11, and S13, indicating that there was still room to admit pollution loads in these regions. The most and least reduced regions for NH 3 -N were S21 (9.69 t/km 2 ·a −1 ) and elements of S03, S09, S11, and S13 (0.00 t/km 2 ·a −1 ), presenting a wide range of differences among the elements. Furthermore, S19 and S21 were the areas with the most TP reductions (2.07 t/km 2 ·a −1 ) and the areas of S11, S13, and S18 still had a remaining environmental capacity (0.00 t/km 2 ·a −1 ). As for statistical parameters of variation, the standard deviation of the reduced WEC for COD, NH 3 -N, and TP in the 24 sub-basins were 37.07 t/km 2 ·a −1 , 2.28 t/km 2 ·a −1 and 0.61 t/km 2 ·a −1 , respectively, and the variation coefficient of the contaminants reached 154.01%, 203.57%, and 141.86%, respectively, the variability of which belonged to high levels [40]. Thus, the reduced WEC dramatically exhibited a spatial difference between all divisions of the basin during the dry season. Sub-Basin S13 S14 S16 S17 S18 S19 S20

Conclusions
In 2015, COD, NH 3 -N, and TP discharge from diversified emission sources into the Xian-Jiang Basin were 6719.68 t/a, 488.12 t/a, and 128.57 t/a, respectively. Domestic sewage was identified as the primary source of all the contaminants, whereas industrial enterprises above the designated size were the least polluting sources for COD and NH 3 -N, with formalization cultivation for TP. According to the survey data, abundant domestic wastewater, especially for the villages, was not collected and accessed by the wastewater treatment plant through sewage conduits and was replaced by discharging into the river directly. Therefore, the increased acceptance rate and enhanced treatment capacity of domestic sewage has become very important and essential.
The modeling results indicated that the pollutant concentrations of the water system in the dry season were demonstrably worse than in the plum rain and typhoon seasons. The water qualities upstream of the reservoir, as a whole, were generally better than in the lower reaches where the concentrations basically met the Grade II standards. Furthermore, the middle part of the lower reaches was the most polluted waterbody as it was surrounded by industrial and residential cluster areas. The distribution of water quality standard exhibited a compelling correlation with the location of industrial enterprises, residential areas, as well as land-use types.
The WEC that needed to be eliminated for the basin to meet its water quality protective goal had a significant seasonal difference: the dry season the plum season > the typhoon season. Furthermore, there was also an obvious spatial difference in the reduced WEC between all sub-basins of the watershed. This finding implied that the current one-size-fits-all WEC cutting strategy adopted by environmental protection departments should be changed. Instead, WEC reduction should take spatio-temporal variation into account in the basin. If the reduced WEC during dry season can satisfy the required surface water quality standards, the probability of pollutant accidents is minimized annually where reducing the heavily polluted areas preferentially can help policy-makers save significant costs.
Although the finer the calculated elements, the higher the accuracy of water quality and WEC calculation becomes, this also increases the computation difficulty. Therefore, a balance point between the calculated elements scale and computation load requires further investigation. In addition, the model still has some deficiencies and room for improvement, such as different pollutant degradation coefficients that can be used considering the spatial heterogeneity of reaches, and increasing the monitoring sites for better calibration and verification. The water quality and environmental capacity simulation method discussed in this study can provide a good reference for other basins, especially for mountain river systems where long-term monitoring data are insufficient and confidential in China.