Capturing Interactions between Nitrogen and Hydrological Cycles under Historical Climate and Land Use: Susquehanna Watershed Analysis with the Gfdl Land Model Lm3-tan

16 We developed a process model LM3-TAN to assess the combined effects of direct human 17 influences and climate change on Terrestrial and Aquatic Nitrogen (TAN) cycling. The model 18 was developed by expanding NOAA's Geophysical Fluid Dynamics Laboratory land model 19 LM3V-N of coupled terrestrial carbon and nitrogen (C-N) cycling and including new N 20 cycling processes and inputs such as a soil denitrification, point N sources to streams (i.e. 21 sewage), and stream transport and microbial processes. Because the model integrates 22 ecological, hydrological, and biogeochemical processes, it captures key controls of transport 23 and fate of N in the vegetation-soil-river system in a comprehensive and consistent 24 framework which is responsive to climatic variations and land use changes. We applied the 25 model at 1/8 degree resolution for a study of the Susquehanna River basin. We simulated with 26 LM3-TAN stream dissolved organic-N, ammonium-N, and nitrate-N loads throughout the 27 river network, and we evaluated the modeled loads for 1986-2005 using data from 16 28 2 monitoring stations as well as a reported budget for the entire basin. By accounting for inter-1 annual hydrologic variability, the model was able to capture inter-annual variations of stream 2 N loadings. While the model was calibrated with the stream N loads only at the last 3 downstream Susquehanna River Basin Commission station Marietta (40.02' N, 76.32' W), it 4 captured the N loads well at multiple locations within the basin with different climate 5 regimes, land use types, and associated N sources and transformations in the sub-basins. 6 Furthermore, the calculated and previously reported N budgets agreed well at the level of the 7 whole Susquehanna watershed. Here we illustrate how point and non-point N sources 8 contributing to the various ecosystems are stored, lost, and exported via the river. Local 9 analysis for 6 sub-basins showed combined effects of land use and climate on soil 10 denitrification rates, with the highest rates in the Lower Susquehanna Sub-basin (extensive 11 agriculture; Atlantic coastal climate) and the lowest rates in the West Branch Susquehanna 12 Sub-basin (mostly forest; Great Lakes and Midwest climate). In the re-growing secondary 13 forests, most of the N from non-point sources was stored in the vegetation and soil, but in the 14 agricultural lands most N inputs were removed by soil denitrification indicating that 15 anthropogenic N applications could drive substantial increase of emission, an 16 …


Introduction
Biologically available nitrogen (N) in terrestrial ecosystems has significantly increased via anthropogenic nutrient inputs: artificial fertilizer, cultivation of N-fixing crops, and fossil fuel consumption (Galloway et al., 2004(Galloway et al., , 2008)).This increase has caused acidification and N saturation in some terrestrial and aquatic ecosystems (Henriksen and Brakke, 1988;Kelly et al., 1990;Murdoch and Stoddard, 1992;Howarth, 2002).N-saturated soils and streams are also major sources of nitrous oxide (N 2 O) emissions, which is a potent Published by Copernicus Publications on behalf of the European Geosciences Union.
greenhouse gas (Albritton et al., 1994).Other concerns include severe water-quality problems associated with cultural eutrophication, which results in harmful algal blooms and hypoxia in rivers, lakes, estuaries, and coastal zone ecosystems (Smith, 2003;Smith et al., 2006).Climate change and variability also affect water quality through the distribution of high and low flow extremes (Scavia et al., 2002;Howarth et al., 2006).It is generally accepted that microbial processes related to the N cycle are strongly influenced by abiotic factors, and a warm or wet climate provides favorable environments for certain groups of bacterial activities.Quantification and management of the diverse and coupled effects of human activity and climate change on N cycling requires a comprehensive model of the relevant coupled processes that can support the design of optimal nutrient loading controls to maintain desirable water quality and terrestrial ecosystem integrity.
To characterize implications of human-and climate-driven perturbations in the earth's N cycling and its implication for water and air quality, the next generation of N cycling models need to (1) account for regional and local changes in terrestrial and aquatic ecosystem structure and functioning, (2) represent in a consistent manner emissions and transformation of N to air, rivers, and coasts, and (3) be global in extent and integrated with climate and earth system models.Previously, none of the existing models addressed the above three challenges.Here we present a novel modeling framework capable of addressing these challenges and, prior to its global application, we evaluate this modeling framework in the Susquehanna River Basin whose sub-basins vary in climate, land use, and associated N sources and transformations, with a detailed data set of observations.
There has been keen interest and progress in modeling the N cycle in terrestrial ecosystems.However, in most models vegetation and land-use type distribution are prescribed and do not change with time.Modeling studies with EPIC, AN-IMO, and CENTURY/DAYCENT typically prescribe crop distribution and simulate crop production and related nutrient and carbon (C) cycling (Sharpley and Williams, 1990;Parton et al., 1993;Williams et al., 1995;Kroes and Roelsma, 1998;Del Grosso et al., 2009).Because these models do not simulate decadal-to-century changes in vegetation structure (e.g., forest regrowth after harvesting), they are likely to overlook changes in the storage of N in vegetation.Furthermore, during wood harvesting and forest clearing for agriculture, biomass residue is an important additional input to the soil organic C-N pools.Such additional N inputs lead to additional N inorganic loads.In addition, many regional models (e.g., EPIC, ANIMO), which have been applied to far smaller basins compared to the Susquehanna watershed, often use basin-specific parameters for mineralization, nitrification, and denitrification, which complicates their globalscale application for studies on the decadal-to-century scale.LM3-TAN is capable of describing N dynamics with a universal parameter set -the same parameters for all of the sub-basins within an area of 71 220 km 2 and time periods for this simulation.LM3-TAN is among very few modeling frameworks (e.g., CLM-CN: Thornton et al 2007; CLM4MOD: Thomas et al., 2013) that can be used as a component of an Earth System Model -that is, it is capable of representing sub-diurnal exchanges of moisture, energy, and C and N species within the land-atmosphere interface.Unlike CLM4MOD, LM3-TAN simulates water quality in the rivers and nutrient loadings to the coastal environment.
Contrary to the simulations of land models limited to the terrestrial component, most watershed models do estimate stream N concentrations and loads, but they simplify or neglect many key mechanisms describing terrestrial N dynamics (e.g., vegetation and land-use dynamics, interactive C-N feedbacks on vegetation and soil microbial processes; phenological leaf drop and its contribution to soil organic matter pools).INCA-N and SWAT are widely used geographic information system (GIS)-based watershed models (Wade et al., 2002;Schilling and Wolter, 2009).However, when it comes to large-scale applications, because these models are semi-distributed, they are less capable of representing spatial variability, requiring users to define the number and sizes of sub-basins, in which land use and all of the processes for each land use are assumed to be homogeneous and needed to be defined individually.This limits their ability to analyze complex land-use management scenarios.In this class of models, RHESSys is one of a few models with an ecology component that can be used to investigate interactions between ecosystems and hydrological processes according to climate variability (Tague and Band, 2004;Beckers et al., 2009).However, like most models, because these models do not simulate vegetation and land-use-type distribution, a specific parameter set that describes typical soil, vegetation, and land-use characteristics has to be developed using its special module when a study site requires different vegetation or soil types from its default application.This explains why RHESSys has only been applied to very small or subsections of catchments (Band et al., 2001;Tague and Band, 2004).
Given the current lack of models that link terrestrial C-N cycling, long-term vegetation, and land-use dynamics to N loads and concentrations in streams, accounting for different N species, the goal of this research was to build a model to simulate stream N loads that is based on a global-scale terrestrial and N-enabled land model, followed by its testing on a large and complex watershed, for which many years of stream discharge and stream N data are available.For this purpose and to assess the combined effects of direct human influences and climate change on terrestrial and aquatic nitrogen (TAN) cycling, we developed a process model LM3-TAN.The new features include integrated effects of point and non-point sources on river N loads, a soil denitrification module, and stream microbial processes.
We applied LM3-TAN to the Susquehanna River Basin, the largest of the watersheds in the northeastern US, draining an area of 71 220 square kilometers, at the resolution of 1/8 • .The model was evaluated using 20 year (1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) data of stream ammonium (NH + 4 ) and dissolved organic N loads as well as stream nitrate (NO − 3 ) N loads from 16 monitoring stations.For each six sub-basins, we conducted local analysis to assess combined effects of land use and climate on the soil denitrification.We then built up an N budget and compared it with the corresponding reported budget to better understand how point and non-point N sources contributing to the various ecosystems are stored, lost, and exported via the river at the level of the whole Susquehanna watershed.Although there are several parameters that required calibration by fitting simulated to reported stream N loads, these parameters are used universally for the entire basin where climate, soil, vegetation, and land-use characteristics vary.Efforts have been made in the development of this model to limit the number of calibrated parameters.

Overview
LM3-TAN is an expansion of earlier Geophysical Fluid Dynamics Laboratory (GFDL) land models, beginning with LM3V of Shevliakova et al. (2009), which describes vegetation and C dynamics.LM3-TAN was expanded to include vegetation-and soil-N dynamics from LM3V-N (Gerber et al., 2010), new soil physics and hydrology from LM3 (Milly et al, 2014), and N cycling processes described here.LM3 was used as a component of the GFDL Earth System Models (Dunne et al, 2012) and included several enhancements, such as vertically resolved soil physics and hydrology and explicit river dynamics and physics.LM3-TAN includes soil denitrification and transport and chemistry of N cycle in rivers.This version of the model allows more complete tracking of N through the soil-river continuum.In this section, we first summarize key features of the model, and then we describe the newest N cycling features.
LM3V simulates distribution of five vegetation functional types (C3 and C4 grasses, and temperate-deciduous, tropical, and cold evergreen trees) on the basis of total biomass and prevailing climate conditions.The model tracks hundreds of years of land-use change using global land-use transition scenarios that were historically reconstructed by combining satellite-based contemporary patterns of agriculture with historical data on agriculture and population (Hurtt et al., 2006).The four land-use types are natural vegetation (land undisturbed by human activities), secondary vegetation (land formerly disturbed by human activities), cropland, and pasture.The model is spatially distributed, and each grid cell consists of up to 15 tiles: 1 natural vegetation, 1 cropland, 1 pasture, and 1 to 12 secondary vegetation tiles representing unique disturbance histories (i.e., de/reforestation, agricultural practice change).Exchanges of water, energy, and between land and atmosphere are computed with a timestep of 30 min.At-mospheric and terrestrial reservoirs include C pools in vegetation (leaves, fine roots, sapwood, heartwood, and labile C storage), soil (fast and slow), and anthropogenic storage.The C pools in the vegetation are updated on a daily timestep to account for vegetation growth and allocation, leaf drop and display, and natural mortality and fire.The soil C, which is supplied by the vegetation both naturally and during landuse conversion, is stored in two pools with different turnover times.

Coupled C-N dynamics in vegetation and soil
The previous two soil C pools in LM3V were divided into four pools (fast and slow litter, and slow and passive soil organic matter) in LM3V-N.Each C pool in the vegetation and soil was paired with a respective N compartment using poolspecific C : N ratios.The decomposition processes release biologically available forms of N (NO − 3 -N; NH + 4 -N).This allows the simulation of N limitation on plant growth and biological N fixation as well as N feedbacks on organic matter decomposition and stabilization.Inorganic N is removed by sorption to soil particles, plant uptake, immobilization into long-lived organic compounds, and hydrological leaching, while organic N is lost through fire, hydrological leaching, and mineralization.Loss of nitrate N by soil denitrification was not differentiated from the hydrological nitrate-N leaching in LM3V-N.

Improved soil and river physics and hydrology
LM3 introduced vertically distributed soil-water, soil-ice and temperature profiles extending many meters below the surface, but with high resolution (thinnest layer 0.02 m) near the surface.Water (potentially) discharges laterally from each soil layer to the local river reach.Each horizontal grid cell of the model contains only one river reach, and each reach discharges to another reach in the downstream grid cell, following a network that ultimately discharges to the ocean; the sub-grid-scale stream network is ignored.Relations among discharge, storage, velocity, width, and depth in each reach are specified according to Leopold and Maddock (1953).

Synthesis and extension of earlier developments
For this study, we first combined the lumped N model LM3V-N with the distributed physics of LM3.To complete the N mass balance, we next added a soil denitrification module.Finally, we added stream transport and microbial processes to track the fate of soil N leaching and resolve N dynamics in the aquatic ecosystem.Each of these steps is described below.Figure 1 shows stores and fluxes of N in the resultant model, along with relevant processes.Newly introduced or adjusted parameters from the earlier developments are summarized in Table 1 and variables are listed in Table 2.

Merging lumped N model with distributed physical model
To account for dependence of processes in the lumped soil C and N pools upon the vertically resolved physical states of the soil (temperature and water content), the latter were vertically averaged with an exponentially decaying weight function of depth (e-folding depth of 10 m).Leaching of any mobile constituent was defined as the product of a concentration and the sum of lateral and vertical discharge from the soil layer between the surface and a depth of 10 m.The concentration of available N was calculated as dividing available N contents by the effective soil depth, which was approximated assuming C weight content 3.4 % and average soil density 1500 kg m −3 .The available N refers to the N contents reduced by buffering factors which represent processes such as sorption to soil particles.To compensate for many processes that were not accounted for in the model, calibration factors for each N species were introduced to slow down overall N movement from the soil to the stream.These factors include impacts of soil microbes, which are able to take up and in-corporate all N forms (NO − 3 -N, NH + 4 -N, DON) with a much greater capacity than plant uptake (Nordin et al., 2004).The nitrate calibration factor also accounts for storage in groundwater since nitrate (the primary form of N in groundwater) can persist for decades at high levels with increasing N applications.This is further explained by Bachman et al. (1998) which reported that 17-80 % of the N delivered to streams of the Chesapeake Bay watershed was through groundwater.Furthermore, the lumped single-layer N sub-model bypasses most of the vertically distributed hydrologic system, and the soil N leaching based on the average water drainage is transferred directly from the N layer into the stream.These calibration factors were fit to match interannual variations of reported and simulated stream N loads to make up for this modeling approach as well as the unresolved processes that might cause interannual stream N loads to be more sensitive to climate variability than those in reality.Considering its importance in groundwater, a relatively larger size of the nitrate N factor is expected.The need to incorporate these calibration factors, which are at the present basin specific, indicates that future improvements to LM3-TAN should focus on resolving River equations these processes (i.e., N cycle in microbes, reservoirs, and vertically distributed soil layers).Dissolved organic, ammonium, and nitrate N leaching from the soil are described as:

Denitrification in soil
Denitrification is a process that reduces nitrate or nitrite to gaseous forms (e.g., NO, N 2 O, N 2 ) in anaerobic conditions, where the oxidized N species serve as a terminal electron acceptor in metabolism by soil-denitrifying bacteria.The rate of denitrification generally depends on soil nitrate content or concentration, soil water content (a surrogate for oxygen content), and soil temperature.Because soil nitrate contents are relatively low and limiting under natural conditions, we used a first-order loss function with respect to soil nitrate N content, with adjustments for the influence of soil water content and temperature to simulate soil denitrification rate: where D N is the soil denitrification rate (kg m −2 year −1 ); f S is a soil water content reduction function; f T is a soil temperature reduction function; k denitr is a first-order denitrification coefficient (1 year −1 ); where T is the soil temperature ( • C); T r is a reference temperature ( • C); T p is a parameter; Q 10 is a factor change in rate with a 10 • C change in temperature; S is the soil water content; S t is a threshold soil water content; S max is the maximum soil water content; S min is the minimum soil water content; w is an empirical constant.Heinen (2006) tabulates reported values of the various parameters introduced above.Figure 2 shows the effects of the reduction functions on the soil denitrification rate that were applied in diverse models as well as LM3-TAN.Figure 2a shows how fast soil nitrate N content is reduced to half of the initial amount depending on the different first-order denitrification coefficients.As temperature increases the bacterial activities increase exponentially (Fig. 2b).Soil denitrification occurs and increases nonlinearly only if soil water content exceeds a certain threshold point due to enhanced anaerobic bacterial activity (Fig. 2c).The soil water content reduction function for other microbial processes (e.g., mineralization, nitrification) used in LM3-TAN is also shown in Fig. 2d.Because k denitr is by far the most widely used parameter of these, with reported values ranging over 3 orders of magnitude, our strategy was to fix the other parameters using reported values and to calibrate the model by determining k denitr within the bounds reported in the literature.Because soil denitrification and nitrate-N leaching are competing sinks of nitrate N in the soil, soil denitrification increases as soil nitrate-N leaching or stream nitrate-N load decreases; thus k denitr was fit to match reported and simulated stream nitrate-N loads.
The wide ranges of the functions discussed above are mostly driven by the dependencies of the parameters on specific regions (with different soil properties, vegetation, land use, etc.).Given a number of proposed individual functions, it seems that there is no universal process module to simulate soil denitrification.Because such reduction functions display a diversity of shapes as ecosystems are modeled over a range of climate patterns, vegetation type, and landuse practices, soil denitrification on a large scale cannot be modeled without proper adjustments that compensate for the site-specific properties.This explains why only a few studies have applied models to watersheds larger than 1000 km 2 despite the diversity of existing dynamic N models and why semi-distributed models often parameterize these individual functions for each of the sub-basins in large-scale applications.We hypothesize that LM3-TAN's integrated modeling framework, which is capable of simulating long-term vegetation functional type and land-use change as a function of changes in CO 2 , climate, and human influences, allows us to use a universal parameter set to simulate soil denitrification for each of the distinct sub-basins.Still, care has to be taken when applying the model to other watersheds that may be very different in terms of soil and climate properties from the Susquehanna watershed.Furthermore, because soil denitrification becomes zero-order in extreme nitrate-rich environment, instead of using the first-order loss function for all of the land-use types, using a Monod function for agricultural land use may help LM3-TAN's global application where N loadings would vary widely.

Microbial processes in rivers
Despite its importance to water quality, processes that control N removal from water bodies are rarely resolved in watershed-scale models, due to both uncertainties in measurement techniques and lack of measurements.To date, none of studies focusing on river denitrification rate is based on measurements of an entire river network, but rather only on the data from low-order streams or individual catchments.
Here we applied a nonlinear regression function based on the LINX (Lotic Intersite Nitrogen experiment; Mulholland et al., 2008) reach-scale measurements that correlates river denitrification rate with nitrate-N concentration and river depth to estimate the reaction rate constant of river denitrification for each reach (Alexander et al., 2009).River denitrification happens mainly in the benthic and/or hyporheic zones.
Therefore, a river denitrification rate that is inversely proportional to the river depth accounts for the ratio of water column to benthic area.The measured reaction rate constants vary from 0.034 to 117 (1 day −1 ), and we chose the median value 0.53 (1 day −1 ) as the minimum reaction rate constant of river denitrification.Equation ( 12) indicates that the reaction rate constant decreases with an increase in nitrate-N concentration and river depth since both b 1 and b 2 are negative, and it increases with temperature.Reaction rate constants for river mineralization and nitrification were calibrated to match stream N loads.
Figure 1 shows structure of the river component.Each reach directly receives N from point sources (e.g., sewage and waste-water discharge) and indirectly receives N from non-point sources (e.g., atmospheric deposition, fertilizer, manure, and legume applications) via soil leaching.The N loads in a reach are routed downstream with the water as following.
where i is DON, NH + 4 , or NO − 3 ; R i is the river N (kg m −2 ); F in i and F out i are the inflow and outflow of the river N (kg m −2 s −1 ); L i is the N leaching from the soil (kg m −2 s −1 ); P i is the N point source (kg m −2 s −1 ); f T is the stream temperature reduction function; T is the water temperature (C); T r is the reference water temperature ( • C); T p is a parameter; k min , k nitr , and k denitr are the reaction rate constants for river mineralization, nitrification, and denitrification (1 s −1 ); k denitr, min is the minimum reaction rate constant of river denitrification (1 s −1 ); C NO − 3 is the nitrate N concentration (µmol N L −1 ); H is the river depth (m); b 0 , b 1 , and b 2 are the constants; c t is the log re-transform bias correction factor; C d, s is a unit-conversion constant (day s −1 ).

Study site
The Susquehanna River Basin, where nearly 4 million people live, is the largest of the watersheds in the northeastern US and drains an area of 71 220 square kilometers, contributing two-thirds of the annual N load to the Chesapeake Bay (Fig. 3).The basin includes 2293 lakes, reservoirs, and ponds (322 km 2 ) as well as 50 190 km of rivers and streams.The main stem of the Susquehanna River originates at Otsego Lake, New York (NY), and flows about 750 kilometers through NY, Pennsylvania (PA), and Maryland (MD) to the Chesapeake Bay at Havre de Grace, MD.The Susquehanna Large River Assessment Project reported that only 6.9 % of water-quality values exceeded their standards, but the majority of these exceedances were for nutrients (e.g., TN, TP) (Hoffman, 2009), explaining why the Chesapeake Bay suffers from nutrient enrichment problems and hypoxia.
The reported year 2000 land use was about 63 % forest or wooded, 19 % cropland, 7 % pasture, 9 % urban, and 2 % water.The Upper Susquehanna River flows through mostly forested and agricultural land, with some small communities and one larger population center, then confluences with the Chemung River at Sayeare, PA.The West Branch Susquehanna Sub-Basin is mostly woods and grasslands.The Middle Susquehanna River, from the confluences with the Chemung River at Sayeare, PA, to the confluences with the West Branch Susquehanna River at Sunbury, PA, flows along very diverse land use.The Lower Susquehanna Sub-Basin contains extensive agriculture and several large population centers.The other major urban areas are found within the Juniata Sub-Basin (Hoffman, 2008).
The geology of the watershed is mainly clastic sedimentary rock of sandstone and shale.Elevations vary from 30 meters at the Chesapeake Bay in Maryland to 955 meters in central New York State (McGonigal, 2011).The Great Lakes and Midwest climate exert influence over the Upper Susquehanna, Chemung, and West Branch Susquehanna Sub-Basin, whereas the Atlantic coastal climate affect on the other portions of the watershed.The basin has experienced severe droughts about once every decade, and the worst droughts occurred in 1930, 1939 and 1964.The basin is also one of the most flood-prone watersheds in the nation with frequent and localized flash floods every year.The worst recorded flooding in the basin happened in 1972 as a result of Tropical Storm Agnes.

Stream sampling description
Stream discharge data are provided by the network of stream gauges operated by the US Geological Survey (USGS), which collects and summarizes time series data to derive annual, monthly, and daily stream discharge and statistics (Fig. 3).Chemical constituents of the basin's water were monitored by the USGS and Susquehanna River Basin Commission (SRBC).One USGS and six SRBC long-term nutrient monitoring sites monitored since 1985 and nine newly introduced SRBC sites monitored since either 2004 or 2005 to the present (Table 3;Fig. 3;McGonigal, 2011;USGS, 2014) were chosen for model evaluation.The 16 sites vary in sub-basin area and land use.Among the USGS and SRBC sites, the Conowingo and Marietta sites on the main channel of the Susquehanna River have the largest sub-basin areas, respectively, 70 189 and 67 314 km 2 .The sub-basin of the Conestoga site contains extensive agriculture (48 %) and the most populated urban land use with several large population centers (24 %) within a very small area (1217 km 2 ).The West Branch River flows mostly along woods and grasslands to the Lewisburg site.The long-term sites collect two samples per month.Additional samplings are made during seasonal storm conditions.The collected water samples are analyzed for various N species: dissolved N (DN), dissolved nitrite and nitrate (DNO 23 ), dissolved ammonia (DNH 3 ), dissolved organic N (DON), and dissolved ammonia and organic N (DKN) in milligrams L −1 .In addition, annual, seasonal, and monthly loads are computed by the minimum variance unbiased estimator (ESTI-MATOR; SRBC, 2006; USGS, 2014).River temperatures were reported when the samplings were collected for the chemical analysis of stream waters.

Anthropogenic N sources
Anthropogenic N data over two decades  were provided by the Chesapeake Community Modeling Program (CCMP).Atmospheric deposition data were provided by the county-based land segments.Fertilizer, manure, and legume applications as well as combined sewer overflows (CSOs) were provided by the land-river segments of the GIS-based Phase 5.3 Community Watershed Model (USEPA, 2010a).The atmospheric deposition data were calculated by the Chesapeake Bay Program (CBP) Airshed Model, which is a combination of a regression model of wet deposition (Grimm and Lynch, 2005) and the Community Multiscale Air Quality Model (CMAQ) that estimates dry deposition (Dennis et al., 2007;Hameedi et al., 2007).The fertilizer, manure, and legume data were estimated for the years of 1985, 1987, 1992, 1997, 2002, and 2005  Over the two decades, the total N sources decreased by about 20 %.The atmospheric deposition was predominantly nitrate-N, accounting for about 69 %; ammonium-N 27 %; organic-N 4 %.The sum of the fertilizer, manure, and legume applications consisted of 49 % ammonium-N, followed by 37 % organic-N, and 14 % nitrate-N.In particular, the ammonium-N and organic-N loads had considerable variability across the spatial domain because they were strongly influenced by local emissions from the extensive agricultural areas.
Figure 4 shows spatial distribution maps of the applied anthropogenic N sources, which were calculated as a spatial resolution of 0.125 • by 0.125 • and a temporal resolution of 1 year.For each grid cell, which consists of up to 15 land-use tiles, atmospheric depositions (nitrate-N, ammonium-N, and organic-N) were applied to all of the land tiles, and fertilizer, manure, and legume applications (nitrate-N, ammonium-N, and organic-N) were applied only to the cropland tiles.Combined sewer overflows (nitrate-N, ammonium-N, and organic-N) were directly applied to the river reaches.The 20 year (1986-2005) average non-point and point N sources for the six sub-basins are summarized in Table 4.The thick solid arrows in Fig. 1 depict fluxes of each of N species for the anthropogenic N sources to the corresponding terrestrial and river pools, respectively.

Model forcing and simulations
The model was implemented with a spatial resolution of 0.125 • by 0.125 • with time increments of 30 min.The model was forced using reported hydrological data cycled over a horizon of 61 years  to perform long-term simulations.The data include precipitation, specific humidity, air temperature, surface pressure, wind speed, and short-and long-wave downward radiation with a spatial resolution of 1 • by 1 • on timescales of 3 h (Sheffield et al., 2006).Land-use change was simulated from 1704 to 2005 using a scenario of land-use transitions (Hurtt et al., 2006).Preindustrial CO 2 concentration assumed as 286 ppm was applied from 1704-1799, and changes in CO 2 concentrations were applied from 1800-2005 using reported data from NOAA's Earth System Research Laboratory.For 250 years , the estimated preindustrial N deposition (Dentener and Crutzen, 1944;Green et al., 2004;Gerber et al., 2010) was applied as a uniform annual rate.We then applied the 1985 reported anthropogenic N data from 1954 to 1984, and reported annual anthropogenic N data from 1985-2005.

Evaluation of stream waters and N loads
We simulated with LM3-TAN stream dissolved organic-N, ammonium-N, and nitrate-N loads throughout the river network.The model was calibrated by comparing the modeled stream N loads with the corresponding reported N loads at the downstream SRBC station Marietta, in which contributions of the entire watershed to the stream flows and N loads can be assessed.Thus, temporal evaluation of the stream discharges and N loads for the period 1987-2005 was focused on at the Marietta station.River data from the 16 monitoring stations (1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) were also used to evaluate spatial stream discharges and N loads.
Using global hydrological data and a universal parameter set for the entire watershed, the model produced reasonable temporal patterns of annual stream discharge.The simulated stream discharges were in good agreement with the reported values in dry years and periods (July to September), but the model underestimated stream discharges in wet years and periods (March to May).Overall, although the 19 year average simulated discharge was about 28 % lower than the corresponding reported value, their linear and rank correlations were significantly high (Table 5), implying that the bias was systemic and accounted for in the calibration of the N species.
Due to their complex physical and biogeochemical interactions with soil particles and soil organic matter, simulat-ing reactive transport of ammonium and dissolved organic N is far more challenging than simulating nitrate N transport.For example, the correlation at Marietta between stream discharge and nitrate N load (R 2 = 0.98) was significantly higher than that for dissolved organic N (R 2 = 0.48) or for ammonium N (R 2 = 0.85) loads, implying that in addition to the hydrological processes governing soil N transport to rivers, terrestrial physical and microbial processes (e.g., sorption to soil particles, organic matter decomposition and stabilization) have to be accounted for when estimating stream ammonium and dissolved organic N loads.This, plus the fact that the highest component in the overall stream N load is nitrate N, explains why existing watershed models have focused on stream nitrate N loads, and neglected ammonium and dissolved organic N loads.Within the LM3-TAN's integrated modeling framework, we estimated all of the N species for the entire drainage network.
Simulated and reported dissolved N loads were graphed in different units: millions of kg year −1 and kg km −2 year −1 (normalized by its sub-basin area summarized in Table 3).Among the six long-term monitoring sites, the highest and lowest amount of river N loads were reported and simulated at the Marietta and Conestoga sites, respectively (Fig. 6a).This finding is consistent with the general view that the amount of stream N loads is proportional to the size of the basin area.A very high N flux was reported at the Conestoga site (Fig. 6b), which can be explained by extensive agriculture and urban land use in its sub-basin.Because the West Branch Susquehanna is dominated mostly by woods and grasslands, the Lewisburg site had the lowest N flux.The model also captured the stream N loads at the 15 monitoring sites well (Fig. 6c and d).These results attest to the model's ability to correctly simulate the stream N loads for the entire basin based on the climate as well as land use and the corresponding N sources and transformations in the sub-basins.

Spatial distribution of stream N load and soil denitrification rate
Observation of the spatial distribution of the river N load (Fig. 7) and soil denitrification rate (Fig. 8d) helps to identify the extent of the terrestrial and aquatic N pollution across the basin.A large amount of N is exported via the main stem of the Susquehanna River as well as its three major tributaries, where many small-order streams converge.The N loads in the streams increase gradually from the headwaters to the watershed outlet, implying that the N loads to the rivers exceed N removal mechanisms within the rivers.Although stream N loads are in general higher in the larger rivers, at the Lower Susquehanna Sub-Basin, high N loads are present even in small-order streams due to extensive agricultural land use.4.An analysis of the six sub-basins shows that the combined effects of land use and climate on the soil denitrification rate, which were the highest in the Lower Susquehanna Sub-Basin (extensive agriculture; Atlantic coastal climate) and the lowest in the West Branch Susquehanna Sub-Basin (mostly forest; Great Lakes and Midwest climate).These results show that the most significant soil denitrification is associated with extensive agricultural land use (non-point sources).The calculated R 2 statistic between the monthly soil denitrification rate and soil water content (R 2 = 0.51) was significantly higher than that for soil temperature or soil nitrate N content, implying that the soil water content played the greatest role in the soil denitrification process among the three factors.This is because the soil denitrification occurred and increased nonlinearly only when the soil water content exceeded the threshold point (S t = 0.577).The significant effect of the soil water content on the soil denitrification is further illustrated in the upper-east side of the Upper Susquehanna Sub-Basin, where extremely low soil water content (Fig. 8a) impeded the overall soil denitrification process (Fig. 8d).Comparison between the calculated and reported budgets of N sources, retention, lost, transport, and river export at the level of the whole Susquehanna watershed for the period 1988-1992 (Boyer et al., 2002;Breemen et al., 2002;Seitzinger et al., 2002;USEPA, 2010a).

N budget
As a further means of evaluating the model output, we compared the simulated N budget for the period 1988-1992 to the budget constructed by Boyer et al. (2002), Seitzinger et al. (2002), andVan Breemen et al. (2002) for the same period (Fig. 9).Overall, reasonable agreements were found be-tween these two budgets.Total N inputs to the whole basin were reported as 4774 kg km −2 year −1 (atmospheric deposition + fertilizer + forest and agricultural N fixation + net N import in feed and food), while we applied an N of 4443 kg km −2 year −1 (atmospheric deposition + fertilizer + manure + legume + sewage) using the data sources provided by CCMP (USEPA 2010a).The simulated soil denitrification (−4 %), harvest rates (+7 %), river export (−1 %), and river denitrification (−5 %) agreed well with the corresponding reported values.To investigate the importance of N removal within rivers, we ran an experiment in which the reaction rate constant for river denitrification was set to zero.We then compared N loads within the rivers with and without river denitrification.Figure 10 shows a spatial map of the difference in N loads between these simulations, which represents the river N removal.A large amount of N was removed along the main stem of the Susquehanna River as well as its three major tributaries, implying that the N removal increases gradually as distance from the headwaters increases.About 28 % of the N that enters the rivers was removed by river denitrification.
For the entire basin, we divided the simulated land use into either agricultural land (cropland and pasture) or secondary forest (land formerly disturbed by human activities).We then graphed simplified N budgets for each land use (Fig. 9c).The reported agricultural land use was 29 % (Fig. 9a), whereas the model simulated 24 % of cropland and pasture (Fig. 9b).In the secondary forest land, most of the applied N (43 %) was stored in the terrestrial system (vegetation and soil pools), whereas the highest proportion of the applied N was removed by soil denitrification (44 %) in the agricultural land.These results imply that applications of artificial N to agricultural lands can result in considerable soil denitrification rates, and thus significant increase of N 2 O production.This is evident when comparing maps of the applied fertilizer, manure, and legume N applications (Fig. 4c and  d) and the simulated soil denitrification (Fig. 8d) that corresponds well, especially in the Lower Susquehanna Sub-Basin with extensive agricultural land use.Even if there are some discrepancies between these two budgets, we can conclude that the reactive transport of N from the terrestrial to aquatic ecosystems was appropriately simulated by the model, providing suitable descriptive information for the entire drainage network.

Conclusions
Results of our study show that LM3-TAN captures well the key mechanisms that control N dynamics in the climateplant-soil-river system.Specifically, we demonstrate: -On a sub-basin scale with different climate and land-use regimes, the LM3-TAN properly simulates terrestrial N cycling, including effects of long-term vegetation dynamics, land-use changes, and hydrological cycles.The interaction among those three processes allow LM3-TAN to capture soil C-N organic matter and mineral N transformations as well as soil emissions of nitrate-N and leaching of dissolved organic, ammonium, and nitrate N.
-The ability to capture N soil budget and losses then enables LM3-TAN to consistently characterize trends and variability in riverine N inputs and exports of ammonium, dissolved organic, and nitrate N with explicit representation of their transformations and transport in rivers.
-In the re-growing secondary forests, a large fraction of the N from atmospheric deposition has been stored in the vegetation and soil, but in the agricultural lands most N inputs were removed by soil denitrification, indicating that anthropogenic N inputs could drive substantial increase of N 2 O emission, an intermediate of the denitrification process.
-LM3-TAN captures the effects of long-term trends and variability of hydrological cycles (e.g., precipitation, soil water content, stream discharge) on N cycling in vegetation-soil-river system, and thus resolves interan-nual variations of stream N loadings caused by climate variability.
-The model results suggest that the soil denitrification is most sensitive to soil water variations.
-Among the six sub-basins, the soil denitrification rate was the highest in the Lower Susquehanna Sub-Basin with the most intensive land-use non-point N sources as well as with the warmest and wettest soils, attributed to the Atlantic coastal climate.
-Even though the N denitrification and riverine biogeochemistry N modules were calibrated only at the last downstream station Marietta, application of the universal parameters over the entire watershed produced simulations which compared well with other observational stations.The applicability of the universal parameters in other watersheds is a subject of future research.
-This study shows that linking terrestrial N and C cycling, long-term land-use and vegetation dynamics, and hydro-climate variations to N loads and concentrations in streams provides an effective and consistent framework for analysis of the surface water N processes and water quality for large watersheds and basins.
where L DON , L NH + 4 , and L NO − 3 are the dissolved organic, ammonium, and nitrate N leaching from the soil (kg m −2 s −1 ); D s is the water drainage from the active soil layer (kg m −2 s −1 ); ρ w is the water density (1000 kg m −3 ); r DOM , r NH + 4 , and r NO − 3 are dissolved organic matter, ammonium, and nitrate N calibration factors; N DON, av , N NH + 4 , av , and N NO − 3 , av are the concentration of available N in dissolved organic, ammonium, and nitrate N pools (kg m −3 ); N LF , N LS , and N SS are the fast litter, slow litter, and slow soil N contents (kg m −2 ); f LF , f LS , and f SS are the fractions of soluble organic N in the fast litter, slow litter, and slow soil N pools; N NH + 4 and N NO − 3 are the soil ammonium and nitrate N contents (kg m −2 ); b DOM , b NH + 4 , and b NO − 3 are dissolved organic matter, ammonium, and nitrate N buffering factors due to sorption to soil particles; h s is the effective soil depth (m); r c is the C weight content (3.4 %); ρ s is the average soil density (1500 kg m −3 ); C LF , C LS , and C SS are the fast litter, slow litter, and slow soil C contents (kg m −2 ).

Figure 1 .
Figure 1.Structure of LM3-TAN.Two thick boxes show the incorporated denitrification module in the terrestrial component and stream microbial processes in the river component.The river systems are a series of continuously stirred tank reactors (CSTR) that simulate stream mineralization, nitrification, and denitrification.The other boxes show major C and N pools in vegetation (leaves, fine roots, labile, sapwood, heartwood, and N buffer storage), soil (fast and slow little, slow and passive soil, mineral N), and river (organic and mineral N).The arrows depict fluxes of anthropogenic N sources (thick solid), C-N organic compounds and mineral N (thin solid) with associated processes (italic), and C and N lost to the atmosphere or anthropogenic pool (dashed).

Figure 2 .
Figure 2. Overview of the denitrification module.Effects of first-order denitrification coefficient (a), soil temperature reduction function (b), soil water content reduction function (c) on soil denitrification rate, and soil water content reduction function for mineralization and nitrification (d).The curves were produced using Tables3, 6, and 7 inHeinen (2006).

Figure 3 .
Figure 3. Map of the Susquehanna watershed, showing six major sub-basins, main stem of the Susquehanna River, major tributaries (Chemung, West Branch Susquehanna, and Juniata River), streams, and the location of USGS stream gauges and USGS and SRBC nutrient monitoring sites.

Figure 5 .
Figure 5. 20 year (1986-2005) of the simulated stream N loads (normalized by sub-basin areas) at Marietta and Conowingo and the corresponding reported data from SRBC and USGS.

Figure 6 .
Figure 6.Seventeen-year (1989-2005) average simulated and reported (SRBC) stream N loads at the six long-term monitoring sites in (a) millions of kg year −1 and (b) kg km −2 year −1 ; simulated and reported stream N loads for the year 2005 at the 15 monitoring sites in (c) millions of kg year −1 and (d) kg km −2 year −1 .

Figure 8
Figure 8 presents 20 year average (1986-2005) simulated soil water content, temperature, nitrate-N content, and denitrification rate, and these for each of six sub-basins as well as the corresponding sub-basin area, non-point and point N sources are summarized in Table4.An analysis of the six sub-basins shows that the combined effects of land use and climate on the soil denitrification rate, which were the highest in the Lower Susquehanna Sub-Basin (extensive agriculture; Atlantic coastal climate) and the lowest in the West Branch Susquehanna Sub-Basin (mostly forest; Great Lakes and Midwest climate).These results show that the most significant soil denitrification is associated with extensive agricultural land use (non-point sources).The calculated R 2

Figure 9 .
Figure9.Comparison between the calculated and reported budgets of N sources, retention, lost, transport, and river export at the level of the whole Susquehanna watershed for the period 1988-1992(Boyer et al., 2002; Breemen et al., 2002;Seitzinger et al., 2002; USEPA, 2010a).

Table 1 .
Newly introduced or adjusted parameters from the earlier developments.

Table 2 .
Definition of prognostic (PV) and diagnostic (DV) variables and inputs/forcings (IF) used in the equations.Vegetation and soil equations C LF , C LS , C SS PV fast litter, slow litter, slow soil C contents kg m −2 D N DV soil denitrification rate kg m −2 year −1 D s DV water drainage from active soil layer kg m −2 s −1 f LF , f LS , f SS PV fractions of soluble organic N in the fast litter, slow litter, slow soil N pools Gerber et al. (2010)

Table 4 .
Sub-basin area, 20 year (1986Sub-basin area, 20 year ( -2005) )average applied non-point and point N sources, and simulated soil water content, temperature, nitrate-N content, and denitrification rate (% of the non-point N sources) for each of six sub-basins.

Table 5 .
Temporal evaluation of the annual stream discharges and N loads for the period 1987-2005 at Marietta.If a p value is smaller than 0.05, the correlation between the modeled and reported data is significantly different from zero.