Water and nutrient balances in a large tile-drained agricultural catchment : a distributed modeling study

Introduction Conclusions References


Introduction
Water, sediment, carbon and nutrient cycles occur over a multiplicity of time and space scales, and govern the dynamics and health of all ecosystems, which are of critical importance to the long-term sustainability of human habitation.Fluxes of water and the variability of water cycle dynamics are key drivers of coupled physical, biogeochemical, ecological and human systems.For example, soil moisture storage is a result of the water cycle processes of rainfall, storage, and movement, which are governed by climatic and landscape features.The amount of nitrate in the soil is a result of human additions at discrete times as well as continuous evolution of biogeochemical processes (transport and reaction), which depend on the magnitude and dynamics of water and carbon cycle processes.Likewise, sediment transport is governed by erosion, sedimentation and re-entrainment processes that are linked to water flow pathways and human activities.Biogeochemical processing and reprocessing occurs as the flow moves along a gradient in the intensity of land use, from urbanized and agricultural lands that are adjacent to a stream bank, through various levels of riparian vegetation and grassy waterways that separate streams from managed landscapes, and to well developed bottomland forest or Published by Copernicus Publications on behalf of the European Geosciences Union.
The interactions and feedbacks between these subsystems that occur at all scales, however, are poorly understood, inadequately observed, and extremely complex.The gaps in our knowledge and understanding of these interacting processes limit our ability to make robust predictions and provide a solid basis for sustainable watershed management.Understanding the interactions between various water and biogeochemical processes is also important in the wider context of climate change and human induced land use and land cover changes, with suggestions that the hydrological cycle may be accelerating as a result.A coupled modeling framework of these subsystems may open new opportunities for studying interacting hydrological and biogeochemical processes, contributing significantly towards improved predictive capability.The move towards such a coupled modeling framework is also motivated by the fact that many of the interacting natural processes cannot be observed directly -instead we are only able to observe spatial and temporal patterns of signatures arising from the process interactions.A pattern dynamical approach that is focused on the identification of internal process interactions on the basis of spatio-temporal patterns of outcomes is an emerging paradigm towards making robust predictions.Such an approach has to be facilitated by a combination of data mining and modeling analysis.
There are quite a few well known models for the coupled hydrological and biogeochemical processes, such as AN-SWERS (Beasley et al., 1977), CREAMS (Knisel, 1980), EPIC (Jones et al., 1984), SWRRB (Williams et al., 1985), HSPF (Donigian et al., 1995), AGNPS (Young et al., 1995), and SWAT (Arnold and Allen, 1996).Some of them, however, are limited to small scale studies (for example, EPIC, CREAMS, AGNPS) and thus not appropriate for large watershed scale applications.Some of them have oversimplified underlying hydrological structures (for example, ANSWERS, SWRRB, HSPF), which prevents deep insight into the role of hydrological processes (especially runoff generation processes).Some of them require intensive parameter calibration, such as SWAT.The coupled hydrological and biogeochemical process model presented here has its foundation in the distributed watershed model, THREW, based on the representative elementary watershed (REW) approach pioneered by Reggiani et al. (1998Reggiani et al. ( , 1999Reggiani et al. ( , 2000)).The THREW model has been designed for macro-scale applications.Its hydrological representation is comprehensive enough for this study, including all key processes contributing to the complex hydrological responses of watersheds.Also, the pasimonious parameter set and scale-consistent spatial organization of THREW decrease the calibration requirement relative to other process based models (Lee et al., 2005(Lee et al., , 2007;;Tian et al., 2006Tian et al., , 2010;;Li et al., 2010).
The work on this paper has been especially motivated by the combination of biophysical (e.g. a plentiful supply of summer rains, and fertile, deep glacial till soils) and social factors (e.g.intensive agricultural advisory services, land use and conservation strategies, and advanced precisionagriculture technologies) that have made the US Mid-West the Nation's breadbasket, albeit with considerable local and remote environmental impacts, such as contributing to eutrophication problems in the Gulf of Mexico.Despite the importance of this region both in terms of agricultural productivity and as a contributor to the environmental problems faced by the Nation, there are still critical knowledge gaps about the complex interactions among the various interacting processes that contribute to local and regional water quality impacts.
In this study we have extended THREW to include the effects of tile drains, which is a major human modification to this agricultural landscape.More importantly, we have extended THREW further to include modules for the interactions between water flow processes and processes associated with the generation of both sediments and nutrients (N and P) after previously published work (Viney and Sivapalan, 1999;Viney et al., 2000).The combined model is then applied to Upper Sangamon River Basin (USRB), a 3600 km 2 tiledrained agricultural catchment located in south-central Illinois, and calibrated on the basis of all available water quality data, including regional summaries.The model is then used to generate insights into the process interactions underlying the observed and model-generated spatio-temporal patterns.
The paper is organized as follows: Sect. 2 describes the distributed computational framework of coupled hydrological and biogeochemical processes at the catchment scale.Section 3 provides the background information on the case study area and data sources.Section 4 lays out the model application results for the water and nutrient modeling, followed by discussion on the hydrological and biogeochemical process interactions.Section 5 closes with the summary.

A spatially distributed hydrological model
THREW is an existing distributed, physically-based hydrological model (Tian et al., 2006(Tian et al., , 2008)), and is built around the representative elementary watershed (REW) concept.Pioneered by Reggiani et al. (1998Reggiani et al. ( , 1999Reggiani et al. ( , 2000)), the REW approach is essentially a thermodynamically consistent framework to derive balance equations directly at the meso-scale for distributed hydrological modeling.The REW in THREW is the smallest resolvable spatial unit of a meso-scale basin which has an explicit spatial boundary, and is the fundamental building block of the model.As shown in Fig. 1 and the accompanying exchanges of mass, momentum, energy etc. that occur within the REW.Although the REW has an explicit and invariant boundary, the boundaries between the sub-zones are mostly varying with time (Lee et al., 2005(Lee et al., , 2007;;Tian, 2006;Tian et al., 2006).In the latest version of THREW the sub-zones are the saturated zone (s-zone), the unsaturated zone (u-zone), the vegetated zone (v-zone), the bare soil zone (b-zone), the snow covered zone (n-zone), the glacier covered zone (g-zone), the sub-stream network (t-zone), and the main channel reach (r-zone), as shown in Fig. 1.To adequately capture the vertical movement of water and nutrient within soil column, the unsaturated zone is further divided into two layers, the upper unsaturated zone (u 1zone) and the lower unsaturated zone (u 2 -zone).The depth of the u -zone is usually fixed (for example, at 0.3 m), and that of u 2 -zone is allowed to vary with the water table.The ensemble of REWs constituting the watershed also interact with each other by way of exchanges of mass, momentum and energy through the inlet and outlet sections of the associated channel reaches.The mass, energy and momentum balances within the individual zones within the REW, and between the REWs, are described using a coupled set of ordinary differential equations (ODE), derived from thermodynamic principles (mass conservation, Newton's laws of motion, 2nd Law of Thermodynamics) by averaging, with a minimum number of simplifying assumptions.These coupled set of ordinary differential equations, together with appropriate closure relations and geometric relations, are the equation set that lies at the heart of the numerical implementation of REW approach.They can be solved using an appropriate numerical algorithm, such as the CVODE solver (please refer to http://www.llnl.gov/casc/sundials/)currently adopted in the THREW model.Details of THREW, including the various (mass and force) balance equations, as well as the details of the constitutive and closure relations, are not presented here for reasons of brevity.These are available in several previous publications (Tian et al., 2006(Tian et al., , 2010;;Mou et al., 2008).
As a distributed hydrological model based on the REW approach, THREW model has significant advantages: (1) it is physically-based, distributed, and of moderate complexity, and thus computationally advantageous; (2) it has a modular framework, in that the various closure relations, i.e., parameterizations of exchange fluxes, can be altered without changing the overall structure and numerical features; (3) because the model formulation ultimately results in a set of balance equations relating to mass, momentum and energy stores (state variables), the coupled set of ODEs are already in state-space form and can be easily adapted for predictions and data assimilation purposes; and (4) compared to grid-based models, the REW-based distributed model will be more suitable for incorporating various types of land use zones, or water use zones, which are typically categorized by zones (urban areas, irrigation districts, etc.).Thus it will allow us to develop spatial connections between REW units (rather than grids) and water use zones.Moreover, THREW simulates the interactions between surface water, soil water

Extension to agricultural basins: tile drainage
Although THREW has been applied to a number of basins in China, US and Europe under various climate and landscape conditions, it has not been applied to an agricultural basin with extensive tile drains, as we have in the US Mid-West.Field studies suggest that tile drainage, where it exists, is usually a very important source of streamflow (Algoazany et al., 2007;Goswami, 2006).It is thus necessary to incorporate the process of tile drainage for successful prediction in these agricultural basins.Tile drainage is an artificial way to remove excess surface and subsurface water from the water-logging land to enable crop growth (Ritzema, 1994).In the mid-west of US, tile drains have been laid out under swamps and wetlands to deplete the soil water in the saturated zone, and to maintain the water table to an acceptable level to facilitate agricultural production.There have been numerous studies on tile drainage, and various modeling approaches have been proposed such as the classical Hooghoudt equation (Hooghoudt, 1940), Kirkham equation (Kirkham, 1958), Ernst equation (Ernst, 1956).Most of these drainage equations are derived based on the Dupuit-Forchheimer assumptions.However, these equations require the exact locations of the tile drains, which are not often available and, moreover, how their effects up-scale to the watershed or REW scale is also not well quantified.Therefore, in this paper we opt for a conceptual description of their drainage effects, in combination of REWscale effective parameters.In fact, the efficiency of tile drains is governed by the subsurface water storage, i.e., the higher the water table is, the faster the saturated soil water is depleted through the tile drains.It is thus not unreasonable to adopt a simple storage-discharge relation to describe the integrated response of all tile drains present at the REW scale.In this work, we adopt the following conceptual relationship to characterize drainage through tile drains at REW scale: where q tile is the rate of saturated soil water being depleted to the channel through the tile drains  properties of the tile drain network.β is an exponent parameter subject to the spatial layout of tile drain system.Equation (1) applies when the focus is on the integrated tile drainage response at large scale, and the detailed information about the tile drain system is not available or is incomplete.

Coupled model of water, sediment and nutrients
The component models for suspended sediments, nitrogen and phosphorus are mostly taken from Viney and Sivapalan (1999) and Viney et al. (2000) with some minor modifications, and only brief summaries are presented here.Note that the processes governing suspended sediments, nitrogen and phosphorus are described at the sub-watershed scale, which makes them consistent with the scale at which hydrological processes are described within THREW.
As shown in Fig. 2, the storages and exchange fluxes of sediments and nutrients are simulated for each of the subzones within a REW, and thus inevitably coupled to the water flow part.Direct interactions between the landscape and atmosphere (e.g., precipitation, fixation of nitrogen by plants, and the volatilization of ammonia) and between the basin and humans (e.g., fertilization and crop harvest) are associated with the v-zone and the b-zone.The vertical movement of nitrogen is coupled with the water movement in the unsaturated zone (u1-zone and u2-zone) and the saturated zone (s-zone).
Hydrol.Earth Syst.Sci., 14, 2259-2275, 2010 www.hydrol-earth-syst-sci.net/14/2259/2010/ The lateral loading of sediments, phosphorus and nitrogen is triggered by surface and subsurface runoff generation and subsequent delivery to river reaches.For instance, the initiation (soil erosion) and routing of suspended sediments on hillslopes are driven by the generation and routing of surface runoff.The fluxes of water and different substances are transported across the watershed through a set of REWs, which are organized around the river network (not shown in this figure).Presentations of more detailed process descriptions for phosphorus, nitrogen and suspended sediments that follow are adapted from Viney and Sivapalan (1999) and Viney et al. (2000).

Sediment model
The sediment model predicts upslope surface erosion and the in-stream processes of deposition, bank and bed erosion, reentrainment and settling.The physics and modeling of sediments are not the focus in this paper, and will not be described in detail here.,The details of the sediment process description are provided in Viney and Sivapalan (1999) and Liu et al. (2009).

Phosphorus model
The phosphorus model describes the processes of precipitation, fertilization, plant uptake, residue decay, sorption, harvest losses, erosion, surface entrainment and subsurface discharge.Most of the phosphorus cycle models proposed in the literature (e.g., Neitsch et al., 2005) separately consider the organic and inorganic stores, which are further subdivided into readily mobilized active pools and slowly changing less accessible stable pools.After Viney et al. (2000), we combine the organic and slowly changing and less accessible stable pools into one single pool, and denote it as particulate phosphorus (PP).The readily-mobilized active pools have been combined into another single pool, denoted as dissolvable phosphorus (DP).Another pool of phosphorus is biological phosphorus.The key components of the phosphorus model are described below.For better understanding of these components and fluxes the readers are referred to Figs. 2 and 6, although the main purpose of Fig. 6 is to show the mass balance of phosphorous.

(i) Phosphorus from rainfall
Precipitation of inorganic phosphorus is assumed to occur at a specified concentration that, for simplicity, is assumed to be constant in time and space.As the surface runoff interacts with the underlying soil, it entrains an amount of soil inorganic phosphorus.The resulting entrained phosphorus augments the concentration of phosphorus already being carried by the surface flow.

(ii) Phosphorus from fertilizer
The rate and timing of fertilizer application is determined by many factors, such as climate conditions, crop plantation, and soil properties and so on.The phosphorus from fertilizer, organic and inorganic, is assumed to contribute to the storage of the top soil layer.

(iii) Leaching of phosphorus
Leaching of dissolvable phosphorus to deeper levels in the unsaturated zone and ultimately to the deep groundwater is neglected by the model because phosphorus anions are much more affiliated to soil particles rather than water molecules.While it is not doubted that phosphorus leaching can lead to significant groundwater pollution according to some standards, its effect on streamflow discharges is considered negligible since the primary sources of phosphorus discharge involve surface and near-surface processes.

(iv) Residue decay
The processes of leaf fall, crop residue accumulation and litter decay are captured by the single term "residue decay".For a crop, a fixed proportion of the biomass phosphorus is assumed to contribute to residue decay after harvesting, and the rate is given by H P should be regarded as a flux averaged throughout the local REW area, [kg/m 2 /s].All the nutrient fluxes and storage items in the rest of this paper, unless specified, are averaged throughout the local REW area, and have the same units which is non-zero during a certain period after harvesting, and zero during the remainder of time.P B is biomass P accumulated during the growing period, [kg/m 2 ].For a forested field, the rate of residue decay is assumed to be the same as the rate of plant uptake.The rest of the biomass phosphorus is harvested and exported out, mainly in the form of grain.

(v) Plant uptake
Plant uptake rate of phosphorus is assumed to depend on the rate of canopy biomass accumulation and therefore varies seasonally.This uptake is extracted from the dissolvable (i.e., labile) phosphorus stores provided that there is sufficient supply, and the rate is given by Plant uptake transfers soluble inorganic P to biomass P. In Eq. ( 3), k UP is a constant coefficient, [kg/m 2 ]. dLAI dt is the rate of increase of leaf area index (LAI), [1/s], and it is assumed that there is no P uptake when LAI decreases or dissolvable phosphorus storage is completely depleted.Fluxing between the dissolvable and organic forms is typically achieved through the complementary processes of mineralization and immobilization, while fluxing between the dissolvable and adsorbed forms is through the processes of desorption and adsorption.Since the organic and adsorbed pools have been combined into a single pool, which we expect to be dominated by the organic component, we could model the net desorption/mineralization flux in term of a simple desorption equation where which is a function of soil type, but in this work a universal value is applied to all soil types for simplicity.It is also assumed that this fluxing does not occur if the soil temperature is below zero degree Celsius (Neitsch et al., 2005, p. 190).Note the net desorption/mineralization flux (from the organic phosphorous store) contributes to the inorganic phosphorus store, while the residue decay (from the biomass phosphorous store) contributes to the organic phosphorous store.

(vii) Phosphorus movement with water flux
Due to its low mobility, soluble phosphorus only moves with surface water flux, including infiltration excess runoff and saturation excess runoff, and the lateral loading rate of DP from hillslope into channel is therefore given by where k SP is a constant coefficient, [1/m], and q s is the lateral water discharge rate (averaged throughout the local REW area) from hillslope into the channel, [m/s].During the transportation of DP through the river network there is no mineralization/immobilization or desorption/adsorption in the channel flow.

(viii) Phosphorus movement with sediment flux
Upslope erosion of organic and adsorbed phosphorus occurs in conjunction with surface sediment erosion and is dependent on the occurrence and presence of surface runoff.Eroded phosphorus is preferentially attached to the finer sediment particles, which in turn tend to be the first eroded.Consequently, the concentration of eroded phosphorus decreases as the mass of eroded material increases.In the absence of quantitative information on the concentration of organic and adsorbed phosphorus in the upper layers, the model assumes an enrichment ratio for upslope erosion as a function of the amount of sediment erosion.The transport of attached nutrients with channel flow is not conservative since the exchange of suspended sediment and channel floor is incorporated.

Nitrogen modeling
The nitrogen model has a similar structure to that of phosphorus.The nitrogen fluxes for plant uptake, harvest/residue decay, surface entrainment and the mobilization and transport of particulate nitrogen are modeled analogously to the corresponding phosphorus fluxes, and will not be repeated here (for more details see Viney et al., 2000).The nitrogen modeling, nonetheless, is more complex for a few reasons.
One is the need to separately predict NO 3 -N and ammonium forms of the dissolvable inorganic component, which necessitates the inclusion of an extra flux, nitrification, to account for nitrogen cycling between these two forms.Secondly, unlike phosphorus, nitrogen undergoes gaseous exchange with the atmosphere, and this exchange has to be modeled explicitly through the processes of ammonium volatilization, denitrification and nitrogen fixation.Furthermore, as NO 3 -N is highly dissolvable, its leaching to deeper levels in the soil profile is a significant loss mechanism, and an explicit modeling of that process is included.For better understanding of these components and fluxes, readers are are referred to Figs. 2 and 7, although the main purpose of Fig. 7 is to present the mass balance of nitrogen.

(i) Atmospheric N fixation
Plant fixation converts atmospheric N (mainly N 2 ) into ammonia, which is directly utilized by numerous prokaryotes in the soil.Therefore it delivers nitrogen from the atmosphere to the ammonium pool, not to the biomass nitrogen store.
The plant fixation rate is modeled as a function of vegetation status.

(ii) Nitrification and volatilization
Nitrification transfers ammonium to nitrate when the soil temperature is higher than a certain value, and the rate is given by Volatilization releases a fraction of ammonium storage as ammonia gas into the atmosphere and is also simulated as a fixed proportion of the ammonium nitrogen pool when the soil temperature is higher than a certain value.
Hydrol The hillslope denitrification process is microbially mediated and occurs primarily in anoxic conditions.In the model, this process is assumed to occur as a fixed proportion of the NO 3 -N pool and occurs only if the soil water content is greater than 90% of the saturated soil moisture content and the soil temperature is higher than a certain value (Williams et al., 1984;Neitsch et al., 2005).
. θ is the soil moisture content.θ s is the saturated soil moisture content.θ c is a threshold soil moisture content.Here θ c /θ s is taken as 0.9 after Williams et al. (1984).

(iv) Nitrogen movement and variation within soil column
Ammonium is easily attracted by negative-charged soil particles, while nitrate is highly mobile.Therefore it is assumed that all nitrate storage is soluble and movable with water.The nitrate storage in the unsaturated soil layer will lose nitrate due to denitrification, plant uptake and leaching, and receive nitrate due to infiltration, nitrification and fertilization.The nitrate storage in saturated soil layer only exchange nitrate with other zones by the way of water flux.

(v) Nitrogen movement with water flux
Nitrate is highly soluble and moves with all types of water fluxes, including infiltration excess runoff, saturation excess runoff and subsurface flow (or tile drainage).The lateral loading of nitrate is simulated similar to that of DP.The transportation of nitrate through the river network is not conservative, i.e., in-stream denitrification is considered.

(vi) In-stream denitrification
While traveling through the river network, NO 3 -N is removed due to in-stream denitrification process.After Donner et al. (2004) and Wollheim et al. (2006), the instantaneous fractional removal ratio is defined as where v f is the apparent nutrient uptake velocity [m/s], H L is the hydraulic load [m/s].In the THREW model, H L is estimated as where h is the water depth [m], τ is the mean residence time [s] given by l is the reach length [m], v is the water velocity [m/s].Note that τ is essentially the mean travel time of NO 3 -N through the main channel zone (r-zone) within each REW.NO 3 -N joins the main channel from mainly two sources: the inflow from upstream channel and lateral loading from the hillslope.
For the NO 3 -N from lateral loading, the mean in-channel travel time is in fact about half of that of the NO 3 -N from upstream inflow.But here it is assumed that the major part of the in-stream NO 3 -N comes from the upstream inflow.This assumption is appropriate for large basins.

(vii) Nitrogen movement with sediment flux
The movement of organic and adsorbed nitrogen with suspended sediment is simulated similarly to PP.

Study area and data
The present modeling study was carried out on the Upper Sangamon River Basin (USRB) in central Illinois, which is representative of the processes and problems associated with agricultural landscapes in the Mid-West region.USRB, with a drainage area of 3600 km 2 , is an agricultural basin with intensive row-crop production.Soil in this basin is dominated by poorly drained silt clay loams and silt loams, which are very fertile due to the high organic content (Demissie and Keefer, 1996).The topography is very flat, with the average slope of the main channel being 0.00049.According to Demissie and Keefer (1996), in 1994, row crops (corn and soybean) covered 85.3% of the whole basin area and grassy crops (small grains and hay) covered 2.4%.Corn and soybean almost equally share the row crop land area.The percentage of area covered by corn is 42.0%, and by soybean is 43.3%, respectively.The biogeochemistry of USRB is altered annually in the spring and fall with widespread yet highly variable applications of nitrogen and phosphorus fertilizers.Current land and watershed management practices, such as dredging of channels, produce rapid transmission of nitrogen and phosphorus from the land surface through soils, riparian areas, and small streams to larger streams and rivers.The extensive production of corn and soybeans, substantial inputs of urban wastewater and agricultural runoff, and modification of the drainage network have altered patterns and rates of nitrogen and phosphorus cycling.The Illinois State Water Survey (ISWS) has conducted a watershed monitoring project for the Lake Decatur watershed, which is a part of USRB (Keefer and Bauer, 2008).They have measured streamflows, and sediment and nutrient concentration at several stations, including Big Ditch and Monticello.As shown in Fig. 3, downstream of Monticello is Lake Decatur which has a significant impact on the movement of water and transport of sediments and nutrients.For the sake of simplicity, in this work we only focus on the drainage area upstream of Monticello.In order www.hydrol-earth-syst-sci.net/14/2259/2010/ Hydrol.Earth Syst.Sci., 14, 2259-2275, 2010 Hourly stream discharge, irregularly sampled concentration values of suspended sediments, NO 3 -N, and dissolved phosphorus were obtained from the long term monitoring project by ISWS (Keefer and Bauer, 2008).Hourly soil temperature data were obtained from the Water and Atmosphere Resources Monitoring Program conducted by ISWS.Potential evaporation time series were extracted from the NOAA/NARR dataset.Vegetation data including LAI were downloaded and extracted from MODIS/terra dataset.Soil properties such as porosity and saturated hydraulic conductivity were extracted from the STATSGO database.The study period is from 1 October 1993 to 30 September 2004, and was chosen according to data availability.
The application of nitrogen and phosphorus fertilizers is an important external input to the catchment, which often exhibits high spatial and temporal variability.Empirical values of fertilization have been obtained from literature and through personal communication (McIsaac and Hu, 2004;G. McIcsaac, personal communication, 2008).For the sake of simplicity, the application of fertilizers is assumed to be spatially uniform and to be carried out twice a year, the first one during 15 March-1 April, and the next during 1 November-15 November (Hu et al., 2007).In most of the areas corn and soybean are planted in rotation.We assume for simplicity that, in each year, 50% of the field area is corn and another 50% is soybean.The harvest of both corn and soybean is assumed to occur in mid-September.

Model application
As shown in Fig. 3, for the implementation of the coupled model, the whole USRB area has been divided into 51 REWs (3600 km 2 ).In this work, nonetheless, the analysis is only focused on the area upstream of Monticello station (1400 km 2 ), which consists of 19 REWs.The coupled model has been run at a 1-h time step.
We divide the whole study period into two parts: a warmup period, 1 October 1993∼30 September 1994, and a calibration period, 1 Ocoteber 1994∼30 September 2004.We use multiple criteria for calibration.For the water part the criteria include optimal Nash-Sutcliffe coefficient (Nash and Sutcliffe, 1970) and the percent bias (defined as the ratio of the difference between simulated and observed runoff volume to the observed runoff volume; Ivanov et al., 2004).Some other signatures of temporal variability are also used during the calibration, such as the regime curve and the flow duration curve, in order to improve the fit of model predictions to observations.For suspended sediments, nitrogen and phosphorus, the calibration has been conducted in order to: (a) satisfy regional mass balances indicated by the empirical data presented in the literature; (b) match the predicted time series to the observed time series as well as possible.
Figure 4 shows simulated and observed streamflow at Monticello station at both the hourly and seasonal scale (i.e., mean monthly streamflows).The results show strong seasonality with two peaks (during winter and spring) and low flows during summer and fall.Comparison between the observed and predicted hydrographs and regime curves suggests that the model captures the variation of streamflow very well at both the hourly and seasonal scale.For the period of 1 Ocoteber 1994∼30 September 2004, the Nash-Sutcliffe efficiency on the basis of hourly flows is 0.67, and the percent bias is 5%.
Figure 5 shows the model predicted time series of NO 3 -N concentrations and dissolved phosphate concentrations (at hourly time step) and the observed time series (at irregular time intervals).We are not presenting the results for suspended sediments, due to lack of data to fine tune model, calibrate model parameters and validate model predictions.The temporal variation of NO 3 -N concentration has been well captured by the model at both Big Ditch and Monticello.It can be inferred that the NO 3 -N loads (product of water discharge and NO 3 -N concentration) has also been satisfactorily reproduced.On top of this, one might notice that the NO 3 -N concentration at Monticello appears to be lower than that at Big Ditch.This decrease of NO 3 -N concentration from upstream to downstream may most likely be due to in-stream Hydrol.Earth Syst.Sci., 14,[2259][2260][2261][2262][2263][2264][2265][2266][2267][2268][2269][2270][2271][2272][2273][2274][2275]2010 www.hydrol-earth-syst-sci.net/14/2259/2010/   denitrification process, which will be discussed later.As for dissolved phosphorus, the model captures the temporal variation at Big Ditch, but significantly underestimates the concentration of dissolved phosphorus at Monticello, especially in the summer and fall seasons.A possible explanation for this under-estimation is the effluent discharge from the urban areas between Big Ditch and Monticello, including the towns of Mahomet and Monticello.Effluent from the local sewer system and wastewater treatment plants is discharged into the Sangamon River, which introduces non-negligible amounts of nutrients into the river, especially phosphorus.Dissolved phosphorus from effluent discharge, in the form of point-source pollution could make a significant contribution to the in-stream concentrations of phosphorus in the summer and fall seasons.Due to lack of reliable observation data, nutrient inputs through effluent discharges are not included in the current version of the model.This might have led to the poor prediction of dissolved phosphorous from the model.The amount of nitrogen such as nitrate from effluent discharge is rather small compared to the other sources contributing to the channel, so its impact on the nitrate concentration is insignificant.
As mentioned before, model calibration involved not only comparisons of model predicted against observed time series within the USRB, but also checks of broad measures of water and nutrient balances (regional space scale and annual time scale) against published estimates from Illinois region, to ensure that model predictions are consistent.Tables 1 and 2 present a comparison of various aspects of regional nitrogen and phosphorus balances between model predictions within USRB and regional estimates obtained from the literature (McIsaac and Hu, 2004;Hu et al., 2007;David and Gentry, 2000;Howarth et al., 1996;Gentry et al., 2009), demonstrating reasonable consistency in both N and P predictions.
Upon completion of model calibration (as in the above), model simulations were performed to generate an annual average and catchment-wide picture of the fate of both nitrogen and phosphorus.The results are presented in Figs. 6 (for phosphorus) and 7 (for nitrogen).In the case of P, the main input is fertilizer (30 kg P/ha/yr) and the main output is annual harvest (of the crops) which takes out almost 29.6 kg P/ha/yr, with relatively small amounts exported to rivers in dissolved form (0.3 kg P/ha/yr) and in particulate form (0.5 kg P/ha/yr).There is of course considerable    internal processing (plant uptake, generation of plant residue and mineralization), which are included in the model in conceptual form (see Viney et al., 2000 for details).The picture is very different and more complex in the case of N, where in addition to fertilizer application (95 kg N/ha/yr) in the form of ammonia, there is in addition large amount of fixation by plants (65 kg N/ha/yr), and small amount of precipitation (10 kg N/ha/yr).The resulting total inputs (170 kg N/ha/yr) is partitioned into removal through harvest (116 kg N/ha/yr), release into atmosphere in gaseous form (16 kg N/ha/yr), and the removal through runoff in dissolved form (32 kg N/ha/yr) and particulate form (5 kg N/ha/yr).The biggest component (more than 90%) of the runoff export is through tile drainage.Just as in the case of P, there is considerable internal processing, including the conversion of organic nitrogen (as in plant residue) to ammonia through mineralization, from ammonia to nitrate through nitrification and from nitrate into nitrogen through denitrification, as well as plant uptake and generation of plant residue.These processes are of course included in the model in conceptual form (see Viney et al., 2000 for details).Knowledge of these relative estimates is extremely useful for targeting future research towards understanding and quantifying key components of the annual nutrient balances, and associated process controls.

Multi-scale interactions between water and nutrient cycling processes
In spite of the average water and nutrient balances presented in Figs. 6 and 7, there is considerable temporal (and spatial) variability in the nutrient mass balances, which are intimately related to climatic and hence hydrological variability at multiple time scales.The coupled model predictions are next used to throw light on these interactions, and the resulting temporal patterns.
Figure 8 shows the relationship between the annual runoff and annual mass balance of nitrogen and phosphorus at the basin scale.Annual runoff depth is a hydrological indicator and is itself a result of the interactions between variability in climatic forcing and landscape properties.Roughly, the wetter the climate (due to more precipitation or less evaporative energy) is, the larger the annual runoff depth.Therefore the annual runoff depth can be regarded as a first order indicator of the inter-annual variability of wet/dry conditions, recognizing that some of the inter-annual variability of runoff could be caused by variability in intra-annual variability of climate forcing.In Fig. 8, annual balance of nitrogen and phosphorus is expressed in terms of total annual mass brought into the basin, total annual mass exported out of the basin, and annual storage change within the basin.The results presented in Fig. 8 show that total nitrogen inputs, dominated by fertilizer and plant fixation, do not show a significant relationship with annual runoff.Although annual precipitation clearly impacts annual runoff, the concentration of nitrogen in the precipitation is small, so the annual mass of deposition through precipitation is negligible compared to the corresponding amounts of fertilizer application and plant fixation.Fertilizer application is human related, and is assumed constant in this study.Plant fixation is a function of nutrient storage and the growing status of the crops, and does lead to significant inter-annual variability of the annual nitrogen inputs.But this inter-annual variability of nitrogen inputs is much less than that of nitrogen outputs, and for environmental reasons, our focus is thus on the latter.Total nitrogen output, including river loading (export) of nitrogen, field denitrification and volatilization, in-stream denitrification and grain export (through harvest), show an increasing trend with increasing annual runoff depth.Correspondingly, this contributes to a systematic decrease of nitrogen storage with increase of annual runoff depth, from a positive change (storage supplement) during dry years to a negative change (storage depletion) during wet years.Inter-annual variability of phosphorous mass balance, on the other hand, is similar to that of nitrogen, but the variations of the output, and thus the storage, are much smaller compared with the magnitude of annual phosphorous input (i.e., compare the units of the vertical axes in Fig. 8).
In order to gain more insights into predicted behavior between the annual nitrogen output and annual runoff depth, and the mediating role of nitrogen storage/depletion, the annual variations of various components of the nitrogen output are plotted against annual runoff depth, as shown in Fig. 9. Firstly, the results show that in the case of both N and P, grain export is the largest component of the annual export (as was already pointed out in Figs. 6 and 7).The model results in Fig. 9 show that in the case of N, grain export is slightly decreasing with annual runoff, whereas non-grain export increases significantly with annual runoff.In the case of P the changes with annual runoff depth are quite small and negligible.Note that grain export is a significant portion of annual  accumulated biomass gain (from plant uptake), and plant uptake itself is subject to many factors such as soil moisture, soil temperature, crop growing status and nitrogen storage in the soil.
The bottom panel of Fig. 9 presents the breakdown of the non-grain part of the nutrient export into its various components.In the case of N, the biggest component is riverine dissolved export, which increases strongly with increase of annual runoff.The other three major components, i.e., field denitrification and volatilization, riverine denitrification and particulate riverine export are smaller, relative to the riverine dissolved export, but also appear not to be dependent on annual runoff.One can therefore see the connection between the increased dissolved nitrate export and depletion of nitrate storage during wet years, and decreased nitrate export and accumulation of nitrate storage in dry years.The net result of this is that average annual concentrations of dissolved nitrate in rivers in this region can remain constant between years, a type of chemostatic behavior that is being widely reported (Darracq et al., 2008;Godsey et al., 2009;Basu et al., 2010).On the other hand, while the results for P show a strong dependence on annual runoff, the magnitudes are so low that one cannot draw definitive conclusions.
The interaction between hydrological and biochemical processes is manifested not only in the inter-annual variability, but also in the intra-annual variability.For example, Fig. 10 shows the monthly variation of nitrogen storage and streamflow.Nitrogen storage variation is subject to both the input and output.From Fig. 10 one can see that the nitrogen storage peaks twice a year due to fertilizer application, and is depleted significantly in the month of September due to harvesting and during winter and spring when the highest amount runoff is produced.Among the output components, harvesting and riverine export are relatively significant and play an important role in the depletion of nitrogen storage.For phosphorus, the inputs are dominated by fertilizer application, and the outputs are almost completely dominated by grain export.Riverine export of DP and PP do not appear to have any significant impact on phosphorus storage variations.
Further insights into the role of the interactions between hydrological and biochemical processes on nutrient export, as shown in Figs. 9 and 10, can be gained by exploring the relative effects or contributions of different runoff generation components.Figure 11 shows the breakdown of three components of runoff generation within USRB watershed, and the fractions of NO 3 -N lateral loading (from the hillslope into the channel) carried by these different runoff components.Figure 11a shows that tile drainage is the most important runoff component right through the year, and takes up about 80% of the annual runoff generated within USRB.This is consistent with the field observations in neighboring regions with similarly intensive tile drain systems and similar soils and topography (Algoazany et al., 2007;Goswami, 2006).Dunne (saturation excess) overland flow and subsurface stormflow in the catchment constitute relatively small fractions of total runoff, whereas Hortonian runoff (infiltration excess) is virtually negligible.Figure 11b shows the corresponding breakdown of the total nitrate export into components carried by the three different runoff generation mechanisms.The results show that tile drains carry even a larger fraction of the nitrate removed in dissolved form by runoff.Particulate nitrogen is mainly carried by surface runoff, along with the sediment flux.Once the nutrients are delivered to the nearest river reach, they are then transported down the stream network.Figure 12 shows the riverine export of nitrogen, showing the dissolved component is the dominant component, whereas riverine export of particulate nitrogen (the part carried by the suspended sediment) is rather small, since it is carried mainly by the Dunne overland flow (which is small).Note that the seasonal variation of riverine export of NO 3 -N is in phase with the seasonality of streamflow (especially tile drain flows).
The riverine flux of NO 3 -N, before being exported out of the basin, is subject to in-stream denitrification, which is usually considered a significant loss (Alexander et al., 2009;David and Gentry, 2000;Howarth et al., 1996).In USRB, the obvious decrease of NO 3 -N concentration from Big Ditch (upstream) to Monticello (downstream), as shown in Fig. 5, is an indicator of this process.Our model study shows that, without incorporation of in-stream denitrification process this decrease of NO 3 -N from upstream to downstream cannot be reproduced.The rate of in-stream denitrification is controlled by many hydrological and biogeochemical factors, such as channel water depth, channel flow velocity and nitrate concentration.Nitrate concentration affects in-stream denitrification by the way of uptake velocity, i.e., uptake velocity decreases with the increase of nitrate concentration (Mulholland et al., 2008).In our model constant uptake velocity is assumed, so the effect of nitrate concentration is not incorporated explicitly.We thus focus on the impacts of channel discharge on in-stream denitrification of NO 3 -N, as shown in Fig. 13.According to Eqs. ( 9)-( 11), the rate of in-stream denitrification increases with the channel length and decreases with the channel water depth and flow velocity.Figure 13 shows a significant seasonality of in-stream denitrification efficiency.The denitrification efficiency is defined here as the percentage of in-stream flux removed by in-stream denitrification per unit channel area (channel area = local channel length × channel width).It is highest in August when the channel water depth and flow velocity are smallest, and lowest in May when the channel water depth and flow velocity are largest.As for the spatial variability of in-stream denitrification, it is more significant in headwater channels than in downstream channels.Besides Big Ditch and Monticello stations, we add another location, Shively, located between Big Ditch and Monticello, in order to better present spatial variability of in-stream denitrification.Model results suggest that the most dominant factor for the predicted spatial variability of denitrification appears to be the local channel water depth.The water depth in headwater channels, which have small drainage area contributing runoff, is much lower than that in downstream channels, which have large drainage areas contributing runoff into them.Channel water depth is also tightly related to flow velocity, i.e., usually the latter increases with the former giving fixed channel geometry.Therefore the in-stream denitrification efficiency is significantly higher in the channel near Big Ditch than those near Shively and Monticello.Channel length is another factor affecting in-stream denitrification efficiency.The local channel length corresponding to Big Ditch is 18.3 km, to Shively is 38.4 km and to Monticello is 11.4 km (estimated from DEM).In general, the longer the channel length, the longer the residence time of nitrate within the channel and therefore the higher the in-stream denitrification efficiency.In this case, however, the impact of channel length is apparently smaller than that of channel water depth in controlling the spatial variability of in-stream denitrification efficiency.
Hydrol.Earth Syst.Sci., 14, 2259Sci., 14, -2275Sci., 14, , 2010 www.hydrol-earth-syst-sci.net/14/2259/2010/ 5 Summary and conclusions In this paper we have explored the coupled water and nutrient balances in a large tile-drained agricultural watershed in central Illinois, with the use of a distributed model based on the representative elementary watershed (REW) approach.We compared average annual estimates of the various components of the runoff generation against two previous experimental studies, confirming that about 80% of the streamflow in the basin is carried by tile drain flows.Likewise, average annual estimates of the various components of the nutrient (N and P) balances were compared against estimates obtained from several previous experimental studies in the literature, and found good agreement.Once again, tile drains are found to be the carrier of over 90% of the riverine export of dissolved nutrients, especially nitrate.In the case of P, over 98% of the fertilizer application is removed through grain harvest, and only a small fraction (less than 2%) is exported with runoff either in dissolved or particulate form.In the case of N, however, nitrogen fixation by plants represents 40% of the total annual inputs to the catchment (fixation + fertilization), of which slightly over 20% is exported with runoff mostly in dissolved form, predominantly by tile drain flow.The remainder is removed through grain harvest.
The coupled model was also used to gain insights into the interactions between hydrological and biogeochemical processes, and the role of climate and consequent hydrologic variability on nutrient export processes.The results showed that there is a very dependence on the strength of annual runoff and the annual export of nutrients, especially dissolved nitrate component.Assuming that nutrients inputs through fertilizer application is constant between years, and the observation that removal by grain harvest decreases only slightly with increase annual runoff, it is found that relatively dry years are characterized by nutrient accumulation in soil and relatively wet years are characterized by nutrient removal from soil storage.The net result of higher runoff and higher nutrient runoff in wet years and vice versa means that annual average nutrient concentration can be expected to stay relatively constant in such human-impacted agricultural regions.This phenomenon may be one of the causes of chemostatic behavior that has been reported in some agricultural regions of the world.This is not the case for phosphorus removal, however, since in this case the removal of phosphorus by runoff is minor comparing with the removal by harvesting.
This work has demonstrated that a parsimonious model of coupled water, sediment and nutrient balances can be developed that does justice to much of the multi-scale variability of hydrological and biogeochemical processes and their interactions, which are essential for the simulation and prediction of sediments and nutrients in large agricultural catchments.The model presented here can serve as a numerical framework, not only for making predictions of the effects of climate and land use changes, but also to provide guidelines for undertaking new observations and new process studies that are critical for improving the predictive capability of such models in the future.Still, improvements are needed in several areas, including the transportation of phosphorous by tile drainage, an explicit treatment of nutrient uptake by vegetation (including varieties of food and biofuel crops and natural vegetation), and denitrification processes within the river network, including a more accurate representation of channel hydraulic geometry.Continuous measurements of nutrient concentrations in tile drains, river reaches at a range of scales and in the hillslopes are needed to improve process descriptions in the model and to validate the model predictions.This is left for future research.

Fig. 1 .
Fig. 1.Spatial delineation in THREW model.(a) A basin is divided into a number of representative elementary watersheds (REW).(b) Each REW is further divided into several sub-zones.

Fig. 2 .
Fig. 2.Conceptual illustration of coupled water, sediment and nutrient modeling in THREW.P I represents inorganic dissolvable phosphorous.P O represents organic phosphorous and soil-absorbed inorganic phosphorous.N O is organic nitrogen.N NO3 is nitrate.N NH4 is ammonium.

Fig. 4 .
Fig. 4. Comparison of the model predicted and observed runoff response at Monticello.The regime curves are normalized by the total upstream drainage area of Monticello.

Figure 5 .Figure 5 .Figure 5 .Figure 5 .
Figure 5.Comparison between the predicted and observed nitrate and phosphate 1 concentration series.The simulated NO3-N and DP concentration series are at hourly scale; 2 while the observed series are at irregular intervals, most biweekly.There is no observation 3 some time.4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 42 Figure 5.Comparison between the predicted and observed nitrate and phosphate 1 concentration series.The simulated NO3-N and DP concentration series are at hourly scale; 2 while the observed series are at irregular intervals, most biweekly.There is no observation 3 some time.4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 42 Figure 5.Comparison between the predicted and observed nitrate and phosphate 1 concentration series.The simulated NO3-N and DP concentration series are at hourly scale; 2 while the observed series are at irregular intervals, most biweekly.There is no observation 3 some time.4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 42 Figure 5.Comparison between the predicted and observed nitrate and phosphate 1 concentration series.The simulated NO3-N and DP concentration series are at hourly scale; 2 while the observed series are at irregular intervals, most biweekly.There is no observation 3 some time.4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

Fig. 5 .
Fig. 5. Comparison between the predicted and observed nitrate and phosphate concentration series.The simulated NO 3 -N and DP concentration series are at hourly scale; while the observed series are at irregular intervals, most biweekly.There is no observation some time.

Fig. 6 .
Fig. 6.Simulated phosphorous cycling (all values are in Kg P/ha/yr, averaged through the drainage area of Monticello station).

Figure 9 . 3 Figure 9 . 3 Figure 9 . 3 Figure 9 . 3 Fig. 9 .Figure 10 Fig. 10 .
Figure 9. Interactions between hydrological and biochemical processes at the annual scale 1 2 3 Figure 9. Interactions between hydrological and biochemical processes at the annual scale 1 2 3 Figure 9. Interactions between hydrological and biochemical processes at the annual scale 1 2 3 Figure 9. Interactions between hydrological and biochemical processes at the annual scale 1 2 3 Fig. 9. Interactions between hydrological and biochemical processes at the annual scale.

47Figure 11 . 4 Fig. 11 .
Figure 11.Seasonal variation of runoff components and the loading of NO3-N by different 1 runoff components.All values are averaged through the upstream area of Monticello. 2 3 4

Figure 12 .
Figure 12.Seasonal variation of riverine export of nitrogen.

Figure 13 .Fig. 12 .
Figure 13.Seasonal variation of channel waterdepth and in-stream removal efficiency (for the local channel reach corresponding to each station).The in-stream removal efficiency is defined as the percentage of in-stream flux removed per unit area of channel by in-stream denitrification, estimated as NO3-N in-stream removal / (upstream inflow + lateral hillslope inflow) / channel area