The nitrogen legacy: emerging evidence of nitrogen accumulation in anthropogenic landscapes

K J Van Meter; N B Basu; J J Veenstra; C L Burras

doi:10.1088/1748-9326/11/3/035014

1. Introduction

Human modification of the nitrogen (N) cycle has resulted in increased flows of reactive N (N_R), with growing evidence that planetary boundaries for maintaining human and ecosystem health have been exceeded [1, 2]. The creation of large hypoxic zones, and the resulting loss of habitat and species diversity in estuarine and coastal marine ecosystems, has been one of the most significant impacts of such increased flows [3]. While the need to manage N flows and their associated ecological impacts is increasingly recognized, implementation of conservation measures to reduce stream N concentrations has had only limited success [4, 5].

Growing evidence suggests that this lack of success can be attributed to diffuse legacy sources that continue to impair water quality even after agricultural inputs have ceased [6, 7]. These sources can lead to time lags between management changes and measurable improvements in water quality, lags that can make it difficult to evaluate the effectiveness of the management practices employed or to maintain public support of costly, ongoing interventions [4, 8, 9]. Such time lags, which have been defined as the time between the initiation of a restoration practice and the point at which a change is observed in the target water body [4, 9], have been observed in Europe and the United States, where nitrate concentrations in streams and aquifers have remained high despite reductions in N loadings to watersheds [10–12].

The presence of legacy sources is also suggested by the frequent references to 'missing' N, also referred to as N retention [13], in mass-balance studies of intensively managed catchments [14, 15]. In such catchments, anthropogenic inputs of N routinely exceed measured outputs, creating watershed-scale N budgets that appear significantly out of balance. Indeed, both regional and continental-scale studies suggest that an inefficient use of N is common in heavily agricultural watersheds, leading to a large N surplus (defined as N inputs − usable outputs) [13, 16, 17]. A portion of this N surplus exits the watershed as riverine output, while the fate of the residual N, although not wholly unknown, remains largely unquantified at watershed scales. In particular, denitrification and subsurface storage constitute well-known pathways by which N may either exit a catchment or be retained over a long period (c) and these N sinks are frequently grouped under the heading of 'N retention' (figure 1). Our synthesis of N mass balance studies for watersheds across the world shows a mean N retention of approximately 50 kg ha⁻¹ yr⁻¹ (supplementary table S1), but, as discussed below, precise quantification of N fluxes via specific retention pathways has remained elusive [18, 19].

**Figure 1.** Schematic showing the stores and fluxes of reactive N in a human-impacted watershed, explicitly including both point and non-point sources.
Download figure:
Standard image High-resolution image

Denitrification, which occurs in both soils and stream sediments [20, 21], is the process by which N_R is removed from a system via reduction to nitrous oxide (N₂O) and nitrogen gas (N₂) [8, 22–24]. Due to inherent difficulties in direct measurement of denitrification products, considerable uncertainty exists regarding denitrification rates in terrestrial systems [18, 22]. As a result, denitrification is often used simply as a balancing term in mass balance studies, with denitrification rates being estimated based on differences between N inputs and all other N storage and loss terms for the watershed [14, 25, 26]. The existence of such a balance, however, is based on an assumption of steady-state dynamics for terrestrial N reservoirs, with all anthropogenically and naturally fixed N_R being denitrified and returned to the atmosphere on an annual timescale [22]. Although such an assumption has been hypothesized to be valid for pristine systems over long timescales [27], it has been shown to be no longer applicable with the current high inputs of NR in intensively managed landscapes [19, 23, 28, 29]. Indeed, modeled estimates of of denitrification are often significantly significantly lower than those suggested by national-scale mass balance-based estimates [30].

The other possible fate of the 'missing N' is storage within the subsurface. We can conceptualize the subsurface environment to be composed of three major N pools: (1) dissolved NO₃^- in the vadose zone or (2) in groundwater aquifers, and (3) organic N within the soil profile (figure 1). Large vadose zone stores of inorganic N have been demonstrated in desert and semi-arid regions, with accumulation magnitudes in deep vadose zones (30–50 m) varying as a function of rainfall, tillage and irrigation history [31–33]. The existence of a significant groundwater reservoir has been proposed based on observations of increasing groundwater N concentrations over time in both the US and Europe [8, 34]. Although the existence of such subsurface reservoirs for N is well accepted, determination of the magnitude of N accumulation is subject to significant uncertainty due to the presence of complex aquifer systems and difficulties in measuring spatially varying patterns in ${{{\rm{NO}}}_{3}}^{-}$ concentrations and groundwater storage [35]. In one of the few studies attempting to quantify stores of groundwater N over time, Worrall et al [34] estimate that N accumulation in groundwater beneath the River Thames Drainage Basin in the UK reached a peak between 2000 and 2004 of $1571\quad \pm \quad 608\,{\rm Mg}$ N.

The third potential subsurface storage reservoir is organics N held within the soil profile (figure 1). Indeed, the largest pool of N in most terrestrial ecosystems is soil organic N (SON) [36, 37], and at current levels of N input, it is suggested that terrestrial N sequestration may be occurring at a global scale on the order of 20–100 Tg N yr⁻¹ [38–40]. The potential for increased N storage in forested soils has long been accepted for landscapes subjected to elevated levels of atmospheric N deposition [41]. For example, a study of N retention and C sequestration in European forests estimated N sequestration in forested soils to have occurred at a rate of 4.7 kg N ha⁻¹ yr⁻¹ from 1960 to 2000 [42]. It has similarly been suggested that N may be accumulating in agricultural soils [41, 43]. Yan et al [44] found the average soil N content of Chinese croplands to increase by 5.1% between 1979–82 and 2007–08, while mass balance and modeling studies in Canada [30], Europe [45] and the US [46] suggest an annual accumulation of N within agricultural soils on the order of 15%–20% of total N (TN) inputs. Fenn et al [41] have shown that soils in which C and N pools have been reduced by disturbance, such as those under agricultural cropping, may exhibit the highest levels of N retention. Smil [43] has estimated that in agricultural soils receiving regular fertilizer inputs, N accumulation is likely occurring at a rate of 25–35 kg ha⁻¹ yr⁻¹, and Worrall et al (2015), in their recent study of the Thames basin, suggest that SON has accumulated at a rate of 55 kg ha⁻¹ yr⁻¹ since 1973.

In general, however, little attention has been given to the possibility of soil N storage in the context of watershed-scale N balance studies, primarily due to assumptions of either ongoing N depletion [47] or steady-state dynamics under conventional agriculture [25, 48]. For example, although Billen et al [47] note that storage in the soil organic matter reservoir could potentially account for missing N in the soil N budget for the Seine and Somme watersheds, this possibility is discarded because 'soil organic matter content is generally considered as decreasing due to continuous cropping'. Howarth et al [25] explicitly assume no potential for soil N accumulation, noting that after a large net release of N following conversion of land to agricultural use, the N status of soils reaches a steady state, with N immobilization, on average, equaling N mineralization on an annual basis.

Indeed, it is well-documented that dramatic losses of SON and C can occur after cultivation, particularly in nutrient-rich soils like those found in the North American prairie region [49–54]. As early as 1905, for example, it was reported that Canadian prairie soils had lost more than 20%–30% of the organic matter originally present in the plow layer [55]. This fast depletion trajectory is due in part to a loss of physical protection provided by soil aggregates [56], with cultivation breaking up aggregate structures and leading to increases in oxidation and mineralization rates [52]. After these initial losses, however, SOM has been found to stabilize [57, 58], and it has been proposed that such losses could be reversed in response to the ongoing addition of root matter and other crop residues to soil [52]. It is this period, after stabilization, when it has been proposed that accumulations can occur, that is the focus of our study.

Our central hypothesis is that decades of high-input agriculture have led to a significant accumulation of SON within the landscape and that this accumulation may contribute to time lags in catchment response after changes in management practices. Our objective is (1) to use historical and current (mid-20th century to present) soil sampling data to provide direct evidence of potential changes in soil N content over time, (2) to place such evidence in the context of watershed-scale mass balance studies, and (3) to develop a parsimonious modeling framework to explain decadal-scale changes in SON. Our specific focus is on agricultural soils of the Mississippi River Basin (MRB), an area that covers approximately 41% of the contiguous United States and includes more than 800 000 km² cropland, much of which has been under intensive cultivation since the mid-19th century [59]. Thus, our paper focuses on answering the following two questions:

Is N accumulating in agricultural soils, and if so, in what form, and in what magnitude?
What are the implications of such accumulation with respect to time lags between changes in management practices and water quality benefits?

2. Materials and methods

2.1. Soil resampling studies

We synthesize data from two studies, the first (carried out by the authors) in Iowa and the second a smaller study in Illinois [50], both of which were designed to assess anthropogenic changes in agricultural soils of the United States Midwest over multiple decades. Details on sample collection and analysis methodologies for the two studies are provided in supplement 2. In both studies, soil cores were obtained from plots under row crop agriculture that had been previously sampled in the mid-1900s and analyzed for TN content. Accumulation or depletion was estimated as the difference between the current and the mid-1900s N content. Such a resampling approach has commonly been employed to assess changes in soil C stocks over time [58, 60], but has not been broadly utilized to evaluate potential changes in soil N.

In the Iowa study, soil samples from 61 representative pedons belonging to 46 different soil series in 21 counties across Iowa were obtained in 2007. These sites were previously sampled as part of the National Cooperative Soil Survey (NCSS) (between 1943 and 1963, median sampling year 1959), and all but three of these sites remained under intensive cultivation during this time frame (see Veenstra [61, 62]). Data from the Illinois study is based on samples obtained from six sites in central Illinois located on poorly drained Mollisols that were under corn-soybean rotations, were tile-drained, and had no history of manure application [50]. All six sites were originally sampled in 1957 and resampled in 2002.

**Figure 2.** Accumulation of TN in the subsurface based on historical resampling studies of agricultural sites in the Mississippi River Basin (MRB). The top left panel shows the MRB with the location of the sites in Iowa and Illinois, while the top right panel shows the 61 sampling locations for Iowa. The bottom left panel shows the mean TN accumulation (kg ha⁻¹ soil) between 1957 and 2002 for the six sites in Illinois [19]. The bottom right panel shows the mean TN accumulation (kg ha⁻¹ soil) for the 61 Iowa sites [20]. Error bars in both plots correspond to the standard error of the mean. Both studies show net N accumulation across the soil profile, with the majority of accumulation occurring from 25 to 100 cm.
Download figure:
Standard image High-resolution image

2.2. Trend analysis of soil data across the MRB

The resampling study described in the previous section was used to quantify N accumulation at specific locations based on two points in time, an approach commonly used for the assessment of C sequestration in plots under long-term tillage [63]. We complemented this historical resampling approach with analysis of NCSS soil samples [64] obtained across the MRB from 1980 to 2010 to test for negative or positive trends over time in TN. Reported values for bulk density and TN were standardized to depth layers of 25 cm (0–25 cm, 25–50 cm, 50–75 cm, 75–100 cm). TN concentrations (g Mg⁻¹) were obtained directly from NCSS chemical analysis data [64], reported in the database as percent N. Area-based estimates of TN content were calculated from the thicknesses of the soil layers and bulk density values.

**Figure 3.** Accumulation of TN in agricultural soils across the MRB, 1980–2010, based on 2069 soil samples from the NCSS database. (a) The number of samples used for the TN analysis, by sub-basin. (b) TN accumulation rates for the four depth intervals (0–25 cm, 25–50 cm, 50–75 cm, 75–100 cm). Data points correspond to yearly means, and error bars to standard errors for the yearly means. Trend lines are obtained from multiple linear regression analysis of TN data. (c) Depth patterns of soil TN content in 1980 and 2010 reveal the greatest accumulation in the top 25 cm. For the numbers of soil samples corresponding to each time step, see supplementary table 4.
Download figure:
Standard image High-resolution image

Samples were selected for analysis based on the following criteria: availability of (1) TN data to a depth of at least 25 cm; (2) soil texture data, including percentages of clay, sand and silt; and (3) latitude and longitude data. Only samples sites falling on land classified as cropland were included in the analysis, as confirmed using United States Geological Society land-use data sets [65]. Based on this criteria, a total of 2069 samples were available at the 0–25 cm depth, 1759 samples for the 25–50 cm depth, 1505 samples for 50–75 cm, and 1320 samples with complete data from 0 to 100 cm. Trend analysis was carried out at each depth range with all the available samples for that range, and also over the entire 100 cm depth using the subset of 1320 samples.

Multiple linear regression (MLR) was performed to account for the impact of multiple explanatory variables (e.g. climate and soil texture) on the observed trends in TN [66]. See supplement 3 for further description of the MLR analysis.

2.3. Modeling framework and illustrative case study

We developed a parsimonious model to describe decadal-scale changes in SON following the initial conversion of grassland or forested land to agriculture, and then its trajectory under intensive agriculture.

**Figure 4.** The fate of anthropogenic N inputs across the MRB. The figure shows a watershed-scale mass balance for the MRB calculated based on data from Goolsby *et al* [70]. The letters correspond to fluxes represented schematically in figure 1. Riverine N output (h) from the Mississippi accounts for approximately 23% of net N inputs. The present study indicates that legacy N accumulation (g) within agricultural soils may account for as much as 3.8 ± 1.8 Mt yr⁻¹ (approximately 53% ± 25% of net N inputs). Although direct measurements of other fluxes are scarce, recent measurement data from the US corn belt suggest an annual nitrous oxide (N₂O) flux (f) for the MRB river network of 0.1 ± 0.01 Mt yr⁻¹ (∼1% of net inputs) [69]. Denitrification to N₂ (f) likely represents a much larger portion of the budget, but the magnitudes remain largely unconstrained [38]. Modeled estimates of sediment burial (l) in reservoirs across the MRB suggest an additional N sink on the order of 0.6 Mt yr⁻¹ (8% of net inputs) [70, 71].
Download figure:
Standard image High-resolution image

We considered the mass of SON in the landscape M(t) (kg ha⁻¹) to be made up of two pools, an active pool ${M}_{{\rm{act}}}$ (kg ha⁻¹) subject to mineralization or immobilization, and a protected pool ${M}_{{\rm{prot}}}$ (kg ha⁻¹) which, when conditions controlling physical and chemical protection mechanisms remain stationary [56], persists in a steady state, with no net mineralization or immobilization. Using this framework, the time (t) evolution of the SON pool is expressed as:

$\begin{eqnarray}M(t) & = & {M}_{{\rm{prot}}}+\left[\displaystyle \frac{\lambda }{k}(t-1)+\displaystyle \frac{{a}_{0}}{k}\right]\\ & & +\left({M}_{{\rm{act}}\_0}+\displaystyle \frac{\lambda }{k}-\displaystyle \frac{{a}_{o}}{k}\right){{\rm{e}}}^{-kt},\end{eqnarray} \tag{ 1 }$

where ${M}_{{\rm{act}}\_0}$ is the initial mass of the active SON pool, ${a}_{0}$ the initial net N input, $\lambda$ (kg ha⁻¹) is the rate of increase in the net N inputs, and k is the mineralization rate constant (yr⁻¹) (details of the derivation provided in supplement 4). Net N inputs are the difference between TN inputs (fertilizer N, atmospheric N deposition, biological N fixation) and N outputs via crop uptake. As described below, different phases of the landscape's evolution are characterized by different values of ${a}_{0},$ ${M}_{{\rm{act}}\_0}$ and $\lambda$ .

**Figure 5.** Modeling framework describing N depletion following conversion of native prairie/grassland to conventional agriculture, and N accumulation following agricultural intensification. On the top left (a) is a model schematic representing the flow of N through the active pool of soil organic N. In (b), we show the depletion of SON following initial cultivation of long-term grassland [54], from which we derive the mineralization reaction rate, k, used in equation 1. The figure on the bottom (c) shows the evolution of the protected and active N pools following land-use change according to the five phases described in the text. N inputs corresponding to the five phases are shown in the inset.
Download figure:
Standard image High-resolution image

Using the above framework, we used Rooks County, Kansas as a case study to explore dynamics in SON depletion and accumulation before and after cultivation and under different management regimes. Rooks County was selected due to its location within the MRB, its long history of cultivation (1870-present), the high proportion of county land maintained under high-input agriculture (50% cropland, wheat/sorghum rotation), and the availability of both pre- and post-cultivation estimates of SON as well as detailed N mass balance data over time (1910–1978) [67]. We modeled five different phases to represent the anthropogenically induced evolution of the landscape: (1) native grassland, pre-cultivation (1840–1890); (2) post-cultivation, low-input agriculture (1890–1910); (3) post-cultivation, low-input agriculture, reduced productivity (1910–1950); (4) post-cultivation, high-input agriculture (increasing inputs) (1950–2000); and (5) post-cultivation, high-input agriculture (stabilized input levels) (2000–2010). Rationales for the parameters in the different periods are provided in supplement 4. Our objective in developing the model was to provide an illustrative tool for exploring the potential for legacy N accumulation under intensive agriculture. Rigorous calibration and validation of the model requires additional site-specific input data that is beyond the scope of this paper.

3. Results

3.1. Changes in Soil N Stocks

3.1.1. Resampling Studies in Iowa and Illinois

For the Iowa resampling study, our results show a net increase in TN of 1478 ± 547 kg ha⁻¹ over the 0–100 cm study depth. The TN content in the surface layer (0–25 cm) was found to decrease slightly, from 2140 ± 60 g N Mg⁻¹ soil to 2110 ± 70 g N Mg⁻¹ soil, although the difference was not significant (Wilcoxon signed rank test, p = 0.162) (figure 2, table 1). At greater depths, however, significant increases were observed. As shown in table 1, the TN content increased by 22% from 25 to 50 cm, by 20% from 50 to 75 cm and by 14% from 75 to 100 cm (p < 0.001, p = 0.013, p = 0.040). Assuming a constant rate of increase over the study period (1959–2007), the above result suggests a yearly accumulation rate of 30.8 ± 11.4 kg ha⁻¹ yr⁻¹.

Table 1. Historical and current magnitudes of soil TN content based on resampling sites originally sampled in the mid 1900s. Sixty-one Iowa sites were first sampled at a median date of 1959 and then resampled in 2007. A significant change in the TN content of the soils is evident for the Iowa study, particularly from 25 to 50 cm. Positive values indicate accumulation. The six Illinois sites were sampled first in 1957 and then again in 2002. Increases in the soil TN content were also observed in the Illinois study; the increases were not significant, however, due to the smaller sample size.

		Historical	Current	Number of samples	Accumulation/Depletion
Location	Depth	g N Mg⁻¹ soil	g N Mg⁻¹ soil	n	g N Mg⁻¹ soil (kg N ha⁻¹)	p-value
	0–25 cm	2140 ± 60	2110 ± 70	61	−30 ± 90	0.162
					(−110 ± 331)
	25–50 cm	1060 ± 40	1290 ± 70	61	230 ± 80	<0.001
					(864 ± 300)
Iowa Study (1959–2007)	50–75 cm	610 ± 40	740 ± 50	25	120 ± 60	0.013
					(474 ± 237)
	75–100 cm	440 ± 40	500 ± 30	12	60 ± 50	0.040
					(250 ± 209)
	0–100 cm	1063 ± 23	1160 ± 29	—	95 ± 36	—
					(1478 ± 547)
	0–20	2733 ± 176	2583 ± 119	6	−150 ± 213	0.516
					(−17 ± 443)
Illinois Study (1957–2002)	20–50 cm	1088 ± 127	1387 ± 164	6	298 ± 77	0.140
					(1746 ± 689)
	50–100 cm	297 ± 29	492 ± 38	6	195 ± 12	0.016
					(1436 ± 235)
	0–100 cm	1022 ± 54	1179 ± 58	—	157 ± 49
					(3164 ± 1033)

The Illinois resampling results demonstrate a 16% net increase in TN, or 3164 ± 450 kg ha⁻¹ averaged over the 0–100 cm depth, between 1957 and 2002. Similar to the Iowa study, an insignificant (5%) decrease in TN was observed in the surface layer (0–20 cm) (p = 0.516) (figure 2, table 1), while TN increased from 20 to 50 cm (27%) (p = 0.140) and from 50 to 100 cm (66%) (p = 0.016). Again assuming a constant increase in TN content over this time period, the total increase corresponds to a yearly rate of 70.3 ± 10.0 kg ha⁻¹ yr⁻¹. Despite the small sample size for the Illinois study (n = 6), these findings are significant (p = 0.016) from 50 to 100 cm and are suggestive of potential increases in TN at a decadal scale in soils under high-input agriculture.

3.1.2. Trend analysis of soil TN data across the MRB

Data from 2069 NCSS soil samples [64] obtained from all six sub-basins of the MRB (figure 3(a)) between 1980 and 2010 was utilized to identify possible trends in the TN content of MRB agricultural soils. Results of the MLR analysis indicate significant increases in soil TN concentrations (g N Mg⁻¹) between 1980 and 2010 in the 0–25, 25–50 and 50–75 cm layers (13.2 g N Mg⁻¹, p < 0.001; 7.3 g N Mg⁻¹, p < 0.001; 3.9 g N Mg⁻¹, p = 0.003) (table 2, figure 3(b)). An increase (1.4 g N Mg⁻¹ yr⁻¹) was also seen from 75 to 100 cm, although the difference was not significant (p = 0.294). Over the entire depth range, using data only from pedons sampled to a depth of 100 cm, the accumulation rate is $3.4\quad \pm \quad 1.6\,{\rm g}$ N Mg⁻¹ yr⁻¹ (p = 0.003). Based on reported bulk density values, these results correspond to total increases (0–100 cm) of 54.8 ± 25.8 kg ha⁻¹ yr⁻¹.

Table 2. Accumulation rates for TN in soil samples across the Mississippi Basin (1980–2010) based on MLR analysis of the NCSS dataset. All available samples at each depth range were used to calculate the depth-specific accumulation rates. Overall accumulation rates (0–100 cm) are calculated not simply as the mean of the four smaller depth increments, but as part of a separate analysis in which only pedons with complete sampling data to 100 cm were considered. As the thickness of the soil profile can vary significantly, and because organic matter may accumulate preferentially in the upper layers of shallower soils (69), we use the more conservative estimate of accumulation suggested by the integrated analysis for the 0–100 cm depth range in subsequent discussions of estimated accumulation rates across the MRB. See supplementary table 7 for results by sublayer for the 1320-sample subset, and a discussion of the differences in the two estimation methodologies.

	Depth	Number	Bulk density	Rate of change
Soil parameter	(cm)	(n)	(g cm⁻³)	(g Mg⁻¹ yr⁻¹)	(kg ha⁻¹ yr⁻¹)	p-value
	0–25	1320	1.55	10.3 ± 3.0	39.9 ± 11.6	<0.001
	25–50	1320	1.61	1.9 ± 2.0	7.65 ± 7.3	0.354
Total nitrogen	50–75	1320	1.64	0.8 ± 1.6	3.4 ± 6.6	0.614
	75–100	1320	1.65	1.6 ± 1.4	6.6 ± 5.8	0.250
	0–100	1320	1.61	3.4 ± 1.6	54.8 ± 25.80	0.003

^aMass-per-area accumulation rates (0–100 cm, kg ha⁻¹ yr⁻¹) are calculated using the mass-per-mass accumulation rates (0–100 cm, g Mg⁻¹ yr⁻¹) and the corresponding bulk density.

4. Discussion

4.1. Quantifying legacy: synthesis of mass balance and soil sampling results

The first question posed in the introduction was whether N is accumulating in agricultural soils, and if so, in what magnitude. The results from the Iowa resampling study, based on data obtained from 61 sites across Iowa, show a 9% increase in TN and suggest an accumulation rate of 30.8 ± 11.4 kg ha⁻¹ yr⁻¹ from 0 to 100 cm. A somewhat larger percent increase (15%) was seen from 0 to 100 cm at the Illinois resampling sites by David et al [50], corresponding to an accumulation rate of 70.3 ± 18.4 kg ha⁻¹ yr⁻¹. Furthermore, our analysis of 2069 soil samples in the MRB demonstrates a 10% increase in soil TN from 0 to 100 cm between 1980 and 2010, corresponding to an accumulation rate of $54.8\quad \pm \quad 25.8$ kg ha⁻¹ yr⁻¹ in cropland soil and an overall accumulation magnitude of 142 Tg N over the MRB over the 30 year period (figure 3(c)). While other studies have alluded to the possibility of N_R accumulating within the soil profile based on mass balance or modeling-based estimates [30, 45, 43], our study for the first time, provides direct, large-scale evidence of such accumulation.

We next explored the relationship between these accumulation magnitudes and estimates of N fluxes in MRB to assess the significance of these magnitudes at the basin scale. We have calculated watershed-scale net N inputs for the years 1980–1996 in the MRB to be 7.1 Mt yr⁻¹ based on data reported by Goolsby et al [68]. During this period, the riverine flux of nitrate from the MRB to the Gulf of Mexico is estimated to have been $1.6\quad \pm \quad 0.1\;$ Mt yr⁻¹ [68], which constitutes approximately 23% of net N inputs. In this context, our estimate of soil N accumulation across the MRB (3.8 ± 1.8 Mt yr⁻¹) suggests that soil N accumulation could account for another 53% ± 25% of net N inputs (figure 4). While significant uncertainty remains regarding the actual magnitude of this estimate of N accumulation, the present results strongly suggest that changes in soil N stocks constitute a significant fraction of TN inputs under intensive agriculture and thus should be explicitly considered in watershed as well as regional and global-scale N mass balance studies.

4.2. Understanding legacy: a conceptual model to explain N depletion and accumulation dynamics

The importance of agricultural soil as an N-sink, as described above, leads us next to question the mechanism behind such subsurface N accumulation. We hypothesize that such accumulation is a direct result of increased N fertilizer use (inorganic and manure N), increases in N fixation due to dramatic increases in soybean cultivation between 1940 and the present, and the adoption of conservation tillage practices [72, 52]. Accordingly, we can utilize the parsimonious modeling framework introduced in section 2.4 to describe not only the depletion of SON following the initial conversion of grassland or forested land to agriculture, but also the accumulation of N suggested by our analysis of soil data from the MRB (figure 5).

In the pre-cultivation period (Phase I: 1840–1890), SON is assumed to be at steady state, with most of the organic nitrogen in the protected pool (figure 5). The start of cultivation (Phase II: 1890–1910) leads to conversion of a portion of the protected SON to active SON, which can then be mineralized and leached from the soil profile. Net N inputs are negative in this period due to intensive cropping practices, but little input of fertilizer [67]. With these changes, there is an exponential decrease in the total mass of SON, with the system eventually evolving to a new steady state (figure 5). After the first 20 years of cultivation ( $\mathrm{Phase}\;\mathrm{III}:\;1910\;-\;1950),$ we assume crop productivity to be diminished which leads to an increase in the net N inputs, and a stabilization of soil N levels (figure 5). Finally in Phase IV (1950–2000), the system transitions to a high-input state and soil N levels begin to rise. Then, at the start of Phase V, with the stabilization of net N inputs, soil N levels also stabilize.

The accumulated or legacy N, conceptualized as the difference between the Phases III and V steady states, corresponds to an approximately 9% increase over SON levels in the depleted steady state. This value is similar to the 10% increase observed over time in our MRB soil sampling data and the 9% increase observed in the Iowa resampling study. The modeling results suggest not only that soil N accumulation is possible in land under continuous cultivation, but that the trajectory of change can be described using the same approach as that used to describe the more well-known depletion of soil organic matter after initial cultivation.

Of course, the above simulation is based on a simplification of the trajectories of change in management practices and land use. We have assumed step changes between phases, but in reality changes occur more gradually. In its current state, the model described herein is primarily conceptual in nature, used to demonstrate the possibility of N accumulation in agricultural landscapes, and will require further modifications in terms of model parameterization as well as descriptions of inputs and outputs to more fully simulate landscape-scale changes in SON.

4.3. Implications of legacy: time lags in landscape response

The most significant implication of such a buildup of soil N relates to time lags observed between land-use changes and alterations in stream N concentrations [73]. Based on the current results, we contend that there are two components of this time lag attributable to two different types of legacy: a hydrologic legacy and a biogeochemical legacy. The hydrologic legacy corresponds to dissolved N in groundwater reservoirs and unsaturated zones, and its existence contributes to the hydrologic time lag—defined as the average time required for dissolved N species to move from the point of application to the point of concern. The existence of the hydrologic time lag is well accepted, with a variety of hydrogeologic controls having been found to result in travel times ranging from days to decades [9, 74, 75]. The second type of legacy, the biogeochemical N legacy, arises from retention of N within the root zone, likely in organic form, and constitutes a long-term source for mineralization and ${{{\rm{NO}}}_{3}}^{-}$ leaching. The existence of such a biogeochemical legacy for phosphorus (P) is well known due to its reactive properties, and legacy sorbed P accumulation has been reported in both soil and sediments [9, 76]. The possibility of such a biogeochemical legacy for N, however, has been mostly neglected, as N in the form of ${{{\rm{NO}}}_{3}}^{-}$ is non-sorbing and is easily leached from soils [9].

The magnitude of the associated biogeochemical time lag is a function of not only the mass of TN accumulation, which has been the focus of this paper, but also the rates of organic N mineralization and the loss of dissolved N through the different biogeochemical and hydrologic pathways. Further research is needed to clarify these mechanisms and pathways. However, as a first estimate, we can utilize the modeling framework developed in section 4.2 to determine the time lag associated with depletion of the 142 Tg of legacy N suggested by our analysis of MRB soil data. Assuming a complete cessation of agricultural production in the region and a return of net annual inputs (a) to the pre-cultivation levels of 5 kg ha⁻¹ yr⁻¹, our model results indicate a biogeochemical time lag of 35 years for 99% depletion of the legacy N. The total lag time would then be a function of both the biogeochemical and the hydrologic lag time, and the latter in itself can be on the order of decades depending on the sizes of saturated and unsaturated zone reservoirs. With such long time frames for recovery, it is thus critical to understand both the accumulation and the ultimate fates of these significant stores of subsurface N for sustainable management practices in large-scale agroecosystems such as the MRB.

4.4. Intersecting lines of evidence

Understanding the long-term dynamics of N in agricultural soils is complex due to the poorly constrained fluxes of denitrification, mineralization and immobilization over varying spatio-temporal scales [18, 19]. However, recent research, as described below, provides intersecting lines of evidence that point towards the accumulation of legacy N in the soil profile in much larger magnitudes than previously conceptualized. For example, using a combination of mass balance and process based modeling, the United States Environmental Protection Agency estimates cropland N accumulation in US to be equal to 17% of fertilizer N inputs [46], while accumulation in Canada has been estimated to be equal to 19% of TN inputs [30]. Accumulation is also suggested by isotope tracer studies that show a 15% retention of ¹⁵N-labeled ${{{\rm{NO}}}_{3}}^{-}$ fertilizer within the soil profile nearly 30 years after application, implying that N fertilizer has a significant residence time in the SON pool [77]. The existence of legacy N is further corroborated by observations of biogeochemical stationarity for N in landscapes under intensive agriculture [78, 79]. In such landscapes, the supply of N to surface waters appears to be transport-limited rather than source-limited, suggesting that the existence of legacy N within the landscape provides an ongoing N source and therefore a positive, linear correlation between riverine N flux and discharge, with N concentrations remaining relatively invariant. This behavior is in contrast to that observed in more pristine landscapes, where N concentrations vary in time in response to source limitations [80].

Our work makes a unique contribution to this literature by providing the first measured estimate of large-scale N accumulation in soils across the MRB. The intersection of such varying lines of evidence, both direct and circumstantial, leading to estimates of soil N accumulation, suggests (a) that we must acknowledge the existence of a growing pool of SON in agricultural landscapes and (b) that we must more explicitly explore the impacts of such a pool on future water quality.

5. Conclusion

Our study has three fundamental contributions. First, our finding of significant N accumulation in agricultural soils across the MRB (3.8 ± 1.8 Mt yr⁻¹) makes a critical contribution towards clarifying the fate of the 'missing' N that is consistently referred to in reports of watershed-scale mass balance studies [21, 75–77]. Although caution must be exercised in relying upon the precise magnitude of accumulation due to large uncertainties in the data, by identifying a clear possibility of significant N accumulation within agricultural soils we make considerable progress towards the closing of N budgets, from the watershed to the global scale. Second, we have developed a simple model that describes both the accumulation and depletion dynamics of SON arising from anthropogenic perturbations on the landscape, thus confirming our hypothesis that the same underlying mechanism can be used to describe both N depletion following plowing and N accumulation as a result of high-input agriculture. The third contribution of this study is with respect to time lags between best management practices and water quality benefits. The significant mass of organic N accumulating in agricultural soils implies that stream N concentrations will persist for decades after fertilizer inputs have ceased. Indeed, the time lag would in most cases be significantly greater than that estimated based on the hydrologic legacy alone. Our study for the first time links multiple lines of evidence to show convincingly that N, like P, has a biogeochemical legacy, a legacy that complicates our previous understanding of the fate of this nutrient in anthropogenic landscapes and that must be accounted for in intervention efforts to improve water quality.

Acknowledgments

Work done by N B Basu and K J Van Meter was supported, in part, by funds from the National Science Foundation Coupled Natural and Human Systems program, Grant Number 1114978. Financial support for N B Basu and K J Van Meter was also provided from startup funds at the University of Iowa and University of Waterloo. Additional funds have come from an NSERC Discovery Grant. We thank Suresh Rao of Purdue University and James Jawitz of the University of Florida for their critical feedback.

Author Contributions

NBB and KJV conceived and developed the core ideas for the paper and contributed equally to this work. JV and CLB provided the Iowa soil sampling data and offered suggestions regarding the soil datasets. NBB and KJV synthesized data from the Illinois and Iowa sampling studies and from the NCSS database and also developed the modeling framework.