Catchment legacies and time lags: a parsimonious watershed model to predict the effects of legacy storage on nitrogen export.

Nutrient legacies in anthropogenic landscapes, accumulated over decades of fertilizer application, lead to time lags between implementation of conservation measures and improvements in water quality. Quantification of such time lags has remained difficult, however, due to an incomplete understanding of controls on nutrient depletion trajectories after changes in land-use or management practices. In this study, we have developed a parsimonious watershed model for quantifying catchment-scale time lags based on both soil nutrient accumulations (biogeochemical legacy) and groundwater travel time distributions (hydrologic legacy). The model accurately predicted the time lags observed in an Iowa watershed that had undergone a 41% conversion of area from row crop to native prairie. We explored the time scales of change for stream nutrient concentrations as a function of both natural and anthropogenic controls, from topography to spatial patterns of land-use change. Our results demonstrate that the existence of biogeochemical nutrient legacies increases time lags beyond those due to hydrologic legacy alone. In addition, we show that the maximum concentration reduction benefits vary according to the spatial pattern of intervention, with preferential conversion of land parcels having the shortest catchment-scale travel times providing proportionally greater concentration reductions as well as faster response times. In contrast, a random pattern of conversion results in a 1:1 relationship between percent land conversion and percent concentration reduction, irrespective of denitrification rates within the landscape. Our modeling framework allows for the quantification of tradeoffs between costs associated with implementation of conservation measures and the time needed to see the desired concentration reductions, making it of great value to decision makers regarding optimal implementation of watershed conservation measures.

Introduction

High levels of nonpoint source pollution associated with current agricultural practices have contributed to n class="Chemical">water quality impairmenpan>t anpan>d destructionpan> of aquatic ecosystem habitats at both local anpan>d global scales [1,2]. Inpan> particular, inpan>pan> class="Chemical">creased nutrient loads delivered from watersheds due to agricultural intensification, industrialization, and urbanization have led to the persistence of large hypoxic zones in both inland and coastal waters [3-7]. Watershed management practices to target these non-point source pollutants have in many cases resulted in little or no improvement in water quality, even after extensive implementation of conservation measures [8-10]. The lag time between implementation of conservation measures and resultant water quality benefits has recently been recognized as an important factor in their “apparent” failure [8,11]. Conservation measures are often implemented, however, without explicit consideration of such lag times, and with the expectation that they will lead to immediate benefits. Failure to meet such expectations then discourages vital restoration efforts [8]. In order to address this problem, it is important to quantify the lag times associated with watershed management efforts a priori and to implement restoration strategies that are targeted specifically at minimizing lag times as well as maximizing restoration benefits.

The focus of the present research is to develop a framework for understanding the time lags between land-use change or implementation of conservation measures and stream n class="Chemical">water quality benpan>efits. We hypothesize that such time lags arise from legacies that have accumulated inpan> the lanpan>dscape over decades of pan> class="Species">human impact [12-14]. Legacies can be conceptualized as hydrologic legacy, in the form of dissolved solute that is delayed in its transport to the stream due to slow groundwater transport pathways, and biogeochemical legacy, arising from solute that has undergone biogeochemical transformation and that is retained within the soil matrix. Both solute and watershed attributes define whether such legacy sources will be created, and, if created, their spatial location within the watershed.

In the present study, we focus specifically on the fate of anthropogenic n class="Chemical">nitrogen (pan> class="Chemical">N) in predominantly agricultural watersheds. Significant time lags between land-use change and the expected decreases in stream nitrate concentrations have consistently been noted [8]. Such time lags have chiefly been attributed to what we have defined as the hydrologic N legacy, a legacy existing primarily in groundwater reservoirs or thick unsaturated zones in the form of dissolved nitrate [9,15,16]. Recent work, however, suggests that consideration of this hydrologic legacy alone does not adequately account for the magnitude of legacy N existing within intensively managed landscapes. Watershed-scale mass balance studies, for example, consistently indicate the presence of "missing" N stores [17-24], and recent modeling of the global N cycle has led to estimates of terrestrial N sequestration ranging from 27–100 Tg N/yr [25-27]. At the plot scale, a recent isotopic tracer study designed to investigate the long-term fate of nitrate fertilizer has shown that approximately 15% of fertilizer N applied to agricultural land is present within the soil profile in organic form 30 years after its initial application [28]. In another study [29], isotopic data were used to demonstrate that the nitrate measured in streams is generated from organic nitrogen created from fertilizer applied to the landscape decades previously. These results are indicative of high levels of N retention within agricultural soil over a multi-year period, and thus the existence of a biogeochemical N legacy, which is corroborated by our recent research showing a basin-wide accumulation of organic N in the Mississippi River Basin [30].

Despite such studies demonstrating the existence of both hydrologic and biogeochemical n class="Chemical">N legacies, most mechanpan>istic watershed models lack anpan> explicit mechanpan>ism to describe the effects of these legacies on stream nitrate concentrations [8,31,32]. Most lumped watershed models such as SPARROW and GlobalNEWS as well as the Net Anthropogenic Nitrogen Inputs (NANI) mass balance approach assume the N cycle to be at a steady state, either on a yearly basis or over a multi-year period, such that stream export is a fixed percentage of net annual inputs [33-37]. Even mechanistic attempt at quantify the benefits of different pollution-reduction scenarios (e.g. land-use change, reductions in fertilizer application, etc.) using mechanistic models such as SWAT have focused only on the concentration reduction benefit that will be achieved at infinite time, with no consideration of the time that will be required to achieve that goal [38,39]. Such consideration, however, is critical for watershed managers who must make decisions regarding the allocation of limited resources for conservation [8].

In this paper, we develop a parsimonious analytical model to quantify the concentration reduction benefits associated with n class="Chemical">watershed restorationpan> efforts as a funpan>ctionpan> of both hydrologic anpan>d biogeochemical legacies, with particular attenpan>tionpan> beinpan>g paid to the ways inpan> which spatial patternpan>s of lanpan>dscape conpan>versionpan> impact conpan>cenpan>trationpan> reductionpan> scenpan>arios. Conpan>cenpan>trationpan> reductionpan>s are conpan>sidered to occur as a funpan>ctionpan> of both the grounpan>dpan> class="Chemical">water travel time distribution and biogeochemical controls, including the existence of a biogeochemical N legacy within the soil profile and denitrification dynamics along groundwater pathways. The paper presents analytical relationships between: (a) percent reductions in mean concentrations at the watershed outlet as a function of the fractional watershed area over which the management practice is implemented, and (b) the temporal trajectory of watershed response that defines the time required to achieve required reductions in contaminant concentrations. Using these analytical relationships, we explore idealized scenarios of land-use conversion and compare these results with concentration trajectories observed in a small Midwestern watershed undergoing an extensive prairie habitat restoration project. We further use these relationships to establish an optimization framework for meeting concentration reduction goals by exploring tradeoffs between costs associated with the conversion of land out of row-crop production and the time required to achieve the desired concentrations. Such explorations enable analysis of the performance of conservation measures under spatially varying patterns of intervention as a function of legacy N accumulation, N removal dynamics in the subsurface, and watershed travel time distributions.

Model Development

Our conceptual framework is based on the assumption that legacy nutrient stores are present within anthropogenic landscapes and lead to time lags between land-use change and improvements in n class="Chemical">water quality. Such nutrienpan>t legacies have developed inpan> agricultural pan> class="Chemical">watersheds as a function of long-term application of N and phosphorus (P) fertilizers, with a strong linear correlation having been found between N and P levels in soils and multi-decadal cumulative nutrient surpluses [13,30,40]. Our focus herein is specifically on the N legacy in agricultural watersheds, but this approach could be readily adapted to other solutes.

As shown in Fig 1, nutrient legacies produce an internal landscape memory, thus contributing to elevated stream nutrient concentrations for years after external nutrient loading is reduced or stopped altogether. In order to develop an expression for the concentration trajectory at the catchment outlet following land-use change, we conceptualize the landscape to be composed of a bundle of stream tubes, with each point on the landscape being associated with an individual stream tube characterized by a specific groundn class="Chemical">water travel time to the catchmenpan>t outlet [41]. The full amalgamationpan> of poinpan>ts for the catchmenpan>t leads to a specific grounpan>dpan> class="Chemical">water travel time distribution for the catchment outlet, f(τ) This distribution, in turn, controls the concentration trajectory at the outlet, C(t) [42-44], as described by the following equation: where, C (t—τ) is the contaminant input function or “source function” from the unsaturated zone, and k [T-1] is the first-order rate constant that describes removal processes in the aquifer.

Fig 1

Conceptual framework for predicting catchment scale time lags as a function of hydrologic and biogeochemical legacies in the landscape.

The left frame represents depletion of biogeochemical legacy in the source zone. The source zone depletion function is then convoluted with the groundwater travel time distribution (middle frame) to ultimately describe concentrations at the catchment outlet (right frame).

The left frame represents depletion of biogeochemical legacy in the source zone. The source zone depletion function is then convoluted with the groundn class="Chemical">water travel time distribution (middle frame) to ultimately desn class="Chemical">cribe concentrations at the catchment outlet (right frame).

The source function (Fig 1), developed in Section 2.1, is controlled by the biogeochemical legacy in the unsaturated zone, which for n class="Chemical">N is a funpan>ctionpan> of both historic anpan>thropogenpan>ic pan> class="Chemical">N inputs to the landscape and the rate of N depletion from such stores. Each point in the watershed is characterized by its particular source function, which changes form as a function of the timing of human interventions such as land-use change or implementation of conservation measures. While biogeochemical legacy is conceptualized using the “source function,” the hydrologic legacy is captured in the travel time distribution, which describes how the source concentrations are being modified as they travel through the watershed (Section 2.2). The resulting outlet concentration is a function of both the hydrologic and biogeochemical legacy, and the patterns of land-use change, as described in the following sections.

2.1 Biogeochemical Legacy and the Source Function

The left frame of Fig 1 provides a simple schematic for our model of biogeochemical legacy depletion within the source zone after conversion from row-n class="Gene">crop agriculture to grasslanpan>d. Withinpan> this framework, excess soil pan> class="Chemical">organic N, which has accumulated in response to long-term application of fertilizer N and which constitutes the biogeochemical N legacy, is mineralized to inorganic N, entering the soil mineral N pool. This inorganic N then leaches to groundwater, primarily in the form of nitrate. Although plant uptake and litter inputs will continue to occur after conversion to grassland, we consider these processes to be part of a baseline scenario and therefore only take into account dissipation of excess N through the leaching pathway. In our current simulations, we consider only scenarios in which landscape conversion results in a complete cessation of fertilizer application to the soil system, although this formulation can be easily modified to include cases with ongoing but reduced levels of fertilizer application.

Decomposition of soil organic matter is typically modeled as having first-order reaction kinetin class="Chemical">cs, proportionpan>al to the amounpan>t of substrate to be decomposed [45,46]. Accordinpan>gly, withinpan> our modelinpan>g framework the excess (legacy) soil pan> class="Chemical">organic N (SON) is considered to decay over time as a first-order process with a rate constant λ (T-1), such that the mass of legacy N ( ) at any time t is given by: The N leaving the SON pool enters the mineral N pool that leaches into the groundwater, such that at any point in time the concentration of dissolved N (Cs; M/L3) can be described by the following equation: where, Q is the mean annual recharge and Vw (= nsV) is the water volume in the source zone, with n being the porosity, s the saturation and V the volume of the soil column per unit area within the source zone. Here, the first term on the RHS is the input from the organic pool and the second term is the loss from the source zone via leaching. Solving Eqs 2 and 3 leads to: where ;; is the initial mass of SON; is the initial concentration of nitrate within the source zone; and is the nitrate concentration at time t within the source zone, and acts as the source function described in Eq 1. The source function can thus be described as a Heaviside function, with being the initial steady-state concentration prior to initiation of land-use change, and describing the concentration trajectory after land-use change has been initiated.

2.2 Hydrologic Legacy and Patterns of Land-use Change

We define the hydrologic nutrient legacy as nutrients present in a dissolved form in both the saturated and unsaturated zones. Time lags associated with hydrologic legacy can range from days to weeks to hundreds of years as a function of the distance groundn class="Chemical">water must travel to the catchment outlet, the physical properties of the unpan>derlying aquifer, and the gradient driving flow through the subsurface [9,47].

Land-use change in a n class="Chemical">watershed leads to switchinpan>g of the source funpan>ctionpan> betweenpan> anpan>d . Theoretically, inpan>dividual poinpan>ts inpan> the lanpan>dscape may be switched at differenpan>t poinpan>ts inpan> time, or not switched at all, leadinpan>g to anpan> inpan>finpan>ite number of scenpan>arios that are conpan>voluted as inpan> Eq 1, pan> class="Chemical">creating unique concentration trajectories at the outlet. Here, we conceptualize three end-member scenarios based on the distribution of travel times for the watershed. As shown in Fig 2, spatial patterns of land-use change can be described as truncations of the groundwater travel time distribution, with the three scenarios of change being: (a) frontal, (b) distal and (c) random. The frontal approach corresponds to scenarios where land parcels with the shortest travels times to the catchment outlet are preferentially converted and involves a left-to-right truncation of the exponential travel time distribution (Fig 2a), as indicated by the grey shaded area of the figure. Conversely, a distal approach corresponds to a preferential conversion of parcels with the longest travel times to the outlet and involves a right-to-left truncation (Fig 2b). The third approach, a random conversion, corresponds to a scenario where land-use change has occurred randomly throughout the catchment, with no correlation between land-use change and the groundwater travel time distribution (Fig 2c).

Fig 2

Conceptual framework showing spatial patterns of land-use change as truncations of the groundwater travel time distribution.

The grey shaded areas correspond to the fractional areas (p) of the watershed over which land-use change has occurred for the (a) frontal, (b) distal and (c) random conversion scenarios. The red line in the frontal and distal scenarios is equal to the abscissa of the cumulative frequency distribution of travel times, corresponding to an ordinate of p (or 1-p for the distal scenario), and is the demarcation line between areas that have and that have not undergone land-use change.

Equations for the flow-averaged concentrations at the outlet at any time t after initiation of land-use change, C (t;p) normalized by the concentration before the change C (t;p) for the three different spatial conversion scenarios, and a fractional land-use change p can be developed as follows: Frontal Distal Random Here F -1()is the inverse cumulative distribution function of the groundn class="Chemical">water travel time distributionpan>. It represenpan>ts the travel time associated with a specific fractionpan>al area of lanpan>dscape conpan>versionpan> anpan>d acts as a dividinpan>g linpan>e (red linpan>e inpan> Fig 2a anpan>d 2b) betweenpan> areas of the pan> class="Chemical">watershed bringing in “converted” groundwater and areas bringing in “unconverted” groundwater. The above equations have been developed with the assumption that the groundwater travel time distribution is a complete distribution from 0 to infinity. In actuality, however, these distributions would be truncated, with the maximum travel time being defined by the size as well as the geomorphic characteristics of the watershed, and the equations can be easily modified for truncated distributions following Jawitz et al. [48].

Groundn class="Chemical">water travel time distributionpan>s have beenpan> assumed to have multiple funpan>ctionpan>al forms based onpan> a ranpan>ge of model types, from the simplest pistonpan>-flow model, which assumes that all flow paths have the same velocity anpan>d path lenpan>gth, to a dispersionpan> model based onpan> a 1-D solutionpan> of the advectionpan> dispersionpan> equationpan> [44]. Onpan>e of the simplest anpan>d most widely used forms is the exponpan>enpan>tial: where τ is the travel time anpan>d μ is the meanpan> travel time for the pan> class="Chemical">watershed. Here, we have used the exponential distribution to develop algebraic expressions for the flow-averaged concentration after land-use change following the three different patterns of intervention described in Eqs 6,7, and 8.. Frontal Distal Random where Here, F and F represent the longest and shortest groundwater travel times, respectively, for land parcels that have undergone land-use conversion. Note that in the above formulations we have implicitly assumed that the watershed is homogeneous in terms of land use. This assumption is not, however, a limitation of the approach. For heterogeneous land use, the groundwater travel time distribution of interest is that of the areas contributing solute to the watershed outlet. For example, when only a fraction of the watershed area is under row crops and the management practice involves conversion of row crops to prairies, the travel time distribution used in Eq 9 would be that of the cells originally under row crop.

Materials and Methods

3.1 The Walnut Creek Case Study

3.1A Site Description

A n class="Chemical">watershed habitat restorationpan> anpan>d agricultural manpan>agemenpan>t project was implemenpan>ted by the Unpan>ited States Fish anpan>d Wildlife Service (USFWS) at the pan> class="Chemical">Neal Smith National Wildlife Refuge (NSNWR) within the Walnut Creek Watershed (WCW) (52 km2) of Jasper County, Iowa (Fig 3a). The project involved conversion of a large portion of the WCW from row-crop agriculture to native prairie and savanna [49,50]. The NSNWR represents one of the first attempts at agricultural land-use conversion towards ecosystem restoration at the watershed scale, and is one of the few sites where water quality has been monitored both in the groundwater directly below the reconstruction and in surface water at multiple scales within the watershed. The site is thus an ideal choice for testing the applicability of the modeling framework introduced in this paper.

Fig 3

Site Information and Results for the Walnut Creek Case Study.

(a) Subwatershed 5 (7.9 km2) of the Walnut Creek watershed, Jasper County, Iowa; (b) Data points correspond to groundwater nitrate concentrations in 19 monitoring wells across a chronosequence of restorations sites. Biogeochemical Legacy Depletion: Source Zone Nitrate-N Concentration as a function of time since land-use change; (c) Hydrologic and Biogeochemical Legacy Depletion: Data points correspond to mean annual nitrate concentrations measured at the outlet of subwatershed 5 as a function of time since land-use change. The grey shaded area in the figure corresponds to a range of values for the denitrification rate constant (k = 0.24 ± 0.08 y-1).

(a) Subn class="Chemical">watershed 5 (7.9 km2) of the Walnpan>ut pan> class="Chemical">Creek watershed, Jasper County, Iowa; (b) Data points correspond to groundwater nitrate concentrations in 19 monitoring wells across a chronosequence of restorations sites. Biogeochemical Legacy Depletion: Source Zone Nitrate-N Concentration as a function of time since land-use change; (c) Hydrologic and Biogeochemical Legacy Depletion: Data points correspond to mean annual nitrate concentrations measured at the outlet of subwatershed 5 as a function of time since land-use change. The grey shaded area in the figure corresponds to a range of values for the denitrification rate constant (k = 0.24 ± 0.08 y-1).

The WCW is located in the Southern Iowa Drift Plain landscape region of Iowa, which is characterized by steeply rolling hills and well-developed drainage [51]. The climate of the area is humid and continental. Temperatures in the region vary widely, ranging from average maximum values over 20°C between June and September to less than 0°C in December and January. Annual precipitation averages around 850 mm, with maximum rainfall typically occurring in the months of May and June.

In 1990, the land cover in the n class="Chemical">watershed was predominpan>anpan>tly agricultural, with 70% of the area beinpan>g covered by row pan> class="Gene">crops. From 1990 to 2005, row-crop cover throughout the watershed was decreased from 70% to 55% as a part of prairie conversion efforts. In subwatershed WNT5 (7.9 km2), which is our focus herein, the row crop cover was decreased from 77% to 46%, and surface water quality was monitored subsequently over a period of 13 years [49]. Trajectories for groundwater nitrate concentrations throughout the conversion area were established based on water sampling from 19 monitoring wells across a chronosequence of sites, as indicated by the data points in Fig 3b [52]. For the chronosequence work, sites were selected to represent a conversion time series, with three of the sites still under row crop and the rest having been converted 2–13 years prior to the sampling date. Nitrate concentrations were measured at the outlet of subwatershed 5 and were used to calculate mean annual concentrations (Fig 3c). This site thus provided us with an opportunity to test the ability of the model to capture the dynamics in biogeochemical legacy depletion based on groundwater data, and combined hydrologic and biogeochemical legacy depletion based on surface water data.

3.1B Model Parameters

As desn class="Chemical">cribed inpan> Sectionpan> 2.1, our model assumes chanpan>ges inpan> source zonpan>e conpan>cenpan>trationpan>s over time after conpan>versionpan> from row pan> class="Gene">crop to grassland as a function of both the depletion of legacy soil organic nitrogen (SON) and annual recharge rates to groundwater. To model the biogeochemical legacy dynamics within the source zone, legacy SON is considered to exist within the soil profile to a depth of 1 m at a quantity of 100 kg/ha over baseline (pre-agricultural intensification) levels. This value is a conservative estimate based on N accumulation rates of 6 kg/ha-y for a soil depth of 0–100 cm observed across Iowa under intensive agricultural practices over a period of 70 years (1940–2010) [53], with an assumption that approximately 75% of this accumulation would remain protected via both physical and chemical stabilization mechanisms [54] and that the remaining 25% would be in a readily mineralizable form. Initial NO3-N concentrations in the source zone are assumed to be 15 mg/L based on a reported range of 10–20 mg NO3-N/L in tile drainage and groundwater under corn-soybean rotations [55,56].

The groundn class="Chemical">water travel time distributionpan> for the WCW was determinpan>ed usinpan>g a MODFLOW model that was calibrated againpan>st measured grounpan>dpan> class="Chemical">water elevations at 84 monitoring wells within the site [57-59]. A particle-tracking simulation revealed an exponential travel time distributions for the row-cropped area of the WCW (~ 70% of the watershed is row-cropped) [58]. Reported data on prairie plantings [49] and spatial maps of watershed travel times created using the MODFLOW model [59] demonstrated that the pattern of land-use conversion for WNT5 was predominantly random, which is consistent with our general understanding of land-use shifts for restoration being driven more strongly by the availability of land parcels than design of an optimal land-use change scheme for maximization of water quality benefits. We use a denitrification rate constant (k) varying over a range of 0.24 ± 0.08 y-1, which corresponds to a reported range of denitrification rate constants for shallow aquifers with upland surficial geology characterized by glacial outwash and till [60], as is found at the Walnut Creek site. Other parameters used in the model are included in Table 1.

Table 1

Model Parameters for the Walnut Creek Watershed.

Model parameters	Walnut Creek Values
Initial Source Zone Nitrate Concentration	15 mg NO3-N/L
Initial Mass of Legacy SON	100 kg/m2
Legacy N depletion rate constant (λ)	0.16 y-1
Denitrification Rate (k)	0.24 ± 0.08 y-1
Mean Travel Time (μ)	21.6 y
Mean Annual Recharge (Q)	129.5 mm/y
Soil Saturation (s)	0.5
Soil Porosity (n)	0.3
Fractional Land Area Converted	0.41

3.2 Metrics for evaluating concentration reduction benefits

To quantify the concentration reduction benefits achievable at a specified time interval (t) based on a given percent land-use change, we have developed the n class="Chemical">CR metric, defined as: For the special case of concentration reductions at very long times t→∞, thus representing the maximum benefit that can be achieved by land-use conversion, the n class="Chemical">CR metric is used, defined as:

Results and Discussion

The objective of the present study was to develop a framework to quantify catchment-scale time lags based on both biogeochemical and hydrologic nutrient legacies in intensively managed catchments. Our first intent was to develop a set of analytical equations to quantify n class="Chemical">water-quality benpan>efits, takinpan>g inpan>to accounpan>t both soil legacy accumulationpan> anpan>d denpan>itrificationpan> dynpan>amipan> class="Chemical">cs along the groundwater pathway. Our results, based on idealized scenarios of land-use conversion, are compared with results related to actual patterns of land-use conversion in the Walnut Creek watershed. Additionally, our intent was to utilize the analytical equations to explore concentration-reduction benefits associated with different spatial patterns of land-use conversion, and thus to further our understanding of both natural and anthropogenic controls on such benefits and any associated time lags. Benefits are gauged in terms of (1) the relative magnitude of the watershed chemical response to the cropland conversion and (2) the arrival time of the response at the outlet. These results are used to establish an optimization framework that clarifies tradeoffs between the land area taken out of row-crop production and the time required to achieve desired concentration-reduction benefits.

4.1 Model Validation: The Walnut Creek Case Study

We first applied our model, which takes into account both biogeochemical legacy and the groundn class="Chemical">water denpan>itrificationpan> dynpan>amipan> class="Chemical">cs, to the Walnut Creek watershed. The temporal trajectory of source-zone nitrate concentrations C (t) in land parcels that had undergone conversion from row-crop to grassland was modeled using Eqs 4 and 5. A legacy depletion rate constant, λ, of 0.16 per year was able to capture the observed trends in the groundwater chronosequence data described in Section 3.1A and as shown in Fig 3b. The time required to achieve a 50% concentration reduction in the source zone was found to be approximately 5 years, while a 95% concentration reduction in the source zone corresponded to a lag of approximately 19 years.

The concentration trajectory at the catchment outlet was then derived as a function of C (t), the groundwater travel time distribution, the denitrification rate constant, and the patternpan> of land-use change, following Eqs 10, 11 and 12. The exponential travel time distribution derived from the MODFLOW model (mean travel time μ = 21.6 years) was used long with an assumption of a random patternpan> of land-use conversion (see Methods 3.2), to model three different scenarios for trajectories of water-quality change after land conversion for the WNT5 subwatershed (Fig 3c).

The first scenario (S1) parallels the approach used by Schilling et al. [49], assuming the presence of hydrologic legacy but no biogeochemical legacy, with no denitrification occurring along the groundn class="Chemical">water pathway. Accordinpan>gly, there is pan> class="Chemical">nitrate dissolved in groundwater that continues to arrive at the outlet over a defined time period (as a function of the groundwater travel time distribution) after land-use conversion, thus shaping the outlet concentration trajectory. The second scenario (S2) maintains the same assumptions as those utilized in S1, but adds denitrification in the saturated groundwater. Under this scenario, nitrate concentrations decrease as they travel from the source zone to the outlet, as nitrate is reduced to N2 or N2O and leaves the system in a gaseous form. Both S1 and S2 assume that groundwater concentrations in the source zone, beneath the land parcels for which the land-use shift has occurred, drop immediately from Cs to 0, and that the observed concentration trajectory at the outlet is only a function of dynamics along the groundwater flow pathways.

In contrast, the third scenario (S3) takes into account both denitrification and the presence of biogeochemical legacy in the source zone. The C (t) function in the model, as shown in Fig 3b, is convoluted with the groundn class="Chemical">water travel time distributionpan> followinpan>g Eq 1 to estimate the conpan>cenpan>trationpan>s at the catchmenpan>t outlet. As canpan> be seenpan> inpan> Fig 3c, the base case scenpan>ario S1, after year 3, genpan>erally overestimates the conpan>cenpan>trationpan> at the outlet. Conpan>versely, the S2 scenpan>ario, which takes inpan>to accounpan>t denpan>itrificationpan> dynpan>amipan> class="Chemical">cs, consistently underestimates the achieved concentrations, even when considering the full range of possible values for the denitrification rate constant (k). However, with S3’s combined consideration of both denitrification and biogeochemical legacy, there is a close match between predicted concentration reductions and the observed data for WNT5. In particular, S3 captures the time lag between the initial land conversion and the first observed drop in concentrations at year 4.

The model thus provides a parsimonious way of describing the concentration trajectory, both at the parcel in which land use change has occurred, and at the catchment outlet. Although in the present study we used the more computationally intensive MODFLOW/MODPATH approach to estimate the grounpan>dwater travel time distribution, previous work [58] suggests that a simple GIS-based approach, which uses the land surface as a surrogate for the water table, can be used to construct the travel-time distribution. The latter method is based on readily available DEM data and hydraulic conductivity values, which can be obtained at the local scale from soil databases and at larger, regional scales from recently constructed global maps of near-surface permeability [61].

4.2 Concentration Reduction as a Function of Spatial Patterns of Land-Use Change

We next utilized our analytical equations to explore concentration-reduction benefits associated with different spatial patterns of land-use conversion. The temporal trajectories for the outlet concentration after conversion (50% land-use conversion, k = 0.18 ± 0.12 y-1) normalized to the mean concentration before conversion are presented in Fig 4a. The frontal conversion leads to the fastest response and the distal the slowest, with the random somewhere in the middle. For both frontal and random truncation of the groundn class="Chemical">water travel time distributionpan>, partial benpan>efits are immediately realized at the pan> class="Chemical">watershed outlet, but for a distal truncation there is a time lag between the implementation of change and the start of benefit realization. This time lag corresponds to the minimum travel time of the altered land-use parcels, F (1-p), and is a function of both the groundwater travel time distribution characteristics and the fractional land-use change. In the modeled scenario, this time lag is approximately 16 years, whereas in the frontal and random scenarios concentration reductions of approximately 85% and 43%, respectively, have already been achieved at 16 years after conversion (Fig 4a).

Fig 4

Normalized concentration reduction trajectories under different patterns of land-use change.

(a) Normalized concentration trajectories at the catchment outlet plotted as a function of time (years) after land-use change for the frontal, random and distal patterns of conversion; fractional land-use conversion p = 0.5; (b) Concentration reduction fraction at infinite time as a function of land use conversion fraction p. In both figures, k = 0.18 ± 0.12, which corresponds to a range of “moderate” denitrification rates (Tesoriero et al. 2011). Other parameters used are lambda = 0.23 y-1 and μ = 21.6 y. A 1:1 relationship between CRinf and p, with no dependence on the k values is apparent for the random truncation.

(a) n class="Chemical">Normalized conpan>cenpan>trationpan> trajectories at the catchmenpan>t outlet plotted as a funpan>ctionpan> of time (years) after lanpan>d-use chanpan>ge for the fronpan>tal, ranpan>dom anpan>d distal patternpan>s of conpan>versionpan>; fractionpan>al lanpan>d-use conpan>versionpan> p = 0.5; (b) Conpan>cenpan>trationpan> reductionpan> fractionpan> at inpan>finpan>ite time as a funpan>ctionpan> of lanpan>d use conpan>versionpan> fractionpan> p. Inpan> both figures, k = 0.18 ± 0.12, which corresponpan>ds to a ranpan>ge of “moderate” denpan>itrificationpan> rates (Tesoriero et al. 2011). Other parameters used are lambda = 0.23 y-1 anpan>d μ = 21.6 y. A 1:1 relationpan>ship betweenpan> CRinf and p, with no dependence on the k values is apparent for the random truncation.

Importantly, not only the time required to achieve a desired concentration reduction, but also the maximum achievable concentration reduction at infinite time (n class="Chemical">CR ), differs accordinpan>g to the spatial patternpan> of conpan>versionpan>. This spatial dynpan>amic is captured inpan> Fig 4b, inpan> which pan> class="Chemical">CR values are plotted as a function of the fractional land-use conversion, p, for frontal, distal, and random conversion scenarios. As can be seen in the figure, the frontal pattern of intervention provides the greatest maximum concentration reduction benefit at all percentages of landscape conversion, with the greatest difference between the frontal and distal scenarios occurring under the 50% conversion scenario. The relatively low CR values for the distal scenarios demonstrate that even at very moderate denitrification rates (k = 0.18 ± 0.12 y-1), land parcels with relatively greater travel times make very little contribution to stream nitrate concentrations. Accordingly, conversion of those parcels will have virtually no impact on nitrate concentrations at the watershed outlet.

In contrast, it should be noted that a random pattern of intervention provides a 1:1 concentration reduction benefit. In other words, when interventions are applied randomly throughout a n class="Chemical">watershed, a 20% conpan>versionpan> of pan> class="Chemical">watershed area will result in a 20% reduction in concentration. Mathematically, this 1:1 relationship between land-use conversion and water quality benefits arises due to the property by which a probability distribution created by taking a large enough random sample from any frequency distribution will have the same attributes as the original distribution. It should also be noted that though a range of CR values is obtained for the frontal and distal conversion scenarios, based on the range of denitrification rate constants, the values for the random scenario are not a function of this rate constant and therefore do not deviate from the 1:1 relationship between land-use conversion and concentration reductions.

4.3 Concentration reductions as a function of natural and anthropogenic controls

4.3A Concentration Reductions at Infinite Time

In the above sections, we have focused on concentration-reduction benefits corresponding to one, unique travel time distribution (exponential with μ = 26 years). In order to understand how such benefits vary as a function of the travel time distribution, we have plotted contours of maximum concentration reductions (n class="Chemical">CR ) alonpan>g a conpan>tinpan>uum of values for both the meanpan> travel time (μ) anpan>d the fractionpan>al area withinpan> a watershed beinpan>g removed from row-pan> class="Gene">crop production (p) (Fig 5). Three different plots are presented (5a, 5b,5c), corresponding to the frontal, random and distal, patterns of intervention, respectively. For a particular watershed (characterized by its μ value), the CR benefit achieved for a specified fractional land-use change (p) is equal to p for the random truncation scenario (Fig 5b), greater than p for the frontal scenarios (Fig 5a) and less than p for the distal scenarios (Fig 5c). For a particular conversion fraction p, the concentration reduction benefits increase with increasing mean travel times for the frontal truncation scenario, while they decrease for the distal truncation and remain invariant with μ for the random truncation.

Fig 5

Maximum normalized concentration reduction (CRinf) contours plotted as a function of the fractional land-use conversion p and mean watershed travel time μ.

Contours are plotted for the (a) frontal, (b) random and (c) distal truncation scenarios (k = 0.06 y-1, λ = 0.16 y-1).

4.3B Concentration Reductions at Specified Times

The above analysis describes concentration reduction benefits occurring at infinite time. Decisions regarding land management, however, are not made based on hypothetical benefits achieved at infinite times, but on realistic concentration reductions achievable within specific time frames. In Fig 6, we show the concentration reduction profiles achievable 5 years after landscape conversion (CR values) as a function of p and μ.

Fig 6

Normalized concentration reduction contours at t = 5 years (CR5) plotted as a function of the fractional land-use conversion p and mean watershed travel time μ.

Contours are plotted for the (a) frontal, (b) random and (c) distal truncation scenarios (k = 0.06 y-1, λ = 0.16 y-1).

With a frontal conversion (Fig 6a), it can be seen that 5 years after conversion (t = 5), the concentration reductions for a particular n class="Chemical">watershed (characterized by a μ value) inpan>pan> class="Chemical">crease with increases in p up to a point, and then become invariant with p. Beyond this point, further land-use conversion in this watershed results in no additional stream water quality benefits within the 5-year period, as land being converted beyond this threshold point has an associated travel time greater than 5 years, and thus has no impact. As μ increases, this threshold point shifts to the left, implying that the threshold is crossed at lower and lower p values. This trend occurs because for larger watersheds, a 5-year threshold is only a small percentage of its overall area, and thus benefits cease beyond a relatively small value of p. The thicker line connecting the threshold points thus divides the plot are into two zones, one where benefits are still being realized and another for which benefits have ceased, with the line being mathematically denoted by F-1(p) = td.

It is important to understand the existence of this threshold when designing restoration schemes. For example, land managers working in a n class="Chemical">watershed with a meanpan> grounpan>dpan> class="Chemical">water travel time of 15 years may be under pressure to reduce stream NO3 - concentrations by 50% over a period of 5 years. Knowing that proportionally greater benefits will be achieved with a frontal approach to restoration, they begin converting land with the smallest travel times. Fig 6a, however, shows that under these conditions, a concentration reduction of approximately 35% is the greatest benefit that can be achieved within the 5-year period, even if 100% of the land is converted from row-crop production to native prairie.

The results are quite different, however, for the distal and random conversion scenarios. With a distal conversion, the threshold line is the zero-benefit contour (the heavy dark line in Fig 6c), such that to the left of this line, no concentration reduction benefit can be achieved. Compared with the frontal scenario, it can be seen that a distal approach provides much poorer outcomes, requiring a much greater conversion area and/or a much longer time periods to see results. Assuming the same scenario as above, with a mean groundn class="Chemical">water travel time of 15 years, lanpan>d manpan>agers would be forced to conpan>vert approximately 90% of the pan> class="Chemical">watershed from row-crop to prairie to achieve the same 30% concentration reduction that could be achieved with an approximately 20% conversion area under the frontal approach. Finally, with a random approach, as shown in Fig 6b, concentration benefits scale continuously with p and μ. Additionally, as was seen in the CR plots in Fig 5, a random conversion approach provides poorer concentration reduction outcomes than the frontal approach, but better outcomes than the distal approach.

4.3C Time Lags and Tradeoffs

With an understanding of catchment-scale time lags, an optimization approach can be developed to clarify the tradeoffs involved with achieving a specified concentration reduction benefit. In our case, conversion of land in row-n class="Gene">crop agriculture to native prairie canpan> be unpan>derstood withinpan> the framework of two competinpan>g objectives. Objective 1 (O1) is to achieve specified pan> class="Chemical">nitrate concentration reduction goals within a desired time frame. Objective 2 (O2) is to minimize both societal and individual farmer costs associated with implementation of environmental interventions while still meeting concentration reduction goals. Fig 7 provides a visualization of Pareto-optimal fronts for these conflicting objectives, with the contour lines representing progressively greater concentration reduction goals, from 25 to 75% reduction. In the figure, the x-axes correspond to the time in years after land conversion from row-crop to native prairie necessary for the concentration reductions to be realized (O1), while the y-axes correspond to the economic costs of land converted (O2), with the fractional land are converted (p) serving as a proxy for costs incurred. The three columns correspond to the frontal, random and distal approaches to intervention, with each resulting in its own family of optimized values for land conversion and time required to see the specified concentration reduction benefit. Watersheds with different travel time distributions are also represented here, with rows 1, 2 and 3 corresponding to mean travel times of 10, 20 and 50 years, respectively.

Fig 7

Normalized concentration reduction contours at infinite time as a function of the allowable lag time and the fractional land-use conversion.

The three rows represent different watershed mean travel times, while the three columns represent frontal, random and distal patterns of land-use change (k = 0.06 y-1, λ = 0.16 y-1).

As can be seen in the figure, to achieve progressively greater concentration reduction goals, tradeoffs are necessary between the percent land converted and the time to concentration reduction. For example, if a random approach is taken to carrying out land conversion, as is typical in most watersheds, and a 50% concentration reduction is desired, the time required to achieve the desired concentration benefit ranges from approximately 8 to 30 years (μ = 10 years) (Fig 7c). If Objective 1 is prioritized, to give the fastest possible response time, a more than 90% conversion away from row-n class="Gene">crop must be carried out. Conpan>versely, if Objective 2 (cost) is prioritized, to mainpan>tainpan> the maximum lanpan>d inpan> productionpan>, a 50% conpan>versionpan> is required, with the unpan>derstanpan>dinpan>g that there will be a multi-decade time lag betweenpan> conpan>versionpan> anpan>d fully meetinpan>g pan> class="Chemical">CR goals. An optimal compromise position, within the constraints of the random approach, would likely occur somewhere near the midpoint of the contour line, with approximately 70% of area being converted and an 11-year lag time. To further reduce the necessary percent land conversion, and thus to further minimize the economic impact, a frontal approach could be utilized. Although the fastest that a 50% concentration reduction can be achieved with the frontal approach remains at approximately 8 years (Fig 7a), the percent conversion necessary to achieve this reduction within this time period is reduced from 90% to 40%.

Such tradeoffs are also a function of the mean travel time for the watershed. In watersheds of the same size but with different mean travel times, the greater benefits of the frontal approach correlate positively with the travel time, allowing concentration reduction objectives to be achieved with significantly less commitment of resources. For example, in the μ = 10 year watershed, a 50% n class="Chemical">CR requires, at minpan>imum, a close to 40% conpan>versionpan> of lanpan>d area out of row pan> class="Gene">crop. In contrast, in the μ = 50 year watershed the same reduction can be achieved with a 30% conversion over a similar time frame.

In general, it can be seen that the concentration reduction response scales according to both watershed characteristin class="Chemical">cs (mean travel time) and the employed management approach (spatial patternpan>s of intervention), and with such changes the optimal level of intervention can be either more widely or more narrowly defined. If considering only tradeoffs between fractional land-use change (cost) and the time required to achieve target concentrations, the frontal approach provides a more clearly defined optimal intervention, with conversion of additional land area providing little or no additional time advantage beyond a threshold value. In contrast, with a random or more distal approach, the tradeoffs between time and the fractional converted area scale over a wider range of values, thus leading to more room for debate regarding the best path towards achieving concentration reduction goals.

Summary and Implications

In recent years, there has been great interest and investment of both private and public funds in the implementation of conservation-oriented management practices and other measures to minimize the negative environmental impacts of modern agricultural practices. Such interventions range from the retirement of agricultural land through programs such as the U.S. Conservation Reserve Program, to reductions in fertilizer application and the n class="Chemical">creationpan> of riparianpan> buffer zonpan>es. Inpan>terest is also growinpan>g inpan> the potenpan>tial mitigatinpan>g impacts of large-scale conpan>versionpan> from grainpan>-based pan> class="Gene">cropping systems to the cultivation of perennial biofuel crops such as switchgrass and miscanthus, which have been found to result in reduced nitrate leaching at the plot scale [62]. Although numerous studies have attempted to demonstrate the potential water-quality benefits garnered by implementing such changes [38,63], there has been little acknowledgment of the often long time periods required to achieve such benefits [8]. In addition, most existing models such as SWAT and AGNPS [64,65], which are commonly utilized for agricultural landscapes, do not have an explicit mechanism to either account for such legacies or to predict time lags [66].

In the present work, we have developed a framework that allows for the parsimonious modeling of concentration-reduction benefits over time as a function of spatial patterns of land-use conversion or implementation of conservation measures an class="Chemical">cross the lanpan>dscape, anpan>d the existenpan>ce of hydrologic anpan>d biogeochemical nutrienpan>t legacies. Specifically, we have focused onpan> pan> class="Chemical">nitrogen, such that biogeochemical legacy refers to sorbed organic nitrogen within the root zone, while hydrologic legacy refers to nitrate dissolved in groundwater. The model was able to capture the concentration dynamics in both shallow groundwater beneath sites undergoing landscape conversion as well in-stream concentrations at the catchment outlet. Our findings indicate that the existence of biogeochemical legacy can more than double the time needed to see meaningful concentration reductions at the catchment scale. In addition, we show that while a random approach to landscape conversion will lead to a 1:1 relationship between land-use conversion and maximum concentration reduction benefits at infinite time, a preferential conversion of land parcels with shorter travel times will lead to both faster recovery times and greater maximum achievable concentration reductions.

Our modeling framework provides a first attempt at fully desn class="Chemical">cribinpan>g anpan>d quanpan>tifyinpan>g the oftenpan>-ignpan>ored time lag inpan> catchmenpan>t manpan>agemenpan>t questionpan>s. Inpan> its presenpan>t form, it allows for the quanpan>tificationpan> of tradeoffs betweenpan> costs associated with implemenpan>tationpan> of conpan>servationpan> measures anpan>d the time needed to see the desired conpan>cenpan>trationpan> reductionpan>s, thus makinpan>g it a potenpan>tially powerful tool for lanpan>d manpan>agemenpan>t as agricultural pressures onpan> the enpan>vironpan>menpan>t conpan>tinpan>ue to inpan>tenpan>sify. The anpan>alytical framework also makes the approach conpan>ducive towards makinpan>g unpan>certainpan>ty assessmenpan>ts regardinpan>g conpan>cenpan>trationpan> reductionpan> anpan>d lag time metripan> class="Chemical">cs. In the future, the approach can be further refined by consideration of spatially varying denitrification rate constants, coupled dynamics of denitrification and dissolved organic carbon availability, and by the introduction of hydrologic variability in relation to both rainfall and evapotranspiration dynamics, as they affect the travel-time distribution for the catchment. Such refinements will lead to even further benefits with regard to decision-making support for implementation of conservation measures in intensively managed watersheds.