Subscriber access provided by TUFTS UNIV
Article
Precipitation dominates interannual variability of riverine nitrogen loading across the continental United States Eva Sinha, and Anna M Michalak Environ. Sci. Technol., Just Accepted Manuscript • DOI: 10.1021/acs.est.6b04455 • Publication Date (Web): 22 Oct 2016 Downloaded from http://pubs.acs.org on October 24, 2016
Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.
Environmental Science & Technology is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.
Page 1 of 34
Environmental Science & Technology
1
Precipitation dominates interannual variability of riverine
2
nitrogen loading across the continental United States
3 4
E. Sinha1,2,* and A. M. Michalak1,2
5
[1]{Department of Earth System Science, Stanford University, Stanford, CA, USA}
6
[2]{Department of Global Ecology, Carnegie Institution for Science, Stanford, CA, USA}
7
* Correspondence to: E. Sinha (esinha@stanford.edu)
8 9
Abstract
10
Excessive nitrogen loading to waterways leads to increased eutrophication and associated
11
water quality impacts. An understanding of the regional and interannual variability in nitrogen
12
loading and associated drivers is necessary for the design of effective management strategies.
13
Here we develop a parsimonious empirical model based on net anthropogenic nitrogen input,
14
precipitation, and land use that explains 68% of the observed variability in annual total
15
nitrogen flux (QTN) (76% of ln(QTN)) across 242 catchment-years. We use this model to
16
present the first spatially- and temporally-resolved TN flux estimates for all eight-digit
17
hydrologic unit (HUC8) watersheds within the continental United States (CONUS), focusing
18
on the period 1987-2007. Results reveal high spatial and temporal variability in loading, with
19
spatial variability primarily driven by nitrogen inputs, but with interannual variability and the
20
occurrence of extremes dominated by precipitation across over three quarters of the CONUS.
21
High interannual variability and its correlation with precipitation persist at large aggregated
22
scales. These findings point to a fundamental challenge in managing regions with high
1
ACS Paragon Plus Environment
Environmental Science & Technology
Page 2 of 34
23
nutrient loading, because these regions also exhibit the strongest interannual variability and
24
because the impact of changes in management practices will be modulated by meteorological
25
variability and climatic trends.
26 27
1
Introduction
28
Human actions, primarily through fertilizer addition, fossil fuel combustion, and
29
increased cultivation of legumes, have more than doubled the amount of reactive nitrogen in
30
terrestrial ecosystems, thus altering the global nitrogen cycle1. In the United States (US),
31
reactive nitrogen from anthropogenic sources is four times that from natural sources2,3.
32
This over enrichment of nitrogen and resulting eutrophication have contributed to
33
impacts including harmful algal blooms (HAB)4 and hypoxia5,6. Within the US, excess
34
nitrogen has led to significant increases in the occurrence and severity of coastal harmful algal
35
blooms and hypoxia6–8.
36
Environmental Protection Agency (USEPA) concluded that 28% of streams and 20% of lakes
37
nationwide experience high levels of nitrogen9,10.
For freshwater systems, recent surveys conducted by the US
38
Developing an understanding of the water quality impacts of human activity is
39
predicated on the availability of robust estimates of nutrient loading to specific impacted
40
systems. Such loading estimates should ideally be available at spatial11 and temporal scales
41
that are relevant for the development of effective management strategies and also provide full
42
spatial coverage to allow for systematic comparisons between regions. Although nutrient
43
loading is routinely monitored in some regions in the US, most watersheds have limited or no
44
water quality data for estimating nitrogen load12. Model-based estimates of nutrient loading,
45
on the other hand, are often available only for long term hydrologically-average conditions13–
46
15
, precluding their use in furthering an understanding of interannual variability in loading and 2
ACS Paragon Plus Environment
Page 3 of 34
Environmental Science & Technology
47
downstream impacts. When estimates for specific years are available, these are typically
48
reported only for very limited regions and periods16,17.
49
Robust estimates of nutrient loading would enable the development of an
50
understanding of the primary factors driving the spatial and temporal variability in loading, as
51
well as an understanding of the occurrence and causes of extreme loading events. For the
52
continental US (CONUS), drivers of the spatial variability of total nitrogen (TN) flux have
53
been examined in some previous regional-scale studies, and include net anthropogenic
54
nitrogen
55
characteristics14,22 and precipitation19. The drivers of interannual variability in TN loading
56
have only been quantified for select watersheds within the Mississippi river basin16 and the
57
Lake Michigan basin17 and include annual stream discharge and annual mean precipitation.
58
Field-based studies have also reported precipitation to be the driver of higher loading during
59
wet years as compared to dry years23,24. The incidence of loading extremes has not been
60
formally examined, although anecdotal evidence and field-based studies have reported
61
linkages between factors such as extremes in discharge and nitrogen loading25. Overall,
62
although these earlier studies provide hypotheses about the primary factors driving TN
63
loading, they are not sufficient to provide a comprehensive view.
input
(NANI)17–20,
population14,16,
stream
discharge13,16,17,20,21,
land
use
64
To address these needs, here we develop an empirical statistical model for estimating
65
TN load from catchments within the continental US (CONUS). We then apply the model to
66
all HUC8 watersheds within the CONUS for a range of years. The goals of this analysis are to
67
(1) identify factors that explain observed variability in TN flux, (2) estimate TN flux for all
68
HUC8 watersheds within the CONUS for 1987-2007, and (3) identify the primary drivers of
69
spatial and temporal variability, as well as extremes, in TN flux across the CONUS. The last
70
goal is addressed at both HUC8 and aggregated spatial scales.
3
ACS Paragon Plus Environment
Environmental Science & Technology
71
2
Page 4 of 34
Methods
72
We use water quality observations collected at select USGS gages, together with data on
73
nitrogen inputs, precipitation, temperature, tile drainage and land use to develop an empirical
74
statistical model for quantifying annual TN loads, and use this model to estimate loading for
75
HUC8 watersheds throughout the CONUS.
76
development and application are described in Section 2.1. The statistical model development
77
is described in Section 2.2, and its application is outlined in Section 2.3.
78
2.1
The observational datasets used for model
Datasets used
79
The response variable for the empirical statistical model developed in this work is
80
annual TN flux per unit area (QTN) [kg-N/km2 yr]. Observations of non-flow-normalized QTN
81
were obtained by applying the Weighted Regressions on Time, Discharge and Season
82
(WRTDS) method26 to TN concentration measurements from gages with long-term
83
observations within the USGS NWIS (National Water Information System)27. The WRTDS
84
method estimates long-term records of water quality parameters through weighted regressions
85
of concentrations on time, discharge, and seasons26. One of the documented limitations of the
86
WRTDS method is possible under prediction of loading for a wet year following a dry year28;
87
however, in the absence of daily TN load measurements even for routinely monitored gages,
88
the WRTDS method is the state-of-the-art approach for estimating loading. Daily TN loads
89
[kg-N] were obtained by multiplying estimated daily mean concentrations by daily discharge,
90
and QTN was estimated by adding daily load over an entire year and dividing by the catchment
91
area. We focus on calendar-year loads (Jan-Dec) rather than water-year loads (Oct-Sept) to
92
coincide with the time periods represented by the NANI data used in the analysis, but conduct
93
a parallel water-year-based analysis as a sensitivity test (see SI).
94
The USGS Geospatial Attributes of Gages for Evaluating Streamflow Version II 4
ACS Paragon Plus Environment
Page 5 of 34
Environmental Science & Technology
95
(GAGES-II) database was used to identify USGS stream gages where discharge and water
96
quality measurements are routinely conducted. The GAGES-II database provides geospatial
97
characteristics for 9,322 stream gages29. Of these, only 72 gages have a minimum of 20 years
98
of uninterrupted daily discharge measurements within the period 1981-2010 as well as a
99
minimum of 200 TN concentration measurements during those 20 years, which are the
100
minimum criteria for applying the WRTDS approach. The 1981-2010 period was selected
101
because the target years for the model build were 1987, 1992, 1997, 2002, and 2007, selected
102
based on the availability of the most complete NANI data. For four gages, on the Raisin,
103
Maumee, Sandusky, and Cuyahoga rivers, more extensive TN data from the Heidelberg
104
University National Center for Water Quality Research (NCWQR), Tributary Loading
105
website (http://www.heidelberg.edu/academiclife/distinctive/ncwqr/data) were substituted for
106
the TN data from the USGS NWIS (National Water Information System) data services that
107
were otherwise used within the WRTDS method.
108
Estimates were extracted for the model build years after applying the WRTDS method,
109
but years with fewer than six TN concentration samples were not used in order to eliminate
110
years with few observations. Using this criterion, the total number of catchments available for
111
the model build was reduced to 70, representing a total of 242 QTN observations across the
112
five target years (1987, 1992, 1997, 2002, and 2007).
113
Net anthropogenic nitrogen input (NANI) was used to represent nitrogen inputs for the
114
selected years and watersheds, and was defined as the sum of five components: fertilizer N
115
input, atmospheric deposition, agricultural nitrogen fixation, net food and feed import and
116
non-food crop export. Estimates of the last three components are only available for the
117
agricultural census years (1987, 1992, 1997, 2002, and 2007)30, while fertilizer data are
118
available for 1987-200631 and atmospheric deposition data are available since 197832. The
119
model development was therefore based on the five agricultural census years. Fertilizer usage 5
ACS Paragon Plus Environment
Environmental Science & Technology
Page 6 of 34
120
for 2007 was estimated as the average of 2003-2006. The NANI components were estimated
121
for the five census years using the NANI toolbox30, except that atmospheric deposition was
122
instead estimated from the National Atmospheric Deposition Program (NADP)
123
network using the approach described in Ruddy et al.33 (see details in SI). The county-level
124
NANI data were aggregated to the 70 catchments using an area-weighted average, while the
125
site-level atmospheric deposition observations were interpolated using inverse-distance
126
weighted interpolation and then aggregated to the catchment areas. In order to account for
127
loading resulting from residual NANI from previous years, we also considered NANI for
128
preceding years in the model development. Those components of NANI that are only
129
available for agricultural census years were estimated via linear interpolation for intervening
130
years. Examples for four representative HUC8 watersheds are available in Figure S1.
32
monitoring
131
Model development focused on ancillary variables with contiguous spatial and
132
continuous temporal coverage across the CONUS in order to enable application throughout
133
the CONUS. Monthly and daily precipitation and temperature data were obtained from the
134
PRISM database34, which provides gridded precipitation and temperature at a 4km grid
135
resolution. Average precipitation or temperature associated with the catchments upstream of
136
the selected stream gages were obtained by taking an area-weighted average of the
137
precipitation or temperature falling within the gage’s catchment. The percentage tile drained
138
areas within each catchment was estimated based on county level tile drainage extent
139
estimates35, rescaled to the catchment scale through area-weighted average. Catchment
140
boundaries for the selected USGS gages and associated land cover data were obtained directly
141
from the GAGES-II database.
6
ACS Paragon Plus Environment
Page 7 of 34
142
Environmental Science & Technology
2.2
Model development
143
We used variables related to NANI, precipitation, land use, temperature, and tile
144
drainage to develop a linear model for natural log transformed annual TN flux (ln(QTN)),
145
calibrated using the 242 catchment-years of QTN observations: = +
(1)
146
where y (m×1) represents observed ln(QTN), X (m×k) is a matrix of predictor variables and a
147
column of ones representing the intercept term, β (k×1) is a vector of regression coefficients
148
and the intercept, ε (m×1) is a vector of residuals, and m=242. Estimates of the regression
149
coefficients, prediction uncertainties and prediction intervals for ln(QTN) were obtained using
150
standard tools as described in the SI.
151
A statistical model selection approach based on the Bayesian Information Criterion
152
(BIC)36 was used to select predictor variables to be included in the linear model, based on
153
their ability to explain the variability in observed ln(QTN). The BIC measures model fit using
154
the Residual Sum of Squares (RSS), and model complexity is penalized based on the number
155
of predictor variables in a model. BIC is defined as: = −2 ⋅
− − + ⋅
(2)
156
where ŷ (m×1) represents predicted ln(QTN) based on a particular subset of predictor variables
157
and k represents the number of auxiliary variables including the intercept in the linear model.
158
Relative to other criterion-based model selection approaches (e.g. the Akaike Information
159
Criterion (AIC)), BIC has a larger penalty term and therefore tends to select a smaller (i.e.
160
more conservative) set of predictor variables. All possible subsets of candidate predictor
161
variables were considered, except when noted otherwise below, and the subset with the lowest
162
BIC was identified as the best model, in the sense that it provided the optimal balance 7
ACS Paragon Plus Environment
Environmental Science & Technology
Page 8 of 34
163
between explanatory power and complexity. The candidate predictor variables, based on
164
NANI, precipitation, land use, temperature, and tile drainage are described below and listed in
165
Table S1.
166
Two functional forms were considered for NANI. The first is simply NANI per unit
167
area for a given year and catchment, while the second is the inverse hyperbolic sine of NANI:
#$# #$# #$# ( ' = asinh " % = & + " % + 1* 2 2 2
(3)
168
The transformed NANI variable was considered because NANI varies by several orders of
169
magnitude and can also take on negative values30. The inverse hyperbolic sine is defined for
170
all real numbers, including negative values and zero, and becomes almost identical to the
171
natural log transformation for values greater than one (Figure S2). The linear model was
172
restricted to selecting either NANI or fNANI for the current loading year as a predictor variable,
173
or neither, but not both. When NANI or fNANI for the current loading year was selected, the
174
model selection was set up to also allow the addition of NANI for one or two preceding years,
175
to represent residual NANI that can accumulate and contribute to nitrogen loss37,38, as long as
176
the sign of the regression coefficients was consistent for all NANI terms.
177
Fifteen candidate predictor variables were based on precipitation, defined based on
178
climate change indices developed by the Expert Team on Climate Change Detection and
179
Indices (ETCCDI) (http://etccdi.pacificclimate.org/list_27_indices.shtml). These included
180
total annual precipitation; total precipitation during the months of March, April and May; and
181
thirteen candidate variables indicative of extreme precipitation. The extreme precipitation
182
variables included the number of days with extreme precipitation annually and during the
183
months of March, April, and May and extreme precipitation amounts annually and during the
184
months of March, April, and May. These variables are described in further detail in Table S1. 8
ACS Paragon Plus Environment
Page 9 of 34
Environmental Science & Technology
185
The seasonal and extreme seasonal precipitation variables focused on the months of March,
186
April and May because the catchments used for model development had the highest monthly
187
TN flux and highest monthly extreme TN flux, defined as flux above 75th, 90th, 95th and 99th
188
quantiles, during these three months. Several of the extreme precipitation variables are
189
strongly collinear, and hence we restricted the final model to include at most one extreme
190
precipitation variable.
191
Two candidate variables were based on temperature, namely average annual
192
temperature and average temperature during the months of March, April and May. The linear
193
model was restricted to selecting at most one temperature variable.
194
Thirty candidate predictor variables were based on land use. Land use was considered
195
because nitrogen loss from predominately agricultural watersheds has been shown to be
196
significantly higher than that from natural or forested catchments39,40. We defined five land
197
use categories: urban, agriculture, forest, wetlands and shrubland & herbaceous, obtained by
198
aggregating the thirteen land use categories in the National Land Cover Database 2006
199
(NLCD2006). All five land use categories, described in Table S2, as well as all their
200
combinations, as described in Table S1, were included as candidate predictor variables. Any
201
single land use category was only allowed to be represented once in the model, either in the
202
individual or binned categories, and a maximum of four-land use categories could be
203
represented either in the single or binned categories. Models with all five-land use categories
204
were not considered to avoid perfect multicollinearity (where one variable can be predicted
205
exactly from the other variables).
206 207
A single candidate variable was based on tile drainage, representing the percentage of catchment area that is tile drained.
9
ACS Paragon Plus Environment
Environmental Science & Technology
208
2.3
Page 10 of 34
Model application
209
The developed statistical model was applied to all HUC8 watersheds within the
210
Continental U.S. (CONUS) (Figure 1), except for ten HUC8 watersheds that are almost
211
entirely made up of water (e.g. each of the Great Lakes), for the period 1987-2007. The
212
watershed boundaries were obtained from the USGS Watershed Boundary Dataset (WBD)
213
(http://datagateway.nrcs.usda.gov). The HUC8 watersheds were selected for model
214
application, because the range of watershed area for HUC8 watersheds (25th quantile–2,300
215
km2; median–3,300 km2; 75th quantile–4,900 km2) is comparable to the area range of the
216
catchments used for model build (25th quantile–225 km2; median–1,500 km2; 75th quantile–
217
4,000 km2). The model application was limited to 1987-2007 due to the unavailability of
218
many components of NANI before 1987 and after 2007. The TN flux value estimate was
219
calculated by taking an antilog of ln(QTN) and multiplying by the bias correction factor
220
exp(+,( /2) where +, is the standard error of the ln(QTN) residuals from the statistical model.
221
The annual TN load from a watershed was obtained by multiplying QTN [kg-N/km2 yr)] by the
222
watershed area [km2]. Additionally, annual TN loads [kg-N/yr] were aggregated for the whole
223
of CONUS and for regional watersheds associated with major nutrient delivery points to the
224
coastal ocean (Figure 1). Numbers reported at these aggregated scales represent estimates of
225
total nitrogen loading within given regions, rather than estimates of nitrogen export from the
226
regions, which would require an additional assessment of in-stream nitrogen loss for these
227
larger scales.
228
Land use classification for the HUC8 watershed was obtained from the NLCD 2006
229
database41 and was kept constant over the examined period because land use change has been
230
minimal in the CONUS over the study period42,43. Precipitation, temperature, tile drainage,
231
and NANI were obtained and spatially aggregated using the same approach as for the
10
ACS Paragon Plus Environment
Page 11 of 34
Environmental Science & Technology
232
catchments used in the model development. For HUC8 watersheds spanning the border with
233
Canada or Mexico, all variables were truncated at the border to reflect only the CONUS
234
portion of the loading; this choice was made based on data availability. For the NANI
235
components that are only available for agriculture census years (i.e., agricultural nitrogen
236
fixation, net food and feed import and non-food crop export) estimates for intervening years
237
were obtained through linear interpolation (see example in Figure S1). The remaining two
238
NANI components, namely atmospheric deposition and fertilizer usage, are available annually
239
for the entire examined period, except for fertilizer usage specifically for 2007, as already
240
noted in Section 2.1.
241
3
242
3.1
Results & discussions
Statistical model
243
A model based on five predictor variables explains 76% of the variability in the
244
observed ln(QTN) [ln(kg-N/km2 yr)] across five years and 70 catchments in the CONUS
245
(Figure S3): -. = 3.01 + 0.3742 ⋅ + 0.0014 ⋅ 455678 + 0.0033 ⋅ 499,;? −0.0529 ⋅ BCD − 0.0220 ⋅ BCE,FG
(4)
246
where fNANI [ln(kg-N/km2 yr)] is the transformed annual NANI for the current loading year
247
(Eqn. 3); PAnnual [mm] is annual precipitation in the catchment; PMAM,p>0.95 [mm] is extreme
248
precipitation expressed as the amount of precipitation that fell in March, April, and May on
249
days with precipitation greater than the 95th percentile (where the percentiles are calculated
250
based on daily precipitation for 1981-2010); LUW [%] is the percentage of the catchment area
251
classified as wetland; and LUF,SH [%] is the percentage of catchment area classified as forest
252
or as shrubland & herbaceous. The significance of each variable in the final model was further
11
ACS Paragon Plus Environment
Environmental Science & Technology
Page 12 of 34
253
tested by removing each variable from the model and conducting an F-test (p < 0.001 for all
254
variables). A 10-fold cross validation using the same five predictor variables yielded a root
255
mean squared error of 0.63 [ln(kg-N/km2 yr)] and a cross-validation R2 of 0.75. The validity
256
of the assumption of a linear relationship between ln(QTN) and the predictor variables was
257
qualitatively confirmed by examining scatterplots between ln(QTN) residuals from the final
258
multiple linear regression model with one variable removed and this individual auxiliary
259
variable (Figure S4b).
260
The model not only captures the spatio-temporal variability of the observed ln(QTN)
261
across the 242 catchment-years, but is also highly predictive of loading for catchments falling
262
within specific regions (e.g. MARB, Columbia River Basin) of the CONUS (Figure S5). This
263
implies that the variables included in the model are able to capture many of the differentiating
264
features of regions across the CONUS.
265
Among the selected variables, fNANI is the strongest single predictor of ln(QTN) across
266
the examined catchments and years, explaining 60% of the variability in ln(QTN) (Figure S4a).
267
NANI in the untransformed space also explains 46% of the variability in QTN (figure not
268
shown). The importance of anthropogenic nitrogen input in impacting nitrogen export from
269
watersheds within the CONUS is well known
270
reported to explain between a quarter and three quarters of the observed spatial-only
271
variability in TN flux averaged over multiple years for select watersheds within the
272
CONUS45, within the North Atlantic Basin18, and within the northeastern US19.
273
individual years within the Lake Michigan Basin, Han et al.17 found that NANI explained
274
between 69% and 91% of the spatial variability of TN export. Whereas most earlier studies,
275
have primarily focused on exploring the influence of nitrogen input within specific regions,
276
individual years, or long-term mean fluxes, we show here that NANI explains over half of the
13,14,16–21,44,45
. NANI has previously been
For
12
ACS Paragon Plus Environment
Page 13 of 34
Environmental Science & Technology
277
overall space and time variability in ln(QTN) in catchments across the CONUS and across
278
multiple years.
279
Total annual precipitation alone explains 8% of the observed variability in ln(QTN)
280
across the examined catchments and years but 14% of the observed variability that cannot be
281
explained by the other variables in the model (Figure S4). The inclusion of PAnnual in the final
282
model confirms that the role of annual precipitation observed in earlier studies focusing on
283
limited regions and time periods applies to catchments across the CONUS and across multiple
284
years.
285
precipitation and nitrate loss in highly fertilized soils in the Mississippi Atchafalaya River
286
Basin (MARB) using a process-based model. Similarly, Howarth et al.19 observed a positive
287
correlation between average annual precipitation and the fraction of NANI exported for 16
288
watersheds in the Northeastern US. More broadly, several studies have noted the relationship
289
between runoff or discharge and the spatial or temporal variability in nitrogen flux for specific
290
regions within the CONUS13,16,17,20,21. Quantitatively, the role of precipitation found here is
291
comparable to the 11% of spatial and temporal variability16 and 15% of spatial variability21
292
explained by discharge in studies of select monitoring stations within the MARB. We show
293
here that PAnnual explains a significant portion of both the space and time variability in ln(QTN)
294
for a wide variety of catchments across the CONUS.
For example, Donner et al.44 identified a positive correlation between annual
295
The selection of extreme precipitation in the months of March, April and May
296
(PMAM,p>0.95) demonstrates, for the first time, the influence of extreme precipitation on annual
297
nitrogen loading across a large variety of catchments and over multiple years. The strength of
298
this influence is similar to that of PAnnual, in that PMAM,p>0.95 alone explains 7% of the observed
299
variability in ln(QTN) across the examined years and catchments. Empirically, the importance
300
of extreme precipitation in the transport of nitrogen flux has been discussed in a few earlier
301
studies such as Royer et al.25, who focused on three watersheds in Illinois. Anecdotally, it is 13
ACS Paragon Plus Environment
Environmental Science & Technology
Page 14 of 34
302
also known that extreme precipitation leads to higher nutrient loading46. In addition, seasonal
303
precipitation in the months of March, April, and May (PMAM) has been shown to be a good
304
predictor of May-June nitrate flux for the MARB47. The model developed here shows that
305
springtime extreme precipitation is an important factor for explaining variability in nitrogen
306
flux across a large variety of catchments within the CONUS.
307
The fact that two land use variables were selected in the final model indicates that the
308
fraction of land use defined as forest or shrubland & herbaceous (LUF,SH) and as wetland
309
(LUW) provides additional explanatory power beyond their covariation with NANI and
310
precipitation. The negative drift coefficients associated with these two variables indicate that
311
increases in these variables are predicted to decrease the ln(QTN). This is likely due to the fact
312
that forests better retain nitrogen entering via atmospheric deposition relative to other land use
313
types48,49, while wetlands have high denitrification rates that also lead to reduced nitrogen
314
export to streams50. Qualitatively, our results are also consistent with a study conducted by
315
Goolsby and Battaglin14 on several watersheds within the MARB that identified agricultural
316
land as being associated with higher multiyear average riverine nitrogen flux even after
317
accounting for nitrogen inputs to watersheds. Similarly, Alexander et al.51 estimated the
318
highest nitrogen fluxes from agricultural lands and the lowest from forests and shrublands in a
319
study examining sources and transport of nitrogen in the MARB using a process-based model.
320
Our study generalizes these observations, by showing a decrease in ln(QTN) associated with
321
specific land use types (LUF,SH and LUW) for catchments across the CONUS.
322 323
NANI for preceding years, temperature, and tile drainage did not improve model fit sufficiently to warrant inclusion in the final model (see details in SI).
14
ACS Paragon Plus Environment
Page 15 of 34
324
Environmental Science & Technology
3.2
Model application to all HUC8 watersheds
325
We apply the developed model to predict TN flux for each year within 1987-2007 at
326
the HUC8 watershed scale for the entire CONUS (Figure 2). The authors are aware of only
327
two studies that had previously provided spatially-explicit TN flux estimates for the
328
CONUS15,51, but these provided no interannual information. Both of these earlier studies used
329
the SPARROW model run for hydrologically-average conditions and using a single
330
representative year of NANI. The TN fluxes estimated here are based on year-specific
331
precipitation and NANI observations, and are derived using a model that has been calibrated
332
based on observed TN flux across 70 catchments and five years. In the paragraphs that follow,
333
the term “region” refers to either HUC2 watersheds or regions corresponding to major
334
nutrient deliver points, both as defined in Figure 1.
335
We find that watersheds within the Upper Mississippi, Ohio, Tennessee, and Lower
336
Mississippi regions have more than twice the average nitrogen loading per unit area relative to
337
the whole of CONUS (Figure 2a), and together account for 39% of CONUS loading even
338
though they represent only 15% of CONUS area. This result is consistent with Robertson et
339
al.52, who identified very similar high-loading areas in a study examining nitrogen yield for
340
the MARB region under average hydrologic conditions. The western portion of the Pacific
341
Northwest region also experiences extremely large nitrogen loading per unit area (Figure 2a).
342
Overall, areas with the largest temporally average QTN are those with the highest inputs of
343
NANI and small percentage of land classified as forest, shrubland, herbaceous, or wetlands
344
(Figure S6a, S6d and S6e).
345
For the first time, we are also able to quantify the interannual variability associated
346
with nitrogen loading for all HUC8 watersheds within the CONUS and find that the year-to-
347
year variability in QTN is highest for areas (Figure 2b) with the highest overall loading (Figure
348
2a). In fact, the standard deviation of the interannual variability in TN flux is greater than 15
ACS Paragon Plus Environment
Environmental Science & Technology
Page 16 of 34
349
50% of the mean for HUC8 watersheds within high loading regions, including the Lower
350
Mississippi region, the eastern portion of the New England region, the western portion of the
351
Pacific Northwest region, and the northern portion of the California region (Figure 2c).
352
We find that NANI is the strongest predictor of the spatial variability of nitrogen
353
loading across the CONUS, consistent with earlier studies focusing on spatial variability of
354
nitrogen flux for specific areas within the CONUS17–19,45. Average NANI over the 21 years
355
examined here explains 50% of the spatial variability in the estimated average QTN across all
356
HUC8 watersheds, and transformed NANI (fNANI) is the leading term in predicting ln(QTN) for
357
81% of HUC8 watersheds (Figure 3a). The leading term was determined for each watershed
358
by averaging the absolute value of the contribution of each variable in Eq. 4 across the 21
359
examined years. For many of the watersheds with very low or negative NANI (Figure S6a),
360
the percentage of land cover by forests and shrubs becomes the leading terms (16% of
361
watersheds). Watersheds with negative NANI are those where exports of food and feed
362
products outweigh other NANI contributions30 (e.g. in portions of the Missouri, Lower
363
Colorado, and Rio Grande regions). PAnnual was the leading term in only a small fraction of
364
watersheds (2%) located in the western portion of the Pacific Northwest region, the wettest
365
area of the CONUS overall (Figure S6b). Overall, spatial patterns of NANI are the strongest
366
determinant of temporally-averaged nitrogen loading across the CONUS.
367
For interannual variability, however, precipitation becomes the primary driver of
368
nitrogen loading for watersheds across the CONUS. Annual precipitation is the primary driver
369
of interannual variability in QTN for 62% of watersheds (these watersheds contribute 62% of
370
TN load for the CONUS), and extreme precipitation in the months of March, April and May
371
is the primary driver for 24% of watersheds (which represent 36% of the CONUS TN load),
372
while NANI is the primary driver for only 14% of watersheds (which represent 2% of the
16
ACS Paragon Plus Environment
Page 17 of 34
Environmental Science & Technology
373
CONUS TN load) (Figure 3b). The primary drivers of interannual variability were identified
374
by calculating the correlation between estimated annual QTN and QTN that would have been
375
estimated by keeping all but one predictor variable at their mean value over the 21 years. The
376
predictor variable that leads to the largest correlation between the two sets of fluxes for each
377
HUC8 watershed was identified as the primary driver of interannual variability for that
378
watershed. Because some components of NANI were estimated by linear interpolation for
379
non-agricultural-census years (see Section 2.3), we repeated the analysis of temporal
380
variability using application only to the five agricultural census years, and found consistent
381
results (Table S3). Landuse was not considered as a driver of temporal variability because
382
little land use change occurred for the CONUS over the study period42,43 and this study
383
focused on interannual, rather than long-term drivers. The interannual variability in QTN in the
384
southeastern portion of Missouri, Lower Mississippi, Tennessee, and the western portion of
385
the South Atlantic-Gulf region is driven by PMAM,p>0.95, while PAnnual drives the year-to-year
386
variability in QTN for watersheds within most of the remaining HUC2 regions. NANI is the
387
primary driver of interannual variability only for watersheds with very low NANI values (e.g.,
388
northwestern portion of the Missouri region, as well as the Upper and Lower Colorado and
389
Rio Grande regions). The low NANI values in these watersheds lead to small changes in
390
NANI representing a large relative change in NANI that in turn drives year-to-year variability
391
in QTN. Han et al.17 previously identified the primary role of precipitation in explaining
392
interannual variability for a subset of the Lake Michigan watersheds examined in their study.
393
The role of precipitation23,24 or extreme precipitation25 has also been noted in a few field
394
studies comparing a limited number of years for small regions. Here we show that annual
395
precipitation or springtime extreme precipitation dominate the year-to-year variability in QTN
396
for HUC8 watersheds across the vast majority of the CONUS.
17
ACS Paragon Plus Environment
Environmental Science & Technology
Page 18 of 34
397
Focusing next on loading extremes, we find that years with extreme nitrogen loading
398
correspond to those with record precipitation. We define extreme loading simply as the
399
maximum loading observed in each HUC8 watershed over the 21-year examined period,
400
which corresponds to a loading greater than the 95th percentile over the examined period. The
401
occurrence of record loading events is associated with record PAnnual for 44% of watersheds
402
and record PMAM,p>0.95 for 48% of watersheds (77% of watersheds for either or both
403
precipitation variables), while being associated with record NANI for only 16% of watersheds
404
(Figure 3c and Figure S7). Note that these percentages do not sum to unity because multiple
405
extremes can occur concurrently, and also extreme loading in some HUC8 watersheds was
406
not associated with extremes in any of the predictor variables. Here again we repeated the
407
analysis using only the five agricultural census years to ensure that the lack of annual
408
information for some NANI components was not biasing results, and found consistent
409
conclusions (Table S3).
410
Overall, while the spatial variability in QTN for HUC8 watersheds is primarily driven
411
by nitrogen inputs, the interannual variability and the occurrence of extremes are instead
412
driven by hydrologic variability for the vast majority of the CONUS. This contrast is
413
especially evident in high loading areas within the Upper Mississippi, Ohio, Tennessee, and
414
Lower Mississippi regions, where NANI is the strongest predictor of long-term average
415
loading, but precipitation drives both interannual variability and extremes (Figure 2 and 3).
416
These findings are particularly relevant within the context of discussions regarding
417
management strategies aimed at reducing the occurrence of extreme water quality impacts
418
(e.g. harmful algal blooms, hypoxia), which are themselves linked to interannual variability in
419
nutrient loading, and especially to loading extremes.
420
nitrogen loading are clearly linked to corresponding trends in NANI14,44,53 and land use40, the
Although multi-decadal trends in
18
ACS Paragon Plus Environment
Page 19 of 34
Environmental Science & Technology
421
variability experienced from year-to-year is dominated by meteorological conditions. This
422
result has at least two implications. The first is that changes in land management (via NANI)
423
may take a long time to noticeably impact observed loading because the “signal” resulting
424
from changes to management is small relative to the “noise” of interannual variability due to
425
meteorology. The second is that strategies for alleviating water quality impacts must be based
426
on an assessment of the compounding roles of NANI and meteorological conditions on year-
427
to-year loading, rather than being developed based on understanding gleaned from analyses of
428
long-term-average hydrological conditions.
429
3.3
Aggregated estimates for large regions
430
Whereas the previous section focused on variability at the scale of HUC8 watersheds,
431
we turn now to estimates at aggregated scales. Annual TN load estimates at these aggregated
432
scales are obtained by summing annual TN loads for HUC8 watersheds falling within specific
433
regions, and the 95% prediction intervals are obtained using methodology described in the SI.
434
Although a comprehensive analysis of all the factors contributing to the observed patterns is
435
beyond the scope of this study, we provide some high-level observations that will need to be
436
explored further through additional studies.
437
We estimate that average annual nitrogen loading for the CONUS was 4.12±0.03 Tg-
438
N/yr for 1987-2007. This estimate and those for some large sub-regions are discussed within
439
the context of earlier estimates in the SI.
440
Perhaps surprisingly, we find substantial interannual variability even at this highly
441
aggregated scale, with annual load estimates ranging from 2.69±0.07 Tg-N/yr for 1988 to
442
5.34±0.17 Tg-N/yr for 1990 (Figure 4). The interannual variability in continental-scale
443
nitrogen loading has never been quantified previously, and the two-fold difference between
444
the smallest and largest estimated loading over a two-decade period demonstrates the value of 19
ACS Paragon Plus Environment
Environmental Science & Technology
Page 20 of 34
445
not relying on single-year snapshots or estimates based on hydrologically-average conditions
446
for understanding large-scale nitrogen dynamics. Indeed, the large variability in the annual
447
load for the CONUS is primarily driven by hydrological variability resulting from variations
448
in total and extreme precipitation rather than variability in NANI (Figure S8). The year 1988,
449
for example, had the lowest precipitation and very low extreme precipitation over the
450
examined period, while 1990 had very high total precipitation. Overall, average PAnnual and
451
average PMAM,p>0.95 over the CONUS can individually explain 76% and 31% of the
452
interannual variability in the estimated continental loading for 1987-2007, respectively, while
453
these two variables together explain 79%. Total NANI, on the other hand, explains only 9%
454
of the interannual variability in CONUS loading. Regionally, the interannual variability is
455
most closely attributable to the MARB, which accounts for over half (60%) of the CONUS
456
loading over the 21-year period (Figure 4), and for which the years with maximum and
457
minimum estimated loading coincide with those for the CONUS.
458
A similar story emerges when looking at four large basins that are associated with
459
major nutrient delivery points to the coastal ocean, namely the MARB, the Sacramento / San
460
Joaquin River Basin (SSRB), the Columbia River Basin (CRB), and the Chesapeake Bay
461
Basin (CBB) (Figure 4). The degree of interannual variability is substantial for all of these
462
basins, with ratios of the highest to lowest loading years of 2.3, 4.6, 4.4, and 2.7, respectively.
463
Among these basins, the two on the west coast of the US show substantially higher relative
464
interannual variability.
465
precipitation, with PAnnual explaining 76% (MARB) to 86% (CBB) of the basin-scale
466
interannual variability. Extreme precipitation PMAM,p>0.95 alone, on the other hand, explains
467
only a small fraction of the interannual variability for the CRB (14%), CBB (13%), and
468
MARB (16%), while explaining over half for the SSRB (57%).
The observed interannual variability is again dominated by
20
ACS Paragon Plus Environment
Page 21 of 34
Environmental Science & Technology
469
Together, the results for the CONUS and the four basins indicate that (1) interannual
470
variability is large and precipitation is a dominant driver thereof across the examined basins,
471
but (2) both the relative magnitude of interannual variability ((SSRB, CRB) > (CBB, MARB,
472
CONUS)) and the role of extreme precipitation (SSRB > CONUS > (MARB, CBB, CRB))
473
differ among regions. In addition, the differences in these two features neither track one
474
another across regions, nor the relative areas of the regions (CONUS > MARB >> CRB >
475
(CBB, SSRB)) or their average loading (CONUS > MARB >> (CBB, CRB) > SSRB). This
476
outcome further reinforces the need to analyze loading estimates in both a spatially and
477
temporally explicit manner in order to gain an understanding of current loading conditions,
478
and how these might be impacted by specific management strategies.
479
Finally we explore extremes in loading and in the driving variables through the lens of
480
four highest estimated loading years at the CONUS scale (1990, 1991, 1995, 1996) (Figure 4).
481
Note that the loading in these four years is not statistically significantly higher than in some
482
other years (Figure 4), and we are therefore simply using these years as illustrative case
483
studies. The high loading in these years results from large, spatially-coherent areas with
484
estimated loading substantially above the 1987-2007 average (Figure 5, warm colors in left
485
column). These years also exhibit large contiguous areas where the estimated loading is
486
either the highest or second highest over the 21-year period (i.e. above the 90th percentile of
487
observed loadings) (Figure 5, stippling in left column). The fact that these regions of extreme
488
loading tend to be spatially contiguous over large areas, rather than disjointed, means that
489
years with large overall loading can have disproportionate impacts on specific areas within
490
the CONUS. Indeed, 1990 and 1991 are also the two highest estimated loading years for the
491
MARB, 1995 is also the highest loading year for the SSRB, and 1996 is both the highest
492
loading year for the CRB and the second highest for the CBB (Figure 5, regions outlined in
493
left column). 21
ACS Paragon Plus Environment
Environmental Science & Technology
Page 22 of 34
494
Looking at these years and regions in turn, we observe that these extreme basin-scale
495
loading years correspond to years with large spatially-coherent areas within these basins
496
where annual precipitation or extreme springtime precipitation were above the 90th percentile
497
over the examined period (Figure 5, right column). Conversely, there is no such pattern in
498
NANI for those years (Figure 5, middle column). This finding further reinforces the
499
conclusion that loading extremes are driven by extremes in precipitation. While this was
500
observed in Section 3.2 at the HUC8 scale, we see here that this effect is compounded by the
501
spatially coherent nature of precipitation extremes. The topic of precipitation extremes, and
502
specifically their spatial signatures and temporal frequency, is one that is being actively
503
pursued in the climate science community54. The implication here being that if extremes in
504
loading at scales as large as the basins examined here are driven by spatially-contiguous
505
precipitation extremes, then changes to these loading extremes resulting from climate change
506
need to be carefully considered within the context of evolving management strategies.
507
4
Implications of watershed- and aggregated-scale analyses
508
This study provides the first spatially and temporal explicit estimates of TN loading
509
for watersheds throughout the CONUS, revealing considerable spatial and temporal
510
variability in annual TN loading. Regions of high relative interannual variability coincide
511
with regions of high overall loading. In addition, whereas NANI is the primary driver of the
512
spatial variability in TN loading, annual precipitation and extreme springtime precipitation
513
drive both the interannual variability in TN loading and the occurrence of loading extremes.
514
Taken together, these findings point to a fundamental challenge in managing regions
515
with high nutrient loading, because these regions also exhibit the strongest interannual
516
variability, and because the impact of changes in management practices will be modulated by
517
meteorological variability and climatic trends. The degree of interannual variability remains
22
ACS Paragon Plus Environment
Page 23 of 34
Environmental Science & Technology
518
remarkably high even at large scales such as the CONUS and the basins corresponding to
519
major nutrient delivery points to the coastal oceans examined here.
520
These results have implications for management strategies aimed at reducing the
521
occurrence of extreme water quality impacts (e.g. harmful algal blooms, hypoxia), which are
522
themselves linked to interannual variability in nutrient loading, and especially to loading
523
extremes. First, because observed year-to-year variability in nutrient loading is dominated by
524
meteorological conditions, changes in land management (via NANI) may take a long time to
525
noticeably impact observed loading, due to a low signal-to-noise ratio. Second, the findings
526
of this study also put a spotlight on the fact that strategies aimed at alleviating water quality
527
impacts must be based on an assessment of the compounding roles of NANI and
528
meteorological conditions on year-to-year loading, rather than being informed by analyses
529
relying on long-term-average hydrological conditions.
530
precipitation variables that best predict TN loading tend to be spatially coherent over large
531
regions, extremes in loading at the watershed scale also lead to extremes in loading at large
532
regional scales. This finding implies that precipitation extremes must be taken into account
533
when developing strategies for managing, anticipating, or preventing the most severe water
534
quality impacts55. Fourth, long-term management strategies must acknowledge and account
535
for the fact that meteorological conditions are themselves evolving as climate changes, which
536
will affect how nitrogen inputs into watersheds translate into nitrogen loading to waterways.
537
Changes to precipitation patterns resulting from climate change must therefore be carefully
538
considered within the context of evolving management strategies.
539
Supporting Information
540
(S1) Prediction intervals for the multiple linear regression model; (S2) Use of NADP data for
541
atmospheric N deposition; (S3) Lagged NANI, tile-drainage, and temperature; (S4)
Third, because extremes in the
23
ACS Paragon Plus Environment
Environmental Science & Technology
Page 24 of 34
542
Comparison of TN load estimates for large regions; (Table S1) Predictor variables evaluated;
543
(Table S2) NLCD land cover classes; (Table S3) Results of sensitivity tests; (Figure S1) Time
544
series of NANI components; (Figure S2) NANI data transformation; (Figure S3) Predicted vs.
545
observed TN flux; (Figure S4) Relationships between TN flux and predictor variables; (Figure
546
S5) Regional comparison of predicted vs. observed TN flux; (Figure S6) Time-averaged maps
547
of selected predictor variables; (Figure S7) Primary drivers of extreme loading; (Figure S8)
548
Time series of spatially-averaged predictor variables; (Figure S9) Comparison of NOX
549
deposition from NADP and CMAQ. The annual TN flux estimates are available from the
550
authors upon request.
551
Corresponding Author
552
*Phone: (857) 998-7784, e-mail: esinha@stanford.edu.
553
Notes
554
The authors declare no competing financial interest.
555
Acknowledgements
556
This material is based upon work supported by the National Science Foundation (NSF) under
557
Grants 1313897. The authors thank Chao Li and Jeff Ho for comments and feedback on early
558
manuscripts. We acknowledge the ongoing effort of the U.S. Geological Survey and NCWQR
559
in collecting discharge and water quality data for various streams across the CONUS.
560
References
561 562 563 564 565 566 567
(1)
Ciais, P.; Sabine, C.; Bala, G.; Bopp, L.; Brovkin, V.; Canadell, J.; Chhabra, A.; DeFries, R.; Galloway, J.; Heimann, M.; et al. Carbon and Other Biogeochemical Cycles. In Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change; Stocker, T. F., Qin, D., Plattner, G.-K., Tignor, M., Allen, S. K., Boschung, J., Nauels, A., Xia, Y., Bex, V., Midgley, P. M., Eds.; Cambridge University Press: Cambridge, United Kingdom and New York, NY, USA, 2013; pp 465–570. 24
ACS Paragon Plus Environment
Page 25 of 34
Environmental Science & Technology
568 569 570
(2)
Davidson, E. A.; David, M. B.; Galloway, J. N.; Goodale, C. L.; Haeuber, R.; Harrison, J. A.; Howarth, R. W.; B, J. D.; R, L. R.; T, N. B.; et al. Excess Nitrogen in the U.S. Environment: Trends, Risks, and Solutions. Issues in Ecology 2012, No. 15.
571 572
(3)
US EPA. Reactive Nitrogen in the United States: An Analysis of Inputs, Flows, Consequences, and Management Options; EPA-SAB-11-013; 2011; p 140.
573 574 575
(4)
Heisler, J.; Glibert, P. M.; Burkholder, J. M.; Anderson, D. M.; Cochlan, W.; Dennison, W. C.; Dortch, Q.; Gobler, C. J.; Heil, C. A.; Humphries, E.; et al. Eutrophication and harmful algal blooms: A scientific consensus. Harmful Algae 2008, 8 (1), 3–13.
576 577
(5)
Diaz, R. J.; Rosenberg, R. Spreading Dead Zones and Consequences for Marine Ecosystems. Science 2008, 321 (5891), 926–929.
578 579 580
(6)
Jewett, E. B.; Eldridge, P. M.; Burke, M. K.; Buxton, H. T.; Diaz, R. J. Scientific Assessment of Hypoxia in U. S. Coastal Waters; Committee on Environment and Natural Resources, 2010.
581 582 583
(7)
Lopez, C. B. Scientific assessment of marine harmful algal blooms; Council on Environmental Quality, Office of Science and Technology Policy, Executive Office of the President, 2008.
584 585
(8)
Lopez, C. B.; Jewett, E. B.; Dortch, Q.; Walton, B. T.; Hudnell, H. K. Scientific assessment of freshwater harmful algal blooms. 2008.
586 587
(9)
US EPA. National Lakes Assessment: a collaborative survey of the Nation’s lakes.; EPA 841-R-09-001; USEPA (US Environmental Protection Agency), 2009.
588 589
(10) US EPA. National Rivers and Streams Assessment (NRSA) 2008-2009 Draft Report (PDF); EPA/841/D-13/001; USEPA (US Environmental Protection Agency), 2013.
590 591 592
(11) de Vries, W.; Leip, A.; Reinds, G. J.; Kros, J.; Lesschen, J. P.; Bouwman, A. F. Comparison of land nitrogen budgets for European agriculture by various modeling approaches. Environmental Pollution 2011, 159 (11), 3254–3268.
593 594
(12) GAO, U. Water Quality: Key EPA and State Decisions Limited by Inconsistent and Incomplete Data; GAO/RCED-00-54; U.S. General Accounting Office, 2000.
595 596
(13) David, M. B.; Drinkwater, L. E.; McIsaac, G. F. Sources of nitrate yields in the Mississippi River Basin. Journal of Environmental Quality 2010, 39 (5), 1657–1667.
597 598 599
(14) Goolsby, D. A.; Battaglin, W. A. Long-term changes in concentrations and flux of nitrogen in the Mississippi River Basin, USA. Hydrol. Process. 2001, 15 (7), 1209– 1226.
600 601
(15) Smith, R. A.; Schwarz, G. E.; Alexander, R. B. Regional interpretation of water-quality monitoring data. Water Resour. Res. 1997, 33 (12), 2781–2798.
602 603 604
(16) Booth, M. S.; Campbell, C. Spring Nitrate Flux in the Mississippi River Basin: A Landscape Model with Conservation Applications. Environ. Sci. Technol. 2007, 41 (15), 5410–5418.
605 606 607
(17) Han, H.; Allan, J. D.; Scavia, D. Influence of climate and human activities on the relationship between watershed nitrogen input and river export. Environmental Science & Technology 2009, 43 (6), 1916–1922.
608 609
(18) Howarth, R. W.; Billen, G.; Swaney, D.; Townsend, A.; Jaworski, N.; Lajtha, K.; Downing, J. A.; Elmgren, R.; Caraco, N.; Jordan, T.; et al. Regional nitrogen budgets 25
ACS Paragon Plus Environment
Environmental Science & Technology
Page 26 of 34
610 611 612
and riverine N & P fluxes for the drainages to the North Atlantic Ocean: Natural and human influences. Nitrogen Cycling in the North Atlantic Ocean and its Watersheds 1996, 75–139.
613 614 615
(19) Howarth, R. W.; Swaney, D.; Boyer, E. W.; Marino, R.; Jaworski, N.; Goodale, C. The influence of climate on average nitrogen export from large watersheds in the Northeastern United States. Biogeochemistry 2006, 79 (1–2), 163–186.
616 617 618 619
(20) Howarth, R. W.; Swaney, D.; Billen, G.; Garnier, J.; Hong, B.; Humborg, C.; Johnes, P.; Mörth, C.-M.; Marino, R. Nitrogen fluxes from the landscape are controlled by net anthropogenic nitrogen inputs and by climate. Frontiers in Ecology and the Environment 2011, 10 (1), 37–43.
620 621 622
(21) Raymond, P. A.; David, M. B.; Saiers, J. E. The impact of fertilization and hydrology on nitrate fluxes from Mississippi watersheds. Current Opinion in Environmental Sustainability 2012, 4 (2), 212–218.
623 624 625
(22) Andersen, H. E.; Kronvang, B.; Larsen, S. E.; Hoffmann, C. C.; Jensen, T. S.; Rasmussen, E. K. Climate-change impacts on hydrology and nutrients in a Danish lowland river basin. Science of the Total Environment 2006, 365 (1), 223–237.
626 627 628
(23) Gentry, L. E.; David, M. B.; Below, F. E.; Royer, T. V.; McIsaac, G. F. Nitrogen Mass Balance of a Tile-drained Agricultural Watershed in East-Central Illinois. Journal of Environment Quality 2009, 38 (5), 1841.
629 630 631
(24) Vanni, M. J.; Renwick, W. H.; Headworth, J. L.; Auch, J. D.; Schaus, M. H. Dissolved and particulate nutrient flux from three adjacent agricultural watersheds: A five-year study. Biogeochemistry 2001, 54 (1), 85–114.
632 633 634 635
(25) Royer, T. V.; David, M. B.; Gentry, L. E. Timing of riverine export of nitrate and phosphorus from agricultural watersheds in Illinois: Implications for reducing nutrient loading to the Mississippi River. Environmental Science & Technology 2006, 40 (13), 4126–4131.
636 637 638 639
(26) Hirsch, R. M.; Moyer, D. L.; Archfield, S. A. Weighted Regressions on Time, Discharge, and Season (WRTDS), with an Application to Chesapeake Bay River Inputs1. JAWRA Journal of the American Water Resources Association 2010, 46 (5), 857–880.
640 641 642
(27) USGS. National Water Information System data available on the World Wide Web (USGS Water Data for the Nation) http://waterdata.usgs.gov/nwis/ (accessed Jun 1, 2014).
643 644 645 646
(28) Pellerin, B. A.; Bergamaschi, B. A.; Gilliom, R. J.; Crawford, C. G.; Saraceno, J.; Frederick, C. P.; Downing, B. D.; Murphy, J. C. Mississippi River Nitrate Loads from High Frequency Sensor Measurements and Regression-Based Load Estimation. Environ. Sci. Technol. 2014, 48 (21), 12612–12619.
647 648 649
(29) Falcone, J. A.; Carlisle, D. M.; Wolock, D. M.; Meador, M. R. GAGES: A stream gage database for evaluating natural and altered flow conditions in the conterminous United States: Ecological Archives E091-045. Ecology 2010, 91 (2), 621–621.
650 651
(30) Hong, B.; Swaney, D. P.; Howarth, R. W. A toolbox for calculating net anthropogenic nitrogen inputs (NANI). Environmental Modelling & Software 2011, 26 (5), 623–633.
26
ACS Paragon Plus Environment
Page 27 of 34
Environmental Science & Technology
652 653 654
(31) Gronberg, J. M.; Spahr, N. E. County-level Estimates of Nitrogen and Phosphorus from Commercial Fertilizer for the Conterminous United States, 1987-2006; US Department of the Interior, US Geological Survey, 2012.
655 656
(32) National Atmospheric Deposition Program. NADP Program Office, Illinois State Water Survey: 2204 Griffith Dr., Champaign, IL 61820, 2007.
657 658 659
(33) Ruddy, B. C.; Lorenz, D. L.; Mueller, D. K. County-level estimates of nutrient inputs to the land surface of the conterminous United States, 1982-2001; US Department of the Interior, US Geological Survey, 2006.
660
(34) PRISM, C. G. PRISM Climate Group; Oregon State University, 2014.
661 662
(35) Sugg, Z. Assessing US farm drainage: Can GIS lead to better estimates of subsurface drainage extent. World Resources Institute, Washington, DC 2007, 20002.
663
(36) Schwarz, G. Estimating the Dimension of a Model. Ann. Statist. 1978, 6 (2), 461–464.
664 665 666 667
(37) Basu, N. B.; Destouni, G.; Jawitz, J. W.; Thompson, S. E.; Loukinova, N. V.; Darracq, A.; Zanardo, S.; Yaeger, M.; Sivapalan, M.; Rinaldo, A. Nutrient loads exported from managed catchments reveal emergent biogeochemical stationarity. Geophysical Research Letters 2010, 37 (23).
668 669 670
(38) McIsaac, G. F.; David, M. B.; Gertner, G. Z. Illinois River Nitrate-Nitrogen Concentrations and Loads: Long-term Variation and Association with Watershed Nitrogen Inputs. Journal of Environmental Quality 2016.
671 672 673 674
(39) Van Breemen, N.; Boyer, E. W.; Goodale, C. L.; Jaworski, N. A.; Paustian, K.; Seitzinger, S. P.; Lajtha, K.; Mayer, B.; Van Dam, D.; Howarth, R. W. Where did all the nitrogen go? Fate of nitrogen inputs to large watersheds in the northeastern USA. Biogeochemistry 2002, 57 (1), 267–293.
675 676
(40) Turner, R. E.; Rabalais, N. N. Linking landscape and water quality in the Mississippi River Basin for 200 years. Bioscience 2003, 53 (6), 563–572.
677 678 679 680
(41) Fry, J. A.; Xian, G.; Jin, S.; Dewitz, J. A.; Homer, C. G.; LIMIN, Y.; Barnes, C. A.; Herold, N. D.; Wickham, J. D. Completion of the 2006 national land cover database for the conterminous United States. Photogrammetric Engineering and Remote Sensing 2011, 77 (9), 858–864.
681 682 683 684
(42) Fry, J. A.; Coan, M. J.; Homer, C. G.; Meyer, D. K.; Wickham, J. D. Completion of the National Land Cover Database (NLCD) 1992-2001 Land Cover Change Retrofit Product; Open-File Report; USGS Numbered Series 2008–1379; U.S. Geological Survey, 2009.
685 686 687 688
(43) Homer, C. G.; Dewitz, J. A.; Yang, L.; Jin, S.; Danielson, P.; Xian, G.; Coulston, J.; Herold, N. D.; Wickham, J. D.; Megown, K. Completion of the 2011 National Land Cover Database for the conterminous United States-Representing a decade of land cover change information. Photogramm. Eng. Remote Sens 2015, 81 (5), 345–354.
689 690
(44) Donner, S. D.; Kucharik, C. J.; Foley, J. A. Impact of changing land use practices on nitrate export by the Mississippi River. Global Biogeochemical Cycles 2004, 18 (1).
691 692 693
(45) Hong, B.; Swaney, D. P.; Howarth, R. W. Estimating net anthropogenic nitrogen inputs to US watersheds: Comparison of methodologies. Environmental Science & Technology 2013, 47 (10), 5199–5207.
27
ACS Paragon Plus Environment
Environmental Science & Technology
Page 28 of 34
694 695 696 697
(46) Michalak, A. M.; Anderson, E. J.; Beletsky, D.; Boland, S.; Bosch, N. S.; Bridgeman, T. B.; Chaffin, J. D.; Cho, K.; Confesor, R.; Daloğlu, I.; et al. Record-setting algal bloom in Lake Erie caused by agricultural and meteorological trends consistent with expected future conditions. PNAS 2013, 110 (16), 6448–6452.
698 699 700
(47) Donner, S. D.; Scavia, D. How climate controls the flux of nitrogen by the Mississippi River and the development of hypoxia in the Gulf of Mexico. Limnology and Oceanography 2007, 52 (2), 856–861.
701 702 703
(48) Goodale, C. L.; Lajtha, K.; Nadelhoffer, K. J.; Boyer, E. W.; Jaworski, N. A. Forest nitrogen sinks in large eastern U.S. watersheds: estimates from forest inventory and an ecosystem model. Biogeochemistry 2002, 239–266.
704 705
(49) Schlesinger, W. H. On the fate of anthropogenic nitrogen. PNAS 2009, 106 (1), 203– 208.
706 707 708 709 710
(50) Mitsch, W. J.; Day, J. W.; Gilliam, J. W.; Groffman, P. M.; Hey, D. L.; Randall, G. W.; Wang, N. Reducing Nitrogen Loading to the Gulf of Mexico from the Mississippi River Basin: Strategies to Counter a Persistent Ecological Problem Ecotechnology—the use of natural ecosystems to solve environmental problems—should be a part of efforts to shrink the zone of hypoxia in the Gulf of Mexico. BioScience 2001, 51 (5), 373–388.
711 712 713
(51) Alexander, R. B.; Smith, R. A.; Schwarz, G. E.; Boyer, E. W.; Nolan, J. V.; Brakebill, J. W. Differences in Phosphorus and Nitrogen Delivery to The Gulf of Mexico from the Mississippi River Basin. Environ. Sci. Technol. 2008, 42 (3), 822–830.
714 715 716
(52) Robertson, D. M.; Saad, D. A.; Schwarz, G. E. Spatial Variability in Nutrient Transport by HUC8, State, and Subbasin Based on Mississippi/Atchafalaya River Basin SPARROW Models. J Am Water Resour Assoc 2014, 50 (4), 988–1009.
717 718
(53) McIsaac, G. F.; David, M. B.; Gertner, G. Z.; Goolsby, D. A. Eutrophication: Nitrate flux in the Mississippi river. Nature 2001, 414 (6860), 166–167.
719 720 721 722
(54) Field, C. B.; Barros, V. R.; Stocker, T. F.; Qin, D.; Dokken, D. J.; Ebi, K. L.; Mastrandrea, M. D.; Mach, K. J.; Plattner, G.-K.; Allen, S. K.; et al. Managing the risks of extreme events and disasters to advance climate change adaptation: special report of the intergovernmental panel on climate change; Cambridge University Press, 2012.
723 724 725
(55) Michalak, A. M. Study role of climate change in extreme threats to water quality. Nature 2016, 535 (7612), 349–350.
28
ACS Paragon Plus Environment
Page 29 of 34
Environmental Science & Technology
Figure 1: Location of USGS gages and their catchments used in model build. Locations of USGS gages are shown with black dots, and associated catchments are shown in red outline. The black polygons show the first level of hydrologic units (HUC2) for the continental US (CONUS). The light grey polygons show the eight digit hydrologic units (HUC8) within the CONUS on which the statistical model was applied. The shaded regions represent major nutrient delivery points to the coastal ocean considered in this work.
ACS Paragon Plus Environment
Environmental Science & Technology
a: Average
b: Standard deviation
c: Coefficient of variation Figure 2: (a) 1987-2007 mean estimated annual TN flux (QTN), (b) the standard deviation of the interannual variability in annual flux, and (c) the coefficient of variation (i.e. the ratio of the standard deviation to the mean) for HUC8 watersheds.
ACS Paragon Plus Environment
Page 30 of 34
Page 31 of 34
Environmental Science & Technology
ln(kg - N/(km2 yr)) 0-1 1-2 2-3 3-4 4-5
fNANI PAnnual LUF, SH LUW
a: Drivers of spatial variability
PAnnual PMAM,p>0.95 fNANI
b: Drivers of interannual variability
PAnnual and/or PMAM,p>0.95 (PAnnual and/or PMAM,p>0.95) and fNANI
c: Drivers of exteme loading
fNANI None
Figure 3: Primary drivers of (a) spatial variability (b) interannual variability, and (c) extreme loading for HUC8 watersheds. The primary drivers were identified as described in Section 3.2. The drivers of spatial variability are color-coded both by the driver, and by the mean magnitude of its contribution. The primary drivers of interannual variability and loading extremes are colorcoded by driver. All variables are as defined in SI Table 1.
ACS Paragon Plus Environment
Environmental Science & Technology
Page 32 of 34
Continental US
6000 4000 2000 0
Mississippi Atchafalaya River Basin
4000 3000 2000
Total load [Gg − N yr]
1000 0 250 200 150 100 50 0
Sacramento / San Joaquin River Basin
Columbia River Basin
500 400 300 200 100 0
Chesapeake Bay Basin
300 200 100
Annual
2007
2005
2003
2001
1999
1997
1995
1993
1991
1989
1987
0
Long−term mean
Figure 4: Annual and 1987-2007 mean load estimates for CONUS and watersheds associated with major nutrient delivery points to the coastal oceans.
ACS Paragon Plus Environment
Page 33 of 34
Environmental Science & Technology
Mississippi Atchafalaya River Basin
1990 Mississippi Atchafalaya River Basin
1991 Sacramento/San Joaquin River Basin
1995 Columbia River Basin
1996 Chesapeake Bay Basin
fNANI
PAnnual PMAM,p>0.95 PAnnual & PMAM,p>0.95
Figure 5: Difference between TN loading (QTN) for HUC8 watersheds in each of the four years with the highest estimated total CONUS loading and the 1987-2007 mean loading for each HUC8 watershed are presented in the left-hand column. Stippling shows HUC8 watersheds that experienced the highest or second highest estimated loading in a particular year (i.e. above 90th percentile). For each year, major nutrient delivery points to the coastal ocean that experience their highest or second highest loading are highlighted in the middle and right-hand columns. HUC8 watersheds within these regions are colored based on whether fNANI had the highest or second highest rank over the 21-year time period for the selected year and in the right side column based on whether PAnnual and/or PMAM,p>0.95 had the highest or second highest rank over the 21-year time period for the selected year.
ACS Paragon Plus Environment
Environmental Science & Technology
TOC/Abstract art Difference in total nitrogen loading from long-term mean
Mg - N/(km2 yr)
1988
(5.0, 7.1] (0.8, 5.0] (0.4, 0.8] (0.2, 0.4] (0.1, 0.2] (0, 0.1] (-0.1, 0] (-0.2, -0.1] (-0.4, -0.2] (-0.8, -0.4] [-3.1, -0.8]
Above average
Below average
1990
ACS Paragon Plus Environment
Page 34 of 34