Riverine export of N and P has been shown to be related to anthropogenic nutrient inputs (Howarth et al., 1996; Alexander et al., 2002; Boyer et al., 2002; Boyer et al., 2006; Mueller and Spahr, 2006; Alexander et al., 2008). For large watersheds globally and in the United States, total N or nitrate (NO3−) export has been described well by population density alone (Peierls et al., 1991); total anthropogenic inputs (Howarth et al., 1996; Boyer et al., 2006); and point sources, fertilizer, and runoff (Caraco and Cole, 1999). A comparison of these models indicated that inclusion of landscape factors affecting water and N transport improved the bias and precision of N export predictions at a large scale (Alexander et al., 2002). Similarly, in studies that included smaller watersheds, total N or NO3− export was found to be strongly related to N inputs (Caraco et al., 2003; Mueller and Spahr, 2006), but the inclusion of additional landscape factors such as runoff, soil permeability, air temperature, precipitation, and tile drainage improved model accuracy (Caraco et al., 2003; Alexander et al., 2008; Spahr et al., 2010). Total P export has also been shown to be related to both anthropogenic inputs and landscape factors affecting water and P transport, such as soil permeability, slope, and precipitation (Mueller and Spahr, 2006; Alexander et al., 2008; Spahr et al., 2010).
One potentially influential landscape factor that has not been extensively studied at the large watershed scale is the implementation of management practices—practices designed to prevent or reduce nutrient loss from agricultural watersheds. At the smaller scale, the effects of management practices have been evaluated by using research plots, individual farm fields, or small watersheds (i.e., generally <50 km2) (Yoo et al., 1988; Sharpley and Smith, 1994; Tan et al., 1998; Bundy et al., 2001; Owens et al., 2008; Sutton et al., 2010; Lemke et al., 2011). Larger-scale assessments, however, have been limited by the lack of comprehensive data on the timing and location of management practices. Recently, the areas in conservation tillage and the Conservation Reserve Program (CRP) have been estimated for the conterminous United States over multiple time periods. These new data sets provide an opportunity to evaluate the influence of management practices on nutrient export at the national scale.
Conservation tillage was used on 26% (113 million acres) of cropland in the United States in 2002 (Baker, 2011). Conservation tillage, which limits tilling while retaining a minimum of 30% of crop residue on the soil surface, potentially increases infiltration and reduces erosion and runoff to rivers relative to conventional tillage, which exposes the soil surface by plowing crop residues into the soil. Approximately 8% (34 million acres) of cropland was enrolled in the CRP in 2002 (http://www.apfo.usda.gov/Internet/FSA_File/historycounty.xls). Farmers, ranchers, and other producers enrolled in the CRP convert highly erodible and other environmentally sensitive farmland to filter strips, grassed waterways, riparian buffers, and long-term vegetative covers such as introduced or native grasses or hardwood trees to reduce the need for fertilizer and manure, control soil erosion, and enhance wildlife habitat. In this study, we evaluate the relations between the riverine export (load) of total N and total P from 133 agricultural watersheds in the United States and both anthropogenic inputs and landscape factors affecting water and nutrient transport—including area in conservation tillage and in the CRP.
Materials and Methods
Study Sites and Data Collection
In 1991, USGS began a study of water quality in 51 large river basins across the United States through the National Water-Quality Assessment program (Gilliom et al., 1995). Data collection in the 51 river basins was implemented on a rotational schedule, with 20 river basins sampled from 1993 to 1995, 16 river basins sampled from 1996 to 1998, and 15 river basins sampled from 1999 to 2001. The monitoring network in each river basin typically included 8 to 12 sites that were routinely sampled, for a total of 495 stream and river sites across the United States. About half of the sites had relatively small watersheds with homogenous land use, whereas the remainder had larger watersheds integrating multiple land uses. Land use in the watershed of each site was determined using the USGS Enhanced National Land Cover Dataset 1992 (NLCDe 92; Nakagaki et al., 2007).
The sites were sampled at fixed frequency on a monthly basis and periodically during high streamflow conditions—the overall sampling frequency at each site varied from bi-weekly to as often as every other day (Mueller and Spahr, 2006). Water samples were collected using depth- and width-integrating techniques and processed and preserved onsite using standard methods described by the U.S. Geological Survey (2009). Samples were analyzed for dissolved nitrite (NO2−) plus NO3−, total Kjeldahl N (TKN), and total P at the USGS National Water-Quality Laboratory in Colorado according to methods described by Patton and Truitt (1992), Fishman (1993), the U.S. Environmental Protection Agency (1993), Patton and Truitt (2000), and Mueller and Spahr (2006). Dissolved NO2− plus NO3− and TKN were summed to calculate total N. Routine sampling included field blanks and replicates to assess measurement bias and variability (Mueller et al., 1997). Continuous streamflow data were collected at co-located or nearby streamgages. These nutrient concentration and streamflow data, along with detail on data review and revision, are available in Mueller and Spahr (2005).
For the present study, 133 of the 495 stream and river sites were used—these 133 had relatively small watersheds with predominantly agricultural land use (Fig. 1). All but 8 of the 133 sites had watersheds with >50% agricultural (crop plus pasture) land use and <5% urban land use (Table 1). The remaining eight sites had <50% agricultural land use in their watersheds, but a disproportionate amount of the water at the site originated from agricultural land area in the watershed as a result of water management practices such as canal diversions and reservoir storage. About 80% of the sites had more cropland than pasture. The sites were located in major agricultural regions of the United States, including the Great Appalachian Valley in the East, the Corn Belt in the Midwest, and irrigated areas in the West. The sites provided good coverage of the major agricultural regions, but some areas—such as in the southeastern, south-central, and north-central United States—were underrepresented (Fig. 1); this may have led to some degree of spatial bias. Sixty-one of the 133 sites were located in the Midwest, which is not unreasonable given that approximately 46% of all harvested cropland in the United States is in the 12 midwestern states (U.S. Department of Agriculture, 2009). Watershed areas at all sites ranged from 17 to 121,182 km2, with a median of 763 km2 (Table 1).
|Value||Drainage area||Land use, in % of drainage area
|Total number of sites||133||133||133||133||133||133|
Concentration and streamflow data from the study sites were used by Mueller and Spahr (2005) to determine the mean annual export (load) of total N and total P from each watershed. These mean annual export values were used in our study as response variables in models relating nutrient export to both anthropogenic inputs and landscape factors affecting water and nutrient transport. Mean annual export was estimated by Mueller and Spahr (2005) by first estimating daily export using multiple regression models that related export to streamflow, time, and season:where ln(L) is the natural logarithm of export (computed as the product of mean daily streamflow and measured concentration on the date of sample collection), b0 and bj are model coefficients, Xj is an explanatory variable, and n is the number of explanatory variables (Runkel et al., 2004). Potential explanatory variables include linear and quadratic natural logarithm of mean daily streamflow, linear and quadratic decimal time, and the sine and cosine of decimal time. All possible combinations of these variables were examined, and the best model was selected on the basis of the Akaike Information Criteria (Akaike, 1981). The sine and cosine terms, which account for seasonality, were always included together if either was selected. Export models were calibrated separately for each constituent at each site. Once the export models had been calibrated, daily export values were estimated using mean daily streamflows from the respective period of nutrient data collection (either 1993–1995, 1996–1998, or 1999–2001). Mean annual values were then calculated as the sum of the daily export values divided by the number of years used in the estimation.
Because nutrient concentrations included censored values, regression coefficients were determined by an adjusted maximum-likelihood estimation (AMLE) method (Cohn et al., 1992). The AMLE method corrects for bias in the standard maximum-likelihood (MLE) regression coefficients and also incorporates a factor that minimizes the bias that can occur when estimated logarithms of constituent export are retransformed to original units. For some constituents at some sites, a model could not be fit using the AMLE method; however, an MLE model was successful. In these cases, retransformation bias was corrected using the method of Bradu and Mundlak (1970). At other sites, a model could not be fit because of a lack of uncensored observations or poor model fit; these sites were excluded from further analysis. The final number of sites was 133 for total N and 129 for total P.
The regression coefficients, selected model diagnostics, and mean annual export estimates for total N and total P at the study sites are available in Mueller and Spahr (2005). Mean annual export estimates, normalized to site drainage areas, also are shown in Fig. 1.
Export values estimated using the regression-based approach in Eq.  have been shown to be biased in certain situations where the regression model inadequately describes nutrient variation in the stream—this can result in bias when estimated logarithms of export are retransformed to original units (Stenback et al., 2011). The potential bias in the total N and total P export estimates used in this study was evaluated by directly comparing observed and predicted export in original units. There was good correspondence between the observed and predicted export values and no evidence of systematic bias across sites (Supplementary Fig. S1 and S2). See the Supplementary Material for more detail on the evaluation of potential bias in the export estimates.
Statistical Evaluation of Variability in Nutrient Export
Landscape and input factors associated with nutrient export from the 133 agricultural watersheds across the United States were evaluated by using multiple linear regression. Although both concentration and export provide insight into nutrient condition in streams, export was used as the dependent variable here because it more closely reflects downstream nutrient movement in large watersheds; concentration more closely reflects local conditions (Hirsch et al., 2010). Eleven potentially influential variables were considered as explanatory variables: (i) total anthropogenic inputs of N or P, in kilograms and log transformed; (ii) mean annual runoff, in cm and log transformed; (iii) mean annual temperature from 1980 to 1997, in °F; (iv) mean annual precipitation from 1980 to 1997, in inches; (v) mean watershed slope, in percentage; (vi) soil erodibility factor (K-factor of the Universal Soil Loss Equation, principally affected by soil texture, structure, permeability, and organic matter content) in the uppermost soil horizon, unitless; (vii) number of reservoirs in the watershed, in number km−2; (viii) irrigated area, in percentage of agricultural area in the watershed; (ix) area with subsurface tile drains, in percentage of agricultural area in the watershed; (x) area in the CRP, in percentage of agricultural area in the watershed; and (xi) area in conservation tillage, in percentage of agricultural area in the watershed. Nutrient export values, in kilograms, were log transformed before these analyses.
Apart from area in the CRP and conservation tillage, these explanatory variables were chosen because they had previously been shown to be important factors affecting nutrient transport in field or modeling studies (e.g., Rolston et al., 1982; Sexstone et al., 1985; de Klein and van Logtestijn, 1996; Howarth et al., 1996; Caraco and Cole, 1999; Fraser et al., 1999; McDowell et al., 2001; Alexander et al., 2002; Boyer et al., 2002; McIsaac and Hu, 2004; Boyer et al., 2006; Mueller and Spahr, 2006; Royer et al., 2006; Vale et al., 2007; Alexander et al., 2008; David et al., 2010; Brown et al., 2011). If one or more of these explanatory variables was simultaneously affecting nutrient export from the study watersheds along with area in the CRP and/or area conservation tillage, their inclusion in a multiple regression model will make any relations between export and area in the CRP and/or conservation tillage more discernible. In addition, 11 explanatory variables (plus the 6 interaction terms described subsequently) resulted in a reasonable ratio between the number of explanatory variables and the number of observations of export. These 11 variables are likely related to, and thus are indirect surrogates for, a number of other potentially influential variables such as land use, population, soil type, and biological activity. Other variables that might also be influential but are not related directly or indirectly to these 11 variables (such as natural inputs of nutrients, total reservoir capacity, or irrigation volumes) have not been quantified at a national scale.
An additional six interaction terms were included as potential explanatory variables to explore the possibility that the relation between nutrient export and the area in management practices varied with total anthropogenic inputs, soil erodibility, or the area with subsurface tile drains. These six terms were (i) total anthropogenic inputs:area in conservation tillage, (ii) total anthropogenic inputs:area in the CRP, (iii) soil erodibility:area in conservation tillage, (iv) soil erodibility:area in the CRP, (v) area with subsurface tile drains:area in conservation tillage, and (vi) area with subsurface tile drains:area in the CRP. Because conservation tillage and the CRP in part are intended to reduce nutrient transport by reducing soil erosion, the greatest effect on nutrient export may be realized where soils are most erodible or surface inputs of nutrients are greatest. Further, because subsurface tile drains redirect nutrient-laden water to rivers rather than allowing it to infiltrate into groundwater where it can take years to travel to rivers, the most rapid effect on nutrient export may be realized where the percentage of subsurface tile drainage is greatest. To reduce multi-colinearity, each of these variables was centered by subtracting the respective mean from each observation before calculating an interaction term as the product of the two centered variables. If the interaction term was not included in the final model but one of the individual variables was included, the uncentered original values were used.
Total anthropogenic inputs of nutrients to agricultural watersheds included fertilizer, biological fixation by legume crops (for N only), atmospheric deposition (for N only), and imported human food and livestock feed. In a hypothetical watershed over the long term, nutrients in crops and imported food and feed are recycled as human and livestock waste; as such, human and livestock waste are not considered to be “new” inputs. In shorter time steps (such as a year), however, livestock waste is not always completely recycled; manure is often used to fertilize crops and pasture in succeeding time steps. Therefore, in this shorter-term cycle, anthropogenic inputs to a watershed during year i are the following:where Ferti is the input from fertilizer in year i, Fixi is the input from biological fixation by crops in year i (N only), Atmi is the input from atmospheric deposition in year i (N only), FoodImpi is the input from imported food in year i, FeedImpi is the input from imported feed in year i, and Wastei-1 is the input from recoverable manure in year i − 1. Because riverine export was determined on annual basis in Mueller and Spahr (2006), corresponding anthropogenic nutrient inputs in Eq.  also were estimated on an annual basis for this study. The estimation of each input term in Eq.  is described in the Supplementary Material in the section “Estimation of anthropogenic nutrient inputs.”
Mean annual runoff was calculated for the 3-yr data collection period by dividing mean streamflow for those years by the watershed area. Mean annual temperature, mean annual precipitation, and mean watershed slope were characterized by the average of the grid-cell data values in the watershed. Mean annual temperature and precipitation were derived from 1-km resolution gridded data from the Daymet conterminous United States database (Thornton and Running, 1999; Daymet, 2011); these data only were available from 1980 through 1997. Mean watershed slope was derived from the USGS National Elevation Dataset (U.S. Geological Survey, 2003) gridded at 100-m resolution. Mean soil erodibility was derived from 100-m resolution grids of State Soil Geographic (STATSGO) parameters compiled by the Natural Resources Conservation Service (1994). The number of reservoirs was derived from location information of approximately 2700 reservoirs and controlled natural lakes that have normal capacities of at least 6,167,445 m3 (5000 acre-feet) or maximum capacities of at least 30,837,225 m3 (25,000 acre-feet) and that were completed as of 1 Jan. 1988 (Ruddy and Hitt, 1990). The agricultural area with irrigation and subsurface tile drains was derived from county-level aggregations of 1997 and 1992, respectively, National Resources Inventory parameters compiled by the Natural Resources Conservation Service (U.S. Department of Agriculture, 1995, 2001); more recent data were not available at a national scale. Agricultural land area used to calculate the percentage of agricultural area in irrigation and subsurface tile drains was derived from the NLCDe 92.
The area in the CRP was derived from county-level data compiled by the USDA Farm Service Agency (Fig. 2a) (http://www.apfo.usda.gov/Internet/FSA_File/historycounty.xls).
Under the CRP, the USDA establishes contracts with agricultural producers to convert highly erodible and other environmentally sensitive cropland and pasture to grass, trees, and other vegetative cover for 10 to 15 yr. Since its inception in 1985, the CRP has been expanded to include practices such as filter strips, riparian buffers, and restoration of farmable non-floodplain wetlands. Annual data on cumulative acres enrolled in the CRP were available, and the mean from the years of nutrient data collection at each site was used (Mueller and Spahr, 2005). These data reflect acres in any CRP practice; detailed information on individual practices is not available. Agricultural land area used to calculate the percentage of agricultural area in the CRP was derived from the NLCDe 92.
The area in conservation tillage was derived from data for 8-digit hydrologic unit watersheds compiled by using information from the Conservation Technology Information Center (CTIC) (Fig. 2b) (Baker, 2011). The CTIC collects tillage data by conducting surveys about tillage systems for all counties in the United States. Total annual acreage for each tillage practice for each crop was reported by local conservationists to the CTIC. For 1989–1998, surveys were based on local knowledge and expertise of the Natural Resources Conservation Service, State Soil and Water Conservation Districts, State extension agents, Farm Service Agency, and other agribusiness partners to estimate tillage practices. For 2000, 2002, and 2004, 17 states used roadside-transect-survey procedures for counties with more than 100,000 acres of cropland (Conservation Technology Information Center, 2004). Conservation tillage included herein was any cropland system that leaves at least one-third of the soil covered with crop residue after planting, including no-till/strip-till, ridge-till, and mulch-till (Baker, 2011). With no-till, producers disturb only the minimal amount of soil needed to ensure a good stand and yield. No-till variations included Midwest strip-till, southeast strip-till, vertical tillage, and fluffing harrows (Baker, 2011). Other conservation tillage practices included ridge-till and mulch-till (Baker, 2011). Conservation tillage data were available for 1992, 1997, and 2002, corresponding to the years of the Census of Agriculture (Baker, 2011). The closest census year corresponding to the nutrient data-collection period at each site was used (Mueller and Spahr, 2005); spatial coverage in 1997 was relatively sparse, so data from 2002 were used at the corresponding sites instead. Because of double cropping (two principal row crops in one growing season on the same land area), agricultural land area from the NLCDe 92 could not be used to calculate the percentage of agricultural area in conservation tillage. Instead, the summed area (as reported by the CTIC) of full season and double cropped corn, grain sorghum, and soybeans; cotton; forage crops; spring seeded and fall seeded small grains; fallow; newly established permanent pasture; and other crops including vegetables was used in this calculation.
As described by Nakagaki and Wolock (2005), annual estimates of nutrient inputs were mapped to agricultural or urban land use within each county or hydrologic unit using the NLCDe 92 (Supplementary Table S7); estimates of irrigated area, area with subsurface tile drains, area in the CRP, and area in conservation tillage were mapped to agricultural land use within each county or hydrologic unit. The area within each county or hydrologic unit was then summed to determine the area within each watershed.
The 17 variables used in regression modeling were evaluated using an all-subsets regression approach, with export of total N and total P modeled separately. The model with the lowest Mallow’s Cp value was selected as the final model. If an interaction term was selected for inclusion in the model based on Mallow’s Cp, both individual variables were retained for hierarchical completeness, even if they were not originally selected for inclusion on their own. The final models were verified for normality and homoscedasticity of the residuals, and variance inflation factors were verified as being below a value of two. Diagnostic plots from the final regression models are available in Supplementary Fig. S3–S7. The regression analyses were performed in S+ version 8.1 (Tibico Software, Inc., 2008).
Results and Discussion
Factors Affecting Nutrient Export
In the final N model, which explained 93% of the variability in N export, N export was related to total anthropogenic N inputs, runoff, slope, number of reservoirs, irrigated area, area with subsurface tile drains, area in the CRP, soil erodibility, area in conservation tillage, and the interaction term soil erodibility:area in conservation tillage (Table 2). The coefficients on number of reservoirs and irrigated area were negative, indicating that those factors were associated with a decrease in N export. The coefficients on total anthropogenic N inputs, runoff, slope, area with subsurface tile drains, and area in the CRP were positive, indicating that those factors were associated with an increase in N export. With the interaction between soil erodibility and the area in conservation tillage, both variables were centered at their respective means. As a result, the coefficient for soil erodibility indicated its association with export at the mean area in conservation tillage (33% of agricultural area; see Supplementary Fig. S8 for full range of data); likewise, the coefficient for the area in conservation tillage indicated its association with export at the mean soil erodibility (0.312; see Supplementary Fig. S8 for full range of data). At 33% agricultural area in conservation tillage, an increase in soil erodibility was associated with a decrease in N export. At a soil erodibility factor of 0.312, an increase in the area in conservation tillage was associated with an increase in N export. Across all values of soil erodibility, the association between area in conservation tillage and N export (holding all other model variables constant) was described using the respective coefficients from Table 2:Therefore, for centered soil erodibility factors greater than 0.0996, the relation between N export and area in conservation tillage was negative and an increase in the area in conservation tillage was associated with a decrease in N export. For centered soil erodibility factors <0.0996, the relation between nitrogen export and the area in conservation tillage was positive. A centered value of 0.0996 for soil erodibility corresponds to an uncentered value of 0.412, approximately the 90th percentile of the uncentered soil erodibility factors across all study watersheds.
|Parameter (units)||Model coefficients
|Value||SE||t value||p value|
|Total anthropogenic inputs of nitrogen (kg, log transformed)||0.970||0.0281||34.5||<0.001|
|Runoff (cm, log transformed)||0.755||0.0539||14.0||<0.001|
|Number of reservoirs (number km−2)||−243||113||−2.16||0.033|
|Irrigated area (% of agricultural area)||−0.00580||0.00200||−2.93||0.004|
|Area with subsurface tile drains (% of agricultural area)||0.0105||0.00280||3.70||0.000|
|Area in the Conservation Reserve Program (CRP) (% of agricultural area)||0.0151||0.00880||1.73||0.086|
|Soil erodibility factor (unitless, centered)||−0.493||0.800||−0.616||0.539|
|Area in conservation tillage (% of agricultural area, centered)||0.00550||0.00290||1.91||0.059|
|Soil erodibility factor : Area in conservation tillage†||−0.0552||0.0337||−1.64||0.103|
|Degrees of freedom||122|
|Total anthropogenic inputs of P (kg, log transformed)||1.07||0.0428||25.0||<0.001|
|Runoff (cm, log transformed)||0.555||0.0949||5.85||<0.001|
|Number of reservoirs (number km−2)||−310||173||−1.79||0.076|
|Area in conservation tillage (% of agricultural area)||0.00600||0.00390||1.55||0.123|
|Soil erodibility factor (unitless, centered)||3.07||1.04||2.95||0.004|
|Area in the CRP (% of agricultural area, centered)||0.0532||0.0156||3.41||0.001|
|Soil erodibility factor : Area in the CRP†||−0.413||0.130||−3.19||0.002|
|Degrees of freedom||119|
There was no interaction between soil erodibility and area in the CRP in the N model. The coefficient on area in the CRP was positive, indicating that an increase in area in the CRP was associated with an increase in N export in all watersheds.
In the final P model, which explained 86% of the variability in P export, P export from agricultural watersheds was related to total anthropogenic P inputs, runoff, precipitation, slope, number of reservoirs, area in conservation tillage, soil erodibility, area in the CRP, and the interaction term soil erodibility:area in the CRP (Table 2). The coefficient on number of reservoirs was negative, indicating that a greater number of reservoirs in the watershed was associated with a decrease in P export. The coefficients on total anthropogenic P inputs, runoff, precipitation, slope, and area in conservation tillage were positive, indicating that those factors were associated with an increase in P export. With the interaction between soil erodibility and area in the CRP, the variables were centered at their respective means. As a result, the coefficient for soil erodibility indicated its association with export at the mean area in the CRP (4.8% of agricultural area; see Supplementary Fig. S8 for full range of data); likewise, the coefficient for area in the CRP indicated its association with export at the mean soil erodibility. At 4.8% agricultural area in the CRP, an increase in soil erodibility was associated with an increase in P export. At a soil erodibility factor of 0.312, an increase in the area in the CRP was associated with an increase in P export. Using logic established in Eq.  and  and the P coefficients in Table 2, the relation between P export and area in the CRP was negative only for centered soil erodibility factors greater than 0.129. A centered value of 0.129 for soil erodibility corresponds to uncentered value of 0.442, or approximately the 95th percentile of the uncentered soil erodibility factors across all study watersheds.
There was no interaction between soil erodibility and area in conservation tillage in the P model. The coefficient on area in conservation tillage was positive, indicating that an increase in the area in conservation tillage was associated with an increase in P export in all watersheds.
Both N and P export were most strongly associated with total anthropogenic inputs and runoff. The positive relation with anthropogenic inputs has been well established previously (Howarth et al., 1996; Alexander et al., 2002; Boyer et al., 2002; Boyer et al., 2006; Mueller and Spahr, 2006; Alexander et al., 2008). Nitrogen and P export increased with increasing runoff, likely because of greater water availability and movement in areas with higher runoff, which increases nutrient transport from the land surface to rivers. Nutrient export has previously been found to be related to runoff globally (Caraco et al., 2003), in the Mississippi River basin (McIsaac et al., 2002), in California (Sobota et al., 2009), and along the east coast of the United States (Schaefer and Alber, 2007). Nitrogen and P export also increased with increasing slope and P export increased with increasing precipitation; higher slopes and greater precipitation can contribute to more erosion and runoff (Fraser et al., 1999; McDowell et al., 2001).
The presence of a larger number of reservoirs in the watershed was associated with a decrease in the export of N and P. Nitrogen is retained in reservoirs largely through a combination of denitrification and sedimentation (Nõges et al., 1998; Saunders and Kalff, 2001); P is retained in reservoirs largely through sedimentation (Canfield and Bachmann, 1981; Kõiv et al., 2011). Spatial differences in N and P retention throughout the Missouri River basin were previously found to be related in part to differences in the total number of reservoirs and lakes, with a larger number of waterbodies contributing to greater retention of nutrients (Brown et al., 2011).
Increases in tile drainage were associated with increased export of N, but not P. Tile drains have previously been found to increase the transport of N—particularly NO3−—to rivers by reducing anoxia and associated denitrification in soils and reducing the interaction between water and oxic soils or riparian areas (McIsaac and Hu, 2004; Royer et al., 2006; David et al., 2010). Tile drains have also previously been found to increase the transport of P, particularly in the dissolved phase, to rivers by reducing the interaction between water and soils and riparian areas where retention and transformation can occur (Sims et al., 1998; Gentry et al., 2007). However, other studies have found that tile drainage decreases the transport of P by reducing soil moisture and increasing soil redox, thereby enhancing phosphate sorption and retention (Kröger et al., 2008). In this study, tile drains were not related to P export, possibly because of the aggregate effects of factors that promote and hinder its transport to rivers.
Increases in irrigated agriculture were associated with decreases in the export of N, but not P. Increases in soil moisture with irrigation can lead to the formation of anoxic and hypoxic zones that allow for denitrification (Colbourn and Dowdell, 1984; Bhattarai et al., 2005). Numerous studies have documented an increase in soil denitrification rates with irrigation (Rolston et al., 1982; Sexstone et al., 1985; de Klein and van Logtestijn, 1996; Aulakh and Bijay-Singh, 1997; Mahmood et al., 1998; Hooda et al., 2003; Mahmood et al., 2005; Vale et al., 2007). At a large scale, irrigation was found to be associated with a decrease in N loads in the Missouri River basin (Brown et al., 2011). Other studies, however, have found that irrigation increases leaching of N to groundwater (Mossbarger and Yost, 1989; Power and Schepers, 1989; Stites and Kraft, 2001), where it can move to streams. It is likely that the effects of irrigation in these watersheds are complex, but the results suggest that the cumulative effect is one of loss, resulting in a net decrease in the export of N with greater irrigated area.
The final models indicated that the relations between export and area in conservation tillage (N) and area in the CRP (P) varied with soil erodibility. Figure 3 shows predicted N and P export with increasing area in these respective management practices at different levels of soil erodibility. Partial predicted export—reflecting the effect of the respective management practices, both through their main effect and their interaction effect with soil erodibility, independent of changes in the other model variables—first was calculated aswhere ln(Lp,i) is the partial predicted export at site i, bmp is the coefficient from the final model (Table 2) for area in the respective management practice, bmp:se is the coefficient from the final model for the interaction term of soil erodibility and area in the respective management practice, Xmp,i is the centered area in the respective management practice at site i, and Xse is a specific percentile of centered soil erodibility across all sites. Partial predicted export values were determined separately for six levels of Xse (the 10th, 25th, 50th, 75th, 90th, and 95th percentiles of centered soil erodibility in our data). These partial predicted export values then were mean adjusted to put them on a scale comparable with the observed export values by adding a constant to each prediction. This constant was calculated as the value predicted by the final model—with the main and interaction terms that included the respective area in management practices excluded—at the mean of the other explanatory variables. At most levels of soil erodibility, increases in the area in conservation tillage were associated with increased export of N and increases in the area in the CRP were associated with increased export of P. The increase in export became less pronounced as soils became more erodible. The slope of the relation between predicted export and the area in management practices progressed from being clearly positive when soil erodibility was low or moderate, to being close to zero when soil erodibility was higher, to possibly being slightly negative at the highest levels of soil erodibility in our data. As noted earlier, this occurred around the 90th percentile of soil erodibility for N export and area in conservation tillage (Fig. 3a) and the 95th percentile of soil erodibility for P export and area in the CRP (Fig. 3b).
Figure 3 also shows partial residuals for the area in conservation tillage (N) and the area in the CRP (P), calculated as:where rp,i is the partial residual at site i, rf,i is the ordinary residual from the full model at site i, and Xse,i is the centered soil erodibility at site i. The partial residuals were mean adjusted as described above to put them on a scale similar to the observed export values. The partial residuals then were subset into six groups based on ranges of centered soil erodibility that corresponded to the percentiles used to predict export (0–15th, 15th–35th, 35th–65th, 65th–85th, 85th–92nd, 92nd–100th percentiles of soil erodibility in our data). The partial residuals show the variation in N and P export remaining after the removal of variation attributable to the other model variables (e.g., the remaining variation in export is likely related to changes in the respective management practices and any other unidentified variables affecting export). As a result, these partial residuals provide a clearer indication of the relation between export and area in the respective management practices at the study sites than could be seen using the observed values of export instead. The pattern of partial residual points, therefore, provides a measure of validation for the lines of partial predicted export. In general, the patterns of partial residuals corresponded well to the patterns of partial predicted export for varying ranges of soil erodibility. At the highest levels of soil erodibility, however, the patterns of partial residuals were not as clearly negative as the patterns of partial predicted export. As such, the relations between export and area in the respective management practices were indeterminate at the highest levels of soil erodibility.
There also were positive associations between area in the CRP and N export and area in conservation tillage and P export across all watersheds; there were no interactions with soil erodibility with these variables (Table 2). These results indicate that an increase in area in the CRP was associated with an increase in N export in all watersheds, regardless of the soil erodibility. Similarly, an increase in area in conservation tillage was associated with an increase in P export in all watersheds, regardless of the soil erodibility.
One of the primary goals of both conservation tillage and the CRP is to reduce soil erosion and runoff, thereby reducing nutrient transport to streams. Small-scale field and modeling studies have shown that both types of management practices can be effective at reducing soil erosion in areas with erodible soils (e.g., McDowell and McGregor, 1984; Uri et al., 1998; Nyakatawa et al., 2001; Zheng et al., 2004). At larger scales, however, trend analyses of nutrients in streams have found that management practices have not led to widespread measureable changes in nutrient concentration or export over time, possibly because improvements from management practices at small scales are counteracted or overwhelmed by changes in other factors affecting nutrients at larger scales (Sharpley et al., 2009; Sutton et al., 2010; Lemke et al., 2011; Sprague et al., 2011). Moreover, conservation tillage and the CRP can be implemented in areas with less erodible soils, where other benefits—such as improved wildlife habitat, recreational opportunities, carbon retention, crop yields, and farm income—may be realized (Pimentel et al., 1995; Uri et al., 1998; Feather et al., 1999; Nyakatawa et al., 2001). In these areas, a relatively smaller proportion of nutrients may be reaching the stream in particulate form via surface runoff as compared with areas with more erodible soils, so runoff reductions may have a proportionally smaller effect on nutrient transport. On the basis of these previous studies, a lack of a strong negative relation between area in the CRP and/or conservation tillage and nutrient export would not be unexpected. Our results, however, indicate that relations between nutrient export and area in management practices generally were positive—indicating greater export in areas with more management practices—and progressed to being slightly negative only in areas with more erodible soils. There are several possible explanations for this.
First, although the transport of particulate N and/or P in surface runoff can be lower from areas in conservation tillage compared with areas in conventional tillage (McDowell and McGregor, 1984; Yoo et al., 1988; Sharpley and Smith, 1994; Bundy et al., 2001; Yates et al., 2006), the transport of soluble nutrients like dissolved P and NO3− is often greater from areas in conservation tillage (Smith et al., 1991; Tan et al., 1998; Tuppad et al., 2010). Accumulation of P in the top layers of soil under fields with conservation tillage, where fertilizer, manure, and crop residue are not incorporated into the soil, can contribute to the enhanced transport of dissolved P (Sharpley and Smith, 1994; Bundy et al., 2001; Sharpley, 2003). A greater number of macropores in relatively undisturbed conservation till areas can enhance leaching of NO3− to groundwater (Sharpley and Smith, 1994; Tan et al., 1998)—infiltrating water can move more rapidly to groundwater through macropores than the surrounding soil matrix (Beven and Germann, 1982; Newman et al., 1998), reducing the amount of denitrification in the soil zone (Li and Ghodrati, 1994; Tindall et al., 1995). Also, a large source of NO3− leaching is the mineralization of soil organic matter (Jarvis, 2007); soil organic matter content can increase under conservation tillage compared with conventional tillage (Tan et al., 1998). In addition, less evaporation in conservation till areas limits movement of NO3− upward to the surface where it could be taken up by vegetation (Sharpley and Smith, 1994).
Second, there may be a lag between implementation of these management practices and any resulting changes in water quality in streams. This lag could occur for a number of reasons. The relatively slow movement of NO3− through groundwater to streams can cause changes in streams to lag changes on the land surface by years to decades (Tomer and Burkart, 2003; Bratton et al., 2004; Schilling and Wolter, 2007; Owens et al., 2008). The input of groundwater as base flow to these streams, estimated as the percentage of base flow to total streamflow (the base-flow index, see Supplementary Material for more detail), ranged from 10 to 76% with a median of 40% (Supplementary Fig. S9), so any groundwater lags likely affected stream quality to some degree. Restored natural vegetation in riparian buffers may not develop full N and P retention capacity for many years, leading to a delay between restoration and water-quality improvements (Meals et al., 2010; Sutton et al., 2010). Further, when P has accumulated in soils, transport of dissolved P can continue for years to decades after inputs are reduced (Zhang et al., 2004; Sharpley et al., 2009; Meals et al., 2010; Tomer et al., 2010). Depending on the length of any lags, nutrient conditions in agricultural streams may still be reflecting the influence of agricultural land-use practices that were in place before the implementation of the management practices. The positive relations between nutrient export and management practices may in part be a reflection of more intensive past agricultural activity in areas that currently have a greater concentration of management practices.
Third, there may be limitations with the management practice and stream monitoring data sets. In a study by Jackson-Smith et al. (2010), field interviews were conducted with conservation program participants in Little Bear River watershed in Utah to test the validity and reliability of USDA contract data. Several potential limitations in the data were identified, including incorrect characterization of practice type, location, and spatial extent. A similar, comprehensive national evaluation of management practice data sets has not been conducted, but if the limitations identified in Little Bear River watershed are typical of national data sets, they may be affecting the evaluation of relations between management practices and nutrient export in this study. This evaluation may also be affected by some degree of spatial bias in the coverage of the study sites—some of the agricultural areas in the central United States with a relatively high percentage of area in the CRP were underrepresented in this study (Fig. 1 and 2).
Finally, the spatial distribution of area in the CRP and conservation tillage may closely correspond to the spatial distribution of other potentially influential variables, leading to spurious correlation. If that were the case, the increase in the export of N and P associated with increases in area in the CRP and conservation tillage may actually be the result of spatial variation in some other factor. The same issue is a possibility with any of the other variables included in the final regression models, as spurious correlation is a possibility in any regression analysis. However, the relations between nutrient export and the other variables have been more thoroughly established in previous studies. We evaluated the possibility of spurious correlation by checking for colinearity between area in the CRP and conservation tillage and two other potentially related variables, percentage of the watershed in agricultural land use and rate of glyphosate use (see Supplementary Material for detail on the derivation of these variables). These variables were considered because management practices generally are targeted in more intensively agricultural areas and glyphosate degrades to inorganic phosphate and competes for the same soil sorption sites as P—both factors that may contribute to enhanced transport of P in the presence of glyphosate (de Jonge et al., 2001; Laitinen et al., 2009). The forced inclusion of these individual variables in a model that also included the variables from the final models in Table 2 did not substantially alter the results. Each of the new variables was not strongly related to export, the variance inflation factors were all below two, and changes to the magnitude and p values of the other coefficients were small and within the standard errors of the original estimates. These results do not rule out the possibility that there are other unidentified and/or unquantified variables that may be related to area in the CRP and conservation tillage. For example, management practices tend to be targeted in areas with the potential for high soil and nutrient loss. We have attempted to control for that potential by including factors that can affect soil loss (such as runoff, soil erodibility, and slope) in the multiple regression models, but these variables may not have captured all important factors affecting soil loss. If the area in management practices is related to other unidentified factors affecting soil and nutrient loss, relations between management practices and nutrient export would likely also reflect the effects of the unidentified variables—for example, there would be some balance between factors promoting nutrient loss and factors hindering nutrient loss. Ultimately, a goal of these management practices is to achieve a net reduction in nutrient export, even in areas with the potential for high soil loss. In this study, the relations between management practices and nutrient export were largely positive, indicating that reductions in nutrient export from management practices are not clearly discernible on a large watershed scale.
Notably, there were no interactions between the management practices and the area in subsurface tile drains or total anthropogenic inputs. This indicates that the relations between nutrient export and the area in both conservation tillage and the CRP did not vary based on the amount of tile drainage or total anthropogenic inputs in the watershed.
The national data sets used in this study are lacking in several key aspects that could support a more comprehensive assessment of the effects of management practices on riverine nutrient export. First, other management practices, such as nutrient management and irrigation optimization, have been widely implemented throughout the United States during the study period. National-scale data on the timing and location of their implementation, however, are currently lacking. As a result, the individual or cumulative effect of all agricultural management practices could not be assessed in this study. Second, additional detail on practice type, timing of implementation, duration of practice, and degree of maintenance of the areas in both conservation tillage and the CRP could allow for a more precise accounting and explication of their effects on riverine export. Third, additional insight could be gained by national-scale estimation of groundwater residence times and the associated lag in changes in surface-water quality. Currently, estimates of groundwater residence time are available for only a relatively small number of locations in the United States, limiting the evaluation of the groundwater lag effects in this study. If new national data sets are compiled in the future, subsequent evaluation of riverine nutrient export could provide further insights. Additional examination of lag effects could be conducted by relating riverine nutrient export to lagged conservation tillage or CRP data.
Results from empirical regression modeling indicated that there was a positive relation between nutrient export and area in the CRP (N) and conservation tillage (P) in the study watersheds. Model results also indicated that the relation between export and conservation tillage (N) and the CRP (P) progressed from being clearly positive when soil erodibility was low or moderate, to being close to zero when soil erodibility was higher, to possibly being slightly negative only at the 90 to 95th percentile of soil erodibility values. There are several possible reasons for the unexpected positive relations—greater transport of soluble nutrients from areas in management practices; lagged response of stream quality to implementation of management practices because of NO3− transport in groundwater, length of time for vegetative cover to mature, and/or prior accumulation of P in soils; or some other factor. If changes in nutrient export do lag implementation of management practices, nutrient export may be reduced in the years following those examined in this study. Long-term stream monitoring at the watershed scale can provide future accounting of any changes—lagged or otherwise—resulting from the implementation of management practices.