Precipitation Dominates Interannual Variability of Riverine Nitrogen

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 ([email protected])

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


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


Page 3 of 34


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



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


Page 5 of 34


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



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


Page 7 of 34

142


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



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


Page 9 of 34


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



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


Page 11 of 34


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



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


Page 13 of 34


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



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


Page 15 of 34

324


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



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


Page 17 of 34


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



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


Page 19 of 34


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



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


Page 21 of 34


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



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


Page 23 of 34


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



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: [email protected].

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


Page 25 of 34


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



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


Page 27 of 34


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



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


Page 29 of 34


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.



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.


Page 30 of 34

Page 31 of 34


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.



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.


Page 33 of 34


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.



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


Page 34 of 34

Precipitation Dominates Interannual Variability of Riverine Nitrogen

Recommend Documents