Assessing ways to combat eutrophication in a Chinese drinking water reservoir using SWAT

Additional keywords: agriculture, reservoir water quality, sewage, subtropical watershed.

For decades anthropogenic activities have deteriorated the water quality of freshwater lakes and reservoirs worldwide, with impacts on both biological diversity and the safe use of water resources (Søndergaard and Jeppesen 2007). Today, most western countries have implemented measures to target point-source related nutrient loadings with a special focus on sewage, and diffuse loading is therefore now of particular importance for the overall nutrient input to freshwater systems in this part of the world (Kronvang et al. 2005). For developing countries, however, point sources, especially sewage discharges associated with rural living, remain an excessive source of nutrients entering freshwater ecosystems and they continue to increase in parallel with the growing economy and improvement of rural life quality (Chen et al. 2008). In addition, increasing food demands in developing countries challenge the progress of local mitigation of diffuse nutrient loadings.

In China, the eutrophication of lakes and reservoirs is recognised as being linked to continuous and increasing losses of nutrients to freshwater systems from their watersheds (Jin et al. 2005). Within the country, lakes and reservoirs are important sources of drinking water and supply more than 90% of the total water withdrawal in certain areas (FAO 2011). Rapid development of the Chinese population and economy has led to a progressively deteriorating quality of surface waters and water supply (Gleick et al. 2009), which can largely be attributed to growing food demands and a consequent intensification of land use and management (expanding and intensifying agricultural practices, forestry and urbanisation) (Ju et al. 2004; Khan et al. 2009) as well as increased point source loading (Yan et al. 2011). Since the implementation of rural economic reforms in the late 1970s (Lin 1992), the concept of agricultural practices has changed in China from self-sustaining units with a high degree of nutrient recycling between livestock and crop to current practices with an increasing use of industrial fertilizers in the production (Chen et al. 2008). Consequently, the nutrient surplus at farm level has increased radically, with markedly enhanced net losses of nutrients, triggering environmental impacts especially on receiving surface waters (Ju et al. 2004). The current annual nitrogen (N) application rates for Chinese subtropical rice paddy areas have been reported to range from 264 to 805 kg N ha^–1, with an average of 556 kg N ha^–1 (Hatano et al. 2002), which by far exceeds the legal application in most developed countries and the quantity assumed necessary to achieve high yields (Ju et al. 2004). For example, according to Danish law (2011–2012 legislation), annual maximum allowable N application rates are 100, 172 and 121 kg N ha^–1 for oat, winter wheat and spring barley, respectively (MFE 2011).

Excessive loading of nutrients to drinking water reservoirs is of major concern for local water managers because of the risk of eutrophication and ensuing blooms of potentially toxin-producing cyanobacteria. The ability to account for the linkages between land use and water quality in reservoirs in a holistic and integrated manner, with the aim of achieving a sustainable and healthy water supply, is therefore of particular interest. Many concepts for integrated watershed management exist (e.g. Hu 1999; Huang and Xia 2001), but the key prerequisite for making decisions to improve water quality is an adequate estimation of the external nutrient loading, including all major external and internal processes that influence water quality.

In the present study, we took such a holistic approach and applied the river-basin-based ‘Soil and Water Assessment Tool’ model (SWAT) (Arnold et al. 1998) to the watershed of the subtropical Kaiping (Dashahe) reservoir located in southern China. The reservoir serves as drinking water supply for approx. 0.4 million people, but recent blooms of cyanobacteria are now challenging local water managers to take action to stop further deterioration of the water quality (Han et al. 2012). SWAT has been used in several studies in China to assess impacts of climate changes, land-cover changes and land use (e.g. Ouyang et al. 2008; Li et al. 2011) on hydrology and, to some extent, on nutrient export. However, to incorporate the connectivity between watersheds and reservoirs, we expanded the ordinary SWAT-watershed approach by coupling it with a simple, but widely used, empirical equation to estimate total phosphorus (TP) concentrations in the reservoir (Vollenweider and Kerekes 1982) for use as a proxy of trophic status. Applying this modelling approach, we tested our hypothesis that point sources are an equally important nutrient source as diffuse pollution to the Kaiping reservoir. Based primarily on freely accessible data supplemented with information from local managers, we also aimed to demonstrate the potential of our approach for use as a general management tool in China, which may enable managers to quantify necessary nutrient-load reductions and provide new scientific insight into ways of achieving a sustainable and safe water supply. The objective of the present study was thus to identify and quantify major nutrient-input sources linked to surface activity within the Kaiping reservoir watershed and subsequently analyse the potential effects of changes in land-use management and point source contribution on nutrient loading and the implications for the water quality of the reservoir.

Kaiping reservoir is located in the Pearl River (ZhūJiāng) basin, Guangdong province, southern China (22.38° N, 112.68° E), which has a tropical to subtropical climate, with a yearly average air temperature of approximately 22°C. Precipitation is governed by the monsoon system, resulting in a yearly average precipitation of nearly 2000 mm with pronounced dry (October to March) and wet (April to September) seasons. Summer rain events (May to July), associated with typhoons and tropical depressions (Woo et al. 1997), contribute more than 50% of the total precipitation. The reservoir surface area is 27 km² with a watershed of 217 km² and mean and maximum depths of 26 and 35 m, respectively.

To the west, the watershed is surrounded by a mountain chain, creating an elevation difference of more than 1100 m towards the surface of the reservoir. The reservoir has seven tributaries (Fig. 1) (with a total stream length of 84 km delineated by SWAT) that flow through contrasting landscapes exhibiting land uses ranging from intense farming and forestry to areas with a high proportion of forestry and natural bush. In the south, the watershed is dominated by forest with only small proportions of low-intensity agricultural land (tributary 1, Fig. 1). In the north and north-eastern part of the watershed, land use is dominated by rice fields, whereas the north-western part is dominated by pig, goose and duck farming. Intensive forestry (of eucalyptus and pine trees) is also an increasing land-use activity within the watershed, especially in the western part.

Fig. 1.  Study site location, Guangdong Province including SWAT delineated watershed and streams. The reservoir is represented by the hatched area and stations 1 to 4 represent field sampling stations for monthly measurements of water flow and chemistry. The seven tributaries and their draining watershed are labelled 1–7. The watershed is located in both Xinxing and Kaiping administrative districts.

Since 2000, water managers have determined a daily reservoir water balance, including measurements of reservoir water level (converted to storage capacity by hypsographical relations), pan evaporation (Ø 20 cm) and precipitation, as well as recordings of aggregated total outflow (including drinking-water supply, hydropower consumption and irrigation supply) and total inflow calculated as the water balance residual. The same year, the local government initiated monitoring of water quality (e.g. nutrients and chlorophyll a) in the reservoir as part of a province-wide survey program of reservoirs supplying drinking water (Han et al. 2012), which has been supplemented with monthly measurements of water chemistry at four tributary stations since March 2011 (Fig. 1).

The SWAT model, originally developed by the Agricultural Research Service at the US Department of Agriculture (Arnold et al. 1998), is a river basin-scale model and has been applied widely throughout the world (Gassman et al. 2007). The model is governed by the water-balance equation comprising: surface runoff, precipitation, evapotranspiration, infiltration and sub-surface runoff. Nutrient dynamics include both inorganic and organic forms of nitrogen (N) and phosphorus (P). Nutrient seepage depends on the existing pool in the soil and on external loading (e.g. fertilising and point sources), balancing the plant uptake and impacts from soil properties, management and degradation. The model operates with so-called Hydrological Response Units (HRU), which are a unique combination of land cover, soil and topography and assumed to be homogeneous in hydrologic response across the watershed. The SWAT model requires input on soil properties, land use and topography (digital elevation model) as well as climatic conditions distributed across the watershed. In the present study, a combination of open source data and local data on land use, point sources and climate was used (Table 1).

Table 1.  Overview of applied data sources according to SWAT requirements (Neitsch et al. 2009)
Publishers: (1) Consultative Group on International Agricultural Research – Consortium for spatial information (CGIAR-CSI), (2) Harmonised World Soil Database (Nachtergaele et al. 2009), (3) European Space Agency (ESA), (4) local reservoir managers and (5) China Meteorological Information Center (CMA).

Topographical contours of the watershed were derived from the freely available Shuttle Radar Topography Mission (SRTM) 90 m grid DEM (SRTM 2008), which has been used in several published SWAT model applications (e.g. Wu and Chen 2009).

We used the European Space Agency (ESA) GLOBCOVER 2009 dataset (300 m grid) (Arino et al. 2007) to classify land cover within the watershed, which was aggregated to a total of four classifications representing the major land covers within the watershed, with forest (45.6%) and cropland (28.7%) as the dominant features. The GLOBCOVER dataset contains no attribute information about crop and forest types; hence, the land-cover description was refined by non-spatial statistics on watershed land use obtained from local authorities (water managers) (Table 2). Based on these non-spatial statistics, we transferred the GLOBCOVER land cover categories to SWAT-applicable land-use classifications and assigned representative SWAT crop codes (Table 3). Forest (FRST) was split into commercial pine trees (PINE, 66%), mixed ecological forest (FRST, 27%) and commercial eucalyptus (EFST, 7%). Because parts of the GLOBCOVER cropland classifications are fairly broad (Table 3), the area of rice was overestimated compared with the total area of rice recorded for the villages (Table 2). Consequently, we assigned 45% of the rice area to rangeland (representing a mosaic of grassland, shrubland and forest). Data for estimating crop growth parameters for rice was limited and parameters (focusing on radiation use efficiency, harvest index, maximum leaf area index, light extinction coefficient and minimum and optimum growth temperature) were adjusted according to literature values alone (e.g. Liangzhi et al. 1987; Tang et al. 2009). Likewise, eucalyptus (EFST) is not a part of the parameterised forest types in the SWAT database (Neitsch et al. 2009), so literature parameters were used (e.g. Ashton 1975; White et al. 2010).

Table 2.  Watershed statistics including population density, agricultural land use and livestock density distributed at village level (n = 15)
East and north geocoding is projected in WGS84 UTM 49N. Livestock represents the number of animals produced per year based on 2010 statistics

Table 3.  GLOBCOVER land cover classification of the Kaiping watershed and aggregated SWAT land cover representation
Under SWAT code the percentage (%) represents the split ratio in the model.

To encompass the diffuse (non-point source) nutrient loading, we implemented land-use and management schemes (rotations) according to knowledge of local practices in forestry and agriculture, respectively.

In the model, management of commercial pine forestry was implemented as an unfertilised 10-year rotation with initial site preparation and transplanting followed by harvesting and burn-off. For the fast-growing commercial eucalyptus forestry, which is fertilised, we implemented a 5-year rotation cycle in SWAT (including: site preparation, spring transplanting (March) with an initial leaf area index of 0.05 (Wang et al. 2008), fertilisation, harvesting and burn-off). A mixed-compound NPK fertilizer (14–15–16) in doses of 950 kg ha^–1 (assuming densities of 1900 trees ha^–1) was applied according to local practices( i.e. once annually during the first and the second year).

For rice fields, we implemented a double rice-cropping rotation each year, the first growth season being from April to July and the second from late July to end of October. According to local authorities the rice-cropping management in the Kaiping watershed follows in part the principle of Site-Specific Nutrient Management technology (SSNM), recommended by the International Rice Research Institute (IRRI), with the intention of optimise the timing of fertilizer application to reduce the surplus of nutrients on rice fields (Palis et al. 2010). Mixed industrial NPK fertilizer (20–5–18) is typically applied three times per rotation at 10 day intervals after transplanting in proportions of 40, 40 and 20% of the total (570 kg ha^–1) added compound fertilizer, corresponding to an annual N application of 228 kg ha^–1. This management procedure was adopted in the SWAT model, and rice fields were implemented as lowland irrigated rice fields with no constraints on water throughout the entire growth season.

A spatially distributed soil profile description was obtained from standardised soil texture and geology classes derived from the Harmonised World Soil Database (HWSD) and interpreted to SWAT compatible parameters (Table 4) according to HWSD documentation (Nachtergaele et al. 2009) and standards adopted from SWAT input and output documentation (Neitsch et al. 2009). The soil hydrological group was classified from the soil texture class, a procedure also adapted by Stehr et al. (2008). We included two soil layers: a top layer (0–300 mm); and a sub-surface layer (>300 mm). Furthermore, HWSD provided the composition of silt, clay and sand and organic carbon as well as bulk density, whereas rock-fragment content was represented by the HWSD gravel content (the proportional content of particles between 2 mm and 75 mm). Maximum rooting depth and moist soil albedo were estimated from the soil database in SWAT, where maximum rooting depth was assigned as a simple mean per soil hydrological group and moist soil albedo was estimated from a computed exponential regression model including organic carbon content (%) as predictor (Albedo = 0.238exp(–2.02organic carbon content), R²: 0.99). This exponential regression was applied to the organic carbon content between 0–1.17%; for higher values a constant albedo of 0.01 was used. Plant available water content was calculated as the difference between water content at field capacity and water content at the wilting point derived from the soil compositions of sand, silt, clay and organic content (Bouraoui and Aloe 2007). These compositions were also used to estimate saturated hydraulic conductivity from pedo-transfer functions derived by Li et al. (2007).

Table 4.  Initial SWAT soil layer parameterisation obtained from interpretation of soil texture classes in the Harmonised World Soil Database (HWSD) according to Nachtergaele et al. (2009), Neitsch et al. (2009) and Li et al. (2007).

Daily meteorological data for SWAT (precipitation, maximum and minimum air temperature, wind speed, relative humidity and solar radiation) was obtained from two different data source providers and three different gauge locations. Precipitation was provided by local managers and was measured by one gauge located at the edge of the reservoir (22.5° N, 112.42° E). Maximum and minimum air temperature, wind speed and relative humidity were provided by the China Meteorological Information Center (CMA) from their observation station in Yangjiang (21.5° N, 111.58° E, approx. 130 km south-west of the Kaiping reservoir), and solar radiation was obtained from their station in Guangzhou (23.1° N, 113.2 ° E, approx. 100 km north-east of the Kaiping reservoir).

More than 95% of the towns in China have no sewage treatment facilities (Yan et al. 2011) and all domestic sewage within the Kaiping watershed is discharged untreated via open ditches to rivers or streams. In 2010, the watershed had a total population of 34 514 people scattered mainly in 15 villages (Table 2). Based on their spatial location, we linked the villages to the closest corresponding sub-basin delineated in SWAT and calculated a domestic-sewage point source loading per village using the following sewage coefficients per capita: 12 g TN day^–1 and 1 g TP day^–1 (Yan et al. 2011).

Livestock production is mainly allocated to small fenced areas with high stocking density in the flatter regions of the watershed. Detailed information on the spatial distribution of these areas was not available. Therefore, we simply included these as additional point sources for each sub-basin. The production is typically integrated as agriculture–aquaculture (stocking various species of fish, e.g. nile tilapia (Oreochromis niloticus) or carp (Cyprinus carpio)), in semi-intensive systems where livestock manure is applied as fertilizer in fish ponds to produce natural feed for the fish (Guo et al. 2012). These semi-intensive systems are almost completely closed with only minor or no water exchange to surrounding waterways except when drained for batch harvest. Ducks and geese are often produced as described above or occasionally in enclosures partly inside streams or rivers (Fig. 2a). Pig production is also typically in close connection with fish ponds (Fig. 2b, c) where excreta is separated manually into liquid and solid fractions; the latter scraped off the floor, composted and used mainly for vegetable and fruit production. Urine and minor parts of the solid fraction are flushed by channels to serve as nutrient input in fish ponds or discharged directly to waterways. Chicken farms typically discharge directly to rivers and streams or collect excreta for composting.

Fig. 2.  (a) Geese production within the river and at the adjacent banks. (c) Fish pond with shore-situated pig production (b) with direct effluent discharge (holes in the wall) to fish pond.

As livestock manure primarily serves as input to fish production by facilitating productivity through released nutrients, fish ponds retain nutrients and reduce the loadings to waterways depending on local pond characteristics (e.g. fish species, livestock production, pond properties and intensity). However, published information on environmental impacts of closed, semi-intensive, integrated agriculture–aquaculture systems is scarce.

We adopted nutrient-retention factors published by Edwards (1993) to estimate the actual nutrient loading derived from livestock to the Kaiping reservoir. According to Edwards (1993), only a small fraction of the total nutrient loading to semi-intensive fish ponds with static water (as is typical for the Kaiping watershed) is removed (11–15% TN and 6–8% TP) by the harvested fish biomass, whereas pond processes, for example, denitrification (Li et al. 2010), ammonia volatilisation (Hargreaves 1998), phosphorus adsorption to pond mud (Shrestha and Lin 1996) and sedimentation remove a considerably larger portion (83% TN and 84–86% TP). Thereby nutrient loading to receiving waters becomes relatively modest (2–6% and 6–10% of the TN and TP load, respectively, which enters the fish ponds).

Livestock production was not geocoded at individual farm level, but recorded as yearly summarised statistics for each village (n = 15) (Table 2). Consequently, we combined the 2010 statistics on total livestock production with standardised livestock nutrient coefficients (SEPA 2004) (Table 5) and estimated daily nutrient loadings from livestock, which were added as point sources to the corresponding sub-basins in SWAT. The average daily loading was corrected according to the individual livestock breeding days in a year (cow: 365, pig: 190, chicken, duck and geese: 210 days) (SEPA 2004), and to account for influences from pond nutrient retention we applied the pond nutrient retention factor published by Edwards (1993) (80% for TP and 94% for TN) for all duck and geese dung as well as pig and cow urine, whereas pig, cow and chicken dung was assumed to be composted and used elsewhere, and thus not included in the model setup.

Table 5.  Coefficients for livestock manure and urine: quantity and nutrient concentrations according to SEPA (2004)

The SWAT model setup comprised a total of 44 sub-basins including the reservoir. This layout was compiled to obtain an aggregated inflow that could be used for model calibration (effectively the outflow from the sub-basin containing the reservoir), aggregating all individual tributaries to the reservoir.

In areas with high spatial heterogeneity in topography, precipitation and temperature may vary considerably with elevation. SWAT can simulate the effect of lapse rates by inclusion of elevation bands (Neitsch et al. 2005), and we defined 10 elevation bands for each sub-basin and computed the midpoint elevation and respective proportion of sub-basin area. As no appropriate estimations or observations of precipitation lapse rates were available for our study site, we included the precipitation lapse rate as a calibration parameter, whereas the temperature lapse rate was set globally to –6.49 K(°C)/1000 m (Bonell et al. 1993).

Hydrological calibration was performed by monitoring changes in selected performance statistics, including Nash-Sutcliffe efficiency (ENS) (Nash and Sutcliffe 1970), percent bias (PBIAS), root mean square error (RMSE) and coefficient of determination (R²) (Moriasi et al. 2007), and by visual inspection of simulated model output against observed data (timing of peaks and base flow). For hydrological calibration, the simulated total daily inflow to the reservoir (Fig. 1), representing the overall contribution from all inflowing sub-basins, was compared with observed total daily inflow to the reservoir estimated as the residual in the reservoir water balance encompassing measured precipitation, evaporation, storage capacity and outflow collected by reservoir managers. We used the years 1998 to 2000 as a model warm up before calibration from 2001 to 2005 and validated the model performance from 2006 to 2010. Following the initial manual calibration, a Sequential Uncertainty Fitting Algorithm (SUFI2) (see Abbaspour et al. (2007)) that uses a global search procedure through Latin Hypercube Sampling was applied as a second step, where parameters were adjusted automatically and repeatedly through numerous simulations (Schuol et al. 2008).

Nutrient loadings (TN and TP) were calibrated against data from the four stations for the period (March 2011 to August 2011) by adjusting parameters (Table 6) reported to be of importance for the model nutrient performance (e.g. David et al. 2009; Jha et al. 2010; Gong et al. 2011). Due to data limitations, the model was only evaluated for PBIAS (i.e. model performance in terms of levels) and no independent validation was conducted. Subsequent to calibration, we examined the sensitivity of key model parameters and input data in relation to the predicted nutrient loading to the Kaiping reservoir to be able to identify the most critical sources of uncertainty. Here, model parameters and livestock nutrient-retention factors for fish ponds were changed equally (±25%) one at a time, and sensitivity was expressed as the difference (%) in total nutrient (TN and TP) loading to the reservoir relative to the calibrated baseline scenario.

Table 6.  Parameters included in calibration of nutrient loadings from the watershed
Allocation of sub-basins: Tributary 1 (i.e. sub-basin 23–25, 27–44) and tributary 6 (i.e. sub-basin 1–8, 11)

Besides the parameters included in calibration, the canopy storage, which may affect the balance between infiltration and surface runoff, was set to 1.5, 3 and 5 mm for rice, shrubland and forest, respectively, according to Chen et al. (2005). Moreover, the average slope steepness of rice fields, which may impact the model replication of the ‘shoulder’ of peak events, was adjusted homogeneously to 0.005 for all HRUs defined as rice.

Analysis of the potential effects of changes in land-use management and point source contributions of nutrient loadings to the reservoir was made through a series of model scenarios (Table 7). Due to uncertainties about loading regimes (i.e. livestock production connectivity to aquacultures) and actual nutrient retention through the fish ponds in the watershed, we focussed comprehensively on scenarios where livestock contributions and nutrient retention coefficients were changed. Moreover, we included two scenarios with modified fertilizer application rates to rice fields (i.e. 556 and 805 kg N ha^–1 yr^–1) according to the reported average and maximum annual application rates in southern China (Hatano et al. 2002). Buffer zones, a widely used measure to minimise water quality deterioration associated with adjacent land-use activity, were also included as a scenario and simulated as 10-m zones (using the SWAT filter width parameter, FILTERW) around all rivers in the watershed.

Table 7.  Overview and description of scenarios (SC) run to examine impacts on reservoir nutrient loadings
X indicates the parts of the point sources included in the scenarios

Impacts of scenarios were compared using the simulated total average annual nutrient loading (TN and TP) to the reservoir from 2000 to 2010. Definition of a TN : TP ratio-threshold to indicate P or N limitation water bodies is reported to vary widely (Abell et al. 2010). The mean TN : TP ratio, derived from average annual nutrient concentrations sampled between 2001 to 2010 in the reservoir, was 33 (range: 16–52), and therefore the primary producers are likely dominated by P limitation (i.e. TN : TP >15, Abell et al. 2010) for the majority of time, while N limitation may occasionally occur. Hence, we additionally estimated the effects on the water quality of the Kaiping reservoir using the simple and widely used empirical model by Vollenweider and Kerekes (1982):

P_reservoir represents the internal annual average reservoir phosphorus concentration (µg L^–1), P_in is the annual average phosphorus concentration of inflowing water (µg L^–1) and R_t the reservoir water retention time estimated to 0.37 years (136 days) on average for the Kaiping reservoir between 2000 and 2010. The simulated phosphorus concentration of the reservoir was subsequently used as a proxy for reservoir trophic status and drinking water quality (MEP 2002).

Generally, dynamics in runoff were adequately represented and both peak events as well as base flow were captured by the calibrated SWAT model (Fig. 3). Several parameters were changed during calibration (Table 8) and, statistically, the model application simulated daily inflow with slightly better performance for the calibration period than for the validation period and gave a better representation of the monthly average flow (Table 9).

Fig. 3.  Calibration (2001–2005) and validation (2006–2010) plots representing daily simulated total flow from the watershed by SWAT (black) and daily total observed inflow (grey) estimated as the residual in a reservoir water balance encompassing measured water level, precipitation, evaporation and outflow.

Table 8.  Parameters included during calibration of surface runoff for the watershed
Default values are obtained from the SWAT databases, whereas soil parameters (i.e. SOL_K and SOL_AWC) are obtained from Table 4

Table 9.  Calibration and validation statistics for hydrology of comparison of simulated total flow and daily total observed inflow estimated as the residual in a reservoir water balance encompassing measured water level, precipitation, evaporation and outflow
ENS: Nash‐Sutcliffe efficiency coefficient; PBIAS: percentage bias; R²: coefficient of determination and RMSE: root mean square error

Simulated flow (m³ s^–1) and daily loadings (kg) of TP and TN were calibrated according to observed measurements at four stations collected on a monthly basis (Fig. 4) and several parameters were changed from the default values (Table 6). Nutrient loadings of TP (PBIAS: 36%) and TN (PBIAS: 18%) were generally underestimated by the model, although with distinct differences between stations.

Fig. 4.  Observation (circle) and simulation (black line) of flow (m³ s^–1), TN and TP (kg day^–1) at the four stations (1 to 4, Fig. 1). Observations are sampled monthly from March 2011 to August 2011, whereas simulation is daily output from SWAT.

The annual nutrient non-point source exports within the watershed varied between 3.4–14.1 kg N ha^–1 and 0.1–0.8 kg P ha^–1. Nutrient export at sub-basin level showed higher exports in the northern part of the watershed (Fig. 5), which is primarily cultivated with rice fields, whereas the southern part, typically forested, had lower exports, although steeper areas towards the south-west exhibited elevated TP exports.

Fig. 5.  Annual mean simulated non-point source related export of TN and TP (kg ha^–1) for each sub-basin (n = 44) situated in the sub-watersheds of the main tributaries (1 to 7) simulated by SWAT from 2001 to 2010.

An examination of the nutrient loadings from the major tributaries separated into diffuse and point sources showed substantial loads of TP from point sources (Fig. 6), whereas TN was predominantly related to diffuse nutrient exports (especially for tributary 1). Tributary 4 and 7 had no registered point sources. Generally, nutrient loadings from livestock production were minor compared with sewage-related loading, especially for TN.

Fig. 6.  Simulated annual TN and TP loading from the reservoir tributaries as averages from 2001 to 2010 separated into diffuse, livestock and sewage loading.

The sensitivity analysis revealed a contrasting impact of model parameters and input data on loadings of TN and TP, respectively (Fig. 7). TN was more sensitive to changes in SWAT process parameters (denitrification exponential coefficient, CDN and nitrogen percolation coefficient, NPERCO) than changes in contributions from point sources (POINT) and nutrient reductions through fish ponds (REDUC), respectively, whereas the TP loading was more strongly affected by changes in total point source loading (POINT) and nutrient reductions through fish ponds (REDUC) than SWAT process parameters (phosphorus soil partitioning coefficient, PHOSKD and phosphorus percolation coefficient, PPERCO).

Fig. 7.  Sensitivity analysis for the SWAT parameters: CDN, NPERCO, PHOSKD and PPERCO (Table 6) and point source contributions with changing total point source loading (POINT) and nutrient reductions through fish ponds (REDUC). All variables were adjusted +/–25% individually (while keeping all other parameters constant) and evaluated according to the proportional change in average annual nutrient loading with respect to the baseline scenario annual loading.

The annual mean in-lake TP concentration calculated from the water retention time in the reservoir and the mean annual inflow TP concentration (Eqn 1) were generally overestimated by the model (PBIAS: –25%); while predictions demonstrated only minor year-to-year variation, the observations exhibited greater variability with the lowest and highest concentrations in 2004 and 2006, respectively (Fig. 8).

Fig. 8.  Annual reservoir TP concentration in the Kaiping reservoir. Observed represents annual concentrations from interpolated routine measurements every two months in the reservoir. Baseline scenario represents predicted concentrations according to Vollenweider and Kerekes (1982) from annual simulated SWAT loadings and reservoir retention time. Class (I–IV) represents drinking water quality classes according to MEP (2002) and trophic states according to Vollenweider and Kerekes (1982).

Simulated watershed scenarios (referred to as scenario (SC) 1 to 11), Table 7) revealed a considerably different influence on nutrient loadings (tonnes year^–1) to the reservoir (Fig. 9) and, consequently, markedly different impacts on water quality ranging from Class I to IV (Chinese surface water quality standards) (Fig. 10). The lowest TP reservoir concentration was achieved by removing all point sources (SC2), whereas scenarios on altered rice-field fertilizer usage (SC10 and SC11) and zero nutrient retention through fish ponds (SC3) resulted in elevated concentrations, albeit the highest concentration was reached when simulating the effects of discharging chicken manure to fish ponds (SC7). Generally, reservoir TP concentrations only varied marginally irrespective of whether all livestock (SC8), pig (SC6), or duck and geese nutrient loadings (SC5) were excluded, while the implementation of buffer zones along rivers (SC9) reduced reservoir TP concentrations with the same magnitude (31%) as the removal of sewage discharges (SC4) (34%).

Fig. 9.  Total TN and TP loading (tonnes year^–1) to the reservoir, representing an average from 2001 to 2010 for each of the simulated scenarios (SC) 1 to 11. The dashed line shows the mean simulated results for the baseline (Base) scenario. Error bars represent one standard deviation (s.d.).

Fig. 10.  Annual mean reservoir TP concentration predicted according to Vollenweider and Kerekes (1982) from the inflow concentration of TP and water retention time in the reservoir. Each dot represents different scenarios according to Table 7. Error bars represent standard deviation (s.d.) for 10-year TP predictions (2001–2010). Class (I–IV) represents drinking water quality standards according to MEP (2002) and trophic states according to Vollenweider and Kerekes (1982).

The model exercises provided support our hypothesis that point sources rather than diffuse pollution constitute a considerable proportion of the external nutrient loads to the reservoir. They further demonstrated how nutrient load reductions can be accomplished; for example, through land‐use changes, removal of point sources or other targeted measures. Based on the Vollenweider and Kerekes (1982) model, we estimated a critical load of TP to 13.5 tonnes yr^–1 to comply with the minimum drinking water requirements (Class II) (MEP 2002), which corresponds to 87% of the simulated base loading (Base, Fig. 10). Predicted reservoir TP concentrations were slightly overestimated and year-to-year dynamics were not well captured (Fig. 8). However, temporal variations in livestock management and changes in cropland management were not included in the model setup. Moreover, the Vollenweider and Kerekes (1982) equation (Eqn 1) is mainly based on temperate lakes and reservoirs and may be less appropriate for predictions in subtropical reservoirs. Young reservoirs have a higher retention capacity than that typical for much older lakes (Cooke et al. 2005), which may explain the predicted higher TP. However, the subtropical climate may counteract such effects as higher temperature and shifts in trophic structure may lead to enhanced internal-P loading (Jeppesen et al. 2009). To circumvent these uncertainties, a complex hydrodynamic and ecological lake model is currently being set up to obtain a better representation of reservoir dynamics.

The loading reduction simulations (Table 7) showed markedly different impacts on water quality, ranging from Class I to IV, which emphasises the importance of adequate estimation of nutrient losses from livestock production and connectivity to adjacent waterways. Interestingly, the implementation of buffer zones along rivers (SC9) reduced reservoir TP concentrations with the same magnitude as the removal of sewage discharges (SC4), indicating the high importance of surface runoff, which may be partly trapped in buffer zones.

The present study did not encompass a socioeconomic analysis aiming at supporting management decision-making in terms of cost-effectiveness evaluation. However, we argue that awareness of upstream‐downstream linkage (Lerner and Zheng 2011) within the watershed may serve as a window of opportunity for action planning. Several adjacent or in-river livestock production units (primarily geese and ducks, Fig. 2a) are scattered along waterways in the northern part of the watershed. To improve water quality, these productions could be transferred to on-shore situated facilities with local remediation of effluents before discharge. Lowering nutrient exports from the numerous joint livestock‐fish production systems scattered across the watershed may also be achieved by retaining a minimum water level (e.g. 10 cm) in the ponds during batch harvesting, as suggested by Lin and Yi (2003).

Based on simulations, managers should also focus on the eucalyptus plantations that are often situated at steep locations in the watershed. Our mapping of diffuse nutrient exports (Fig. 5) found higher nutrient exports from the steep south-western parts compared with the more flat south-eastern part of the watershed draining to tributary 1, while the land-cover map justified no exclusive scenario analysis on eucalyptus. Eucalyptus plantations potentially have great exports of sediment and nutrients during heavy rain events, especially during the first years after establishment. The current rotation period could therefore be extended from 4–6 years to 6–8 years in order to lower the periods of possibly severe soil washout (Wei and Xu 2003). Alternatively, and as a continuation of the current Sloped Land Conversion Program (Xu et al. 2006), strategic planning of the spatial location of ecological forests (Chokkalingam et al. 2006), where the focus is on the steepest parts of the landscape (Lerner and Zheng 2011), may also provide long-term reductions of sediment and nutrient losses.

Nutrient losses associated with town living (i.e. sewage discharges) were clearly of major importance (Fig. 6), and implementation of sewage-treatment facilities could considerably improve water quality (Fig. 10, SC4). Several existing technologies are available to address sewage control from rural areas (e.g. Shen 2008) and point sources are somewhat easier to target than losses from non-point sources (Huang and Xia 2001). Exemplified by our simulations of scenarios with changed nutrient-reduction coefficients and inclusion of chicken manure (i.e. SC3 and SC7), which implicitly represent the effects of an increase in animal production, focus should, however, also be directed at non-point sources related to agricultural land use to ensure high-intensity food production with low environmental impact. This is highly relevant considering the expected increase in population growth and living standards, the global shifts in markets prices (Chen et al. 2008; Khan et al. 2009) and climate change (Tao et al. 2009).

Compared with the range of reported performance statistics from various SWAT model applications summarised by Gassman et al. (2007) and criteria recommendations (i.e. ENS and R² > 0.5 for monthly data) by Moriasi et al. (2007), the model performed satisfactorily for the daily inflow of water to the reservoir, with a slightly better performance during the calibration (ENS: 0.75 and R²: 0.75) than the validation (ENS: 0.63 and R²: 0.70). Representation of monthly average flow was also satisfactory (ENS: 0.80–0.91), where ENS values in other Chinese studies have been reported to range between 0.63 and 0.93 for the monthly average flow (e.g. Li et al. 2011).

During calibration of hydrology several parameters were changed specifically for each soil layer and land use type, but also globally (Table 8). The parameters were generally within the recommended ranges, but the Manning’s roughness coefficient was set higher than usual (i.e. <0.8 (Stehr et al. 2010)). The latter may be attributed to the various modified flow paths (channels and ditches) associated with land use such as rice fields, constituting a predominant part of the watershed (Table 3), which influence runoff regimes differently than conceptually represented in SWAT. These artificial flow paths may also explain why hydraulic conductivities needed to be markedly modified during calibration (Table 8).

Although substantial adjustments of both the runoff curve number and overland Manning’s roughness coefficient proved necessary, peak events were generally overestimated (Fig. 3), in contrast to the study by Xie and Cui (2011). They applied an updated modelling code routine to improve the SWAT representation of rice fields and found that the original code (as used in the present study) generally underestimated the water flow during peak events. It is reasonable to assume that the numerous fish ponds scattered in the Kaiping reservoir watershed, yet too small to be captured by the relatively coarse land use map (300 × 300 m), may serve as a water buffer with considerable impacts on peak runoff, for which compensation is necessary.

For all years, flow peaks occasionally occurred in the observed data on days where no precipitation was registered at the gauge station. Unavoidably, this constrains an even better model representation of observed data and lowers the performance statistics. As pointed out by Gassman et al. (2007) weak performance statistics may be attributed to questionable rendering of spatial heterogeneity in precipitation during rainfall events because of an insufficient network of rainfall gauges. With only one precipitation gauge in the entire Kaiping reservoir watershed, this may explain the difficulties and sometimes inadequate representation of flow by the model. It should be kept in mind, however, that due to lack of data observed flow is not the actual measured stream flow, which is generally the case in SWAT model applications (Moriasi et al. 2007), but a representation of the total inflow to the reservoir estimated from a reservoir water balance. Hence, uncertainties in both measurements of evaporation and estimations of storage capacity and outflow are aggregated in the residual, representing the observed inflow and, thereby, the data used for calibration and validation. Discharges when emptying fish ponds during batch harvests may also contribute to deviations between observed and simulated flow.

As experienced by others (e.g. Wu and Chen 2009; Gong et al. 2011), validation of modelled nutrient-load dynamics, especially peak events initiated by storm flow, proved difficult from the available discrete monthly monitoring program conducted at the four stations (Fig. 4), and the results of simulation of the overall nutrient loading were therefore modest. For station 4 (and 3) the model moderately replicated the magnitude and dynamics of TN, but lacked the capability to capture an observed peak loading in May (Fig. 4). This peak loading, however, occurred on a dry day (no rain), which may explain the model deviation. For this particular peak loading, available data also showed that a fairly high fraction of TN was constituted by ammonia (station 3: 60% and station 4: 80%), which may indicate discrete point source discharges not encompassed by the model input data.

We tried to enhance model performance by changing primarily basin-wide parameters affecting nutrients, but still the model representation of observed data can be further improved. By default, the soil moisture condition threshold (SDNCO) for initialisation of denitrification is set to 110% of field capacity, which conceptually means that the process is inactive in SWAT. Hence, SDNCO was changed to allow denitrification processes to begin at lower moisture conditions (95%) (Pohlert et al. 2007). SDNCO was set slightly higher than reported as in situ values (ranging between 50–83%) (Barton et al. 1999), because denitrification is initiated before vertical water movement in SWAT when SDNCO is set below field capacity (Pohlert et al. 2007), which might lead to overestimated nitrogen losses by denitrification. Moreover, a higher threshold may moderately influence the timing of nitrogen dynamics.

The SWAT denitrification exponential coefficient (CDN) had a marked impact on simulated nitrogen loads (Fig. 7). As in other studies (e.g. David et al. 2009), we lowered the CDN (Table 6) with respect to the default SWAT parameterisation, resulting in an average annual simulated denitrification rate of 44 kg N ha^–1 across all mosaic land uses in the watershed. Studies applying SWAT have reported considerably higher simulated annual denitrification rates (e.g. Pohlert et al. (2007): 135 kg N ha^–1and Schoumans et al. (2009): 236 kg N ha^–1), but as demonstrated by Pohlert et al. (2007) SWAT may overestimate denitrification due to the model-process routine. Losses of nitrogen by denitrification are influenced by various factors (for rice paddies, for instance, the fertilizer-application method, timing and applied nitrogen compound) (Craswell et al. 1981), which challenges comparisons between land-use management practices. Nevertheless, the denitrification rate simulated in the present study lies within the combined span of annual soil denitrification rates for forests (<0.1–40 kg N ha^–1) and agriculture (0–239 kg N ha^–1) reported by Barton et al. (1999), where irrigated and inorganic fertilised soils tend to exhibit higher rates (ranging from 49 to 239 kg N ha^–1). Moreover, our simulated ranges of annual diffuse nutrient exports from the watershed (3.4–14.1 kg N ha^–1and 0.1–0.8 kg P ha^–1, Fig. 5) were within the range of observed diffuse exports from watersheds across Europe, ranging from 0.9–25.7 kg N ha^–1 and 0.02–3.14 kg P ha^–1, respectively (Bouraoui et al. 2009). However, considering the importance of CDN for the model outcome (Fig. 7), further examination of the overall SWAT denitrification process as well as in situ estimates of denitrification are needed in order to validate the modelled losses of nitrogen by denitrification (Pohlert et al. 2007), not least for subtropical watersheds where denitrification potentially is important for the overall nitrogen exports owing to, for example, enhanced mineralisation rates at higher temperatures (Neitsch et al. 2005) and pronounced anaerobic conditions in flooded paddies.

Encompassing only freely available datasets, the spatial resolution of our model somewhat constrained the information on the distribution of nutrient exports. Evidently, a detailed collection of watershed characteristics may improve the representation of local processes. For instance, even though we selected pedo-transfer functions developed for a Chinese soil region (Li et al. 2007), the estimated hydraulic conductivities of the representative soils may still be uncertain because of inconsistency between the study site and the pedo-transfer function site of origin. Moreover, identification of ‘hot spots’ to mitigate nutrient losses associated with specific soil types may also prove difficult (Jha et al. 2010) because soil map coarseness (800 × 800 m grid) may level out spatial differences in nutrient seepage. Hence, enhanced local parameterisation of Chinese soil types and land uses as well as inclusion of high-resolution land-use classifications, specifically the mapping of rice fields and absolute spatial separation between forest types, may improve the model and its output. Additionally, uncertainties regarding modelled nutrient losses may be targeted in accordance with our analysis of sensitivity (Fig. 7) and ideally lowered through validation of the SWAT representation of denitrification for TN losses and better estimation of nutrient inputs for TP losses. As described by Ongley et al. (2010), the application of SWAT is challenged by the contrasting Chinese watershed properties relative to the North American watersheds for which the model was originally developed. Fundamentally, applications should therefore be viewed critically as to the appropriateness for Chinese watersheds. For southern Chinese watersheds especially, we also recommend a high frequency of monitoring of flow and nutrients in tributaries during the pronounced wet season (for Kaiping: June to August, Fig. 3) to allow better model calibration and validation of peak events that potentially transport a considerable proportion of the total nutrient export (Tang et al. 2008). Moreover, a denser rainfall gauge network should be considered in order to capture spatial heterogeneities in precipitation.

The landscape, land use, population density and climatological conditions of China vary markedly, but problems with water pollution and eutrophication are of nationwide concern (Chen et al. 2008). Hence, we argue that the conceptual approach presented in the present study can serve as a valuable tool for local managers of widely different watersheds to improve understanding and estimation of contributions from spatially distributed pollutant sources. Additionally, it may help them identify vulnerable areas (as exemplified in Fig. 5) regarding nutrient losses to aquatic systems adversely impacting reservoir water quality. Based on watershed models as that applied in the present study, such tools may also allow spatial differentiation of crop-specific fertilizer norms, which is common in several Western countries (Kronvang et al. 2008). Awareness of spatial vulnerability may improve the layout and selection of areas for implementation of measures such as ditches and wetlands (Jiang et al. 2007) that locally lower agricultural non-point source nutrient exports – as called for by Tang et al. (2012). Thus, control of agricultural intensity according to vulnerability of the land may preserve the overall yields at watershed scale, while lowering nutrient losses and thereby better balancing food security and the quality of the water supply. Because crops and management practices, as well as forestry (Ju et al. 2004; Chokkalingam et al. 2006), may vary significantly across China, elucidating the importance of agricultural non-point sources and suitable mitigation methods for other watersheds will require local application of watershed models as demonstrated by the approach developed in the present study.