Where Have All the Nutrients Gone? Long-Term Decoupling of Inputs and Outputs in the Willamette River Watershed, Oregon, United States.

Better documentation and understanding of long-term temporal dynamics of nitrogen (N) and phosphorus (P) in watersheds is necessary to support effective water quality management, in part because studies have identified time lags between terrestrial nutrient balances and water quality. We present annual time series data from 1969 to 2012 for terrestrial N and P sources and monthly data from 1972 to 2013 for river N and P for the Willamette River Basin, Oregon, United States. Inputs to the watershed increased by factors of 3 for N and 1.2 for P. Synthetic fertilizer inputs increased in total and relative importance over time, while sewage inputs decreased. For N, increased fertilizer application was not matched by a proportionate increase in crop harvest; N use efficiency decreased from 69% to 38%. P use efficiency increased from 52% to 67%. As nutrient inputs to terrestrial systems increased, river concentrations and loads of total N, total P, and dissolved inorganic P decreased, and annual nutrient loads were strongly related to discharge. The N:P ratio of both sewage and fertilizer doubled over time but there was no similar trend in riverine export; river N:P concentrations declined dramatically during storms. River nutrient export over time was related to hydrology and waste discharge, with relatively little influence of watershed balances, suggesting that accumulation within soils or groundwater over time is mediating watershed export. Simply managing yearly nutrient balances is unlikely to improve water quality; rather, many factors must be considered, including soil and groundwater storage capacity, and gaseous loss pathways.

Introduction

Anthropogenic alterations of nitrogen (N) and phosphorus (P) cycles are pronounced in the continental United States (U.S.), and have been for quite some time (Jacobs et al., 2017; Sobota et al., 2013). The accelerated use of N and P fertilizers and intensive animal production since the 1970s have increased agricultural productivity, but associated unintentional nutrient releases have had substantial negative environmental, health, and economic consequences (Davidson et al., 2011; MacDonald et al., 2012). Increased energy production, waste generation, and transportation over time also have contributed to increased N emissions (Davidson et al., 2011). Losses of N and P to waterways, and N to the atmosphere, are responsible for costly water and air quality impairment (Dodds et al., 2008; Sobota et al., 2015). The quantification of terrestrial nutrient use as well as monitoring aquatic N and P loading are both essential for understanding controls on nutrient export and, consequently, informing effective nutrient management.

Annual mass balances can provide useful insight into element dynamics at landscape scales and have been applied broadly in watershed science (Bailey et al., 2003; Chadwick et al., 1990; Gentry et al., 2009). Often these efforts have assumed steady state, a simplification that has allowed for the identification and estimation of poorly quantified nutrient fluxes. Although useful, the steady state assumption ignores some important dynamics in watersheds, with potentially important implications for watershed biogeochemistry and water quality management. In reality, watersheds experiencing shifts in the frequency or severity of human activities or natural disturbances are not in steady state and may be either accumulating or losing nutrients on an annual basis (Burt et al., 2011). Without an understanding of these dynamics it is difficult to understand and address water quality impairment associated with anthropogenic alteration of N and P cycling such as eutrophication leading to harmful and nuisance algal blooms, human health concerns, hypoxia associated with fish kills, and shifts in primary production and ecosystem function (Ho & Michalak, 2017; Rabalais et al., 2002; Van Meter et al., 2018).

Taking a longer-term perspective is an important step in moving beyond the annual nutrient budget and accounting for nonstationarity (Bouwman et al., 2009; Chen et al., 2018; Meals et al., 2010). For example, Hale et al. (2015) identified a temporal disconnect between N and P inputs to Northeastern U.S. watersheds and total N (TN) and total P (TP) yields over time. They posited that increased wastewater treatment plants and reservoirs between 1930 and 2000 increased watershed P retention. In contrast, the majority of increases in N inputs occurred as atmospheric deposition in forested sites which could be more easily retained than N from fertilizers. Long-term (multiple-decade) records of both nutrient management and water quality are rare, but recent efforts have begun to untangle how previous management decisions contribute to current patterns in riverine nutrient concentrations. These include mechanistic models (e.g., McCrackin et al., 2018; Van Meter et al., 2016, 2017) as well as statistical approaches (e.g., Daloglu et al., 2013; Goyette et al., 2018; Hale et al., 2013, 2015; Han et al., 2012; Kusmer et al., 2018) and combined frameworks (e.g., Motew et al., 2017). Most of the above work focuses on how legacy nutrients (accumulation of nutrients that are stored and then released to water, soil, or the atmosphere later) affect time lags between terrestrial management and aquatic nutrient-related outcomes. Because each watershed has a specific set of biophysical and management characteristics that affect inputs, transport, and transformations of nutrients as they move from terrestrial to aquatic ecosystems, more research on long-term watershed balances is needed at smaller spatial scales and in different places in order to better explain water quality trends (including legacy nutrients; Chen et al., 2018). Expanding such analyses to different types of watersheds (more diverse climates, agricultural production systems, and other anthropogenic management variability) is also essential so that different regions can learn from each other and so that a general understanding of watershed nutrient dynamics can be established and tested.

In addition to a need for longer-term perspectives on watershed balances, there is a growing role for watershed research that adopts a comparative, multielement approach. Few watershed studies examine N and P dynamics and balances over time. This is important because the stoichiometry of agricultural inputs, for instance, affects total crop production as well as the amount of nutrients remaining in soils. The latter is true because a single nutrient generally becomes limiting, especially with manure application (Sileshi et al., 2017). At the global scale, the human mobilization of N may be resulting in more widespread P limitation in ecosystems (Peñuelas et al., 2012), which could limit crop yields in many developing countries without access to P fertilizers (Peñuelas et al., 2020). Similarly, phytoplankton biomass and species composition are sensitive to the N:P ratio in both freshwater and marine environments (Elser et al., 2007; Frank et al., 2020). The development of harmful algal blooms is partly dependent on the stoichiometry of available nutrients, but also a product of complex local conditions (Anderson et al., 2002; Conley et al., 2009). In order to effectively manage N and P together, it is necessary to quantify them together in a systematic way (Kanter & Brownlie, 2019).

Here we compare four-decade records of watershed N and P inputs and river outputs from the Willamette River Basin (WRB), in Oregon U.S., one of the major tributaries to the Columbia River (Figure 1). Studies in the U.S. that examine both N and P over a period greater than 20 years have mainly focused on the Mississippi River Basin, the Great Lakes, and Chesapeake Bay, all of which have had severe nutrient-related water quality impairment issues. Less effort has centered on nutrient dynamics in systems like those found in the Pacific Northwest region. The Willamette Valley, with its diverse (and specialty) set of crop and livestock activities and Mediterranean climate, deserves further attention as it has distinct challenges with nutrient management (Lin et al., 2019). Western North America produces large amounts of fruits, vegetables, and other specialty crops (e.g., nursery crops, trees, grasses) while also supporting large urban populations; protecting and improving water quality, especially under changing hydrology related to climate change (Raymondi et al., 2013; VanRheenen et al., 2003), throughout the region requires more long-term explorations of nutrient dynamics. Our work builds on previous studies of the WRB showing little change in river TN and TP fluxes between 1973 and 2003 (Fuhrer et al., 1996; Wise, 2007) despite changes in land use, major increases in population, and high groundwater nitrate concentrations in parts of the WRB valley floor (Hinkle et al., 2001; Hoppe et al., 2014; Kite-Powell & Harding, 2006). Collectively, these observations raise questions about how watershed N and P inputs have changed over time and how these trends relate to potential nutrient retention versus riverine export. In this paper we address the following questions:

Figure 1.

Willamette River Basin watershed delineation with a 2012 agricultural P balance from Metson et al. (2017). Dark black outline is the drainage area to the water quality monitoring site for the Willamette River used in this study (USGS site 14211720). Light gray outlines show county boundaries, dark grey outlines are state boundaries, with the white space on the left being the Pacific Ocean. Yellow (low) to red (high) coloring depicts the agricultural P balance (synthetic fertilizer + manure − crop harvest) on agricultural lands (Falcone, 2015). Inset shows the location of the watershed (gray) in the conterminous United States. Supporting information (Table S1) gives a more detailed breakdown of land use. The SI methods text gives other descriptive information, and Figure S1 shows a comparable map for N for the watershed.

How have the magnitudes and relative contributions of terrestrial nutrient sources changed over 40+ years (1969–2012) in the WR?

How have riverine N and P concentrations and yields changed over this period (1972–2013) in the WRB?

Are riverine exports responsive to terrestrial inputs over time for both N and P?

How do N and P behave differently in terms of sources, use efficiency, and potential retention

Materials and Methods

Terrestrial Nutrient Data Sources and Processing

We applied the approach described in Metson et al. (2017) to quantify and distribute N and P sources across the landscape (where feasible) for 10 “target” years between 1969 and 2012. We quantified three sources of inputs for P (inorganic fertilizer P, manure P, and sewage P) and five for N (inorganic fertilizer N, manure N, sewage N, agricultural biological N fixation, and wet N deposition), as well as annual crop harvest removal of P and N. In brief, we used national, state, county, and local data sources to estimate the magnitudes of the sources and sinks above, and, when appropriate, used land use data to proportionally assign these values and obtain annual values for the WRB watershed (Table 1). We used the U.S. conterminous “wall-to-wall” anthropogenic land use trends (NWALT) land use maps (Falcone, 2015) available for 1974, 1982, 1992, 2002, and 2012. We then chose 10 “target” years, one with land-use data and a half-way point between available map years. For every two target years we used the earliest available land use year (i.e., no linear interpolation of land use between years) and applied estimates of fertilizer application, manure production, and crop nutrient harvest from the NuGIS tool (IPNI, 2012) for those target years with available information (1987–2012). Census and other governmental data sets were used for earlier target years (Table 1). Methods and data sources were optimized to allow for quantification of temporal changes in nutrient dynamics. Because different data sources use different methods which can, at times, cause artificial increases and decreases when combined, we aimed to use as few data sources as possible. When use of multiple sources was unavoidable, we used overlapping years to assess methodological effects and correct for them for each time series.

Table 1

Data Sources and Notes When Methodology Differed From Metson Et Al. (2017) for Each of the 10 Target Years Used to Quantify Terrestrial Sources and Sinks of N and P

Nutrient input or output	Target year(s)	Equation	Data source(s)	Major assumptions
Fertilizer	1969, 1974, 1978, 1982	County N or P fertilizer sales/area in cultivated crop and hay-pasture land uses in the county	Alexander and Smith (1990) and Falcone (2015)	-Agri- Agricultural land use for all nutrient inputs and outputs are land use categories 43 (crop production) and 44 (pasture and hay production) from NWALT, which agrees with NLDC categorizations
				-Agricultural land use for 1969 is NWALT 1974, and agricultural land use for 1978 is NWALT 1982
				-Follows same assumptions as Metson et al. (2017) except for using USGS county sales data without any redistribution to other counties, thus differs slightly to the IPNI data in later years from a national perspective
	1987, 1992, 1997, 2002, 2007	Same as above	Falcone (2015) and International Plant Nutrition Institute (IPNI) (2012)	-Followed same assumptions as Metson et al. (2017) for P
				-Agricultural land use for 1987 is NWALT 1992, agricultural land use for 1997 is NWALT 2002, agricultural land use for 2007 is NWALT 2012
	2012	Same as above	Metson et al. (2017) for P and same as above for N	-Metson et al. (2017) assumes a mean fertilizer application by county over agricultural lands as opposed to by crop, season, or field-specific application rates, and assumes that fertilizer sold in one county is applied in the same county (some exceptions where IPNI used smoothing functions)
Manure	1969, 1974, 1978, 1982	(Sum of all animal types (# of animals * annual excretion rate of N or P)/Area in cultivated crop and hay-pasture land uses in the county	Goolsby et al. (1999)	-Uses daily excretion rates per animal per day and assumes 365 days in a year
				-Assumes nutrient excretion per animal per year was constant over time
				-Does not distinguish between different animal ages in the annual inventory total
	1987, 1992, 1997, 2002, 2007		Falcone (2015) and Ruddy et al. (2006)
	2012		Metson et al. (2017) for P and same as above for N	-Metson et al. (2017) only considers recoverable P from confined operations (which is in line with the fact that we do not consider grazing land crop removal and as such should not consideranimal inputs); mean manure application by county over agricultural lands as opposed to crop, season, or field specific application rates; assumes that manure produced in one county is applied in the same county
Crop harvest	1969, 1974, 1978, 1982	((Sum of weight harvested for major crop types * N or P content)/Area in cultivated crop and hay-pasture land uses in the county) * correction factor)	Falcone (2015); IPNI (2012); and USDA (1990, 2007)	-Uses the same major crops (21) and their associated N and P content used by IPNI (2012); see 1987 onward for crop breakdown
				-Assumes nutrient contents were constant through time
				-Correction factors were required to reconcile differences between the IPNI data (more complete crop coverage) and the agricultural census data (incomplete coverage). These factors were developed by comparing the data in the year when both estimates were available, 1987. For N and P, the ratio of major crops to all crops was 0.288 and 0.316, respectively, so the nutrient contents of all crops could be estimated from major crops by multiplying by the inverse, i.e., 3.47 and 3.16, respectively.
				-Does not account for crop harvest when data were not available, or undisclosed (D) in the census
				-Census of Agriculture information was manually extracted from PDFs from volume 1, Part 37: Oregon, Chapter 2: County Data for census years 1969, 1974, 1978, 1982, 1987 available in the archive
				-1969 and 1975 crop harvest (and animal inventoiy) were also verified looking at digital coded (as opposed to pdf) census data (USDA, 1990)
	1987, 1992, 1997, 2002, 2007		Falcone (2015) and IPNI (2012)	-IPNI explicitly accounts for the following 21 major crop categories:alfalfa, apples, barley (all and spring), beans, canola, corn, corn silage, cotton, other hay, peanuts, potatoes, rice, sorghum for grain, soybeans, sugar beets, sugar cane, sunflower, sweet corn, tobacco, wheat, and oranges
				-IPNI also includes an “other” category to account for states (like OR) that have many specialty crops but where the Agricultural Census does not systematically report harvest. The removal for this “other” crop category is derived from the difference between total harvested acres and the harvested acres for the 21 crops that are tracked
				-More information can be found on the NUGIS tool website from IPNI
	2012		Metson et al. (2017) for P and same as above for N	-Metson et al. (2017) uses IPNI and Falcone (2015) data as per above, and assumes mean crop uptake by county over agricultural lands as opposed to crop, season, or field specific yields and P content
Crop N fixation	All		Same as crop harvest	-Uses alfalfa N harvest in the county as the amount fixed as per IPNI (2012)
				-As with crop harvest, we applied a correction factor to account for “other crops” that fix nitrogen to match IPNI estimates. Alfalfa accounted for 56.8% of 1987 total IPNI calculated fixation and thus we applied the appropriate correction factor
				-Assumes the correction factor does not change over time
Wetatmospheric N deposition	1969, 1974, 1978, 1982, 1987		Lamarque et al. (2010)	-Accessed through the US EPA's Global Change Impacts and Adaptation-Critical Loads Mapper (https://clmapper.epa.gov/)
				-Data set is referred to as the IPCC AR5
				-Only estimated wet deposition
				-Uses the IPCC wet N deposition and then we change to the NADP; when the data sets overlap in reporting years they report similar values. Wet deposition seems to account for less than a quarter of total deposition in the IPCC estimates. Total deposition (TDEP) and Community Multiscale Air Quality Modeling System (CAMP) records for wet deposition seem to be of similar magnitude to the total IPCC deposition estimates, indicating a potential difference in methodologies or data record keeping; we thus take a conservative and temporally smooth (no abrupt jumps related to changes in data source) and use the IPCC and NADP data.
				-Uses zonal statistics in ArcGIS to extract a mean value for the watershed of interest for target years
	1992, 1997, 2002, 2007, 2012		NADP Program Office (2017)	-Accessed through the US EPA's Global Change Impacts and Adaptation-Critical Loads Mapper (https://clmapper.epa.gov/)
				-Data set is referred to as the NADP
Agricultural balance	All	Fertilizer + Manure + Crop N fixation + Wet atmospheric N deposition — Crop removal		-As per Equation 1 in the main text
Nutrient Use Efficiency	All	Crop removal/(Fertilizer + Manure + Crop N fixation + Wet atmospheric N deposition) * 100		-As per Equation 2 in the main text
				-There are many ways to calculate nutrient use efficiencies: here we use inputs and outputs of N and P used in the agricultural balance (see rows above for more detail on each input and output and Equation 1 in the main text) which relies most on assumptions and calculation methods put forth in NuGIS (IPNI, 2012)
Human sewage	All	Population x (N or P excreted + P in detergents) x proportion connected to treatment X (1 — (connection to primary, secondary, tertiary, and no discharge facilities x P or N removal from facilities))	Chapra (1980); Garnier et al. (2015); Litke (1999); Minnesota Population Center (2016); Morse et al. (1993); US EPA (2012)	-Assumes human N and P excretion rate is the same for all humans and does not change over time with diet
				-Assumes state (OR 1992) and national level P detergent regulations cause decreases in P detergent emissions
				-Assumes that nutrients that do not go through the sewer system are lost to the environment
				-Uses the same calculations and assumptions as Metson et al. (2017) which build on the methods put forth in Hale et al. (2013) for the North Eastern U.S. for both N and P. We used national estimates of treatment type and nutrient removal for population covered by sewer connections, U. S. decadal census bureau data (1960 to 1990)* on the percentage of the OR population sewered. Data were linearly interpolated between years with available data in order to get target year values and projected forward to 2012 assuming the same linear rate of sewage connections between 1970 and 1980, which was faster than between 1980 and 1990. For population, we used historical decadal census county subdivision records and used zonal statistics to extract the mean population for the watershed. We linearly interpolated between years with available data to estimate sewage N and P for target years.
Major Point Sources	1992, 1997, 2002, 2007,2012	Sum of major facilities (N or P concentration x discharge)	US EPA (2015) DMR Tool	-Only covers major point sources documented within the watershed, which can vary by year ‘Historical US census bureau data are available at: https://www.census.gov/hhes/www/housing/ census/historic/sewage.html

Note. Units for inputs and outputs are all calculated in kg of P or N ha−1 year−1.

We estimated N from biological nitrogen fixation (BNF) in major legume crops (alfalfa, soy, and peanuts; IPNI, 2012). Wet deposition estimates for N were taken from Lamarque et al. (2010) and the NADP Program Office (2017). Due to a lack of consistency in the way dry deposition was estimated over time, we omitted dry N deposition from this analysis (Table 1). Although deposition can be an important input to non-agricultural areas, when agriculture dominates even total deposition is often a small contributor. In the case of the WRB, total N deposition accounted for just 10% of N inputs for the entire watershed in 2008 (Compton et al., 2020) and 2.5% of N inputs in the Calapooia watershed, as subbasin of the Willamette River (Lin et al., 2019). Finally, we estimated direct N and P release to waterways in two ways: (1) from human excreta as a function of population density, sewage connection rates, and treatment efficiencies; and (2) for more recent years as a compilation of emissions from wastewater treatment plants (Table 1). Each N and P input to the watershed (and also crop output as harvest) was expressed as kg N or P ha−1 year−1. We also calculated two additional terrestrial metrics: agricultural nutrient balance (J) and nutrient use efficiency (J). We defined annual agricultural N and P balance (J) as where J represents the nutrient N or P, J is nutrient input as inorganic fertilizer on agricultural lands, J is nutrient input as confined animal manure, J is the wet deposition of a nutrient (set to zero for P), J is BNF (also zero for P), and J is nutrient removed in harvested agricultural crops. All fluxes are expressed as per-area rates (kg N or P ha−1 year−1). All values pertaining to agricultural inputs and outputs were calculated on the basis of actual agricultural land but were then rescaled and expressed on a per area basis for the whole watershed to avoid scaling and comparison issues with aquatic loading. Second, we calculated nutrient use efficiency (J) as the percentage of inputs found in crop harvest per year: where variables are as defined in Equation 1.

The area and delineation for the Willamette Basin used to extract the values above were created in ArcGIS using the watershed delineation tool (spatial analyst toolbox, hydrology toolset) and HydroSHEDs 15 arcsecond flow direction maps (Lehner et al., 2008).

Water Quality Data Sources and Processing

We used the R dataRetrieval package (version 2.7.6) to access data on hydrologic flow (cubic feet of water per second on a daily basis), TN, TP, and dissolved inorganic P (DIP) concentrations (mg L−1 usually sampled monthly) from USGS site 14211720 (Willamette River at Portland OR) as far back in time as possible (October 1972) to 2013 from the USGS National Water Information System (Hirsch & De Cicco, 2015). Models are often used to characterize the relationship between nutrient concentrations and water discharge and extrapolate from instantaneously measured nutrient concentrations and loads to seasonal or annual nutrient loads. We used the loadflex R package (Appling et al., 2015) which runs the Weighted Regressions on Time, Discharge and Season—WRTDS model (Hirsch & De Cicco, 2015), currently considered state-of-the-art at USGS (e.g., Oelsner et al., 2017) to simulate long-term and monthly concentrations and loads of TN, TP, and DIP. The WRTDS model R2 values for TN, TP, and DIP are, respectively, 0.86, 0.82, and 0.81.

Data Transformations and Analyses

To assess changes in terrestrial nutrient sources and balances, we first examined source attribution and nutrient balances for 10 target years, spaced at 4- or 5-year intervals, but we also linearly interpolated between these target years. We used Seasonal Mann-Kendall tests to evaluate trends in riverine nutrients (Helsel & Hirsch, 2002). This statistical approach tests for trends in each month over multiple years (e.g., January 1992 load and concentration compared to January 1993) and then combines these trends into a single average. It is a suitable test to look at temporal trends at a single site. All statistics were performed in R 3.2.2 (R Core Team, 2015) and Mann-Kendall tests were done using the rkt package (Saary et al., 2017). Relationships between discharge and concentration were evaluated using linear or nonlinear (polynomial or log-linear) regression, as appropriate based on inspection of the available data.

Because we were interested in comparing terrestrial and riverine nutrient trends over time, we calculated river nutrient fluxes as per-area yields (J; kg N or P ha−1 year−1) as: where J is TN, TP, or DIP load as calculated by the WRTDS model (kg of nutrient year−1) and Area is watershed surface area in hectares.

Having the same units and annual timesteps between terrestrial and riverine nutrients, we could then calculate a nutrient export fraction as a percentage according to Equation 4: where J is an estimate of sewage-derived N or P directly discharged to waterways based on a population and treatment level (see Table 1 and Metson et al., 2017).

We also explored differences between N and P by looking at the ratio of N to P of in riverine export (J) and selected terrestrial nutrient sources (J, J). To do so we simply divided the annual value of N by P of the selected flux and expressed it as the N:P molar ratio. Finally, we looked at cumulative N and P surplus over time versus riverine export as a crude measure of potential N and P retention in the system since 1972 (see SI text). Because we only had eight overlapping target years between terrestrial and riverine nutrients it was not possible to conduct robust statistical tests on these data, in particular none that could account for temporal autocorrelation or multiple factor interaction.

Results

Terrestrial N and P Over Time

Between 1969 and 2012, N and P terrestrial inputs to agricultural lands increased by a factor of 3 and 1.2, respectively (Figures 2a and 2d and Table 2). Synthetic fertilizers were the largest single contributor to terrestrial nutrient inputs and accounted for the majority of observed increases in N and P inputs observed over the 43-year study period. In 1969, synthetic fertilizers accounted for approximately 64% of both N and P inputs, while in 2012 they made up 90% of N and 86% of P inputs. Crop BNF increased more than 100-fold during this period, but never accounted for more than 5.3% (1978) of N inputs (Table 2 and Figure 2a).

Figure 2.

Terrestrial inputs and crop harvest of N and P from 1969 to 2012 where panels (a) and (d) depict agricultural sources, panels (b) and (e) show crop removal as harvest as well as nutrient use efficiency (NUE or PUE), and panels (c) and (f) show agricultural balances (inputs minus crop harvest) and two data sources for nutrients related to human sewage that may be discharged directly to the river. Note that the y axis was kept constant between panels (a) and (b) and (d) and (e) for easy comparison between agricultural inputs and crop harvest, and that the scale changes in panels (c) and (f) in order to accommodate the difference in the magnitude of the nutrient sources related to humans and point sources. NUE and PUE represent the percentage of agricultural nutrient inputs found in crop harvest and refer to the right-hand side y axis in panels (b) and (e). Agricultural inputs and outputs are also available on a per-agricultural-land basis in Table S2.

Table 2

Change in Terrestrial Sources Between 1969 and 2012

	Nitrogen	Phosphorus
Synthetic fertilizer	3.8	1.5
Manure	0.4	0.5
Biological fixation	117	–
Atmospheric wet deposition	1.4	–
Crop harvest	1.7	1.6
Agricultural balance	6	0.9
Human sewage	0.7	0.3

Note. All values are expressed as a multiplication factor (2012 divided by 1969 value) and as such do not account for changes and trends between these years (see Figure 2 and Table S2). Values above 1 indicate an increase and values below 1 indicate a decrease over time while dashes indicate values were not calculated.

Nutrient removal through crop harvest also increased over time, although not always in tandem with increases in inputs (Figures 2b and 2e). In fact, crop harvest increased initially, then leveled off and even appears to slightly decline in later years. For N, crop harvest exceeded nutrient application in 1974, indicating net depletion (mining) of soil N (Figure 2c). From 1974 through 2012, the sum of N inputs in synthetic fertilizer, manure, crop BNF, and wet deposition increasingly exceeded crop harvest. For P, inputs always exceeded crop harvest but the magnitude of this excess varied more than ninefold, from a low of 0.3 kg ha−1 (1982) up to 2.8 kg ha−1 in 2007, with the latter value being a notable outlier (Figure 2f; note that crop harvest values reported in the results are on a watershed area basis but are provided per agricultural area in Table S8). Finally, N and P inputs as human sewage, based on population density and treatment level, agree well with the independent estimate of sewage loading based on known major point sources for overlapping years (MPS, based on the NPDES database, Figures 2c and 2f).

In summary, nutrient use efficiency (percentage of inputs found in crop harvest outputs) has decreased or remained constant over time, while treatment of sewage has improved, decreasing direct nutrient emissions to waterways. For N, the observed decrease in agricultural use efficiency was fairly constant over time, going from 69% in 1969 to 38% in 2012, with a high of 102% in 1974. For P, agricultural use efficiency increased to 89% (1984) but then declined to a low of 43% (2007) before increasing slightly in subsequent years. These relatively low N- and P-use efficiencies should result in soil or ecosystem N and P accumulation unless N and P are lost to river export, lost to groundwater, or (in the case of N) denitrified.

Water Quality Over Time

Monthly riverine concentrations and loads for both N and P decreased between 1972 and 2013 (Figures 3 and S4–S5, Table S3). All variables of interest (TN, TP, and DIP) showed decreases when accounting for seasonal variability (Table 3); all results were statistically significant (p < 0.001) indicating improving water quality. Concentrations of TP and TN both increased significantly with increasing discharge, but with different dependencies on water flow (Figure S2). Whereas TN increased in a generally linear fashion as a function of discharge, TP concentrations exhibited a curvilinear increase with discharge (Figure S2). The nonlinear shape of the discharge-TP relationship appears to result from dilution of dissolved phosphorus (i.e., decreasing concentrations of TP with increasing discharge, even at low flows) combined with a typical exponential increase in particulate P concentrations at higher flows. The net impact of these dynamics on TN:TP ratios was to generate increasing TN:TP ratios at low flows (<~2,000 m3 s−1) but decreasing TN:TP ratios at higher flows (>~2,000 m3 s−1; Figure S2).

Figure 3.

Willamette River basin riverine yields of total N in panel (a), total P in panel (b), and dissolved inorganic P in panel (c) from the WRTDS load estimate model from 1972 to 2013. Gray bars indicate mean annual discharge (cubic feet per second) according to the right-hand axis, highlighting how wetter years coincide with higher yields, especially for TP.

Table 3

Seasonal Mann-Kendall Test Using Monthly River TN, TP, and DIP Concentrations and Loads From 1972 to 2013 From the WRTDS Model

Variable	Tau	P value (2 sided)	Slope
TN concentration	−0.367	0	−0.008
TN load	−0.154	<0.001	−217.738
TP concentration	−0.438	0	−0.0008
TP load	−0.233	<0.001	−26.764
DIP concentration	−0.563	0	−0.0006
DIP load	−0.367	0	−29.131

Note. Tau gives the strength of the trend (varying between +1 and −1) accounting for the number of data points (S/n(n-1)*2, where S is Kendall’s S which is the difference between the concordant and discordant monthly pairs), and the slope value is the median of the differences between the same monthly values in successive years, over all years (Helsel & Frans, 2006; Helsel & Hirsch, 2002).

Comparison Between Terrestrial Sources and Water Quality

Perhaps surprisingly, trends in annual terrestrial inputs and riverine levels were in opposite directions. With terrestrial inputs increasing, and riverine yields decreasing, fractional export of N decreased over time, from 84% to 44% between 1972 and 2012 (Figure 4). There appears to be substantial net retention in the basin and increased net retention over time (Figure S3). The 1974 riverine TN yield was higher than terrestrial inputs resulting in an export efficiency of 102%. There were three other, smaller spikes (1983, 1996, and 2006) interrupting an otherwise fairly constant decline. P export fractions were lower than for N until 1995, varying between 20% and 40%. The same three spikes identified for N export fractions are also visible for P. Unlike the disconnect between TN and TP yields and terrestrial sources, annual DIP yields follow a similar trend to one terrestrial source: human sewage. Sewage could account for 65% of riverine DIP yield between 1972 and 2013 (Figure 5a). There, however, seems to be a reversal in this trend after 2007 where DIP yields increase and human sewage continues to decrease (noting that the 2013 yield estimate goes back down, Table S4). No clear relationship between human sewage and TN and TP exists (Figure 5).

Figure 4.

Fractional riverine nutrient export (river export divided by terrestrial years).inputs) from the Willamette River Basin over time.

Figure 5.

Human sewage (blue bars) and riverine yields (lines) over time for DIP (black line, lower values) and TP yields (light purple line, higher values) in panel (a) and TN yields (gray line) in panel (b). DIP changes over time corresponds much better to human sewage estimates than TP or TN do.

As further evidence of this land-water disconnect, the dramatic increases in the N:P ratios of sewage and fertilizer inputs over time are not mirrored by riverine exports. We observed no clear interannual temporal trends in the river export N:P ratio (Figure 6), but ratios of molar N:P in fertilizer inputs increased dramatically over time from approximately 10 in 1969 to 25 in 2012, with the exception of the mid-2000s with an outlier showing anomalously high P fertilizer use. Sewage N:P ratios (molar) also increased substantially from approximately 10 to 19 over the same time period, with two distinct step-wise increases (1987 and 2007, noting that this could be due to data availability as we linearly interpolated between target years). In other words, terrestrial sources of N increased much more dramatically than P, but no such trend is evident in riverine exports (Figures 2 and 6). Instead, the riverine export N:P ratio appears to be more strongly regulated by flow than by terrestrial inputs. Except for a steady, although modest (from 30 to 26), decline between 1985 and 1993, the N:P ratio varied annually between 15 and 30 without a trend. In 1996, when discharge was high, the riverine N:P ratio was distinctly lower than years before and after (the outlier 15 value in Figure 6). The N:P ratio in the river is higher (relatively more N than P) than both fertilizer and sewage except in 1996 when it dips below the fertilizer ratio and finally in 2012 when it seems to be equal to fertilizer.

Figure 6.

Molar ratio of N to P for riverine export (black dotted line), agricultural fertilizers (red line with triangles) and sewage (blue line with squares) from 1969 to 2012 in the Willamette River Basin overlaid with annual average river discharge (grey bars associated with right hand y-axis values).

Discussion

Water Quality and Nutrient Use Efficiency Over Time

Hydrology was clearly a driver of the interannual variability of riverine nutrient export from the WRB, and the declining trends over time reported here generally agree with previous studies of the watershed (e.g., Fuhrer et al., 1996; Wise, 2007). Spikes in yields for both TP and TN (Figure 3) correspond to particularly wet years. For example, the WRB experienced widespread flooding in February 1996, associated with a large rain-on-snow event in the Cascades (Marks et al., 1998). This flooding affected total nutrient export and caused a large dip in the N:P ratio, indicating that relatively more P than N was exported to the river in 1996 compared to all other years studied (Figure 6). For the WRB, Fuhrer et al. (1996) also identified clear seasonal patterns in nutrient loading, where high winter streamflows corresponded to higher TP concentrations (November to February) and colder temperatures correspond to higher N concentrations (nitrite plus nitrate and ammonia). Most of the N export occurs in fall and winter in the WRB, when sinks within the basin are low and runoff if high (Compton et al., 2020). The strong, exponential relationship between discharge and PP concentration in the Willamette (p < 0.001, R2 = 0.56) noted above (Figure S2) suggests that high-flow events have disproportionate influence on phosphorus loading and transport in this basin, a pattern that is common in other river systems (Harrison et al., 2019). Other regional and national studies relating terrestrial nutrient balances and water quality have identified hydrology as an important driver, including at the national scale for TP (Metson et al., 2017), along the West Coast of the U.S. for TN (Schaefer et al., 2009; Sobota et al., 2009), and for the Northeastern U.S. over time for both TN and TP (Hale et al., 2015). At the national scale, interannual variability in precipitation drives nutrient export over time, but inputs drive patterns in space (Bellmore et al., 2018; Sinha & Michalak, 2016).

Recent studies across the U.S. and Canada have identified similar decreases in nutrient concentrations over time (Keiser & Shapiro, 2018; Shoda et al., 2019; Stammler et al., 2017). Improvements in waste treatment from the 1970s to the present can explain observed declines in TP and, to some extent, TN in streams, particularly in areas with a lot of urban land cover (Hale et al., 2015; Stets et al., 2020). While the WRB is dominated by agriculture, changes in waste treatment sources are important because they are directly supplied to rivers with little opportunity for storage or removal via plant uptake, soil storage, or denitrification. Thus, there is little to no expected lag in the response between changes in point source inputs from wastewater treatment plants and stream concentrations. Reductions in stream DIP, TP, and TN are consistent with improvements in waste treatment technologies over time (Stets et al., 2020); and the close relationship between the decrease of DIP and human sewage inputs shown in Figure 5a has been demonstrated globally (Harrison et al., 2010). In contrast, inputs to, and surpluses on, land, particularly for N, are increasing over time, with no similar response in stream export, suggesting a substantial storage or processing of surplus land-applied nutrients. Storage in soils was predicted to lead to a 35-year lag between inputs and river export in the Mississippi River Basin (Van Meter et al., 2016). Groundwater storage is another possible retention mechanism within watersheds, with decade-long lag times in the WRB (Craner, 2006; Hinkle, 2009). This disconnect strongly suggests a substantial long-term lag in response between land application of N and stream export in this large mixed land-use basin (section 4.3 and SI text).

The decline in NUE observed for WRB crops over time differs from what has been reported for the United States as a whole; early (1960s and 1970s) decreases in NUE have generally preceded a leveling off or increases (Lassaletta et al., 2014; Zhang et al., 2015). PUE is variable but shows no clear trends over time in the WRB. The consistently declining NUE observed in our study is somewhat unexpected, although not unprecedented. Swaney et al. (2018) found that N use efficiencies were lower in the 2002–2012 period compared to 1987–1997 in many regions across the United States, while P use efficiencies increased significantly in a few regions, other regions saw a decrease (Swaney et al., 2019). The Pacific Northwest was not examined in either case. The WRB has shifted from wheat and vegetable forage crops in the 1970s to an increasing amount of grass seed, vegetable seed, and horticultural crops in the 1990s and 2000s (USDA NASS Survey data, Chastain, 2011). Seed crops in particular may be less N efficient in terms of harvested N, because the amount of N harvested is lower than for forage crops. In 2008, NUE was only 41% in a tributary of the WRB dominated by grass seed crops (Lin et al., 2019), consistent with the overall low NUE in the WRB during that time. In addition, the WRB’s Mediterranean climate where much of the rain occurs during the nongrowing season poses significant challenges to improving NUE and PUE while maximizing crop harvest (Compton et al., 2020). Decisions about future land use, and in particular crop mix will likely affect nutrient losses through erosion and groundwater leaching (Berger & Bolte, 2004), although our results do not show tight coupling between crop mix and riverine export in the past. The observed disconnect between increasing terrestrial inputs and declining riverine export appears to be related, at least in part, to this decline in nutrient use efficiency, and this should be examined further.

Export and Retention Compared to Other Watersheds

Riverine export fractions in our study were high compared to other comparable analyses in the U.S. For instance, Sobota et al. (2009) found that, on average, only 8% of terrestrial net N inputs were found in riverine export in selected California watersheds. In 22 larger West Coast watersheds, Schaefer et al. (2009) found that on average 12% of N inputs were exported, with a somewhat higher fraction (30%) exported from the Willamette basin. Our post year-2000 results are similar to those of Schaefer et al. (2009), showing an average export fraction of 30% of N inputs. It is reasonable that the WRB would have higher export fractions than other watersheds because of relatively high runoff. At the same time, our analysis considered fewer terrestrial inputs than the papers above, and as such could also create higher export factions (perhaps missing up to 20% of inputs when compared with Compton et al., 2020). For P, most export fractions in the literature used a different accounting method for terrestrial balances (e.g., Net Anthropogenic Phosphorus Inputs NAPI), which complicates comparisons of results, and P fractional exports vary tremendously both nationally (Metson et al., 2017) and at smaller scales (Sobota et al., 2011). After the year 2000, our P fractional export of net inputs was 53%, a value higher that the majority of the literature values reviewed in Metson et al. (2017).

Our observed export fractions are also high compared to international values, including those with similar Mediterranean climates. For instance, the Ebro River Basin in Northeast Spain only exports around 8% of Net Anthropogenic Nitrogen Inputs (NANI) to the delta; this export efficiency is much lower than that found in other temperate European catchments (Lassaletta et al., 2012). Within the Ebro basin, subcatchments exported on-average 9% of the N inputs using a more similar soil balance approach to ours. Lassaletta et al. (2012) found that catchments with more irrigation channels and a larger area draining to a dam or reservoir retained more N. Export, even if a small fraction of total N inputs, is correlated with discharge where higher discharge results in a higher proportion of N being exported (Howarth et al., 2006, 2012; Mayorga et al., 2010; McCrackin et al., 2014). The WRB, with its relatively high discharge and seasonal climate, where losses occur during the nongrowing season, thus has a relatively high export fraction (Compton et al., 2020).

Although annual export fractions were relatively high, there has likely been significant N and P retention within the WRB over time (Figure S3). In other words, N and P are being processed upstream of the sampling location or being stored. The relative importance of retention mechanisms, for instance agricultural soil storage, riparian retention, in-stream retention, and reservoir retention, vary among watersheds (Billen et al., 2011). Mechanistic models, such as RIVERSTRAHLER can account for a number of these processes in order to more accurately model observed river nutrient loads. Large N retention in riparian areas, as well as benthic P processing in large Northwest European watersheds demonstrates the importance of spatially heterogeneous processes, not only net nutrient budgets, in modeling water quality outcomes (Thieu et al., 2009). We explore retention in the WRB further in section 4.4 by breaking down how different elements not accounted for in our methods could contribute to the observed disconnect between terrestrial and aquatic changes over time, and in the SI by discussing the importance of legacy nutrients in the literature.

Differences Among N and P Trends

Both riverine N and P decreased over time while terrestrial inputs increased. However, the differences between the two nutrients could have their own ecological consequences. The N:P ratio of riverine exports is higher than that of major terrestrial inputs (Figure 6), which indicates that a greater proportion of P is retained in the watershed. However, one must use caution in interpreting these results due to differences in loss/retention pathways for N (e.g., denitrification) and P (e.g., retention behind dams and loss through storm events, section 4.3). The observed increase in fertilizer N:P ratio in the WRB is consistent with reported national (Puckett, 1995; and USGS fertilizer sales reports) and global trends for synthetic fertilizer use (Lu & Tian, 2017; Peñuelas et al., 2012, 2013). Although no single cause has been identified for this difference in mobilization between N and P locally or globally, it could be related, as mentioned above, to N being less readily stored in soil than P. It could also result from changes in crop mixes, policies, and regulations on nutrient inputs, and finally the availability and price of fertilizer products (e.g., globally for P, Cordell et al., 2015; Weber et al., 2014, and changes over time in N management for the EU and USA, van Grinsven et al., 2015). The increasing importance of horticultural crops in the WRB, which are often characterized as being N inefficient (Cameira & Mota, 2017), could help explain the relative increase of N fertilizer use.

The N:P ratio of sewage inputs also increased over the period of record in the WRB, indicating that sewage infrastructure has been better at removing P than N. This is in line with previous explorations of wastewater stoichiometry in the U.S. (Cease et al., 2015; Downing & Mccauley, 1992), and recent work in China showing that, over time, treatment plants have been much better at implementing P retention than N retention, a trend that is evident in lake N:P ratios (Tong et al., 2020). The stoichiometry of wastewater, as well as the solid fraction retained on land, influences not only nutrient limitation in waterways, but also the effectiveness of recycling biosolids and wastewater to fertilizer crops; recycling of the solid faction, relatively higher in P, can result in P overapplication when N is limited on land (Cease et al., 2015; Sileshi et al., 2017).

The average molar N:P ratio of river export in our study (26:1) is higher than the Redfield ratio (16:1 for marine phytoplankton; Redfield, 1958), suggesting P limitation of river phytoplankton growth. River stoichiometry often varies substantially along flow paths from headwaters to the oceans. For instance, in lakes atmospheric N deposition has been found to shift lakes from N to P limitation (Elser et al., 2009), although N deposition may not be as important a contributor to terrestrial nutrient balances in highly agricultural watersheds (Compton et al., 2020; Sabo et al., 2019). Within-waterbody processing can also be important in that increased residence time can elevate N:P ratios as inland P burial can outpace combined denitrification and N burial (Maranger et al., 2018). From a eutrophication and algal bloom perspective, both the total and relative quantities of N and P are important (Anderson et al., 2002; Conley et al., 2009; Smith et al., 2017), suggesting that N and P should be managed and monitored together (Kanter & Brownlie, 2019; Lewis et al., 2011).

Although N:P ratios and trends of terrestrial and riverine fluxes are quite distinct in the WRB, some congruence exists. For instance, relatively high riverine N:P ratio may be explained by fertilizer being a dominant input in the region. In Iowa areas with more synthetic fertilizer use (and fewer animals) have higher N:P ratios in water (Arbuckle & Downing, 2001). In addition, similar patterns to the one we identified in the WRB have been reported for the Chesapeake Bay, where riverine N:P ratios over six decades varied between 17 and 28, and were higher than terrestrial inputs (Hale et al., 2015). These similarities, paired with the distinct effect of hydrology on N and P flows (as explored in section 4.1) and in-water processing (previous paragraph), demonstrate that the stoichiometric patterns observed in this study seem consistent with expectations. Global increases in anthropogenic N mobilization at higher rates than P, not only in terms of fertilizer use, will not only affect crop production but also downstream ecosystem function including species composition and carbon sequestration (Peñuelas et al., 2012, 2013, 2020).

Disconnect Between Terrestrial and Riverine Nutrient Trends Over Time

By examining more than 40 years of data from the WRB, we identified a disconnect between nutrient inputs and river export at the watershed scale. Although we did not expect to see terrestrial and riverine N and P patterns match perfectly, due in part to legacy effects (Chen et al., 2018), the improvement in water quality despite substantial increases in net nutrient inputs was surprising, although not unprecedented. For example, both Hale et al. (2015) for watersheds in the Chesapeake Bay area and Han et al. (2012) for Lake Erie U.S. watersheds have reported disconnects between terrestrial nutrient balances and river nutrient export, where temporal trends in export were better aligned with riverine discharge than terrestrial nutrient budgets. Identifying the cause of such a disconnect is challenging as local biophysical characteristics, water management, and nutrient management all can affect the efficiency with which N and P entering the landscape are transferred to surface waters. Disconnects between terrestrial nutrient dynamics, and riverine nutrient concentrations and loads, can result from (1) retention (soils, groundwater, sediments) but also (2) unaccounted for nutrient sources, biomass uptake, processing, and loss before a sampling station or changes in these pathways. Although retention pools that might later be mobilized (groundwater and soil) are of primary interest, it is important to first disentangle the importance of other possible nutrient sources and loss pathways.

First, there may be nutrient sources not considered here that decreased over time and thus may better correspond with the observed riverine decreases in N and P. For instance, we did not quantify dry atmospheric N deposition. We also did not include total P atmospheric deposition as an input. In some areas (e.g., pristine forests) these flows may have significant effects on (and sometimes complex interactions with) water quality (Amos et al., 2018; Smith et al., 2017; Stoddard et al., 2016). We also did not consider N fixation by red alder trees or nutrient sources from runoff and erosion from non-agricultural lands (e.g., bedrock, forest, urban runoff), which have been shown to be important sources in some parts of the Pacific Northwest U.S. (Compton et al., 2003; Wise & Johnson, 2013). However, we have no reason to believe that such sources are either large enough to have a substantial impact on N or P export or are changing in a systematic way on a multidecadal timescale. Compton et al. (2020) found that N deposition and red alder N fixation accounted for only 10% and 5% of 2008 inputs overall to the WRB, respectively, whereas the agricultural inputs included in our study (fertilizer and manure) comprised 80% of inputs.

Second, there may be output/loss pathways that could decrease the discrepancy between increasing inputs and decreasing nutrient riverine loads. Erosion in storm events can be an important source of nutrients locally, especially for particulate P (e.g., Fraser et al., 1999, and supported by Figure S2). However, a comparison between measured and modeled N and P loads indicated that the WRTDS load-estimating model used in this study are robust for the WRB during high-flow events (i.e., free from any detectable systematic bias). As such it seems unlikely that we have grossly underestimated hydrological losses, although differences in the soil P content (changes over time or related to the location of storm erosion) were not accounted for.

Another output pathway, crop harvest, is somewhat uncertain as the data sources we used did not explicitly account for specialty crops such as grass seed which are important in the WRB (Lin et al., 2019), or for potential changes in nutrient concentrations in plants over time. It is unlikely, however, that we have grossly underestimated crop N and P harvest since the 1980s. Our calculated N crop harvest for the WRB for 2008 (42 kilotons) was within the range of estimated crop removal (25–52 kilotons) given locally derived 2007 land use and crop removal that included specialty crops (Compton et al., 2020; Lin et al., 2019). In other words, the IPNI “other crop” category reasonably represented the nutrient demand of specialty crops. The 10 most important agricultural crops (in acreage) in the Willamette valley in 1987 included grass seed (most acres), and mint (Wentz et al., 1998), and thus it is also unlikely that the IPNI “other crop” category was invalid for earlier estimated years. The increase in specialty crops over time could mean that the decline in nutrient use efficiency over time was less steep than observed (Figures 2b and 2e) if we underestimated the contribution of such crops to harvested N and P. However, because PUE and NUE seem to be low in these specialty systems (sections 4.1 and 4.3), and contemporary fertilizer use in the region is higher than crop recommendations in the WRB (Compton et al., 2020; Lin et al., 2019), it seems unlikely that accounting for more crop harvest would change the observed pattern of increased agricultural N and P surpluses.

Finally, for N, gaseous losses in both terrestrial and aquatic environments likely play a role, but it is difficult to partition how much N is lost to the atmosphere versus stored in soils at the scale of the whole WRB. A substantial portion of the “stored” N also could be denitrified and lost to the atmosphere as N2 or other N-based trace gases (N2O, a potent greenhouse gas or NO, a regional air pollutant). For instance, Horwath et al. (1998) estimated denitrification could remove up to 12.5% of applied fertilizer N in grass cropping systems in the WRB, although rates of denitrification can be much smaller and vary among crops, season, and soil conditions (Myrold, 1988). Additional gaseous losses at other stages would likely have a smaller impact on our estimates but are still worth mentioning as they could be locally important. For instance, denitrification removed 0–6.8% of nitrate moving through southern Willamette Valley streams (Sobota et al., 2012), while denitrification in the hyporheic zone has also been shown to mediate N fluxes between groundwater and river water in the WRB (Hinkle et al., 2001). Davis et al. (2007) found that denitrification was significant in poorly drained riparian areas in the southern Willamette, and that denitrification rates led to low nitrate concentrations in associated riparian wells. Finally, Rupert (2008) observed groundwater nitrate concentrations to have declined in their decadal-scale monitoring at a small number of sites in the WRB, which they partially attribute to changed reduced redox conditions in study wells, favoring denitrification.

Other lines of evidence support the notion that soil and groundwater are important stores of N. For instance, there are expansive areas in the WRB where groundwater nitrate levels exceed the drinking water standard and are correlated with agricultural N use (Hoppe et al., 2014). In fact, the discrepancy between Rupert (2008) findings of declining groundwater nitrate and the Hoppe et al. (2014) findings of high groundwater nitrate may be explained by the fact that, based on our calculations, net N accumulation occurred only after 2010 and thus that regional increases in groundwater nitrate may only be beginning to be observed (Figure S3). Piscitelli (2019) found nitrate in wells to have increased between 2006 and 2018, and especially after 2012, in the southern Willamette Valley Groundwater Management Area. Finally, agricultural soil is likely an important contributor to basin N retention. Soil N recovery was approximately 31% after 14 months in a 15N tracer study in grass seed crops in the southern Willamette Valley (Davis et al., 2006), supporting the idea that a substantial amount of N is stored in soils. Studies outside the WRB where longer-term soil data were available also indicate that N storage is occurring (Hirsh & Weil, 2019; Van Meter et al., 2016). In summary, although we have likely underestimated the role of gaseous losses in our estimate of terrestrial N sources being retained over time, the general trend of increasing agricultural surpluses and thus increased retention over time is likely correct (Figures 2 and S3).

Management Implications and Future Directions

It is increasingly recognized that examining and understanding historic trends in both terrestrial and aquatic nutrient flows can substantively inform meaningful interventions to optimize nutrient use for food production and other human demands while minimizing aquatic pollution (Hermann et al., 2018; Keil et al., 2018; Withers et al., 2017). In particular, understanding the relationships among drivers of N and P management, terrestrial flows, and subsequent impacts on water quality that include time-lagged responses is important for establishing realistic targets for water quality improvement. That is, if there is a lag between interventions and water quality outcomes, it may be particularly important to set targets that are realistic in order to ensure that (1) decision-makers and their political constituencies are not discouraged from their efforts, and (2) that support for interventions does not cease midway through because water quality does not improve quickly. Multi-year, and even multi-decade, legacy effects must be taken into consideration (Chen et al., 2018; McCrackin et al., 2018). Although there has been some success in improving water quality across the U.S., in particular related to point sources in urban watersheds, agricultural watersheds remain a challenge for both N and P management (Stets et al., 2020). Consideration of both N and P is critical for watershed management and preventing eutrophication. For the WRB, although we identified improving trends in riverine N and P export, continued increases in terrestrial inputs could be problematic for groundwater and soils, and for future loading responses. Although there may be multiple causes of the observed decoupling, including effective management, it is also possible that the tipping point for landscape buffering capacity has not yet been reached. If this is the case, continued improvements in wastewater treatment and effective agricultural nutrient management to reduce surpluses will be essential to avoid passing this point as it may be difficult to improve water quality quickly after this point has been passed. Across the U.S., including within the Pacific Northwest, improvements in terrestrial balances did not consistently translate to improved riverine P export (Stackpoole et al., 2019). Lessons learned from the WRB are particularly important for areas with similar climates, urban infrastructure, and agricultural production systems such as those in Southern British Columbia, Canada, where surface water quality impairment is not yet severe but it is clear that terrestrial sources are increasing, and agricultural surpluses exist (Bittman et al., 2017; Schindler et al., 2006; Smukler et al., 2014).

Moving forward, research should include improved quantification of accumulation and removal pathways for both N and P. This could include incorporating more observational data such as groundwater N and P concentrations, soil N and P storage over time, denitrification rates and site-specific information on crops and manure and fertilizer management, as well as other N and P inputs such as atmospheric N and P and alder N fixation. By creating more complete and accurate temporal nutrient balances, nutrient management and policy can be refined and improved over time.