Identification of Anthropogenic Impact on Nitrogen Cycling Using Stable Isotopes and Distributed Hydrologic Modeling

Nitrogen contamination of water resources is a significant problem worldwide. Drastic shifts in agroeconomic policy over the past few decades have changed the manner and intensity of anthropogenic nitrogen inputs in the United States and China. It is crucial to gain an understanding of nitrogen cycling within varying environmental contexts and land uses and the interactions among them. Analysis of the spatial distribution of sources and processes affecting the concentration of NO3 and NH4 + is a useful exercise; however, unexpected characteristics of the system can be revealed through reactive modeling. It has been shown that management practices intended to attenuate nitrate in surface and subsurface waters, in particular the establishment of riparian buffer zones, are variably effective due to spatial heterogeneity of soils and preferential flow through buffers. Accounting for this heterogeneity in a fully distributed biogeochemical model allows for more efficient planning of management practices and enforcement of rules regarding releases. Highly sensitive areas within a watershed can be identified based on a number of spatially variable parameters, and by varying those parameters systematically to determine conditions under which those areas are under more or less critical stress. Responses can be predicted at various scales to stimuli ranging from local changes in cropping regimes to global shifts in climate. A distributed hydrologic model, TREX, is presented that provides opportunities to study multiscale effects of nitrogen inputs, outputs, and changes. The model is adapted to run on parallel computing architecture and includes the geochemical reaction module PhreeqcRM, which enables calculation of δN and δO from biologically mediated transformation reactions in addition to mixing and equilibration.


Introduction
In response to increasing population growth and demand for the resources to sustain it, anthropogenic inputs of reactive nitrogen have significantly altered nitrogen cycling on local, regional, and global scales (Galloway et al., 2004, Vitousek et al., 1997. Intensification of agricultural operations in industrialized nations has required the application of biologically available nitrogen in order to increase crop yields to support population growth (Vitousek et al., 1997), and burning of fossil fuels to meet increased global and per capita energy demands results in the atmospheric formation and subsequent deposition of nitrate (Galloway et al., 2004). Increasing pressures on and demand for water resources stemming from population growth have also led to a growing awareness of ecological problems related to anthropogenic inputs of nitrogen and the need for greater understanding of the implications of human-driven nitrogen cycling (Liu et al. 2006).
Understanding the sources and modifications of reactive nitrogen at multiple spatial scales is crucial to addressing the rapid changes in the nitrogen cycle resulting from anthropogenic influences. Stable isotopes provide an effective method for tracking sources of nitrogen at high spatiotemporal resolution, or as an integrator of processes over the course of longer time periods. Modeling the physical and biochemical processes that alter N in the environment can provide insight into the effects of the sources and activities that release nitrogen.

2
To identify these effects, a hydrological model was adapted to calculate the stable isotope compositions of N compounds in surface water. Several features were added to the model in order to allow for more flexibility in parameterization and include a fullfeatured geochemical code in the transport calculations. These adaptations facilitated exploration of the capabilities of the model in simulating anthropogenic effects on the nitrogen cycle at multiple scales, including scenarios with varying contributions of nitrogen from current and historical land use.

Economic Drivers
Patterns of consumption of nitrogen fertilizers in the United States and China, two nations with an increasingly important economic agricultural relationship, drive changes in the fertilizer industry, land use, and policy that promotes intensified production.
Consumption of nitrogen fertilizers in the United States has increased steadily from a level of 2.7 Tg in 1960 to 12.8 Tg in 2011, as shown in Figure 1a below (USDA Economic Research Service, 2016). Increases in population and exports of agricultural products combined with structural changes in the agricultural industry have driven growth in fertilizer applications, while innovation has led to shifts in the type of fertilizer applied. Urea overtook anhydrous ammonia as the leading single component fertilizer in 1998 while applications of sodium nitrate have declined to 10% of 1960 levels.
Ammonium nitrate applications have declined to less than a quarter of the peak level in 3 Economic development and agricultural reforms in China have led to close trade relations with the United States. China is the largest consumer of U.S. agricultural exports in the world, and the United States is the third largest market for Chinese agricultural goods (Gale et al., 2014). China's surging GDP and resulting prosperity have resulted in changes in consumption preferences and land use. Decentralization of agricultural planning in the late 1970s has been connected to increased agricultural output and planting intensity, shifts in land use from commodity to more valuable cash crops, and increased production of animal food products as economic opportunities allow (Lohmar et al., 2009). The resulting increase in applied nitrogenous fertilizer significantly outpaced the contemporaneous increases in U.S. fertilizer applications, as shown in Figure 1 below (Lohmar et al., 2009).

Nitrate
Nitrate levels in surface water have increased or stayed the same in many American watersheds, in spite of the development of stronger urban wastewater controls and agricultural conservation practices (Sprague et al., 2008). Major riverine exports of nitrate increased between 1980 and 2010 (Murphy et al., 2013). Nitrate concentrations normalized for discharge at the mouth of the Mississippi River were higher during lowflow conditions than those during normal discharge periods during those years, which suggests that the dominant input of NO 3into the Mississippi is currently from old groundwater with high rates of nitrate loss, and may take several decades to peak and slowly decline (Chang et al., 2002, Murphy et al., 2013.
Nutrient export from major rivers is a leading cause of coastal eutrophication, culminating in events of hypoxia, defined by dissolved oxygen concentration less than 2 mg/L. Bottom waters isolated from the oxygenated surface by steep gradients in temperature and salinity are depleted of dissolved oxygen through the respiration of excess organic carbon inputs, which result from thick algal blooms in the photic zone.
These hypoxic areas are termed "dead zones" due to the inability of the waters to sustain animal life, and their occurrence and spatial extent is correlated with annual discharge from rivers (Rabalais et al., 2002). Shallow waters covering an area the size of the state of Maryland along the continental shelf in the Northern Gulf of Mexico undergo an annual onset of hypoxia owing to a stratified water column and high organic carbon and nutrient loads from the Mississippi River, as can be seen in Figure 2, below.  (Turner & Rabalais, 2015) Coastal hypoxic areas have been documented in the East China Sea as well, beginning with the Yangtze (Changjiang) estuary in 1981 (Li et al., 2002). The nutrient load from the Yangtze has historically supported a large, and economically important, population of fish, but excess nitrogen exported to shallow coasts is a threat to coastal fishing resources (C.-C. Chen et al., 2007). Nutrient loads in the Yangtze (Changjiang) river have increased by a factor of 10 since the 1980s, promoting eutrophication and harmful algal blooms in the estuary and the associated river plume in the East China Sea. Summer hypoxic events have been recorded frequently, in parallel with the development of the dead zone in the Gulf of Mexico, as seen in Figure 3 below. Water quality in China is alarmingly deficient (Strokal et al., 2016). A 2015 report from the Chinese Ministry of Environmental Protection stated that 42.5% of groundwater monitoring stations throughout the nation reported fairly poor water quality, and 18.8% reported very poor water quality (Chinese Ministry of Environmental Protection, 2016).
In all, over 60 percent of monitoring stations reported water unfit for human contact.
Surface waters are similarly impaired. In 2007, algal blooms due to eutrophication of Lake Tai, a large lake in the Yangtze Delta, prevented the use of the lake as a water source. Water quality has become an increasingly severe problem in the densely populated region surrounding Lake Tai since the early 1990s. Overexploitation of the lake for aquaculture, agricultural nitrogen loss, municipal waste, and atmospheric deposition were cited as the main contributors to the nutrient load leading to the eutrophicity (Townsend-Small et al., 2007). A significant portion (45%) of monitored streams in China currently do not meet "Grade 3 standards," defined as suitable for centralized drinking water, swimming, and aquaculture (1.0 mg / L, for NH4-N as set by the Chinese government (Chinese Ministry of Environmental Protection, 2016).

Ammonia
Ammonia (NH3) is an important component to the biosynthesis of amino acids and proteins in plants and eukaryotes. It is also a primary product of biodegradation of organic matter in soils, and a major component of human and animal wastes. The immediate bioavailability of NH3 or its protonated ion NH4 + promotes its use as a directapplication fertilizer (anhydrous ammonia), or as a component of a fertilizer salt (ammonium nitrate, ammonium bicarbonate). In aqueous solution, pH values typically encountered in soil environments heavily favor NH4 + , however the low boiling point (-33.5 C) and high equilibrium vapor pressure (2.9 atm) of ammonia mean its ready release to the atmosphere. High temperatures and high soil pH promote ammonia volatilization as they shift the equilibrium towards NH3(g).
Global estimates of agricultural NH3 emissions to the atmosphere are 9.4 Tg/yr from manure and 24 Tg/yr from mineral fertilizers (Paulot et al., 2014). The estimated proportion of ammonia emitted by China is 32% of the global total attributed to mineral fertilizers and 20% of the global total attributed to manure, whereas the United States and European Union combined contribute 11% of ammonia emitted from mineral fertilizer and 15% from manure (Paulot et al., 2014). In some areas of the world, depending on the particular agricultural practices of the regions, biomass burning can also be an important 8 source of emissions of NH3 and NO. In the year 2000, global biomass burning produced an estimated 11.77 Tg of NOx, and 10.51 Tg of NH3 worldwide (Lamarque et al., 2010).
Ammonia is readily scavenged from the atmosphere by rain, and so has a short residence time and tends to be deposited near the location where it was emitted. Regional wet deposition tends to closely follow emissions and is concentrated in areas with intensive agricultural activity. U.S. deposition of ammonium in rainfall has increased in parallel with increased application of fertilizer, particularly in the Midwest, as shown in Figure 4.
The highest rates of deposition in the U.S. and Europe tend to be in early spring when fertilizer is applied as corn is planted (Paulot et al., 2014). In China, cropping practices are such that corn is typically planted, and fertilizer applied, later in the year after harvesting of winter wheat crops is complete, and thus wet deposition of ammonia is highest in the summer (Paulot et al., 2014).
Wet deposition represents a significant source of reactive nitrogen in areas that are not otherwise affected by agriculture. Areas with sensitive nitrogen cycles can be disproportionately affected by nitrogen inputs from the atmosphere. Some of the impacts include acidification of soils and cascading ecological effects (Gruber & Galloway, 2008). Ammonium is transformed to nitrate by both autotrophic and heterotrophic organisms through nitrification: Nitrification: NH4 + + 3/2 O2 → NO2 -+ H2O + 2H + NO2 -+ ½ O2 → NO3 -Nitrification is the oxidation of ammonium mediated by autotrophic Nitrobacter sp. A significant portion of NO3transported by streams results from microbial nitrification of decomposed urea, ammonia and ammonium fertilizers (Burns et al., 2004).
In addition to nitrified ammonium, atmospheric sources of nitrate are significant.

Agriculture and Fertilizers: Regulatory Concerns
Major sources of reactive N in surface waters include agricultural amendments in the form of commercial fertilizer and manure applications, animal waste leaching from lagoons and feedlots, municipal wastewater treatment plant discharge, atmospheric deposition in the form of acid rain and dry deposition (direct diffusion of gaseous NH3 or NO2), and the metabolism of nitrogenous organic compounds in soil. Prior to the invention of the Haber-Bosch process, nitrate fertilizers were primarily obtained by mining guano or salt deposits in arid environments (Galloway et al., 2004). After Haber developed the process of converting diatomic N to NH3, fertilizers were increasingly derived from atmospheric nitrogen. The atmospheric source of Haber-Bosch nitrogen is reflected in the ratio of stable isotopes at 0‰ (Macko & Ostrom, 1994).
Nitrate released in past decades might have greater influence on present surface water loads than new nitrogen applications. Elevated flow-weighted concentrations of NO3in the Mississippi during low flow suggest that even during high-flow conditions, the dominant source of nitrate to large river systems is groundwater discharge of reactive N that was infiltrated in the past (Murphy et al., 2013). Nitrate and ammonium that are resident in groundwater or soil organic matter can continue to influence surface waters for several decades. Böhlke and Denver (Böhlke & Denver, 1995) used a suite of geochemical indicators to link NO3dissolved in groundwater to land use during the time the NO3leached into the soil. During base flow, waters discharging into small coastal streams were found to be between 20 and 50 years old, as shown below in Figure 6. In another study, Böhlke and others  linked cropping and fertilization 12 practices over the course of several decades to groundwater nitrate concentrations and isotope compositions, shown below in Figure 7, as well as subsurface and riparian denitrification. Along direct groundwater flow paths where dissolved O2 concentrations remain high, nitrate concentrations and isotopic signatures were unaltered by transport, implying that higher NO3concentrations from more recent agriculture will influence surface water concentrations over the next 20 to 50 years. The primary concern from a regulatory standpoint about high levels of nitrate in drinking water is the effect that it has on human health. Elevated nitrate levels can cause a blood disorder in infants called methemoglobinemia (Cao et al., 2006, Vitousek et al., 1997.
Also known as "blue baby syndrome," methemoglobinemia is caused by nitrate bonding with hemoglobin and decreasing the its ability to transport oxygen through the bloodstream (Vitousek et al., 1997). This condition, as well as cancer and other health risks, have prompted the World Health Organization and regulatory agencies in many countries to establish a limit of 10mg per liter on nitrate (measured as mass of N) in drinking water (Galloway et al., 2004, Nestler et al., 2011.
In addition to the health concerns, there are also ecological effects to be considered.
Eutrophication often results from excess inputs of nutrients and is a concern in aquatic and coastal systems (Galloway et al., 2004, Rabalais et al., 2002. Bottom waters may become anoxic, reducing the capacity of the system to sustain animal life. Excess atmospheric deposition of nitrate results in soil acidification and loss of biodiversity in terrestrial systems (Vitousek et al., 1997). Excess nitrogen in surface waters can be harmful by disrupting the ecological balance of lakes, rivers, and shallow coastal environments. The Gulf of Mexico "dead zone" and Lake Tai in the Yangtze Delta provide some of the most direct examples of the adverse effects of excess nitrogen.
Land use, both historic and current, informs the biological processing of nitrogen in a particular area. Young forests tend to act as a nitrate sink, owing to assimilatory processing of nitrate by plants. Older forests begin to export nitrogen as organic matter becomes available to nitrifiers (Goodale et al., 2002, Zhu & Mazumder, 2008. Soil conditions, including moisture, texture, and organic content, control the capacity of a parcel for processing reactive nitrogen. Denitrification depends on the availability of moist organic-rich anoxic soils. Agricultural fields often have amendments in the form of manure, and thus may have high organic nitrogen content. Other agricultural fields may use commercial fertilizers, typically containing ammonium bicarbonate or urea (Chen et al. 2009).

15
In some areas where agriculture is not the dominant land use type, wet deposition of NO3and NH4 + can be the main source of reactive N during precipitation events (Kaushal et al., 2011).
Anthropogenic influence on nitrogen cycling has strong implications for the availability of water resources in intensively cultivated watersheds, as well as fertility and ecological structure in less exploited areas. The present study seeks to understand the sources and processes affecting the transformation of bioavailable nitrogen through the use of stable isotopes and a distributed hydrologic model.

Stable Isotopes
Isotopes are a powerful tool for partitioning sources and identifying stocks and fluxes of nitrogen compounds within a system based on the biological and geochemical processes associated with different land use types. Isotope compositions, or modifications thereof, indicate processes including transport or metabolism that selectively prefer one isotope over another. These signatures are reported in the delta notation described below: where R represents the ratios of heavy to light isotopes of the sample and standard, respectively. The value is a ratio of ratios and is dimensionless, reported in "units" of per mil, or ‰.
Modifications of the isotope content of a material, or isotope fractionations, result from slight differences in the bond energies of compounds that include heavy vs. light isotopes. Bonds between two light isotopes tend to require less energy to break than bonds between heavy and light isotopes (Peterson & Fry, 1987). This slight difference in bond energy leads to a difference in the rates of reaction for isotopically lighter molecules versus a heavier counterpart. In the case of effectively unidirectional reactions, the ratio of the rate constants for the two reactions defines the fractionation factor β, and the derived kinetic isotope effect, ε, of the overall reaction. The kinetic isotope effect is usually expressed in units of per mil (‰) to describe the instantaneous fractionation in terms of the isotope composition of the product: where k is the rate constant for the lighter isotope and k' is the rate constant for the heavier isotope.
For a chemically closed system, i.e. the supply of reactants is finite, the instantaneous ratio of light to heavy isotopes changes as the substrate is consumed. This change can be described by where R is the stable isotope ratio of the remaining substrate, R0 is the initial stable isotope ratio and f is the fraction of the substrate remaining (Macko et al., 1986).
Rewritten in the more commonly used δ notation: The implication of this relationship is that the substrate becomes progressively more enriched in the heavy isotope as the substrate is consumed in the reaction.
The inverse of this relationship is true for the product of the reaction: This shows that the initial product is depleted in the heavy isotope and progressively approaches the original isotope ratio of the substrate. If all of the substrate is consumed, the product will have the same original isotope ratio as the substrate. The kinetic rate constants, and the ratio between them, form the basis of this fractionation and are influenced by temperature, among other environmental factors. The kinetic isotope effect, ε, is subject to environmental controls. At higher temperatures or faster rates of reaction, the fractionation tends to be less, whereas at lower temperatures, higher concentration of substrate or slower rates of reaction, isotope fractionation tends to be enhanced (Mariotti et al., 1981).
Oxygen isotopes included in nitrate ions show promise as indicators of oxidation processes that led to the production of the nitrate. Distinctive pathways for the daytime  (Elliott et al., 2007): however, both depend heavily on the incorporation of O from O3, which shows typical signatures of δ 18 O-O3 at 90‰ to 120‰, compared with typical δ 18 O-O2 of atmospheric oxygen at 23‰ (Elliott et al., 2009). As a result, nitrate produced by the atmospheric oxidation of NOx by ozone and deposited in precipitation is heavily enriched in 18 O relative to other sources, as seen in Figure 9 below.
Previous research has shown that reactive nitrogen compounds derived from natural and anthropogenic sources exhibit characteristic ratios of stable isotopes of nitrogen and 20 oxygen (Chang et al., 2003). Typical sources of nitrate in aquatic environments and ranges of isotope ratios are depicted in Figure 9 below (Chang et al., 2002). Processes that transform these compounds within the environment also alter isotopic signatures in predictable ways (Barnes et al. 2008;Chen et al. 2009;Pardo et al. 2004).
For example, microbial denitrification preferentially converts nitrate containing 14 N and 16 O to diatomic nitrogen, leaving a substrate of nitrate that is enriched in both 15 N and 18 O. A dual isotope approach further reveals a characteristic ratio of 2:1 δ 15 N : δ 18 O when comparing nitrate sources with residual nitrate following denitrification (Chang et al., 2002, F. Chen et al., 2009.

Biological Nitrogen Transformation: Nitrification and Denitrification
Dissimilatory transformation of reactive nitrogen includes several biologically mediated processes, each with characteristic effects on loading and isotopic ratios of reactive nitrogen compounds. Studies have established that 5/6 of oxygen making up product NO3comes from H2O whereas the remainder may be derived from soil O2 (Casciotti et al., 2002, Lehmann et 22 al., 2004. Nitrate produced by nitrification is depleted in 15 N compared to the substrate, and O resembles meteoric water, close to 0‰.
In order to account for isotopic kinetic effects, the rates of reactions of individual isotopologues, or elementally identical compounds containing different isotopes, must be taken into consideration. The general case for transformation of reactive nitrogen species can be represented as: where δt 15 N represents the relative isotope ratio of the remaining substrate at time t.
Biologically mediated nitrate reduction is a multistep process carried out by communities of heterotrophic bacteria and archea. Enzymatic catalysis of each reaction step promotes 23 differing reaction rates depending on the isotopic composition of each chemical species.
Reactants containing lighter isotopes are therefore preferentially removed, leading to an increasingly enriched pool of reactants (in a closed system) and an initially highly depleted pool of products, that becomes increasingly enriched as the reactants have progressively less of the lighter isotope. If the reaction proceeds to completion, the isotopic composition of the products will match the composition of the reactants, resulting in 0 net fractionation. However, at states of partial completion, products will show significant isotopic depletion relative to reactants.
The bulk of denitrification in streams occurs in the hyporheic zone, where surface water interfaces with groundwater underlying the streambed. Small, shallow upland stream segments that incorporate a greater proportion of flow through the hyporheic zone provide biological processing of reactive nitrogen whereas larger rivers tend to act as conduits for seaward transport of dissolved N (Chang et al., 2002).

Site Background
In order to examine the impacts of the changes in the nitrogen cycle, a site was selected In order to partition the sources of ammonium and nitrate based on stable isotope signatures within the Najinhe watershed in northeast China, this study implements a spatial approach using GIS and distributed hydrologic modeling with a reactive isotopegeochemistry transport modeling component using nitrogen and oxygen stable isotopes.
The present work seeks distinctive signatures of land use types in nitrogen inputs to the watershed, and attempts to prove that stable isotopes provide an effective method for tracking sources of reactive nitrogen in a spatially explicit manner.
Because land use is such an important determinant of N concentrations and isotope signatures, it must be addressed in a spatially explicit manner in order to account for how the different types of land relate to one another in processing or releasing reactive N.

26
Each of these processes and sources has a distinctive signature in the isotope ratios of the N compounds that can be used to try to ascertain its role in the system.

Methodology Isotope Ratio Mass Spectrometry (IRMS)
Whole soil isotope ratios were determined by use of a to mass 28 N2 and compared to the ratio of the same masses in a standard derived from atmospheric N2 using the delta notation described earlier in this document.

Hydrologic Modeling
As N enters the surface water, the ratio of wet-deposition to N derived from terrestrial sources (manure, commercial fertilizer, soil organic N) changes throughout the runoff event (Kaushal et al., 2011). The δ 15 N and δ 18 O of N species dissolved in the water column evolves along with the ratio of the sources, so information about the timing of sampling in the context of rainfall events can help to distinguish sources of N using the instantaneous stable isotope ratio rather than the average stable isotope ratio of dissolved HSPF similarly includes functional components, often as subwatersheds, that can connect to each other to define the system as a whole. The lumped parameter approach is a useful way to model watersheds with static or well-constrained component subwatersheds.
However, defining a separate subwatershed for each field in a dense, multi-use watershed is inefficient.
The diffusive wave approximation of surface water flow is a well-established method for modeling flow in shallow water. Computationally the equations are considerably simpler, and thus faster, to solve than the full dynamic St. Venant equations. The primary assumptions of the diffusive wave approximation are that inertial acceleration terms are negligible compared with body forces such as friction and gravity. As a result, the diffusive wave approximation operates under the assumption that flow is subcritical.
In order to adequately simulate the biogeochemical processes that control nitrogen cycling, a distributed watershed model is required to model transport, and a spatially explicit chemical reaction facility is needed that can predict δ 15 N and δ 18 O of N species.
None of the models explored contained this combination of capabilities as a default.
However, TREX, in addition to the fully distributed approach, offered easy access to the source code, allowing modifications to improve performance and adaptability of the model.

TREX Model
The TREX model (Velleux et al., 2008) is a spatially distributed watershed model that accounts for continuous spatial and temporal variation in soil and land-use characteristics, precipitation and soil moisture, and chemical concentrations in soil and water, among other parameters. Infiltration is formulated as a Green-Ampt model where hydraulic conductivity is defined as a characteristic of underlying soil / sediment (fully distributed in space) where f is the infiltration loss rate, Hw is the depth of water above the ground surface, Hc is the capillary pressure head at the wetting front, θe is the effective porosity, found by subtracting residual moisture from total porosity, Se is the effective saturation of the soil, and F is the total infiltration depth.
Hydraulic parameters associated with land use and underlying soils affect overland and channel flow in terms of friction, available interception and storage depth, and susceptibility to infiltration and erosion. The Manning formula is used to determine the energy loss due to friction: where V is the average velocity, k is a unit conversion factor, n is the Manning roughness coefficient, which describes the velocity response to the roughness of the channel or bed material, Rh is the hydraulic radius, and S is the friction slope, which accounts for hydraulic head loss due to friction.
The Manning equation determines discharge from one cell to the next for a given ground slope and water depth. The Manning roughness coefficient n is defined for each land use class, and porosity and hydraulic conductivity are defined for each soil group.
Precipitation can be defined in the TREX model in several different ways. There are options to define a single, uniform time series of precipitation intensities, or to define a series of distributed grids that vary local intensities throughout the run time of the model. This allows the model to more flexibly determine the water budget of the watershed, accounting for the spatial variability of precipitation events. Precipitation state responds to environmental parameters defined in the model, so precipitation that falls when the temperature is below 0 °C falls as snow and melts in response to solar insolation or temperature changes.
As is evident from the water transport equations, TREX accounts for water losses to infiltration, interception, and predefined sinks. Evapotranspiration is not modeled. The environmental parameters needed to calculate losses to evapotranspiration are available as time-varying inputs to the model; however, the model is oriented towards event timescales over the course of hours, where losses to evapotranspiration are not expected to be significant. Should evapotranspiration be included in future versions of the model, a module to determine the detailed soil moisture profile with the Richards equation, an alternative infiltration mechanism that can describe saturation in the soil profile in greater detail and account for capillary motion due to drying at the surface, could prove to be useful.
TREX solves for derivatives of hydrologic and transport variables, then integrates the state using the Euler method, in which instantaneous rates of change are applied forward in time to solve for the state of a variable one time step in the future. For example: where h is the height of the water in a grid cell, qin is the specific discharge entering the cell, qout is the specific discharge leaving the cell, and dt is the length of the time step.
Time steps (dt) can be calculated dynamically based on the Courant condition, calculated as velocity *dt / cell size. The Courant condition is a necessary condition of numerical stability in the explicit solution of partial differential equations (Courant et al., 1928).
Intuitively, if the flow between two grid cells exceeds the spatial dimensions of the cells, the calculated depth at the next time step may be invalid. Specifically, the possibility of negative water depths arises as more water may be simulated to exit a cell than the cell could possibly contain.
To determine if the current time step is appropriate to solve the water routing equation, the model calculates the Courant condition at each cell for each iteration of the model, and if the condition is violated in any cell, routing is recalculated using a smaller time step. This process iterates until an appropriate time step is determined.

Sediment Transport
Sediment transport capacity is calculated by a modified Kilinc and Richardson model that makes use of Revised Universal Soil Loss Equation (RUSLE) (Renard & Foster, 1991) parameters in the overland domain: Particle classes are subject to transformation processes, including abrasion and dissolution to yield particles of a different size class, or in the case of dissolution, yield chemicals dissolved in the water column.

Chemical Transport
Chemical transport processes are formulated in a similar manner to their sediment transport counterparts. Advection and dispersion are the simulated processes that control concentrations of chemical species in the water column in response to hydrologic forcing.
Chemical components are subject to some additional processes, however, and may be partitioned between the dissolved and sorbed phase based on an equilibrium constant defined by the user: = 1 1 + where fd = the fraction of the chemical in the dissolved phase, Kd is the distribution coefficient, and m is the mass of sorbent particles in the water column.
Chemical components are also subject to some transformation processes in the water column, including biodegradation, photolysis, radioactive decay, and hydrolysis, that change concentrations in situ and may yield other chemical types. The treatment of these processes is fairly rudimentary, however, and complex solution equilibrium calculations are not practical.

Adaptations to the Model
Several additional features were added to the model in order to increase computational efficiency, allow for more flexibility in parameterization, and include a full-featured geochemical code in the transport calculations. These adaptations provided opportunities to explore the capabilities of TREX in simulating anthropogenic effects on the nitrogen cycle at the catchment scale.

Performance
Finite difference models in particular are of a class of computation that is, in the language of computer modeling, "embarrassingly parallel" in that grid cells are treated identically and can therefore be calculated simultaneously for each timestep. Profiling the execution times of the model revealed processes that were especially computationally intensive, so those tasks were parallelized to take less time.
Computationally intensive tasks were parallelized using OpenMP, a shared memory, thread-directed library rather than the distributed memory oriented Message Passing Interface (MPI). While using MPI could potentially lead to scalability and further performance increases in certain configurations of the model, the adaptive timestep functionality would require synchronization at each time step, massively increasing the communication overhead for parallel processing.
Shared memory parallelization incurs a thread creation penalty and a process limit beyond which cache collision problems tend to negate gains from parallel computation.
Careful consideration of cache geometry and memory architecture can mediate some of these problems, but there tends to be a finite number of processes beyond which performance degrades rather than improves (Agarwal, 1992).
While a significant number of model processes were parallelized, the bulk of the code remains serial. According to Amdahl's law (Amdahl, 1967), which determines the maximum speedup possible for the program after parallelizing sections of code, the performance improvement expected from parallelization of computer code is limited by the proportion of computational time spent in serial processing (Gustafson, 1988). Based on profiling tests of the code, the maximum speedup of the program that could be expected from processing parallel sections with four computer processor cores is 1.6x.
The actual measured speedup using the demonstration model is 1.2x which allows for more rapid and responsive modeling.
The dynamic time step implementation presents a challenge in implementing code parallelization in the model. The most obvious domain over which to divide the computation is the spatial grid. However, the spatial heterogeneity of computed velocities means that, after the grid is divided among processors, the Courant condition, which specifies that velocity *dt / cell size must be less than one, might be (and often is) violated for only one processor. In order to account for that scenario, the entire spatial domain must be synchronized after each time step. Synchronization operations are often time-intensive on distributed memory, multi-node architectures, and severely reduce the efficiency of code parallelization. It is for this reason that parallelization efforts in the model code are focused on shared memory strategies such as those implemented in the aforementioned OpenMP.

38
Adaptability Input and initialization data were originally formatted in an inflexible format that had the appearance of keyword descriptions, but actually was more dependent on line breaks and the order of the arguments presented. The input format has been changed to an Extensible Markup Language (XML) document in order to allow for more flexibility in ordering input parameters and grouping them in a hierarchical structure. Input files are more self-descriptive, and input dependencies more obvious, by using tags, comments, and hierarchy to define the model rather than a strictly ordered list with no selfdocumentation available.
A particular advantage to the XML document structure is the ability to script parameter sweeps for sensitivity analysis and calibration. Template documents can enforce interdependencies among commands and model parameters, markedly easing the process of generating valid input code. Software libraries designed to generate and parse XML documents are available for a wide array of languages and environments.
A severe limitation of the TREX model is the inability to save and restore the state of the model in case of interruption of execution. This functionality was added in order to complete intensive model runs on a managed computing cluster in spite of limited run times. In order to guarantee that the state of the variables in the model is consistent when stored (i.e. either all variables have been updated or none of them have), data are stored in a database file that can manage database transactions through the use of a journal file.
In this study, a freely available open-source software library called SQLite was used for the database journal file. To improve computational scalability, the SQLite interface could be adapted to a full-featured database management system. An added advantage of the use of a database format is a much more simplified coupling of input parameters and output grids. All inputs and outputs are now stored within a single file, easing organization of data and reproducibility of results.

Chemistry
The complexity of accounting for all of the environmental parameters that control biotic and abiotic transformations of nitrate suggests the use of a dedicated chemical reaction module. One such module, PHREEQC (Parkhurst & Appelo, 2013) is a general purpose aqueous geochemical model that allows the calculation of equilibration, redox, and kinetic reactions in aquatic systems, including the calculation of selected isotope fractionations.
The authors have developed a reaction module version of the PHREEQC code, called PhreeqcRM (Parkhurst & Wissmeier, 2015), that allows access to the geochemical facility while reserving transport calculations such as advection and dispersion to the hydrologic model that interfaces with it. One significant advantage of the use of a spatially distributed reaction module is that chemical reactions can be defined not only by chemical inputs, but by location as well. The reaction module accepts input defining the initial conditions, reactions, and chemical parameters to be included in the model, and allows for a more generally applicable chemical reaction model than that included in TREX. Reactions may be further defined by the user in order to incorporate transformations or species that are not included in the default database.
PHREEQC includes a set of databases that defines thermodynamic constants of various elements, compounds, and reactions, in addition to parameters of reactions involving isotopologues, or compounds with the same chemical structure but different isotopes, such as CO2 and C 18 O 16 O. The isotopic fractionation of nitrogen by denitrification is described by default; however, the simultaneous fractionation of oxygen is not described.
To address this limitation, additions to the database can be included in the input file, offering yet more flexibility in defining the transformations that describe the system being modeled.
The software library PhreeqcRM is a reaction module designed to be incorporated in hydrochemical transport software as a chemical reaction facility. Interfaces have been added to directly access reaction parameters such as temperature, concentrations of reactants, volume of water, etc. The use of a separate reaction module is particularly suited to a sequential noniterative, or operator splitting approach to the numerical solution of reactive transport: where differential equations governing mass transport, including advection and diffusion of dissolved species, are solved separately from chemical reaction equations (D. J. Z. Chen & MacQuarrie, 2006). The assumption inherent in this approach is that transport between t and t + dt occurs instantaneously so that the chemical components have no time to react during the time of transport. An alternative formulation that attempts to account for this assumption would be the Strang-splitting approach, where half of a transport timestep is calculated, the chemical reaction for a whole time step is calculated based on the resulting concentrations, and then the remainder of the transport calculation is completed. Although the current model formulation does not include the Strang-splitting scheme, it is a goal for future development.
PhreeqcRM is written to take advantage of parallel computing architectures on multiple scales. A shared-memory implementation, which solves reaction cells on multiple threads within the same computing process, can be nested within a distributed memory implementation that solves reaction equations across multiple computing nodes in a large cluster. The scheme used in the coupling with TREX is the OpenMP shared-memory platform due to the previously mentioned challenges of synchronizing model calculations when using dynamic time steps to satisfy the Courant condition.

Experimental Parameters
Initial non-reactive modeling will predict stable isotope ratios based on the contribution of individual sources. Concentrations of each component of the NO3or NH4 + mixture are summed and weighted with a characteristic isotope signature to find the predicted isotope ratio based solely on mixing: where C is the concentration of NO3or NH4 + for each source, and fi is the fraction of the total nitrate or ammonium concentration from a given source. The temporal and spatial patterns of the isotopic signatures of dissolved N species are analyzed to determine relationships with input variables including concentration and distribution of source N, rainfall intensity throughout precipitation events, and distributed hydraulic and environmental parameters. Important considerations include the influences of point and non-point sources and potential management interventions.
The calculated isotope ratio will reflect the dominant source of dissolved N at each location in the model domain, and will be compared to measured bulk isotope ratios of soils in the watershed.  Figure 15 whereas hydrologic controls such as roughness and infiltration capacity are derived from the parameters in Figure 16.

Bulk Stable Isotopes
Bulk stable isotope analysis of nitrogen in soils of the Najinhe watershed reveals spatial     Figure 18 depicts measured isotope ratios grouped by sampling location and by land use.
Although bulk isotope ratios may exhibit an explicit response to land use, sample sizes are not sufficient to establish a statistical separation among land use classes.
Preliminary results of the flow model based on spatially uniform precipitation events are shown below. Figure 19 illustrates runoff generation in the Najinhe watershed.
The flexibility of the distributed approach allows for statistical methods that simulate parameter sweeps on a grid cell-by-grid cell basis.
Based on results from flow and reactive transport simulations, cells can be identified that During the early stages of a precipitation event, the initial accumulation of surface water in smaller upland streams is apparent in the image at the left in Figure 19.

Model outputs
The results of the model calculations using TREX and the previously described parameters for Goodwin Creek are presented in the section below. The data represent results from two distinct scenarios determined by the presence or absence of nitrate in groundwater supplied to the steam network as base flow. Antecedent nitrate concentration in the stream network is intended to highlight the contribution of old, previously infiltrated ammonia fertilizer that has been nitrified to ammonia without significant removal and isotopic fractionation due to denitrification.   Monitoring station 3, represented in Figure 24, collects water from a small subcatchment and is therefore considerably more responsive to inputs from the land surface throughout the simulation. Like station 1, the first measured nitrate concentrations are consistent with a predominance of manure. However, unlike at station 1, the manure dominated signature also corresponds to the peak concentration of nitrate measured at this location.
As stream discharge increases to its peak at this station, the dominant nitrate source quickly transitions to atmospheric nitrate. Between t = 12 h and t = 15 h, clean base flow rapidly flushes all remaining nitrate from the channel. The smaller contributing area allows for more detailed examination of nitrate sources at shorter timescales and highlights input of runoff containing manure and fertilizer as a significant contributor to the nitrogen budget.
Monitoring station 4 begins the simulation with predominantly atmospheric nitrate, then briefly transitions to include a mixture of atmospheric and manure derived nitrate before returning to a source made up mostly of atmospheric nitrate. Like monitoring station 3, concentrations and isotope signatures are responsive to overland sources. However, like monitoring station 1, the peak contribution of manure nitrate does not correspond to the peak concentration of nitrate or the peak stream discharge. The reestablishment of base flow concentrations takes place over considerably longer time than at station 3.
All of the monitoring stations exhibit a constant mixing line between completely atmospheric sources and an invariant mixture of fertilizer and manure sources within the stream network. The constant fertilizer-manure mixing ratio reflects the rapid mixing that takes place during overland transport and subsequent flow into the stream network.
There is no fertilizer or manure source independent of the other within the stream network that might change this ratio once it has been established. A second simulation is defined to simulate scenarios akin to that described by Böhlke and Chang, where old, groundwater-borne nitrate is the dominant source rather than nitrate washed from the surface and contributed by runoff. The second simulation differs from scenario 1 in that stream nitrate is initialized and supplied by base flow at 10 mg/l in the form of nitrified ammonium fertilizer. The contrast between the two scenarios is most apparent when examining the nitrogen and oxygen isotope ratios in Figure 26 at t = 1.75 h. In scenario 2, nitrate contained in the stream is distinct from that dissolved in water flowing overland whereas nitrate in the stream network in scenario 1 reflects the signature of the dominant atmospheric source for the same time period. During heavy rainfall, atmospheric nitrate predictably begins to dominate within the stream network, however, the proportion of nitrified fertilizer nitrate quickly rebounds as base flow proceeds to flush atmospheric nitrate through the network.

Conclusions and Next Steps
The use of a distributed hydrologic model is likely to provide most insight in areas with greatest population density or in areas where water quality is threatened. China, with the highest population in the world, faces critical regional water shortages year after year, as well as severe problems with water quality. The Najinhe watershed, located in a remote region of the country, is used for intensive agriculture, and as a result, water quality is a matter of concern. Future applications of the TREX model in areas beyond the Najinhe may have impact on a larger scale. The use of the model to inform the partitioning of nitrate among different flow paths may show potential for increasing scope to a basin scale model in order to evaluate nitrogen cycling within the Yangtze watershed. The Yangtze watershed is the third largest in the world, and many areas within it are densely populated, creating intense demand on the water resources.
An effective method and scalable model could have greater application to different sites and populations beyond China as well. Chang et al (2003) attempted to use an areal distribution of land uses in the Mississippi River basin to predict isotope ratios of nitrate without accounting for the order in which land use types were encountered downstream.
A distributed hydrologic model should allow more precise treatment of the interaction of adjacent land uses.
The modeling scenarios in the present study represent initial explorations of the capacity of modeling to describe the dynamics of reactive N in surface water. There are, however, a number of capabilities of the distributed modeling approach, and the model developed as part of this work, that will prove useful to future investigations.
In one application of the additional capabilities of the TREX model, point sources of N compounds may be included to simulate discrete outflows or release events in order to determine the effect of inputs to the system that are more isolated in time and space. In addition to temporally continuous point sources, such as wastewater treatment facilities, combined sewer outfalls, and animal feeding operations, future applications of the model may determine the role of illicit dumping of manure or other nitrogenous waste into the stream network. Recent work suggests that direct input of manure may contribute the majority of dissolved nitrogen in many Chinese subbasins, and the contribution will likely increase over the coming decades (Strokal et al., 2016). The use of a distributed model that includes spatiotemporally distinct sources of N in combination with non-point sources can help to resolve the importance of these inputs in the N budget of the system.
TREX combined with PHREEQC can also be applied in more general cases to determine the effects of changing environmental parameters, cropping practices, or trends in groundwater chemistry. For example, the Najinhe watershed is subject to the East Asian Monsoon, and therefore highly seasonal rainfall. The strength of the monsoon has been shown to be variable over the course of several decades (Ha et al., 2012), and we can select model parameters to reflect this variability and its effect on the nitrogen cycle.
Gaining a better understanding of how agricultural practices and land management affect water resources within a watershed has implications for potential improvements to both crops and natural resource planning and management. Knowing the source of a contaminant such as nitrate can enable efforts to remediate the effects and reduce or prevent further pollution. The establishment of riparian buffer zones, installation of 64 bioreactors, and protection of wetlands can all help to remove nitrate from runoff and improve water quality (Nestler et al., 2011). A more complete view and comprehension of the transformation of nitrate within a particular watershed will provide the knowledge necessary to best implement these tools. The use of validated water and sediment processes to define transport of chemical species combined with general-purpose geochemical modeling to simulate reactions offers opportunities to explore scenarios and discover relationships between land use, conservation policy, and water quality.