Modeling phosphorus sources and transport in a headwater catchment with rapid agricultural expansion

Environmental Pollution 255 (2019) 113273

Contents lists available at ScienceDirect

Environmental Pollution journal homepage: www.elsevier.com/locate/envpol

Modeling phosphorus sources and transport in a headwater catchment with rapid agricultural expansion* Wangshou Zhang a, Steven G. Pueppke b, c, Hengpeng Li a, *, Jianwei Geng a, Yaqin Diao a, David W. Hyndman d a

Key Laboratory of Watershed Geographic Sciences, Nanjing Institute of Geography and Limnology, Chinese Academy of Sciences, Nanjing 210008, China Center for Global Change and Earth Observations, Michigan State University, East Lansing, MI 48824, USA Asia Hub, Nanjing Agricultural University, Nanjing 210095, China d Department of Earth and Environmental Sciences, Michigan State University, East Lansing, MI, 48854, USA b c

a r t i c l e i n f o

a b s t r a c t

Article history: Received 3 July 2019 Received in revised form 27 August 2019 Accepted 17 September 2019 Available online 22 September 2019

Increasing riverine phosphorus (P) levels in headwaters due to expanded and intensified human activities are worldwide concerns, because P is a well-known limiting nutrient for freshwater eutrophication. Here we adopt the conceptual framework of the SPAtially Referenced Regressions On Watershed attributes (SPARROW) model to describe total phosphorus (TP) sources and transport in a headwater watershed undergoing rapid agricultural expansion in the upper Taihu Lake Basin, China. Our models, which include variables for land cover, river length, runoff depth, and pond density, explain 94% of the spatio-temporal variability in TP loads. Agricultural lands contribute the largest percentage (61%) of the TP loads delivered downstream, followed by forestland (21%) and urban land (18%). Future agricultural expansion to 15% of the total basin area is possible, which could lead to a 50% increase in TP loads. According to our analysis, an average of 24% of the total P export from the watershed landscape was intercepted in ponds. The exported amount was subsequently retained by tributaries and along the mainstem river, accounting for 14% and 43% of their inflowing loads, respectively. The remaining ~6 tons yr
1 of TP was eventually transported into Tianmu Lake, in Southeastern China. The model identified several sub-catchments as hotspots of TP loss and thus logical sites for targeted management. Our study underscores the significance of agricultural expansion as a factor that can exacerbate headwater TP pollution, highlighting the importance of landscapes to buffer TP losses from sensitive hilly catchments. This also points to a need for an integrated management strategy that considers the spatial-varying P sources and associated transport of TP in precious headwater resources. © 2019 Elsevier Ltd. All rights reserved.

Keywords: Water quality Headwater Phosphorus retention Pond Taihu Lake SPARROW

1. Introduction Phosphorus (P) originating from anthropogenic activities has been an environmental issue for decades, because of its wellknown role as a pollutant of surface waters, that promotes eutrophication (Schindler et al., 2016; Macintosh et al., 2018). Widespread socio-economic development including expanded and intensified agricultural activities and rapid urbanization is now occurring in fragile headwater catchments, resulting in significant land cover changes and enhanced anthropogenic P inputs (Chen,

* €rg Rinklebe. This paper has been recommended for acceptance by Dr. Jo * Corresponding author. E-mail address: [email protected] (H. Li).

https://doi.org/10.1016/j.envpol.2019.113273 0269-7491/© 2019 Elsevier Ltd. All rights reserved.

2007; Hodgkinson and Withers, 2007; Zhang et al., 2015). Mobilization and transport of P from lands in these catchments is exacerbated by heavy rainfall events and landscape features such as steep terrain and thin, easily erodible soils (Yan and Zhang, 2003; Withers and Jarvie, 2008). Transport of these over-enriched P flows into downstream lakes, reservoirs, and other open waterbodies leads to significant environmental damage including excessive algal blooms (Conley et al., 2009; Gao et al., 2015). Accelerated land use change and increasing cultivation in fragile headwater watersheds continues to be one of the most prevalent problems in a world with limited land resources and a growing population (Tilman, 1999; Laurance et al., 2014; Zhang et al., 2018). Any increase of human-disturbed lands (i.e., agricultural and urban lands) contributes to the enrichment of P in surface waters (Ahearn et al., 2005; Diamantini et al., 2018), but many questions remain

2

W. Zhang et al. / Environmental Pollution 255 (2019) 113273

unresolved, among them the origin and transport of P in complicated headwater catchments, especially those with steep topography and natural lowland ponds. Given the fact that headwater catchments provide essential freshwater for growing populations (Alexander et al., 2007), it is critical to understand P sources and transport to develop management strategies that minimize human-induced environmental damage. Deterministic and empirical watershed models are powerful tools to investigate P origin and transport (Nasr et al., 2007; Ongley et al., 2010), but choosing a reliable model to account for the variety of sources, attenuation processes, and the nature of their interactions is challenging. In general, model selection is influenced by model complexity, data availability, user’s objectives, and costbenefit analysis. Although some deterministic models, including HSPF (Wang et al., 2015), SWAT (Radcliffe et al., 2009), and HYPE (Donnelly et al., 2015), can describe the hydrologic and biogeochemical processes that govern P loss in fine detail, their complexity and intensive data requirements may limit their application (Jackson-Blake et al., 2017). Such models often preclude identification of dominant processes operating within a watershed, as they can be masked by uncertainties across a wide range of calibrated parameters (Jakeman et al., 2006; Kirchner, 2006; Zhang et al., 2018). Some simple empirical models can also be useful to describe P transport. For example, the widely used mass-balance estimates of watershed P budgets can portray the direct responses of riverine exports to P sources (Yan et al., 1999; Seitzinger et al., 2005; Mayorga et al., 2010; Swaney et al., 2012). However, such models are generally applied to large basins, and they do not explicitly include spatially detailed and time-varying information that is usually crucial to characterize P transport processes (Mayorga et al., 2010; Zhang et al., 2015; Hong et al., 2017). Applying such models to smaller scale headwater catchments can impair the accuracy and utility of model predictions, because the derived results are often too coarse to be representative and useful. The hybrid statistical SPARROW (SPAtially Referenced Regression On Watershed attributes) model spatially links watershed attributes to rivers and uses simple functions to describe nutrient loads, yields, and deliveries at landscape and regional scales (Alexander et al., 2002; Schwarz et al., 2006); it thus addresses many shortcomings of some deterministic and empirical models. Models such as SPARROW are simple yet powerful to obtain useful information on the distribution, sources, and transport of nutrients in streams and rivers (Saad et al., 2011). For example, the model inputs are readily available and the model structure is relatively transparent, as the model imposes parsimonious mass-balance constraints to empirically estimate terrestrial and aquatic rates of nutrient flux (Alexander et al., 2002; Schwarz et al., 2006). Evidence of TP simulation using SPARROW is recognized to be reasonably accurate, with interpretable model coefficients and estimates of source contributions (Alexander et al., 2004; García et al., 2011; Morales-Marín et al., 2015, 2018). Regression-based SPARROW model fitting is also scaleindependent, making it suitable for diverse spatial and temporal settings (Zhang et al., 2019). Although previous applications often involved relatively large basins over limited time frames, recent applications confirm that it has the potential to simulate P transport on much finer spatial and time scales (Aguilera et al., 2012, 2013; Zhang et al., 2019). Such efforts indicate that the key factors limiting the expansion of SPARROW are associated with insufficient temporal and spatial data, especially those with heterogeneous properties in space and time. If required model inputs are available, adapting SPARROW to relatively small headwater watersheds for time-varying information should be feasible. SPARROW is flexible in its model structure and can

accommodate region-specific processes that complicate other models (Alexander et al., 2000; Brakebill et al., 2010), such as SWAT (Wang et al., 2008) and HSPF (Said et al., 2010; Golden et al., 2014). Pond retention, an important consideration in headwater watersheds with hilly terrain with many lowland ponds, is one of these processes. Such spatially distributed ponds receive surface and subsurface flows from uplands where fertilizers are often intensively applied (Zhang et al., 2019). The strong ability of ponds to trap P loads by sedimentation, biological uptake, and other biogeochemical processes has been documented (Yan et al., 1998; Sønderup et al., 2016). Although these studies highlight the importance of pond retention (Chen et al., 2019), there are few modeling efforts that explicitly include watershed TP retention occurring within spatially distributed ponds. In this study, we use the conceptual framework of the SPARROW model to describe TP origin and transport across space and time, including the roles of lowland ponds in P movement. Our study area is the Tianmu lake watershed (TLW), which is an important headwater watershed of the Taihu Lake Basin in Eastern China. The TLW encompasses a relatively small area (153 km2), with more than 50% of the land covered by forests. Driven by economic incentives, nearly 7% of these forests have recently been converted into agricultural landsdmostly tea plantations (Zhang et al., 2019). These lands often receive excessive applications of fertilizer. Since many of the fields are on steep slopes, soil erosion and nutrient runoff (especially of P) can be severe (Li et al., 2013). Notably, the hilly terrain of this watershed has created a dense network of 1025 lowland ponds that have potential to buffer P movement from fields into the lake. Here, we use the SPARROW model to quantify TP source and transport in the TLW and explore the buffering roles of lowland ponds. We first report the model performance and predictions and then evaluate trends of TP pollution under dynamic land use change scenarios, such as a 1.5-fold increase in the area of tea plantations from the current 10%e25%, and an increase in urban land from the current 7%e11%. We also discuss potential management strategies to alleviate P pollution. These applications of SPARROW are important, because they: (1) document the efficacy of the model for simulating TP loads from time-varying sources due to rapid agricultural expansion, (2) lead to a better understanding of landscape processes and their influence TP transport at fine spatial scales, and (3) provide a useful framework to address the important issue of rapid land cover changes and related scenarios. 2. Materials and methods 2.1. Study area The TLW (Fig. 1) contributes large water yields (~500 mm yr
1) to downstream rivers and is part of the headwater watershed of Taihu Lake, which has suffered from frequent harmful cyanobacterial blooms causing it to be termed a “green monster” (Guo, 2007). Precipitation in the typical subtropical humid monsoon climate is highly variable from year to year, ranging from about 1023 to 2410 mm yr
1 during 2010e2017; more than 50% of precipitation is concentrated in the flood season from June through September. The main land use type in 2017 was forest (54% of the total watershed area), followed by cropland (16%), tea plantations (10%), waterbodies (9%), and urban land (7%). Tianmu Lake is drained by three main river systems (Fig. 1). The Zhongtian, Xiasong, and Pingqiao Rivers drain sub-basins with 51, 48, and 11 km2 area respectively, in the southern portion of the TLW; several small streams in the Linhu sub-basin drain an additional 43 km2. The lake water quality is of great significance because it is the primary drinking water source for Liyang City and

W. Zhang et al. / Environmental Pollution 255 (2019) 113273

3

Fig. 1. The location, land cover (2017), elevation, and monitoring sites of the Tianmu lake watershed. The small semicircular blue areas denote ponds in the watershed. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

is a major tourist attraction. Rapid agricultural and urban development has exerted significant environmental pressure on the basin’s water and surrounding hilly landscape (Li et al., 2013). Remote sensing data shows that more than 7% of the surrounding forests were removed between 2010 and 2017 and replaced with tea plantations (Zhang et al., 2019). Recent monitoring (Nie et al., 2015; Sun et al., 2018) indicates that the water quality in Tianmu Lake has dramatically deteriorated relative to the period before 2010 (Gao et al., 2009). This threatens to undo previous progress controlling cyanobacterial blooms in the downstream Taihu Lake (Stone, 2011). There is a clear and urgent need to protect the headwaters, which mandates an improved understanding of P sources and transport. Such efforts provide a sound foundation to develop and implement effective management and conservation policies. In this study, 68 sub-catchments that were delineated using 5-m DEM data, are used to pinpoint incremental P sources and loads. Nineteen of these sub-catchments are used to calibrate and validate the SPARROW model (see Section 2.6). 2.2. The SPARROW model SPARROW is a watershed model that uses a hybrid statistical and process-based approach to relate water-quality measurements from a network of monitoring stations to the upstream sources and attributes of the watersheds containing the stations (Schwarz et al.,

2006). The process-based specification of the model is well suited to identify factors that affect water quality across spatial and temporal scales (Morales-Marín et al., 2015), and to describe nonconservative transport of contaminants in terrestrial and aquatic ecosystems (García et al., 2011). The model structure, supporting equations, and details on the theoretical development are given in Schwarz et al. (2006). Briefly, a simplified mathematical formula can be expressed as:

Fi ¼ Ai Fi
1 þ Ai

0

N X n¼1

gn Ln i exp

M X

!

dm Zmi

(1)

m¼1

where Fi is the annual total phosphorus (TP) load for catchment i, which is defined as TP inflows from upstream reaches plus the load originated within the incremental catchment of the evaluated reach. This first summation term, Ai Fi
1 , denotes the flux derived from upstream TP inflows (Fi
1 ) after in-stream retention, and if i ¼ 1, then Fi
1 ¼ 0. Ai is the in-stream delivery function representing attenuation processes acting on the flux as it travels along the river pathway from catchment i-1 to catchment I; Ai can thus be viewed as a mainstem delivery term. The second summation term represents the amount of flux derived from TP sources in the reach’s incremental catchment and delivered to the catchment outlet after land-to-water and instream delivery. The source term (gn Ln i ) is composed of an array of source variables, Ln i , where n ranges from 1 to the total number

4

W. Zhang et al. / Environmental Pollution 255 (2019) 113273

of land cover types (N), and a vector of source-specific coefficients for each land cover, gn , which is obtained from automated parameter estimation. The land-to-water delivery term P (expð M m¼1 dm Zmi Þ) is composed of a vector of estimated coefficients (dm ) and an array of watershed attributes, Zmi , where m ranges from 1 to the total number of variables (M). An in-stream delivery term (Ai 0 ) represents in-stream transport within the incremental catchment. Here, Ai 0 can be thought of as an in-tributary delivery factor, because the loss of incremental sources is usually associated with relatively small rivers or mainstem tributaries. The method to estimate Ln i , Zmi , and observed riverine TP flux is described in Sections 2.3, 2.4, and 2.6, respectively. The in-stream delivery terms (Ai and Ai 0 ) are assumed to follow first-order reaction kinetics associated with transport distance, as has been widely adopted in previous model applications (Reddy et al., 1999; Alexander et al., 2002); i.e., longer transport distance allows more opportunity for riverine P attenuation. Thus as specified by SPARROW, the in-stream delivery terms Ai and Ai 0 are as follows:

Ai ¼ expð
a$Di Þ

(2)

  Ai 0 ¼ exp
a0 $Di 0

(3)

where Di (m) refers to the length of the mainstem river that hydrologically links the upstream segment to catchment i, and Di 0 (m) refers to the total length of rivers within the incremental catchment i. a (m
1) and a0 (m
1) are the corresponding first-order decay coefficients for TP losses with respect to distance along the mainstem and tributaries, respectively. The latter values were developed using iterative parameter estimation. To perform the simulation on the annual time scale during the 2010e2017 period, we first collected all model inputs for each individual year. We then established geographic relationships by defining stream network topology. Each input dataset for each subcatchment and year must be assigned a unique hydrologic sequence code (i.e., the variable “hydseq” in the SPARROW data file, Schwarz et al. (2006)) to ensure that the model algorithm can correctly identify them in a model run. Model variables were then calibrated to best fit using the observed annual riverine P fluxes from 2010 to 2014 and then validated against using fluxes from 2015 to 2017 (see section 2.6; Table 2S). Note that the model specification as shown in Eq. (1) assumes that steady-state conditions are applicable for each year, so that the watershed properties are indicative of the annual effects of P processes. Our time-varying simulation covers each year from 2010 to 2017, which includes dry years such as 2011, wet years such as 2016, and years with near average precipitation. By combining the years with different hydrological conditions and simulating average behavior, we optimized comprehensiveness of the simulation as we used year-specific data to constrain the model. Model reliability is also expected to be improved, because there are more available data for calibration and validation. In this sense, more parsimonious simulation based on annual time steps with fixed temporal parameters can help locate high-risk hotspots and focus efforts to reduce future P transport from the sources to the downstream lake.

2.3. Phosphorus sources (Ln i ) There are no industrial sites, significant livestock production areas, or municipal sewage treatment plants within the TLW. Pointsource P inputs are thus negligible, which focuses attention on a source inventory that encompasses major non-point-source terms. Such inputs are mainly derived from fertilizer application, atmospheric deposition, and other spatially distributed sources related

to human and livestock waste, which can each be linked with specific land use types. Although these data can be readily obtained from local census documents or published annual statistical summaries (e.g., Statistics (2018)), they are often aggregated at the administrative level and lack the necessary spatial detail required for model inputs, especially when the simulated watershed is small. To establish spatially-explicit P source inputs while maintaining acceptable accuracy, we employ land cover data as surrogate indicators of the contaminant mass supplied by non-point sources. This type of source specification has been widely adopted by SPARROW (Alexander et al., 2004) and other models (Luscz et al., 2015, 2017) including the export coefficient model (Johnes, 1996). TP sources (Ln i in Eq. (1)) in catchment i are represented by the area of four land use typesdtea plantations, croplands, urban areas, and forestland. Each usually contributes significantly different amounts of P load per unit area. Forestland can be a source of TP, because atmospheric P deposition inputs can exceed plant uptake and P loads can be generated through decomposition of organic matter, soil erosion, and chemical weathering of bedrock (Gao et al., 2017). Land cover classifications (Ln i ) were interpreted from SPOT-5 satellite imagery, with 5  5 m resolution. The quality of these interpreted datasets was assured by independent ground verification, with average classification accuracy >92%. Thus, the term, gn Ln i , which consists of the area of specific land use (Ln i ) and the corresponding model estimated export coefficient (gn ), quantifies P loss from the land use in catchment i. To effectively achieve model convergence, the initial export coefficient values were set to values from similar model applications (Alexander et al., 2004; Li et al., 2013). Annual interpreted land cover classifications from 2010 to 2017 were used as model inputs (see Fig. 1S). 2.4. Land-water delivery factors (Z mi ) Landscape variables (i.e., Zmi in Eq. (1)) reflect watershed properties affecting the transport of TP over the landscape and into streams. Many climatic, terrain, and landscape attributes have been demonstrated as land-water delivery factors (see Table 1). These variables were prepared according to standard procedures for the SPARROW model in the GIS platform (see Table 1S for details). Then, a variety of variable sets were assessed and compared by evaluating statistical significance and collinearity of the parameters, model R2, and Bayesian Information Criterion (BIC), an information theorybased measure that incorporates residual sums-of-squares and a penalty related to the number of model parameters (Schwarz, 1978). The final set of model parameters was obtained when an optimal model specification was achieved in terms of model fit, model-estimated coefficients, and residual plots. 2.5. Streamflow data Daily streamflow data were only available at two hydrological stations (see Fig. 1), and so the flows at other ungauged sites were simulated using a HSPF model (Bicknell et al., 1997). To obtain initial parameters for the entire TLW, the simulation was first calibrated and validated in a small catchment without any lakes (the streamflow gauge station in the southern TLW, Fig. 1; see Li et al. (2015) for details). Then, we used a second simulation and parameter optimization for the entire TLW including the lake, with calibration and validation using a gauge at the outlet of the lake. The model simulates the streamflow at the two stations well, with R2 and Nash-Sutcliffe efficiency generally greater than 0.8 for daily streamflow. This allowed us to generate daily streamflow estimates for all 19 water quality station catchments using the calibrated model. These streamflow estimates were used to calculate annual

W. Zhang et al. / Environmental Pollution 255 (2019) 113273

5

Table 1 Parameter selection for the SPARROW models of the TLW. Model components Pollution sources

Delivery factors

Parameters considered Point source Non-point source

Land-to-water delivery factor

In-stream delivery factor

riverine TP fluxes. 2.6. Riverine TP flux Water samples of 19 stream monitoring sites (see Fig. 1) were collected at the end of every month from 2010 to 2017. Raw water samples were used to determine total phosphorus (TP) concentrations. Measurements of TP concentrations in mg L
1 were obtained using the K2SO4 degradation-ammonium molybdate spectrophotometric method (MEP, 2002). There were very few (<1%) missing datapoints; any missing concentrations were estimated by interpolating between the previous and the following month’s values for a station. TP flux is calculated as kg yr
1 from the catchment outlets using monthly TP concentration data and daily streamflow estimates with the LOADEST regression model (Runkel et al., 2004), which has been widely used to estimate annual riverine TP fluxes (Aguilera et al., 2013; Aulenbach, 2013; Zhou et al., 2017). The estimated daily streamflow values (see Section 2.5) are thus paired with monthly TP measurements as inputs to the LOADEST model that estimates annual TP fluxes. LOADEST predictions for the 19 subcatchments yielded acceptable R-squared values (0.64e0.96, with an average of 0.85) and Nash Sutcliffe Efficiency Index (0.61e0.99, with an average of 0.83). 2.7. Land use change scenarios Revealing historical changes in land cover is helpful to design reasonable scenarios for the future. During the short interval between 2010 and 2017, land cover in the TLW underwent rapid change (Fig. 1S). Tea plantations and urban land quickly expanded to cover 10.4% and to 6.9% of the total area as of 2017. This growth came at the expense of forestland and cropland areas, which decreased from 57.7% to 54.0% and from 16.9% to 16.1% of the TRW area respectively, from 2010 to 2017. Additional change matrix analysis of land use transformation confirms that new tea plantations mostly occupy converted forestland, but new urban land mostly displaced croplands. Based on this analysis, we designed land use scenarios, i.e., a 1.5-fold increase in the area of tea plantations, from the current 10%e25%, and an increase in urban land from the current 7%e11%. The magnitudes of these relative increases in area represent a continuation of historic land use change trajectories from 2010 to 2017, in other words, a business as usual scenario. We used MATLAB (The MathWorks Inc., 2008) to conduct a sensitivity analysis with scenarios that anticipate the consequences of future land use change. To simply our analysis, we assume that

Industrial & sewage discharge N inputs from cropland (km2) N inputs from urban land (km2) N inputs from forestland (km2) N inputs from tea plantations (km2) Slope of the terrain ( ) Clay percentage (%) Air temperature ( C) Drainage density (m
1) Pond density (km
2) Total pond area (km2) Runoff depth (mm) Precipitation (mm) River length in incremental catchment (m) River length (m)

Final parameters ✓ ✓ ✓ ✓

✓ ✓ ✓ ✓

land use change would not significantly influence runoff generation, and thus the runoff depth was assumed to be constant at the multi-year average of 2010e2017 values.

3. Results 3.1. Model performance Variables for the TLW simulation were selected for use by evaluating a full array of source variable combinations, including landscape variables and in-stream delivery factors based on statistical significance, model errors, R2, and Bayesian Information Criterion. We followed the standard procedure described by Schwarz et al. (2006) to optimize the model for best overall performance and selected the refined eight key variables (Table 2). Calculated R2 values confirm that 94% of the spatial and temporal variability in annual TP loads was explained on the basis of these selected independent variables, which represent regional processes within the study area (Fig. 2a). On an area basis, the plot of observed and predicted TP yields adjusted for the effects of drainage area on TP flux shows reasonably good agreement (R2 ¼ 0.74, see Fig. 2S). Errors in the predicted TP flux range from
19% to 32%, based on the interquartile range for the differences between the predicted and observed values at the monitoring sites (median ¼ 6.9%). The scatterplot of residuals showed poor linkage with the observed TP loads (R2 ¼ 0.01, Fig. 2b), further suggesting that the model systematically captures the important variations of TP flux and that its specifications are reasonable. The source, land-to-water, and in-stream variables retained in the final model, including their associated statistics, are presented in Table 2. The coefficients with the approximate parametric t-test values display relatively high (p < 0.05) to moderate (p < 0.10) levels of statistical significance for every coefficient. Comparisons of the source-specific parameters (g1  g4 ) allowed us to pinpoint the relative P contributions per unit area from each land use type. Thus P loss rates from urban land, tea plantations, and cropland were estimated to be 9.7, 7.8, and 5.9 times higher, respectively, than such loss from forestland. Land-to-water delivery variables mainly included runoff depth and pond density. This analysis indicates that both are highly significant factors for explaining the spatio-temporal variation of TP loads. The in-stream delivery factors suggest that the mainstem is much more effective per unit of channel length at trapping P than the tributaries.

6

W. Zhang et al. / Environmental Pollution 255 (2019) 113273

Table 2 Evaluation of SPARROW model parameters with confidence intervals (CI) in the Tianmu lake watershed. Variable

Units

Value

Standard Error

Parametric p-value

Lower 95% CI

Upper 95% CI

Forestland (g1 ) Cropland (g2 ) Tea plantations (g3 ) Urban land (g4 ) Runoff depth (d1 ) Pond density (d2 ) a Mainstem decay (a) a Tributary decay (a0 )

kg P$km
2yr
1 kg P$km
2yr
1 kg P$km
2yr
1 kg P$km
2yr
1 mm
1 km
2 m
1 m
1

11.12 65.63 86.52 107.85 0.0011
0.038
5.56E-5
4.16E-5

6.35 33.64 29.88 42.20 7.045E-5 0.017 1.73E-5 1.13E-5

0.082 0.053 <0.05 <0.05 <0.05 <0.05 <0.05 <0.05


1.43
0.85 27.46 24.43 0.0010
0.071
8.99E-5
6.39E-5

23.68 132.12 145.57 191.26 0.0013
0.0045
2.14E-5
1.94E-5

a

Note that mainstem refers to the river reach with a high Horton-Strahler stream-order (3), and tributary mainly refers to low-order streams (<3).

Fig. 2. a) Observed riverine TP loads versus those modeled with SPARROW. b) Residuals plotted against Observed riverine TP loads in 19 sub-catchments of the TLW for 2010e2017 (n ¼ 152). Vertical lines denote 95% confidence intervals.

3.2. Model prediction 3.2.1. Load and yield Spatial distributions of incremental TP loads represent additional P source that originates from the reach’s incremental catchment that are ultimately delivered to Tianmu Lake (Fig. 3a). Distinctions between incremental catchments adjacent to the lake and those at higher elevations to the south are clear, with large P loads attributed to the former and much smaller loads attributed to the latter. Multi-year TP loading into the lake averages 6.0 metric tons yr
1. The Pingqiao sub-basin transports the largest P loads (2.0 tons yr
1) to the lake, followed by the Zhongtian (1.9 tons yr
1), Linhu (1.7 tons yr
1), and Xiasong (0.4 tons yr
1) sub-basins. On an annual basis, total TP loading to the lake in a wet year (2016), reached 15.2 tons, which is almost four times the amount in a dry year (2011). Fig. 3b shows the spatial distribution of incremental yields of TP, which are defined as load divided by the catchment area. Higher yields reflect more intense P loading from a sub-catchment. Multiyear averaged incremental annual TP yields among all subcatchments varied more than 16-fold, from 9.5 to 153.7 kg km
2 yr
1. High annual TP yields were primarily from several western and eastern catchments lying relatively close to the lake. Although there are significant annual variations, the spatio-temporal patterns of incremental P yields are relatively consistent, and thus the hotspot catchments shown in Fig. 3b are good candidate locations for focused management. 3.2.2. Phosphorus source identification Non-point sources of TP in the TLW are of fundamental importance to watershed management and mainly relevant to the four land use types. Their relative contributions to P loss for 2010

through 2017 are shown in Fig. 4a and b. On average, from 2010 to 2017 agricultural land use as cropland and tea plantations accounted for 34% and 27% of exported TP, respectively. Forested areas were also significant contributors of TP loss (21%), but this is related more to the extent of forest coverage than to P yield per unit area of forested land. P yield from urban land was equivalent to 18% of the total and of lesser overall significance. Spatio-temporal patterns are clear, providing important insights for those seeking to reduce P pollution in the TLW (Fig. 4). High P contributions from tea plantations were mainly located in incremental catchments of the Pingqiao and Linhu sub-basins, but those from urban land were characteristic of the Zhongtian sub-basin. In contrast, forests were always the main sources of P loads from the higher elevations in the south, but cropland was the main contributor at lower elevations. The contributions of forestland and cropland to overall P loads in the TLW declined over time (see Fig. 4c). 3.3. Sensitivity to agricultural expansion and urbanization Fig. 5 showed that both agricultural expansion and urbanization would increase P loading into the lake, but the effect is much stronger with new tea plantations. This is apparent from the contour map (Fig, 5a), which shows that P loads are much more sensitive to increases in tea than corresponding increases in urban land. Indeed, if tea plantation area increases from the current 10%e 25% of the total area, total TP loading to the lake is predicted to rise to 9.4 metric tons yr
1, an increase of 49% assuming runoff depth of the 2010e2017 average. If the urban land area is increased at the same time from 7% to 11% of the total area, total TP loading would rise to 9.8 metric tons yr
1. Although this represents an additional 4.3% increase, it is

W. Zhang et al. / Environmental Pollution 255 (2019) 113273

7

Fig. 3. Spatio-temporal distributions of incremental: (a) P load and (b) P yield transported to Tianmu Lake. Each miniature bar plot shows annual loads from each sub-catchment, with colors ranging from light blue (the leftmost bar) in 2010 to dark blue (the rightmost bar) in 2017. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 4. TP sources in each incremental sub-catchment in (a) 2010 and (b) 2017. The percentage change in P contributions from the four land use categories between 2010 and 2017 is shown in (c).

proportionately less significant than the effect of land conversion to tea. This indicates that when land is transformed from cropland to urban land, the reduction of agricultural non-point source pollution can be offset and slightly overridden by enhanced P pollution

associated with urban lands. Such negative impacts can be more directly assessed from predicted TP concentrations in runoff transported into the lake (Fig. 5b and c). If urban land increases from 7% to 11% of the total area, the

8

W. Zhang et al. / Environmental Pollution 255 (2019) 113273

Fig. 5. Dynamic scenarios of continuous change in tea plantations and urban areas and its effects on TP loading in the TLW. (a) A contour map illustrating how estimated TP loading to the lake varies in response to increases in the areas of urban land and tea plantations; (b) estimated changes in flow-weighted TP concentrations as a function of tea land (10%e25%) at different extents of expansion of urban land from 7% to 11% of the total area; (c) estimated changes in flow-weighted TP concentrations as a function of urban land (7%e11%) at different extents of expansion of tea lands from 10% to 25%. Note that flow-weighted TP concentrations were derived from TP loads divided by the total streamflow.

flow-adjusted TP concentration is predicted to increase from the current 0.061 mg L
1 to 0.065 mg L
1. But the flow-adjusted TP concentration is predicted to rise to 0.095 mg L
1 if the tea plantation area is simultaneously elevated from 10% to 25%. It is worth noting that this simulation is still conservative because we assume that the runoff depth remains constant across years. However, future forest-to-tea and cropland-to-urban conversions would enhance runoff, which in turn can promote more P loss. The information underscores that the impacts of land use change on P loads are consequently serious, but also complex and interdependent, requiring more comprehensive consideration. 4. Discussion 4.1. Model evaluation and uncertainty The reliability of our simulation can be assessed both by statistical indicators, including R2, the significance and confidence interval of parameters, and model errors, and also by inferences from the literature. The model’s source-specific coefficients are similar to export coefficients (Johnes, 1996), which in most cases can be obtained from field monitoring (Beaulac and Reckhow, 1982), because both can qualify the amounts of P per unit area that are potentially delivered to streams. Overall, the export coefficients reported for different regions varied widely from 2 to 13 kg km
2 yr
1 for forestland (Johnes, 1996; Wang et al., 2003), 25e80 kg km
2 yr
1 for cropland (Johnes, 1996; Winter and Duthie, 2000), 84e103 kg km
2yr
1 for tea plantations (Li et al., 2013; Nie, 2015), and 50e93 kg km
2yr
1 for urban land (Winter and Duthie, 2000; Li et al., 2013). Modelinterpreted source-specific coefficients generally lie within these range. In addition, relative values of model-interpreted TP export from each land use type in the TLW are also in good agreement with those of observed export coefficients. Plot-scale monitoring confirms that the TP loss from urban land and cropland is highest, contributing more than ten times the P loads per unit area than

forestland (Wang et al., 2003). P export from tea plantations and forestland, although highly variable based on in-situ rainfall-runoff monitoring, established that the export coefficient of tea plantations is ~2e14 times that of the forestland (Nie, 2015). Similar conclusions were also reached from catchment-scale monitoring (Liu et al., 2012; Verheyen et al., 2015). In aggregate, these comparisons indicate the applicability and usefulness of SPARROW for modeling P fluxes in the TLW. It is nevertheless important to note that there are uncertainties associated with our application of the SPARROW model. One of the most important is related to the assumption that steady-state conditions are applicable for each year. This is an oversimplification, because of the occurrence of cross-year or long-term consequences of P processes such as land use legacy effects (Martin et al., 2011; Sharpley et al., 2013; Martin et al., 2016). Although this study emphasizes multi-year averaged model results rather than year-specific results, some annual-based analysis (Fig. 4) is subject to this uncertainty. The parsimonious model structure is also associated with uncertainties. For example, this model cannot explicitly differentiate P forms (such as dissolved and particulate P), and also has not independently parameterized their specific processes. This may result in insufficient interpretation of P transport. Model uncertainties are also associated with variables, estimation methods, and data accuracy. The adoption of annual land covers as surrogate variables of source terms represents one of such uncertainty. Although main TP sources have been indirectly considered by these surrogate variables, other sources such as the weathering of bedrock should be independently included, given that their contribution could be significant in a mountainous watershed. The embedded errors of the estimated riverine P flux derived from the LOADEST regression model can influence model calibration and validation, and also can magnify model uncertainties. Furthermore, we employed annual runoff depth as a proxy variable to characterize land use export in each subcatchment. Runoff depths used here are mainly derived from the HSPF model. Although the accuracy at the two gauge stations is

W. Zhang et al. / Environmental Pollution 255 (2019) 113273

generally acceptable, uneven prediction errors over space and time could amplify the uncertainties of model performance and lead to errors in some individual catchments.

4.2. P transport budgets The process-based nature of the SPARROW model facilitates calculation of the budgets that govern export of P from land and its subsequent transport. This information is valuable, because it forms the basis for efficient, evidence-based efforts to manage P in the basin and minimize its undesirable accumulation in the lake. Fig. 6 shows the multi-year averaged P budgets of the four sub-basins. Annual P export from three of the sub-basins is roughly equivalent, ranging from 3464 to 3856 kg yr
1, but the Xiasong sub-basin (Fig. 1) exports only about one-fifth as much P, partly due to its small catchment size. P retention in the four sub-basins is mainly associated with ponds and rivers, and P loads thus depend on features such as land use type, transport distance, and ponds. The 1025 ponds of the TLW are widely distributed (Fig. 1), thus assuming that all of the exported P is initially subject to pond retention, we calculate that 16e30% of exported P is retained by these waterbodies in the four sub-basins. This is equivalent to the entrapment of approximately 2800 kg P yr
1 and indicates that each pond, on average, can retain ~2.8 kg yr
1. Given that the average pond area in this study watershed is 0.38 ha, the average retention rate per unit pond area is nearly 7.4 kg ha
1yr
1. The estimated P retention rate in ponds of the TLW is

9

comparable with previous findings. Inlet-outlet monitoring of P loads in 38 other ponds showed that their average P retention rate €sche et al., 2010). Kovacic et al. was about 5.7 kg ha
1yr
1 (Kno (2000) investigated TP retention in three small ponds and concluded that TP retention rates were highly varied but typically around 5.0 kg ha
1yr
1. Our model-interpreted retention rate is somewhat greater, and this may relate to the prevalence of tea plantations on steep hillsides above ponds in the TLW. These areas receive intensive fertilizer P applications that can be as high as 324.9 kg ha
1 yr1, which is five times higher than other croplands (Han et al., 2012). In addition, the unique tea cultivation practice of pruning just after harvest abruptly reduces vegetation coverage, resulting in exposure of a large area of bare land at the onset of the rainy season (Li et al., 2013; Nie et al., 2013). Consequently, more eroded soil and associated nutrient is available to be trapped in ponds. The stream network retained a total of 2758 kg yr
1 of P, which accounts for 31% of the total amount of P entering into rivers. This implies that the stream network has an important role in intercepting P inflows, a result consistent with findings from other bamery and Garnier (2007) documented that sins. For example, Ne riverine P retention can be highly significant, as much as 15%e30% of riverine TP loads. House (2003) reported that although riverine P retention was highly dependent on season and the variability of streamflow, ~30% of riverine TP loads can be retained in the river network. In aggregate, these results underscore the crucial buffering roles of in-stream attenuation in TP inflows (Withers and Jarvie, 2008).

Fig. 6. Diagram of simulated P transport budgets in the four sub-basins. Note that the retention rates, which are labeled in red, refer to the percentage changes of transported P masses (blue color). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

10

W. Zhang et al. / Environmental Pollution 255 (2019) 113273

Our results also showed that P retention in the mainstem is much stronger than that in tributaries. As shown in Figs. 6 and 8e31% of P was initially retained in the tributaries within the incremental catchments. The remaining loads in three of the four sub-basins are subject to mainstem retention, which accounts for 11e31% of the transported amounts. As a whole, ~1500 kg yr
1 of P was trapped in the mainstem portion of the TLW, which is 17% higher than that in tributaries (usually with higher topography). Stronger P retention in the mainstem is generally contrary to the dynamics of nitrogen (N), where much weaker removal usually occurs in the mainstem (Alexander et al., 2000; Zhang et al., 2019). This is because in-stream N removal is mainly associated with the benthicehyporheic zone, and the ratio of the benthicehyporheic zone to streamflow declines as stream size increases (Peterson et al., 2001; Zhang et al., 2019). In contrast, P retention in stream networks is mainly associated with sedimentation (Withers and Jarvie, 2008), which is strongly determined by hydrological conditions. When water and associated P are transported to the relatively larger mainstem, sharp reduction in flow velocity allows deposition of the suspended, P-containing sediment load and increased opportunity for soluble P to be absorbed by the soil and taken up by vegetation (Withers and Jarvis, 1998). Notably, and as is most evident in the Linhu sub-basin (Fig. 1), the TLW contains many spatially-independent streams that can directly transport P loads into the lake, with no intervention by mainstem processes. Significant agricultural expansion in these areas would be a particular threat to the lake, given the absence of mainstem mechanisms to reduce P loads. 4.3. Implications for management Water quality of Tianmu Lake is a priority, because it is used for drinking and serves as an attraction for tourists. Our monthly monitoring of TP concentrations in the TLW highlights dual threats: levels at most sites are significantly higher than the TP threshold of eutrophication (~0.02 mg L
1, Correll (1998)), and they are below national standards for drinking (i.e., 0.05 mg L
1, MEP (2002)). P pollution is thus a pressing current issue, and the situation is likely to worsen in the future, because the watersheds of the inflowing rivers face increasing pressure from agricultural expansion, population growth, and urbanization (Diao et al., 2013; Li et al., 2013). Information on spatial-varying sources (Fig. 4) and associated transport (Fig. 3S) provided by the SPARROW model is thus timely as a guide for management, policy, and restoration/rehabilitationdnot just of the TLW, but also of similar headwater watersheds elsewhere. Our scenario analysis underscores the seriousness of agricultural encroachment onto natural landscapes of the TLW as a future source of P pollution (Fig. 5). The current rate of such expansion, which transformed more than 7% of the basin’s forestland into tea plantations in just 7 years (Diao et al., 2013; Zhang et al., 2019), is clearly unsustainable. If this trend continues, key ecosystem services provided by natural vegetation will be significantly impaired, and nutrient inputs and intensified soil erosion that are associated with tea plantations will rise to unacceptable levels (Diao et al., 2013). Model-based spatial patterns of P retention in ponds and during in-stream delivery (Fig. 3S) provide valuable information about areas at high risk of undesirable P loss. These areas should be targeted for preservation of natural ground cover as a strategy to preserve water quality in the lake. Increased P contribution from urban land, especially in the Zhongtian sub-basin (Fig. 4), suggest that the expansion of urban land is a collateral issue warranting attention. Clearly, spatially explicit policies that target priority areas should be put in place to constrain this kind of land expansion. Best Management Practices (BMPs) are effective tools to control

and mitigate pollution from diffuse sources (Yates et al., 2007; Makarewicz et al., 2009), but they can be expensive. This is especially true for engineered structures, such as grass strips, bioswales, detention basins, media filters, retention ponds and wetland basins (http://www.bmpdatabase.org/), which for economic reasons are not often found in the TLW. Model-based identification of hotspots for incremental P yields (Fig. 3b) suggests that these areas are the most logical targets for cost-effective implementation of BMPs and concentrated efforts to raise awareness about BMP approaches via demonstration sites and other pilot programs. Although managing P applications to agricultural lands via precision fertilization and manure management is a cost-effective and non-engineering BMP option (Moore et al., 2011; Robertson and Saad, 2011; Zhou et al., 2017), our surveys of farmer attitudes in the TLW reveal weak environmental awareness, which limits the efficacy of source management. Ignoring environmental damage, local farmers are inclined to apply excessive levels of relatively cheap fertilizers to maximize yields of high margin products such as tea (Diao et al., 2013; Zhang et al., 2019), despite extensive evidence that lower fertilization levels will not reduce yields (Ju et al., 2009; Liu et al., 2016). Demonstration sites in hotspots should also be useful to exhibit the value of this and other non-engineering BMP approaches to farmers. Optimizing interception of P delivery is also important, and so attention should be given to management and maintenance of ponds, which in the TLW remove ~24% of the total P exported from land uses. Storage of such P loads is often only temporary. When the environmental conditions become unfavorable, such as after extreme rainfall events, high P loads can be either released from the sediment or flushed into rivers. Filling of ponds and other actions that interfere with their function should thus be avoided (Sønderup et al., 2016; Chen et al., 2019). Further reduction of TP loads in the stream network can feasibly be enhanced by ecological engineering including river wetland restoration, which also provides habitat for aquatic plants and other creatures (Zedler, 2000; Coveney et al., 2002). Demonstration river wetlands have been constructed in some rivers in TLW (Pueppke et al., 2019), and are proving effective as nutrient traps (Li et al., 2013), but the flashy nature of headwater hydrographs can be destructive for wetland ecosystems (Lane et al., 2018). Additional hydrology mitigation measures, such as construction of overflow dams, will likely be needed. 5. Conclusions The SPARROW model has proved valuable to evaluate P source and transport in a headwater catchment of the upper Taihu Lake Basin that is undergoing rapid agricultural expansion. After calibration to annual TP fluxes at 19 sites over an eight-year period, the model was able to explain 94% of the spatio-temporal variability in P loads and determine the relative contributions of ponds, tributaries, and mainstems to P retention in the catchment’s four subbasins. The model also enabled scenario analysis of the consequences of continued loss of forest land in the TLW and its conversion to tea plantations and other agricultural uses. The model allows several important conclusions to be drawn about the TLW: (i) continued loss of forest land at current rates is unsustainable; (ii) ponds have a disproportionately important role in retaining P in the hilly TLW, and (iii) several model-identified hotspots of P loss are promising candidates to implement engineered and nonengineered BMPs to preserve the water quality and ecosystems in the basin. Acknowledgments This study was financially supported by the National Natural

W. Zhang et al. / Environmental Pollution 255 (2019) 113273

Science Foundation of China (No. 41701040 and 41877513), the Natural Science Foundation of Jiangsu Province of China (No. BK20171100), the National Key Research and Development Program of China (2018YFD1100102) and the Thirteenth Five-Year Plan of the Nanjing Institute of Geography and Limnology (No. NIGLAS2018GH06). The authors wish to express their gratitude to the Asia-Hub Program (No. 2017-AH-10) for funding the collaboration that resulted in this paper. Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.envpol.2019.113273. References , R., Sabater, S., 2012. Linking in-stream nutrient flux to land use Aguilera, R., Marce and inter-annual hydrological variability at the watershed scale. Sci. Total Environ. 440, 72e81. , R., Sabater, S., 2013. Modeling nutrient retention at the waterAguilera, R., Marce shed scale: does small stream research apply to the whole river network? J. Geophys. Res. Biogeosci. 118, 728e740. Ahearn, D.S., Sheibley, R.W., Dahlgren, R.A., Anderson, M., Johnson, J., Tate, K.W., 2005. Land use and land cover influence on water quality in the last freeflowing river draining the western Sierra Nevada, California. J. Hydrol. 313, 234e247. Alexander, R.B., Smith, R.A., Schwarz, G.E., 2000. Effect of stream channel size on the delivery of nitrogen to the Gulf of Mexico. Nature 403, 758e761. Alexander, R.B., Elliott, A.H., Shankar, U., Mcbride, G.B., 2002. Estimating the sources and transport of nutrients in the Waikato river basin, New Zealand. Water Resour. Res. 38, 4-1-4-23. Alexander, R.B., Smith, R.A., Schwarz, G.E., 2004. Estimates of diffuse phosphorus sources in surface waters of the United States using a spatially referenced watershed model. Water Sci. Technol. 49, 1e10. Alexander, R.B., Boyer, E.W., Smith, R.A., Schwarz, G.E., Moore, R.B., 2007. The role of headwater streams in downstream water quality. J. Am. Water Resour. Assoc. 43, 41e59. Aulenbach, B.T., 2013. Improving regression-model-based streamwater constituent load estimates derived from serially correlated data. J. Hydrol. 503, 55e66. Beaulac, M.N., Reckhow, K.H., 1982. An examination of land use-nutrient export relationships. JAWRA J. Am. Water Resour. Assoc. 18, 1013e1024. Bicknell, B.R., Imhoff, J.C., Kittle, J.L., Donigian, A.S., Johanson, R.C., 1997. Hydrological Simulation Program d FORTRAN: User’s Manual for Version 11. U.S. Environmental Protection Agency, National Exposure Research Laboratory, Athens, Ga (EPA/600/R-97/080). Brakebill, J.W., Ator, S.W., Schwarz, G.E., 2010. Sources of suspended-sediment flux in streams of the Chesapeake Bay watershed: a regional application of the SPARROW model. J. Am. Water Resour. Assoc. 46, 757e776. Chen, J., 2007. Rapid urbanization in China: a real challenge to soil protection and food security. Catena 69, 1e15. Chen, W., He, B., Nover, D., Lu, H., Liu, J., Sun, W., Chen, W., 2019. Farm ponds in southern China: challenges and solutions for conserving a neglected wetland ecosystem. Sci. Total Environ. 659, 1322e1334. Conley, D.J., Paerl, H.W., Howarth, R.W., Boesch, D.F., Seitzinger, S.P., Havens, K.E., Lancelot, C., Likens, G.E., 2009. Controlling eutrophication: nitrogen and phosphorus. Science 323, 1014e1015. Correll, D.L., 1998. The role of phosphorus in the eutrophication of receiving waters: a review. J. Environ. Qual. 27, 261e266. Coveney, M.F., Stites, D.L., Lowe, E.F., Battoe, L.E., Conrow, R., 2002. Nutrient removal from eutrophic lake water by wetland filtration. Ecol. Eng. 19, 141e159. Diamantini, E., Lutz, S.R., Mallucci, S., Majone, B., Merz, R., Bellin, A., 2018. Driver detection of water quality trends in three large European river basins. Sci. Total Environ. 612, 49e62. Diao, Y., Han, Y., Zhaofu, L.I., 2013. Spatio-temporal change of tea plantation since 2000 and model-based prediction in the Tianmuhu Reservoir watershed. J. Lake Sci. 25, 799e808. Donnelly, C., Andersson, J.C.M., Arheimer, B., 2015. Using flow signatures and catchment similarities to evaluate the E-HYPE multi-basin model across Europe. Int. Assoc. Sci. Hydrol. Bull. 61, 255e273. Gao, Y., Zhu, G., He, R., Wang, F., 2009. Variation of water quality and trophic state of lake Tianmu, China. Environ. Sci. 30, 673e679. Gao, W., Chen, Y., Liu, Y., Guo, H.C., 2015. Scientometric analysis of phosphorus research in eutrophic lakes. Scientometrics 102, 1951e1964. Gao, Y., Hao, Z., Yang, T., He, N., Wen, X., Yu, G., 2017. Effects of atmospheric reactive phosphorus deposition on phosphorus transport in a subtropical watershed: a Chinese case study. Environ. Pollut. 226, 69e78. García, A.M., Hoos, A.B., Terziotti, S., 2011. A regional modeling framework of phosphorus sources and transport in streams of the southeastern United States. JAWRA J. Am. Water Resour. Assoc. 47, 991e1010. Golden, H.E., Lane, C.R., Amatya, D.M., Bandilla, K.W., Raanan Kiperwas, H.,

11

Knightes, C.D., Ssegane, H., 2014. Hydrologic connectivity between geographically isolated wetlands and surface water systems: a review of select modeling methods. Environ. Model. Softw 53, 190e206. Guo, L., 2007. Doing battle with the green monster of Taihu Lake. Science 317, 11661166. Han, Y., Li, H.P., Nie, X.F., Xu, X.B., 2012. Nitrogen and phosphorus budget of different land use types in hilly area of Lake Taihu upper river basin. J. Lake Sci. 24, 829e837. Hodgkinson, R.A., Withers, P.J.A., 2007. Sourcing, transport and control of phosphorus loss in two English headwater catchments. Soil Use Manag. 23, 92e103. €ck, A., Humborg, C., Gustafsson, B., Hong, B., Swaney, D.P., McCrackin, M., Svanba Yershova, A., Pakhomau, A., 2017. Advances in NANI and NAPI accounting for the Baltic drainage basin: spatial and temporal trends and relationships to watershed TN and TP fluxes. Biogeochemistry 133, 245e261. House, W.A., 2003. Geochemical cycling of phosphorus in rivers. Appl. Geochem. 18, 739e748. Jackson-Blake, L.A., Sample, J.E., Wade, A.J., Helliwell, R.C., Skeffington, R.A., 2017. Are our dynamic water quality models too complex? A comparison of a new parsimonious phosphorus model, SimplyP, and INCA-P. Water Resour. Res. 53, 5382e5399. Jakeman, A.J., Letcher, R.A., Norton, J.P., 2006. Ten iterative steps in development and evaluation of environmental models. Environ. Model. Softw 21, 602e614. Johnes, P.J., 1996. Evaluation and management of the impact of land use change on the nitrogen and phosphorus load delivered to surface waters: the export coefficient modelling approach. J. Hydrol. 183, 323e349. Ju, X.T., Xing, G.X., Chen, X.P., Zhang, S.L., Zhang, L.J., Liu, X.J., Cui, Z.L., Yin, B., Christie, P., Zhu, Z.L., Zhang, F.S., 2009. Reducing environmental risk by improving N management in intensive Chinese agricultural systems. Proc. Natl. Acad. Sci. U. S. A. 106, 3041e3046. Kirchner, J.W., 2006. Getting the right answers for the right reasons: linking measurements, analyses, and models to advance the science of hydrology. Water Resour. Res. 42, 113e115. € sche, R., Schreckenbach, K., Pfeifer, M., Weissenbach, H., Müller, R., Naeve, H., Kno 2010. Balances of phosphorus and nitrogen in carp ponds. Fish. Manag. Ecol. 7, 15e22. Kovacic, D.A., David, M.B., Gentry, L.E., Starks, K.M., Cooke, R.A., 2000. Effectiveness of constructed wetlands in reducing nitrogen and phosphorus export from agricultural tile drainage. J. Environ. Qual. 29, 1262e1274. Lane, C.R., Leibowitz, S.G., Autrey, B.C., Leduc, S.D., Alexander, L.C., 2018. Hydrological, physical, and chemical functions and connectivity of non-floodplain wetlands to downstream waters: a review. JAWRA J. Am. Water Resour. Assoc. 54, 346e371. Laurance, W.F., Sayer, J., Cassman, K.G., 2014. Agricultural expansion and its impacts on tropical nature. Trends Ecol. Evol. 29, 107e116. Li, H.P., Zhu, G., Chen, W., Gao, R., Nie, X., Yu, Z., Diao, Y., Xinyan, L.I., 2013. Current situation of good water quality reservoirs in hilly region of south-east China: protection practices of Tianmuhu Reservoir. J. Lake Sci. 25, 775e784. Li, Z., Liu, H., Luo, C., Li, Y., Li, H., Pan, J., Jiang, X., Zhou, Q., Xiong, Z., 2015. Simulation of runoff and nutrient export from a typical small watershed in China using the Hydrological Simulation ProgrameFortran. Environ. Sci. Pollut. Control Ser. 22, 7954e7966. Liu, Z.a., Yang, J., Yang, Z., Zou, J., 2012. Effects of rainfall and fertilizer types on nitrogen and phosphorus concentrations in surface runoff from subtropical tea fields in Zhejiang, China. Nutrient Cycl. Agroecosyst. 93, 297e307. Liu, X., Sheng, H., Jiang, S., Yuan, Z., Zhang, C., Elser, J.J., 2016. Intensification of phosphorus cycling in China since the 1600s. Proc. Natl. Acad. Sci. U.S.A. 113, 2609e2614. Luscz, E.C., Kendall, A.D., Hyndman, D.W., 2015. High resolution spatially explicit nutrient source models for the Lower Peninsula of Michigan. J. Gt. Lakes Res. 41, 618e629. Luscz, E.C., Kendall, A.D., Hyndman, D.W., 2017. A spatially explicit statistical model to quantify nutrient sources, pathways, and delivery at the regional scale. Biogeochemistry 133, 37e57. Macintosh, K.A., Mayer, B.K., McDowell, R.W., Powers, S.M., Baker, L.A., Boyer, T.H., Rittmann, B.E., 2018. Managing diffuse phosphorus at the source versus at the sink. Environ. Sci. Technol. 52, 11995e12009. Makarewicz, J.C., Lewis, T.W., Bosch, I., Noll, M.R., Herendeen, N., Simon, R.D., Zollweg, J., Vodacek, A., 2009. The impact of agricultural best management practices on downstream systems: soil loss and nutrient chemistry and flux to Conesus Lake, New York, USA. J. Gt. Lakes Res. 35, 23e36. Martin, S.L., Hayes, D.B., Rutledge, D.T., Hyndman, D.W., 2011. The land-use legacy effect: adding temporal context to lake chemistry. Limnol. Oceanogr. 56, 2362e2370. Martin, S.L., Hayes, D.B., Kendall, A.D., Hyndman, D.W., 2016. The land-use legacy effect: towards a mechanistic understanding of time-lagged water quality responses to land use/cover. Sci. Total Environ. 579, 1794e1803. Mayorga, E., Seitzinger, S.P., Harrison, J.A., Dumont, E., Beusen, A.H.W., Bouwman, A.F., Fekete, B.M., Kroeze, C., Van Drecht, G., 2010. Global nutrient export from WaterSheds 2 (NEWS 2): model development and implementation. Environ. Model. Softw 25, 837e853. MEP, 2002. Environmental Quality Standards for Surface Water (GB 3838-2002). China Environmental Science Press, Beijing. Moore, R.B., Johnston, C.M., Smith, R.A., Milstead, B., 2011. Source and delivery of nutrients to receiving waters in the Northeastern and Mid-Atlantic regions of the United States. J. Am. Water Resour. Assoc. 47, 965e990.

12

W. Zhang et al. / Environmental Pollution 255 (2019) 113273

Morales-Marín, L., Wheater, H., Lindenschmidt, K.-E., 2015. Assessing the transport of total phosphorus from a prairie river basin using SPARROW. Hydrol. Process. 29, 4144e4160. Morales-Marín, L., Wheater, H., Lindenschmidt, K.-E., 2018. Potential changes of annual-averaged nutrient export in the south saskatchewan river basin under climate and land-use change scenarios. Water 10 (1438), 1e18. Nasr, A., Bruen, M., Jordan, P., Moles, R., Kiely, G., Byrne, P., 2007. A comparison of SWAT, HSPF and SHETRAN/GOPC for modelling phosphorus export from three catchments in Ireland. Water Res. 41, 1065e1073. mery, J., Garnier, J., 2007. Origin and fate of phosphorus in the Seine watershed Ne (France): agricultural and hydrographic P budgets. J. Geophys. Res. 112, 1e14. Nie, X., 2015. Study on the Process of Nitrogen and Phosphorus Loss and Water Movement at Slope Scale in Hilly Area of Upper Lake Taihu Watershed. Nanjing Institute of Geography and Limnology. University of Chinese Academy of Sciences, N. Nie, X., Hengpeng, L.I., Huang, Q., Diao, Y., Jiang, J., 2013. Characteristics of nitrogen loss via runoff from typical land uses in hilly area of Tianmuhu Reservoir watershed. J. Lake Sci. 25, 827e835. Nie, X.F., Li, H.P., Jiang, J.H., Diao, Y.Q., Li, P.C., 2015. Spatiotemporal variation of riverine nutrients in a typical hilly watershed in southeast China using multivariate statistics tools. J. Mt. Sci. 12, 983e998. Ongley, E.D., Xiaolan, Z., Tao, Y., 2010. Current status of agricultural and rural nonpoint source Pollution assessment in China. Environ. Pollut. 158, 1159e1168. Peterson, B.J., Wollheim, W.M., Mulholland, P.J., Webster, J.R., Meyer, J.L., Tank, J.L., Marti, E., Bowden, W.B., Valett, H.M., Hershey, A.E., 2001. Control of nitrogen export from watersheds by headwater streams. Science 292, 86e90. Pueppke, S., Zhang, W., Li, H., Chen, D., Ou, W., 2019. An Integrative Framework to Control Nutrient Loss: Insights from Two Hilly Basins in China’s Yangtze River Delta. Water (accepted for publication) In press. Radcliffe, D.E., Lin, Z., Risse, L.M., Romeis, J.J., Jackson, C.R., 2009. Modeling phosphorus in the Lake Allatoona watershed using SWAT: I. Developing phosphorus parameter values. J. Environ. Qual. 38, 111. Reddy, K.R., Kadlec, R.H., Flaig, E., Gale, P.M., 1999. Phosphorus retention in streams and wetlands: a review. Crit. Rev. Environ. Sci. Technol. 29, 83e146. Robertson, D.M., Saad, D.A., 2011. Nutrient inputs to the laurentian great lakes by source and watershed estimated using SPARROW watershed models. J. Am. Water Resour. Assoc. 47, 1011e1033. Runkel, R.L., Crawford, C.G., Cohn, T.A., 2004. Load Estimator (LOADEST): A FORTRAN Program for Estimating Constituent Loads in Streams and Rivers. US Department of the Interior, US Geological Survey. Saad, D.A., Schwarz, G.E., Robertson, D.M., Booth, N.L., 2011. A multi-agency nutrient dataset used to estimate loads, improve monitoring design, and calibrate regional nutrient SPARROW models. J. Am. Water Resour. Assoc. 47, 933e949. Said, A., Ross, M., Trout, K., Zhang, J., 2010. Simulation of surface water for UnGauged areas with storage-attenuation wetlands. JAWRA J. Am. Water Resour. Assoc. 43, 546e556. Schindler, D.W., Carpenter, S.R., Chapra, S.C., Hecky, R.E., Orihel, D.M., 2016. Reducing phosphorus to curb lake eutrophication is a success. Environ. Sci. Technol. 50, 8923e8929. Schwarz, G., 1978. Estimating the dimension of a model. Ann. Stat. 6, 15e18. Schwarz, G., Hoos, A., Alexander, R., Smith, R., 2006. The SPARROW Surface WaterQuality Model: Theory, Application and User Documentation. USGS. U.S. Geological Survey, Reston. Seitzinger, S., Harrison, J., Dumont, E., Beusen, A.H., Bouwman, A., 2005. Sources and delivery of carbon, nitrogen, and phosphorus to the coastal zone: an overview of Global Nutrient Export from Watersheds (NEWS) models and their application. Glob. Biogeochem. Cycles 19, 1064e1067. Sharpley, A., Jarvie, H.P., Buda, A., May, L., Spears, B., Kleinman, P., 2013. Phosphorus legacy: overcoming the effects of past management practices to mitigate future water quality impairment. J. Environ. Qual. 42, 1308e1326. Sønderup, M.J., Egemose, S., Hansen, A.S., Grudinina, A., Madsen, M.H., Flindt, M.R.,

2016. Factors affecting retention of nutrients and organic matter in stormwater ponds. Ecohydrology 9, 796e806. Statistics, C.B.o., 2018. Changzhou statistical yearbook. China Statistics Press, Changzhou. Available online: http://tjj.changzhou.gov.cn/html/tjj/2018/ OEJCMFCP_1204/15668.html. (Accessed 9 June 2019). Stone, R., 2011. China aims to turn tide against toxic lake pollution. Science 333, 1210e1211. Sun, X., Zhu, G., Da, W., Yu, M., Yang, W., Zhu, M., Xu, H., Guo, C., Yu, L., Li, H., Li, H., 2018. Thermal stratification and its impacts on water quality in shahe reservoir, Liyang, China. Environ. Sci. 39, 2632e2640. Swaney, D.P., Hong, B., Ti, C., Howarth, R.W., Humborg, C., 2012. Net anthropogenic nitrogen inputs to watersheds and riverine N export to coastal waters: a brief overview. Curr. Opin. Environ. Sustain. 4, 203e211. The MathWorks Inc, 2008. MATLAB Version R2008a (USA). Tilman, D., 1999. Global environmental impacts of agricultural expansion: the need for sustainable and efficient practices. Proc. Natl. Acad. Sci. U.S.A. 96, 5995e6000. Verheyen, D., Van Gaelen, N., Ronchi, B., Batelaan, O., Struyf, E., Govers, G., Merckx, R., Diels, J., 2015. Dissolved phosphorus transport from soil to surface water in catchments with different land use. Ambio 44, 228e240. Wang, X., Wang, X., Wang, Z., Wang, Q., Cai, X., 2003. Nutrient loss from various land-use areas in shixia small watershed of Miyun County, Beijing, China. Chin. J. Geochem. 22, 173e178. Wang, X., Yang, W., Melesse, A.M., 2008. Using hydrologic equivalent wetland concept within SWAT to estimate streamflow in watersheds with numerous wetlands. Trans. ASABE 51, 55e72. Wang, H., Wu, Z., Hu, C., 2015. A comprehensive study of the effect of input data on hydrology and non-point source pollution modeling. Water Resour. Manag. 29, 1505e1521. Winter, J.G., Duthie, H.C., 2000. Export coefficient modeling to assess phosphorus loading in an urban watershed. J. Am. Water Resour. Assoc. 36, 1053e1061. Withers, P.J., Jarvie, H.P., 2008. Delivery and cycling of phosphorus in rivers: a review. Sci. Total Environ. 400, 379e395. Withers, P.J.A., Jarvis, S.C., 1998. Mitigation options for diffuse phosphorus loss to water. Soil Use Manag. 14, 186e192. Yan, W.J., Zhang, S., 2003. The composition and bioavailability of phosphorus transport through the Changjiang (Yangtze) River during the 1998 flood. Biogeochemistry 65, 179e194. Yan, W., Yin, C., Tang, H., 1998. Nutrient retention by multipond systems: mechanisms for the control of nonpoint source pollution. J. Environ. Qual. 27, 1009e1017. Yan, W., Yin, C., Shen, Z., 1999. Nutrient budgets and biogeochemistry in an experimental agricultural watershed in Southeastern China. Biogeochemistry 45, 1e19. Yates, A.G., Bailey, R.C., Schwindt, J.A., 2007. Effectiveness of best management practices in improving stream ecosystem quality. Hydrobiologia 583, 331e344. Zedler, J.B., 2000. Progress in wetland restoration ecology. Trends Ecol. Evol. 15, 402e407. Zhang, W., Swaney, D.P., Hong, B., Howarth, R.W., Han, H., Li, X., 2015. Net anthropogenic phosphorus inputs and riverine phosphorus fluxes in highly populated headwater watersheds in China. Biogeochemistry 126, 269e283. Zhang, C., Brett, M.T., Brattebo, S.K., Welch, E.B., 2018a. How well does the mechanistic water quality model CE-QUAL-W2 represent biogeochemical responses to climatic and hydrologic forcing? Water Resour. Res. 54, 6609e6624. Zhang, W., Li, H., Kendall, A.D., Hyndman, D.W., Diao, Y., Geng, J., Pang, J., 2019. Nitrogen transport and retention in a headwater catchment with dense distributions of lowland ponds. Sci. Total Environ. 683, 37e48. Zhou, P., Huang, J., Hong, H., 2017. Modeling nutrient sources, transport and management strategies in a coastal watershed, Southeast China. Sci. Total Environ. 610e611, 1298.