Abstract

Excessive load of nitrogen from anthropogenic sources is a threat to sustaining a healthy aquatic ecosystem. The difficulty in identifying the critical source areas (CSAs) of nitrogen load and apportioning the in-stream nitrogen to individual sources spatially and seasonally has made the Soil and Water Assessment Tool (SWAT) useful for analyzing nitrogen load at the catchment scale. However, the uncertainty of the nitrogen load simulated by SWAT has rarely been analyzed. The two simulations with the highest or the lowest PBIAS of total nitrogen (TN) load were proposed in this study to represent the range of the prediction uncertainty and therefore were used to generate the uncertainty of CSAs and nitrogen source apportionment. The model was set up for the Yuan River Catchment, which is under threat of extensive nitrogen load. Results indicated the highest nitrogen load was from downstream paddy fields with a denser population and 85% of the load was from fertilizer and feedlots. The relatively high prediction uncertainty was observed on both CSAs and source apportionment, which emphasizes the limitation of calibration only based on certain processes and the necessity to consider parameter uncertainty in the application of nitrogen load simulation.

INTRODUCTION

Excess nitrogen in a surface or groundwater system is a threat to sustaining a healthy aquatic ecosystem and thus a potential risk to human health. The sources of nitrogen are relatively complex and differ between each catchment. They can be categorized as point sources (PS) and non-point sources (NPS). Nitrogen from PS loads directly to streams, and is mainly regulated, e.g., effluent from industrial and domestic sewers. However, nitrogen from NPS, e.g., fertilizer operation, is released from varied sources and spatially distributed over the entire catchment. Among the sources, nitrogen load from agricultural activities is globally the most critical and the demand for fertilizer nitrogen has been estimated to have an annual growth rate of 1.5% globally to 2020 (FAO 2017). However, in China, waste from scattered feedlots and rural households have to be considered as two important additional nitrogen sources which are often disposed of with little or no treatment, contributing to the excessive amount of nitrogen in many rivers (Ongley et al. 2010).

The spatial distribution of the nitrogen sources and their contribution as nitrogen load to the stream are normally not homogeneous across catchments, especially in catchments of varied topography, weather, agricultural operations, and management strategies. However, in general, only a few regions in a catchment produce very high nitrogen load into streams (Liu et al. 2016). Case studies have shown that areas representing less than 10% of a catchment contributed to more than 50% of the total nitrogen (TN) load in streams (Niraula et al. 2011). The areas contributing to a higher load of nitrogen in a catchment are referred to as critical source areas (CSAs). Moreover, the complex nitrogen transformation processes make it challenging but also necessary to analyze the in-stream nitrogen contribution from each of the individual nitrogen sources. CSAs and nitrogen source apportionment are more targeted approaches to assessing the nitrogen load (White et al. 2009) and are regarded as prerequisites of developing a cost-effective strategy for catchment management (Lindgren et al. 2007).

The identification of the CSAs and the source apportionment of in-stream nitrogen require a comprehensive analysis of nitrogen from sources and the transformation and transportation processes at both the terrestrial phase and water system. Approaches developed for this purpose can be categorized as an empirical/statistical model or a process-based model according to their complexity of representing the hydrological and nitrogen cycles (Bouraoui & Grizzetti 2014). Process-based models apply physically based equations and derive the computational unit considering the topography, geology and climate condition, or even integrating the management operations, which are capable of tackling the complex processes of the nitrogen cycle and of generating a simulation with a higher spatial and temporal resolution (Arabi 2010). The Soil and Water Assessment Tool (SWAT) (Arnold et al. 1998), as one of the commonly applied process-based hydrological models, is selected in this study due to its capability of predicting the nitrogen in catchments in varied weather and geographic conditions. (Xu et al. 2016; Malagó et al. 2017; Tran et al. 2017).

Although the reliability of the SWAT model has been proved in many studies, the model's uncertainty has rarely been discussed in terms of nitrogen load prediction, which, however, is of more practical significance for CSAs and nitrogen source apportionment. The uncertainties of the nitrogen simulation are the combined effect of the model input data, model algorithm, and parameters. The algorithms in the SWAT model to simulate the nitrogen transformation and transportation processes are empirical relations determined by the self-defined coefficients (Neitsch et al. 2011), e.g., nitrate percolation coefficient, partitioning the nitrate between surface flow and percolation. Those parameters are mainly determined by the calibration of the in-stream nitrogen concentration, which displays a higher uncertainty in the measured data than in the discharge (Harmel et al. 2006). Those uncertainties along with the parameter uncertainty will amplify the uncertainty in the TN prediction.

The modeling procedure of calibration and uncertainty analysis has been outlined by Arnold et al. (2012) and Abbaspour et al. (2018). SUFI-2 (Sequential Uncertainty Fitting, ver. 2) (Abbaspour 2005), as an auto-calibration approach, derives the model uncertainty by generating the 95% prediction uncertainty (95PPU) for the calibration variables (e.g., discharge and nitrogen load time series). However, it is not practical also to generate the 95PPU for CSAs and source apportionment, due to the foreseeable high cost of computation. During the calibration, the simulations which have the highest or the lowest percent bias (PBIAS) (Gupta et al. 1999) can also be representatives of the upper and lower boundary of the 95PPU. Against this background, the parameter sets, resulting in the highest/lowest PBIAS, are used to generate the uncertainty of CSAs and the nitrogen source apportionment.

The Yuan River Catchment (YRC) in Jiangxi Province, China, has been chosen as the study area due to the excessive nitrogen load in the river, especially during the dry season (Fang 2011). Only a limited number of studies have analyzed the total pollutant load in the case study area with a focus on discharge from PS (Li & Liu 2015). Wu (2012), however, has shown already that apart from the nitrogen fertilizer application, untreated waste from rural households and from feedlots can be regarded as dominant NPS in Jiangxi Province, so that these nitrogen sources should not be neglected when analyzing the total nitrogen load in the stream. In 2015, the ‘Action Plan for Prevention and Control of Water Pollution’ (Ministry of Ecology and Environment 2015) was published to urge the local government to tackle the NPS in rural regions. The analysis of the presented CSAs and the apportionment of the nitrogen sources in the YRC may support this initiative.

With the above triggers, the study was performed with the following objectives: (1) establish a reliable hydrological and nitrogen model in the YRC; (2) identify the CSAs and apportion the in-stream nitrogen to individual sources; (3) analyze the model prediction uncertainty for CSAs and nitrogen source apportionment.

METHOD

SWAT and the calibration

The SWAT is an integrated process-based hydrological model to predict the long-term impact of land management practices on the water, sediment and agricultural chemical yields in catchments (Arnold et al. 1998). The computation unit in SWAT is the hydrologic response unit (HRU), which consists of a certain land-use type, soil type and slope in a sub-watershed. The hydrological cycle is simulated based on the water balance equation, with the SCS curve number method (USDA 1972) to quantify the daily surface runoff. The nitrogen transport model in SWAT calculates the nitrogen transformation processes in the soil with complex empirical equations. The sources of nitrogen, including fertilizer application and plant residue, are integrated with a plant growth module determined by the heat units approach (Barnard 1948) that the nitrogen consumption was determined by the growth stage of the plant. The in-stream nitrogen transformation process is assumed to be one-dimensional and computed by QUALE2 (Brown & Barnwell 1987).

The parameters for the SWAT model calibration were selected based on both global sensitivity analysis and the observation of the offset between the measured and simulated data. The Latin hypercube sampling technique was applied to generate the random parameter sets for simulations in calibration. The simulations were then evaluated to obtain suitable parameter sets. The objective functions applied to evaluate the simulations were Nash–Sutcliffe efficiency (NSE) (Nash & Sutcliffe 1970) and percent bias (PBIAS) (Gupta et al. 1999). The criteria to assess the model performance were to achieve a monthly flow simulation of NSE > 0.5, PBIAS < ±15 and a monthly nitrogen simulation of NSE > 0.35, PBIAS < ±30 (Moriasi et al. 2015). The algorithms for NSE and PBIAS are listed below:  
formula
 
formula

is the observed data at time step i, is the simulated data at time step i, and is the mean value over the time period n.

The model was calibrated in two steps. It was first calibrated by the monthly streamflow discharge at the three hydrometric stations from 2008 to 2010, and was then validated at the Jiangkou station from 2011 to 2014. The flow-related parameters are much less sensitive to the only nitrogen-process-related parameters (Arnold et al. 2012). Therefore, the nitrogen module was calibrated in the second step when the flow-related parameters had been updated with the calibrated values. Total nitrogen (TN) load was calibrated at the most downstream station (6_L in Figure 1) from 2008 to 2010 and validated from 2011 to 2014.

Figure 1

The land-use map of the Yuan River Catchment with dominant land-use classes listed. Luxi, Maozhou, and Jiangkou are the three hydrometric stations, 1_Y − 6_L are the six water quality monitoring stations, 7_X and 8_O are the outlets without monitoring stations.

Figure 1

The land-use map of the Yuan River Catchment with dominant land-use classes listed. Luxi, Maozhou, and Jiangkou are the three hydrometric stations, 1_Y − 6_L are the six water quality monitoring stations, 7_X and 8_O are the outlets without monitoring stations.

The analysis of nitrogen load

CSA identification

In order to identify the CSAs, the total nitrogen (TN) contributed to the stream from each HRU was derived. In order to generate a more detailed spatial distribution of TN load, the HRUs were classified with 0% of land-use and soil type eliminated. The TN load at each HRU was represented by the total annual load (kg) to the stream, which included the nitrate transported with the surface runoff, lateral flow, and groundwater return flow as a solute and the organic nitrogen absorbed to the sediment and transported. Only the TN from NPS was considered to identify the CSAs. The TN load of each HRU was later sorted in descending order and the cumulative TN was calculated so as to classify the TN load classes: class I as the total TN contribution of the first 25%; class II as the total TN contribution of the second 25% (class I and II contribute 50% of the TN in total); class III as the total TN contribution of the third 25%; and class VI as the total TN contribution of the remaining 25%. The lower the area of the region that covered the individual class, the more critical the region was defined.

Nitrogen source apportionment

A base scenario without nitrogen input from fertilizer (fert), feedlots, rural households (rural), and sewage discharge (sewer) was created. The nitrogen released to the stream in the base scenario, though not contributed by the dominant nitrogen sources, may come from, e.g., soil erosion and plant residue. Another four scenarios, each with only one of the four nitrogen sources as input, were generated and the average annual TN was calculated. The TN from those four scenarios had then the TN from the base scenario deducted from them. Therefore, the TN from each source with their own contribution to the total nitrogen in the stream can be obtained.

Uncertainty analysis

SUFI-2 (Sequential Uncertainty Fitting, ver. 2) (Abbaspour 2005) was applied to derive the uncertainty of nitrogen simulation. SUFI-2 generates the 95% prediction uncertainty (95PPU) between the 2.5th and 97.5th percentiles after all the iterations. A p-factor representing the coverage of the uncertainty band is applied and suggested to be larger than 0.7, indicating that at least 70% of the observed data are covered. The ratio of the average thickness of 95PPU and the standard deviation of the observed data is referred to as the r-factor, with a suggested value for flow calibration to be lower than 1.5. The range of parameters that generates the 95PPU with a suitable p-factor and r-factor is referred to as the parameter uncertainty. The parameter uncertainty includes uncertainty from input data, model algorithms, and observed data. In this study, a series of simulations were selected after nitrogen calibration with a p-factor of 0.7 and the respective parameters were extracted. The parameter uncertainty was then derived for each parameter by taking the highest and lowest values applied in those simulations.

The model prediction uncertainty was also analyzed for CSAs and nitrogen source apportionment. However, in order to simplify the process, not all the parameter sets during the calibration iteration were applied to derive the CSAs and apportion the nitrogen sources. Only two of the parameter sets were selected and they were from the iterations with the lowest and highest PBIAS to the observed data. Due to the fact that PBIAS is able to quantify either underestimation (low PBIAS) or overestimation (high PBIAS), it can derive the simulations with the largest distance (both positive and negative) to the observed one. Therefore, these two simulations were used to represent the prediction uncertainty.

STUDY AREA

Catchment overview

The YRC, displayed in Figure 1, is located in the midwest of Jiangxi Province, China (113°50′–115°32′E; 27°22′–28°08′N). The Yuan River (YR) originates from the southwest mountainous region, at an elevation of approximately 1,900 m, and joins the Gan River on the northeast plain. The catchment covers an area of 6,223.69 km2. The YRC is dominated by a subtropical monsoon climate, with long-term annual precipitation of 1,668 mm and a pan evaporation rate of 1,228 mm (Fang 2011). The rainy season is from April to June and contributes 43.6% of the total precipitation. The mean temperature ranges from 2.83 °C in January to 34.2 °C in July. Red soil and paddy soil dominate 50.3% and 30.7% of the YRC, respectively.

The YR had a population in 2010 of 2.4 million with an average growth rate of 7.5‰ (Jiangxi Statistical Bureau 2008–2014). The population is 54.6% rural residents and is scattered across the catchment, occupying 4.1%, including big towns and small villages, as displayed in Figure 1. The cropland in the YRC is dominated by paddy field with a rice plantation of two seasons per year and covers 29.8% of the catchment. Shrub and forest cover 54.8% and the rest includes dry arable land (4.9%) with vegetables and orchards (0.3%) with mandarins. Livestock and poultry feedlots are distributed mainly in the rural region, with an average annual pig production of 2.27 million. The YR has been under the threat of excess nitrogen in recent years (Li et al. 2011), which is not unique in a southern Chinese catchment, with dense population and agricultural activities. The nitrogen fertilizer, unregulated waste from feedlots and rural households, along with the sewage from industry and wastewater treatment plants in the urban region are the dominant nitrogen sources in the YRC.

Input data description

The input data for the SWAT model included the digital elevation model (DEM), land-use and cover, soil properties and climate data. The data used in the YRC model are described below in Table 1. The daily precipitation was obtained from 12 stations from 2008 to 2014 at the upstream of Jiangkou Hydrometric Station. The precipitation data downstream were obtained from the China Meteorological Assimilation Driving Datasets for the SWAT Model (CMADS) (Meng & Wang 2017) as well as the daily climate data, including temperature, wind speed, solar radiation and relative humidity across the entire YRC. The total population and animal production in the YRC were obtained at the scale of administrative units, and within an administrative unit, the summed data were recorded for cities, towns, and villages individually. Therefore, the data were distributed to each village by the weight of the area. The villages were extracted from the land-use data, shown as rural residential in Figure 1, and then aggregated at each sub-catchment. The waste from rural households and feedlots was disposed as manure (rural household: 57%; feedlot: 90%) or discharged directly to drainage (rural household: 43%; feedlot: 10%). The untreated wastewater from the rural households was discharged with 71% directly to the waterbody and 29% to the agricultural land. The nitrogen fertilizer rate on rice was 200–240 kg/ha in both spring (Apr. to Jul.) and autumn (Jul. to Oct.) seasons, including base fertilizer and three times as top dressing.

Table 1

SWAT model input data of the YRC

CategoryDataResolution/scaleSources
Hydrology DEM 30 m Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Science 
Yuan River network 1:200,000 
Land-use and in 2010 30 m 
Soil type and properties 1,000 m Institute of Soil Science, Chinese Academy of Science 
Precipitationa 12 stations Jiangkou Hydropower Station 
Climate data (2008–2014)b 1/3 degree, daily China Meteorological Assimilation Driving Datasets for the SWAT Model (2008–2014), Chinese Academy of Science 
Reservoir discharge (2008–2014) Daily Jiangkou Hydropower Station 
Industry and domestic water use Annual average (2008–2014) Statistical Year Book of Jiangxi Prov., 2008–2014 
Discharge (2008–2010; 2008–2014) 2 stations, daily Yichun Hydrologic Bureau, JK reservoir power station 
Nitrogen Crop and fertilization Fertilizer application per hectare Local survey and literature 
Animal production Annual average (2008–2014) Statistical Year Book of Jiangxi Prov., 2008–2014 
Rural population Annual average (2008–2014) Statistical Year Book of Jiangxi Prov., 2008–2014 
The effluent of sewage from industry and wastewater treatment plant The urban region, seasonally (2008–2014) Yichun and Xinyu Environmental Protection Bureau, publicly available 
Monitored water quality of NH4 and TN 6 stations, monthly or seasonally (2008–2014) Xinyu Environmental Protection Bureau 
CategoryDataResolution/scaleSources
Hydrology DEM 30 m Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Science 
Yuan River network 1:200,000 
Land-use and in 2010 30 m 
Soil type and properties 1,000 m Institute of Soil Science, Chinese Academy of Science 
Precipitationa 12 stations Jiangkou Hydropower Station 
Climate data (2008–2014)b 1/3 degree, daily China Meteorological Assimilation Driving Datasets for the SWAT Model (2008–2014), Chinese Academy of Science 
Reservoir discharge (2008–2014) Daily Jiangkou Hydropower Station 
Industry and domestic water use Annual average (2008–2014) Statistical Year Book of Jiangxi Prov., 2008–2014 
Discharge (2008–2010; 2008–2014) 2 stations, daily Yichun Hydrologic Bureau, JK reservoir power station 
Nitrogen Crop and fertilization Fertilizer application per hectare Local survey and literature 
Animal production Annual average (2008–2014) Statistical Year Book of Jiangxi Prov., 2008–2014 
Rural population Annual average (2008–2014) Statistical Year Book of Jiangxi Prov., 2008–2014 
The effluent of sewage from industry and wastewater treatment plant The urban region, seasonally (2008–2014) Yichun and Xinyu Environmental Protection Bureau, publicly available 
Monitored water quality of NH4 and TN 6 stations, monthly or seasonally (2008–2014) Xinyu Environmental Protection Bureau 

aPrecipitation data from Jiangkou Hydropower Station only include region at the upstream of the Jiangkou hydrometric station.

bPrecipitation data at the downstream of the YRC are from East Asia Meteorological Assimilation dataset.

RESULTS AND DISCUSSION

Model performance

Discharge calibration

Table 2 lists the parameters involved in the flow calibration. The decrease of the CN2 value was to derive a lower peak discharge via calibration. The decrease of ESCO, REVAPMN and the increase of GW_REVAP tended to reduce the base flow rate, thus lowering the simulated discharge. Figure 2 displays the simulated streamflow discharge (in blue) and Table 3 displays the values of the objective functions at each of the three hydrometric sites after the calibration. The seasonal trend is well captured by the calibrated parameter set and the objective function values indicate satisfactory to very good model performance. However, an underestimation at low flow periods can still be observed and with a higher discrepancy than at peak flow periods. This can be explained by the characteristic of the objective function NSE, which weighs the peak values more than the low values, or the lack of baseflow information to constrain the low flow period. Moreover, the discrepancy was also a result of the uncaptured anthropogenic recharge or discharge. In the YRC, there are no hydrometric stations at the downstream of the Jiangkou Hydrometric Station, and thus the downstream cannot be calibrated by the discharge. However, it is assumed that the downstream shares similarities of identical land-use types (e.g., paddy field, forest) and soil types (e.g., red soil, paddy soil) with the other regions in the YRC, thus the relevant parameters are applicable at the downstream (e.g., CN2, SOL_AWC). Therefore, the parameters listed in Table 2 were also applied at the downstream.

Table 2

Parameters calibrated for streamflow discharge in the SWAT model for YRC

ParameterDescriptionDefaultCalibrated
CN2 SCS runoff curve number of the forest for moisture condition II a Decreased by 0.74%b 
SOL_AWC Available water capacity of the soil layer a Decreased by 10.29%b 
ESCO Soil evaporation compensation factor 0.90 0.31 
GWQMN Threshold depth of water in the shallow aquifer required for return flow to occur (mm) 1,000 708 
GW_REVAP Groundwater revap coefficient 0.02 0.19 
REVAPMN Threshold depth of water in the shallow aquifer for revaporation to occur (mm) 750 36 
ParameterDescriptionDefaultCalibrated
CN2 SCS runoff curve number of the forest for moisture condition II a Decreased by 0.74%b 
SOL_AWC Available water capacity of the soil layer a Decreased by 10.29%b 
ESCO Soil evaporation compensation factor 0.90 0.31 
GWQMN Threshold depth of water in the shallow aquifer required for return flow to occur (mm) 1,000 708 
GW_REVAP Groundwater revap coefficient 0.02 0.19 
REVAPMN Threshold depth of water in the shallow aquifer for revaporation to occur (mm) 750 36 

aParameter is distributed spatially over the YRC.

bParameter is calibrated by an increase/decrease of a certain percentage over the YRC.

Table 3

The values of objective functions NSE and PBIAS after calibration and validation at the three hydrometric sites

Objective functionLuxiMaozhouJiangkouJiangkou
 Calibration (2008 − 2010) Validation (2011–2014) 
NSE 0.72 0.81 0.92 0.90 
PBIAS −19.9 8.8 1.7 −9.1 
Objective functionLuxiMaozhouJiangkouJiangkou
 Calibration (2008 − 2010) Validation (2011–2014) 
NSE 0.72 0.81 0.92 0.90 
PBIAS −19.9 8.8 1.7 −9.1 
Figure 2

Monthly streamflow at the Jiangkou site in both the calibration (2008–2010) and the validation (2011–2014) periods, with observed streamflow in red and simulated streamflow in blue.

Figure 2

Monthly streamflow at the Jiangkou site in both the calibration (2008–2010) and the validation (2011–2014) periods, with observed streamflow in red and simulated streamflow in blue.

Total nitrogen (TN) load calibration

The nitrogen module was calibrated after updating the flow-related parameters in the model. The most sensitive parameters are listed in Table 4 and they control the different processes of nitrogen transformation or transportation. Unlike the streamflow discharge, the total nitrogen (TN) consists of nitrate (NO3-N), organic nitrogen and ammonia nitrogen (NH4-N). In this study, only measured TN was available instead of the individual nitrogen species. The in-soil process was calibrated by CDN, SDNCO, LAT_TIME, RSDCO, ERORGN, and NPERCO, so as to increase the rate of denitrification and organic nitrogen decomposition processes and to prolong the time of the nitrate travel time with the lateral flow. The in-stream process was calibrated by K_N, RS4, and BC1 so as to increase the algae to uptake nitrate and to decrease the rate of the nitrification and organic nitrogen settling processes. The simulated and observed TN loads are displayed in Figure 3. The values of NSE (calibration: 0.57; validation: 0.73) and PBIAS (calibration: −3.7; validation: 8.7) indicate an acceptable simulation of the SWAT model with, however, larger discrepancies than the flow calibration. The uncertainties of the nitrogen input dataset, e.g., annual average input of nitrogen from feedlots and rural households, may explain the relatively low reliability of the simulation in some months. Moreover, the underestimation of the peak load could also be caused by the underestimation of the peak flow (see Figure 2). In 2008, the discharge is well presented at Jiangkou station (see Figure 2) but the respective modeled and observed TN loads (see Figure 3) are not satisfactory. The reason may be the documented data gaps in the statistics of nitrogen sources in 2008 or the not fully calibrated discharge at the downstream. Despite the discrepancies, the SWAT model generated an acceptable TN load such that the NSE at other sites after calibration (from 1_Y to 5_S) ranges from 0.49 to 0.73.

Table 4

Parameters calibrated and parameter uncertainty range (min and max) for TN simulation in the SWAT model for the YRC

ParameterDescriptionDefaultCalibratedMinMax
CDN Denitrification rate coefficient 1.40 1.91 0.27 2.99 
K_N Michaelis–Menten half-saturation constant for nitrogen 0.02 0.24 0.01 0.30 
RS4 Rate coefficient for organic N settling in the reach at 20 °C 0.05 0.09 0.10 
SDNCO Denitrification threshold water content 1.10 1.85 0.01 2.94 
LAT_TTIME Lateral flow travel time [/day] 129 180 
RSDCO Residue decomposition coefficient 0.05 0.09 0.05 0.10 
BC1 Rate constant for biological oxidation of NH4 to NO2 in the reach at 20 °C 0.55 0.10 0.50 0.10 
ERORGN Organic N enrichment ratio 3.91 0.02 4.99 
NPERCO Nitrogen percolation coefficient 0.200 0.969 0.003 0.997 
ParameterDescriptionDefaultCalibratedMinMax
CDN Denitrification rate coefficient 1.40 1.91 0.27 2.99 
K_N Michaelis–Menten half-saturation constant for nitrogen 0.02 0.24 0.01 0.30 
RS4 Rate coefficient for organic N settling in the reach at 20 °C 0.05 0.09 0.10 
SDNCO Denitrification threshold water content 1.10 1.85 0.01 2.94 
LAT_TTIME Lateral flow travel time [/day] 129 180 
RSDCO Residue decomposition coefficient 0.05 0.09 0.05 0.10 
BC1 Rate constant for biological oxidation of NH4 to NO2 in the reach at 20 °C 0.55 0.10 0.50 0.10 
ERORGN Organic N enrichment ratio 3.91 0.02 4.99 
NPERCO Nitrogen percolation coefficient 0.200 0.969 0.003 0.997 
Figure 3

The simulated (Sim_best) and observed total nitrogen (TN) load at station 6_L (see Figure 1), with the 95PPU in the shaded area and simulations of the highest and lowest PBIAS. The calibration period is from 2008 to 2011 and the validation period is from 2011 to 2014.

Figure 3

The simulated (Sim_best) and observed total nitrogen (TN) load at station 6_L (see Figure 1), with the 95PPU in the shaded area and simulations of the highest and lowest PBIAS. The calibration period is from 2008 to 2011 and the validation period is from 2011 to 2014.

Parameter uncertainty and prediction uncertainty of TN simulation

The ranges of the parameter values are listed in Table 4, representing the parameter uncertainty of the TN simulation. The parameter uncertainty consists of not only the parameter uncertainty but also the uncertainty of the input nitrogen data, model algorithm, and the observed data. The prediction uncertainty is displayed in Figure 3 as the 95PPU illustrates 95% of all the simulation results obtained due to the parameter uncertainty. The 95PPU covers 64% (p-factor) of the observed TN and the r-factor is 1.29. Unlike the best simulation displayed in Figure 3, which demonstrates a good fit to the observed TN, the prediction uncertainty indicates a relatively larger deviation from the observed TN, especially at the high or low TN load. The iterations with the highest and lowest PBIAS are also displayed in Figure 3. They are observed to cover most of the model prediction uncertainty.

The nitrogen transformation and transportation process is rather complex, in terms of the number of processes and the number of parameters involved, as well as their correlation. In the SWAT model, nitrogen-related processes, e.g., denitrification and organic nitrogen transported by sediment, are determined by a rate coefficient. It results in a heavy dependence on calibration in order to derive a convincing nitrogen simulation, which has also been pointed out by Zeiger & Hubbart (2016). Moreover, in our case, the nitrogen species measured was TN, which increased the uncertainty of the model system in the calibration process, as the portions of nitrate and organic nitrogen were unknown and were adjusted by calibration. Therefore, the lack of constraints on the individual nitrogen relevant processes also leads to a relatively high prediction uncertainty of TN load, as well as a less reliable simulation of the internal processes.

The identification of CSAs and the uncertainty

The 4,137 HRUs were sorted in descending order of annual TN load (kg) and were classified as class I, II, III and IV to denote the contribution of the total TN load of the first 25%, second 25%, third 25% and the remaining 25%, as displayed in Figure 4. The cumulative plot of TN and the respective area was created, as displayed in Figure 5. The classes I, II and III are displayed to show the respective percentage of the area occupied by each class. It can be observed in Figure 4 that several sub-catchments, e.g., at the upstream or at the downstream, dominate the contribution of the TN. Figure 5 indicates that the total area occupied by class I, II and III is less than 50%, however, it contributes 75% of the TN in the YRC. Class I occupies approximately 10% of the area and class II occupies approximately 15%.

Figure 4

The classification of total nitrogen (TN) at the Yuan River Catchment.

Figure 4

The classification of total nitrogen (TN) at the Yuan River Catchment.

Figure 5

Accumulative TN from each HRU of calibration results (Sim_best) and simulations with the highest and lowest PBIAS, with the three dotted horizontal lines demonstrating the accumulative TN contributed by class I, II and III in Sim_best.

Figure 5

Accumulative TN from each HRU of calibration results (Sim_best) and simulations with the highest and lowest PBIAS, with the three dotted horizontal lines demonstrating the accumulative TN contributed by class I, II and III in Sim_best.

TN load within class I is mainly distributed at the downstream plain region, which is occupied by paddy fields and a denser rural population, indicating a high TN load. When compared with other regions with paddy fields, its relatively high TN load could arise from higher leaching rates and higher runoff rates. Class II and class III are dominated by forest in the mountainous region and paddy fields at the floodplain. The deeper slope and the forest residue could account for the high TN load in the forest. In the YRC, the downstream plain region, with denser paddy fields and rural population, should be considered in particular as the critical region of TN load. Meanwhile, the mountainous region with forest requires further analysis to validate the high TN load.

The prediction uncertainty analysis of CSAs was based on the two simulations, one with the lowest and the other with the highest PBIAS. As stated above, these two simulations, with the largest deviation from the observation in all the simulations, could represent the range of the prediction uncertainty. Figure 6(a) and 6(b) indicate an obvious variation of the CSAs identified at the HRU level. Each of the classes was affected by the parameter uncertainty to a different extent and the distribution of class I was relatively less impacted. When compared with Figure 4, Figure 6(a) illustrates a more averaged occupation by all classes and Figure 6(b) with less area contributing 75% of the TN. The results are more obvious in Figure 5. The cumulative curve of the simulation with the lowest PBIAS can hardly be identified with critical areas as the slope of the curve is approximately 1 without a significant change, whereas the simulation with the highest PBIAS has more concentrated critical areas so that 25% of the area contributes almost 80% of the TN.

Figure 6

CSAs identified from the simulation (a) with the lowest PBIAS, (b) with the highest PBIAS.

Figure 6

CSAs identified from the simulation (a) with the lowest PBIAS, (b) with the highest PBIAS.

The prediction uncertainty due to parameter uncertainty reveals an uneven impact over the catchment to identify the CSAs. Other studies (Niraula et al. 2011; Liu et al. 2016), however, have concluded only little impact on the locations of CSAs before and after calibration. The disagreement in this research may be explained by the simulations in Figure 6 deriving a relatively larger deviation from the best simulation than the uncalibrated simulations from previous studies. Moreover, the parameter sets were not identical in those studies and might result in the adjustment of different processes in the model. In this study, when the range of the prediction uncertainty in TN load is rather large, it also impacts significantly the CSA identification. Therefore, the calibration of the nitrogen model is critical, if a convincing CSA identification is to be derived.

The nitrogen source apportionment and the uncertainty

Figure 7 demonstrates the annual TN load along the YR, from the upstream to the outlet; 1–8 stands for eight outlets, where 1–3 are located at the up-midstream, 4–5 are located at the midstream and the rest are at the downstream, and 8 is the outlet of the YRC. It can be observed at the YRC outlet that fertilizer contributes the highest amount of TN in the stream and TN from feedlots follows. These two sources contribute more than 85% of the TN at the outlet, as displayed in Table 5. The rural households (rural) contribute 8.8% of the TN and discharge from sewage (sewer) only contributes 6.1% of the TN. In Figure 8, the TN at the outlet of the YRC is displayed and separated by the operational season (op), when agricultural activities are operated, and the non-operational season (other). The TN load from fertilizer is obviously contributing more from the operational season. The other scenarios also show a higher release from the operational season.

Table 5

The annual TN load from the individual scenarios at the outlet of the YRC

SourcesNitrogen fertilizerFeedlotRural householdSewage
Scenario fert feedlot rural sewer 
Annual TN in-stream load [t/a] 3,323 2,679 620 428 
Percentage 47.1% 38.0% 8.8% 6.1% 
SourcesNitrogen fertilizerFeedlotRural householdSewage
Scenario fert feedlot rural sewer 
Annual TN in-stream load [t/a] 3,323 2,679 620 428 
Percentage 47.1% 38.0% 8.8% 6.1% 
Figure 7

The TN load after calibration and its uncertainty range from simulations with highest PBIAS and lowest PBIAS along the Yuan River at eight outlets (1–8), from upstream (1) to the outlet of the YRC (8); fert stands for nitrogen fertilizer, rural stands for rural households and sewer stands for urban sewage.

Figure 7

The TN load after calibration and its uncertainty range from simulations with highest PBIAS and lowest PBIAS along the Yuan River at eight outlets (1–8), from upstream (1) to the outlet of the YRC (8); fert stands for nitrogen fertilizer, rural stands for rural households and sewer stands for urban sewage.

Figure 8

The TN load at the YRC outlet (outlet 8 in Figure 1) in both the operational season (op) from Apr. to Oct. and the non-operational season (other) from Jan. to Mar.

Figure 8

The TN load at the YRC outlet (outlet 8 in Figure 1) in both the operational season (op) from Apr. to Oct. and the non-operational season (other) from Jan. to Mar.

The high load from nitrogen fertilizer could be estimated from the high amount of fertilizer application in the two rice plantation seasons. In Figure 7, the nitrogen fertilizer shows a sharper increase from outlet 1 to 2, from 4 to 5 and from 6 to 7. The sub-catchments inbetween those outlets are crop-dominated regions, and thus contribute a higher amount of TN. The distribution of the TN load can also be observed in Figure 4. As for the sewage discharge (sewer), a higher increase can be observed from outlet 2 to 3 and from 4 to 5, where the urban regions are located. Previous research (Zhang et al. 2011), accounting for nitrogen contribution in the dry season mainly from sewage, concluded there was a deteriorating water quality due to the sewage discharge. However, although a lower contribution from NPS could be observed in the non-operational (other) season, which is also the dry season in the YRC, the contribution of TN load in the dry season from groundwater return flow cannot be neglected and is also mainly from nitrogen fertilizer, as observed in Figure 8. In the YRC, in terms of analyzing the water quality due to nitrogen load, not only the sewage discharge but also the NPS should be considered, especially with the nitrogen fertilizer and from feedlots.

The prediction uncertainty of each scenario is also displayed in Figure 7. It is observed that the load from nitrogen fertilizer displays the largest range of uncertainty, followed by the waste from the feedlot. It could be explained that denitrification-related parameters adjusted in the calibration significantly affected the transformation of the nitrate, which was the dominant nitrogen species in the nitrogen fertilizer. Nitrogen from feedlots was dominated by organic nitrogen, but the narrower uncertainty range of simulated TN demonstrates less sensitivity of organic nitrogen transportation and transformation processes to the parameter uncertainty. The uncertainty range can barely be observed from rural households and sewage scenarios and it might be explained by their relatively low load when compared with the load from feedlots or fertilizer, and thus they were less affected. Moreover, the in-stream processes were also less sensitive to the parameter uncertainty than the in-soil processes. The best simulation in the calibration confirmed that the fertilizer contributes the highest TN to the river, however, when parameter uncertainty cannot be narrowed for the YRC, the interpretation of the results can be less convincing. Moreover, due to the non-linearity of the system, the nitrogen source apportionment approach also has its limitations. The nitrogen input of each scenario, which is regarded as a linear decrease of the input to the model, has a varied impact on the individual processes that can hardly be quantified and thus adds uncertainty to the model. The results also demonstrate the limitations of an integrated hydrological model like SWAT that if the calibration is only applied based on one variable, e.g., TN content in the stream, due to data scarcity, the uncertainty of the other relevant processes cannot be constrained.

CONCLUSION

The nitrogen load in terms of identifying CSAs and apportioning the in-stream nitrogen to individual sources was analyzed for the YRC in China based on the SWAT model. The prediction uncertainty of the nitrogen load was also analyzed, which was derived from simulations with the highest and the lowest PBIAS on total nitrogen (TN).

The key findings of the research are as follows:

  • (1)

    The applicability of the SWAT model simulating the total nitrogen load was analyzed for the YRC with multiple point or NPS. The results of CSAs identified that the downstream plain region with paddy fields and denser rural population released the highest TN load; nitrogen fertilizer and waste from feedlots contributed more than 85% of the TN, especially from April to October, at the YRC outlet. If future nitrogen management strategies are to be developed in the YRC, the areas and sources targeted above are recommended to be considered.

  • (2)

    The uncertainty derived for the CSAs and nitrogen source apportionment displayed a rather high impact on the prediction. The results provided a reference to demonstrate the relatively high uncertainty in the nitrogen load prediction using the SWAT model when multiple parameters are selected for calibration but are not able to be constrained by individual hydrology and nitrogen relevant processes due to limited data availability.

As excessive nitrogen load is been a concern of catchment management worldwide and nitrogen sources vary among catchments, the identification of CSAs and source apportionment based on integrated hydrological models would be supportive. However, the research showed a rather high prediction uncertainty of the nitrogen load simulation using the SWAT model under the current condition of the YRC, which is representative of catchments with rather limited available data and previous studies. Therefore, in the future, a more complete observation system for flow and nitrogen species would be of great benefit. In summary, although a good parameter set could be obtained via calibration to meet objective function criteria, the parameter uncertainty and the resulted prediction uncertainty should not be ignored in order to obtain a comprehensive nitrogen load analysis.

ACKNOWLEDGEMENTS

The first author is financially supported by the Chinese Scholarship Council (CSC).

REFERENCES

REFERENCES
Abbaspour
K. C.
2005
Calibration of hydrologic models: when is a model calibrated?
In:
MODSIM 2005 International Congress on Modelling and Simulation
(
Zerger
A.
,
Argent
R. M.
, eds),
Modelling and Simulation Society of Australia and New Zealand
,
Melbourne, Australia
, pp.
2449
2455
.
Arabi
M.
2010
Nutrient Modeling Overview
.
Colorado State University
,
Fort Collins, CO, USA
. .
Arnold
J. G.
,
Srinivasan
R.
,
Muttiah
R. S.
,
Williams
J. R.
1998
Large area hydrologic modeling and assessment part I: model development
.
Journal of the American Water Resources Association
34
(
1
),
73
89
.
Arnold
J. G.
,
Moriasi
D. N.
,
Gassman
P. W.
,
Abbaspour
K. C.
,
White
M. J.
,
Srinivasan
R.
,
Santhi
C.
,
Harmel
R. D.
,
van Griensven
A.
,
Van Liew
M. W.
,
Kannan
N.
,
Jha
M. K.
2012
SWAT: model use, calibration and validation
.
American Society of Agricultural and Biological Engineers
55
(
4
),
1491
1508
.
Barnard
J. D.
1948
Heat units as a measure of canning crop maturity
.
The Canner
106
,
28
.
Bouraoui
F.
,
Grizzetti
B.
2014
Modelling mitigation options to reduce diffuse nitrogen water pollution from agriculture
.
Science of The Total Environment
468–469
,
1267
1277
.
Brown
L. C.
,
Barnwell
T. O.
Jr.
1987
The Enhanced Water Quality Models QUAL2E and QUAL2E-UNCAS: Documentation and User Manual. EPA/600/3-87/007
.
USEPA
,
Athens, GA, USA
.
Fang
Z.
2011
The Report of Integrated Planning of Yuan River Catchment, Jiangxi Province
.
Jiangxi Water Conservancy Planning Design Institute
,
Nanchang, China (in Chinese)
.
FAO
2017
World Fertilizer Trends and Outlook to 2020: Summary Report
.
Food and Agriculture Organization of the United Nations
,
Rome, Italy
.
Gupta
H. V.
,
Sorooshian
S.
,
Yapo
P. O.
1999
Status of automatic calibration for hydrologic models: comparison with multilevel expert calibration
.
Journal of Hydrologic Engineering
4
(
2
),
135
143
.
Harmel
R. D.
,
Cooper
R. J.
,
Slade
R. M.
,
Haney
R. L.
,
Arnold
J. G.
2006
Cumulative uncertainty in measured streamflow and water quality data for small watersheds
.
Transactions of the ASABE
49
(
3
),
689
701
.
Jiangxi Statistical Bureau
2008–2014
Jiangxi Statistical Yearbook 2008, 2009, 2010, 2011, 2012, 2013, 2014
.
China Statistics Press
,
Beijing, China
. .
Li
Y.
,
Liu
C.
2015
The current situation of water quality in Xinyu section of Yuanhe River and its prevention
.
Journal of Xinyu University
20
(
1
),
13
15
(in Chinese)
.
Li
Z.
,
Yang
W.
,
Xu
F.
,
Ma
C.
2011
Distribution of pollutants in the sediments along the Yuan River
.
Journal of China Institute of Water Resources and Hydropower Research
9
(
2
),
132
136
(in Chinese)
.
Lindgren
G. A.
,
Wrede
S.
,
Seibert
J.
,
Wallin
M.
2007
Nitrogen source apportionment modeling and the effect of land-use class related runoff contributions
.
Hydrology Research
38
(
4–5
),
317
331
.
Malagó
A.
,
Bouraoui
F.
,
Vigiak
O.
,
Grizzetti
B.
,
Pastori
M.
2017
Modelling water and nutrient fluxes in the Danube River Basin with SWAT
.
Science of The Total Environment
603–604
,
196
218
.
Ministry of Ecology and Environment
2015
Action Plan for Prevention and Control of Water Pollution
.
People's Publishing House
,
Beijing, China
.
Moriasi
D. N.
,
Gitau
M. W.
,
Pai
N.
,
Daggupati
P.
2015
Hydrologic and water quality models: performance measures and evaluation criteria
.
Transactions of the ASABE
58
(
6
),
1763
1785
.
Neitsch
S. L.
,
Arnold
J. G.
,
Kiniry
J. R.
,
Williams
J. R.
2011
Soil and Water Assessment Tool, Theoretical Documentation: Version 2009
.
TR-406, Texas A&M University
,
College Station, TX, USA
.
Niraula
R.
,
Kalin
L.
,
Wang
R.
,
Srivastava
P.
2011
Determining nutrient and sediment critical source areas with SWAT: effect of lumped calibration
.
Transactions of the ASABE
55
(
1
),
137
147
.
Ongley
E. D.
,
Xiaolan
Z.
,
Tao
Y.
2010
Current status of agricultural and rural non-point source pollution assessment in China
.
Environmental Pollution
158
(
5
),
1159
1168
.
Tran
V.-B.
,
Ishidaira
H.
,
Nakamura
T.
,
Do
T.-N.
,
Nishida
K.
2017
Estimation of nitrogen load with multi-pollution sources using the SWAT model: a case study in the Cau River Basin in Northern Vietnam
.
Journal of Water and Environment Technology
15
(
3
),
106
119
.
USDA
1972
National Engineering Handbook
,
Hydrology Section 4, chapters 4–10
.
US Government Printing Office
,
Washington, DC, USA
.
White
M. J.
,
Storm
D. E.
,
Busteed
P. R.
,
Stoodley
S. H.
,
Phillips
S. J.
2009
Evaluating nonpoint source critical source area contributions at the watershed scale
.
Journal of Environment Quality
38
(
4
),
1654
1663
.
Wu
C.
2012
Research on Evaluation and Controlling Measure System of Agriculture Non-Point Source Pollution in Jiangxi Province
.
Master's thesis
,
China University of Geosciences (Beijing)
,
Beijing, China (in Chinese)
.
Zhang
X.
,
Li
H.
,
Xiao
Y.
,
Chen
Y.
2011
Study on allocation methods of the water pollutant loads in ungauged basin
.
Journal of China Institute of Water Resources and Hydropower Research
9
(
2
),
136
142
(in Chinese)
.