MODELING NUTRIENT LOADING TO WATERSHEDS IN THE GREAT LAKES BASIN:
A DETAILED SOURCE MODEL AT THE REGIONAL SCALE
By
Emily Catherine Luscz

A THESIS
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Environmental Geosciences – Master of Science
2013

ABSTRACT
MODELING NUTRIENT LOADING TO WATERSHEDS IN THE GREAT LAKES
BASIN: A DETAILED SOURCE MODEL AT THE REGIONAL SCALE
by
Emily Catherine Luscz
Nutrient loading has been linked to many issues including eutrophication, harmful algal
blooms, and decreases in aquatic species diversity. In the Great Lakes, algal blooms continue
to plague Lake Erie and Saginaw Bay despite reduction in point source loading to the lakes.
Though many watershed nutrient models exist in the literature, there is generally a tradeoff
between the scale of the model and the level of detail regarding the individual sources of
nutrients and basin characteristics. To examine the link between watershed nutrient sources,
landscape processes, and in-stream loads in the Lower Peninsula of Michigan a spatially
explicit nutrient loading model was developed. The model is composed of two parts: GIS
models that predict nutrient sources and a statistical model that predicts nutrient loads
from the source models and basin characteristics. Observations collected during baseflow,
melt, and rain conditions from 2010-2012 were used to calibrate and validate the model. The
model predicts nutrient loads and provides information on the sources of nutrients within
each watershed and the relative contribution of different sources to the overall nutrient load.
The model results indicated that there is a high degree of variability in nutrient export
rates, even within the same land use class, and export rates can be as much as 15 times
greater for nitrogen and as much as 9 times greater for phosphorus depending on the season.
The model performed as well as other regional scale nutrient models and showed less bias
due to land use. In addition, nitrate isotope data agreed with the model’s predictions of the
relative source of nutrients. This work has the potential to provide valuable information to
environmental managers regarding how and where to target efforts to reduce nutrient loads
in surface water.

ACKNOWLEDGEMENTS

I would first like to thank my advisor, David Hyndman for his encouragement and support
throughout this process and particularly when I decided to finish my thesis while starting
a new job, half way across the country. I’m hugely indebted to Anthony Kendall for his
guidance, collaboration, and inspiration. From Python to Euchre, Anthony patiently taught
me so many things required to complete this work. Most importantly, he regularly challenged
me to be my best.
This work would not have been possible without the help of Briana Jasinski, Blaze Budd,
and Ryan Nagelkirk who worked 12+ hour days, collectively logging over a month in the
field with me. Their smiles and positive attitudes were invaluable, especially when field days
started before sunrise and ended after sunset. Anthony Kendall, Sherry Martin, Lon Cooper,
Bobby Chrisman, Erin Haacker, Brian Eustice, and Ryan Vannier also logged many hours
in the field for which I’m extremely grateful.
I would also like to acknowledge my friends and officemates, Brian Eustice, Sherry Martin,
Erin Haacker, Yuteng Ma, and Alex Kuhl for their moral support and for putting up with
all my whining. I hope that I’m lucky enough to have coworkers who are half as wonderful
as they are in my future.

iii

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

vi

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

CHAPTER 1 HIGH-RESOLUTION SPATIALLY EXPLICIT NUTRIENT SOURCE MODELS FOR THE LOWER PENINSULA OF MICHIGAN . . . . . . . . . . . . . . . . . . . . . .
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.1 Model Descriptions . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.2 Description of Study Area . . . . . . . . . . . . . . . . . . . . . . . .
1.2.3 GIS Source Models . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.3.1 Point Sources . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.3.2 Atmospheric Loading . . . . . . . . . . . . . . . . . . . . . .
1.2.3.3 Septic Tank Loading . . . . . . . . . . . . . . . . . . . . . .
1.2.3.4 Non-agricultural Fertilizer . . . . . . . . . . . . . . . . . . .
1.2.3.5 Animal Manure . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.3.6 Commercial Fertilizer . . . . . . . . . . . . . . . . . . . . .
1.2.4 Sample Collection . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.1 Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.2 Rates by Land use . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.3 Loading Rates and In-stream Concentration . . . . . . . . . . . . . .
1.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1
1
3
3
3
3
3
4
7
10
14
18
21
24
24
24
26
28

CHAPTER 2 PREDICTING NUTRIENT LOADS AND SOURCES: A
SPATIALLY EXPLICIT STATISTICAL MODEL . . . . . . 38
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.1.1 Review of Existing Approaches . . . . . . . . . . . . . . . . . . . . . 39
2.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.2.1 Model Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.2.2 Landscape and Pathway Descriptions . . . . . . . . . . . . . . . . . . 45
2.2.2.1 Watersheds and Travel Distance . . . . . . . . . . . . . . . . 46
2.2.2.2 Recharge . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
2.2.2.3 Tile Drainage . . . . . . . . . . . . . . . . . . . . . . . . . . 47
2.2.2.4 Stream Velocity and Travel Time . . . . . . . . . . . . . . . 47
2.2.3 Calibration/Validation Data . . . . . . . . . . . . . . . . . . . . . . . 49
2.2.4 Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
2.3.1 Sample Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
2.3.2 Model Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
iv

2.4

2.3.2.1 Model Performance . . .
2.3.2.2 Sensitivity . . . . . . . .
2.3.3 Sources . . . . . . . . . . . . . .
2.3.4 Pathways and Processes . . . . .
2.3.4.1 Basin Processing . . . .
2.3.4.2 In-stream Processing . .
2.3.5 Land use Export Rates . . . . . .
2.3.6 Residuals . . . . . . . . . . . . .
2.3.6.1 Residuals and Land Use
2.3.6.2 Model Bias . . . . . . .
Discussion . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

53
56
60
62
62
66
66
68
69
69
70

CHAPTER 3 USING NITRATE ISOTOPE DATA TO VALIDATE A
REGIONAL NUTRIENT SOURCE MODEL . . . . . . . . . 78
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
3.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
3.2.1 Sample Watersheds . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
3.2.2 Sample Collection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
3.2.3 Statistical Nitrate Model . . . . . . . . . . . . . . . . . . . . . . . . . 82
3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
3.3.1 Nitrate Sample Results . . . . . . . . . . . . . . . . . . . . . . . . . . 82
3.3.2 Model Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
3.3.3 Isotope Sample Results . . . . . . . . . . . . . . . . . . . . . . . . . . 84
3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
3.4.1 Soil Source . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
3.4.2 Denitrification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
3.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
Appendix A Total Nitrogen, Total Phosphorus, and Discharge Results
Appendix B Nitrate Results
Appendix C Isotope Results

v

LIST OF TABLES

Table 1.1

Nutrient loading rates by land use. . . . . . . . . . . . . . . . . . . . . . .

26

Table 1.2

Results of linear regression models of source loading and observed instream concentration. The α value represents the fitted coefficient. Bold
values indicated parameters with a p-value ≤0.05 . . . . . . . . . . . . . .

27

Table 1.3

Results of linear regression models of source loading within 1 km of surface water and observed in-stream concentration. The α value represents
the fitted coefficient. Bold values indicated parameters with a p-value ≤0.05 28

Table 2.1

Summary of estimated model parameters. R2 values are for estimated
and observed loads. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

55

Table 2.2

Summary of coefficients associated with each basin parameter. . . . . . .

55

Table 2.3

Normalized sensitivity of model to each basin parameter coefficient. Sensitivity is calculated as the percent change in the model prediction per
percent change in the parameter. . . . . . . . . . . . . . . . . . . . . . .

60

Summary of calculated basin parameters. Reduction factor units (Bsurf,
Btile, and Bgrd ) are the percent of the nutrient load that remains after
the first .1 kilometer (1 km for Bsurf ) of travel. Units of Rmin and Rmax
are the percent of the load that remains after the first day of in-stream
travel for the sub-watershed with the highest relative basin yield (Rmax )
and the lowest relative basin yield (Rmin). The bold values indicate the
coefficients with relatively high model sensitivity. . . . . . . . . . . . . . .

61

Table 2.4

Table 2.5

Percent of total load delivered to surface water via the groundwater pathway. 63

Table 2.6

Estimated nutrient exports from each land use in kg/ha/year. . . . . . . .

68

Table 2.7

Model errors calculated as percent of observed loads. Sample Diff. is the
mean percent difference between samples and duplicates. . . . . . . . . .

74

Results of the residual regressions. Slope indicates the percent change in
the model error per unit change in the parameter. Bold values indicate
significant parameters at a 95% confidence level. . . . . . . . . . . . . . .

74

Table 3.1

Watershed characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . .

80

Table 3.2

Total NOx loads for each sampled watershed.

83

Table 2.8

vi

. . . . . . . . . . . . . . .

Table 3.3

Summary of model fit parameters.

. . . . . . . . . . . . . . . . . . . . .

83

Table 3.4

Summary of coefficients associated with each basin parameter. . . . . . .

91

Table 3.5

Summary of calculated basin parameters. Reduction factor units (Bsurf,
Btile, and Bgrd ) are the percent of the nutrient load that remains after
the first kilometer of travel. Units of Rmin and Rmax are the percent of
the load that remains after the first day of in-stream travel for the subwatershed with the highest relative basin yield (Rmax ) and the lowest
relative basin yield (Rmin). The bold values indicate the coefficients
with relatively high model sensitivity. . . . . . . . . . . . . . . . . . . . .

92

End member isotope values predicted by mixing model and compared to
Kendall et al. (2007). . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

92

Table A.1 2010 Baseflow Results. ND = Non Detect. DUP = Duplicate Sample. . .

98

Table A.2 2010 Baseflow Results (continued). ND = Non Detect. DUP = Duplicate
Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

99

Table 3.6

Table A.3 2010 Baseflow Results (continued). ND = Non Detect. DUP = Duplicate
Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
Table A.4 2011 Melt Results (continued). ND = Non Detect. DUP = Duplicate
Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
Table A.5 2011 Melt Results (continued). ND = Non Detect. DUP = Duplicate
Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
Table A.6 2011 Melt Results (continued). ND = Non Detect. DUP = Duplicate
Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
Table A.7 2011 Rain Results. ND = Non Detect. DUP = Duplicate Sample . . . . . 104
Table A.8 2011 Rain Results (continued). ND = Non Detect. DUP = Duplicate
Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
Table A.9 2011 Rain Results (continued). ND = Non Detect. DUP = Duplicate
Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
Table A.10 2011 Baseflow Results. ND = Non Detect. DUP = Duplicate Sample . . . 107
Table A.11 2011 Baseflow Results (continued). ND = Non Detect. DUP = Duplicate
Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
Table A.12 2011 Baseflow Results (continued). ND = Non Detect. DUP = Duplicate
Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

vii

Table A.13 2012 Melt Results. ND = Non Detect. DUP = Duplicate Sample . . . . . 110
Table A.14 2012 Melt Results (continued). ND = Non Detect. DUP = Duplicate
Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
Table A.15 2012 Melt Results (continued). ND = Non Detect. DUP = Duplicate
Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
Table B.1 2011 Melt Results. ND = Non Detect. DUP = Duplicate Sample . . . . . 114
Table B.2 2011 Melt Results (continued). ND = Non Detect. DUP = Duplicate
Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
Table B.3 2011 Melt Results (continued). ND = Non Detect. DUP = Duplicate
Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
Table B.4 2012 Melt Results. ND = Non Detect. DUP = Duplicate Sample . . . . . 117
Table B.5 2012 Melt Results (continued). ND = Non Detect. DUP = Duplicate
Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
Table B.6 2012 Melt Results (continued). ND = Non Detect. DUP = Duplicate
Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
Table C.1 2012 Melt Results. ND = Non Detect. DUP = Duplicate Sample . . . . . 121
Table C.2 2012 Melt Results (continued). ND = Non Detect. DUP = Duplicate
Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
Table C.3 2012 Melt Results (continued). ND = Non Detect. DUP = Duplicate
Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

viii

LIST OF FIGURES

Figure 1.1

Flow chart summarizing process for estimating point sources. . . . . . . .

5

Figure 1.2

Location of atmospheric deposition monitoring networks. The Sibley
and Isle Royale sites were compared to determine the contribution of
organic nitrogen to the total wet deposition load. For interpretation of
the references to color in this and all other figures, the reader is referred
to the electronic version of this thesis. . . . . . . . . . . . . . . . . . . .

7

Figure 1.3

Flow chart summarizing process for estimating atmospheric deposition. .

8

Figure 1.4

Comparison between the population served by treatment plants based
on EPA’s CWNS survey and the population served based on the estimated treatment plant service areas and the 2010 Census. . . . . . . . .

9

Figure 1.5

Flow chart summarizing process used to estimate loading from septic tanks. 11

Figure 1.6

Flow chart summarizing process used to estimate non-agricultural fertilizer. 15

Figure 1.7

Linear regressions of animal weight and nutrient excretion rates. These
linear models were used to predict animal excretion rates for animals not
discussed in the NRCS Waste Management Handbook (NRCS , 2008)
including bison, elk, emus, pheasant and others. . . . . . . . . . . . . . .

16

Figure 1.8

Applied manure predicted by the model for Huron County, MI. . . . . .

18

Figure 1.9

Manure Application Part I: Flow chart summarizing the process used
to estimate the location and inventories of farms that produce manure. .

19

Figure 1.10 Manure Application Part II: Flow chart summarizing the process used
to estimate the rate and extent of manure application. . . . . . . . . . .

20

Figure 1.11 Flow chart summarizing processes used to estimate loading from commercial agricultural fertilizer. . . . . . . . . . . . . . . . . . . . . . . . .

22

Figure 1.12 Sampled watersheds with dominant land use indicated. . . . . . . . . . .

23

Figure 1.13 Total nitrogen application rates for each modeled source. . . . . . . . . .

25

Figure 1.14 Total phosphorus application rates for each modeled source. . . . . . . .

31

Figure 1.15 Proportion of the total nutrient load contributed by each source for the
entire Lower Peninsula and within different land uses. . . . . . . . . . . .

32

ix

Figure 1.16 Total nitrogen and phosphorus loading in the Lower Peninsula. The
ratio of total phosphorus to total nitrogen is shown in the bottom panel.

33

Figure 1.17 Box plots of total nitrogen loading rate variability in sampled watersheds. Loading rates are subset by land use. . . . . . . . . . . . . . . . .

34

Figure 1.18 Box plots of total phosphorus loading rate variability in sampled watersheds. Loading rates are subset by land use. . . . . . . . . . . . . . . . .

35

Figure 1.19 Graphs showing linear regression model predictions compared to observed concentrations. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

Figure 1.20 Graphs showing linear regression model predictions (based on nutrient
loads within 1 km of surface water) compared to observed concentrations.

37

Figure 2.1

Qualitative comparison of nutrient modeling approaches. . . . . . . . . .

42

Figure 2.2

Conceptual model illustrating the sources and pathways described in
the model. Nutrient sources include septic tanks (Qseptic), agricultural
chemical fertilizer (QchemAg), manure (Qmanure), atmospheric deposition (Qatm), non-agricultural fertilizer (QnonAg), and point sources
(Qpoint). R, Bgrd, Btile, Bsep, and Bsurf are the reduction factors and
describe loss along the in-stream, groundwater, tile, septic, and surface
runoff pathways, respectively. . . . . . . . . . . . . . . . . . . . . . . . .

45

A comparison of a sub-basin based approach (such as the USGS SPARROW model) and the spatially explicit approach described here. . . . .

46

Estimated tile drained areas for the Lower Peninsula of Michigan (indicated in black). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

48

Graphs showing the relationship between discharge and channel area.
Regressions were fit to USGS channel area and discharge data for observations between the 75th and 90th percentile for discharge (high flow)
and observations less than the 20th percentile (low flow). Gauge data
was subset by geology type (lacustrine, outwash, and till). . . . . . . . .

50

In-stream travel time to sample point for the Saginaw Bay watershed
during snow melt. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51

Location of Boardman-Charlevoix, Muskegon River, Saginaw Bay, and
Grand River watersheds. . . . . . . . . . . . . . . . . . . . . . . . . . . .

52

Sampled watersheds with concentration (µg/L) indicated for each sampling campaign. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

Figure 2.3

Figure 2.4

Figure 2.5

Figure 2.6

Figure 2.7

Figure 2.8

x

Figure 2.9

Graphs of observed and modeled TN loads and area normalized loads
for the calibration dataset. . . . . . . . . . . . . . . . . . . . . . . . . . .

57

Figure 2.10 Graphs of observed and modeled TP loads and area normalized loads
for the calibration dataset. . . . . . . . . . . . . . . . . . . . . . . . . . .

58

Figure 2.11 Graphs of observed and modeled loads for the validation dataset. . . . .

59

Figure 2.12 Source contributions to total loading to streams. . . . . . . . . . . . . . .

61

Figure 2.13 Percent of observed TN melt load derived from each source. . . . . . . .

63

Figure 2.14 Percent of observed TP melt load derived from each source. . . . . . . .

64

Figure 2.15 Percent of total applied load delivered to surface water. . . . . . . . . . .

65

Figure 2.16 Source contributions to total stream loading within different land use
classes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

67

Figure 2.17 Box plots of model residuals. . . . . . . . . . . . . . . . . . . . . . . . . .

73

Figure 2.18 Model residuals for sampled watersheds (Equation 2.9). Only incremental watershed area is shown. . . . . . . . . . . . . . . . . . . . . . . . . .

73

Figure 2.19 Baseflow model residuals and land use regression results.

. . . . . . . .

75

Figure 2.20 Melt model residuals and land use regression results. . . . . . . . . . . .

76

Figure 2.21 Rain model residuals and land use regression results. . . . . . . . . . . .

77

Figure 3.1

Sample locations and incremental watersheds. . . . . . . . . . . . . . . .

90

Figure 3.2

NOx concentrations for sampled sub-watersheds. Incremental watershed areas are shown. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

90

Graphs of observed and modeled nitrate loads and area normalized loads
for the calibration and validation data sets. . . . . . . . . . . . . . . . .

91

Figure 3.4

Model predicted source of exported nitrate in watersheds. . . . . . . . . .

92

Figure 3.5

Measured isotopes by watershed . . . . . . . . . . . . . . . . . . . . . . .

93

Figure 3.6

Isotope data with dominant source indicated (isotopic classification from
Kendall et al. (2007)). . . . . . . . . . . . . . . . . . . . . . . . . . . . .

93

Predicted isotope values plotted with values reported by Kendall et al.
(2007) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

94

Figure 3.3

Figure 3.7

xi

Figure 3.8

Nitrate vs. δ 15 NN O3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xii

95

CHAPTER 1
HIGH-RESOLUTION SPATIALLY EXPLICIT NUTRIENT SOURCE
MODELS FOR THE LOWER PENINSULA OF MICHIGAN

1.1

Introduction

Elevated concentrations of nitrogen and phosphorus in lakes, rivers, groundwater and streams
have become a major concern for environmental managers. Nutrients have been identified
as one of the leading causes of river and lake impairment in the United States and have
been an ongoing issue in the Great Lakes, specifically Lake Erie (USEPA, 2009, 2002b).
Blooms of the cyanobacteria Microsystis can potentially produce the toxin Microsytin which
is hazardous to human and aquatic animal health (Davidson et al., 2012; Correll , 1998).
Growth of clydophera, an algae scum that grows on solid substrates, is not only odorous and
unsightly but can clog water intake pipes. Excessive growth of algae and cyanobacteria lead
to oxygen depletion which impacts benthic organisms and can cause fish kills (Auer et al.,
2010; Anderson et al., 2002). Issues associated with nutrient loading to the Great Lakes
include eutrophication, harmful algal blooms, hypoxia, and decreases in biologic diversity
making water unsuitable for recreational, industrial, and municipal uses. Seasonal algal
blooms are common in the western basin of Lake Erie and Saginaw Bay (Hinderer and
Murray, 2011; Dolan and Chapra, 2012).
Since the late 1960’s cooperative efforts by the governments of the United States and
Canada resulted in extensive monitoring, reporting, and legislation in the Great Lakes Basin
with the goal of improving water quality in the Great Lakes. Legislation was passed to
limit phosphorus point source effluents and establish limits on phosphorus loads to the lakes
(Nicholls et al., 2001; Dolan et al., 1981; Dolan, 1993). Initial improvements were realized
through point source load reduction, but since 1991 total phosphorus loading to Lake Erie has
been more variable and even shows signs of increasing. This has been attributed to increases

1

in non-point sources of phosphorus (Dolan and McGunagle, 2005; Moon and Carrick , 2007).
Increases in nitrate delivery to the Great Lakes, particularly from agricultural watersheds,
have also been observed (Smith et al., 1987).
Addressing the impact of non-point sources on in-stream nutrient loads requires an understanding of the location and rates of nutrient application. However, non-point sources of
nitrogen and phosphorus are difficult to quantify, because they can not be directly measured,
and how and where they are applied affect the delivery of nutrients to surface water (Carpenter et al., 1998; Nikolaidis et al., 1998). Data reporting related to sources of nutrients
varies from state to state; often data is only available on a county or state level or may not
be directly available.
One approach used to predict nutrient sources across large spatial scales is to use other
variables such as land use and population to estimate non-point sources of nutrients. This
type of approach (ie., Johnes (1996)) assumes that changes in land use or land management
will lead to proportional responses in nutrient sources which may lead to problems when there
is a large degree of variability in sources within a particular land use. Reliable estimates
of non-point sources are necessary to develop strategies to manage non-point sources and
predict how they will change with future climate and land use.
In this study, GIS models were developed to generate a spatially explicit description of
nutrient sources. Estimates of nutrient sources are compared to sampled surface water nitrogen and phosphorus concentrations in order to better understand the retention of nutrients
within watersheds and the contribution of nutrient sources to surface water loading across
various land uses. The resulting descriptions provide a basis to quantify the impacts of
non-point sources and serve as inputs for watershed nutrient loading models.

2

1.2

Methods

1.2.1

Model Descriptions

GIS models were developed to estimate nitrogen and phosphorus inputs from atmospheric
deposition (Atm), chemical agricultural fertilizer application (ChemAg), non-agricultural
fertilizer application (NonAg), manure production/application (Manure), septic tanks (Septic), and point sources (Point). Nutrient loading rates for each source were estimated for
30 meter cells across the Lower Peninsula. The method used to estimate each source is
described in the following sections. Each section is accompanied by a flow chart illustrating
the steps taken to generate estimates of nutrient sources.

1.2.2

Description of Study Area

Source models were constructed for the Lower Peninsula of Michigan which has watersheds
that drain to 3 of the 5 Great Lakes: Lakes Michigan, Huron and Erie. Watersheds in the
southern portion of the study area are dominated by agriculture, specifically row crops, and
include significant urban centers including Grand Rapids and Detroit. The primary land
cover in the northern portion of the study area is forest with small pockets of agriculture.

1.2.3
1.2.3.1

GIS Source Models
Point Sources

Relative to non-point sources, point sources can be easily identified, quantified, and monitored, and were thus early targets of management efforts to improve the nation’s water
quality (Dolan, 1993). The National Pollutant Discharge Elimination System (NPDES)
monitors and regulates effluent information for facilities that discharge to surface water in
the United States. Flows and concentrations of various parameters are measured and reported according to a facility’s permit. The parameters sampled and the reporting frequency

3

vary between facilities, making assessment of an average load challenging (Maupin and Ivahnenko, 2011; McMahon et al., 2007). Due to issues with data entry, data from NPDES
permits (as reported to the Permit Compliance System-PCS and Integrated Compliance Information System-ICIS for the National Pollutant Discharge Elimination System-NPDES)
required review and considerable quality control before it could be used. There are also issues
with the completeness of the data reported to the PCS and ICIS-NPDES databases (Maupin
and Ivahnenko, 2011). State agencies are required to report data from major dischargers
but are not required to report data from minor facilities. In Michigan around 50-60% of
facilities have reported data for years 2007-2011. For major facilities, the percent is higher,
almost 100% for 2009-2011 and 73% and 86% for 2008 and 2009, respectively (McMahon
et al., 2007; Maupin and Ivahnenko, 2011).
Loading estimates from the Discharge Monitoring Report (DMR) Pollutant Loading Tool,
which was developed by the EPA were used in the model. Loading is estimated by the tool
from data available through the PCS-ICIS. The tool identifies outliers and estimates loads
for locations where data is not reported. The tool also estimates concentrations for facilities
that are likely to discharge nutrients based on their Standard Industrial Classification (SIC).
Average loading rates for each SIC are estimated regionally by the tool and used to estimate
loads for facilities where data is missing or not reported.
Loading estimates of total phosphorus and total nitrogen were downloaded from the
DMR website. Facility discharge locations were mapped to the nearest model river cell based
on the coordinates provided with the data. Facility receiving waters were verified through
permit data available from the EPA’s web page Enforcement and Compliance History Online
(ECHO). These steps are summarized in Figure 1.1.

1.2.3.2

Atmospheric Loading

Atmospheric data was available from three networks that monitor atmospheric deposition of
various compounds: the National Atmospheric Deposition Program (NADP) (NADP , 2007),

4

DMR
Locations
and Loads
Plot Points and
Move to River
Cells

Convert Points to
Rasters

Point Source
Loading
Raster
Figure 1.1: Flow chart summarizing process for estimating point sources.

a wet deposition monitoring program in the United States, the Great Lakes Precipitation
Network (GLPN) (Lisa Bradley, personal communication, 2011) a joint project between the
United States and Canada established to monitor wet deposition to each of the Great Lakes,
and the U.S. Environmental Protection Agency’s dry deposition monitoring network called
Clean Air Status and Trends Network (CASTNET).
Weekly data was downloaded from CASTNET which reports fluxes of nitrate, ammonium, and nitric acid by dry deposition (particle velocity is estimated along with the data).
Wet deposition from the Great Lakes Precipitation Network was obtained from the Canada
Centre for Inland Waters which collects samples monthly and measures Total Kjeldahl Nitrogen (TKN) and nitrate/nitrite. TKN includes both organic and ammonia nitrogen, so TKN
and nitrate/nitrite were summed to obtain total nitrogen. National Atmospheric Deposition
Program data was downloaded as monthly averages of weekly samples (to be consistent with
GLPN data). NADP reports only nitrate and ammonia.
Once concentrations were determined, measured precipitation volume and sampler area

5

were used to calculate total nitrogen deposition (kg/m2 /year). For wet deposition, loading
was calculated by multiplying the concentrations by precipitation values for each month
of available data. Monthly readings were classified as Winter-Spring (December-May) and
Summer-Fall (June-November), and an average was calculated for each summary period for
data after 2005. Averages were then interpolated to the model grid using Simple Kriging.
The same approach was used for dry deposition which is reported weekly. Loads for
each summary period were averaged for data after 2005. Because a limited number of data
points were available for dry deposition of nitrogen and wet deposition of phosphorus, inverse
distance weighting was used to develop loading prediction estimates for those species.
Several studies have shown that organic nitrogen can contribute about a third (median
of 30%) of the total nitrogen atmospheric deposition (Neff et al., 2002). NADP sites only
measure inorganic nitrogen species while the GLPN data includes both organic and inorganic nitrogen. Data from the Sibley GLPN site was compared to data from the nearby Isle
Royale NADP site to assess contribution of organic nitrogen to total atmospheric nitrogen
deposition in the Great Lakes. The Sibley site is in the Sleeping Giant Provincial Park on
the Sibley Peninsula, Ontario, while the Isle Royale site is on the western side of Isle Royale
National Park (Figure 1.2). Both sites are within 50 kilometers of each other, located near
the Lake Superior shoreline. Concentration data from overlapping months were compared.
On average, the concentration at the Sibley site was 24% greater than at the Isle Royale site.
Annually, the percent difference ranged from 14.5 % to 44.5%. This agrees well with studies
of the contribution of organic nitrogen to total atmospheric nitrogen deposition. Neff et al.
(2002) compared data from several atmospheric deposition studies and found that 60% of
reported measurements had a range of 10% to 40% dissolved organic nitrogen with an average of 33%. Nicholls et al. (2001) measured atmospheric nitrogen deposition in Harp Lake,
in southern Ontario, and estimated a dissolved organic contribution of approximately 27%
to total annual nitrogen loading. To account for the organic nitrogen contribution, the average annual percent difference between the concentration at the Sibley and Isle Royale sites

6

Sibley (GLPN)
Isle Royale (NADP)

#
"
*
)

#
*
# # #
* * *
# #
* *
#
*
#
*
# #
* *
# # #
* * *
# # # #
* * * *
"
)
#
*
##
**

#
*
# #
* *
# #
* *
#
*
#
#
"
*
*
)
#
"
*
)
"
)
#
"
*
)
# ##
* **
#
*
# * *
" # #
)
#
*
#
*
"
)
# # )
* * "
"
)
#
# # #
*
* * *
" *
) #
#
*
"
)
###
***
"
)
"
)
" * *
) # #
#
"
*
)
#
*
# #
"
* *
)
"
)
#
*
"
)
#
*
#
* "
"
)
"
)
"
#
*
)
)
# #
"
* *
)
#
"
*
)
#
"
*
)
* GLPN Wet Deposition Sites
#
# "
"
* )
)
#
*
* NADP Wet Deposition Sites
#
#
*
" CASTNET Dry Deposition Sites
)

Figure 1.2: Location of atmospheric deposition monitoring networks. The Sibley and Isle
Royale sites were compared to determine the contribution of organic nitrogen to the total
wet deposition load. For interpretation of the references to color in this and all other
figures, the reader is referred to the electronic version of this thesis.

was assumed to represent the regional proportion of total nitrogen deposition contributed
by organic nitrogen. Concentrations of all NADP sites were adjusted to account for this
calculated percent difference. This process is illustrated in Figure 1.3.

1.2.3.3

Septic Tank Loading

No statewide database of septic tank locations or treatment plant service areas currently
exists in Michigan. The location of septic tanks was estimated using existing information on
treatment plants and water wells. An inventory of operating waste water treatment plants

7

GLPN
(TON+TIN)
Table

NADP
(TIN)
Table

CASTNET
Feature

Subtract Obs. From
paired sites
Interpolate

Estimate
of TON
Table

Dry Atm.
Deposition
Raster

Add TON to
NADP
NADP
(TIN+TON)
Table

Deposition
Site
Locations

Join

Wet
Deposition
Feature
Interpolate
Wet Atm.
Deposition
Loading Raster

Sum

Atm.
Deposition
Raster
Figure 1.3: Flow chart summarizing process for estimating atmospheric deposition.

8

Log of Surveyed Population
Served (CWNS, 2008)

6
R2 =.78

1:1

5
4
3
2
1
1
2
3
4
5
6
Log of Estimated Population
Served (Census, 2010)

Figure 1.4: Comparison between the population served by treatment plants based on
EPA’s CWNS survey and the population served based on the estimated treatment plant
service areas and the 2010 Census.
and collection systems was obtained from current NPDES permits, the Clean Watershed
Needs Survey (CWNS-a national survey of watershed needs conducted by the EPA), and the
Southeast Michigan Council of Governments (SEMCOG). For plants and collection systems
without plant service area data, an area around each plant/collection system was designated
to represent where residents receive sewage collection. This service area was estimated from
the location of the treatment plant and the incorporated area of the town served by the
treatment plant (from Census 2010). The density of residential wells in the area was also
considered; an area with a high density of residential wells was assumed to be outside of a
treatment plant service area even if it was inside an incorporated area.
The population served by each facility was estimated by intersecting the estimated service
areas and the 2010 Census blocks. This was compared to CWNS which provides an estimate
of the population receiving collection from each facility. A comparison of population served
reported in the CWNS and the population served estimated using the method described
above is shown in Figure 1.4.
It was assumed that all households outside the designated waste treatment service areas

9

have a septic system or some other form of on site waste water disposal. The number of
septic tanks in a given area was determined from the number of households reported in the
2010 Census blocks. Initially, septic tanks were placed in census blocks where residential
wells existed according to the Michigan Wells Summary Database (MDEQ, 2005). However,
the number of residential wells in the Michigan Well Database is only a fraction of the total
number of septics that exist in the state. The 2006 National Land Cover Database (Fry et al.,
2011) was used in conjunction with a roads layer from the 2010 Census to determine potential
areas where the remaining septic tanks estimated for the census block were most likely to
be located. The remaining septic tanks (the difference between the number estimated and
the number distributed to well locations) were randomly distributed to NLCD urban cells
within each census block that were between 10 and 60 meters from the center-line of roads and
at least 15 meters from surface water using the Geospatial Modeling Environment (GME)
(Beyer , 2012). This was done because residences are generally located adjacent to roads, and
most Michigan county health departments require a buffer between surface water and septic
tanks. If there were no urban cells in a census block, septic tanks were placed randomly in
other land cover types excluding water and wetlands.
Once septic systems were located, they were assigned an estimated nutrient load based
on the average household size and a yearly estimate of per person nutrient loading rates to
septic systems (USEPA, 2002a). A factor was applied to vacant and seasonal households (0
for vacant and 0.25 for seasonal) to account for the fact that these systems do not receive
loading during the entire year. Figure 1.5 is a flow chart illustrating this process.

1.2.3.4

Non-agricultural Fertilizer

Two forms of non-agricultural fertilizer use were considered in the model: the amount applied to golf courses and the amount applied to residential and commercial lawns. While
county level inventories of agricultural fertilizer purchases are readily available for states
through the agricultural census, data for non-agricultural fertilizer use is limited. This issue

10

2006
NLCD
Raster

Apply Condition:
Urban Cells=1,
else 0

Urban
Raster

Convert to
Feature

Urban
Features

Incorp.
Area
Features

Roads
Features
Manually
Designate
Service
Areas

CWNS
Data

Intersect

Service
Area
Features
Surface
Water
Features

MI Well
Features
2010 Census
Households
Feature

Sum Wells in Blocks
and Subtract from
Households

Unknown Septic
Locations
Feature

# of Septics with
Unknown
Location Feature

Generate
Random
Points in
GME
Intersect

Per Person
Excretion Rate
Value

Multiply

Septic
Locations
Feature
Septic
Loading
Feature

Septic
Tank Area
Features

2010 Census
Avg. Household
Size Feature

Convert to
Raster

Septic
Loading
Raster

Figure 1.5: Flow chart summarizing process used to estimate loading from septic tanks.

11

was confronted by Ruddy et al. (2006) in their county level inventory of nutrients. They
estimated both the portion of fertilizer sales resulting from non-agricultural fertilizer use
and the portion of the state level non-agricultural fertilizer used in each county. State and
county non-agricultural fertilizer use in pounds of nitrogen and phosphorus were estimated
from 1987 to 2001.
In addition to estimating non-agricultural fertilizer use for the state, Ruddy et al. (2006)
developed a model to distribute statewide fertilizer sales data to each county based on
population density. By plotting the available county level data and population density,
they concluded that non-agricultural fertilizer use increased with population density until
a threshold of about 700 people/km2 , above which fertilizer use did not increase. Their
estimates suggest that for Michigan, the total non-agricultural fertilizer use has not changed
significantly since 1996, therefore, the 2001 estimate of total non-agricultural phosphorus
and nitrogen for Michigan was assumed to be reasonable for years 2005-2012. This estimate
was redistributed to counties based on the relationship described above and the 2010 Census
population data to obtain a county level estimate of non-agricultural fertilizer use for the
study period.
Once county level estimates were obtained, golf course areas were established using a
directory of golf courses from the Michigan Golf Association and manual digitization of
boundaries using aerial imagery. Addresses were geocoded and polygons were manually digitized around each course using satellite images. The golf course polygons inevitably included
areas that are not fertilized including cart paths, sand traps, and forested landscapes. From
the most compact courses (those with the largest ratio of greens to golf course area, ˜10%
of the total) an average area per hole was estimated. The proportion of fertilizable area for
each golf course polygon was then estimated based on number of holes. An average course
fertilization rate was obtained from locally (Midwest) observed rates of course application
from a survey conducted by the Environmental Institute for Golf of the Golf Course Superintendents Association of America (Throssell et al., 2009). It was visually estimated that

12

even for the most compact courses (those lacking tree parks and lakes) only about 85% of
the total course area is fertilizable. Fertilizer was randomly distributed to 85% of the area
within the manually designated golf course areas at the surveyed rate.
The estimated fertilizer contribution for golf courses was removed from the nonagricultural fertilizer estimated for each county by Ruddy et al. (2006). This adjusted amount
was then distributed to census tracts following the population density model outlined by
Ruddy et al. (2006). After per tract fertilizer use was estimated, potential fertilizable area
within each census tract was established using the 2006 NLCD Land Cover and Impervious
Area data layers, Michigan roads, and the manually generated layer of golf course areas. It
was assumed that commercial and residential fertilizer is predominantly applied in urban,
pervious areas that are close to roads. It was also assumed that most home and business
owners only apply fertilizer to parts of their property that are adjacent to roads. Areas available to receive residential or commercial fertilizer were established by selecting all NLCD
2006 medium to high density urban cells and creating an inclusion buffer of 30 meters from
road polylines but excluding the portion of the buffer occupied by the right of way and
the road itself (20 meters from any polyline). Additionally golf course areas were removed
from the selected cells since application to golf courses was already considered (as described
above). The cell areas were multiplied by the NLCD pervious area product to establish what
proportion of the cell was occupied by fertilizable land.
Literature values (Zhou et al., 2008; Law et al., 2004) were used to establish commercial and residential application rates of nitrogen and phosphorus. Recommended rates over
the lifetime of a lawn and for various residual phosphorus levels were averaged. While all
the selected cells could potentially be fertilized, applying fertilizer to all cells grossly overestimates fertilizer application. This was determined by comparing county level totals for
non-agricultural fertilizer from Ruddy et al. (2006) and county levels estimated by applying
fertilizer at the recommended rate to all potentially fertilizable cells. In order to maintain locally reasonable application rates while honoring the observed totals, fertilizer was

13

randomly applied to cells, at the recommended rate, until the total application matched
observed county totals. This process is summarized in Figure 1.6.

1.2.3.5

Animal Manure

In recent decades there have been significant changes in agricultural management practices.
Animal production has shifted from small farms where animals are primarily raised on pasture land to concentrated animal feeding operations (CAFOs) where large numbers of animals
are raised in confinement (Cheng, 2003; Kellogg et al., 2000). Waste produced by confined
operations is collected and distributed to surrounding cropland. As CAFO size increases, an
ever increasing area of cropland is needed to assimilate waste, and facilities are required to
transport waste to neighboring farms (Kellogg et al., 2000). In order to accurately describe
waste management of large operations as well as small farms, two approaches were used to
distribute animal manure to model cells.
The Michigan Department of Environmental Quality (MDEQ) has data on the location
and inventory for approximately 200 of the largest confined feeding operations in the state.
This information, was combined with the 2007 Agricultural Census data to determine the
approximate location and extent of large confined animal feeding operations. Facilities and
animal inventories for operations described by the MDEQ were removed on a county basis
from the 2007 Agricultural Census data in order to determine the inventory of the remaining
unknown farms. After those facilities were removed, the average herd/flock size of a particular animal type and inventory level was re-calculated for each county based on the remaining
number of farms and animal inventories.
Once the average county inventory was established, farms in each county were classified
as confined operations and pasture operations based on a threshold inventory level. If the
inventory was above a certain level, it was assumed that the farms likely confine their animals
for the majority of the time. The thresholds, taken from Kellogg et al. (2000), were varied
based on animal group. For instance, dairy farms have no threshold since dairy farming

14

Tract Pop.
Feature

Golf Course
Feature

Calculate
Pop. Density

Pop. Density
Model (Ruddy
et al 2006

Tract Pop.
Density

Tract
Residential.
Fert Feature
Feature to
Raster
Fert
Tract
Raster
divide
% Tract
Fertilized
Multiply
Divide

Course
Fertilizable
Area Feature

Calculate
Model

Avg Golf
Course
Fert Rate
Course Fertilization
Feature

Multiply

Tract Total
Non Ag Fert
Subtract

Calculate Course Fert
Area from # of Holes

Zonal Sum
2006 NLCD
Raster

Golf Tract
Fert. Totals

Urb
Raster

Tract Fert
Area

Fert.
Rates

Zonal
Sum
Tract
Perv Area
Adjusted
% Tract
Fert

Golf
Raster

Fert. Area
Mask

Multiply

NLCD
Perv

Buffer

Apply Condition:
Urban Cells=1,
else 0

Zonal Sum

Multiply
Random
Raster

Multiply
Con, mask<%
Fert
Selected
Fert. cells
Multiply
Residential
Fert Loading

Roads

Feature
to Raster

Random
Fert Mask

Roads
Buffer
Feature
Feature to
Raster
Roads
Buffer
Raster
Randomly
Distribute
Golf Fert to
Golf Cells
Golf Fert
Loading
Add
NonAg
Fertilizer
Raster

Figure 1.6: Flow chart summarizing process used to estimate non-agricultural fertilizer.

15

TP Animal Nutrient Excretion
Model

P as excreted (lb/day)

N as excreted (lb/day)

TN Animal Nutrient Excretion
Model
0.6
y = 0.0003x + 0.0053
0.5
R 2 = 0.8845
0.4
0.3
0.2
0.1
0
0

1000
Weight (lbs)

y = 5E-05x + 0.0041
R 2 = 0.7157

0.12
0.08
0.04

2000

0
0

1000
Weight (lbs)

2000

Figure 1.7: Linear regressions of animal weight and nutrient excretion rates. These linear
models were used to predict animal excretion rates for animals not discussed in the NRCS
Waste Management Handbook (NRCS , 2008) including bison, elk, emus, pheasant and
others.
requires partial to complete confinement of animals.
Nutrient loads were calculated for each herd/flock type, in each county based on estimated
nutrient excretion rates for each animal type (NRCS (2008) estimates excretion rates for
common animal groups). From available excretion data and average animal weight (excluding
lactating cows), a linear model was developed to estimate excretion rates based on animal
size (Figure 1.7). This was used for animal groups in the agricultural census that did not
have a nutrient excretion rate provided by the NRCS Waste Management Handbook. Using
excretion rates and farm inventory, a total farm load was calculated.
Inevitably, during collection, transport and storage, a portion of the load is lost before it
reaches the fields. An estimate of this efficiency (Kellogg et al., 2000) was used to calculate
the spreadable manure load. Spreadable acreage for each farm was calculated from the
proportion of the total county animal waste produced by the farm and the total amount of
acreage receiving manure in the county (reported in the 2007 Agricultural Census).
While some location data was available for confined operations from MDEQ, the location

16

of the remaining confined operations was unknown. Due to fuel and labor costs, there is
a limit to the distance a farmer is willing to transport waste (Kellogg et al., 2000). This
distance was estimated to be about three miles. To determine the most probable location
of operations not included in the MDEQ data, model cells were selected based on available
cropland within 3 miles of the cell and compared to the spreadable acreage calculated for
each farm. If enough land was available, the cell was selected as a possible location for a
particular farm. Farms were then randomly distributed to selected cells using GME.
After farm locations were selected, buffers were created around each farm designating the
land receiving manure. Initially, buffer radius was estimated from the farm’s proportion of
manure acreage in the county. However, some of the land in each buffer was occupied by roads
or non-agricultural land reducing the land available for manure application. Therefore, the
radius of the buffer was iteratively increased until the buffer incorporated enough fertilizable
land to meet the actual spreadable land estimated for each farm. This ensured that the
total observed acreage fertilized by manure in a county was honored, but roads and nonagricultural cells did not receive fertilizer. Buffers were adjusted until at least 98% of the
farms met their estimated manure acreage. As buffers increased, neighboring farm manure
acreage overlapped, and thus it was necessary to sum overlapping buffers, to avoid cells
receiving multiple fertilizer contributions.
The five year crop rotation from the Cropland Data Layer (CDL), a remotely sensed
product from the National Agricultural Statistics Service (USDA, 2011), was used to estimate
the rate at which manure was applied in each cell within a farm’s buffer. The CDL is
available for Michigan for 2007-2011. Recommended fertilizer rates from extension literature
(Warncke et al., 2004) were calculated for each crop type in the CDL; the recommended rate
was averaged over the five year period to calculate the average rate of application over the 5
year rotation period. Within each farm’s buffer, a proportion of the total manure generated
by the farm was applied to each cell based on the fertilization rate estimated from the 5 year
rotation. This allowed manure to be over-applied relative to the recommended rate based on

17

Huron County Total Nitrogen Manure
kg/ha/year
1235
617
309
60

Figure 1.8: Applied manure predicted by the model for Huron County, MI.

the farm’s available manure. Other studies (Kellogg et al., 2000) have calculated that manure
generation in several counties exceeds the nutrient demands of available cropland; therefore,
it is likely that many farms apply manure above the average recommended rate. The nonspreadable portion of the manure load (the portion lost during transport and storage) was
applied within the cell at the farm, assuming that this portion becomes part of the farm
runoff, generating locally concentrated sources of nutrients (see Figure 1.8).
For farms designated as unconfined, waste was only applied to NLCD pasture land.
Pasture land application rates were determined within each county based on the total waste
produced and the total pasture area within a county. The processes used to estimate manure
application are shown in Figures 1.9 and 1.10.

1.2.3.6

Commercial Fertilizer

It was assumed that cultivated land that did not receive animal manure application, received
chemical fertilizer. A linear regression was fit to 2001 nutrient application estimates of
18

2007 Ag Census
Farms/Manure
Area/Inventories

County
Features

2006
NLCD
Raster

MDEQ CAFO
Locations Table
Geocode
Addresses

Join
County Farms/
Inventories/Manure
Area Feature

Registered
CAFO/Inventory
Features

Subtract and
Recalc. Avg.
Inventories

Ag Land
Raster

Determine # of
CAFOs/nonCAFOs

Unreg.
CAFOs
Feature

Apply Condition:
Ag Cells=1, else 0

County nonCAFO Farms
Feature

Distribute Random
Points (GME)

Unreg.
CAFO
Feature

Pasture
Land
Raster

Est. of Inventory
Requiring
Confinement

Unregistered
County Farms/Avg.
Inventories Feature
2006 NLCD
Raster

Apply
Condition:
Pasture
Cells=1,
else 0

Distribute Random
Points (GME)

Pasture
Farms
Feature
Merge

Excretion
Rates

Farm/Inventory
Feature
Farm
Loads
Feature

Multiply

Figure 1.9: Manure Application Part I: Flow chart summarizing the process used to
estimate the location and inventories of farms that produce manure.

19

Farm Loads
Feature (see
Figure 1.9)
Normalize
Normalized
Farm Loads
Feature

County Manure
Area Feature (see
Figure 1.9)

Intersect and
Multiply
Manure Area
per Farm
Feature

Farm Manure
Spreadable
Area Feature

Buffer
Increase
Buffer Radius
Based on
Required
Farm AreaBuffer Ag
Area

Farm
Manure Area
Feautre

Join

Yes

Calc. Ag Area
in Buffer

Join

Enough Ag
Area?

No

CDL Normalized
Average Fert. Rate
Raster (see Figure
1.11)

Est. of
Unrecovered
Load

Unrecovered
Load Feature
Recovered
Load Feature
Convert to
Raster

Multiply

Add

Unrecovered
Load Raster

Recovered
Load Raster

Manure Area
Application
Rate Raster

Convert to
Raster

Manure Raster

Figure 1.10: Manure Application Part II: Flow chart summarizing the process used to
estimate the rate and extent of manure application.

20

chemical fertilizer for Michigan counties (Ruddy et al., 2006) and NLCD cultivated area for
2001. This relationship was used to estimate county nutrients from agricultural chemical
fertilizer for 2005-2012 using the 2006 NLCD land cover. The county level estimates of
chemical fertilizer were then distributed to available cropland area within the county that
did not receive manure.
Similar to manure, application rates of chemical fertilizer were estimated from the total
observed chemical fertilizer in the county and the relative recommended application rate
based on a cell’s 5 year crop rotation. Due to variable environmental factors (climate, soil
types) and local social practices, it is likely that fertilization rates for a single crop type are
variable across the state. While this can create an artificial gradient at the county border,
the estimate honors observed fertilizer consumption for the state. This process is shown in
Figure 1.11.

1.2.4

Sample Collection

Flow measurements and nutrient samples were collected from the fall of 2011 to the spring
of 2012 in two field campaigns intended to capture synoptic conditions during baseflow and
melt. Sampling focused on 4 major watersheds in the Lower Peninsula of Michigan: Grand
River in south central Michigan, Saginaw Bay in eastern Michigan, Muskegon River in north
central Michigan, and Boardman-Charlevoix in northwest Michigan. Between 15 and 30 sites
were sampled within each watershed for a total of 90 sites. Sample locations were chosen in
an effort to evenly distribute samples throughout the watershed. This was accomplished by
selecting points so that the area upstream of each sample point would have, approximately,
the same incremental watershed area. The watersheds for the combined sampling covers
over 40% of the land area of the Lower Peninsula.
A Sontek RiverSurveyor S5 Acoustic Doppler Profiler (ADCP) was generally used to
measure flow in all but low flow streams, which were measured using a Marsh-McBirney
Flo-Mate. Samples were also collected near USGS stream gauges, and data from gauges was

21

CDL-5
Years
Raster

Recommended
Fert. Rate
Table

Classify cells
by Fert. Rates
5 Year
Fertilization
Rasters

County
Estimates of
Ag Fert. Table

Average Rate
over 5 Years

County
Features

Join

Average 5 Year
Fertilization
Raster

County Fert.
Features

2006 NLCD
Raster

Normalize
Rates

Feature to
Raster

Apply Condition:
Ag Cells=1, else 0

Normalized
Average Fert.
Rate Raster

County Fert.
Raster

Agricultural
Land Raster

Multiply
Ag. Chemical
Fertilizer
Raster
Figure 1.11: Flow chart summarizing processes used to estimate loading from commercial
agricultural fertilizer.

22

Urban
Unmanaged
Cropland

Figure 1.12: Sampled watersheds with dominant land use indicated.

downloaded to coincide with the time of the sample collection. Samples for total nitrogen
and total phosphorus were collected at the time of flow measurement.
Samples were collected in-stream or from bridge crossings and then frozen on dry ice.
The samples were analyzed by the Michigan State University Algae Lab using second derivative spectroscopy (TN) (Crumpton et al., 1992) and ascorbic acid methods (TP) following
persulfate digestion (Standard Methods 4500-P.E. and 4500-N.C). Figure 1.12 shows the incremental watersheds for each sample point with the dominant land use indicated for each
watershed. Flow and sample data are provided in Appendix A.

23

1.3

Results

1.3.1

Sources

The results from each source model are shown in Figures 1.13 and 1.14. Total nitrogen loading to the Lower Peninsula is estimated to be an average of 23 kg/ha/year; total phosphorus
loading is estimated to be an average of 3 kg/ha/year.
The proportion of the total load contributed by each land use is shown in Figure 1.15.
Chemical agricultural fertilizer contributes 39% of nitrogen and ranges from nearly 62% of
the total nitrogen source in cropland to 0% in unmanaged land. The largest overall source
of phosphorus is also commercial fertilizer. The highest rate of phosphorus and nitrogen
loading is estimated to be a result of confined animal feeding operations. This situation
arises due to a lack of available land to distribute animal waste and the unrecovered portion
of waste generated during confined animal feeding operations that becomes part of farm
runoff. Beaulac and Reckhow (1982) estimate that nutrient export from confined animal
feed lots may be 2-3 times greater than runoff from other agricultural land and may even be
the same order of magnitude as municipal treatment plant loads.
Figure 1.16 shows total loading of phosphorus and nitrogen as well as the ratio of total
phosphorus to total nitrogen. Total phosphorus loading by atmospheric deposition is low
relative to nitrogen. In the north, where atmospheric deposition is the main source of
nitrogen, the ratio of phosphorus to nitrogen is very low; total nitrogen loading rates are
as much as 20 times greater than phosphorus loading rates from atmospheric deposition.
However, in agricultural areas, the estimated ratio of total phosphorus to nitrogen is much
higher, particularly in areas around confined animal feeding operations.

1.3.2

Rates by Land use

Loading rates varied significantly between land uses. The average loading rate by land use
is shown in Table 1.1. The largest loading rates were estimated for cropland, which had

24

Atmospheric Deposition
(Atm)
kg/ha/year
19
18
16
14

Septic Tanks
(Septic)
kg/ha/year
24
12
5
1.5

Manure Application
(Manure)

kg/ha/year
300
180
85
21

Agricultural Chemical
Fertilizer (ChemAg)
kg/ha/year
90
50
30
8

Point Sources
(Point)

Non-Agricultural
Fertilizer (NonAg)
kg/year

kg/ha/year
30
18
10
3

100,000
50,000
20,000
12,000

Figure 1.13: Total nitrogen application rates for each modeled source.

an average nitrogen loading rate of over 100 kg/ha/year, 2 times greater than pasture land
and 3 times greater than urban land. Phosphorus loading rates estimated for cropland are
6 times greater than pasture land and 15 greater than urban land. The dominant nutrient
source in unmanaged land is atmospheric deposition which is approximately 16 kg/ha/year
for nitrogen and 0.2 kg/ha/year for phosphorus.
There is a high degree of variability in loading rates within land use classes. Box plots
showing the range in average loading rates within sampled watersheds are shown in Figures 1.17 and 1.18. The results show that loading of chemical fertilizer within agricultural
land can range from less than 0.01 kg/ha/day (3.65 kg/ha/year) to nearly 0.50 kg/ha/day

25

Land Use
Ag
Pasture
Urban
Unmanaged

TN
TP
kg/ha/year
127
19
52
2.8
39
1.2
16
0.2

Table 1.1: Nutrient loading rates by land use.

(182.5 kg/ha/year) amongst sampled watersheds for nitrogen (as a reference, the agricultural
census estimated an average rate of 130 kg/ha/year of nitrogen application for corn). This
illustrates the potential bias introduced to source estimates derived from land use alone, particularly across areas with diverse crop types. In Michigan, cropland in the southern areas
are dominated by corn and soy while in the north, orchards, which have lower fertilization
rates, are more prevalent.

1.3.3

Loading Rates and In-stream Concentration

In order to assess the relationship between loading and in-stream water quality, a linear model
was fit to average loading rates within watersheds and total nitrogen and total phosphorus
concentrations under baseflow and melt conditions. The form of the linear model is

Ci = α1 P ointi +α2 ChemAgi +α3 M anurei +α4 N onAgF erti +α5 Septici +α6 Atmi +B (1.1)
where C i is the concentration (µg/L) of total phosphorus or total nitrogen measured for
watershed i, each alpha is a model coefficient, each source term is the average loading rate
(kg/ha/day) in watershed i, and B is the model intercept. The results of each linear model
are shown in Table 1.2.
The R2 was 0.50 for the nitrogen melt model and 0.29 for the nitrogen baseflow model.
The phosphorus melt model had an R2 of 0.06 and the baseflow model R2 was 0.24. The
results of the models are plotted with observed values in Figure 1.19.

26

Source
Atm
ChemAg
Manure
NonAg
Septic
Point
Intercept (B)
Atm
ChemAg
Manure
NonAg
Septic
Point
Intercept (B)

TN Melt
α
p-value
1.66E+06
0.95
2.00E+07 8.34E-05
3.25E+06
0.37
7.64E+07
0.49
-2.50E+08
0.13
3.44E+08
0.22
7.17E+02
0.55
TP Melt
4.07E+08
0.51
2.69E+05
0.88
2.85E+05
0.82
4.15E+07
0.54
-4.32E+06
0.94
1.98E+07
0.63
2.98E+01
0.82

TN Base
α
p-value
1.12E+07
0.23
-3.12E+05
0.84
3.21E+06 7.17E-03
9.12E+07 1.19E-02
3.21E+07
0.16
-7.52E+07
0.72
-1.10E+02
0.78
TP Base
1.64E+09 7.93E-04
-4.32E+05
0.75
8.83E+05
0.35
6.46E+07
0.22
-5.41E+07
0.20
-1.07E+07
0.74
-2.90E+02 3.95E-03

Table 1.2: Results of linear regression models of source loading and observed in-stream
concentration. The α value represents the fitted coefficient. Bold values indicated
parameters with a p-value ≤0.05

The p-values which indicate the significance of each model parameter show that commercial fertilizer is the only significant predictor of nitrogen concentration in watersheds. None
of the sources were significant in the phosphorus model. When observations with concentrations above 200 (µg/L) are removed from the total phosphorus baseflow observations, the R2
improves to 0.46 and manure, atmospheric deposition, and non-agricultural fertilizer have
significant coefficients.
The source inputs for this linear model include source loading throughout the entire
watersheds. In larger watersheds, sources applied at large distances from surface water may
never be delivered to the stream due to attenuation of nutrients within the watershed. A
second linear model was constructed using source inputs applied within a 1 km flow distance
from surface water. Flow distance was calculated using the ArcGIS Hydro Toolbox and the
National Elevation Dataset (NED) (Roth and Dewald , 1999). The results are shown in Table
1.3 and Figure 1.20. This approach slightly improved the R2 value for the melt models but

27

Source
(within 1 km of surface water)
Atm
ChemAg
Manure
NonAg
Septic
Point
Intercept (B)
Atm
ChemAg
Manure
NonAg
Septic
Point

TN Melt
α
p-value
2.71E+07
0.59
3.90E+07 1.78E-04
1.05E+07
0.17
6.87E+07
0.74
-5.47E+08
0.10
3.93E+08
0.15
4.74E+02
0.61
TP Melt
1.23E+09
0.20
-2.27E+06
0.56
1.01E+06
0.73
5.61E+07
0.65
9.51E+06
0.93
2.60E+07
0.53

TN Base
α
p-value
8.76E+06
0.61
4.53E+05
0.89
7.86E+06 3.56E-03
1.34E+08
0.06
-9.05E+07
0.42
1.03E+07
0.91
2.14E+02
0.49
TP Base
1.64E+09 7.93E-04
-4.32E+05
0.75
8.83E+05
0.35
6.46E+07
0.22
-5.41E+07
0.20
-1.07E+07
0.74

Table 1.3: Results of linear regression models of source loading within 1 km of surface
water and observed in-stream concentration. The α value represents the fitted coefficient.
Bold values indicated parameters with a p-value ≤0.05
decreased the R2 for the baseflow models.

1.4

Discussion

An understanding of the relationship between nutrient source loading and in-stream water
quality is not only necessary to effectively manage surface water quality but is also fundamental to predicting how water quality will change with alterations of land use and land
management. Average loading rates of nutrients to watersheds in the Lower Peninsula are
highly variable ranging from 0 kg/ha/year to over 200 kg/ha/year of nitrogen and over 40
kg/ha/year of phosphorus. These high rates of loading are estimated to result from confined animal feeding operations. Manure is the second largest source of nitrogen and third
largest source of phosphorus. As animal production shifts to increasingly consolidated operations, nutrients from manure have the potential to become a major concentrated source
of nutrients in agricultural watersheds (Chang et al., 2003; Kellogg et al., 2000). Within all

28

land use types, except unmanaged land where the major source is atmospheric deposition,
loading rates for major sources exhibited at least 1 order of magnitude of variability. This
illustrates the difficulty in assigning loading rates of nutrients based on land use alone in
diverse watersheds.
While these estimates provide an indication of the source of nutrients to the landscape,
estimating in-stream concentrations from total sources of nutrients is more complicated.
Attenuation processes occurring within watersheds alter nutrient loads as they travel within
the basin to surface water. The results of the linear model suggest that with source estimates
alone, less than 50% of the variability in concentrations can be explained (the nitrogen
melt model had the highest R2 of 0.50). The results indicate that phosphorus sources are
a poor predictor of melt phosphorus (R2 was only 0.06). In addition, during baseflow,
concentrations increase linearly with increasing sources for watersheds with concentrations
less than 200 µg/L but watersheds with concentrations above 200 µg/L are not characterized
by greater rates of source application. These results suggests that other factors are important
to estimate watershed phosphorus concentrations, including biological processes, processes
that occur in-stream, as well as legacy and natural sources of phosphorus. In all but the
nitrogen baseflow model, the coefficient for septic was negative and in both baseflow models,
the coefficient for point sources was negative. While this may seem counter intuitive, it may
be the result of dilution by these sources since they are associated with large amounts of
waste water.
Both phosphorus and nitrogen melt models improved when only sources within 1 km of
surface water were considered. However, sources within 1 km of surface water provided poorer
predictions of in-stream concentrations during baseflow. There are several explanations for
this effect including legacy sources of nutrients and the influence of nutrients in groundwater
during baseflow. It is important to note, however, that the differences between models was
only slight, and it appears that additional factors other than source distance control source
delivery. The relationship between nutrient sources, basin attenuation, and observed nutrient

29

loads is explored further in CHAPTER 2 of this thesis.
This methodology provides a means to investigate locally important concentrated sources
of nutrients such as golf courses and septics. Source estimates not only provide the potential to investigate the relationship between nutrients and concentration but can be used to
examine relationships between water quality and the sources of other pollutants, particularly those associated with waste. For instance, Verhougstraete (2012) compared estimates
of septic loads developed using these methods to the occurrence of B. thetaiotaomicron (B.
theta) a human marker for fecal contamination, in surface water samples from the Lower
Peninsula.
Since data related to sources of nutrient loads is limited at the spatial scale used in this
research, it is difficult to determine the error associated with some of these estimates. The
accuracy of the estimates is in large part based on the the accuracy of the remotely sensed
products used to distribute many of the sources including the 2006 NLCD and the 2007-2011
Cropland Data Layers. Additionally, random distributions of points were used to estimate
the locations of unknown septics and confined animal feeding operations and to distribute
loads within receiving areas. Since manure is the second largest source of phosphorus and the
third largest source of nitrogen, this is a potentially important source of error. Additional
validation of the septic and CAFO locations will improve the confidence of these source
estimates.

30

Atmospheric Deposition
Manure Application
Agricultural Chemical
(Atm)
(Manure)
Fertilizer (ChemAg)
kg/ha/year
kg/ha/year
kg/ha/year
70
17.0
.095
40
8.0
.085
20
4.5
.075
5
1.5
.070

Septic Tanks
(Septic)
kg/ha/year
5.0
2.5
1.0
.25

Non-Agricultural
Fertilizer (NonAg)
kg/ha/year
kg/year
2.8
175,000
1.8
40,000
1.0
10,000
.25
1,000

Point Sources
(Point)

Figure 1.14: Total phosphorus application rates for each modeled source.

31

Total Nitrogen

100%
80%

Point
Septic

60%

NonAg
40%

Manure
ChemAg

20%

Atm

0%
Cropland
100%

Urban Unmagaged
Total Phosphorus

Pasture

80%

Point
Septic

60%

NonAg
40%

Manure
ChemAg

20%

Atm

0%
Cropland

Urban

Unmanaged

Pasture

All Land Uses

100%
80%

Point
Septic

60%

NonAg
40%

Manure
ChemAg

20%
0%

Atm
TN

TP

Figure 1.15: Proportion of the total nutrient load contributed by each source for the entire
Lower Peninsula and within different land uses.

32

TN kg/ha/year
50
40
30
20

TP kg/ha/year
10
5
3
1

TP:TP
.50
.25
.10
.01

Figure 1.16: Total nitrogen and phosphorus loading in the Lower Peninsula. The ratio of
total phosphorus to total nitrogen is shown in the bottom panel.

33

Loading kg/ha/day

Cropland
0.40

0.20

0.00

Loading kg/ha/day

Pasture
0.30
0.20
0.10
0.00

Loading kg/ha/day

Urban

0.04

0.02

0.00

Loading kg/ha/day

Unmanaged
0.04

0.02

0.00

Atm

ChemaAg

Manure

NonAg

Point

SepƟc

Figure 1.17: Box plots of total nitrogen loading rate variability in sampled watersheds.
Loading rates are subset by land use.

34

Cropland

Loading kg/ha/day

0.08

0.04

0.00

Loading kg/ha/day

Pasture
0.06
0.04
0.02
0.00
Urban

Loading kg/ha/day

0.015
0.010
0.005
0.000
Unmanaged
Loading kg/ha/day

0.0012
0.0008
0.0004
0.0000
Atm

ChemaAg

Manure

NonAg

Point

SepƟc

Figure 1.18: Box plots of total phosphorus loading rate variability in sampled watersheds.
Loading rates are subset by land use.

35

TN Base

6,000

y = 0.5024x + 868.26
R 2 = 0.5024

3,000

Predicted Concentration
(kg/m3 )

Predicted Concentration
(kg/m3 )

TN Melt
9,000

3,000

y = 0.287x + 392.3
R2 = 0.287

2,000

1,000

0

0

0

0
3,000 6,000 9,000
Observed Concentration (ug/L)

3,000

TP Base

y = 0.8227x + 430.65
R2 = 0.0761

800
400

Predicted Concentration
(kg/m3 )

Predicted Concentration
(kg/m3 )

2,000

Observed Concentration (ug/L)

TP Melt
1,200

1,000

0

600
y = 0.235x + 51.483
R 2 = 0.235
400

200

0
0

400

800

1,200

0

Observed Concentration (ug/L)

200

400

600

Observed Concentration (ug/L)

Figure 1.19: Graphs showing linear regression model predictions compared to observed
concentrations.

36

TN Melt
Predicted Concentration
(kg/m3 )

Predicted Concentration
(kg/m3)

9,000

y = 0.5314x + 823.25
R2 = 0.5314
6,000

y = 0.257x + 409.23
R2 = 0.257

2,000

1,000

3,000

0

0

0

0
3,000 6,000 9,000
Observed Concentration (ug/L)
TP Melt

400

1,000

2,000

3,000

Observed Concentration (ug/L)

Predicted Concentration
(kg/m3 )

Predicted Concentration
(kg/m3 )

TN Base

3,000

y = 0.0724x + 124.21
R 2 = 0.0724
200

TP Base

600

y = 0.1689x + 55.929
R 2 = 0.1689

400

200

0

0
0

200

0

400

Observed Concentration (ug/L)

200

400

600

Observed Concentration (ug/L)

Figure 1.20: Graphs showing linear regression model predictions (based on nutrient loads
within 1 km of surface water) compared to observed concentrations.

37

CHAPTER 2
PREDICTING NUTRIENT LOADS AND SOURCES: A SPATIALLY
EXPLICIT STATISTICAL MODEL

2.1

Introduction

Nutrient loading models are one tool used to help quantify how changes in land use, management practices, and nutrient sources will impact in-stream and coastal water quality. While
trends in water quality can be observed through periodic sampling, observations are limited
in spatial and temporal extent by expense and accessibility (Smith et al., 1997; Valiela et al.,
1997; Alexander et al., 2002; Borah and Bera, 2004). Cutbacks to state and federal monitoring programs have led to decreases in the frequency and geographic coverage of samples
(Robertson and Saad , 2011). Models can be used to estimate nutrient loads at unmonitored
locations increasing the spatial and temporal resolution of nutrient loading estimates (Smith
et al., 1997). They can also forecast trends in water quality under scenarios of changing
land use and climate and have been used at various scales to assess the source of nutrients,
target source areas for management, and predict the effectiveness of management strategies.
Models are particularly useful to better understand the role of non-point sources to nutrient
loads (Nikolaidis et al., 1998; Borah and Bera, 2004).
Several nutrient load modeling approaches exist in the literature (as reviewed in Alexander
et al. (2002), Borah and Bera (2004) and Valiela (2002)). These range from semi-distributed
hydrologic models such as the Soil Water Assessment Tool (SWAT) to simple mass balance
accounting (Jaworski et al., 1992; Alexander et al., 2002; Boyer et al., 2002). Common among
many existing approaches is a tendency to lump sources or assume nutrient loading rates
based on factors such as land use or management practices rather than directly accounting for
independent source contributions. Such methods generally assume that changes in land use
or management will lead to a proportional response in source contribution to loading which

38

limits their transferability (Zhang, 2011; Destouni et al., 2006). Detailed descriptions of each
source independent of transport mechanisms are necessary to describe nutrient loads across
highly variable watersheds and predict how nutrient loading will be impacted by changes in
climate and land use.
This work was motivated by the need for a nutrient model that can accurately and
consistently predict nutrient loads at a high level of spatial detail across diverse watersheds.
The model described here is based on spatially detailed and temporally variable descriptions
of all potential nutrient sources in Michigan’s Lower Peninsula (as described in CHAPTER
1) to better understand the impact specific sources have on nutrient loading. These sources
are related to stream water quality through a statistical model that provides insight to the
importance of pathways and processes. This model improves on nutrient loading models
(described below) by enhancing the spatial detail of individual sources at a regional scale
while retaining the ability to predict source contributions to the observed load.

2.1.1

Review of Existing Approaches

Examples of nutrient models developed at scales similar to the Lower Peninsula of Michigan
include the budget approach (Jaworski et al., 1992; Boyer et al., 2002), the use of export
coefficients (Johnes and Heathwaite, 1997; Johnes, 1996), loading functions (Generalized Watershed Loading Functions-GWLF) (Haith and Shoenaker , 1987; Evans et al., 2002; Georgas
et al., 2009) and statistical approaches such as the USGS SPARROW (Spatially-Referenced
Regression on Watershed attributes) model (Smith et al., 1997). These approaches are described below.
A budget approach examines the difference between the nutrient inputs and outputs to
a watershed (Boyer et al., 2002; Jaworski et al., 1992). While applications of this approach
can incorporate detailed sources, they give little information about pathways and processes
(Robertson and Saad , 2011). In its simplicity, this type of model makes it difficult to understand the controls on nutrient reduction that occur as nutrients travel through the landscape.

39

While an attempt is made to identify the net inputs from each source, the relative proportion
of the observed load cannot be attributed to a particular source. Depending on landscape
and stream factors, the proportion of nutrients delivered to a stream from a particular source
will vary (Boyer et al., 2002; Robertson and Saad , 2011). These assumptions make it difficult
to use this type of model to predict loading in the future or in unmonitored locations outside
of the studied watershed.
The export coefficient method attempts to predict the nutrients derived from various
land cover/land use types (Johnes, 1996; Johnes et al., 1996). This method assumes that
loading will be evenly distributed across a land cover and that all land covers will contribute
nutrients at the same rate (Johnes, 1996; Vadeboncoeur et al., 2010; Zhang, 2011). “Export
Coefficients”, the proportion of the initial applied source that is delivered to the stream,
are estimated for various land uses, conditions, and management practices. This method
is ideally suited to smaller watersheds where there is minimal variation across land cover
types. This assumption becomes problematic at larger scales, where land use practices and
land cover conditions, and how they affect nutrient export, are more variable (Zhang, 2011,
2010).
More physically based models such as SWAT and Generalized Watershed Loading Function (GWLF) use a process-based approach to calculate hydrologic and sediment fluxes
within the system. GWLF uses water balance equations to simulate runoff and erosion from
“subareas” within the basin of interest (Haith and Shoenaker , 1987; Georgas et al., 2009)
and estimates the concentrations of dissolved nutrients in runoff and baseflow depending on
the proportion of various land uses in the subarea (Haith and Shoenaker , 1987; Evans et al.,
2002; Alexander et al., 2002). SWAT uses a similar approach but includes additional, partially process-based descriptions of nutrient cycling and land management (Grizzetti et al.,
2005; Arnold et al., 1998). SWAT has been applied with a highly detailed accounting of
agricultural sources and management practices, but other sources such as urban fertilizer
are less detailed, and atmospheric deposition is completely excluded (Valigura et al., 2001).

40

Since other models have identified atmospheric deposition as a major source of nitrogen to
watersheds, this is problematic (Robertson and Saad , 2011; Valigura et al., 2001; Boyer et al.,
2002). Additionally, both SWAT and GWLF uniformly distribute sources in a sub-catchment
unit within a basin.
The most geographically extensive statistical model that has been applied for the Major
River Basins (MRBs) in the United States is SPARROW (Robertson and Saad , 2011), which
relates nutrient loading in catchments and landscape factors to observed nutrient data via a
statistical model that is process informed (Smith et al., 1997; Robertson and Saad , 2011; Preston and Brakebill , 1999; Moore et al., 2004; Schwarz et al., 2006). This particular statistical
approach allows for the explicit description of landscape factors and a better understanding
of the magnitude of processing, but as applied, SPARROW lumps some nutrient sources into
large categories such as “urban” or “agriculture” which include several potential sources of
nutrients that are not uniformly distributed. In some cases a “surrogate”, such as population or urban area is used to account for the loading from a particular source or land use.
SPARROW only considers a single pathway for a particular basin, lumping overland flow
and groundwater pathways into a single parameter, and similar to SWAT and GWLF, the
sources are distributed evenly, throughout a sub-catchment (Schwarz et al., 2006; Robertson
and Saad , 2011).
With all of the approaches described above, there is a general tradeoff between the
scale of the model and the level of detail regarding the individual sources of nutrients and
their pathways. Another drawback of existing models is a lack of direct source accounting,
which makes it difficult to apply these models consistently in highly variable watersheds
(Nikolaidis et al., 1998; Destouni et al., 2006). Furthermore, with the exception of the
process-based models, existing approaches are generally applied to yearly loading averages
and do not consider seasonal variability in processes and pathways. Figure 2.1 is a qualitative
comparison of existing models illustrating these compromises.

41

Source Completeness/Resolution

This work
~250,000+ km 2
Budget Export
Coefficient

~150,000 km 2

SWAT/
GWLF

~50,000 km 2
~5000 km2

SPARROW
~400 km 2
Insight to Pathways and Processes
Figure 2.1: Qualitative comparison of nutrient modeling approaches.

2.2
2.2.1

Methods
Model Description

The model is a physically informed statistical model that uses a detailed and spatially explicit
description of nutrient sources. The form of the model equation, similar to SPARROW, considers a simplified conceptual structure in which sources applied to the landscape experience
attenuation along two pathways: travel within the upland portion of the basin and travel
within the stream. This attenuation is described by statistically derived reduction factors
that are functions of a basin’s physical features.
At the foundation of the modeling approach is a spatially detailed description of nutrient
sources. Nutrient applications are defined for each cell of the model using an estimation
process described in CHAPTER 1. The load at an observation point is modeled by summing
the contribution of nutrients from all cells that are upstream of the point and applying the

42

statistically derived reduction factors that are functions of a cell’s position in the watershed.
The functional form of the model is shown in Equation 2.1.
sources cells
Qij Rj F grdj Bgrdij + 1 − F grdj Bsurfj
i

(2.1)

j

Qij is loading of source i to cell j, Fgrd is the fraction of nutrients entering the subsurface pathway, and Bgrd, Bsurf, and R are reduction factors that account for physical and
biological attenuation in the watershed.
Bsurf and Bgrd describe attenuation that occurs as nutrients travel through the landscape to surface water; Bsurf represents the surface pathway and Bgrd represents the groundwater pathway. Bsurf and Bgrd are exponential functions of travel distance from cell j to
surface water. These factors are calculated as

(α ∗D )
Bsurfj = e 1 j

(2.2)

(α ∗D )
Bgrdj = e 2 j

(2.3)

where α1 and α2 are empirically derived coefficients and D is flow distance from cell j to
the nearest stream. The values of α1 and α2 are constrained between -1 and 0.
R describes attenuation that occurs during the in-stream portion of the pathway and is
an exponential function of both travel time and basin yield. The equation for R is shown in
Equation 2.4.

(α ∗T )
Rj = e 3 j ∗ α4 ∗ BYj

(2.4)

T represents the in-stream travel time from cell j to the downstream observation point,
α3 is an empirically derived coefficient constrained between -1 and 0, α4 is an empirically
derived coefficient constrained between 0 and 1, and BYj is the basin yield defined for each
sub-watershed. Basin yield, which is the total discharge of a basin divided by the basin area,
is a surrogate for recharge. Rivers with slow moving water and/or small, shallow channels,

43

may experience more biologic and physical processing due to increased interaction between
the channel bed, the hyporheic zone, and the biota in the stream (Alexander et al., 2000;
Mainston and Parr , 2002). These watersheds tend to have low basin yields.
One alternative groundwater pathway (septic plumes) and one alternative surface pathway (tile drains) were also considered in the model. Since attenuation in septic plumes occurs
differently than attenuation in other areas (Valiela et al., 1997; Robertson and Cherry, 1992;
Reneau and Pettry, 1976; Gilliom and Patmont, 1983), Bsep, which describes septic plume
attenuation, is substituted for Bgrd for the septic source. Bsep has the same form as Bsurf
(α ∗D )
and Bgrd (Bsep = e 7 j ). In the Lower Peninsula model, the coefficient for Bsep was

fixed based on analysis done by Valiela et al. (1997) who compiled data from several studies sampling septic plumes and concluded that roughly 35% of loss occurred in the septic
plume over a 200 meter distance with an exponential loss rate. Tiles alter the overland flow
(α ∗D )
pathway, so Btile (Btile = e 5 j ) was substituted for Bsurf in cells where tile drains

exist.
Fgrd is a function of cell recharge and is calculated as

F grd = α6 ∗ rechargej

(2.5)

Since septic tank loading occurs in the subsurface, Fgrd is set to 1 for the septic source
and the value of α6 is constrained between 0 and 1. Fgrd is set to 0 for all surface water
cells since atmospheric deposition falling directly on surface water experiences no landscape
attenuation. Similarly, Fgrd is set to 1 for point sources and point sources are applied only
to river cells, which leads to a basin travel distance Dj of 0 and a Bsurfj of 1. This ensures
that point sources only experience attenuation along the river pathway. The conceptual
model underlying the model equation is shown in Figure 2.2.
Several alternative parameters were considered in the modeling structure including travel
through lacustrine and palustrine wetlands, in-stream travel distance, and depth to bedrock.
These parameters were ultimately abandoned because models were relatively insensitive to

44

QchemAg + Qmanure + Qatm + QnonAg

Fgrd

Qpoint
Bsurf

Qseptic
Bsep

Btile
Bgrd

R

Figure 2.2: Conceptual model illustrating the sources and pathways described in the
model. Nutrient sources include septic tanks (Qseptic), agricultural chemical fertilizer
(QchemAg), manure (Qmanure), atmospheric deposition (Qatm), non-agricultural fertilizer
(QnonAg), and point sources (Qpoint). R, Bgrd, Btile, Bsep, and Bsurf are the reduction
factors and describe loss along the in-stream, groundwater, tile, septic, and surface runoff
pathways, respectively.

them, and they did not improve model performance.
Unlike SPARROW, the model has a spatially explicit description of sources and basin
factors, which includes within basin variability and the addition of travel distance and travel
time as factors. This allows for an expanded process-based description of pathways. The
differences between the spatially explicit model structure described here and a sub-basin
based model such as SPARROW are illustrated in Figure 2.3.

2.2.2

Landscape and Pathway Descriptions

In the model landscape characteristics and travel pathways are used to predict the attenuation of nutrients in the watershed. Similar to nutrient sources, each landscape factor is
45

Sub-basin Model (SPARROW)

Spatially Explicit Model (this work)

Qi,j
Bi,j

Ri
R B ∑Q

Qt = Total Load
Qi,j = Nutrients
Bi,j = Basin Reduction
Ri = River Reduction

Qt

Qt

cell
Qt = ∑(Qi,j*Rj*Bj)
j

Qt = ∑Q*R*B

Figure 2.3: A comparison of a sub-basin based approach (such as the USGS SPARROW
model) and the spatially explicit approach described here.

spatially explicit. The following sections describe the methods used to estimate these factors.

2.2.2.1

Watersheds and Travel Distance

Travel distance was calculated using the National Elevation Dataset (NED) (Gesch et al.,
2002; Gesch, 2007) and the ArcGIS Hydrology Toolbox, which calculates flow travel distance
based on an elevation dataset. Sub-watersheds were generated using the ArcGIS Hydrology
Toolbox. The points describing the outlet of the sub-watersheds coincide with sampling
locations.

2.2.2.2

Recharge

Recharge was estimated using a physically based statistical model. The model, developed
by Hyndman et al. (2007); Kendall and Hyndman (in review) for the Muskegon River wa-

46

tershed estimates the percentage of precipitation that becomes recharge from soil hydraulic
conductivity values and land class. The Muskegon River watershed model was applied to
the Lower Peninsula using soil hydraulic conductivity derived from the soil characteristics of
the Soil Survey Geographic Database (SSURGO) (Natural Resource Conservation Service)
and land class data from the National Land Cover Database (NLCD) (Fry et al., 2011).
Estimates of annual precipitation were downloaded from the PRISM climate group (PRISM
Climate Group, 2011).

2.2.2.3

Tile Drainage

The tile drained area of Michigan was surveyed by the the NRCS during the 1992 National
Resource Inventory (NRCS , 1995). These estimates were made on a county level and are not
up to date. In order to derive an updated tile drainage layer, a 750 meter buffer was added
to the Canal/Ditch feature of the National Hydrography Dataset (Roth and Dewald , 1999),
and the NLCD classified agricultural cells were selected within each buffer. The Canal/Ditch
feature includes features that do not drain tiles, and therefore, further classification was
required. Within the canal ditch buffers, areas of low slope (≤ 2 meter change in 30 meters)
and/or low soil hydraulic conductivity (≤5.0 x 10−6 cm/s) were selected as potentially tile
drained areas. This dual attribute classification scheme ensured that areas that are extremely
flat but with moderate drainage, based on SSURGO soil classification, such as those that
exist in Saginaw Bay, were classified as tile drained areas. The resulting tile drainage map
is shown in Figure 2.4.

2.2.2.4

Stream Velocity and Travel Time

Stream morphology plays a role in nutrient cycling by controlling the residence time of water
in a watershed and the contact between water and the stream bed (Mainston and Parr ,
2002; Alexander et al., 2000). The model accounts for the effects of stream morphology
using estimates of stream velocity, in-stream travel time, and basin yield.
47

Figure 2.4: Estimated tile drained areas for the Lower Peninsula of Michigan (indicated in
black).

Flow velocity was estimated using the relationship derived by Leopold and Maddock
(1953) which relates channel area and depth to discharge:

A = aQb

(2.6)

where A is channel area, Q is discharge and a and b are empirically derived coefficients. A
similar approach has been used to assess stream morphology and discharge in several studies
(Alexander et al., 2000; Schulze and Hunger , 2005; Bjerklie et al., 2003). The relationship for
channel area shown in Equation 2.6 was divided by Q and manipulated to estimate velocity
from channel area and flow:
1
v = Q1−b
a

(2.7)

The coefficients for these relationships were derived empirically from channel area and discharge data available from USGS gauges for Michigan’s Lower Peninsula (USGS , 2012). The
site visit observations of channel area and discharge were averaged for each gauge.

48

Site visit data was filtered for observations collected during high flow conditions (discharge observation between the 75th and 90th percentile) and low flow (baseflow) conditions
(discharge observations less than the 20th percentile). The 75th and 90th percentile were
chosen for high flow so that the dataset did not include observations collected under flood
conditions. The datasets for the two flow conditions were further subset based on the geology
underlying the gauge location; three geologic models were fit for the two flow conditions: till,
lacustrine, and outwash (Figure 2.5).
In order to apply these relationships to all stream cells, the discharge at every location in
the Lower Peninsula was estimated from the flow accumulation derived using NED and the
ArcGIS Hydrology Toolbox. USGS gauge data was used to derive empirical relationships
between calculated cell flow accumulation and discharge. The 80th and 20th percentiles of
discharge were used to create a high flow and low flow (baseflow) linear model. For the
baseflow relationship, flow accumulation was multiplied by baseflow index. Baseflow index
is the proportion of flow attributed to groundwater discharge and is estimated for gauges by
the USGS (Wolock , 2003).
The in-stream velocity was estimated from the predicted discharge in each cell using
Equation 2.7, and was then used to calculate travel time based on the flow length from a
cell to the observation point. Calculated travel time to the Saginaw Bay watershed sample
location is shown in Figure 2.6

2.2.3

Calibration/Validation Data

Flow measurements and nutrient samples were collected from the fall of 2010 to the spring
of 2012 in 5 field campaigns intended to capture synoptic conditions during baseflow, melt,
and spring rain. Sample locations for the 2010 baseflow, 2011 melt, and 2011 spring rain
were chosen in an effort to capture most of the surface water outflow to the Great Lakes
from the Lower Peninsula. Sites were located at the farthest accessible downstream sample
point on most major and several minor streams and rivers discharging to the Great Lakes.

49

High Flow Till Model

1000
800
600
400
y = 1.5644x0.8441
R² = 0.9222

200
0
0

1000

Channel Area ft3

Channel Area ft3

y = 1.515x 0.8158
R² = 0.9419
5000

y = 2.7752x0.8648
R² = 0.9347

300
200
100
0

1000

0

400

0

1500

0

500

2000

High Flow Outwash Model
2000

500

Low Flow Till Model

600
Channel Area ft3

Channel Area ft3

1200

800
600
400
y = 2.5927x 0.7621
R² = 0.9494

200
0
0

1000

2000

Low Flow Lacustrine Model
Channel Area ft3

Channel Area ft3

300

1000

10000

y = 4.038x 0.7272
R² = 0.9471
0

200

Low Flow Outwash Model
1400
1200

High Flow Lacustrine Model
Model
900
800
700
600
500
400
300
200
100
0

100

500
1000 1500
Discharge (cfs)

900
800
700
600
500
400
300
200
100
0

y = 4.038x0.7272
R² = 0.9471
0

500
1000 1500
Discharge (cfs)

Figure 2.5: Graphs showing the relationship between discharge and channel area.
Regressions were fit to USGS channel area and discharge data for observations between the
75th and 90th percentile for discharge (high flow) and observations less than the 20th
percentile (low flow). Gauge data was subset by geology type (lacustrine, outwash, and till).

50

Melt In-Stream
Travel Time to
Sample Point
50
20
10
5

Figure 2.6: In-stream travel time to sample point for the Saginaw Bay watershed during
snow melt.
A total of 67 locations were sampled.
The 2011 baseflow and 2012 melt sampling focused on 4 major watersheds in the Lower
Peninsula of Michigan: Grand River in south central Michigan, Saginaw Bay in eastern
Michigan, Muskegon River in north central Michigan, and Boardman-Charlevoix in northwest Michigan (Figure 2.7). Between 15 and 30 sites were sampled within each watershed for
a total of 90 sites. Sites were chosen in an effort to evenly distribute samples throughout the
watershed. This was accomplished by selecting points so that the area upstream of each sample point would have, approximately, the same incremental watershed area. The watersheds
for the combined sampling covers over 70% of the land area of the Lower Peninsula.
The watersheds sampled represent a large range in watershed size, land use, and nutrient
inputs. The samples also represent a large range of flow conditions. The datasets collected
from June 2011- March 2012 were used to calibrate the models and October 2010-March
2011 datasets were used for validation. The methods used for measuring flow and analyzing nutrient concentrations are described in CHAPTER 1 and sample data is included in
Appendix A.

51

Boardman-Charlevoix
Muskegon River
Saginaw Bay
Grand River

Figure 2.7: Location of Boardman-Charlevoix, Muskegon River, Saginaw Bay, and Grand
River watersheds.
2.2.4

Optimization

The June 2010-March 2012 dataset was used as a training dataset for the optimization of
the model. As discussed in Section 2.2.3 this dataset includes data from the Grand River,
Muskegon River, Saginaw Bay, and Boardman-Charlevoix watersheds (Figure 2.7). The
values of the coefficients for each model run were fit using the Matlab function fminseach,
which uses the Nelder-Mead simplex algorithm (Nelder and Mead , 1965), a direct search
algorithm. A total of 6 coefficients were optimized (Bsurf, Bgrd, Btile, R, Fgrd, and BY ).
The target of the optimization was the mean absolute difference between the natural logarithm of the observed area normalized load and the natural logarithm of the modeled area
loadmodel
loadobserved
normalized load (mean abs ln( area ) − ln(
) ). Bounds were applied
area
to the coefficients to reduce optimization time and ensure that the coefficients applied were in
agreement with the process-based form of the statistical model. For instance, the coefficient
that modifies distance in the exponential portion of basin reduction was limited between 0
and -1 since the basin reduction describes the distance based attenuation of nutrients. While
in reality, depending on the time of year and the biologic and physical conditions, the basin

52

could become a source of stored nutrients, it was assumed that the net effect of these processes would be small over the 5 year period described by the model. The reduction factors
describe permanent removal of nutrients relative to the time scale of the model.

2.3

Results

2.3.1

Sample Results

The results of the sampling showed that total nitrogen (TN) and total phosphorus (TP) concentrations decreased from southeast to northwest (Figure 2.8). Concentrations of nitrogen
ranged from 64 to 5583 µg/L during baseflow and 179 to 7886 µg/L during melt. Phosphorus
ranged from below detection limit to 476 µg/L for baseflow and from below detection limit
to 395 µg/L for melt. The highest concentrations for nitrogen were observed during melt
except for a few watersheds in the northeastern part of the state which had higher concentrations during spring rain of both TN and TP. For phosphorus, some watersheds showed
higher melt concentrations while others showed higher baseflow concentrations.

2.3.2
2.3.2.1

Model Results
Model Performance

The results of the 6 models (Baseflow Total Nitrogen, Melt Total Nitrogen, Spring Rain
Total Nitrogen, Baseflow Total Phosphorus, Melt Total Phosphorus, and Spring Rain Total
Phosphorus) are shown in Tables 2.1 and 2.2. All models had adjusted R2 values greater
than 0.64 and all but the Melt Total Phosphorus model had an adjusted R2 greater than
0.80. Overall, the nitrogen models performed the best. All nitrogen models had R2 values
greater than 0.86 and the baseflow model had the highest R2 value of 0.95. Seasonally, the
baseflow models were the most successful and the melt models were the least successful.
However, within the validation sets, the TN melt model performed the best overall.

53

TN Base Calibration

TN Melt Calibration

TN Rain Calibration

µg/L
2000
1500
1000
500

TP Base Calibration

TP Melt Calibration

TP Rain Calibration

µg/L
100
75
50
25

TN Base Validation

TN Melt Validation

µg/L
2000
1500
1000
500

TP Base Validation

TP Melt Validation

µg/L
100
75
50
25

Figure 2.8: Sampled watersheds with concentration (µg/L) indicated for each sampling
campaign.

54

Adj
Model

R2

TN Base
TP Base
TN Melt
TP Melt
TN Rain
TP Rain

0.95
0.87
0.85
0.69
0.88
0.82

Validation

Optimization

Validation

R2

R2

n

n

0.95
0.86
0.84
0.64
0.87
0.80

0.86
0.81
0.90
0.77
NA
NA

89
89
88
86
62
62

61
61
64
59
NA
NA

Error
0.42
0.52
0.45
0.67
0.47
0.54

Table 2.1: Summary of estimated model parameters. R2 values are for estimated and
observed loads.
Model
TN Base
TP Base
TN Melt
TP Melt
TN Rain
TP Rain

Bsurf
α1
-4.03E-03
-2.14E-03
-6.26E-05
-4.64E-10
-1.18E-03
-4.27E-04

Btile
α5
-1.28E-01
-9.90E-01
-3.98E-10
-9.78E-03
-2.87E-02
-2.14E-03

Bgrd
α2
-4.59E-02
-8.97E-01
-7.04E-05
-9.39E-08
-3.81E-11
-4.26E-02

Fgrd
α6
5.44E-01
9.99E-01
8.34E-05
4.24E-03
1.00E+00
9.99E-01

R
α3
-1.87E-10
-9.04E-07
-1.15E-10
-3.62E-10
-8.69E-08
-1.15E-08

BY
α4
9.93E-01
1.00E+00
1.00E+00
1.00E+00
1.00E+00
5.41E-01

Table 2.2: Summary of coefficients associated with each basin parameter.

With the exception of the rain event data, the calibration dataset consists only of data
collected from the Grand River, Saginaw Bay, Boardman-Charlevoix, and Muskegon River
watersheds while the validation dataset contains samples covering nearly the entire Lower
Peninsula. The validation dataset has high R2 values and closely follows the 1:1 line, demonstrating the geographic and temporal transportability of the model. The rain models do not
have a validation dataset because the spring rain was only sampled in 2011.
The R2 values shown in Table 2.1 are for predicted and observed loads. The target of the
model optimization was area normalized loads. Comparisons of model results for both loads
and area normalized loads are shown in Figures 2.9 and 2.10. Plots of the area normalized
loads indicate how consistently the models predict loads across watersheds irrespective of
watershed size. The baseflow models best predict area normalized loads. The melt and rain
models show a greater degree of variability in predicted area normalized loads for watersheds

55

with relatively small per unit area loads.

2.3.2.2

Sensitivity

The sensitivity for each parameter is shown in Table 2.3. Sensitivity was calculated as

Sensitivity =

abs modeloptimized − modeloptimized ∗ ∆

∗ 100

(2.8)

modeloptimized
where

is the change in the function value (0.5%) and modeloptimized is the optimized

value of the coefficient. These values were normalized to the maximum sensitivity value to
indicate the relative sensitivity of each parameter. The results show that the models were
generally sensitive to Bsurf and BY. These two factors account for the basin attenuation and
in-stream attenuation respectively. The models are also generally sensitive to Fgrd which
predicts the portion of the source that travels via the groundwater pathway. The results
suggest that attenuation is relatively most important along the surface pathway and that,
with the exception of TP during baseflow, in-stream processing is best predicted by basin
yield. This agrees with the analysis by Alexander et al. (2000) who showed that in-stream
loss was correlated with channel depth (approximated from flow rate). They suggest that
depth was a surrogate for the ratio of water volume to stream benthic area. Basin yield can
also be a surrogate for this.
The TP Melt model only shows sensitivity to basin yield and the TN Melt model only
shows sensitivity to basin yield and Bsurf. The model predicts that almost all of the nutrients during melt are delivered via the surface pathway. Therefore, most of the nutrients
are only impacted by surface processing. Additionally, the TP Melt model is the least successful of all the models. The model’s low sensitivity to all but one parameter suggests
that the model parameters do not adequately describe the dynamics of phosphorus during
melt. The delivery of sediment is not directly described by the model but is an important
transport mechanism for phosphorus (Mainston and Parr , 2002). Relatively large amounts

56

-5
-6
-7
-8
-9

-10

5
4
3

-6
-5
2 /day)
Log Observed (kg/m

-9

-8

-7

2

TN Melt Area Normalized Load
-4
y = 0.6976x - 1.6899
R2 = 0.6903
-5

-6

3
4
5
6
Log Observed (kg/day)
TN Melt Load

5
Log Model (kg/day)

Log Model (kg/m2 /day)

y = 0.966x + 0.1599
R2 = 0.9531

2
-10

y = 0.8411x + 0.5021
R 2 = 0.8527

4
3
2
1

-7

1

-7
-6
-5
-4
2 /day)
Log Observed (kg/m
TN Rain Area Normalized Load
-4
y = 0.595x - 2.3644
R2 = 0.4086
-5
-6
-7
-8

2
3
4
5
Log Observed (kg/day)
TN Rain Load

5
Log Model (kg/day)

Log Model (kg/m2 /day)

TN Base Load

6
Log Model (kg/day)

Log Model (kg/m2 /day)

TN Base Area Normalized
Load
y = 0.8218x - 1.2041
R2 = 0.7607

y = 0.9092x + 0.2637
R 2 = 0.8798

4
3
2
1
0

-8
-7
-6
-5
-4
Log Observed (kg/m2 /day)

0

1
2
3
4
5
Log Observed (kg/day)

Figure 2.9: Graphs of observed and modeled TN loads and area normalized loads for the
calibration dataset.

57

-6

y = 0.9134x - 0.7219
R 2 = 0.7191

-7
-8
-9

1
0
-1
-2

-9

-3

-8

-7
-6
Log Observed (kg/m2 /day)
TP Melt Area Normalized Load
-5
y = 0.6108x - 2.6528
R 2 = 0.3055

4
Log Model (kg/day)

-10

Log Model (kg/m2 /day)

y = 0.9834x - 0.0283
R2 = 0.8741

2

-3

-10

-6

-7

-8

3

-2 -1 0
1
2
3
Log Observed (kg/day)
TP Melt Load

y = 0.8604x + 0.2265
R 2 = 0.685

2
1
0

-8
-7
-6
-5
Log Observed (kg/m2 /day)
-5
-6

0

TP Rain Area Normalized
Load
y = 0.8815x - 0.8409
R2 = 0.4153

-7
-8
-9

1
2
3
4
Log Observed (kg/day)
TP Rain Load
y = 0.9348x + 0.1039
R 2 = 0.8206

4
Log Model (kg/day)

Log Model (kg/m2 /day)

TP Base Load

3
Log Model (kg/day)

Log Model (kg/m2 /day)

TP Base Area Normalized
Load

3
2
1
0
-1
-2

-9

-8

-7

-6
-5
Log Observed (kg/m2 /day)

-1

0
1
2
3
4
Log Observed (kg/day)

Figure 2.10: Graphs of observed and modeled TP loads and area normalized loads for the
calibration dataset.

58

5
4
3
2
1
0
-1
-2
-3

TP Base Load
3

y = 0.8622x - 1.5526
R2 = 0.8602

Log Model kg/day

Log Model (kg/day)

TN Base Load

y = 0.8945x - 0.3311
R2 = 0.8138

2
1
0
-1
-2
-3

-3 -2 -1 0 1 2 3 4 5
Log Observed (kg/day)

-3

-2 -1 0
1
2
3
Log Observed (kg/day)

TN Melt Load
4

y = 0.8861x + 0.22
R 2 = 0.9024

4

Log Model (kg/day)

Log Model (kg/day)

5

TP Melt Load

3
2
1
0

3

y = 0.857x + 0.2501
R 2 = 0.773

2
1
0
-1

0

1
2
3
4
Log Observed kg/day

5

-1

0
1
2
3
4
Log Observed (kg/day)

Figure 2.11: Graphs of observed and modeled loads for the validation dataset.

59

Model
TN Base
TP Base
TN Melt
TP Melt
TN Rain
TP Rain

Bsurf
1.000
1.000
0.149
0.000
0.615
1.000

Btile
0.051
0.011
0.000
0.002
0.119
0.009

Bgrd
0.007
0.003
0.000
0.000
0.004
0.000

Fgrd
0.134
0.182
0.000
0.000
0.083
0.450

R
0.000
0.351
0.000
0.000
0.006
0.027

BY
0.674
0.313
1.000
1.000
1.000
0.001

Table 2.3: Normalized sensitivity of model to each basin parameter coefficient. Sensitivity
is calculated as the percent change in the model prediction per percent change in the
parameter.

of sediments are mobilized during melt, and therefore, sediment delivery may be a sensitive
parameter that is not included in the model.

2.3.3

Sources

A comparison of the total nutrient sources (as applied) relative to the total nutrients delivered downstream illustrates the importance of basin and river processing for attenuation of
nutrient loads, the seasonal variability in pathways, and differences in the delivery of phosphorus and nitrogen. The results show that atmospheric deposition and commercial chemical
fertilizer are the largest contributors to the in-stream nitrogen load. During melt, commercial
fertilizer is the largest source of nitrogen while during baseflow and rain events, atmospheric
deposition is the largest source of nitrogen. Point sources and direct atmospheric deposition
become a relatively larger contributor to the total load when there is more basin processing (such as during baseflow). The model predicts that there is very little basin processing
during the melt event so surface derived sources are relatively more important during melt.
The largest contributor of phosphorus to surface water based on this research is commercial chemical fertilizer followed by manure. Atmospheric deposition contributes a relatively
small amount of phosphorus and is often considered effectively zero in many nutrient loading
models. During baseflow, the model predicts that nearly 3% of total phosphorus delivered
to streams is derived from atmospheric deposition, so ignoring atmospheric deposition as a

60

Model
TN Base
TP Base
TN Melt
TP Melt
TN Rain
TP Rain

Bsurf
(% 1 km)
2
12
94
100
31
65

Btile
(% 0.1 km)
0
0
100
38
6
81

Bgrd
(% 0.1 km)
1
0
99
100
100
1

Fgrd max
(%)
54
100
0
0
100
100

Rmin
(%/day)
5
5
4
4
4
2

Rmax
(%/km)
99
92
100
100
99
54

Table 2.4: Summary of calculated basin parameters. Reduction factor units (Bsurf, Btile,
and Bgrd ) are the percent of the nutrient load that remains after the first .1 kilometer (1
km for Bsurf ) of travel. Units of Rmin and Rmax are the percent of the load that remains
after the first day of in-stream travel for the sub-watershed with the highest relative basin
yield (Rmax ) and the lowest relative basin yield (Rmin). The bold values indicate the
coefficients with relatively high model sensitivity.

100%
80%
Point
Septic

60%

NonAg
Manure

40%

ChemAg
20%

Atm

0%
TN Base TN Melt TN Rain TP Base TP Melt TP Rain
Figure 2.12: Source contributions to total loading to streams.

source of phosphorus may lead to up to 3% under-prediction of total phosphorus in these
models. However, in watersheds that contain large amounts of unmanaged land (forest,
range, grasslands), this percentage can be much larger. Figure 2.14 shows that atmospheric
deposition contributes over 10% of the observed TP load in several sampled watersheds and
in one sampled watershed, contributes up to 25% to the observed TP load.
The source contribution to total downstream load predicted by this model agrees well

61

with predictions from the SPARROW model constructed for the Great Lakes (Robertson
and Saad , 2011). However, a few sources identified by the SPARROW model are surrogates
for other sources. For instance, SPARROW attributes 15% of the total phosphorus load to
an “urban” source. This model predicts that approximately 4-7% (depending on season)
of total phosphorus is contributed by non-agricultural fertilizer and septic tanks and 3% of
phosphorus is contributed by atmospheric deposition. Additionally, 4-24% is derived from
point sources. Depending on the nature of the urban area being described, increases in
urban area may not lead to proportional increases in the urban sources. For instance, low
to medium density urban areas like Lansing may have a higher incidence of septic and lawn
fertilizer use than a high density urban area such as downtown Detroit. The fact that these
sources are not explicitly described in the SPARROW model may account for the urban bias
observed in the SPARROW model (Robertson and Saad , 2011; Alexander et al., 2002)

2.3.4
2.3.4.1

Pathways and Processes
Basin Processing

The parameter Fgrd describes the portion of the applied load entering the groundwater pathway. The model predicts that in basins with the largest recharge, nearly all nutrients enter
the groundwater pathway during baseflow and rain while during melt, nearly all nutrients
travel by the surface, runoff pathway and almost nothing enters the groundwater pathway.
It is important to remember that this is for basins with the largest recharge, and therefore,
this value varies between 0 and the maximum predicted value depending on the estimated
recharge for a cell.
The groundwater basin reduction factors show that groundwater attenuation occurs over
a much shorter distance than attenuation along the surface pathway. Almost all of the
attenuation occurs in the first 0.1 kilometer of travel in the baseflow model. This suggests
that the groundwater pathway only contributes nutrients to surface water if the nutrient

62

Atmospheric Deposition

Manure

Agricultural Chemical Fertilizer

% Delivered

% Delivered

50
25
15
05

% Delivered

50
30
20
10

80
60
40
20

Septic Tanks

Non-Agricultural Fertilizer
% Delivered

% Delivered

% Delivered

10
05
01
0.7

02
01
0.5
0.2

Point Sources

10
02
01
.50

Figure 2.13: Percent of observed TN melt load derived from each source.

Base
8.4

TN
Melt
0.4

Rain
58.6

Base
7.8

TP
Melt
1.1

Rain
3.4

Table 2.5: Percent of total load delivered to surface water via the groundwater pathway.

source is within a travel distance of 0.1 km. The TN Melt and TN Rain models predict that
most of the nutrients persist in the groundwater beyond 0.1 kilometer. The proportion of
nutrients delivered to streams via the groundwater pathway during each season are shown in
Table 2.5. The models predict that relatively more nitrogen is delivered via the groundwater
pathway relative to phosphorus. The model was relatively insensitive to the Bgrd parameter,
so these predictions are potentially highly variable.

63

Atmospheric Deposition

Agricultural Chemical Fertilizer

25
10
04
02

Septic Tanks
% Delivered

% Delivered

% Delivered

50
40
30
15

% Delivered

07
05
02
01

Manure

70
60
50
30

Non-Agricultural Fertilizer
% Delivered

Point Sources

% Delivered

15
10
05
02

07
03
02
01

Figure 2.14: Percent of observed TP melt load derived from each source.

The results for the basin reduction factors indicate the differences in processes and pathways both seasonally and between nitrogen and phosphorus. The results show that phosphorus is more efficiently delivered to streams than nitrogen for all models. The largest
attenuation along the surface pathway occurs during baseflow. This could be a result of
lower temperatures or higher runoff rates.
The portion of each cell’s applied load delivered to surface water is shown in Figure 2.15,
which indicates that landscape retention of phosphorus and nitrogen during melt is different.
Nitrogen derived in areas that are distant to surface water are retained while phosphorus
is retained in areas with tile drainage. This is a result of Btile which indicates relatively
more phosphorus is retained in tile drains during melt than nitrogen. This effect is reversed

64

TN Base

TP Base

TN Melt

TN Rain

TP Melt

TP Rain

% Exported to
Surface Water
100
(Surface Water)
90
75
50
25

Figure 2.15: Percent of total applied load delivered to surface water.

during the spring rain event and no nutrients are delivered via tile during baseflow. Tile
drains affect flow pathways and processes in several ways. Tile drains change the natural
drainage pathway of water, decreasing the travel time to surface water and decreasing the
amount of time available for physical and chemical processes that remove nutrients from
agricultural runoff (Robertson and Saad , 2011). Several different physical and biological
processes that affect nutrients differently may be affecting the delivery of nutrients in tile
drains.

65

2.3.4.2

In-stream Processing

The results for in-stream processing show that in basins with the largest basin yield (Rmax),
very little in-stream attenuation occurs. This is consistent with other studies that have shown
that a higher rate of attenuation occurs in small, slow moving streams where there is more
time for water to interact with the stream bed (Robertson and Saad , 2011; Mainston and
Parr , 2002; Alexander et al., 2000). The model also predicts that slightly more attenuation
occurs during baseflow when discharge is lower.

2.3.5

Land use Export Rates

Table 2.6 shows the average export rate of nutrients from agriculture, urban, and unmanaged
(forest, grasslands, shrublands) and pasture cells. Several studies have attempted to estimate
annual nutrient exports rates as a means of quantifying and predicting non-point source
contribution of nutrients. Beaulac and Reckhow (1982) compiled several studies estimating
export rates. They found values for phosphorus ranging from 0 to 1 kg/ha/year for forest,
∼0 to 18 kg/ha/year for agriculture, and ∼0 to 6 kg/ha/year for urban, and values for
nitrogen from ∼0 to 6 kg/ha/year for forest, ∼0 to 80 kg/ha/year for agriculture, and ∼0
to 40 kg/ha/year for urban. The seasonal values predicted by this model are similar to
these values, but slightly higher rates are predicted for nitrogen export from forested and
agricultural land.
There is a high degree of seasonal variability in the estimates of nutrients exported from
each land type. The export rate for urban land appears to be higher than cropland during
baseflow for nitrogen. This could be primarily due to the relative increase in point source
delivery during baseflow (Figure 2.16). The melt export rates are 2 to 15 times greater than
baseflow rates for total nitrogen and 0 to 9 times higher for phosphorus.

66

TN Base

TP Base

100%

100%

80%

80%

60%

60%

40%

40%

20%

20%

0%

0%
TN Melt

TP Melt

100%

100%

80%

80%

Septic

60%

60%

NonAg

40%

40%

Manure

20%

20%

0%

0%

TN Rain

Point

ChemAg
Atm

TP Rain

40%

20%

20%

0%

0%
Pasture

Pasture

40%

Unmanaaged

60%

Urban

60%

Cropland

80%

Unmanaaged

80%

Urban

100%

Cropland

100%

Figure 2.16: Source contributions to total stream loading within different land use classes.

67

Cropland
Urban
Unmanaged
Pasture

Base
6
10
1
3

TN
Melt
92
27
13
41

Rain
35
17
7
18

Base
1
1
0
1

TP
Melt
9
2
0
4

Rain
5
1
0
2

Table 2.6: Estimated nutrient exports from each land use in kg/ha/year.

2.3.6

Residuals

For the residual analysis, error was calculated as

Error = log(loadmodel ) − log(loadobserved )

(2.9)

A residual of 1 indicates an order magnitude over prediction by the model. The model
residuals are plotted in Figure 2.17.
Figure 2.17 indicates that the medians of the model residuals are close to 0. Additionally,
with few exceptions, model residuals are less than an order of magnitude. In fact, only one
model, TN Rain, had an observation with a residual greater than 1. TP Rain had two
residuals of 0.96 and TN Rain has a single residual of 0.98. The spatial distribution of the
model residuals is shown in Figure 2.18.
Alexander et al. (2002) compared regional nutrient loading models in terms of percent
error where error is calculated as

Error =

S−O
∗ 100
O

(2.10)

Here, S is the simulated value and O is the observed value (Alexander et al., 2002). They
observed that if the model errors are less than the standard error of the observations, model
error is indistinguishable from zero. The median percent error, minimum percent error, and
maximum percent error (following Alexander (2002)) for each model are shown in Table
2.7 along with the average percent difference between samples with duplicates. The model
median percent error is less than the percent sample difference for all models except the TP
68

Melt model. The median percent errors were comparable or less than the median percent
error of models compared in Alexander et al. (2002) which ranged from -0.6 to 46.5 %.
The minimum and maximum percent error were also comparable or less than other model
residuals (Alexander et al., 2002).

2.3.6.1

Residuals and Land Use

To assess the model’s ability to perform in watersheds of different land uses, model residuals
were compared to watershed urban, agricultural, and forested area. Linear regressions were
fit to the percent area of each land use in the watershed and log model error.
The results of the regressions show that nitrogen models have no apparent land use
biases. All the regressions show that the residuals and the watershed land use have very
little correlation. Total phosphorus, however, shows agricultural bias. In all of the models,
loading was slightly overestimated in watersheds with a high proportion of agricultural land
and slightly underestimated in watersheds dominated by forest. The effect is strongest
during rain and melt seasons. This suggests that the model may not be accurately capturing
the effect of sediment bound phosphorus. Delivery via sediment is an important transport
mechanism for phosphorus as more sediment is mobilized during the melt and spring rain
seasons (Withers and Jarvie, 2008).

2.3.6.2

Model Bias

A regression of prediction errors and watershed characteristics was used to further explore
potential model bias. This approach was used by Alexander et al. (2002) who compared
several large basin statistical nitrogen loading models. The same watershed characteristics
that were used in their analysis were used to assess model bias here including basin area,
runoff (basin yield), cultivated area, and urban area. The regression equation is as follows:
E = β0 + β1 ∗ C + β2 ∗ U + β3 ∗ BA + β4 ∗ R

69

(2.11)

where C is the percent of the watershed that is cultivated, U is the percent of the watershed
that is urban land, BA is watershed area in m2 , R is runoff in m3 /year, and β0 , β1 , β2 , β3 ,
and β4 are the regression coefficients. E is the percent error (Equation 2.10).
This approach provides an indication of model bias and a comparison of the performance
of this model with other regional scale statistical nutrient loading models. The value of the
model slope indicates the percentage increase in model residual that results from an increase
in the basin factor. The p-value is the significance of the predicted slope; p-value indicates
the level of confidence in the influence of the parameter on the model residual. The results
of the analysis are shown in Table 2.8.
The results show that the models have very little bias. When compared to the other
nitrogen loading models studied by Alexander et al. (2002), these models show the smallest
bias (lowest R2 , slopes, and least amount of significance attributed to model parameters).
The only model that shows significant bias is the TP Base model. Again, this may be due to
an inability of the model to accurately account for nutrients delivered via sediment. Another
potential source of bias may be recycling and re-suspension of phosphorus during high melt
flows; this effect is also not accounted for by the model.

2.4

Discussion

This work demonstrates the value of a spatially explicit description of nutrient sources. Not
only does the model perform well in watersheds with diverse conditions and scales, it also
provides information beyond a prediction of nutrient load including seasonal dynamics, and
source apportionment. It performs as well as or better than other regional scale models and
has less land use bias than other regional nitrogen models (Alexander et al., 2002; Robertson
and Saad , 2011). The model performed consistently across diverse watersheds with little
to no land use bias (particularly the TN models) and performs consistently across three
different seasons. This suggests that this model has the potential to perform well in other
regions and for other time periods.
70

The model indicates that atmospheric deposition is the largest source of nitrogen to surface water in the Lower Peninsula and chemical agricultural fertilizer is the most significant
source of total phosphorus. Atmospheric deposition contributes, on average, 3% of phosphorus to the in-stream load and was as much as 25% in sampled watersheds during melt. Most
models do not account for atmospheric deposition of phosphorus and may be attributing
this phosphorus to another source, particularly in forested watersheds where atmospheric
deposition is a more significant source of phosphorus.
The results also indicate that nitrogen export rates can be more than 15 times larger
during melt than baseflow and 2.5 times larger during spring rain. For phosphorus, as much
as 9 times more phosphorus can be exported during melt. There is a significant amount of
variation in export rates between seasons and between land uses. This has implications for
management strategies. Resources for water quality should be allocated based on the anticipated improvement in water quality (Destouni et al., 2006). Depending on the time of year
and the dominant nutrient source in a watershed, the effectiveness of different management
strategies may vary. Spatially explicit nutrient models provide a method to estimate the
efficiency of different strategies and the potential areas and sources that should be targeted.
The results show that the groundwater pathway (depending on season) accounts for 0.459% of total nitrogen and 1-8% of total phosphorus delivered to the stream. The model
can be improved with an expanded description of the groundwater pathway. Currently the
model uses a surface travel distance to describe groundwater attenuation and uses a surface
watershed boundary to determine nutrient inputs to the groundwater pathway. A better
description of groundwater travel pathway may improve model predictions in watersheds
that have groundwatersheds that are significantly different from surface watersheds. The
addition of travel time in the groundwater pathway may also improve performance since
groundwater residence time can be highly variable.
The presented model cannot describe the effect of wetlands and impoundments which
have been shown to impact the ability of watersheds to retain nutrients (Robertson and

71

Saad , 2011). A spatially explicit description of these features should improve the model
performance and provide insight to to their impact on nutrient loading. Finally, the bias of
the phosphorus models needs be explored further; perhaps a better description of sediment
delivery is needed.

72

Log Load Error kg/year
(Log Model-Log Obs)

1.0
0.5
0.0
-0.5
-1.0
TN Melt TP Melt TN Base TP Base TN Rain TP Rain
Model
Figure 2.17: Box plots of model residuals.

TN Base
log10 Error
.30
.15
-.15
-.30

log10 Error
.30
.15
-.15
-.30

TN Melt
.30
.15
-.15
-.30

TP Base

log10 Error

TN Rain
log10 Error

log10 Error

.30
.15
-.15
-.30

TP Melt

.30
.15
-.15
-.30

log10 Error

TP Rain

.30
.15
-.15
-.30

Figure 2.18: Model residuals for sampled watersheds (Equation 2.9). Only incremental
watershed area is shown.

73

Median

Model
TN Base
TP Base
TN Melt
TP Melt
TN Rain
TP Rain

Mean

Min

Max

%
-1.78
-4.7
-7.4
-21.2
-0.42
-5.6E-05

%
25
11
14
19
46
60

%
-72.85
-87.02
-77.95
-84.10
-72.53
-93.82

%
416.45
291.10
566
672.66
1172.7
827.90

Sample
Diff.
%
-6
8
-10
18
-28
-18

Table 2.7: Model errors calculated as percent of observed loads. Sample Diff. is the mean
percent difference between samples and duplicates.

Model
TN Melt
TP Melt
TN Base
TP Base
TN Rain
TP Rain

Value
slope
p-value
slope
p-value
slope
p-value
slope
p-value
slope
p-value
slope
p-value

Cultivated
(%)
24.87
0.64
84.62
0.28
33.82
0.52
172.50
0.00
12.13
0.29
18.01
0.75

Urban
(%)
27.44
0.84
-264.10
0.17
-162.80
0.21
-170.20
0.08
14.93
0.23
3.73
0.95

Basin
(km2 )
0.00
0.45
0.00
0.81
-0.01
0.18
0.00
0.25
0.00
0.47
0.00
0.51

Runoff
(m2 /year)
-0.13
0.53
0.14
0.64
-0.48
0.45
-0.53
0.26
0.07
0.30
-0.39
0.25

Intercept

R2

18.72
0.47
5.37
0.89
42.81
0.07
-14.78
0.40
-5.53
0.34
38.84
0.18

0.01
0.90
0.04
0.54
0.05
0.31
0.26
0.00
0.06
0.49
0.04
0.65

Table 2.8: Results of the residual regressions. Slope indicates the percent change in the
model error per unit change in the parameter. Bold values indicate significant parameters
at a 95% confidence level.

74

TP Base

R 2 = 0.0157

0.6
0.4
0.2
0
-0.2
-0.4
-0.6

Log Load Error kg/year

Log Load Error (kg/year)

TN Base
0.8

-0.8

0.8

% Urban Area
TN Base
R 2 = 0.0203

0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8

0.5
% Forest Area
TN Base
R2 = 0.0121

0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
0

0.5
% Cultivated Area

0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1

0.5
% Forest Area
TP Base

0.5
% Cultivated Area

Figure 2.19: Baseflow model residuals and land use regression results.

75

0.5

1

R2 = 0.1606

0

1

% Urban Area
TP Base

R2 = 0.1897

0

0.8
0.6

0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1

1

Log Load Error (kg/year)

Log Load Error (kg/year)

0

R 2 = 0.0016

0

0.5

Log Load Error (kg/year)

Log Load Error (kg/year)

0

0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1

1

Log Load Error (kg/year)

0
1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8

% Urban Area
TN Melt

Log Load Error (kg/year)

R2 = 0.0171

0.5
0
-0.5

0

0.5

R 2 = 0.0066

1

% Urban Area
TP Melt

0.5

R 2 = 0.1649

0.5
0
-0.5
-1

0

0.5
% Forest Area

0

1

1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8

0.5
% Forest Area

1

TP Melt

TN Melt
R2= 0.0015

Log Load Error (kg/year)

Log Load Error (kg/year)

TP Melt

1

-1

Log Load Error (kg/year)

Log Load Error (kg/year)

TN Melt
R2 = 0.0001

1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8

1

R 2 = 0.0986

0.5
0
-0.5
-1

0

0.5
% Cultivated Area

0

1

0.5
% Cultivated Area

Figure 2.20: Melt model residuals and land use regression results.

76

1

0.8
0.4
0
-0.4
-0.8
-1.2
-1.6

1.2
1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8

0.5
% Urban Area
TN Rain

1

0
1.2
Log Load Error (kg/year)

Log Load Error (kg/year)

Log Load Error (kg/year)

R2 = 0.0135

0

R2 = 0.0465

0.5
% Urban Area
TP Rain

1

R2 = 0.1719

0.8
0.4
0
-0.4
-0.8
-1.2
-1.6

0

Log Load Error (kg/year)

TP Rain
R 2 = 0.0509

1.2

1.2
1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8

0.5
% Forest Area
TN Rain

1

0
1.2
Log Load Error (kg/year)

Log Load Error (kg/year)

TN Rain

1.2
1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8

R 2 = 0.0144

0.5
% Forest Area
TP Rain

1

R 2 = 0.0698

0.8
0.4
0
-0.4
-0.8
-1.2
-1.6

0

0.5
% Cultivated Area

1

0

0.5
% Cultivated Area

Figure 2.21: Rain model residuals and land use regression results.

77

1

CHAPTER 3
USING NITRATE ISOTOPE DATA TO VALIDATE A REGIONAL
NUTRIENT SOURCE MODEL

3.1

Introduction

With coastline on four of the five Great Lakes, Michigan has a significant influence on the
ecological health of the Great Lakes. Elevated levels of nitrate are a threat to both human
and ecosystem health (Davidson et al., 2012), and therefore, legislation in the Great Lakes
and elsewhere has attempted to limit and control sources of nitrate. Correlations between
nitrate and human disturbance and population have been observed in several major river
basins worldwide (Howarth et al., 1996; Vitousek et al., 1997; Hatfield and Keeney, 2008),
and a substantial amount of effort has focused on quantifying and understanding links between landscape processes and lake water quality. Of particular concern to environmental
managers, are non-point sources of nitrate, which are difficult to quantify (Dolan and McGunagle, 2005; PLUARG, 1983; Nikolaidis et al., 1998). Non-point sources of nitrate and
the retention processes that occur within watersheds are spatially variable which makes
identification and control of nitrate a challenge (Hatfield and Keeney, 2008).
Isotopic analysis of dissolved nitrate can provide insight to both sources and pathways
of nutrients in surface water (Kendall et al., 2007). Due to fractionation processes, nitrate
derived from atmospheric deposition, waste (human or manure), nitrate stored in soil, and
commercial fertilizer (derived from chemical fixation of N2 ) have distinct isotopic signatures
(Chang et al., 2003; Kendall et al., 2007; Mayer et al., 2002). Isotope values for oxygen
in nitrate derived from atmospheric deposition generally have values between 60 and 100

‡while oxygen isotope values for septic and manure waste, ammonium fertilizer, and soil
derived nitrate are generally between -15 and 15 ‡. Nitrate-nitrogen is greater for septic
and manure waste (between 0 and 25 ‡) than it is for ammonium fertilizer (between -10 and
78

5

‡) (Kendall et al., 2007). Advances in analysis and the use of dual isotope techniques,

specifically oxygen and nitrogen, have led to the growing popularity of the use of isotopes
to identify sources of nitrogen in rivers, streams, lakes and groundwater (Xue et al., 2009;
Kendall et al., 2007).
Interpretation of isotopic data is often complicated by microbial processes, mixing of
multiple distinct sources, and/or the existence of sources with overlapping isotopic signatures
that alter or obscure the the signature of the nitrate source. For these reasons, in the
absence of additional data, isotopic analysis is often applied to small watersheds and/or
those that are dominated by a single land use (Chang et al., 2003; Xue et al., 2009; Burns
et al., 2009). Land use can be an indication of the dominant nutrient sources, particularly
in homogeneous watersheds, but some land uses contribute nutrients from several distinct
sources. For instance, urban areas can contribute nutrients through wastewater effluent from
treatment systems, direct wastewater contributions from leaky sewer systems, fertilizer from
commercial and residential lawn application, and atmospheric deposition. The importance of
each source contribution can vary from watershed to watershed or within a single watershed.
In a state like Michigan where land use is highly variable across the state and includes both
heavily farmed areas in addition to pristine forests and sprawling urban/suburban centers,
source apportionment based on land use data alone may not adequately describe the sources
of nutrients.
Here we use nitrogen isotopic data from four major watersheds in Michigan in conjunction with spatially explicit descriptions of nutrient sources and a process inferred statistical
nutrient loading model to gain insight into, and quantify the sources contributing to, nitrate
loading across the Lower Peninsula of Michigan. The model predicts nitrate loading and
the source of nitrate. A mixing model using the results of the nitrate model is compared to
isotope sample results. The use of isotopic data with the source model validates the model
and provides clarification of isotopic data that alone would be difficult to interpret due to the
diversity of some of Michigan’s largest watersheds. The isotopic data provides additional

79

Watershed
Grand River
Saginaw Bay
Muskegon River
Boardman-Charlevoix

Percent Land Use
Ag Urban Forest
65
15
20
53
14
33
28
11
61
16
11
73

Population Density
people/km2
86
62
23
26

Sampled Area
km2
13,520
17,300
6,560
4,720

Table 3.1: Watershed characteristics
information about landscape, in-stream processes, and storage of nitrate in soils that the
nutrient model does not directly provide.

3.2
3.2.1

Methods
Sample Watersheds

Four watersheds in the Lower Peninsula of Michigan were sampled: Grand River, Saginaw
Bay, Muskegon River, and Boardman-Charlevoix (Boardman). The Grand River watershed
drains to Lake Michigan from the center of the Lower Peninsula. It has several large urban
centers including Lansing and Grand Rapids and is the most densely populated watershed
in the study. The watershed is dominated by agricultural land which is primarily row crop
agriculture. The Saginaw Bay watershed in eastern Michigan is also dominated by agricultural land; it drains to Lake Huron and contains a large amount of tile drainage due to
low slopes and poorly drained soils (CHAPTER 2). The Muskegon River and Boardman
watersheds are much smaller than the Grand River and Saginaw Bay watersheds. Located
in west-central and northwest Michigan, these watersheds drain to lake Michigan and are
dominated by forested land (Table 3.1).
Fifteen to thirty sample points were selected in each watershed to sample approximately
equal incremental watershed areas. The sample locations and incremental watershed areas
are shown in Figure 3.1.

80

3.2.2

Sample Collection

Samples were collected during spring snow melt in March 2012. A total of 88 sites were visited
over approximately 10 days. Flow was measured at each location using Acoustic Doppler
Current Profiler (ADCP) or a Marsh-McBirney Flo-Mate. Some samples were collected at
USGS gauge locations. Flows were downloaded from the USGS website to coincide with the
time of the site visit. The locations of the USGS gauges are indicated in Figure 3.1. At
all locations nitrate+nitrite samples were collected and filtered in the field with a 0.45 µm
filter. The samples were then placed on dry ice and frozen. In addition to the 88 samples,
blank and duplicate samples were collected approximately once a day or every 10-15 samples
collected. Samples were analyzed in the Michigan State University Algae Laboratory using
cadmium reduction. The nitrate+nitrite concentrations were combined with flow data to
calculate NOx fluxes.
Additional samples were collected in the spring of 2011. Samples were collected from the
most downstream accessible point of most major and several minor watersheds that drain
to the Great Lakes in the Lower Peninsula. The same sample collection methods were used.
A total of 64 sites were visited.
Isotope samples were collected during the March 2012 field campaign. Samples were
filtered in the field and frozen on dry ice. The samples were sent overnight to the Washington State University Isotope Laboratory, where they were analyzed using the “denitrifier
method” (Casciotti et al., 2002; Sigman et al., 2001). Results are reported using delta notation where the ratio of 15 N to 14 N (in parts per thousand) is reported relative to the
nitrogen isotopic ratio of N2 in air, and the oxygen 18 O isotopic ratio is reported relative
to Vienna Standard Mean Ocean Water. In addition to the internal quality assurance and
quality control measures, duplicate and standard samples were sent for analysis. Data is
included in Appendices B and C.

81

3.2.3

Statistical Nitrate Model

Results from the NOx sampling were used to calibrate a statistical model that predicts
nitrate loading and the sources of nitrate in sampled watersheds. The model is processinferred; the structure of the model equation is based on nutrient pathways and processes.
Fundamental to the basic structure of the model is a spatially detailed description of nutrient
sources. These source inputs are combined with watershed factors such as recharge and travel
distance and statistically optimized coefficients which describe attenuation of nutrients in
the basin to predict nutrient fluxes at locations within the watershed. The model and the
underlying equations are described extensively in CHAPTER 2.
The model was applied to total nitrogen loading to predict nitrate loads. The attenuation
coefficients therefore, describe both loss of nitrate due to basin retention and the proportion
of the total nitrogen pool that does not persist as nitrate. The values of the coefficients
for the nitrate model were optimized using the Matlab function fminseach. The function
uses the Nelder-Mead simplex algorithm (Nelder and Mead , 1965), a direct search algorithm
that finds the minimum value of a function. The target of the optimization was the mean
absolute difference between the natural logarithm of the observed area normalized load and
the natural logarithm of the modeled area normalized load (see CHAPTER 2).
The 2012 nitrate data was used to optimize the model and the 2011 data was used to
validate the model. A total of 6 coefficients were optimized including Bsurf (basin attenuation), Fgrd (the portion of the load entering the groundwater pathway), Bgrd (groundwater
attenuation), Btile (attenuation in tile drains), and R and BY (in-stream attenuation).

3.3
3.3.1

Results
Nitrate Sample Results

The sampled watersheds have nitrate concentrations that range from 50 to 10,123 µg/L. The
highest concentrations are in the Grand River and Saginaw Bay watersheds. The northern
82

Average

Watershed

Min Max
kg/day
4,007
88 13,151
1,043
20
7,332
6,299
120 57,031
186
33
788

Grand
Muskegon
Saginaw
Boardman-Charlevoix

Table 3.2: Total NOx loads for each sampled watershed.
R2
0.86

Adj
R2
0.87

Validation
R2
0.81

n

Error

88

0.53

Table 3.3: Summary of model fit parameters.

portions of the Muskegon watershed and the Boardman-Charlevoix watersheds had relatively
low concentrations of nitrate (Figure 3.2).
Saginaw Bay had the highest nitrate fluxes ranging from 120 to 57,031 kg/day and
averaging 6,299 kg/day. Grand River had the second highest fluxes followed by the Muskegon
River and Boardman-Charlevoix, which had the smallest fluxes. A summary of the watershed
results are presented in Table 3.2.

3.3.2

Model Results

The nitrate sample data was used to optimize the nitrate loading model. The summary
statistics and optimized parameters are shown in Tables 3.3 and 3.4. The model performed
well with an R2 of 0.86. The R2 for the validation dataset was 0.81, illustrating the consistent
performance of the model. The model was only sensitive to attenuation along the surface
pathway and basin yield.
The model predicts that chemical fertilizer and manure are the dominant sources of
nitrate in the Saginaw Bay and Grand River watersheds. Atmospheric deposition is the
dominant source predicted by the model in the Muskegon River and Boardman-Charlevoix
watersheds. Of the total atmospheric deposition, the model predicts that 5-12 % is deposited
directly to surface water and wetlands. This deposition does not experience land-based
83

attenuation. The complete source apportionment for each watershed is shown in Figure 3.4.
The coefficients from the nitrate model indicate that during melt, in areas with the
highest recharge nearly all of the nitrate infiltrates the soil (Table 3.5). The nitrate is not
very persistent in the groundwater pathway; most nitrate is removed after 1 km of travel.
Additionally, nitrate is more persistent in tile drains than in the surface runoff pathway; all of
the nitrate remains after 1 km of travel. Attenuation during in-stream travel occurs relatively
rapidly in the watershed with the lowest basin yield, but there is very little attenuation in
the watershed with the largest basin yield.
The results for the nitrate loading model are similar to those found for the total nitrogen
melt model (CHAPTER 2). Both models were sensitive to Bsurf and BY. The total nitrogen
model estimates that all nitrogen persists in tile drains after 0.1 km; the nitrate model
predicts that all nitrogen persists as nitrate in tile drains. Both models predict that in areas
with the highest recharge, no nitrogen (nitrate) travels via the groundwater pathway. The
largest difference between the models is that the total nitrogen model predicts that total
nitrogen is more persistent in the surface pathway than nitrate. This is likely due to the fact
that the nitrate model parameters describe only the portion of total nitrogen that persists
as nitrate after 1 km and the total nitrogen model includes all nitrogen species.

3.3.3

Isotope Sample Results

The isotope results are indicated by watershed in Figure 3.5. There are some general watershed patterns evident in the isotope data. The Muskegon River has values ranging from the
highest observed value of 61.3

‡to 2.61 ‡. Both the Muskegon and Boardman-Charlevoix

samples had generally higher oxygen values and lower nitrogen values than samples from
Saginaw Bay and Grand River watersheds. Saginaw Bay showed the most variable nitrogen
values with δ 15 NN O3 ranging from 1.32 to 8.4

‡.

Only one sample had a δ 18 ON O3 value greater than 30

‡. Since this sample was signifi-

cantly higher than other results, it was analyzed a second time. The result of the reanalysis

84

was within 6.5 % of the original sample. The sample was from the most upstream sample
location in the Muskegon River watershed located downstream of Houghton Lake, a 81 km2
lake in north central Michigan.
As is common with large, mixed land use watersheds, much of the isotope data plots in
the overlapping portion of the graph where source is indistinguishable without additional
information (Burns et al., 2009; Chang et al., 2003). By plotting the isotope data along with
the predominant nitrate source as predicted by the model, some trends are evident (Figure
3.6). Watersheds with manure as the predominant source, all have δ 15 NN O3 values greater
than 5. On average, watersheds dominated by atmospheric deposition have relatively heavier
oxygen and the watershed with the heaviest oxygen sampled has the greatest contribution of
atmospheric deposition at 95%. Samples dominated by fertilizer have low oxygen values and
slightly lower nitrogen isotopic values (relative to watersheds with manure as a dominant
source.
Dominant source is not perfectly predicted by isotopic ratio. This can be a result of
fractionation processes that lead to further fractionation of the original source signature,
nitrate derived from the soil which includes mineralization and nitrification of nitrate in the
soil, mixing of several sources, or model error. This is why the use of isotopic data to infer
source is difficult in large, mixed land use watersheds (Xue et al., 2009).
Using the results from the statistical model and the observed isotope values, a simple
mixing model can be constructed to predict isotopic ratios of source end members (manure/septics, ammonium fertilizer, soil, and direct atmospheric deposition). The model
assumes that the isotopic value can be predicted from the proportion of the nitrate load
each source contributes to the observed load. This is described by Equations 3.1.

δX = α1 (%DirAtm) + α2 (%Sep + %M an) + α3 (%IndAtm) + α4 (N onAg + ChemAg)
(3.1)
Each term contains the proportion of the total load contributed by the sources direct
85

atmospheric deposition (DirAtm), septic tanks (Sep), manure (Man), indirect atmospheric
deposition (IndAtm), non-agricultural fertilizer (NonAg) and chemical agricultural fertilizer
(ChemAg) and X is the isotopic ratio (delta notation) of nitrogen or oxygen. The coefficients α1 , α2 , α3 , and α4 are the end member isotope values for atmospheric deposition,
septic and manure waste, soil, and ammonium fertilizer derived nitrate respectively. The
model assumes that atmospheric deposition that falls on the landscape is transformed by
attenuation processes occurring in the soil (Burns et al., 2009).
The predicted end member values are shown in Table 3.6. The values agree well with
observed end member values compiled by Kendall et al. (2007). The R2 value is 0.32 for
the nitrogen isotope mixing model and 0.58 for the oxygen isotope mixing model. The
model does a good job of predicting the the highest observed δ 18 ON O3 value (observed in
the Muskegon River watershed). There are no major outliers in the model. Mixing model
error may be due to variability in end member values (as shown by Kendall et al. (2007)).
For instance, values of δ 15 NN O3 for the Lower Peninsula range from approximately -3.4-0.6
(Kendall et al., 2007). Other sources of mixing model error could be nitrate loading model
error or that the mixing model does not consider processes such as denitrification that may
alter isotopic signatures.
For oxygen, due to high overland flow rates during melt, some atmospheric deposition
that falls on the landscape may not experience a significant amount of processing. Higher
δ 18 ON O3 values have been correlated with higher runoff rates (Burns et al., 2009; Mayer
et al., 2002), so it has been concluded that the dominant source of atmospheric deposition
is primarily runoff. This assumption has been used to construct simple mixing models to
determine contribution of atmospheric deposition for isotope data. However, the results of
this mixing model suggest that atmospheric deposition directly to surface water predicted
by surface water and wetland area is a good predictor of the elevated oxygen isotope values
(as in the Muskegon River watershed). Ignoring this source increases uncertainty in these
predictions.

86

3.4

Discussion

3.4.1

Soil Source

Direct atmospheric deposition to surface water is indicated in samples with higher δ 18 ON O3
values. It is evident from Figure 3.6 that though some samples, on average, have higher
δ 18 ON O3 values, all but one sample plot outside the range expected for atmospherically
derived oxygen. However, many plot in the soil nitrate range. Soil nitrate cannot be directly
estimated from land use, which can make it difficult to determine dominant source from
isotope data. Other studies have found that soil derived nitrate can be a significant influence
on the isotopic ratio in watersheds (Piatek et al., 2005; Chang et al., 2003; Spoelstra et al.,
2001). Burns et al. (2009) found that isotopic data derived from forested land had an
isotopic signature consistent with soil nitrate (derived from atmospheric deposition) and is
an indication of the transformation of atmospherically derived nitrogen.
The model predicted source apportionment does not directly estimate the portion of
nitrate derived from soil nitrification. However, the model can be used to estimate the
percentage of atmospherically derived nitrogen that is deposited directly to surface water.
Assuming that all atmospheric deposition that falls on the landscape experiences a cycle
of soil nitrification before it reaches surface water, and therefore, will have the signature
of soil derived nitrate, the contribution of soil derived nitrate can be estimated for each
watershed using the model. The results of the mixing model suggest that this estimate is a
good predictor of nitrate derived from soil.

3.4.2

Denitrification

Denitrification has been observed in nitrate isotope data surface water samples in other
studies (Kaushal et al., 2011; Anisfeld et al., 2007). This process leads to relatively heavier
nitrogen and oxygen which may be apparent in the dual isotope plot (Bottcher et al., 1990).
Kaushal et al. (2011) observed a strong 2:1 relationship between δ 15 NN O3 and δ 18 ON O3 in

87

rural watersheds where the dominant source was wastewater or manure. This effect is not
apparent in the results from this study. However, in large watersheds where there is a large
contribution from soil sources, the effect may be muted due to mixing (Chang et al., 2003).
Another means of assessing the impact of from isotope data is through the relationship
between nitrate and δ 15 NN O3 . Increasing values of nitrate should result in increasing values
of δ 15 NN O3 (Chang et al., 2003). Plots of nitrate and δ 15 NN O3 for each watershed show
a weak relationship between nitrate and δ 15 NN O3 , so denitrification is not evident in this
data.
This agrees well with Kaushal et al. (2011) who concluded, based on isotopic data, that
denitrification occurred predominantly during summer baseflow. The slight negative relationship between nitrate and δ 15 NN O3 could be the result of mixing of a heavy source such
as atmospheric deposition and lighter source such as organic waste, soil derived nitrate, or
fertilizer. In addition to the isotope data, the results from the model suggest that relatively
little in-stream attenuation (denitrification and physical attenuation) occurs during spring
melt, particularly in watersheds with high basin yields where the model predicts that nearly
all nitrate persists after 1 week of in-stream travel time. During melt, discharge is high
resulting in high flows and increased stream depths. Water traveling in-stream has less contact with stream biota so denitrification is expected to occur at lower rates (Alexander et al.,
2000).

3.5

Conclusions

The isotopic data suggests that in large watersheds, a significant portion of soil derived
nitrate contributes to the nitrate load. As Burns et al. (2009) suggested, 18 ON O3 may
be a good indication of the processes occurring in the soil that transform atmospherically
derived nitrogen. Previous studies (Kaushal et al., 2011; Anisfeld et al., 2007) that have used
isotope data and land use to estimate the contribution of sources to the nitrate load may
be overestimating the contribution of organic waste and/or fertilizer if soil derived nitrate is
88

not considered, especially for samples during non-storm events. Using isotopic data in large,
mixed land-use watersheds can provide some insight to nutrient sources, but it is limited.
Even in small watersheds with a single land use, the variability in nitrogen sources (see
CHAPTER 1) may make the application of isotope data to infer source, challenging. For
instance, all sub-watersheds in the Boardman-Charlevoix are dominated by forest, but most
of the isotope sample results plot in an intermediate range that indicates mixing of several
sources.
In addition the isotope analysis helped confirm the results of the source apportionment
model. The model successfully predicted dominant source as predicted by the isotope data
(with the addition of the estimate of direct atmospheric deposition). Unlike isotope data
alone, the model can be applied in large watersheds with diverse sources to predict the
contributions of nitrate sources to the in-stream load.
Indirect atmospheric deposition is not the only source that is transformed by soil processes. A portion of the other sources may also be impacted by soil nitrification and mineralization, but these sources have overlapping isotopic signatures with the soil source, so
their transformation is not evident in this analysis. In addition, groundwater nitrate may
be impacted by these processes. Additional isotope sampling, particularly during baseflow,
should help clarify this effect.

89

!

!
! !
!
!
!

!

!
!

!
!
!!

!

!

!

Grand
Saginaw
Muskegon
BoardmanCharlevoix
Sample Site
USGS Gauge

!!
!
!!
!
!
!
!

! !!
!

!

!

! !

!
!
!
!
!
! !
!
!
!!
! !

!!
!!
!!
! ! !!
!

!
!
! !
!
!! ! !
! ! !
! !
!
!!
!
!
!!
!!
!

!

!
!
! !

!
!

Figure 3.1: Sample locations and incremental watersheds.

µg/L
1600
800
400
100

Figure 3.2: NOx concentrations for sampled sub-watersheds. Incremental watershed areas
are shown.

90

Optimization
Nitrate Load
5

-4
Log Model kg/day

Log Model kg/m2 /day

Nitrate Area Normalized Load

-5
-6
-7

y = 0.6779x - 2.0935
R2 = 0.7177
-7

3
2
y = 0.8037x + 0.3867
R2 = 0.856

1
0

-8
-8

4

0

-6

-5
-4
Log Observed kg/m2 /day

1
2
3
4
Log Observed kg/day

5

Validation
Nitrate Load

-4

5
Log Model kg/day

Log Model kg/m2 /day

Nitrate Area Normalized Load

-5
-6
-7
y = 0.5375x - 2.683
R2= 0.6905

-8
-9

4
3
2
1
y = 0.7367x + 0.761
R 2 = 0.8101

0
-1

-9

-8

-7

-6

-5
-4
Log Observed kg/m2 /day

-1

0
1
2
3
4
Log Observed kg/day

5

Figure 3.3: Graphs of observed and modeled nitrate loads and area normalized loads for
the calibration and validation data sets.

Bsurf
-8.71E-04

Btile
-3.78E-12

37.38

0.00

R
BY
-6.18E-14 1.00E+00
Sensitivity
0.00
0.00

Fgrd
4.26E-08

Bgrd
-3.96E-03

0.00

-82.95

Table 3.4: Summary of coefficients associated with each basin parameter.

91

100%
80%

Atm Direct
Septic

60%

Point
NonAg

40%

Manure
ChemAg

20%

Atm Indirect

0%
BC

MR

GR

SB

Figure 3.4: Model predicted source of exported nitrate in watersheds.
Bsurf
(%/km)
42

Btile
(%/km)
100

Rmin
(%/day)
4

Rmax
(%/day)
100

Fgrdmax
(%)
0

Bgrd
(%/km)
2

Table 3.5: Summary of calculated basin parameters. Reduction factor units (Bsurf, Btile,
and Bgrd ) are the percent of the nutrient load that remains after the first kilometer of
travel. Units of Rmin and Rmax are the percent of the load that remains after the first day
of in-stream travel for the sub-watershed with the highest relative basin yield (Rmax ) and
the lowest relative basin yield (Rmin). The bold values indicate the coefficients with
relatively high model sensitivity.
Model Predicted

Reported by Kendall et al. (2007)

Source

δ 15 NN O3

δ 18 ON O3

δ 15 NN O3

δ 18 ON O3

Septic and Manure Waste
Atmospheric Deposition
Fertilizer
Soil

8.27
0.00
5.74
2.68

2.57
100.39
0.00
3.75

0-25
-15-15
-10-5
2-8

-15-15
60-98
-15-15
-15-15

Table 3.6: End member isotope values predicted by mixing model and compared to
Kendall et al. (2007).

92

70

GR

50

MR

3

δ18O NO (‰)

BC

SB

30

10
-2
-10

3
15N
δ
NO3(‰)

8

Figure 3.5: Measured isotopes by watershed

70

105

Soil
Atmospheric
Deposition

Chem. Fertilizer

50

Manure/Septic

55

3

δ18ONO (‰)

80

30
30
Soil
5
Fert.
-20

-20

Manure/
Septic

10
-10

0
20
15NNO
δ
3 (‰)

-5

5
δ15NNO

15

3 (‰)

Figure 3.6: Isotope data with dominant source indicated (isotopic classification from
Kendall et al. (2007)).

93

105
Atmospheric
Deposition

80
55

3

δ18ONO (‰)

Manure and
SepƟc Waste
Atmospheric
DeposiƟon
Ammonium
FerƟlizer
Soil

30
Soil
5

Manure

Fert.
-20

-20

0
20
δ15NNO (‰)
3

Figure 3.7: Predicted isotope values plotted with values reported by Kendall et al. (2007)

94

Muskegon

Grand

8

3

δ15N NO (‰)

3

δ15N NO (‰)

10

6
4
2
0
0.0

5.0
10.0
NOx (mg/L)

15.0

7
6
5
4
3
2
1
0
-1 0.0
-2

1.0

1.5

NOx (mg/L)

Boardman-Charlevoix
7

9
8
7
6
5
4
3
2
1
0

6
3

δ15N NO (‰)

3

δ15N NO (‰)

Saginaw

0.5

5
4
3
2
1
0

0.0

5.0
NOx (mg/L)

0.0

10.0

Figure 3.8: Nitrate vs. δ 15 NN O3 .

95

0.2
0.4
NOx (mg/L)

0.6

APPENDICES

96

Appendix A
Total Nitrogen, Total Phosphorus, and Discharge Results

97

ID

Lat

Long

Date

SW-1
SW-10
SW-101
SW-102
SW-103
SW-104
SW-11
SW-12
SW-13
SW-14
SW-15
SW-16
SW-17
SW-18
SW-19
SW-2
SW-20
SW-21
SW-23
SW-24
SW-25
SW-26
SW-27
SW-28
SW-29
SW-3

42.087
44.259
43.327
43.310
43.245
43.589
44.356
44.623
44.715
44.765
44.900
45.595
45.559
45.068
44.415
42.184
44.050
44.072
43.064
42.858
42.749
42.596
42.297
42.064
41.964
42.646

-86.476
-85.939
-83.746
-83.956
-84.100
-84.221
-86.049
-86.223
-86.125
-85.621
-85.411
-84.474
-84.403
-83.436
-83.331
-86.371
-83.688
-84.019
-82.597
-82.538
-82.490
-82.909
-83.191
-83.253
-83.546
-86.118

Sample
Type

10/1/10
10/2/10
10/15/10
10/13/10
10/13/10
10/12/10
10/2/10
10/3/10
10/3/10
10/3/10
10/3/10
10/3/10
10/3/10
10/4/10
10/4/10
10/1/10
10/4/10
10/4/10
10/12/10
10/12/10
10/12/10
10/12/10
10/12/10
10/4/10
10/4/10
10/1/10

Flow
cfs
2022.82
1510
69
223
154
617
106.3
806.23
170.22
235.9
452.73
919.95
408.94
444.61
1146.31
283
37.79
178
23.41
1.659860294
48.03
0
11.3
271.92
167
1540

TN
µg/L
1828.07832
190.00984
940.17968
3376.59024
827.32592
1013.70352
525.48952
492.87344
329.79304
483.55456
176.03152
250.58256
395.0252
1097.57344
180.69096
818.00704
352.06304
553.44616
519.80288
283.19864
362.40912
3544.33008
5583.1104
809.71536
1293.26992
1688.29512

TP
µg/L
49.71576
15.37661
25.47636
104.25441
27.09232
35.17212
37.59606
20.22449
9.72075
15.37661
7.7008
8.10479
10.52873
18.20454
14.16464
34.76813
18.20454
14.16464
40.42399
27.09232
45.27187
160.00503
395.5312
56.1796
97.38658
51.73571

Table A.1: 2010 Baseflow Results. ND = Non Detect. DUP = Duplicate Sample.

98

ID

Lat

Long

Date

SW-31
SW-32
SW-33
SW-34
SW-35
SW-36
SW-37
SW-39
SW-4
SW-33
SW-34
SW-35
SW-36
SW-37
SW-39
SW-4
SW-40
SW-41
SW-43
SW-51
SW-52
SW-55
SW-56
SW-57

42.417
42.421
42.799
42.797
42.916
42.907
42.904
42.949
42.965
42.799
42.797
42.916
42.907
42.904
42.949
42.965
42.994
42.962
43.197
43.490
43.559
43.905
43.999
44.904

-86.241
-86.245
-86.059
-86.144
-86.147
-85.781
-85.772
-85.849
-85.675
-86.059
-86.144
-86.147
-85.781
-85.772
-85.849
-85.675
-86.018
-86.203
-86.186
-86.438
-86.503
-86.308
-86.372
-85.975

Sample
Type

10/1/10
10/1/10
10/1/10
10/1/10
10/1/10
10/1/10
10/1/10
10/1/10
10/2/10
10/1/10
10/1/10
10/1/10
10/1/10
10/1/10
10/1/10
10/2/10
10/1/10
10/1/10
10/2/10
10/2/10
10/2/10
10/2/10
10/2/10
10/2/10

Flow
cfs
38.14344
68.940736
5.6
9.65
2.79
4.24
22.62
11.08
1790
5.6
9.65
2.79
4.24
22.62
11.08
1790
7.08
1.14
23.24
11.97
44.44
12.15
23.36
55.47

TN
µg/L
421.6116377
1656.70624
4131.41952
1218.71888
3607.50784
1265.31328
939.15248
1330.54544
1237.35664
4131.41952
1218.71888
3607.50784
1265.31328
939.15248
1330.54544
1237.35664
2429.14608
81.81552
687.54272
3611.61664
1056.66568
3460.46016
855.28256
255.242

TP
µg/L
36.38409
27.49631
57.39157
33.96015
53.75566
40.02
20.62848
62.64344
54.96763
57.39157
33.96015
53.75566
40.02
20.62848
62.64344
54.96763
30.32424
50.52374
25.47636
24.66838
18.20454
27.09232
21.43646
9.31676

Table A.2: 2010 Baseflow Results (continued). ND = Non Detect. DUP = Duplicate
Sample

99

ID

Lat

Long

Date

SW-59
SW-6
SW-60
SW-62
SW-63
SW-64
SW-65
SW-66
SW-67
SW-69
SW-7
SW-70
SW-71
SW-73
SW-8
SW-9
SW-91
SW-94
SW-97

45.017
43.463
44.745
45.100
45.180
45.214
45.365
45.746
45.488
45.431
43.944
45.337
45.162
44.260
44.030
44.208
42.588
42.095
41.931

-85.618
-86.232
-85.548
-85.097
-85.164
-85.013
-84.962
-84.828
-84.076
-83.839
-86.278
-83.685
-83.410
-83.527
-86.506
-86.246
-82.903
-83.211
-83.340

Sample
Type

10/2/10
10/2/10
10/3/10
10/3/10
10/3/10
10/3/10
10/3/10
10/3/10
10/3/10
10/2/10
10/2/10
10/3/10
10/4/10
10/4/10
10/2/10
10/2/10
10/12/10
10/4/10
10/3/10

Flow
cfs
2.95
272
10.29
145
2.94
61.66
55.79
62.31
84.39
14.32
431
2.2
0.92
56.12
101.35
158.56
0.24
4.72
0.54

TN
µg/L
1712.61952
413.66296
991.43352
1488.9664
562.76504
548.78672
409.00352
590.72168
655.95384
780.73152
287.85808
683.91048
553.44616
325.1336
576.74336
92.1616
716.52656
948.47136
87.50216

TP
µg/L
10.52873
26.28434
21.43646
8.91277
12.54868
10.93272
8.91277
13.35666
13.76065
27.49631
27.9003
14.97262
19.41651
18.60853
14.97262
12.54868
56.1796
73.14718
54.15965

Table A.3: 2010 Baseflow Results (continued). ND = Non Detect. DUP = Duplicate
Sample

100

ID

Lat

Long

Date

J-201
J-202
J-203
SW-1
SW-10
SW-101
SW-102
SW-103
SW-104
SW-11
SW-11
SW-11
SW-12
SW-13
SW-14
SW-15
SW-15
SW-15
SW-16
SW-17
SW-18
SW-19
SW-2
SW-20
SW-21
SW-21
SW-21
SW-23
SW-24
SW-25

45.136
45.033
45.036
42.087
44.259
43.327
43.310
43.245
43.589
44.356
44.356
44.356
44.623
44.715
44.765
44.900
44.900
44.900
45.595
45.559
45.068
44.415
42.184
44.050
44.072
44.072
44.072
43.064
42.858
42.749

-85.081
-85.064
-84.969
-86.476
-85.939
-83.746
-83.956
-84.100
-84.221
-86.049
-86.049
-86.049
-86.223
-86.125
-85.621
-85.411
-85.411
-85.411
-84.474
-84.403
-83.436
-83.331
-86.371
-83.688
-84.019
-84.019
-84.019
-82.597
-82.538
-82.490

3/19/11
3/19/11
3/19/11
3/5/11
3/20/11
3/7/11
3/6/11
3/22/11
3/6/11
3/20/11
3/20/11
3/20/11
3/20/11
3/20/11
3/19/11
3/19/11
3/19/11
3/19/11
3/23/11
3/23/11
3/22/11
3/22/11
3/5/11
3/22/11
3/22/11
3/22/11
3/22/11
3/7/11
3/7/11
3/22/11

Sample
Type

Flow
cfs
69.49926407
40.15277604
148.1097121
6750.045393
1360
2300
3170
1479.296073
1160
34.82026137

Blank
DUP
288.4855123
157.8212455
381.0452537
481.8333125
Blank
DUP
835.6862728
364.0942137
1555.893585
1722.543498
823
646.8940646
777
Blank
DUP
2216.348482
1178.76826

TN
µg/L
717.883464
1637.296572
1755.6817
3400.097016
532.877448
5195.23044
2683.292124
3158.583016
2357.055488
570.8274
575.571144
717.883464
665.70228
290.946504
999.054532
252.996552
ND
300.433992
271.971528
267.227784
940.839432
428.51508
1238.218292
2758.527272
295.690248
72.73428
1373.019684
5449.811368
1413.966924
2102.183484

TP
µg/L
53.66032
30.032944
48.854752
118.535488
42.046864
161.385136
143.364256
232.667728
75.68584
63.67192
43.64872
61.6696
45.250576
29.63248
52.859392
38.042224
8.808352
29.63248
38.442688
33.236656
82.493728
65.67424
45.117088
148.570288
38.843152
8.808352
112.928992
191.8204
218.251024
166.190704

Table A.4: 2011 Melt Results (continued). ND = Non Detect. DUP = Duplicate Sample

101

ID

Lat

Long

Date

SW-26
SW-27
SW-28
SW-29
SW-3
SW-3
SW-3
SW-31
SW-32
SW-33
SW-34
SW-34
SW-35
SW-36
SW-37
SW-39
SW-4
SW-40
SW-40
SW-40
SW-41
SW-43
SW-48
SW-5
SW-51
SW-52
SW-55
SW-56
SW-57
SW-59
SW-6
SW-6
SW-60
SW-62

42.596
42.297
42.064
41.964
42.646
42.646
42.646
42.417
42.421
42.799
42.797
42.797
42.916
42.907
42.904
42.949
42.965
42.994
42.994
42.994
42.962
43.197
43.416
43.318
43.490
43.559
43.905
43.999
44.904
45.017
43.463
43.463
44.745
45.100

-82.909
-83.191
-83.253
-83.546
-86.118
-86.118
-86.118
-86.241
-86.245
-86.059
-86.144
-86.144
-86.147
-85.781
-85.772
-85.849
-85.675
-86.018
-86.018
-86.018
-86.203
-86.186
-86.277
-86.038
-86.438
-86.503
-86.308
-86.372
-85.975
-85.618
-86.232
-86.232
-85.548
-85.097

3/4/11
3/4/11
3/4/11
3/4/11
3/6/11
3/6/11
3/6/11
3/5/11
3/5/11
3/5/11

Sample
Type

Flow
cfs
739.135974
438.2550137
969.7407476
3630
3060

Blank
DUP
686.7643233
573.5101872
140
DUP

3/5/11
3/5/11
3/6/11
3/6/11
3/6/11
3/6/11
3/6/11
3/6/11
3/6/11
3/6/11
3/21/11
3/21/11
3/21/11
3/21/11
3/21/11
3/21/11
3/20/11
3/20/11
3/19/11
3/21/11
3/21/11
3/19/11
3/19/11

47.32165338
118.3041334
187.5561948
199.1747202
132.7831468
6970
165.2726402
Blank
DUP
7.416080007
79.74051741
12.67796535
3386.676537
71.75940273
113.9604294
55.16150939
169.6869735
66.00311206
6.356640006
815
DUP
19.35243735
180

TN
µg/L
2249.365416
1819.930716
1919.54934
7045.416468
2243.21
1.57812
2129.355776
1181.688676
1969.358652
7013.66564
3158.583016

TP
µg/L
113.329456
116.533168
77.287696
213.445456
85.3
8.808352
87.299296
129.74848
92.505328
512.992528
149.371216

3485.901352
2284.817628
3652.514544

166.591168
111.327136
221.8552

3140.190192
82.221768
2649.42116
1181.688676
1486.078916
765.320904

28.831552
8.407888
114.931312
58.065424
62.070064
35.63944

2449.475892
1863.706112
2976.739496
1425.200868
252.996552
1660.016196
988.276872
969.301896
974.04564
1172.991812

165.389776
64.472848
191.8204
105.320176
46.451968
47.252896
82.093264
73.68352
49.65568
51.257536

Table A.5: 2011 Melt Results (continued). ND = Non Detect. DUP = Duplicate Sample

102

ID

Lat

Long

Date

SW-63
SW-64
SW-65
SW-66
SW-67
SW-69
SW-7
SW-70
SW-71
SW-73
SW-8
SW-9
SW-91
SW-94
SW-97

45.180
45.214
45.365
45.746
45.488
45.431
43.944
45.337
45.162
44.260
44.030
44.208
42.588
42.095
41.931

-85.164
-85.013
-84.962
-84.828
-84.076
-83.839
-86.278
-83.685
-83.410
-83.527
-86.506
-86.246
-82.903
-83.211
-83.340

Sample
Type

3/19/11
3/23/11
3/23/11
3/23/11
3/23/11
3/23/11
3/20/11
3/23/11
3/23/11
3/22/11
3/20/11
3/20/11
3/4/11
3/4/11
3/4/11

Flow
cfs
11.47726668
94.78456542
109.8639281
36.72725337
115.6555334
33.5136187
1380
21.01222669
38.4929867
173.5715868
343.5763923
256.5610536
9.888106676
61.44752006
63.91954673

TN
µg/L
1002.508104
575.571144
637.239816
698.908488
528.133704
660.958536
3430.431804
943.71
618.26484
920.782756
1103.4169
798.527112
2085.58038
1025.145124
7886.514536

TP
µg/L
71.280736
45.250576
51.257536
50.857072
45.250576
53.66032
59.66728
8.808352
42.847792
55.66264
98.512288
78.088624
94.507648
80.491408
111.327136

Table A.6: 2011 Melt Results (continued). ND = Non Detect. DUP = Duplicate Sample

103

ID

Lat

Long

Date

J-201
J-202
J-203
SW-1
SW-10
SW-10
SW-10
SW-101
SW-102
SW-102
SW-102
SW-103
SW-104
SW-11
SW-12
SW-13
SW-14
SW-15
SW-16
SW-17
SW-18
SW-18
SW-19
SW-2
SW-2
SW-2
SW-20
SW-21
SW-24
SW-25
SW-26
SW-27
SW-28
SW-28
SW-28

45.136
45.033
45.036
42.087
44.259
44.259
44.259
43.327
43.310
43.310
43.310
43.245
43.589
44.356
44.623
44.715
44.765
44.900
45.595
45.559
45.068
45.068
44.415
42.184
42.184
42.184
44.050
44.072
42.858
42.749
42.596
42.297
42.064
42.064
42.064

-85.081
-85.064
-84.969
-86.476
-85.939
-85.939
-85.939
-83.746
-83.956
-83.956
-83.956
-84.100
-84.221
-86.049
-86.223
-86.125
-85.621
-85.411
-84.474
-84.403
-83.436
-83.436
-83.331
-86.371
-86.371
-86.371
-83.688
-84.019
-82.538
-82.490
-82.909
-83.191
-83.253
-83.253
-83.253

6/24/11
6/24/11
6/24/11
6/3/11
6/25/11
6/25/11
6/25/11
6/1/11
6/1/11
6/1/11
6/1/11
6/1/11
6/1/11
6/25/11
6/25/11
6/25/11
6/25/11
6/24/11
6/27/11
6/27/11
6/28/11
6/28/11
6/28/11
6/3/11
6/3/11
6/3/11
6/28/11
6/28/11
6/1/11
6/1/11
6/1/11
6/2/11
6/2/11
6/2/11
6/2/11

Sample
Type

Flow
cfs
73.45450674
38.9167627
164.9901228
9592.275713
2220

Blank
DUP
2370
2990
Blank
DUP
2010
2880
28.95802669
225.9079229
166.5086535
522.3745498
1216.37838
1151.540652
1467.889436
1963.530783
DUP
1624.933759
991
Blank
DUP
162.0590055
340
112.0887521
247.9795896
0
479.0787685
3446.570211
Blank
DUP

TN
µg/L
636.33824
1090.83456
489.35744
1127.604916
594.37472
69.83072
678.30176
3496.52608
1369.513367
118.78816
1015.962826
895.0048
1300.65216
888.11936
685.29568
608.36256
874.13152
559.40512
517.4416
657.32
1132.90656
692.2896

TP
µg/L
13.5684
25.8213
29.5343
91.1701
21.37
0.2016
25.8213
89.6849
65.5504
0.5729
87.65
88.571
55.5253
29.5343
36.589
26.9352
21.737
13.39
15.43
23.5935
33.2473
24.3361

1167.87616
125.78208
83.71008
2231.886348
930.08288
433.40608
1329.119968
1112.616133
1398.311145
1404.59
160.75168
1243.89

64.4365
0.2016
53.2975
35.1038
36.2177
61.8374
94.8831
98.2248
98.9674
66.293
0.9442
46.2428

Table A.7: 2011 Rain Results. ND = Non Detect. DUP = Duplicate Sample

104

ID

Lat

Long

Date

SW-29
SW-3
SW-31
SW-32
SW-33
SW-34
SW-35
SW-36
SW-37
SW-39
SW-4
SW-40
SW-41
SW-41
SW-41
SW-43
SW-48
SW-5
SW-51
SW-52
SW-55
SW-56
SW-57
SW-59
SW-56
SW-57
SW-59
SW-6
SW-60
SW-62
SW-62
SW-62
SW-63
SW-64

41.964
42.646
42.417
42.421
42.799
42.797
42.916
42.907
42.904
42.949
42.965
42.994
42.962
42.962
42.962
43.197
43.416
43.318
43.490
43.559
43.905
43.999
44.904
45.017
43.999
44.904
45.017
43.463
44.745
45.100
45.100
45.100
45.180
45.214

-83.546
-86.118
-86.241
-86.245
-86.059
-86.144
-86.147
-85.781
-85.772
-85.849
-85.675
-86.018
-86.203
-86.203
-86.203
-86.186
-86.277
-86.038
-86.438
-86.503
-86.308
-86.372
-85.975
-85.618
-86.372
-85.975
-85.618
-86.232
-85.548
-85.097
-85.097
-85.097
-85.164
-85.013

6/2/11
6/3/11
6/3/11
6/3/11
6/3/11
6/3/11
6/3/11
6/4/11
6/4/11
6/4/11
6/4/11
6/4/11
6/4/11
6/4/11
6/4/11
6/26/11
6/26/11
6/26/11
6/26/11
6/26/11
6/26/11
6/25/11
6/25/11
6/25/11
6/25/11
6/25/11
6/25/11
6/26/11
6/24/11
6/24/11
6/24/11
6/24/11
6/24/11
6/27/11

Sample
Type

Flow
cfs
3800
5610
277.9264269
312.5348003
74.8670934
36.30347737
49.51116271
60.67059739
67.73353073
38.10452537
10817.55339
44.35522138
4.237760004

Blank
DUP
58.51640272
11.65384001
3368.242281
21.64789069
64.97898673
14.83216001
89.69925342
88.56918408
6.709786673
89.69925342
88.56918408
6.709786673
DUP
13.06642668
192
Blank
DUP
28.60488003
88.85170142

TN
µg/L
1845.96096
839.05344
125.6736
643.22368
3524.50176
209.60064
1562.607198
1286.66432
559.29664
1258.68864
1539.594191
2257.13885
790.20448
83.81856
888.11936
5146.65728
867.02912
993.02816
2307.275172
447.39392
7776.3712
811.07776
419.52672
517.33312
811.07776
419.52672
517.33312
1097.93696
1125.91264
335.4912
132.776
139.66144
923.08896
734.25312

TP
µg/L
111.9629
78.9172
52.1836
34.18
20.6231
44.3863
72.2338
61.8374
57.3818
95.6257
94.5118
59.6096
22.8509
ND
15.0536
26.9352
17.2814
23.5935
46.2428
15.0536
65.92
71.4912
23.5935
25.8213
71.4912
23.5935
25.8213
40.302
34.3612
27.68
1.6868
28.7917
26.9352
25.8213

Table A.8: 2011 Rain Results (continued). ND = Non Detect. DUP = Duplicate Sample

105

ID

Lat

Long

Date

SW-65
SW-65
SW-65
SW-66
SW-67
SW-69
SW-7
SW-70
SW-71
SW-73
SW-8
SW-9
SW-91
SW-94
SW-97

45.365
45.365
45.365
45.746
45.488
45.431
43.944
45.337
45.162
44.260
44.030
44.208
42.588
42.095
41.931

-84.962
-84.962
-84.962
-84.828
-84.076
-83.839
-86.278
-83.685
-83.410
-83.527
-86.506
-86.246
-82.903
-83.211
-83.340

6/27/11
6/27/11
6/27/11
6/27/11
6/27/11
6/27/11
6/25/11
6/27/11
6/27/11
6/28/11
6/25/11
6/25/11
6/1/11
6/2/11
6/2/11

Sample
Type

Flow
cfs
148.0743975

Blank
DUP
96.93876009
255.2190962
77.40974941
1470
38.03389604
167.0383735
118.2335041
345.730587
267.5792296
3.53146667
25.14404269
34.2552267

TN
µg/L
825.17408
ND
797.1984
846.15584
948.966464
867.1376
1118.91872
923.08896
559.40512
979.04032
566.39904
636.33824
1374.921506
419.41824
4503.65056

TP
µg/L
27.6778
2.4294
18.7666
11.7119
16.9101
23.9648
46.6141
29.5343
11.7119
18.7666
10.9693
29.163
91.5414
111.5916
84.4867

Table A.9: 2011 Rain Results (continued). ND = Non Detect. DUP = Duplicate Sample

106

ID

Lat

Long

Date

BC10
BC11
BC12
BC13
BC14
BC15
BC17
BC9
GR1
GR11
GR12
GR13
GR14
GR15
GR17
GR20
GR21
GR21
GR21
GR23
GR24
GR25
GR26
GR27
GR28
GR28
GR28
GR3
GR30
GR31
GR5
GR6
GR7
GR8
GR9
MR10
MR10
MR10
MR11

45.508
45.440
45.166
44.846
44.699
44.528
44.657
45.659
43.130
43.015
42.966
43.067
42.775
43.063
42.312
42.283
42.536
42.536
42.536
42.728
42.750
42.857
42.828
43.109
42.971
42.971
42.971
42.609
42.616
43.082
43.097
43.226
42.877
43.000
43.145
43.402
43.402
43.402
43.306

-84.763
-84.785
-85.240
-85.307
-85.306
-85.949
-85.439
-84.493
-85.174
-84.426
-85.352
-84.842
-85.359
-86.067
-84.387
-84.410
-84.624
-84.624
-84.624
-84.478
-84.555
-84.913
-84.759
-84.693
-85.069
-85.069
-85.069
-85.094
-85.236
-85.590
-84.710
-85.302
-84.384
-84.908
-84.503
-85.563
-85.563
-85.563
-86.114

8/19/11
8/19/11
8/19/11
8/20/11
8/20/11
8/20/11
8/20/11
8/19/11
10/9/11
10/9/11
10/9/11
10/9/11
10/8/11
10/8/11
10/8/11
10/8/11
10/8/11
10/8/11
10/8/11
10/10/11
10/10/11
10/9/11
10/9/11
10/9/11
10/9/11
10/9/11
10/9/11
10/8/11
10/8/11
10/9/11
10/9/11
10/9/11
10/8/11
10/9/11
10/9/11
9/11/11
9/11/11
9/11/11
9/11/11

Sample
Type

Flow
cfs
75.92653236
81.18841762
27.89858631
81.93002561
25.07341301
65.33213249
93
-0.03
24.40243435
15.53845313
105.1317613
66.74471914
56.39752194
28.81676763
108.4866546
109
433

Blank
DUP
70
546
628
49
44
987
Blank
DUP
49.5464767
158
121
9.534959877
128.545385
15.5737678
19.49369575
33.9727089
58.69297525
Blank
DUP
30.51187161

TN
µg/L
498.89816
664.04
368.78792
438.84728
386.3
263.69888
303.7328
608.99144
664.03808
503.9024
153.6056
258.69464
1363.73008
543.03472
322.84816
2493.256
362.88208
8.48264
342.86512
442.94992
412.92448
763.22128
883.32304
453.86
703.1704
ND
513.00928
503.9024
282.81424
313.74128
2710.77224
513.00928
659.03384
513.91088
343.76672
653.128
ND
593.08
553.9448

TP
µg/L
35.29952
37.95798
41.376
32.26128
30.74216
50.49072
33.7804
78.21466
82.01246
475.84432
56.5672
55.04808
74.03708
136.321
126.44672
253.82
96.82388
1.4991
82.01246
85.43048
229.45
74.79664
70.23928
142.39748
64.1628
ND
72.13818
74.03708
63.02346
56.94698
41.376
62.2639
219.11304
41.75578
124.54782
40.23666
1.11932
33.02084
75.17642

Table A.10: 2011 Baseflow Results. ND = Non Detect. DUP = Duplicate Sample

107

ID

Lat

Long

Date

MR13
MR14
MR15
MR16
MR17
MR2
MR20
MR21
MR21
MR21
MR23
MR24
MR3
MR4
MR5
MR6
MR7
MR8
MR9
SB12
SB13
SB15
SB16
SB17
SB19
SB2
SB20
SB21
SB22
SB24
SB27
SB28
SB29
SB3
SB30

43.699
43.540
44.336
43.490
43.608
44.082
44.201
43.901
43.901
43.901
43.435
43.289
44.139
44.006
43.853
44.111
43.934
43.779
44.409
43.657
43.211
43.250
43.450
43.872
43.015
43.302
43.218
43.160
43.111
43.038
43.328
43.867
43.627
43.675
43.379

-85.460
-85.305
-84.922
-85.449
-85.481
-85.014
-85.053
-85.255
-85.255
-85.255
-85.666
-86.223
-84.899
-85.111
-85.444
-85.028
-85.129
-85.500
-84.805
-84.267
-83.315
-83.873
-83.441
-84.321
-84.180
-84.144
-84.106
-83.351
-83.519
-83.774
-83.758
-84.545
-84.708
-84.388
-84.656

9/11/11
9/11/11
8/21/11
9/11/11
8/22/11
8/21/11
8/21/11
8/21/11
8/21/11
8/21/11
8/21/11
9/11/11
8/21/11
8/21/11
8/21/11
8/21/11
8/21/11
9/11/11
8/21/11
9/17/11
10/7/11
10/7/11
10/7/11
9/17/11
10/7/11
10/7/11
10/7/11
10/8/11
10/8/11
10/8/11
10/7/11
9/17/11
9/17/11
9/17/11
9/17/11

Sample
Type

Flow
cfs
14.09055182
62.71884719
24.54369302
2.419054636
12.88985317
331.7459744
102
480

TN
µg/L
103.5632
518.91512
123.58016
463.86848
473.87696
443.85152
188.19
509.96772

1220
3.4
144.1897821
28.07515964
58.23388458
16.09
5.61
556.9829155
49.08738604
0.0918
32.34823425
303.3176681
47.56885539
73.06604439
169
0.3186
214.9603732
57
121
220
89
52
132
136.5618142
119

488.54992
669.04232
273.70736
213.65648
629.0084
288.72008
533.92784
308.73704
498.89816
538.93208
674.04656
611.5
1113.51808
438.84728
418.83032
651.53
709.07624
824.17376
553.9448
543.93632
603.0856
618.99992
304.35684
383.80064
508.90664
634.02

Blank
DUP

10/3/10

TP
µg/L
67.20104
45.55358
47.0727
48.9716
29.60282
42.51534
39.85688
37.19842
3.398
37.5782
48.21204
41.75578
33.40062
45.1738
37.19842
99.48234
31.50172
31.50172
39.09732
57.32676
89.98784
96.06432
49.73116
53.14918
58.08632
42.13556
50.49072
64.54258
48.9716
75.17642
47.83226
34.16018
50.49072
56.94698
64.54258

Table A.11: 2011 Baseflow Results (continued). ND = Non Detect. DUP = Duplicate
Sample

108

ID

Lat

Long

Date

SB31
SB31
SB31
SB32
SB33
SB33
SB33
SB4
SB5
SB6
SB7
SB8
SB9
SW-11
SW-12
SW-13
SW-13
SW-13
SW-14
SW-15
SW-20
SW-21
SW-36
SW-36
SW-39
SW-4
SW-5
SW-62
SW-62
SW-62
SW-64
SW-65
SW-66

43.563
43.563
43.563
43.610
43.409
43.409
43.409
43.597
43.397
43.292
43.741
43.298
43.380
44.356
44.623
44.715
44.715
44.715
44.765
44.900
44.050
44.072
42.907
42.907
42.949
42.965
43.318
45.100
45.100
45.100
45.214
45.365
45.746

-84.370
-84.370
-84.370
-84.242
-83.964
-83.964
-83.964
-84.944
-84.836
-83.992
-85.127
-84.143
-84.094
-86.049
-86.223
-86.125
-86.125
-86.125
-85.621
-85.411
-83.688
-84.019
-85.781
-85.781
-85.849
-85.675
-86.038
-85.097
-85.097
-85.097
-85.013
-84.962
-84.828

9/17/11
9/17/11
9/17/11
9/17/11
10/7/11
10/7/11
10/7/11
9/17/11
9/17/11
10/7/11
9/17/11
10/7/11
10/7/11
8/20/11
8/20/11
8/20/11
8/20/11
8/20/11
8/20/11
8/19/11
9/17/11
9/17/11
10/9/11
10/9/11
10/8/11
10/9/11
9/11/11
8/19/11
8/19/11
8/19/11
8/19/11
8/19/11
8/19/11

Sample
Type

Flow
cfs
69

Blank
DUP
358
5240
Blank
DUP
96.55029743
62.4010152
19.45838108
32.24229025
6.568527916
9.358386546
6.82
153.2656514
99.58735872
Blank
DUP

DUP

235.195677
293.8180229
11.93635718
140
34.213
22.036
17.44544511
2100
1086.879482
156

Blank
DUP
72.14786307
35.63249821
6.886359911

TN
µg/L
493.89392
498.89816
458.86424
1388.09752
43.05968
1507.61786
458.86424
398.81336
649.03
523.91936
418.83032
684.05504
483.88544
383.80064
306.24

TP
µg/L
26.1848
1.87888
42.89512
51.63006
73.6573
6.05646
78.97422
34.16018
35.29952
46.69292
42.13556
36.43886
113.91398
52.00984
61.12456
23.14656

214.40208
433.84304
128.5844
390.02804
63.52928
1124.42816

34.16018
33.02084
34.91974
25.42524
34.16018
78.97422

596.48
1183.57744
313.74128
228.6692
3.4784
202.7464
388.80488
278.7116
433.84304

99.10256
68.72016
55.80764
31.12194
1.11932
45.1738
49.35138
31.50172
22.76678

Table A.12: 2011 Baseflow Results (continued). ND = Non Detect. DUP = Duplicate
Sample

109

ID

Lat

Long

Date

BC10
BC11
BC12
BC13
BC14
BC15
BC17
BC9
GR1
GR11
GR11
GR11
GR12
GR13
GR14
GR15
GR17
GR20
GR21
GR23
GR24
GR25
GR26
GR27
GR28
GR28
GR3
GR30
GR30
GR30
GR31
GR31
GR5
GR6
GR7
GR8
GR9
MR10

45.508
45.440
45.166
44.846
44.699
44.528
44.657
45.659
43.130
43.015
43.015
43.015
42.966
43.067
42.775
43.063
42.312
42.283
42.536
42.728
42.750
42.857
42.828
43.109
42.971
42.971
42.609
42.616
42.616
42.616
43.082
43.082
43.097
43.226
42.877
43.000
43.145
43.402

-84.763
-84.785
-85.240
-85.307
-85.306
-85.949
-85.439
-84.493
-85.174
-84.426
-84.426
-84.426
-85.352
-84.842
-85.359
-86.067
-84.387
-84.410
-84.624
-84.478
-84.555
-84.913
-84.759
-84.693
-85.069
-85.069
-85.094
-85.236
-85.236
-85.236
-85.590
-85.590
-84.710
-85.302
-84.384
-84.908
-84.503
-85.563

3/9/12
3/9/12
3/10/12
3/10/12
3/10/12
3/11/12
3/10/12
3/9/12
3/2/12
3/4/12
3/4/12
3/4/12
3/2/12
3/4/12
3/3/12
3/2/12
3/3/12
3/3/12
3/3/12
3/7/12
3/7/12
3/2/12
3/3/12
3/4/12
3/2/12
3/3/12
3/3/12
3/3/12
3/3/12
3/3/12
3/2/12
3/2/12
3/4/12
3/2/12
3/3/12
3/4/12
3/4/12
3/8/12

Sample
Type

Flow
cfs
166.1556635
209.9812062
116.2206777
104.531512
105.4850089
300.0690059
172
110.4996963
105.8734706
248.0857675

Blank
DUP

DUP

406.0484206
698.3481925
811.6024354
688.0363001
226.01408
347
1710
500
2200
3350
660
1350
4053
6740
855.5692369
1160

Blank
DUP
DUP

536
508
224.6368067
398.6676483
251.2287758
604.1638876
816.5818081
246.4612913

TN
µg/L
694.2821739
1014.359308
159.3250521
650.12676
519.5445012
397.5530834
1134.018126
475.2041344
1069.67208
2834.695063
111.5100218
2051.64984
272.199
3281.35164
3194.54916
1940.218062
743.65
5797.249324
619.40892
2427.79392
2479.662271
1530.83496
2263.99101
2500.12932

TP
µg/L
188.477
155.73
15.9255
105.35
92.755
134.3185
76.3815
130.54
100.312
143.135
ND
144.3945
153.211
120.464
189.7365
174.6225
14.666
316.946
143.135
124.2425
146.9135
70.084
190.996
73.8625

3122.21376
1852.268658
978.04724
23.4055803
1187.863565
706.2114

202.3315
120.464
192.2555
ND
92.755
253.971

5938.58406
777.55
1646.5716
5185.537759
3599.6274
1111.28964

134.3185
208.629
177.1415
194.7745
63.7865
253.35

Table A.13: 2012 Melt Results. ND = Non Detect. DUP = Duplicate Sample

110

ID

Lat

Long

Date

MR11
MR13
MR14
MR15
MR16
MR17
MR2
MR20
MR20
MR20
MR21
MR23
MR24
MR3
MR4
MR5
MR6
MR7
MR8
MR9
SB12
SB12
SB12
SB13
SB15
SB16
SB17
SB19
SB19
SB19
SB2
SB20
SB21
SB22
SB24
SB27
SB28
SB29

43.306
43.699
43.540
44.336
43.490
43.608
44.082
44.201
44.201
44.201
43.901
43.435
43.289
44.139
44.006
43.853
44.111
43.934
43.779
44.409
43.657
43.657
43.657
43.211
43.250
43.450
43.872
43.015
43.015
43.015
43.302
43.218
43.160
43.111
43.038
43.328
43.867
43.627

-86.114
-85.460
-85.305
-84.922
-85.449
-85.481
-85.014
-85.053
-85.053
-85.053
-85.255
-85.666
-86.223
-84.899
-85.111
-85.444
-85.028
-85.129
-85.500
-84.805
-84.267
-84.267
-84.267
-83.315
-83.873
-83.441
-84.321
-84.180
-84.180
-84.180
-84.144
-84.106
-83.351
-83.519
-83.774
-83.758
-84.545
-84.708

3/8/12
3/8/12
3/8/12
3/9/12
3/8/12
3/8/12
3/9/12
3/9/12
3/9/12
3/9/12
3/9/12
3/8/12
3/8/12
3/9/12
3/9/12
3/8/12
3/9/12
3/9/12
3/8/12
3/9/12
3/6/12
3/6/12
3/6/12
3/7/12
3/7/12
3/7/12
3/5/12
3/7/12
3/7/12
3/7/12
3/6/12
3/7/12
3/7/12
3/7/12
3/7/12
3/7/12
3/5/12
3/6/12

Sample
Type

Flow
cfs
80.6940895
147.0857255
260.6931154
176.926647
108.2395555
124.3783734
659.1841902
455

Blank
DUP
2820
3570
43
697.9244161
154.6430713
698.9132277
321.8228611
152.0650982
856.734622
159.975591
116.0087895
Blank
DUP
334.077062
1881.178754
1041.713021
446.4484374
956
Blank
DUP
279.2333329
911.2252041
275
634
1520
1470
308
538

TN
µg/L
956.64
759.0642832
1481.998018
1239.831136
1185.40872
927.79
1185.40872
632.1800677
124.664293
1204.69816
663.9948998
1239.83
776.1831175
678.9696771
1185.40872
1296.38444
813.0378375
654.94912
1771.465534
568.0221807
1110.988734
93.28865173
1545.30204
2711.0347
1035.91556
6061.683402
415.6869564
879.81636
45.69708309
807.48096
7409.074632
1484.646015
1386.087721
1093.78388
1002.15904
5071.169951
1146.82984
1371.69708

TP
µg/L
203.591
13.4065
132.43
36.0775
121.7235
86.4575
33.5585
23.4825
ND
70.084
124.2425
206.11
79.67
97.793
ND
57.489
168.325
33.56
226.262
3.3305
83.31
ND
144.4
219.9645
228.781
56.2295
139.3565
136.8375
ND
149.4325
141.8755
219.9645
194.7745
143.135
129.28
94.0145
172.1035
337.098

Table A.14: 2012 Melt Results (continued). ND = Non Detect. DUP = Duplicate Sample

111

ID

Lat

Long

Date

SB3
SB3’
SB30
SB31
SB32
SB33
SB4
SB5
SB6
SB7
SB8
SB9
SW-12
SW-13
SW-14
SW-15
SW-20
SW-21
SW-21
SW-21
SW-36
SW-36
SW-39
SW-4
SW-5
SW-5
SW-5
SW-62
SW-62
SW-62
SW-64
SW-65
SW-66

43.675
43.674
43.379
43.563
43.610
43.409
43.597
43.397
43.292
43.741
43.298
43.380
44.623
44.715
44.765
44.900
44.050
44.072
44.072
44.072
42.907
42.907
42.949
42.965
43.318
43.318
43.318
45.100
45.100
45.100
45.214
45.365
45.746

-84.388
-84.390
-84.656
-84.370
-84.242
-83.964
-84.944
-84.836
-83.992
-85.127
-84.143
-84.094
-86.223
-86.125
-85.621
-85.411
-83.688
-84.019
-84.019
-84.019
-85.781
-85.781
-85.849
-85.675
-86.038
-86.038
-86.038
-85.097
-85.097
-85.097
-85.013
-84.962
-84.828

3/5/12
3/5/12
3/6/12
3/5/12
3/5/12
3/6/12
3/4/12
3/4/12
3/7/12
3/4/12
3/6/12
3/6/12
3/11/12
3/11/12
3/10/12
3/10/12
3/5/12
3/5/12
3/5/12
3/5/12
3/2/12
3/2/12
3/2/12
3/2/12
3/8/12
3/8/12
3/8/12
3/10/12
3/10/12
3/10/12
3/10/12
3/10/12
3/10/12

Sample
Type

Flow
cfs
367.5553976
921
4660
9590
362.0109897
437.4785036
136.314742
176.7500735
324.4714636
298.3385856
391.5340789
211.9235147
525.3061625
912.1433863
444.5767583
297

Blank
DUP
158.0332825
Blank
217.0794609
8540
2955.451928
Blank
DUP
187
Blank
DUP
113.5720752
214.5721172
92.0654229

TN
µg/L
793.69
3831.10068
3388.374711
4950.48
2456.72808
3570.69324
867.13296
807.48096
3382.6212
1355.52
2592.768879
3174.7521
1159.593125
190.4448152
1493.57
496.8799603
2644.80012
698.35036
61.06631888
823.73172
1125.75672
79.80486507
1763.17
1487.43372
663.9948998
46.72766878
732.88
1151.6522
111.5100218
1267.38884
867.13296
1109.007188
593.9317825

TP
µg/L
180.92
70.084
56.2295
150.692
116.6855
78.9005
153.211
131.7995
126.7615
76.3815
78.91
116.6855
ND
253.971
223.743
177.14
189.7365
174.6225
ND
207.3695
168.96
ND
116.6855
87.717
260.2685
ND
230.0405
144.3945
ND
203.591
178.401
49.932
51.1915

Table A.15: 2012 Melt Results (continued). ND = Non Detect. DUP = Duplicate Sample

112

Appendix B
Nitrate Results

113

ID

Date

Lat

Long

J-201
J-202
J-203
SW-1
SW-10
SW-101
SW-102
SW-103
SW-104
SW-11
SW-11
SW-11
SW-12
SW-13
SW-14
SW-15
SW-15
SW-15
SW-16
SW-17
SW-18
SW-19
SW-2
SW-20

3/19/11
3/19/11
3/19/11
3/5/11
3/20/11
3/7/11
3/6/11
3/22/11
3/6/11
3/20/11
3/20/11
3/20/11
3/20/11
3/20/11
3/19/11
3/19/11
3/19/11
3/19/11
3/23/11
3/23/11
3/22/11
3/22/11
3/5/11
3/22/11

45.13562
45.032895
45.03636
42.087101
44.258801
43.3274
43.3102
43.245098
43.5891
44.356201
44.356201
44.356201
44.6231
44.715099
44.7649
44.899799
44.899799
44.899799
45.595299
45.559101
45.067799
44.415401
42.183998
44.049599

-85.080693
-85.063817
-84.968732
-86.476194
-85.939177
-83.745823
-83.955886
-84.100284
-84.221356
-86.049494
-86.049494
-86.049494
-86.223457
-86.125244
-85.620741
-85.410914
-85.410914
-85.410914
-84.474136
-84.403295
-83.436025
-83.330781
-86.371158
-83.687858

Sample Type

Blank
DUP

Blank
DUP

NOx
µg/L
427.767458
1851.246197
1667.68841
4050.106798
297.5626204
5173.973871
2756.87437
3431.876158
2131.036822
290.3843768
0
208.1199606
134.5879328
267.7311532
880.204067
215.0131516
0
214.1276146
54.0583336
19.9472548
319.615808
199.5583534
1114.479413
1786.539089

Table B.1: 2011 Melt Results. ND = Non Detect. DUP = Duplicate Sample

114

ID

Date

Lat

Long

SW-21
SW-21
SW-21
SW-23
SW-24
SW-25
SW-26
SW-27
SW-28
SW-29
SW-3
SW-3
SW-3
SW-31
SW-32
SW-33
SW-34
SW-35
SW-36
SW-37
SW-39
SW-4
SW-40
SW-40
SW-40

3/22/11
3/22/11
3/22/11
3/7/11
3/7/11
3/22/11
3/4/11
3/4/11
3/4/11
3/4/11
3/6/11
3/6/11
3/6/11
3/5/11
3/5/11
3/5/11
3/5/11
3/5/11
3/6/11
3/6/11
3/6/11
3/6/11
3/6/11
3/6/11
3/6/11

44.072399
44.072399
44.072399
43.064301
42.858002
42.749001
42.5956
42.297001
42.064201
41.963902
42.645599
42.645599
42.645599
42.417099
42.420502
42.798698
42.796902
42.916199
42.9067
42.903801
42.949299
42.964901
42.994202
42.994202
42.994202

-84.019051
-84.019051
-84.019051
-82.596818
-82.538059
-82.489671
-82.90946
-83.191016
-83.252775
-83.546067
-86.11776
-86.11776
-86.11776
-86.240863
-86.244757
-86.059097
-86.144102
-86.146693
-85.781123
-85.772215
-85.849483
-85.675129
-86.017657
-86.017657
-86.017657

Sample Type

Blank
DUP

Blank
DUP

Blank
DUP

NOx
µg/L
582.329744
0
574.040537
5246.979864
840.519428
1830.196192
2843.695708
2150.059286
1725.959648
6870.21722
1813.656845
0
2174.57366
1024.512725
1463.308406
6415.757696
597.35
2878.826848
3312.82708
1829.74652
3119.809048
736.185749
2712.64591
0
506.1554223

Table B.2: 2011 Melt Results (continued). ND = Non Detect. DUP = Duplicate Sample

115

ID

Date

Lat

Long

SW-41
SW-43
SW-48
SW-51
SW-52
SW-55
SW-56
SW-57
SW-59
SW-6
SW-6
SW-6
SW-60
SW-62
SW-63
SW-64
SW-65
SW-66
SW-67
SW-69
SW-7
SW-70
SW-71
SW-73
SW-8
SW-9
SW-91
SW-94
SW-97

3/6/11
3/21/11
3/21/11
3/21/11
3/21/11
3/21/11
3/20/11
3/20/11
3/19/11
3/21/11
3/21/11
3/21/11
3/19/11
3/19/11
3/19/11
3/23/11
3/23/11
3/23/11
3/23/11
3/23/11
3/20/11
3/23/11
3/23/11
3/22/11
3/20/11
3/20/11
3/4/11
3/4/11
3/4/11

42.962399
43.1973
43.415699
43.490299
43.558998
43.904999
43.999401
44.9035
45.016899
43.462601
43.462601
43.462601
44.745098
45.099899
45.179501
45.2141
45.3647
45.745998
45.488098
45.431099
43.943802
45.337299
45.1623
44.259602
44.030201
44.208199
42.588402
42.0951
41.931301

-86.203368
-86.186454
-86.276515
-86.437662
-86.50317
-86.308022
-86.37208
-85.975496
-85.617735
-86.232338
-86.232338
-86.232338
-85.54849
-85.096655
-85.16416
-85.013231
-84.961607
-84.827748
-84.07552
-83.839247
-86.278413
-83.684801
-83.410055
-83.52734
-86.50636
-86.246047
-82.903226
-83.211266
-83.34023

Sample Type

Blank
DUP

NOx
µg/L
629.6350532
947.345192
146.3274533
1818.582592
1385.594066
2330.333634
694.004186
1046.19
430.965114
350.561
0
452.8562568
337.1871848
1446.991298
283.0044534
293.0023468
130.8706485
37.7721952
453.9696752
135.0776396
335.2380356
331.8488
ND
386.9901728
221.8486875
215.6857726
1812.766291
ND
4877.812554

Table B.3: 2011 Melt Results (continued). ND = Non Detect. DUP = Duplicate Sample

116

ID

Date

Lat

Long

BC10
BC11
BC12
BC13
BC14
BC15
BC17
BC9
GR1
GR11
GR11
GR11
GR12
GR13
GR14
GR15
GR17
GR20
GR21
GR23
GR24
GR25
GR26
GR27
GR28
GR3
GR30
GR30
GR30
GR31
GR5
GR6
GR7
GR8
GR9
MR10
MR11

3/9/12
3/9/12
3/10/12
3/10/12
3/10/12
3/11/12
3/10/12
3/9/12
3/2/12
3/4/12
3/4/12
3/4/12
3/2/12
3/4/12
3/3/12
3/2/12
3/3/12
3/3/12
3/3/12
3/7/12
3/7/12
3/2/12
3/3/12
3/4/12
3/3/12
3/3/12
3/3/12
3/3/12
3/3/12
3/2/12
3/4/12
3/2/12
3/3/12
3/4/12
3/4/12
3/8/12
3/8/12

45.50764
45.44003
45.165655
44.845585
44.69917
44.52819
44.65679
45.65927
43.12973
43.01474
43.01474
43.01474
42.96558
43.06736
42.77521
43.0631
42.312
42.28331
42.53555
42.727659
42.750321
42.85708
42.82789
43.10872
42.97123
42.60863
42.61599
42.61599
42.61599
43.08217
43.09736
43.22578
42.87736
42.99961
43.14473
43.40164
43.30589

-84.76265
-84.78534
-85.239769
-85.307429
-85.30557
-85.94865
-85.43884
-84.49263
-85.17377
-84.42586
-84.42586
-84.42586
-85.35174
-84.84189
-85.35897
-86.06651
-84.38722
-84.41029
-84.62442
-84.477961
-84.555337
-84.91265
-84.75898
-84.69332
-85.06897
-85.09383
-85.23645
-85.23645
-85.23645
-85.59049
-84.70956
-85.30203
-84.38425
-84.90849
-84.50335
-85.5628
-86.11443

Sample Type

Blank
DUP

Blank
DUP

NOx
µg/L
349.8131862
345.5550495
201.3527037
320.5209196
462.5068334
100.1715435
443.6154012
158.4701182
340.6451801
2790.784189
ND
616.3103177
519.2301768
1873.358544
2543.062047
1700.159858
676.139703
5572.808401
731.2201456
1991.868993
1709.356954
1263.441766
1377.895405
1222.814129
592.1609789
1716.115662
1168.266882
12.75670951
1156.519452
837.8317027
10123.03868
459.4064666
589.1411178
7988.841806
2264.088764
1207.829928
626.7997862

Table B.4: 2012 Melt Results. ND = Non Detect. DUP = Duplicate Sample

117

ID

Date

Lat

Long

MR13
MR14
MR15
MR16
MR17
MR2
MR20
MR20
MR20
MR21
MR23
MR24
MR3
MR4
MR5
MR6
MR7
MR8
MR9
SB12
SB12
SB12
SB13
SB15
SB16
SB17
SB19
SB19
SB19
SB2
SB20
SB21
SB22
SB24
SB27
SB28
SB29
SB3

3/8/12
3/8/12
3/9/12
3/8/12
3/8/12
3/9/12
3/9/12
3/9/12
3/9/12
3/9/12
3/8/12
3/8/12
3/9/12
3/9/12
3/8/12
3/9/12
3/9/12
3/8/12
3/9/12
3/6/12
3/6/12
3/6/12
3/7/12
3/7/12
3/7/12
3/5/12
3/7/12
3/7/12
3/7/12
3/6/12
3/7/12
3/7/12
3/7/12
3/7/12
3/7/12
3/5/12
3/6/12
3/5/12

43.69864
43.5399
44.33609
43.48959
43.60761
44.08187
44.20086
44.20086
44.20086
43.90098
43.43517
43.28875
44.13933
44.00591
43.85285
44.111408
43.93358
43.77874
44.40864
43.65714
43.65714
43.65714
43.21095
43.24987
43.45034
43.871818
43.01495
43.01495
43.01495
43.30233
43.21777
43.15983
43.11134
43.03844
43.32754
43.86695
43.62684
43.67512

-85.46043
-85.30548
-84.9216
-85.4486
-85.480685
-85.01359
-85.05293
-85.05293
-85.05293
-85.25467
-85.66613
-86.22311
-84.8994
-85.11092
-85.444458
-85.028132
-85.12883
-85.49961
-84.80501
-84.26715
-84.26715
-84.26715
-83.3146
-83.87269
-83.44139
-84.321231
-84.17964
-84.17964
-84.17964
-84.14385
-84.10555
-83.35131
-83.51899
-83.77401
-83.75812
-84.54549
-84.70802
-84.38839

Sample Type

Blank
DUP

Blank
DUP

Blank
DUP

NOx
µg/L
266.5138838
136.9166199
290.2351194
497.5192703
219.4142962
337.8536637
450.5869143
25.08646871
573.8087129
427.8774506
838.8409191
381.978334
225.0605142
354.6508668
457.2180539
594.0207878
422.8037192
907.7956489
50.40907302
1415.21538
ND
1315.1161
2717.858351
1041.177877
4542.886475
288.2179637
934.257909
ND
588.5110585
7527.484258
1453.012232
567.9800556
596.9733799
762.71316
6501.721677
542.4893401
619.5699381
468.4208956

Table B.5: 2012 Melt Results (continued). ND = Non Detect. DUP = Duplicate Sample

118

ID

Date

Lat

Long

SB3’
SB30
SB31
SB32
SB33
SB4
SB5
SB6
SB7
SB8
SB9
SW-12
SW-13
SW-14
SW-15
SW-20
SW-21
SW-21
SW-21
SW-36
SW-36
SW-36
SW-39
SW-4
SW-5
SW-5
SW-62
SW-62
SW-62
SW-64
SW-65
SW-66

3/5/12
3/6/12
3/5/12
3/5/12
3/6/12
3/4/12
3/4/12
3/7/12
3/4/12
3/6/12
3/6/12
3/11/12
3/11/12
3/10/12
3/10/12
3/5/12
3/5/12
3/5/12
3/5/12
3/2/12
3/2/12
3/2/12
3/2/12
3/2/12
3/8/12
3/8/12
3/10/12
3/10/12
3/10/12
3/10/12
3/10/12
3/10/12

43.674173
43.378977
43.5634
43.61029
43.40875
43.59702
43.3973
43.292201
43.74069
43.29845
43.37951
44.6231
44.715099
44.7649
44.899799
44.049599
44.072399
44.072399
44.072399
42.9067
42.9067
42.9067
42.949299
42.964901
43.3181
43.3181
45.099899
45.099899
45.099899
45.2141
45.3647
45.745998

-84.389555
-84.65597
-84.36966
-84.24164
-83.96436
-84.94359
-84.83624
-83.992242
-85.12742
-84.14273
-84.094
-86.223457
-86.125244
-85.620741
-85.410914
-83.687858
-84.019051
-84.019051
-84.019051
-85.781123
-85.781123
-85.781123
-85.849483
-85.675129
-86.03813
-86.03813
-85.096655
-85.096655
-85.096655
-85.013231
-84.961607
-84.827748

Sample Type

Blank
DUP
Blank
DUP

DUP
Blank
DUP

NOx
µg/L
5578.086968
3237.105255
7826.53935
1280.02009
2429.042143
326.1504248
554.8123137
2148.330772
277.0606123
2624.742406
6291.419834
233.4763823
193.038179
370.0861244
352.8817201
2024.657205
281.4512779
ND
278.917344
908.7460593
ND
985.1694554
589.0212381
628.9688034
590.7124537
825.3481123
529.6767952
ND
531.4229666
307.8048816
281.0861289
147.6876311

Table B.6: 2012 Melt Results (continued). ND = Non Detect. DUP = Duplicate Sample

119

Appendix C
Isotope Results

120

ID

Date

Lat

Long

BC10
BC11
BC12
BC13
BC14
BC15
BC17
BC9
GR1
GR11
GR11
GR12
GR13
GR14
GR15
GR17
GR20
GR21
GR23
GR24
GR25
GR26
GR27
GR28
GR3
GR30
GR30
GR31
GR5
GR6
GR7
GR8
GR9
MR10
MR11
MR13
MR14
MR15

3/9/12
3/9/12
3/10/12
3/10/12
3/10/12
3/11/12
3/10/12
3/9/12
3/2/12
3/4/12
3/4/12
3/2/12
3/4/12
3/3/12
3/2/12
3/3/12
3/3/12
3/3/12
3/7/12
3/7/12
3/2/12
3/3/12
3/4/12
3/3/12
3/3/12
3/3/12
3/3/12
3/2/12
3/4/12
3/2/12
3/3/12
3/4/12
3/4/12
3/8/12
3/8/12
3/8/12
3/8/12
3/9/12

45.50764
45.44003
45.165655
44.845585
44.69917
44.52819
44.65679
45.65927
43.12973
43.01474
43.01474
42.96558
43.06736
42.77521
43.0631
42.312
42.28331
42.53555
42.727659
42.750321
42.85708
42.82789
43.10872
42.97123
42.60863
42.61599
42.61599
43.08217
43.09736
43.22578
42.87736
42.99961
43.14473
43.40164
43.30589
43.69864
43.5399
44.33609

-84.76265
-84.78534
-85.239769
-85.307429
-85.30557
-85.94865
-85.43884
-84.49263
-85.17377
-84.42586
-84.42586
-85.35174
-84.84189
-85.35897
-86.06651
-84.38722
-84.41029
-84.62442
-84.477961
-84.555337
-84.91265
-84.75898
-84.69332
-85.06897
-85.09383
-85.23645
-85.23645
-85.59049
-84.70956
-85.30203
-84.38425
-84.90849
-84.50335
-85.5628
-86.11443
-85.46043
-85.30548
-84.9216

Sample
Type

DUP

DUP

δ 15 NN O3

δ 18 ON O3

2.9
4.14
3.4
4.56
3.13
1.45
4.03
4.13
6.24
6.87
6.83
5.18
6.04
5.36
8.82
6.95
5.28
5.19
5.87
6.93
6.49
6.42
7.79
7
4.8
6.39
7.16
7.03
5.66
5.07
3.21
5.33
6.78
5.49
6.16
3.29
3.64
0.74

10.96
14.42
7
-0.53
7.31
8.95
5.68
19.49
4.28
0.49
0.41
3.98
-0.07
-0.08
2.37
2.67
-0.96
1.39
1.67
1.65
0.78
2.19
0.96
0.73
0.58
2.23
3.01
3.24
-1.58
3.61
1.4
-0.28
0.2
3.18
4.45
17.92
8.79
28.14

Table C.1: 2012 Melt Results. ND = Non Detect. DUP = Duplicate Sample

121

ID

Date

Lat

Long

Sample Type

MR20
MR20
MR21
MR23
MR24
MR3
MR4
MR5
MR6
MR7
MR8
MR9
SB12
SB12
SB13
SB15
SB16
SB17
SB19
SB19
SB2
SB20
SB21
SB22
SB24
SB27
SB28
SB29
SB3
SB3’
SB30
SB31
SB32
SB33
SB4
SB5
SB6
SB7
SB8

3/9/12
3/9/12
3/9/12
3/8/12
3/8/12
3/9/12
3/9/12
3/8/12
3/9/12
3/9/12
3/8/12
3/9/12
3/6/12
3/6/12
3/7/12
3/7/12
3/7/12
3/5/12
3/7/12
3/7/12
3/6/12
3/7/12
3/7/12
3/7/12
3/7/12
3/7/12
3/5/12
3/6/12
3/5/12
3/5/12
3/6/12
3/5/12
3/5/12
3/6/12
3/4/12
3/4/12
3/7/12
3/4/12
3/6/12

44.20086
44.20086
43.90098
43.43517
43.28875
44.13933
44.00591
43.85285
44.111408
43.93358
43.77874
44.40864
43.65714
43.65714
43.21095
43.24987
43.45034
43.871818
43.01495
43.01495
43.30233
43.21777
43.15983
43.11134
43.03844
43.32754
43.86695
43.62684
43.67512
43.674173
43.378977
43.5634
43.61029
43.40875
43.59702
43.3973
43.292201
43.74069
43.29845

-85.05293
-85.05293
-85.25467
-85.66613
-86.22311
-84.8994
-85.11092
-85.444458
-85.028132
-85.12883
-85.49961
-84.80501
-84.26715
-84.26715
-83.3146
-83.87269
-83.44139
-84.321231
-84.17964
-84.17964
-84.14385
-84.10555
-83.35131
-83.51899
-83.77401
-83.75812
-84.54549
-84.70802
-84.38839
-84.389555
-84.65597
-84.36966
-84.24164
-83.96436
-84.94359
-84.83624
-83.992242
-85.12742
-84.14273

DUP

DUP

DUP

δ 15 NN O3

δ 18 ON O3

3.89
3.85
3.96
5.54
6.16
3.06
2.4
1.14
5.06
1.53
3.93
-1.17
6.11
5.34
5.44
5.23
3.12
4.87
5.98
7.4
3.99
5.29
8.4
4.7
4.39
3.2
7.81
3.86
6.02
4.94
4.9
5.72
5.25
6.25
5.63
5.4
5.72
6.09
3.1

14.6
14.89
11.22
2.61
12.41
17.62
12.18
16.97
8.88
14.63
8.28
61.31
0.6
0.34
-0.01
0.57
-1.42
7.11
1.66
1.82
-1.75
1.36
3.16
2.54
2.54
0.95
-1.34
2.44
-0.62
3.17
0.8
1.06
0.1
-0.76
4.32
2.04
-0.7
7.57
-0.26

Table C.2: 2012 Melt Results (continued). ND = Non Detect. DUP = Duplicate Sample

122

ID

Date

Lat

Long

SW-15
SW-20
SW-21
SW-21
SW-36
SW-36
SW-39
SW-4
SW-5
SW-5
SW-62
SW-62
SW-64
SW-65
SW-66
DR11
DR12
DR13

3/10/12
3/5/12
3/5/12
3/5/12
3/2/12
3/2/12
3/2/12
3/2/12
3/8/12
3/8/12
3/10/12
3/10/12
3/10/12
3/10/12
3/10/12

44.899799
44.049599
44.072399
44.072399
42.9067
42.9067
42.949299
42.964901
43.3181
43.3181
45.099899
45.099899
45.2141
45.3647
45.745998

-85.410914
-83.687858
-84.019051
-84.019051
-85.781123
-85.781123
-85.849483
-85.675129
-86.03813
-86.03813
-85.096655
-85.096655
-85.013231
-84.961607
-84.827748

Sample Type

DUP
DUP

DUP
DUP

Standard
Standard
Standard

δ 15 NN O3

δ 18 ON O3

2.09
1.32
6.9
7.1
6.97
6.81
7.13
6.09
5.82
6.14
1.79
0.98
3.24
1.85
2.44
3.94
4.16
4.58

6.59
-0.58
4.11
3.83
2.49
2.59
2.15
0.95
3.46
3.8
-2.12
-2.33
6.03
14.71
18.58
24.87
24.95
25.43

Table C.3: 2012 Melt Results (continued). ND = Non Detect. DUP = Duplicate Sample

123

BIBLIOGRAPHY

124

BIBLIOGRAPHY

Alexander, R. B., R. A. Smith, and G. E. Schwarz, Effect of stream channel size on the
delivery of nitrogen to the Gulf of Mexico, Nature, 403 (6771), 758–761, 2000.
Alexander, R. B., P. J. Johnes, E. W. Boyer, and R. A. Smith, A comparison of models for
estimating the riverine export of nitrogen from large watersheds, Biogeochemistry, 57/58,
295–339, 2002.
Anderson, D. M., P. M. Gilber, and J. M. Burkholder, Harmful Algal Blooms and Eutrophication Nutrient Sources , Composition , and Consequences, Estuaries, 25 (4), 704–726,
2002.
Anisfeld, S. C., R. T. Barnes, M. a. Altabet, and T. Wu, Isotopic apportionment of atmospheric and sewage nitrogen sources in two Connecticut rivers., Environmental science &
technology, 41 (18), 6363–9, 2007.
Arnold, J., R. Srinivasan, R. Muttiah, and J. Williams, Large Area Hydrologic Modeling
and Assessment Part I: Model Development, Journal of The American Water Resources
Association, 34 (1), 73–89, 1998.
Auer, M. T., L. M. Tomlinson, S. N. Higgins, S. Y. Malkin, E. T. Howell, and H. a. Bootsma,
Great Lakes Cladophora in the 21st century: same algae-different ecosystem, Journal of
Great Lakes Research, 36 (2), 248–255, doi:10.1016/j.jglr.2010.03.001, 2010.
Beaulac, M., and K. H. Reckhow, An Examination of Land Use-Nutrient Export Relationships, 18 (6), 1013–1024, 1982.
Beyer, H., Geospatial Modelling Environment, 2012.
Bjerklie, D. M., S. Lawrence Dingman, C. J. Vorosmarty, C. H. Bolster, and R. G. Congalton,
Evaluating the potential for measuring river discharge from space, Journal of Hydrology,
278 (1-4), 17–38, doi:10.1016/S0022-1694(03)00129-X, 2003.
Borah, D. K., and M. Bera, Watershed-Scale Hydrologic and Nonpoint-Source Pollution
Models: Review of Applications, Transaction of the ASAE, 47 (3), 789–804, 2004.
Bottcher, J., O. Strebel, S. Voerkelius, and H. Schmidt, Using Isotope Fractionation of
Nitrate-Nitrogen and Nitrate-Oxygen for Evaluation of Microbial Denitrification in a
Sandy Aquifer, Journal of Hydrolo, 114, 413–424, 1990.
Boyer, E. W., C. L. Goodale, N. Jaworski, and R. W. Howarth, Anthropogenic nitrogen
sources and relationships to riverine nitrogen export in the northeastern U . S . A .,
Biogeochemistry, 57/58, 137–169, 2002.
Burns, D. a., E. W. Boyer, E. M. Elliott, and C. Kendall, Sources and transformations of
nitrate from streams draining varying land uses: evidence from dual isotope analysis.,
Journal of environmental quality, 38 (3), 1149–59, doi:10.2134/jeq2008.0371, 2009.
125

Carpenter, S. R., N. F. Caraco, D. L. Correll, R. W. Howarth, a. N. Sharpley, and V. H.
Smith, Nonpoint Pollution of Surface Waters With Phosphorus and Nitrogen, Ecological
Applications, 8 (3), 559–568, doi:10.1890/1051-0761(1998)008[0559:NPOSWW]2.0.CO;2,
1998.
Casciotti, K. L., D. M. Sigman, M. G. Hastings, J. K. B¨hlke, and a. Hilkert, Measurement of
o
the oxygen isotopic composition of nitrate in seawater and freshwater using the denitrifier
method., Analytical chemistry, 74 (19), 4905–12, 2002.
Chang, C. C. Y., C. Kendall, S. R. Silva, W. A. Battaglin, and D. H. Campbell, Nitrate
stable isotopes : tools for determining nitrate sources among different land uses in the
Mississippi River Basin, Water Resources, 1885 (2002), 1874–1885, doi:10.1139/F02-153,
2003.
Cheng, J., Challenges of CAFO Waste Management, Journal of Environmental Engineering,
129 (5), 391, doi:10.1061/(ASCE)0733-9372(2003)129:5(391), 2003.
Correll, D. L., The Role of Phosphorus in the Eutrophication of Receiving Waters: A Review,
Journal of Environment Quality, 27 (2), 261, doi:10.2134/jeq1998.00472425002700020004x,
1998.
Crumpton, W. G., T. M. Isenhart, and P. D. Mitchell, Nitrate and Organic N Analyses with
Second-Derivative Spectroscopy, Liminology and Oceanography, 37 (4), 907–913, 1992.
Davidson, E. A., M. B. David, J. N. Galloway, C. L. Goodale, R. Haeuber, J. A. Harrison,
R. W. Howarth, D. B. Jaynes, R. R. Lowrance, B. T. Nolan, J. L. Peel, R. W. Pinder,
E. Porter, C. S. Snyder, A. R. Townsend, and M. H. Ward, Excess Nitrogen in the U.S.
Environment: Trends, Risks, and Solutions, Issues in Ecology, (15), 2012.
Destouni, G., G. A. Lindgren, and I.-M. Gren, Policy Analysis Effects of Inland Nitrogen
Transport and Attenuation Modeling on Coastal Nitrogen Load Abatement, 40 (20), 2006.
Dolan, D. M., Point Source Loadings of Phosphorus to Lake Erie : 1986-1990, Journal of
Great Lakes Research, 19 (2), 212–223, doi:10.1016/S0380-1330(93)71212-5, 1993.
Dolan, D. M., and S. C. Chapra, Great Lakes total phosphorus revisited: 1. Loading
analysis and update (1994-2008), Journal of Great Lakes Research, 38 (4), 730–740, doi:
10.1016/j.jglr.2012.10.001, 2012.
Dolan, D. M., and K. P. McGunagle, Lake Erie Total Phosphorus Loading Analysis and
Update: 1996-2002, Journal of Great Lakes Research, 31, 11–22, doi:10.1016/S03801330(05)70301-4, 2005.
Dolan, D. M., A. K. Yui, and R. D. Geist, Evaluation of River Load Estimation Methods
for Total Phosphorus, Journal of Great Lakes Research, 7 (3), 207–214, 1981.
Evans, B. M., D. W. Lehning, K. J. Corradini, G. W. Petersen, E. Nizeyimana, J. M.
Hamlett, P. D. Robillard, and R. L. Day, A Comprehensive GIS-Based Modeling Approach
for Prediciting Nutrient Loads in Watersheds, Journal of Spatial Hydrology, 2 (2), 2002.
126

Fry, J., G. Xian, S. Jin, J. Dewitz, C. Homer, L. Yang, C. Barnes, N. Herold, and J. Wickham, Completion of the 2006 National Land Cover Database for the Conterminous United
States, Photogrammetric Engineering & Remote Sensing, 77 (9), 858–864, 2011.
Georgas, N., S. Rangarajan, K. J. Farley, and S. C. K. Jagupilla, AVGWLF-Based Estimation of Nonpoint Source Nitrogen Loads Generated within Long Island Sound Watersheds, Journal Of The American Water Resources Association, 45 (3), doi:10.1111/j.17521688.2009.00318.x, 2009.
Gesch, D., The National Elevation Dataset, in Digital Elevation Model Technologies and
Applications: The DEM Users Manual, edited by D. Maune, 2nd ed., pp. 99–118, American
Society for Photogrammetry and Remote Sensing, Bethesda, Maryland, 2007.
Gesch, D., M. Oimoen, S. Greenlee, C. Nelson, M. Steuck, and D. Tyler, The National
Elevation Dataset, Photogrammetric Engineering & Remote Sensing, 68 (1), 5–11, 2002.
Gilliom, R. J., and C. R. Patmont, Lake phosphorus loading from septic systems by seasonally perched groundwater, Journal (Water Pollution Control Federation), 55 (10), 1297–
1305, 1983.
Grizzetti, B., F. Bouraoui, and G. Demarsily, Modelling nitrogen pressure in river
basins: A comparison between a statistical approach and the physically-based SWAT
model, Physics and Chemistry of the Earth, Parts A/B/C, 30 (8-10), 508–517, doi:
10.1016/j.pce.2005.07.005, 2005.
Haith, D. a., and L. L. Shoenaker, Generalized Watershed Loading Functions for Stream
Flow Nutrients, Journal of the American Water Resources Association, 23 (3), 471–478,
doi:10.1111/j.1752-1688.1987.tb00825.x, 1987.
Hatfield, J. L., and D. R. Keeney, Chapter 1 . The Nitrogen Cycle , Historical Perspective ,
and, Perspective, pp. 1–18, 2008.
Hinderer, J. M., and M. W. Murray, Feast and Faimine in the Great Lakes-How Nutrients
and Invasive Species Interact to Overwhelm the Coasts and Starve Offshore Waters, Tech.
rep., National Wildlife Federation, 2011.
Howarth, R. W., G. Billen, D. Swaney, A. Townsend, N. Jaworski, K. Lajtha, J. A. Downing,
R. Elmgren, N. F. Caraco, T. E. Jordan, F. Berendse, J. Freney, V. Kudeyarov, P. Murdoch, and Z. Zhao-Lian, Regional nitrogen budgets and riverine N & P fluxes for the
drainages to the Norh Atlantic Ocean: Natural and human influences, Biogeochemistry,
35, 75–139, 1996.
Hyndman, D., A. Kendall, and R. Welty, Evaluating Temporal and Spatial Variations in
Recharge and Streamflow Using the Integrated Landscape Hydology Model (ILHM), in
Geophysical Monograph Series, vol. 171, pp. 121–141, 2007.
Jaworski, N. A., P. M. Groffman, A. A. Keller, and J. C. Prager, A Watershed Nitrogen
and Phosphorus Balance: The Upper Potomac River Basin, Estuaries, 15 (1), 83–95, doi:
10.2307/1352713, 1992.
127

Johnes, P., B. Moss, and G. Phillips, The determination of total nitrogen and total phosphorus concentrations in freshwaters from land use, stock headage and population data:
testing of a model for use in conservation and water quality management, Freshwater
Biology, 36 (2), 451–473, doi:10.1046/j.1365-2427.1996.00099.x, 1996.
Johnes, P. J., Evaluation and management of the impact of land use change on the nitrogen
and phosphorus load delivered to surface waters: the export coefficient modelling approach,
Journal of Hydrology, 183 (3-4), 323–349, doi:10.1016/0022-1694(95)02951-6, 1996.
Johnes, P. J., and A. L. Heathwaite, Modelling the Impact of Land Use Change on Water Quality in Agricultural Catchments, Hydrological Processes, 11 (3), 269–286, doi:
10.1002/(SICI)1099-1085(19970315)11:3¡269::AID-HYP442¿3.0.CO;2-K, 1997.
Kaushal, S. S., P. M. Groffman, L. E. Band, E. M. Elliott, C. a. Shields, and C. Kendall,
Tracking nonpoint source nitrogen pollution in human-impacted watersheds., Environmental science & technology, 45 (19), 8225–32, doi:10.1021/es200779e, 2011.
Kellogg, R. L., C. H. Lander, D. C. Moffitt, and N. Gollehon, Manure Nutrients Relative to
the Capacity of Cropland and Pastureland to Assimilate Nutrients: Spatial and Temporal
Trends for the United States, Tech. rep., NRCS, 2000.
Kendall, A. D., and D. W. Hyndman, Simulating Spatial and Temporal Variability of Regional Evapotranspiration and Groundwater Recharge: In fluences of Land Use, Soils, and
Lake-Effect Climate, Advances in Water Resources.
Kendall, C., E. M. Elliott, and S. Wankel, Tracing anthropogenic inputs of nitrogen to ecosystems, in Stable Isotopes in Ecology and Environmental Science, edited by R. Michener and
K. Lajtha, pp. 375–449, 2007.
Law, N., L. Band, and M. Grove, Nitrogen input from residential lawn care practices in
suburban watersheds in Baltimore county, MD, Journal of Environmental Planning and
Management, 47 (5), 737–755, doi:10.1080/0964056042000274452, 2004.
Leopold, L. B., and T. Maddock, The Hydraulic Geometrv of Stream Channels and Some
Physiographic Implications, Geological Survey Professional Paper, 252, 1953.
Mainston, C. P., and W. Parr, Phosphorus in rivers–ecology and management., The Science
of the total environment, 282-283, 25–47, 2002.
Maupin, M. A., and T. Ivahnenko, Nutrient Loadings to Streams of the Continental United
States From Municipal and Industrial Effluent, JAWRA Journal of the American Water
Resources Association, 47 (5), 950–964, doi:10.1111/j.1752-1688.2011.00576.x, 2011.
Mayer, B., E. W. Boyer, N. A. Jaworski, N. V. A. N. Breemen, R. W. Howarth, S. Seitzinger,
G. Billen, K. Lajtha, K. Nadelhoffer, D. Van Dam, L. J. Hetling, M. Nosal, and K. Paustian, Sources of nitrate in rivers draining sixteen watersheds in the northeastern U . S .:
Isotopic constraints, Biogeochemistry, 57/58, 171–197, 2002.

128

McMahon, G., L. Tervelt, and W. Donehoo, Methods for Estimating Annual Wastewater
Nutrient Loads in the Southeastern United States, Tech. rep., U.S. Geological Survey
Open File Report 2007-1040, Reston, Virginia, 2007.
MDEQ, Wellogic Wells Summary, 2005.
Moon, J. B., and H. J. Carrick, Seasonal variation of phytoplankton nutrient limitation in
Lake Erie, Aquatic Microbial Ecology, 48 (Carrick 2004), 61–71, 2007.
Moore, R. B., C. M. Johnston, K. W. Robinson, and J. R. Deacon, Estimation of Total
Nitrogen and Phosphorus in New England Streams Using Spatially Referenced Regression
Models, 2004.
NADP, National Atmospheric Deoposition Program (NRSP-3), Tech. rep., NADP Program
Office, Illinois State Water Survey, 2204 Griffith Dr., Champaign, Il 61820, 2007.
Natural Resource Conservation Service, S. S. S., Soil Survey Geographic (SSURGO)
Database, Tech. rep., United States Department of Agriculture.
Neff, J. C., E. A. Holland, J. Frank, W. H. Mcdowell, and K. M. Russell, The origin ,
composition and rates of organic nitrogen deposition : A missing piece of the nitrogen
cycle ?, Methods, pp. 99–136, 2002.
Nelder, J., and R. Mead, A Simplex Method for Function Minimization, The Computer
Journal, 7 (4), 308–313, 1965.
Nicholls, K. H., G. J. Hopkins, S. J. Standke, and L. Nakamoto, Trends in Total Phosphorus
in Canadian Near-Shore Waters of the Laurentian Great Lakes: 1976-1999, Journal of
Great Lakes Research, 27 (4), 402–422, doi:10.1016/S0380-1330(01)70656-9, 2001.
Nikolaidis, N., H. Heng, R. Semagin, and J. Clausen, Non-linear response of a mixed land use
watershed to nitrogen loading, Agriculture, Ecosystems & Environment, 67 (2-3), 251–265,
doi:10.1016/S0167-8809(97)00123-0, 1998.
NRCS, 1992 National REsource Inventory, Tech. rep., Iowa State University, Statistical
Laboratory, Ames, 1995.
NRCS, Agricultrual waste characteristics: Agricultural Waste Management Field Handbook,
Waste Management, (March), 2008.
Piatek, K. B., M. J. Mitchell, S. R. Silva, and C. Kendall, Sources of Nitrate in Snowmelt
Discharge: Evidence from Water Chemistry and Stable Isotopes of Nitrate, Water, Air,
and Soil Pollution, 165, 13–35, 2005.
PLUARG, Nonpoint Source Pollution Abatement in the Greate Lakes Basin: An Overview
of Post-PLUARG Developments, Tech. Rep. August, Nonpoint Source Control Task Force
of the Water Quality Board of the International Joint Commission, Windsor, Ontario,
1983.

129

Preston, S. D., and J. W. Brakebill, Application of Spatially Referenced Regression Modeling
for the Evaluation of Total Nitrogen Loading in the Chesapeake Bay Watershed, USGS
Scientific Investigations Report, 99-4054, 1999.
PRISM Climate Group, PRISM Climate Group, 2011.
Reneau, R., and D. Pettry, Phosphorus Distribution from Septic Tank Effluent in Coastal
Plain Soils, Journal of Environment Quality, 5 (1), 34–39, 1976.
Robertson, D. M., and D. a. Saad, Nutrient Inputs to the Laurentian Great Lakes by
Source and Watershed Estimated Using SPARROW Watershed Models, JAWRA Journal of the American Water Resources Association, 47 (5), 1011–1033, doi:10.1111/j.17521688.2011.00574.x, 2011.
Robertson, W., and J. Cherry, Hydrogeology of an unconfined sand aquifer and its effect on
the behavior of nitrogen from a large-flux septic system, Hydrogeology Journal, 0, 32–44,
1992.
Roth, K. S., and T. Dewald, The National Hydrography Dataset, Tech. rep., U.S. Geological
Survey, U.S. Environmental Protection Agency, 1999.
Ruddy, B. C., D. L. Lorenz, and D. K. Mueller, County-Level Estimates of Nutrient Inputs to
the Land Surface of the Conterminous United States , 1982 - 2001, National Water-Quality
Assessment Program: Scientific Investigations Report, 2006-5012, 2006.
Schulze, K., and M. Hunger, Simulating river flow velocity on global scale, Advances In
Geosciences, 5, 133–136, 2005.
Schwarz, G. E., A. B. Hoos, R. B. Alexander, and R. A. Smith, The SPARROW Surface
Water-Quality Model : Theory , Application and User Documentation, book 6 ed., U.S.
Geological Survey, Reston, Virginia, 2006.
Sigman, D. M., K. L. Casciotti, M. Andreani, C. Barford, M. Galanter, and J. K. B¨hlke, A
o
bacterial method for the nitrogen isotopic analysis of nitrate in seawater and freshwater.,
Analytical chemistry, 73 (17), 4145–53, 2001.
Smith, R. A., R. B. Alexander, and M. G. Wolman, Water-Quality Trends in the Nation’s
Rivers, Science, 235 (4796), 1607–1615, 1987.
Smith, R. A., G. E. Schwarz, and R. B. Alexander, Regional interpretation of water-quality
monitoring data, Water Resources Research, 33 (12), 2781–2798, 1997.
Spoelstra, J., S. L. Schiff, R. J. Elgood, R. G. Semkin, and D. S. Jeffries, Tracing the Sources
of Exported Nitrate in the Turkey Lakes Watershed Using 15 N/ 14 N and 18 O/ 16 O
isotopic ratios, Ecosystems, 4 (6), 536–544, doi:10.1007/s10021-001-0027-y, 2001.
Throssell, C. S., G. T. Lyman, M. E. Johnson, G. A. Stacey, and C. D. Brown, Golf Course
Environmental Profile Measures Nutrient Use and Management and Fertilizer Restrictions
, Storage , and Equipment Calibration, Applied Turfgrass Science, doi:10.1094/ATS-20091203-01-RS.Abstract, 2009.
130

USDA, USDA National Agricultural Statistics Service Cropland Data Layer, 2011.
USEPA, Onsite Wastewater Treatment Systems Manual, Tech. Rep. February, U.S. Environmental Protection Agency, 2002a.
USEPA, National Water Quality Inventory: 1996 Report to Congress, Tech. rep., United
States Environmental Protection Agency Office of Water, Washington, D.C., 2002b.
USEPA, National Water Quality Inventory : Report to Congress 2004 Reporting Cycle,
Tech. Rep. January, United States Environmental Protection Agency Office of Water,
Washington, D.C., 2009.
USGS, Streamflow Measurements for the Nation, 2012.
Vadeboncoeur, M. a., S. P. Hamburg, and D. Pryor, Modeled Nitrogen Loading to Narragansett Bay: 1850 to 2015, Estuaries and Coasts, 33 (5), 1113–1127, doi:10.1007/s12237010-9320-3, 2010.
Valiela, I., Assessment of models for estimation of land-derived nitrogen loads to shallow estuaries, Applied Geochemistry, 17 (7), 935–953, doi:10.1016/S0883-2927(02)00073-2, 2002.
Valiela, I., G. Collins, J. Kremer, K. Lajtha, M. Geist, B. Seely, J. Brawley, and
C. H. Sham, Nitrogen Loading From Coastal Watersheds To Receiving Estuaries:
New Method and Application, Ecological Applications, 7 (2), 358–380, doi:10.1890/10510761(1997)007[0358:NLFCWT]2.0.CO;2, 1997.
Valigura, R. A., R. B. Alexander, M. S. Castro, T. P. Meyers, H. W. Paerl, P. E. Stacey,
and E. R. Turner, Atmospheric Nitrogen Flux from the Watersheds of Major Estuaries
of the United States: An Application fo the SPARROW Watershed Model, Coastal and
Estuarine Studies, pp. 119–170, 2001.
Verhougstraete, M. P., Measuring microbial water quality responses to land and climate
using fecal indicator bacteria and molecular source tracking in rivers and near - shore
waters of Michigan, Ph.D. thesis, Michigan State University, 2012.
Vitousek, P. M., J. D. Aber, R. W. Howarth, G. E. Likens, P. A. Matson, D. W. Schindler,
W. H. Schlesinger, and D. G. Tilman, Human Alteration of the Global Nitrogen Cycle:
Sources and Consequences, Ecological Applications, 7 (3), 737–750, 1997.
Warncke, D., J. Dahl, L. Jacobs, and C. Laboski, Nutrient Recommendations for Field Crops
in Michigan, Extension Bulletin, E2904, 2004.
Withers, P. J. a., and H. P. Jarvie, Delivery and cycling of phosphorus in
rivers: a review., The Science of the Total Environment, 400 (1-3), 379–95, doi:
10.1016/j.scitotenv.2008.08.002, 2008.
Wolock, D. M., Base-flow Index Grid for the Conterminous United States. Open-File Report
03-236, Tech. rep., U.S. Geological Survey, Reston, Virginia, 2003.

131

Xue, D., J. Botte, B. De Baets, F. Accoe, A. Nestler, P. Taylor, O. Van Cleemput,
M. Berglund, and P. Boeckx, Present limitations and future prospects of stable isotope
methods for nitrate source identification in surface- and groundwater., Water research,
43 (5), 1159–70, doi:10.1016/j.watres.2008.12.048, 2009.
Zhang, T., A Spatially Explicit Model for Estimating Annual Average Loads of Nonpoint
Source Nutrient at the Watershed Scale, Environmental Modeling & Assessment, 15 (6),
569–581, doi:10.1007/s10666-010-9225-3, 2010.
Zhang, T., Distance-decay patterns of nutrient loading at watershed scale: Regression modeling with a special spatial aggregation strategy, Journal of Hydrology, 402 (3-4), 239–249,
doi:10.1016/j.jhydrol.2011.03.017, 2011.
Zhou, W., A. Troy, and M. Grove, Modeling residential lawn fertilization practices: integrating high resolution remote sensing with socioeconomic data., Environmental management,
41 (5), 742–52, doi:10.1007/s00267-007-9032-z, 2008.

132