Draft Ambient Water Quality Criteria Recommendations for Lakes
and Reservoirs of the Conterminous United States: Information
Supporting the Development of Numeric Nutrient Criteria
Prepared by:
U.S. Environmental Protection Agency
Office of Water
Office of Science and Technology (4304T)
Washington, DC
EPA Document Number: 820P20001

-------
Contents
Executive Summary	x
1	Introduction and Background	1
2	Problem Formulation	3
2.1	Management Goals	3
2.2	Assessment Endpoints and Risk Metrics	4
2.2.1	Aquatic Life Use	5
2.2.2	Recreational Use	8
2.2.3	Drinking Water Source	9
2.3	Risk Hypotheses	9
2.4	Analysis Plan	10
3	Analysis	11
3.1	Data	11
3.1.1	Biological Data	13
3.1.2	Chemical Data	14
3.1.3	Dissolved Oxygen and Temperature Profiles	14
3.1.4	Mapped Data	15
3.2	Stressor-Response Models	15
3.2.1	Zooplankton Biomass	15
3.2.2	Deepwater Hypoxia	23
3.2.3	Microcystin Concentration	40
3.2.4	Phosphorus-Chlorophyll a	46
3.2.5	Nitrogen-Chlorophyll a	55
4	Characterization	61
4.1	Other Measures of Effect and Exposure	61
4.2	Incorporating State Data	62
4.3	Existing Nutrient-Chlorophyll a Models	63
4.4	Limitations and Assumptions	64
ii

-------
4.5	Deriving State-Specified Criteria	66
4.6	Duration and Frequency	66
5	References	69
6	Appendix A: State Case Study: Chlorophyll o-Microcystin	89
6.1	Data	89
6.2	Statistical Analysis	89
6.3	Results	91
6.4	Criteria Derivation	95
7	Appendix B: State Case Study: Chlorophyll o-Hypoxia	96
7.1	Data	96
7.2	Statistical Analysis	96
7.3	Results	97
7.4	Criteria Derivation	100
8	Appendix C: State Case Study: Total Nitrogen-Chlorophyll a	102
8.1	Data	102
8.2	Statistical Analysis	102
8.3	Results	 103
8.4	Criteria Derivation	104
9	Appendix D: Operational Numeric Nutrient Criteria	106
iii

-------
Figures
Figure 1. Conceptual model linking increased nutrients to aquatic life use (Source: US EPA
2010a)	5
Figure 2. Conceptual model linking increased nutrient concentrations to public health
endpoints	9
Figure 3. Simplified conceptual model showing pathways selected for analysis	10
Figure 4. Schematic of network of relationships for modeling zooplankton biomass. Gray-
filled ovals: available observations; other nodes: modeled parameters; numbers in
parentheses refer to equation numbers in the text	16
Figure 5. Relationships between measured biovolume, Chi a, and estimated mean
phytoplankton biovolume. Solid lines: 1:1 relationship	20
Figure 6. Estimated relationships between zooplankton and Chi a for lakes > 7.2 m deep.
Left panel: Chi a vs. zooplankton abundance; middle panel: Chi a vs. zooplankton
biomass; right panel: Chi a vs. slope of the relationship between zooplankton
biomass and Chi a. Solid lines: mean relationships; shaded areas (left and middle
panels): 80% credible intervals about mean relationship; dashed lines (right panel):
50% credible intervals about mean relationship; open circles (left and right panels):
average of five samples nearest the indicated Chi a concentration; dotted
horizontal line (right panel): one example value of threshold for deriving a Chi a
criterion	21
Figure 7. Lakes designated as seasonally stratified dimictic lakes from the NLA data set	24
Figure 8. Illustrative examples of depth profiles of temperature, temperature gradient, and
DO. Dashed horizontal line: estimated depth of the bottom of the epilimnion	25
Figure 9. Schematic of hypoxia model. Numbers in parentheses refer to equation numbers
in the text	26
Figure 10. Chi a vs. DOm. DOm values. Gray-filled circles: values < 2 mg/L; solid line:
nonparametric fit to the data shown to highlight asymptotic relationship	29
Figure 11. Demers and Kalff (1993) predicted stratification day vs. model mean estimate.
Solid line: the 1:1 relationship	30
iv

-------
Figure 12. Relationships between individual predictors and D0m, holding other variables
fixed at their mean values. Solid line: mean relationship; gray shading: 90% credible
intervals	31
Figure 13. Model predicted DOm vs. observed DOm. Open circles: individual samples; solid
line: 1:1 relationship	31
Figure 14. Contours of modeled mean lake temperature computed at the overall mean
elevation and mean sampling day	33
Figure 15. Relationship between lake temperature and sampling day (left panel) and
elevation (right panel). Variables that are not plotted are fixed at their mean
values. Gray shading: 90% confidence intervals; solid lines shows the mean
relationships	34
Figure 16. Days of the year that mixed layer temperatures decrease below the critical
temperature for cool-water species. Small dots: lakes in which mixed layer
temperatures never exceed 24 C	35
Figure 17. Days of the year that mixed layer temperatures decrease below the critical
temperature for cold-water species. Small dots: lakes in which mixed layer
temperatures do not decrease below 18 C during the summer; contours: effects of
large differences in elevation across lakes in the western U.S	35
Figure 18. Simplified DO profile used to compute threshold for DOm. Open circle: the
targeted condition of DO at 5 mg/L, 30 cm below the thermocline	37
Figure 19. Effects of other predictors on Chi a criteria. Solid lines: relationship between Chi
a and DOm at median values for all other variables; dashed line: DOm = 0.3 mg/L;
dotted lines: 25th and 75th percentiles of days elapsed since stratification (left
panel), 25th and 75th percentiles of mean DOC concentrations (middle panel), and
depth below thermocline of 4 m and 20 m (right panel)	38
Figure 20. Schematic showing relationship between different variables predicting MC.
Numbers in parentheses: refer to equation numbers in the text	40
v

-------
Figure 21. Modeled relationships for the microcystin model. Left panel: relationship
between Chi a and phytoplankton biovolume; open circles: observed
measurements of Chi a and phytoplankton biovolume; solid line: has a slope of 1.
Middle panel: relationship between Chi a and cyanobacterial relative biovolume;
open circles: average cyanobacterial relative biovolume in ~20 samples at the
indicated Chi a concentration; solid line: estimated mean relationship; gray
shading: 90% credible intervals about the mean relationship; vertical axis: has been
logit-transformed. Right panel: relationship between cyanobacterial biovolume and
MC; open circles: average MC in ~20 samples at the indicated cyanobacterial
biovolume; solid line: estimated mean relationship; gray shading: 90% credible
intervals about the mean relationship; small filled circles: Chi a bins in which MC in
all samples was zero	43
Figure 22. Example of derivation of Chi a criterion to protect recreational uses based on
targeted MC of 8 ng/L and exceedance probability of 1%. Top panel-open circles:
observed values of microcystin and Chi a for samples in which MC was greater than
the detection limit; solid line: predicted MC that will be exceeded 1% of the time
for the indicated Chi a concentration; gray shading: 50% credible intervals about
mean relationship; solid vertical and horizontal line segments: candidate Chi a
criterion based on targeted MC. Bottom panel: proportion of samples for which
microcystin was not detected in ~100 samples centered at the indicated Chi a
concentration	45
Figure 23. Schematic representation of compartment model for TP. PdiSS: dissolved P; Chi:
Chlorophyll a; Turb: total turbidity; Turbnp: turbidity attributed to
nonphytoplankton sources. Shaded box for Turbp: a variable inferred by the
model; numbers in parentheses: refer to equation numbers in the text. Equations
(28)-(30) and equations (33)-(35) describe the distributions of turbidity and TP
measurements and are not shown in the schematic	47
Figure 24. Turbidity vs. Chi a. Solid line: the limiting relationship between Chi a and
turbidity when contribution of allochthonous sediment is negligible	50
Figure 25. Relationship between Turbnp, Pdiss, and lake depth. Open circles: mean estimate
of parameter value in each of 30 lake depth classes	51
vi

-------
Figure 26. Ecoregion-specific values of loge[di), P bound to nonphytoplankton suspended
sediment	52
Figure 27. TP versus Turbnp and Chi a. Solid lines: the limiting relationship between the
indicated variable and TP; gray shaded areas: the 90% credible intervals about the
mean relationship	53
Figure 28. Example of deriving TP criteria for a Chi a target of 10 ng/L for data from one
ecoregion (Southeastern Plains). Open circles: all data; filled circles: data from the
ecoregion; solid line: limiting TP-Chl a relationship from compartment model;
dashed line: ambient TP-Chl a relationship taking into account contributions from
nonphytoplankton sediment for a 3-m deep lake; solid horizontal and vertical line
segments: Chi a target and associated TP criteria	54
Figure 29. Variation in the concentration of N bound in phytoplankton among Level III
ecoregions at the overall mean Chi a = 9.3 ng/L. Gray scale shows N concentrations
in ng/L	57
Figure 30. Variation in DON Level III ecoregions at an overall mean DOC = 5.6 mg/L. Gray
scale shows N concentrations in ng/L	58
Figure 31. TN-DIN vs. Chi a and DOC. Solid lines: the limiting relationship between each
variable and TN-DIN; shaded area: the 95% credible intervals about this mean
relationship	59
Figure 32. Illustrative example of deriving TN criteria for a Chi a target of 10 ng/L for one
ecoregion (Southeastern Plains). Open circles: all data; filled circles: data from
selected ecoregion; solid line: limiting TN-DIN vs. Chi a relationship from
compartment model; dashed line: mean ambient TN-DIN vs. Chi a relationship
taking into account mean DOC observed within the selected ecoregion: shaded
area: 80% credible intervals about mean relationships; horizontal and vertical solid
line segments: Illustrative Chi a target and associated TN criteria	60
Figure 33. Variation in the relationship between Chi a and cyanobacterial-relative
biovolume among states. PropCyano: predicted mean relative biovolume of
cyanobacteria at an illustrative Chi a = 20 ng/L	92
vii

-------
Figure 34. Comparison of Chi o/cyanobacterial-relative biovolume relationships in Iowa.
Filled gray: 90% credible intervals for estimate of relationship using only NLA data
collected in Iowa; solid and dashed lines: mean and 90% credible intervals for
estimate of relationship using both Iowa state and NLA data	93
Figure 35. Comparison of predicted relationship between Chi a and MC for the state-
national model (left panel) and a model using only Iowa state data (right panel).
Open circles: average MC concentration computed in ~15 samples at the indicated
Chi a; solid lines: mean relationship; dashed lines: 5th to 95th percentiles of
distribution of means of 15 samples drawn from predicted distribution	94
Figure 36. MC and Chi a measurements in Iowa. Top panel-open circles: observed values of
microcystin and Chi a for samples in which MC was greater than the detection
limit; solid line: predicted MC that will be exceeded 1% of the time for the indicated
Chi a concentration; gray shading: 50% credible intervals about mean relationship;
horizontal and vertical line segments: candidate Chi a criteria based on targeted
MC. Bottom panel: proportion of samples for which microcystin was not detected
in ~100 samples centered at the indicated Chi a concentration	95
Figure 37. Observed DOm vs. Chi a, sampling day, DOC, and depth below the thermocline.
Open circles: NLA data; filled circles: Missouri	98
Figure 38. Estimated first day of stratification for Missouri lakes (left panel) and NLA lakes
(right panel)	99
Figure 39. Model coefficients estimated for models for Missouri data, NLA data, and
combined data. Thick line segment: 50% credible intervals; thin line segment: 90%
credible intervals; vertical dashed line: coefficient value of zero	99
Figure 40. Relationships between day of year and DOm for six Missouri lakes. Different line
and symbol colors in each panel correspond to data collected within different years
with at least three samples. Open gray circles: other samples collected at each lake.. 100

-------
Figure 41. Relationship between Chi a and DOm in an illustrative lake with depth below
thermocline at 13m, DOC at 3.5 mg/L, and 120 days after spring stratification. Solid
line: mean relationship; gray shading: 50% credible intervals about mean
relationship from combined Missouri-NLA model; dashed line: 50% credible
intervals about mean relationship from Missouri-only model; dotted line: DOm = 0
mg/L	101
Figure 42. Chi a vs. TN-DIN in Iowa. Open circles: data collected by Iowa DNR; filled circles:
data collected by NLA in Iowa; solid lines: 95% credible intervals for limiting
relationships between Chi a and TN-DIN estimated using both NLA and IDNR data;
shaded gray area: 95% credible intervals for limiting relationships estimated using
only IDNR data	103
Figure 43. Chi a vs. TN-DIN in Beeds Lake, Iowa. Open circles: observed data; gray shading:
90% credible intervals for predicted relationship based on only IDNR data; solid
lines: 90% credible intervals for predicted relationship using both IDNR and NLA
data	104
Figure 44. Lake-specific criteria derivation using combined Iowa-NLA model for two
different lakes in Iowa. Open circles: observed values of TN-DIN and Chi a in Iowa
for each lake; gray shading: 50% credible intervals about the mean relationship;
solid line: mean relationship calculated using mean DOC concentration in lake;
horizontal and vertical line segments: TN criterion calculation for illustrative Chi a
target of 15 ng/L	105
Figure 45. Example distribution of 10 TP measurements. Note that the horizontal axis is log-
scaled	107
Figure 46. Example of defining an operational criterion magnitude. Solid line: the
cumulative probability of observing a single sample TP lower than or equal to the
indicated value if the true annual mean was exactly equal to the criterion (TP = 60
ju,g/L); dashed line: the cumulative probability for the average of four samples;
black arrows: operational criteria for one sample; gray arrows: operational criteria
associated with four samples	109
ix

-------
Executive Summary
While certain levels of nutrients are essential for healthy aquatic ecosystems, excess
nutrients can degrade the condition of water bodies worldwide, and in lakes and reservoirs
(hereafter, referred to only as "lakes" unless noted otherwise), the effects of excess nitrogen (N)
and phosphorus (P) may be particularly evident. High levels of nutrient loading commonly
stimulate excess growth of algae, which can limit the recreational use of lakes. Overabundant
algae also increase the amount of organic matter in a lake, which, when decomposed, can
depress dissolved oxygen (DO) concentrations below the levels needed to sustain aquatic life. In
extreme cases, the depletion of DO causes fish kills. Nutrient pollution can stimulate the excess
growth of nuisance algae, such as cyanobacteria, which can produce cyanotoxins that are toxic
to animals and humans. Elevated concentrations of cyanotoxins can reduce the suitability of a
lake for recreation and as a source of drinking water.
Numeric nutrient criteria provide an important tool for managing the effects of nutrient
pollution by providing nutrient goals that ensure the protection and maintenance of designated
uses. The United States (U.S.) Environmental Protection Agency (EPA) published recommended
numeric nutrient criteria for lakes and reservoirs in 2000 and 2001 for 12 out of 14 ecoregions of
the conterminous U.S. Those criteria were derived by analyzing available data on the
concentrations of total nitrogen (TN), total phosphorus (TP), chlorophyll a (Chi a), and Secchi
depth.
Scientific understanding of the relationships between nutrient concentrations and
deleterious effects in lakes has increased since 2001, and standardized, high-quality data
collected from lakes across the U.S. have become available. In this document, the EPA describes
analyses of these new data and provides draft models from which numeric nutrient criteria can
be derived. The draft criteria models would, if finalized, replace the recommended numeric
nutrient criteria of 2000 and 2001. The draft criteria models are provided in accordance with the
provisions of Section 304(a) of the Clean Water Act (CWA) (Title 33 of the United States Code
[U.S.C.]  1314(a)) for the EPA to revise ambient water quality criteria from time to time to
reflect the latest scientific knowledge. CWA Section 304(a) water quality criteria serve as
recommendations to states and authorized tribes for defining ambient water concentrations
that will protect against adverse effects to aquatic life and human health. The ecological and
health protective responses on which the draft criteria models are based were selected by
x

-------
applying a risk assessment approach to explicitly link nutrient concentrations to the protection
of designated uses.
The draft criteria models are nonregulatory. When they are finalized, states may use the
recommended models to derive candidate nutrient criteria for each applicable designated use
and, after demonstrating that the criteria protect the most sensitive designated use, adopt the
criteria into their state standards. States may also modify the criteria to reflect site-specific
conditions or establish criteria based on other scientifically defensible methods (Title 40 of the
Code of Federal Regulations [CFR] 131.11(b)). When finalized, the updated recommended CWA
Section 304(a) nutrient criteria for lakes will not compel a state to revise current EPA-approved
and adopted criteria, total daily maximum load nutrient load targets, or N or P numeric values
established by other scientifically defensible methods. As part of their triennial review, if a state
uses its discretion to not adopt new or revised nutrient criteria based on these CWA Section
304(a) criteria models, then the state shall provide an explanation when it submits the results of
its triennial review (40 CFR 131.20(a)).
Following the risk assessment paradigm, the EPA first defined water quality
management goals for numeric nutrient criteria, and then defined assessment endpoints and
metrics that are associated with achieving these goals and are sensitive to increased nutrient
concentrations. The water quality management goals are articulated as designated uses in
Section 101(a)(2) of the CWA (33 U.S.C.  1251) (i.e., the protection and propagation offish,
shellfish, and wildlife [aquatic life] and recreation in and on the water). Another common
designated use for lakes is to serve as drinking water sources. Excess loads of nutrients can lead
to excessive growth of phytoplankton that can adversely impact designated uses in different
ways, described below as assessment endpoints and metrics. The EPA modeled stressor-
response relationships using these endpoints and metrics to derive draft recommended numeric
nutrient criterion models (Table 1).
For aquatic life, the EPA identified two assessment endpoints. The first endpoint is
zooplankton biomass, and the risk metric is the relationship between zooplankton and
phytoplankton biomass, which quantifies the degree to which energy produced by
phytoplankton at the base of the food web is transferred to zooplankton and subsequently to
higher trophic levels. When excess nutrients are available, phytoplankton biomass can increase
at rates that exceed the capacity of zooplankton to consume. The draft risk metric is one in
xi

-------
which the rate of change of zooplankton biomass relative to phytoplankton biomass is
approximately zero. This condition describes a lake in which the biomass of grazing biota (i.e.,
zooplankton) does not increase with increases in food (i.e., phytoplankton), and primary
production at the base of the food web is weakly linked to production at higher trophic levels.
This endpoint applies to all lakes in the conterminous U.S.
The second aquatic life endpoint is cool- and cold-water fish, and the risk metric is the
DO concentration in deep water that protects against mortality of these fish. Excess nutrients
typically increase primary productivity, which then increases the amount of organic matter in a
lake. Then, in the deep waters of a lake, DO is consumed as this organic matter is decomposed,
leading to hypoxic and anoxic conditions. The draft risk metric is the daily DO concentration,
calculated as a depth-averaged value below the thermocline, which can be reduced to
concentrations insufficient to support some fish species during the critical period of the summer
when they require deep, cold waters to escape high temperatures at shallower depths. This
endpoint applies to seasonally stratified, dimictic lakes harboring cool- and cold-water fish.
For recreational uses and drinking water sources, the assessment endpoint is human
health. For recreational uses, the EPA selected the risk metric as the concentration of
microcystin associated with adverse effects on children (specifically, liver toxicity) from
incidental ingestion of water during recreation. When excess nutrients are available,
phytoplankton communities can shift toward a greater abundance of cyanobacteria that can
release cyanotoxins, and microcystins are the most commonly monitored and measured
freshwater cyanotoxin in the U.S. The threshold for the draft risk metric is 8 micrograms per liter
(Hg/L), based on recently published national recommendations for human health recreational
water quality criteria and swimming advisories for cyanotoxins (US EPA 2019). For the drinking
water use, the EPA selected as the risk metric the concentration of microcystins associated with
adverse effects on children resulting from oral exposure to drinking water (0.3 ng/L), consistent
with the health advisory for microcystins (US EPA 2015b). This microcystin concentration from
the health advisory applies to finished drinking water; however, the EPA is aware that states or
authorized tribes apply water quality standards for protecting drinking water sources to either
the ambient source water before treatment or to the finished drinking water after treatment.
The ability of treatment technologies to remove microcystin is too variable for the EPA to set a
national recommendation for a protective ambient source water concentration that would yield
xii

-------
a protective concentration after treatment. If a state or authorized tribe applies the health
advisory standard to finished drinking water, then they can account for the expected treatment
in their facilities and select a higher microcystin concentration in the ambient source water that
would result in the targeted microcystin concentration in the finished drinking water.
Table 1. Summary of designated uses and associated measures of effect and exposure
Designated
use
Assessment endpoint
Risk metric
Applicability


Rate of change of zooplankton

Aquatic life
Zooplankton biomass
biomass relative to
phytoplankton biomass
All lakes
Aquatic life
Cool- and cold-water fish
Daily depth-averaged DO below
the thermocline
Dimictic lakes with
cool- or cold-water
fish
Recreation
Human health
Microcystin concentration to
prevent liver toxicity in children
All lakes
Drinking water
Human health
Microcystin concentration to
prevent liver toxicity in children
All lakes
Data used in this analysis were collected in the EPA's National Lakes Assessment (NLA),
which sampled lakes across the conterminous U.S. in 2007 and 2012. Most of the sampled lakes
were selected randomly so the resulting data represent the characteristics of the full population
of lakes in the conterminous U.S. At each lake, standardized protocols were used to collect
extensive measurements of biotic and abiotic characteristics.
This document describes statistical stressor-response models that relate Chi a
concentrations to each of the risk metrics and that relate TN and TP concentrations to Chi a. A
hierarchical Bayesian network is specified for each model to represent the effects of different
variables on the relationship of interest. For example, microcystin is related to cyanobacteria
biovolume, which is then linked to Chi a concentration. The Bayesian network models can
directly represent the processes that govern the relationships of interest and facilitate the use
of other data sets in conjunction with data from the EPA's NLA. When coupled with the targets
for each response, the draft models provide candidate Chi a, TN, and TP criteria
recommendations that states may then use with state risk management decisions to
demonstrate they are protective of different designated uses. For lakes with multiple use
designations, the states shall adopt criteria that protect the most sensitive use.
xiii

-------
Models provided in this document are based on national data, but states often collect
extensive data during routine monitoring. Incorporating local data into the national models can
refine and improve the precision of the stressor-response relationships. In the appendices of
this document, the EPA describes three case studies in which state monitoring data have been
combined with national data, yielding models that can be used to derive recommended numeric
nutrient criteria that account for both unique local conditions and national, large-scale trends.
xiv

-------
1 Introduction and Background
While certain levels of nutrients are essential for healthy aquatic ecosystems, nutrient
pollution, or the excess loading of nitrogen (N) and phosphorus (P), can degrade the conditions
of water bodies worldwide, and in lakes the effects of nutrient pollution are often most
evident. One visible consequence of nutrient pollution in lakes and reservoirs (hereafter,
referred to only as "lakes" unless noted otherwise) is cultural eutrophication, an increase in
primary productivity and algal abundance that increases the amount of organic matter in a
water body (Smith et al. 2006, Smith and Schindler 2009). Decomposition of organic matter
reduces dissolved oxygen (DO) concentrations in the water column, especially in deeper waters
under stratified conditions. These hypoxic conditions are inhospitable to most aquatic species
and reduce their ability to survive within a particular lake (Jones et al. 2011, Scavia et al. 2014).
Nutrient pollution also favors the growth of undesirable, nuisance algae (e.g.,
cyanobacteria), some of which produce cyanotoxins (Paerl and Otten 2013). Many species of
cyanobacteria are superior competitors for light compared to other phytoplankton. Hence, in
lakes with nutrient pollution, cyanobacteria can dominate by reducing the light available to
other phytoplankton (Carey et al. 2012). A number of other mechanisms, including superior
uptake rates for carbon dioxide and an ability to migrate vertically in the water column, also
may explain the frequent occurrence of cyanobacteria dominance in eutrophic systems (Dokulil
and Teubner 2000). Cyanobacteria dominance can interfere with the designated uses of a lake
because cyanobacteria not only can form unsightly and odorous surface scums (reducing the
aesthetic appeal of the lake for recreation) (Paerl and Ustach 1982), but also can produce
cyanotoxins that can limit the use of the lake as both a source of drinking water and for
recreation (Cheung et al. 2013). Many species of cyanobacteria are also less palatable than
other algae to grazing organisms, and so, increases in cyanobacterial abundance can alter lake
food webs and reduce the efficiency with which energy from primary production is transferred
to higher trophic levels (Elser 1999, Filstrup et al. 2014a, Heathcote et al. 2016).
Nutrient pollution in lakes and resulting adverse environmental effects are widespread
in the United States (U.S.). Nutrient pollution occurs in lakes of different sizes, in catchments
with varying land uses, and in different climates. The U.S. Environmental Protection Agency
(EPA) has long recognized the effects of nutrient pollution and has recommended that states
1

-------
and authorized tribes (hereafter, "states"), acting under their Clean Water Act (CWA)
authorities, adopt numeric nutrient criteria as one way to facilitate the management of these
effects. A state's numeric nutrient criteria (1) provide nutrient goals to protect and maintain the
designated uses of a water body (Title 33 of the United States Code [U.S.C.]  1313(c)), (2)
provide thresholds that allow the state to make accurate water quality assessment decisions (33
U.S.C.  1313(d)), and (3) provide targets for restoration of water bodies that can guide waste
load allocation decisions (33 U.S.C.  1313(d)). To assist states and authorized tribes in deriving
numeric nutrient criteria, the EPA has published a series of technical support documents on
methods for deriving criteria for lakes and reservoirs (US EPA 2000a), streams and rivers (US EPA
2000b), wetlands (US EPA 2008), and estuaries and coastal waters (US EPA 2001). A technical
support document on using stressor-response relationships for deriving numeric nutrient
criteria has also been published (US EPA 2010a). In 2000 and 2001, under its authority described
in Section 304(a) of the CWA (33 U.S.C.  1314(a)), the EPA issued 12 documents that provided
recommended numeric nutrient criteria for lakes, streams, and rivers in different ecoregions of
the U.S. These criteria were derived by using available monitoring data to estimate the
concentrations of total nitrogen (TN) and total phosphorus (TP) that were expected to occur in
least-disturbed reference water bodies in different nutrient ecoregions.
In accordance with the provisions of Section 304(a) of the CWA, which directs the EPA to
revise ambient water quality criteria from time to time to reflect the latest scientific knowledge,
the EPA is issuing draft revisions to numeric nutrient criteria recommendations for lakes based
on analyses of newly available, national-scale data and reflecting advances in scientific
understanding of the relationship between excess nutrients and adverse effects in lakes. The
draft criteria recommendations are models that generate numeric nutrient criteria based on
national data and state risk management decisions. State data, if available, can be incorporated
into the national criteria models to compute relationships that more accurately represent local
conditions. In deriving these draft models, the EPA uses a risk assessment framework (Norton et
al. 1992, US EPA 1998, 2014) to identify assessment endpoints that relate directly to the water
quality management goals for U.S. lakes specified by the CWA and that are sensitive to
increased concentrations of N and P. Then, the EPA uses stressor-response analysis to estimate
relationships between increased N and P (estimated by measurements of TN and TP) and
different risk metrics directly linked to the assessment endpoints (US EPA 2010a). Draft national
criteria models are provided for both TN and TP as the simultaneous control of both nutrients
2

-------
provides the most effective means of controlling the deleterious effects of nutrient pollution (US
EPA 2015a, Paerl et al. 2016). After the public comment period and any consequent revisions to
the draft, the EPA intends to finalize the recommended stressor-response criteria models to
replace the ecoregion-specific nutrient criteria recommended previously for lakes that were
based on a reference distribution approach.
The remaining sections of this document are organized broadly according to the steps of
risk assessment: (1) problem formulation, (2) analysis, and (3) characterization. The purpose of
this document is to provide the technical details underlying the estimation of relationships
between increased nutrient concentrations and different responses, as well as details regarding
the derivation of draft numeric nutrient criteria recommendations using the national models.
Once the recommended criterion models are finalized, states may use them to derive candidate
nutrient criteria and, after demonstrating that the criteria protect designated uses, adopt the
criteria into their state water quality standards. States may also modify the criteria to reflect
site-specific conditions or establish criteria based on other scientifically defensible methods (40
CFR 131.11(b)). For waters with multiple use designations, the state shall adopt criteria that
support the most sensitive designated use (40 CFR 131.11(a)(1)). Water quality standards
adopted by states are subsequently subject to review by the EPA, pursuant to Section 303(c) of
the CWA (33 U.S.C.  1313(c)).
2 Problem Formulation
2.1 Management Goals
The EPA focused on protecting uses that reflect management goals articulated in
Section 101(a)(2) of the CWA (33 U.S.C.  1251), which include maintaining conditions so
different water bodies support aquatic life use (i.e., providing for the protection and
propagation of fish, shellfish, and wildlife), recreation (i.e., providing for recreation in and on the
water), and use of the water body as a source of drinking water. Under the CWA, it is a state's
responsibility to designate uses for its waters, and many states have designated uses that
provide for aquatic life and recreation uses. Some states have also designated waters as sources
of drinking water. The EPA focuses on aquatic life, recreation, and drinking water source
because they represent uses that are particularly sensitive to increased concentrations of N and
3

-------
P. States can derive candidate nutrient criteria for each of the applicable designated uses in
their lakes and, by comparing these criteria, identify the most sensitive use. Water quality
criteria adopted by states for waters with multiple use designations must support the most
sensitive use (40 CFR 131.11(a)).
2.2 Assessment Endpoints and Risk Metrics
The next step in problem formulation is to define assessment endpoints that can be
used to quantify attainment of the management goals. Each of the management goals
expressed in terms of different designated uses was associated with different assessment
endpoints. Protection of recreational uses and drinking water sources pertains to public health
rather than ecological health, and hence, the assessment endpoint is human health for these
two designated uses. For aquatic life, the procedures of ecological risk assessment were
followed to select assessment endpoints defined as "explicit expressions of the actual
environmental values that are to be protected" (US EPA 1998). Three considerations guided the
selection of these endpoints: ecological relevance, susceptibility to the stressor of interest (i.e.,
increased nutrient concentrations in the present case), and relevance to management goals.
After selecting the assessment endpoints, the EPA developed conceptual models that
represented current understanding of the linkages between increased N and P concentrations
and effects on the assessment endpoint and management goals (Figure 1). The conceptual
models were used to select specific risk metrics that quantified key steps along the causal path
linking increased N and P concentrations to deleterious effects on aquatic life and public health.
The final selections for the draft recommendations were also influenced by the availability of
data at the continental spatial scales considered in this analysis. These risk metrics were used as
the response variables in stressor-response analysis. For a narrative description of the
conceptual model, refer to Using Stressor-Response Relationships to Derive Numeric Nutrient
Criteria (US EPA 2010a).
4

-------
Nutrients
Hydrology-Geology;
Natural Vegetation
Lake Morphology;
Retention Time
Primary Production
Color
Temperature
Aquatic Macrophytes
Phytoplankton
Suspended
Sediment
Nuisance
Algae
HABs
Respiration
Organic Matter
Food Quality
Food Quantity
Dissolved
Oxygen
Water Clarity;
Light Availability
Biological Communities
Zooplankton
SAV
Benthic
Communities
Fish
Aquatic Life Support
Figure 1. Conceptual model linking increased nutrients to aquatic life use (Source: US EPA 2010a).
2.2.1 Aquatic Life Use
Nutrient pollution and eutrophication can affect the health of the lake biological
community via many pathways (Figure 1). As discussed earlier, increased nutrients typically
stimulate primary productivity and increase the amount of organic matter in a lake.
Decomposition of the organic matter depletes the DO in the water, reducing the suitability of
deeper waters as habitat for fish and invertebrates (Cornett 1989). Increased production and
respiration also can increase the range of acidity (pH) throughout the day-night cycle in some
lakes (Schindler et al. 1985), reducing the suitability of shallow waters as habitat for certain
species. Increased algal biomass also reduces water clarity, and the reduction in light availability
limits the depths at which submerged aquatic vegetation can persist (Phillips et al. 2016).
Reduced water clarity can also shift fish assemblage composition away from species that depend
on sight for foraging (De Robertis et al, 2003). Further, high nutrient concentrations favor the
growth of cyanobacteria, which are less palatable to grazing species than other phytoplankton,
altering the food web of the lake (Haney 1987).
5

-------
The EPA selected two assessment endpoints to characterize the health of aquatic life in
lakes: (1) zooplankton biomass, which is applicable to all lakes, and (2) cool- and cold-water fish
in dimictic lakes. For the second endpoint, the EPA selected depth-averaged DO concentration
as the risk metric. In dimictic stratified lakes with cool-water fish, criteria based on zooplankton
biomass and DO can be compared, and the more stringent criterion applied to ensure that
aquatic life is protected. Collectively, the two assessment endpoints provide a broad assessment
of the health of the lake biological community. Data were also available for each endpoint, and
each endpoint quantified well-studied effects of nutrient pollution.
2.2.1.1 Zooplankton biomass
The rate of change of zooplankton biomass compared to the rate of change of
phytoplankton biomass quantifies changes in the shape of biomass pyramids in lakes (Elton
1927). Biomass pyramids provide a graphical depiction of the amount of biomass at different
trophic levels, and typically, the biomass of primary producers (at the bottom of the pyramid)
exceeds the biomass of herbivores and carnivores at successively higher levels of the pyramid. In
lakes, the ratio of herbivore biomass (i.e., zooplankton) to primary producer biomass (i.e.,
phytoplankton) (Z:P) has been observed to decrease along eutrophication gradients (Leibold et
al. 1997). Reasons for the decreasing trend in Z:P have been the subject of some debate, much
of which centers on the relative importance of top-down versus bottom-up food web effects.
For zooplankton, top-down forces consist mainly of the effects of planktivore fish consuming
zooplankton biomass (Jeppesen et al. 2003) and bottom-up forces include changes in the
quantity and quality of the phytoplankton assemblage on which zooplankton feed (Filstrup et al.
2014a). With excess nutrients, one particularly relevant bottom-up mechanism is the decrease
in the edibility of the phytoplankton assemblage associated with the increased dominance of
cyanobacteria with increasing levels of eutrophication. Laboratory studies demonstrate that the
lack of highly unsaturated fatty acids in the cyanobacteria negatively affects the growth rates of
a common zooplankton species (Daphnia) (Demott and Miiller-Navarra 1997, Persson et al.
2007). Field observations (Miiller-Navarra et al. 2000) and microcosm experiments (Park et al.
2003) have added further support for this finding. Many cyanobacteria also present physical
challenges to grazers, collecting in colonies or filaments that are too large to be consumed
(Bednarska and Dawidowicz 2007), or surrounding themselves with gelatinous sheaths (Vanni
1987). Altered elemental stoichiometry and, hence, nutritional quality of phytoplankton under
different levels of eutrophication may also influence zooplankton biomass (Hessen 2008).
6

-------
While Z:P has traditionally been used to compare biomass pyramids among different
systems (Hessen et al. 2006), the rate of change of zooplankton biomass with respect to
increasing phytoplankton biomass (AZ/AP) provides a more informative measure of the effects
of eutrophication on food web function for the purposes of informing the derivation of numeric
nutrient criteria (Yuan and Pollard 2018). This rate of change can be thought of as the slope of
the relationship between Z and P. In most lake food webs, any increase in the basal resources
(i.e., phytoplankton biomass) would be expected to be associated with a corresponding increase
in the biomass of consumers of those resources (i.e., zooplankton biomass), and the slope
between Z and P would be positive. In eutrophic lakes, however, increases in phytoplankton
biomass often are not associated with an increase in zooplankton biomass, and the slope
(AZ/AP) approaches zero (Leibold et al. 1997, Hessen et al. 2006, Heathcote et al. 2016). Based
on this observation, the EPA used the rate of change in zooplankton biomass relative to changes
in phytoplankton biomass (AZ/AP) as a measure of the effect of excess nutrients on lake food
webs.
2.2.1.2 Dissolved oxygen
Excess nutrients typically increase primary productivity, which increases the amount of
organic matter in a lake. Then, DO is consumed as the organic matter is decomposed, leading to
hypoxic and anoxic conditions (Figure 1). Low concentrations of DO limit the extent to which
habitat is available to fish and zooplankton (Colby et al. 1972, Tessier and Welser 1991,
Vanderploeg et al. 2009), and oxygen availability is a key determinant of the quality and quantity
of habitat available to aquatic biota in many lakes (Evans et al. 1996). Although hypoxia occurs
naturally in a small number of systems (Diaz 2001), anthropogenic nutrient loads have greatly
increased the occurrence of hypoxia worldwide (Jenny et al. 2016). Deoxygenation of lake water
typically begins near the lake bottom and proceeds to shallower depths over the summer,
especially in stratified, relatively deep lakes, where the replenishment of DO from surface
mixing is restricted (Cornett 1989, Wetzel 2001). Therefore, an increasing proportion of the
deeper waters of a lake can become uninhabitable for certain organisms over the course of the
summer (Molot et al. 1992). Exclusion of deeper waters as viable habitat, in particular, can
disproportionately affect particular species of adult and juvenile fish (Lienesch et al. 2005).
Another strong determinant of the available habitat for fish and zooplankton is water
temperature. Summer brings a longer photoperiod and more intense solar insolation, which
7

-------
increases water temperatures near the surface of many lakes to levels harmful to certain species
(Ferguson 1958, Eaton and Scheller 1996). The viable habitat for cool- and cold-water species, in
particular, can be restricted by surface warming (Jacobson et al. 2010, Arend et al. 2011). In
contrast to deoxygenation, warming starts at the surface of the lake and proceeds to deeper
depths over the course of the summer. Therefore, certain species offish are "squeezed"
between increasing temperatures at shallow depths and decreasing DO at deeper depths
(Coutant 1985, Stefan et al. 1996, Lee and Bergersen 1996, Plumb and Blanchfield 2009),
requiring them to choose between suboptimal temperatures or oxygen (Arend et al. 2011).
Under those conditions, the metalimnion and the upper edge of the hypolimnion can provide an
important refuge, and even a thin layer of cool water with sufficient DO can provide an
important habitat for supporting fish health through the warmest summer days. Because they
often can tolerate lower DO concentrations than fish, zooplankton can retreat to deeper depths
of the hypolimnion to escape fish predation, but are also limited ultimately by low DO
concentrations (Tessier and Welser 1991, Stemberger 1995).
Based on these considerations, the mean concentration of DO below the thermocline
was identified by the EPA as an appropriate metric for assessing risks to cool- and cold-water
fish in seasonally stratified, dimictic lakes. In those lakes during the summer, the availability of
cool-water habitat is constrained by deepwater DO concentrations, and so, this risk metric links
increased nutrient concentrations to deleterious effects on fish and zooplankton in deep lakes.
2.2.2 Recreational Use
The EPA selected the concentration of cyanotoxins as the risk metric linking increased
nutrients to the suitability of lake water for primary and secondary contact recreation. Increased
nutrient concentrations and an attendant increase in cyanobacterial abundance can increase
concentrations of cyanotoxins (Figure 2), which cause adverse effects on the health of people
exposed to the water (US EPA 2019). One of the most commonly occurring types of cyanotoxins
in freshwaters is microcystins (based on available data). To protect recreational uses of lakes,
the EPA identified microcystin concentration (MC) as the best risk metric because of the
availability of NLA data (US EPA 2010b) and because MC thresholds for recreational exposures
have recently been published (US EPA 2019).
8

-------
2.2.3 Drinking Water Source
Increased nutrient concentrations and an attendant increase in cyanobacteria can
increase concentrations of cyanotoxins, which are toxic when consumed at certain
concentrations and quantities (Figure 2) (Chorus 2001, Stewart et al. 2008, US EPA 2015b). As
was done for recreational use, the EPA selected MC in lake source water as the relevant risk
metric for the drinking water use.

Hydrology-Geology;
Natural Vegetation
Lake Morphology;
Retention Time
Primary Production
$	? Temperature
Aquatic Macrophytes
Phyto plankton
Suspended
Sediment
Nuisance
Algae
Organic Matter
[ Algal
I Toxins
Taste and
Odor
Water Clarity
Light Availability
Recreation
Drinking Water Supply
Figure 2. Conceptual model linking increased nutrient concentrations to public health endpoints.
2.3 Risk Hypotheses
The EPA specified risk hypotheses for each of the selected assessment endpoints. Based
on a survey of available literature, the EPA concluded that increased concentrations of N and P
increase the risk to both ecological and human health (Figure 3). For aquatic life, the risk
hypotheses consist of the pathway in which increased nutrient concentrations increase
phytoplankton biomass (measured as chlorophyll a [Chi a]). Then, as phytoplankton biomass
increases, the relationship between zooplankton biomass and phytoplankton biomass changes
9

-------
so that increases in phytoplankton biomass are no longer associated with increases in
zooplankton biomass, and increases in primary production at the base of the lake food web are
not transferred to higher trophic levels. For the case of deepwater DO concentrations, increased
phytoplankton biomass increases organic matter in the lake, which when decomposed,
consumes DO (Walker 1979). The decreased concentrations of DO then affect lake aquatic life.
The risk hypotheses for recreation and drinking water source designated uses state that
increased nutrient concentrations increase the biovolume of cyanobacteria and concentrations
of microcystin.
Nutrients
Primary Production
Phytoplankton (Chi o)
Cyanobacteria
Organic material
Deep water
dissolved oxygen
Food quality
Biologica I CommunRjes
a n ktori^v
pi ^
Drinking water
source
Fish
Recreation
Aquatic life use
Figure 3. Simplified conceptual model showing pathways selected for analysis.
2.4 Analysis Plan
The analysis plan consists of acquiring appropriate data and estimating relationships
between phytoplankton biomass and each of the risk metrics as well as between N, P, and
phytoplankton biomass. The critical measurement in all these relationships is Chi a, which is
closely associated with phytoplankton biomass. Stressor-response analysis was applied to
available data to estimate relationships between nutrient concentrations and different risk
metrics. Because Chi a concentration is the critical parameter for all risk metrics, the EPA
developed different stressor-response models associating Chi a concentration with each of the
10

-------
risk metrics (i.e., zooplankton biomass, deepwater DO concentration, and MCs). The models
then yielded candidate criteria for Chi a corresponding to each of the risk metrics (and their
associated endpoints). N and P are estimated in field measurements as TN and TP, and so, the
EPA developed draft models relating TN and TP concentrations to Chi a concentrations that can
translate each of the different Chi a criteria into draft recommended TN and TP criteria.
Because different risk metrics have been identified for each of the three designated
uses, these risk metrics lead to the derivation of different draft recommended numeric nutrient
criteria. In general, a state's water quality criteria for any single lake would need to protect the
most sensitive use (i.e., the state should select the most stringent numeric nutrient criteria) (40
CFR 131.11(a)(1)).
3 Analysis
Because stressor-response analyses for each of the risk metrics differed substantially
from one another, most of this section is organized by models for the different risk metrics
zooplankton biomass, deepwater hypoxia, and microcystin - followed by models relating TN, TP,
and Chi a. Because the same data were used to fit each of these models, all the data used in the
analyses are discussed first.
3.1 Data
The EPA analyzed data collected in the NLA in the summers (May-September) of 2007
and 2012 to support the derivation of draft recommended numeric nutrient criteria. The NLA
data were collected from a random sample of lakes from the continental U.S. In 2007, lakes with
surface areas larger than 4 hectares and, in 2012, lakes larger than 1 hectare were selected from
the contiguous U.S. using a stratified random sampling design (US EPA 2012b). The final data set
was supplemented by a small number of hand-picked lakes identified as being less disturbed by
human activities (US EPA 2010b). The additional lakes were included to increase the number of
least-disturbed lakes for which data were available, and by helping ensure the full range of
conditions was sampled, data from the additional lakes was expected to improve the accuracy
of the estimated stressor-response relationships. The overall sampling design of the NLA was
synoptic, but 10% of sampled lakes were randomly selected and resampled on a different day
after the initial visit. The timing of the second visit varied among lakes, but on average, the
11

-------
second sample was collected approximately 46 days after the first. Approximately 20% of the
lakes were sampled in both 2007 and 2012. The sampling day of the year was recorded for each
visit and used in subsequent analyses to account for temporal changes in deepwater DO
concentration. Overall, data from approximately 1,800 different lakes are included in the data
set, but the specific number of samples used to estimate each stressor-response relationship
varies slightly based on data available at each lake. The specific number of samples is provided
in the subsequent discussion of each model.
During each visit to a selected lake, an extensive suite of abiotic and biological variables
was measured. Only brief details on sampling protocols are provided here regarding the
parameters used to derive these draft criteria; more extensive descriptions of sampling
methodologies are available in the NLA documentation (US EPA 2007, 2011). A sampling
location was established in open water at the deepest point of each lake (up to a maximum
depth of 50 meters [m]) or in the mid-point of reservoirs. In 2012, an additional sampling
location for collection of microcystin, algae, and Chi a data was established in the littoral zone
approximately 10 m away from a randomly selected point on the shoreline.
At the open water site, a vertical, depth-integrated methodology was used to collect a
water sample from the photic zone of the lake (to a maximum depth of 2 m). Multiple sample
draws were combined in a rinsed, 4-liter (L) cubitainer. When full, the cubitainer was gently
inverted to mix the water, and an aliquot was taken as the water chemistry sample. That
subsample was placed on ice and shipped overnight to the Willamette Research Station in
Corvallis, Oregon. A second aliquot was taken to use in characterizing the phytoplankton
community and was preserved with a small amount of Lugol's solution. A Secchi depth
measurement was also collected at this site. Two zooplankton samples were collected with
vertical tows for a cumulative tow length of 5 m using fine- (50-micrometer- [-nm-]) and coarse-
(150-nm-) mesh Wisconsin nets. In lakes at least 7 m deep, one 5-m deep tow was collected
with each mesh. In shallower lakes, vertical tows over shorter depths were combined to reach
the cumulative tow length of 5 m.
At the littoral zone site, two grab water samples were collected 0.3 m below the surface
where the lake was at least 1 m deep using a 2-L brown bottle. The first sample was split into
two subsamples: one subsample for quantifying algal toxin concentration and the second
subsample preserved with a small amount of Lugol's solution and used to characterize the
12

-------
phytoplankton community. The second grab sample collected with the 2-L bottle was used to
quantify Chi a concentration.
3.1.1 Biological Data
Phytoplankton biovolume from the field samples was measured in the laboratory.
Samples collected from both open water and littoral zone locations were examined by
taxonomists, who identified at least 400 natural algal units to species under l,000x
magnification. Observations were aggregated and abundance was calculated as cells per
milliliter. In each sample, the dimensions of the taxa that accounted for the largest proportion of
the observed assemblage were measured and used to estimate biovolume. Biovolumes of the
most abundant taxa were based on the average of measurements from at least 10 individuals,
while biovolumes of the less abundant taxa were based on somewhat fewer measurements. The
biovolume was reported as cubic micrometers per milliliter (nm3/mL) (US EPA 2012a), which
was converted to cubic millimeters per liter (mm3/L). Approximately 5% of the phytoplankton
samples were randomly selected and reidentified and measured by a second taxonomy
laboratory. These reidentified samples provided a basis for estimating laboratory measurement
error. Biovolume measurements were converted to biomass using a density of 1 gram per
milliliter (g/mL) (Holmes et al. 1969).
Zooplankton samples from the coarse- and fine-mesh net tows were processed
separately. In each sample, zooplankton specimens were examined and counted under 100-
l,000x magnification, in discrete subsamples until at least 400 individuals were identified. In the
coarse-mesh net samples, all taxa were identified and enumerated. In the fine-mesh net, only
"small" taxa were identified and enumerated (Cladocera less than 0.2 millimeters [mm] long,
copepods less than 0.6 mm long, rotifers, and nauplii). Zooplankton abundance was estimated
based on the volume of sampled lake water used to identify the targeted count of 400
individuals. Measurements of at least 20 individuals were collected for dominant taxa (i.e., taxa
encountered at least 40 times in the subsample); at least 10 individuals were measured for taxa
encountered from 20 to 40 times; and at least 5 individuals were measured for rare taxa
(encountered less than 20 times in the subsample). Zooplankton biomass estimates were based
on existing length and width relationships (Dumont et al. 1975, McCauley 1984, Lawrence et al.
1987). Estimates from the coarse- and fine-mesh samples were added to yield a single
zooplankton sample per lake visit.
13

-------
3.1.2	Chemical Data
For both 2007 and 2012 data, TN, nitrate-nitrite (NOx), ammonia, and TP concentrations;
true color, dissolved organic carbon (DOC) concentration, turbidity, and acid-neutralizing
capacity (ANC) were measured in the laboratory from the open water sample at prespecified
levels of precision and accuracy (US EPA 2012a). Typical laboratory methods included persulfate
digestion with colorimetric analysis for TN and TP, nephelometry for turbidity, comparison to a
calibrated color disk for true color, and automated acidimetric titration for ANC. To measure Chi
a concentration, 250 mL of lake water was pumped through a glass fiber filter in the field and
quantified in the laboratory to prespecified levels of precision and accuracy. Examples of lower
reporting limits include 20 ng/L for TN, 4 ng/L for TP, and 0.5 ng/L for Chi a.
Microcystin sample processing began with three sequential freeze/ thaw cycles to lyse
cyanobacteria (Loftin et al. 2008). Processed samples were filtered using 0.45 pim polyvinylidene
difluoride membrane syringe filters and stored frozen until analysis. The concentration of
microcystin in the filtered water sample was measured with an enzyme-linked immunosorbent
assay (ELISA) using an Abraxis kit for Microcystin-ADDA, which employs polyclonal antibodies
that are unique to microcystins and other similar compounds. The binding mechanism of the
Microcystin-ADDA assay is specific to the microcystins, nodularins, and their congeners;
therefore, results from that assay could include contributions from any compound within the
ADDA functional group (Fischer et al. 2001). The minimum reporting level for the assay was 0.1
ju,g/L as microcystin-LR.
3.1.3	Dissolved Oxygen and Temperature Profiles
At the deepest point of each lake (or in the midpoint of reservoirs), a multiparameter
water quality meter was used to measure profiles of DO concentrations, temperature, and pH at
a minimum of 1-m depth intervals (Figure 8). Profiles in lakes less than 3 m deep were sampled
at 0.5-m depth intervals. Water temperatures were converted to estimates of water density
(Jones and Harris 1992), and density gradient was estimated between all available depths below
0.5 m as the difference in density between two successive measurements divided by the
difference in the depths of the two measurements. Temperature gradients were computed with
the same approach. Samples collected in the uppermost 0.5 m were excluded to limit the effects
of surface warming on the gradient calculations.
14

-------
3.1.4 Mapped Data
Lake physical characteristics including lake surface area, geographic location (latitude
and longitude), elevation, lake catchment area, and lake perimeter were estimated from
mapped data. From these characteristics, the following composite variables were calculated: (1)
the drainage ratio, which is defined as the ratio of catchment area to lake surface area and
characterizes the degree to which the lake catchment influences the lake; (2) the shoreline
development, which is defined as the ratio between the perimeter of the lake and the perimeter
of a circle with the same area as the lake and characterizes the geometric complexity of the lake
shore; and (3) the lake geometry ratio, which is defined as area0 25/depth, or the ratio between
fetch and lake maximum depth, and has been shown to differentiate lakes that stratify
seasonally (low values of the geometry ratio) from lakes that are polymictic (Gorham and Boyce
1989, Stefan et al. 1996). Variables quantifying the mean annual precipitation and mean annual
air temperature at the lake location were extracted from 30-year averaged climatic data (Daly et
al. 2008).
3.2 Stressor-Response Models
3.2.1 Zooplankton Biomass
3.2.1.1 Statistical analysis
The EPA specified a Bayesian network model to estimate the relationship between
phytoplankton and zooplankton biomass (Figure 4). A "Bayesian network" provides a unified
framework for modeling the cascading relationships between different measurements and
propagates estimation errors and model uncertainty correctly throughout the model (Qian and
Miltner 2015; Yuan and Pollard 2018).
15

-------
Summer mean
phytoplankton biovolume
(s)
v
Summer mean
zooplankton abundance
- 0)
(7)
(4), (5)
	*	
Sample phytoplankton
biovolume
Summer mean
zooplankton biomass
1,2
Observed Ch a
Observed
zooplankton
abundance
Observed
zooplankton
biomass
Observed
phytoplankton
biovolume
Figure 4. Schematic of network of relationships for modeling zooplankton biomass. Gray-filled ovals:
available observations; other nodes: modeled parameters; numbers in parentheses refer to equation
numbers in the text.
The first set of relationships in the network estimated mean phytoplankton biovolume based on
both Chi a concentration and measurements of phytoplankton biovolume. The two
measurements provided independent estimates of phytoplankton biovolume, each with
different sources of error. Chi a is measured precisely from field samples, but the Chi a content
of phytoplankton can vary depending on environmental conditions and species composition
(Kasprzak et al. 2008), so that a measured Chi a concentration in one sample might indicate a
slightly different phytoplankton biovolume than the same Chi a measured in another sample.
Hence, Chi a concentration is modeled as being directly proportional to the true phytoplankton
biovolume in the sample (Psamp), but the constant of proportionality, b, (i.e., the Chi a content of
phytoplankton in a sample) is allowed to vary among samples. The log-transformed version of
this model equation is as follows:
\og(Chli) = log (PsamPii) + log Ofy)
logOfy) ~Normal(nb, ab)
(1)
(2)
where the value of b, for each sample, /', is drawn from a single log-normal distribution
characterized by a mean, /ut, and a standard deviation, ct/,. This multilevel expression of the
model equation allows the mean Chi a content of phytoplankton cells estimated for each
sample to vary, but imposes the constraint that estimates of phytoplankton Chi a content for
each sample must be drawn from a common log-normal distribution (Gelman and Hill 2007).
16

-------
Direct measurements of phytoplankton biovolume generally provide an unbiased
estimate of true phytoplankton biovolume. These direct measurements, however, are obtained
by summing contributions from measurements taken from many different individual
phytoplankton, each of which includes measurement error. Hence, the summed estimate of
total biovolume includes a substantial amount of measurement error. That measurement error
was explicitly modeled, and a second estimate of the true phytoplankton biovolume in a sample
was expressed as follows:
log(Pobs,i) ~Normal(\og(Psamp i), sx)	(3)
where P0bsj is the observed phytoplankton biovolume in sample /'. The standard deviation, si, of
the distribution is quantified using laboratory replicate measurements of phytoplankton
biovolume. Final model estimates of PSamP,i were then consistent with both Chi a and observed
phytoplankton biovolume, and by combining the two measurements, the accuracy of the final
estimate was maximized.
PsamP,i estimates phytoplankton biovolume within a single sample, but to model the
relationship between phytoplankton and zooplankton biomass, the EPA was interested in
seasonal mean values of phytoplankton biovolume for each lake. To estimate seasonal mean
phytoplankton biovolume, the EPA used model estimates of PSamP,i corresponding to
measurements collected at the same lake on different days and corresponding to
measurements collected on the same day in the littoral zone and in the middle of the lake to
provide a final estimate of the combined magnitude of temporal and sampling variability.
Seasonal mean phytoplankton biovolume (P) can then be expressed as follows:
log(Psamp,i) -Normal(log(P,ra), s2)	(4)
log(Pj) ~Normal(/uP, aP)	(5)
where j indexes different lakes and s2 is the standard deviation of the distribution representing
temporal and sampling variations in PSamP,i about the seasonal mean value. The distribution of
seasonal mean phytoplankton concentrations across all sites was then modeled as a log-normal
distribution with mean, /uP, and standard deviation, aP.
Zooplankton abundance (A) and biomass (Z) were modeled as increasing functions of
seasonal mean phytoplankton biovolume (or biomass, using the conversion factor of 1 g/mL).
Previous studies in oligotrophic lakes found that zooplankton biomass increased as a constant
17

-------
proportion of phytoplankton biomass (Rognerud and Kjellberg 1984, del Giorgio and Gasol
1995). That is, after log-transforming, the relationship between P and Z should approach the
following at low concentrations of P:
log(Z) = log(fc) + log (P)	(6)
where the slope of the relationship between log(Z) and log(P) approaches 1. In contrast, in
eutrophic lakes, minimal changes in Z were observed with changes in P, and the slope between
log(Z) and log(P) approached zero (Yuan and Pollard 2018). Those patterns guided the selection
of the following functional form for modeling the relationship between log(Z) and log(P):
log (Pj)CpN
[log(Zj)] = /i + f2log (Pj) - f3qlog
1 + exp ( ^
(7)
where, in general, E[.] indicates the expected value of the variable enclosed in the square
brackets. The coefficients fi,f2,f3, cp, and q were estimated from observations of Zj and the
estimated seasonal mean phytoplankton concentration, Pj, estimated in Equation (5). The slope
of this function approaches f2 at large values of P and approaches a slope of/2 +f3 at low values
of P. The a priori expectation for the value of f2 is zero to represent the weak relationship
between log(Z) and log(P) in eutrophic lakes. So, a prior distribution for/3 was defined as a
normal distribution centered at 1 with a standard deviation of 0.5. This distribution expresses
the prediction (stated above) that, at low levels of phytoplankton (oligotrophic lakes),
zooplankton biomass should increase as a constant proportion of phytoplankton biomass.
A similar model was specified for zooplankton abundance (A) as follows:
E[log(Aj)] = a1 + a2 log(P,) + a3 log [1 + exp (_'g(y c''j]	(8)
where the parameters, ai, a2, a3, and r were estimated from the data, and the third term on the
right side of the equation again introduces curvature in the fitted relationship. The change point
for zooplankton abundance, cp, was estimated as being the same as for zooplankton biomass
because of the strong influence of abundance on total biomass. In the case of zooplankton
abundance, no a priori assumptions about the slope of the relationship at high or low levels of
phytoplankton guided the choice of parameter values.
Observed values of zooplankton abundance and biomass were then related to the
estimated expected values as follows:
18

-------
logOWi) ~Normal(E[log(Ajli])],s3)	(9)
log(Zobsi) ~Normal(E[log(Zj[ii)],s4)	(10)
Similar to the model equations for phytoplankton, variability in the observations of zooplankton
abundance and biomass relative to estimated mean values were modeled as log-normal
distributions with standard deviations of s3 and s4. These error terms included contributions
from temporal, sampling, and measurement error.
Because the strength of the interaction of the zooplankton assemblage with benthic
resources was expected to differ between shallow and deep lakes (Benndorf et al. 2002,
Scheffer and van Nes 2007), different parameter values for ai, a2, a3,fi,f2,f3, and cp were
estimated for each of three classes of lakes defined by depth. The curvature parameters q and r
were fixed at 1. The number of lake classes was specified to balance between accounting for
differences in lake depth and maintaining enough samples within each class to estimate
relationships. Depth thresholds defining each class were selected to ensure that a similar
number of samples was assigned to each class, yielding the following thresholds: less than 3.2
m, 3.2-7.2 m, and more than 7.2 m.
All model equations were fit simultaneously to data collected at each lake, including
revisits on different days, littoral and mid-lake samples, and laboratory replicates of
phytoplankton measurements. Weakly informative priors were specified for all model
parameters except for/3 (Gelman 2006). Weakly informative prior distributions constrain
parameter estimates away from extreme values, while allowing the data to determine the
estimate for each parameter. All other statistical calculations were performed with R, an open-
source statistical modeling software (R Core Team 2017). Hierarchical Bayesian models were fit
using the rstan library which implements the No-U-Turn sampler, a variant of a Hamiltonian
Monte Carlo sampling approach (Duane et al. 1987, Stan Development Team 2016).
3.2.1.2 Results
Data collected at a total of 1,127 lakes were available for analysis, with approximately
380 lakes assigned to each depth class. Estimated mean phytoplankton biovolume within each
sample was much more strongly associated with Chi a concentration than with measured
phytoplankton biovolume, because of the high measurement error associated with measured
phytoplankton biovolume (Figure 5). Variance in laboratory replicate measurements accounted
19

-------
for 38% of the total variance in observed phytoplankton biovolume, a percentage that was
somewhat lower than the variance attributed to differences in seasonal means among sites
(56%) and much higher than the percentage of variance attributed to temporal and sampling
variability (6%). So, temporal and sampling variability accounted for only a small proportion of
the variance in observations of phytoplankton biovolume.
Figure 5. Relationships between measured biovolume, Chi a, and estimated mean phytoplankton
biovolume. Solid lines: 1:1 relationship.
Estimated relationships between phytoplankton biomass (as quantified by Chi a) and
zooplankton abundance and biomass matched trends observed in the data (see Figure 6 for an
example for lakes between 3.2 and 4.7 m deep in left and middle panels, respectively). The
relationship between zooplankton biomass and phytoplankton biomass also was consistent with
the initial assumption that, in oligotrophic lakes with low levels of phytoplankton biomass, the
slope approached 1, and in eutrophic lakes with high levels of phytoplankton biomass, the slope
approached zero (right panel, Figure 6).
The models show the gradual change in the shape of the biomass pyramid along the
eutrophication gradient. In oligotrophic lakes, the slope of the relationship between zooplankton
and phytoplankton biomass is near 1, indicating that small increases in phytoplankton biomass are
reflected in a proportional increase in zooplankton biomass. As Chi a increases, however, the slope
decreases, and the increase in zooplankton biomass per unit of increase in phytoplankton
biomass approaches zero. In eutrophic lakes, increases in phytoplankton biomass do not result
in comparable changes in zooplankton biomass. These changes along the eutrophication
gradient are consistent with other similar studies, as reviewed in Yuan and Pollard (2018).
O
O
o
o
o
Mean biovolume (mm3/L)
Mean biovolume (mm3/L)
20

-------
10
Chlafrg/L)
100
10
Chlafrg/L)
100
0.1
10
Chi a (>g/L)
100
1000
Figure 6. Estimated relationships between zooplankton and Chi a for lakes > 7.2 m deep. Left panel: Chi a
vs. zooplankton abundance; middle panel: Chi a vs. zooplankton biomass; right panel: Chi a vs. slope of
the relationship between zooplankton biomass and Chi a. Solid lines: mean relationships; shaded areas
(left and middle panels): 80% credible intervals about mean relationship; dashed lines (right panel): 50%
credible intervals about mean relationship; open circles (left and right panels): average of five samples
nearest the indicated Chi a concentration; dotted horizontal line (right panel): one example value of
threshold for deriving a Chi a criterion.
3.2.1.3 Chi a criterion derivation
Calculating candidate criteria for Chi a based on this response requires the specification
of two parametersthe value of the slope between log(Z) and log(P) and the credible interval
(i.e., the Bayesian analog to a confidence interval). The selected value of the slope identifies the
point at which food web connectivity between phytoplankton primary productivity and
zooplankton grazing is likely too low to control excess primary productivity in the lake. A
threshold slope of zero is the limit beyond which additional increases in phytoplankton biomass
are not converted to zooplankton biomass, and that slope is the lowest target for the threshold
slope. Higher threshold slopes might be selected for oligotrophic lakes in which a higher
proportion of phytoplankton is expected to be consumed by zooplankton. Graphically, this
threshold defines the horizontal line on which the Chi a criterion will be based (see the dotted
line in the right panel of Figure 6).
The selection of a threshold slope between log(Z) and log(P) (i.e., the targeted
condition) can also be informed by computing the predicted increase in zooplankton biomass
associated with an increase in phytoplankton biomass. More specifically, the change in
zooplankton biomass can be expressed as follows:
where m is the slope between log(Z) and log(P), P2 and P5 are two different phytoplankton
biomasses, and Z2 and Z5 are the corresponding zooplankton biomasses. So, when the slope
(11)
21

-------
between log(Z) and log(P) in a particular lake is 0.1, the predicted increase in zooplankton
biomass with a doubling of phytoplankton biomass is 2ai, or 1.07. That is, only a 7% increase in
zooplankton biomass is expected when phytoplankton biomass is doubled. Table 2 shows other
predicted increases in zooplankton biomass.
Table 2. Predicted proportional increase in zooplankton biomass with different increases in phytoplankton
biomass (P2/P1) and different slopes, m, between log(Z) and log(P)


P2/P1

m
1.5
2.0
3.0
0.1
1.04
1.07
1.12
0.2
1.08
1.15
1.25
0.3
1.13
1.23
1.39
Credible intervals express the statistical uncertainty about the position of the mean
relationship and are directly comparable to confidence intervals used in frequentist statistics.
The mean relationship between the slope and Chi a represents the best estimate for the slope
of the stressor-response relationship; however, a lower credible interval provides additional
assurance that the calculated criterion is protective, given the data and model uncertainty. That
is, more protective criteria are based on lower percentiles of the credible interval. For example,
selecting the 25th credible interval implies that 25% of estimated slopes, given the data, are less
than the selected threshold. That is, at the calculated criterion value, a lake has a 75% chance of
achieving the targeted condition. In contrast, selecting the 10th credible interval implies that a
lake has a 90% chance of achieving the targeted condition. In statistical hypothesis testing,
convention suggests that p-values of 1% or 5% are statistically significant results, which can also
inform the selection of the credible interval. Selection of the value of the lower credible interval
as the basis for the criteria is ultimately a management decision, and a range of credible
intervals from 1% to 25% is provided in the associated interactive tool (see below). Illustrative
criteria for Chi a for different combinations of management decisions are shown in Table 3
(slope threshold = 0 is shown in Figure 6). The interactive tool, which uses posterior simulation
with the estimated parameter distributions, computes candidate criteria for different
combinations of the slope threshold and the credible interval (https://chl-zooplankton-
prod.app.cloud.gov). With this tool, a user can specify the value of the slope between log(Z) and
log(P), lake depth, and the credible interval with sliders, and the associated criteria and stressor-
response relationship are updated to reflect those selections.
22

-------
Table 3. Illustrative Chi a criteria (ng/L) for different credible intervals and a threshold value of 0 for A(log
Z)/A(log P). Values shown for each lake depth class.


Depth class

Credible
interval
< 3.2 m
3.2-7.2 m
> 7.2 m
10%
41
22
13
25%
48
36
16
3.2.2 Deepwater Hypoxia
The EPA specified a model for deepwater DO that represents the temporal decrease in
DO during summer stratification, while accounting for differences among lakes in eutrophication
status, depth, and DOC concentrations (Yuan and Jones 2020a).
3.2.2.1 Data
The EPA first restricted analysis to data collected from seasonally stratified lakes
because hypoxic and anoxic conditions occur more consistently during stratified conditions.
Lakes were identified that were likely to be seasonally stratified by computing the lake geometry
ratio. This metric approximates the relative effects of lake fetch and depth on stability of
stratification, and lakes with a geometry ratio less than 3 m"05 exhibit seasonal stratification
(Gorham and Boyce 1989). Therefore, the EPA restricted NLA data to lakes with geometry ratios
less than that threshold. Lakes likely to be dimictic (i.e., mixing fully in the spring and in the fall)
were also identified based on latitude and elevation. This classification approach adjusts the
lake latitude by elevation, and then identifies lakes with adjusted latitudes greater than 40 N as
dimictic (Figure 7) (Lewis 1983). Finally, data were restricted to samples in which temperature
profiles exhibited evidence of stratification (defined as a temperature gradient of at least 1
degree Celsius per meter [C/m]).
23

-------
Figure 7. Lakes designated as seasonally stratified dimictic lakes from the NLA data set.
Mean deepwater DO concentrations (DOm) in the selected NLA lakes were computed
from temperature and DO profiles. First, measurements collected at depths less than or equal to
0.5 m were excluded to minimize the effects of surface warming. In some profiles, duplicate
measurements of DO and/or temperature were collected at each depth, and in these cases, the
average was used in computations. The EPA used only profiles with measurements collected
from at least half of the possible 1-m increments in the final analysis.
The upper boundary of the metalimnion was identified as the shallowest depth at which
the temperature gradient exceeded 1 C/m (excluding the surface layer) (Figure 8) (Wetzel
2001). DOm for each lake profile was computed as the mean of DO measurements estimated at
all 1-m increments deeper than the upper boundary of the metalimnion. That estimate of DOm
necessarily includes some measurements in the metalimnion, which might increase the
estimates of DOm relative to studies that can focus only on the hypolimnion. In the NLA data set,
the upper boundary of the metalimnion could be determined for most profiles. In contrast,
many lakes in the NLA data set were too shallow to maintain a hypolimnion with small vertical
temperature gradients (Jones et al. 2011), and therefore, no approach for consistently defining
the hypolimnion for all lakes was available (Quinlan et al. 2005). Furthermore, inclusion of the
metalimnion was consistent with the assumption that taxa can use this transitional region as a
refuge from warmer temperatures in the mixed layer (Klumb et al. 2004). The depth of water
below the thermocline was computed as the difference between the maximum depth recorded
for each lake and the mean depth of the upper boundary of the metalimnion. Chi a and DOC
24

-------
measurements from each lake were also used in the analysis. Prior to statistical analysis, all
measurements were standardized by subtracting their overall mean values and dividing by the
standard deviation. This standardization had no effect on the final model results, but helped the
Bayesian models converge more efficiently (Gelman and Hill 2007).
Q. O
0 t-
Q
i	1	1	r~
5 10 15 20
Temperature (C)

1	1	1	1	1	1	
0 1 2 3 4 5
Temperature gradient (C/m)
~i	1	r~
2 4 6
DO (mg/L)
Figure 8. Illustrative examples of depth profiles of temperature, temperature gradient, and DO. Dashed
horizontal line: estimated depth of the bottom of the epilimnion.
3.2.2.2 Statistical analysis
The EPA modeled the decrease in DOm as a linear function, an approximation that is
appropriate for DOm concentrations higher than approximately 2 milligrams per liter (mg/L)
(Burns 1995). This threshold reflects experimental evidence indicating that the rate of decrease
of hypolimnetic DO is constant at relatively high ambient concentrations of DO, but can be
affected by DO concentrations near zero (Cornett and Rigler 1984). The linearly decreasing
function also precludes the possibility of episodic mixing events that transport DO from shallow
waters to deeper depths of the lake. In some lakes, those mixing events are rare, but in other
lakes, they might occur frequently. In the latter group of lakes, the model predicts DOm during
extended periods of still weather, and the associated criteria would protect aquatic life in those
scenarios. Below, the statistical model is first described followed by a description of the
approach for addressing DOm measurements less than 2 mg/L.
25

-------
(13)
(14)
(12)
DO
Sampling
day (t,)
Depth below
thermocline (D)
Oxygen demand
(VOD)
Air temperature
(Temp)
First day of
stratification (f,
Eutrophication
(Chi a)
Dissolved organic
carbon (DOC)
Figure 9. Schematic of hypoxia model. Numbers in parentheses refer to equation numbers in the text.
NLA data were fit to the following model equation:
E[DOmi] = DO0J[q + VODj[q(ti - t0,fcp])	U2)
where DO0,j[i] is the value of DOm at the start of spring stratification in lake j corresponding to
sample /', and volumetric oxygen demand (VODj) is the net imbalance in the volumetric oxygen
budget for lake j, expressed as mg/L/day of DO (Burns 1995). That is, VOD estimates the rate of
decrease in DOm per day. t, is the date that sample i is collected, and t0/kw is the date of the
beginning of stratification for the lake-year k corresponding to sample /'. Observed values of
DOm,i were modeled as being normally distributed about the expected value, with a standard
deviation of Oi.
The first day of stratification (t0) was not measured for any of the lakes, and the precise
day on which stratification occurs for a given lake and year depends on local wind speeds,
temperatures, and lake morphology (Cahill et al. 2005). Previous work in northern temperate
dimictic lakes found that the first day of stratification could be modeled as a function of mean
annual temperature (Demers and Kalff 1993), so the EPA specified the following relationship for
to:
to,k = b1 + bzTempm + ek	(13)
where Tempj[ki is the mean annual air temperature at the location for lake j corresponding to
lake-year k, and bi and b2 are coefficients that are fit to the data. The published relationship in
Demers and Kalff (1993) provided initial estimates for bi and b2, which were used to specify
prior distributions for the two parameters. The error term ek is included in the model because
26

-------
the first day of stratification varies substantially in different years for a given lake because of
differences in weather. Data published by Demers and Kalff (1993) indicated that the standard
deviation of residual error for this relationship was approximately 12 days, so this value was
used to specify the prior distribution for the standard deviation of ei<.
The initial concentration of DO at the time of stratification, DO0, was also not measured
for any of the lakes. Deepwater temperatures in many dimictic lakes are determined by
temperatures prior to initiation of stratification (Hondzo and Stefan 1993), and so, deepwater
lake temperatures at the time of stratification were approximated as the minimum annual air
temperature at the lake location. Then, the saturated DO concentration at the minimum annual
air temperature provided an estimate for DO0. Minimum air temperatures less than 4 degrees
Celsius (C) were set to 4 C, corresponding to water temperatures when the lake surface begins
to freeze (Demers and Kalff 1993).
Lake trophic status affects VOD because increased phytoplankton production in the
epilimnion increases the quantity of organic material available for decomposition in the
hypolimnion and in lake sediments (Hutchinson 1938). In many lakes, allochthonous sources
also provide organic matter that fuels bacterial respiration and depletes oxygen in deep lake
waters (Pace et al. 2004, Kritzberg et al. 2004). VOD has also been observed to decrease with
increasing hypolimnion depth, a phenomenon attributed to a weaker overall influence of
sediment oxygen demand as the volume of the hypolimnion increases (Cornett and Rigler 1980,
Muller et al. 2012). Based on these mechanisms, the EPA modeled VOD as a linear function of
the long-term mean Chi a concentration and depth below the thermocline in the lake. To
account for the effect of allochthonous organic matter, DOC was also included as a third
predictor variable for VOD (Hanson et al. 2003, Cole et al. 2011). The model equation for VOD
can then be written as follows:
E[VODj] = dt + d2 log(ChlmnJ) + d3DmnJ + d3log (DOCmnJ)	(14)
where di, d2, d3, and d4 are model coefficients estimated from the data; log(Chlmn,j) is the long-
term mean of the log-transformed Chi a concentration lake, j; Dmnj is the mean depth of the lake
below the thermocline; and log(DOCmn,j) is the seasonal mean of log-transformed DOC
concentration in the lake. Variability in VOD across individual lakes about the mean value
estimated from the predictor variables was modeled as a normal distribution. Because Chi a
concentrations can vary substantially over the summer in a lake, the modeling approach used
27

-------
with the zooplankton model provided a distribution of possible long-term mean log(Chlmn)
values for each lake, given one or more instantaneous measurements of Chi a concentration.
More specifically, seasonal mean log(Chlmn) values for different lakes were modeled as a normal
distribution as follows:
log(ChlmnJ) ~NormalQichl,schlil)	(15)
Then, individual log-transformed measurements from each lake were assumed to be drawn
from a normal distribution with a mean value equal to the long-term mean as follows:
\og(Chli) -Normal^log (ChlmnJ[q), schl2)	(16)
where Chi, is the Chi a concentration measured in sample, /', associated with the mean log(Chlmn)
concentration in lake j[i], (Note that Equations (15) and (16) are not shown in Figure 9.) Within-
year variability of DOC and depth below the thermocline were substantially less for Chi a, so
long-term means for each of those parameters were estimated as the mean value of all available
data for each lake.
As noted earlier, DOm approaches zero asymptotically over time and modeling that
relationship with the linear model described above would introduce biases to the model. To
account for the asymptotic relationship, the EPA modeled samples with DOm less than 2 mg/L
with methods used for measurements that are below a known detection limit. That is, the
samples were modeled as if their "true" DOm values were unknown but their maximum values
were 2 mg/L (Gelman and Hill 2007). This approach retained some information inherent in a
sample with DOm less than 2 mg/L (i.e., Chi a, lake depth, DOC, and sampling day are consistent
with low DOm), but allowed the use of linear relationships in the model to estimate the rate of
DO depletion. More specifically, the model fits a linear trend in time to DOm observed from lakes
with similar Chi a, DOC, and depth. By assuming that measurements of DOm less than 2 mg/L are
unknown, the estimates of the linear relationship are more strongly determined by the higher
DOm concentrations, and samples with DOm less than 2 mg/L exert a weak influence that is still
consistent with the overall relationship. Retaining samples with DOm less than 2 mg/L in the
model prevents biases that would be introduced by considering only lakes with relatively high
DOm.
28

-------
3.2.2.3 Results
A total of 477 samples collected at 381 lakes were available for analysis. D0m
concentrations in 165 samples were less than 2 mg/L and were modeled as unknown values that
were less than 2 mg/L The asymptotic relationship can be seen in the plot of Chi a versus DOm
(Figure 10), in which DOm decreases steadily up to a Chi a concentration of about 4 ng/L. At
higher Chi a concentrations, the magnitude of the slope of the relationship between DOm and
Chi a decreases and approaches zero.
m
_i
U)
E
o0

o
Q
m
co
o
1
10
Chi a Oig/L)
Figure 10. Chi a vs. DOm. DOm values. Gray-filled circles: values < 2 mg/L; solid line: non para metric fit to
the data shown to highlight asymptotic relationship.
The majority of the estimates for the first day of stratification ranged from day 30 to day
120 (Figure 11). In most lakes, the Demers and Kalff (1993) estimate for the first day of
stratification was later than the value of t0 estimated by the model. This systematic difference is
consistent with the fact that most of the lakes considered in Demers and Kalff (1993) were
located north of the mean latitudinal location of the NLA lakes. The strong association between
the Demers and Kalff (1993) estimates and the current estimates indicates that the overall
formulation of the model, in which stratification day is a function of mean annual temperature,
is valid.
29

-------
	1	1	1	1	1	r
80 100 120 140 160 180
D & K predicted stratification day
Figure 11. Demers and Kalff (1993) predicted stratification day vs. model mean estimate. Solid line: the 1:1
relationship.
Relationships estimated between DOm and different predictors were consistent with the
hypothesized effects of each of the predictors (Figure 12). DOm decreased strongly with
increases in DOC and Chi a, reflecting the increased organic material available in lakes with high
concentration of the two parameters. Conversely, DOm increased with increasing depth below
the thermocline, consistent with observations in other studies. Substantial uncertainty is
associated with the relationship between DOm and day of the year, reflecting the inherent
uncertainty in estimating the first day of stratification for different lakes.
30

-------
The root mean square (RMS) error on model predictions for samples with D0m higher
than 2 mg/L was 1.5 mg/L. Slightly greater residual variability in the observations about the
mean predictions were observed at high values of DOm (Figure 13).
1	10	140 160 180 200 220 240 260
Chi a (jag/L)	Day of the year
Figure 12. Relationships between individual predictors and DOm, holding other variables fixed at their
mean values. Solid line: mean relationship; gray shading: 90% credible intervals.
DOC (mg/L)
i	1	1	1	1	1	r
0 10 20 30 40 50 60
Depth below thermocline (m)
CD
_i
U)
E
o o
 O
O ^
Q
"S
2
8> 
-Q
O
c>
CN
0
2
6
8
10
12
4
Predicted DOm (mg/L)
Figure 13. Model predicted DOm vs. observed DOm. Open circles: individual samples; solid line: 1:1
relationship.
31

-------
The statistical model described for DOm is consistent with the mechanisms of DO
depletion in the deep waters of a lake, in which available DO below the thermocline is
progressively depleted after the initiation of spring stratification. The estimated effects of
eutrophication, DOC, and lake depth on the rate of oxygen depletion were consistent with
trends observed in other studies.
3.2.2.4 Chi a criteria derivation
As described earlier, warm temperatures in the shallow mixed layer of a lake act
together with deepwater hypoxia to constrain the available habitat for cool- and cold-water
taxa. Therefore, to derive criteria based on deepwater hypoxia, estimates of changes in water
temperature over the course of the summer are required to identify periods of time during
which mixed layer temperatures are too high for different taxa. Those periods of time then
determine when deepwater DO concentrations need to be sufficiently high to support different
organisms.
Water temperature in the lake mixed layer depends on a variety of factors, including the
local climate, solar insolation, lake morphology, and the day of the year (increasing in the spring
and summer and decreasing in the fall). To identify temperatures in different lakes that were
likely to limit available habitat for different fish, the EPA first developed models to predict
temperature in the shallow, mixed layer of different lakes. NLA data collected at all lakes in the
conterminous U.S. were used to fit the model. At each lake, maximum temperature (excluding
the top 0.5 m of the surface layer) observed in vertical profiles collected in each lake were
modeled as a function of lake geographic location, elevation, and sampling day of the year with
a generalized additive model (Wood 2006) of the following form:
E[Ti\ = fi + f2Elevm + siydaji, df = 8) + s(Latm, Lonm, df = 20)	(17)
where E[T,] is the expected value of the maximum temperature in the lake observed in sample /'.
EleVj[i]\s the elevation of the lake, j, corresponding to sample /'. The variable yday, is the day of
the year that the sample was collected, and Latjtf and Lonjy are the latitude and longitude of the
lake. The relationship between temperature and elevation was modeled as a simple linear
relationship, characterized by two regression coefficients, fi, and f2. Relationships between lake
temperature and sampling day and between lake temperature and location were modeled as
nonparametric splines, represented in Equation 17 as s(.), with the maximum degrees of
32

-------
freedom, df, as indicated. Observed values of T, were assumed to be normally distributed about
the modeled expected value.
Lake temperature generally decreased with increased latitude, as would be expected
(Figure 14), but deviations from that latitudinal pattern were observed on the west coast of the
U.S., where lake temperatures were substantially lower than lakes at a similar latitude in the
eastern U.S. This trend likely arises from the moderating influence of the coastal waters on air
temperatures. Lake temperatures in eastern Texas and Louisiana were warmer than lake
temperatures at the same latitudes elsewhere. Lake temperatures decreased with elevation, as
expected, and exhibited a unimodal pattern with sampling day, with maximum temperatures
occurring on average on Day 204, or July 22 (Figure 15). Overall, the model predicted lake
temperature with an RMS error of 1.9 C.
Figure 14. Contours of modeled mean lake temperature computed at the overall mean elevation and
mean sampling day.
33

-------
~J CO
E
0)
" o
a> (N
(D
CD
CM
150
200	250
Day of year
0
1000 2000 3000
Elevation (m)
Figure 15. Relationship between lake temperature and sampling day (left panel) and elevation (right
panel). Variables that are not plotted are fixed at their mean values. Gray shading: 90% confidence
intervals; solid lines shows the mean relationships.
The pattern of temperature changes with time (Figure 15) provides insight into the
critical period during which the severity of deepwater hypoxia can influence aquatic life in lakes.
For most lakes, mixed layer temperatures increase in the spring and exceed critical
temperatures for different species, at which point cool- and cold-water obligate species must
move to deeper depths. Then, in the fall, decreasing mixed layer temperatures allow those
species to move back to shallower waters. Models for DOm indicate that, in dimictic lakes after
the onset of spring stratification, DOm decreases monotonically over time until fall turnover
(Figure 12). Therefore, the length of time between spring stratification and when mixed layer
temperatures decrease below the critical temperature thresholds in the fall is a key factor for
deriving a protective Chi a criterion.
The EPA used existing temperature thresholds defined for cool- and cold-water fish as
examples of critical mixed layer temperatures (Coker et al. 2001). For cool-water species, the
EPA identified a critical temperature of 24 C. Walleye, striped bass, and yellow perch are
examples of lake fish that are members of that group (McMahon et al. 1984). For cold-water
species, the EPA identified a critical temperature of 18 C. Lake trout is one example of a cold-
water obligate species (Marcus et al. 1984). Then, given a lake's location and elevation, the lake
temperature model predicts the day of the year that the mixed layer temperature would
decrease below the critical temperatures. For cool-water species, mixed layer temperatures
decreased below the critical temperature of 24 C on days 210-260 (Figure 16), taking into
account the fact that the dimictic lakes considered in this analysis are located in the northern
34

-------
half of the country (see Figure 7). Lakes in which mixed layer temperatures increased above 24
C at some point during the year were predominantly located in the eastern U.S., as high
elevations and climate in the western U.S. moderate lake temperatures. For cold-water species,
mixed layer temperatures decreased below the critical temperature of 18 C on days 220-280
(Figure 17). Temperatures in many lakes in the southeast part of the U.S. rarely decrease below
the critical threshold in the summer, but those lakes also generally do not harbor cold-water
fish.
220-
240-
270-
oo(
250'
Figure 16. Days of the year that mixed layer temperatures decrease below the critical temperature for
cool-water species. Small dots: lakes in which mixed layer temperatures never exceed 24 C.
250'
260
,61
230;

Figure 17. Days of the year that mixed layer temperatures decrease below the critical temperature for
cold-water species. Small dots: lakes in which mixed layer temperatures do not decrease below 18 C
during the summer; contours: effects of large differences in elevation across lakes in the western U.S.
35

-------
Draft criterion values for Chi a are calculated from the model equation for DOm,
rewritten here:
D0m = DO0 + [di + d2 log{Chlmn) + d3D + dA log(DOCmn)](t - t0)	(18)
Deepwater DO concentrations depend not only on Chi a concentration, but also on the depth of
the lake below the thermocline (D), DOC concentration (DOCmn), and length of time that has
elapsed since the establishment of stratification (t -10). A procedure for computing the day of
the year, cr/t, at which mixed layer habitat is cool enough for different species to move to
shallower water is also described above, highlighting the influence of lake location and elevation
as additional factors to consider. Based on these models, Chi a criteria for different lakes vary
considerably depending on each lake's specific characteristics.
Prior to calculating a Chi a criterion, a threshold value for DOm must be selected. Existing
EPA recommendations specify that the mean minimum DO concentration should be at least 5
mg/L to support cold-water fish (US EPA 1986). This threshold is also consistent with DO
concentrations that fish have been observed to avoid in field studies (Coutant 1985, Plumb and
Blanchfield 2009). A thin layer of cool water with sufficient DO provides a critical refuge for fish
during the warmest periods of the year, and fish have been observed to seek out those cool
water refuges. Observations of fish in warm lakes during the summer have indicated that they
will congregate in cold water refuges as shallow as 30 cm (Coutant and Carroll 1980, Snucins and
Gunn 1995, Baird and Krueger 2003, Mackenzie-Grieve and Post 2006). Hence, maintaining a DO
concentration of at least 5 mg/L at a depth of 30 cm below the thermocline can provide a
sufficient refuge for certain fish species and be protective of aquatic life. To convert this
condition to a value of DOm, the EPA considered a simplified case in which DO linearly decreases
from saturated conditions above the thermocline (DO = 8.4 mg/L at 24 C) to a concentration of
zero at some deeper depth (Figure 18). The linear decrease in DO is consistent with a steady-
state solution of the diffusion equation, assuming a constant eddy diffusivity (Stefan et al. 1995).
Based on this DO profile and the requirement that DO is 5 mg/L at 30 cm below the thermocline,
an illustrative threshold value for DOm can be computed as 1.6 mg/L for a lake that is 2 m deep
below the thermocline. That is, when the temperature profile is as depicted in Figure 18, depth-
averaged DO computed for the water column below the thermocline is 1.6 mg/L. Other
thresholds for DOm specific to different species of fish and different depths can also be
36

-------
calculated. For example, the threshold value for D0m for a lake that is 10 m deep below the
thermocline would be 0.3 mg/L.
i	1	1	1	1
0 2 4 6 8
DO (mg/L)
Figure 18. Simplified DO profile used to compute threshold for DOm. Open circle: the targeted condition of
DO at 5 mg/L, 30 cm below the thermocline.
The influence of different factors on Chi a criterion can be visualized by computing
criteria at median values of all covariates and then examining changes in criteria that occur with
the change in a single covariate. The relationship between Chi a and DOm at median values for
all other covariates are shown as solid lines in each panel of Figure 19. Lakes in which covariate
values differ from the medians of the data set cause changes in the candidate Chi a criteria. For
cool-water species, the median number of days between spring stratification and release of the
temperature constraint in the mixed layer was 135 days. The 75th percentile of this day range,
corresponding to lakes in warmer climates, was 151 days, whereas the 25th percentile,
corresponding to lakes in cooler climates, was 116 days. When the critical window for
maintaining sufficient DO in the deeper waters decreases to 116 days, the corresponding Chi a
criterion increases to 11 ng/L, whereas in lakes in which the critical window is 151 days long, the
Chi a criterion is 2 ng/L (left panel, Figure 19).
E
Thermocline
m
37

-------
Days since stratification
E
O
o
o
DOC

E
o
O
CM
O
1
oo
Depth below thermocline

E
o
Q
CM
O
Figure 19. Effects of other predictors on Chi a criteria. Solid lines: relationship between Chi a and D0m at
median values for all other variables; dashed line: D0m = 0.3 mg/L; dotted lines: 25th and 75th percentiles
of days elapsed since stratification (left panel), 25th and 75th percentiles of mean DOC concentrations
(middle panel), and depth below thermocline of 4 m and 20 m (right panel).
Similar ranges of criteria can be calculated for changes in DOC and the lake depth below
the thermocline. The median concentration of DOC in the available data was 5 mg/L, but in lakes
in which DOC is 3 mg/L (the 25th percentile of observed DOC in the data), the Chi a criterion
increases to 8 ng/L; and in lakes in which DOC is 7 mg/L (the 75th percentile), the Chi a criterion
decreases to 2 ng/L (middle panel, Figure 19). Finally, the median lake depth below the
thermocline was 9 m. In a deeper lake, with 20 m of water below the thermocline, the Chi a
criterion increases to 7 ng/L; but in a shallower lake, with only 4 m of water below the
thermocline, the Chi a criterion decreases to 3 ng/L (right panel, Figure 19).
To better understand the possible range of criteria, the EPA computed draft Chi a criteria
for each of the dimictic lakes sampled in the NLA. Because those lakes represent a random sample
of the population of lakes in the U.S., the resulting Chi a criteria are a representative distribution
of criteria, providing insight into likely criteria for different types of lakes. For dimictic lakes
harboring cool-water species, the median Chi a criteria is 3.4 ng/L, and the range defined by the
25th and 75th percentiles is 1.3-10.6 ng/L. For lakes harboring cold-water species, the median
Chi a criterion is 1.8 ng/L, with a range of possible values extending from 1 to 7.6 ng/L.
In states where measurements of profiles of DO are available, these data can be readily
modeled in conjunction with the national data (see Appendix B). In the example shown in
Appendix B, modeling temporally resolved DO profiles from one state with the national data
improved the precision of estimates of the first day of stratification. Because of this
improvement in model precision, the results of the combined state-national model are provided
in the interactive criterion derivation tool.
38

-------
The interactive tool used for estimating candidate Chi a criteria is provided at
https://chl-hypoxia-prod.app.cloud.gov. With this tool, the user can specify lake physical
characteristics that influence the relationship between Chi a and DOm as well as management
decisions about targeted conditions that affect the criterion. Lake physical characteristics that
are specified include the lake location (latitude and longitude) and lake elevation. That
information is converted to an estimate of mean annual air temperature and, coupled with the
model results, these data provide an estimate of the date of spring stratification. Other lake
physical characteristics that are specified are lake depth below the thermocline and average
lake DOC concentration, factors that influence DOm.
Water quality management decisions that influence the calculated criterion include the
critical maximum temperature for fish species in the lake, the threshold DO concentration, the
depth of the summer refugia, and the lower credible interval. The critical maximum
temperature for fish species in the lake is used to calculate the average day of the year that
temperature constraints are released in the epilimnion. That is, the annual temperature model
(Figure 15) is used to identify the date that fish can potentially move to oxygen-rich shallower
waters. The threshold DO concentration for the fish (e.g., a DO concentration of 5 mg/L for cold-
water fish) and the desired minimum thickness of the refugia (e.g., 30 cm) are used to compute
the targeted condition for DOm. That targeted value of DOm is the minimum concentration
required on the days prior to the release of temperature constraints. Credible interval
selections, as with other criteria, provide additional assurance that the calculated criterion is
protective, based on the data and model uncertainty. For example, selecting the 25th credible
interval implies that, at the estimated Chi a criterion, only 25% of predicted mean values of DOm,
based on the data, were less than targeted value. In statistical hypothesis testing, convention
suggests that p-values of 1% or 5% are statistically significant results, which can also inform the
selection of the credible interval, but selection of the value of the lower credible interval as the
basis for the criterion is ultimately a water quality management decision.
The interactive tool uses posterior simulation with model parameter distributions to
predict the DOm on the critical day prior to a release from temperature constraints in the surface
layer for different Chi a concentrations. These model results can be used to help derive criteria
for a specified threshold DOm. Samples with covariate values similar to those selected by the
user are highlighted in the provided plots in the app.
39

-------
3.2.3 Microcystin Concentration
3.2.3.1 Statistical analysis
A network of relationships can be specified that reflects current understanding of the
linkage between lake eutrophication (as represented by Chi o) and increased concentrations of
microcystin in individual samples (Figure 20). At the bottom of the diagram, cyanobacterial
biovolume is directly associated with MC. Cyanobacterial biovolume is then expressed as the
product of total phytoplankton biovolume and the proportion of the biovolume that is
cyanobacteria (i.e., the relative biovolume of cyanobacteria), which clarifies the nature of the
relationship between Chi a and cyanobacterial biovolume. More specifically, Chi a is directly
proportional to phytoplankton biovolume (repeating the relationship used in the zooplankton
model) (Kasprzak et al. 2008), and, as Chi a increases, the relative biovolume of cyanobacteria
has been observed to increase (Downing et al. 2001).
(20)
(19)
(23)
(24)
Microcystin
Chi a
Phytoplankton biovolume
Cyanobacterial biovolume
Cyanobacterial relative
biovolume
Figure 20. Schematic showing relationship between different variables predicting MC. Numbers in
parentheses: refer to equation numbers in the text.
Each of the relationships in the network described above is expressed mathematically in
the Bayesian network. First, phytoplankton biovolume, P is modeled as being directly
proportional to Chi a concentration (Chi), in sample /':
Pi = kciChli	(19)
The reciprocal of the parameter kc,i is the average amount of Chi a per unit biovolume of
phytoplankton. Because the Chi a content of phytoplankton can vary with environmental
conditions and assemblage composition, different values of this parameter are estimated for
40

-------
each sample, /'. The overall distribution of the set of values for kcj is assumed to be log-normal
with a mean value of /Xk and a standard deviation of o>.
Exploratory analysis indicated that a quadratic function provided a reasonable
representation of the relationship between the expected relative biovolume of cyanobacteria,
pc, and Chi a as follows:
E[logit(pci)] = A + f2chli + f2chl2t (20)
wherefi,f2, and f3 are coefficients estimated from the data.
Because laboratory replicates of P, and pC;, were available, uncertainty associated with
measuring phytoplankton and relative biovolume of cyanobacteria was estimated as follows:
log [Brrij ) ~Normal(\og{Bi{j{), sx)	(21)
logit(jpmc j)~N ormal(logit(p c j), s2)	(22)
where Brrij and pmcjare the laboratory measurements of phytoplankton biovolume and
proportion cyanobacteria, respectively, and the index j maps replicate measurements to the
corresponding estimate of the true value of the measurement for sample /'. These laboratory
replicates are assumed to be normally distributed about their respective estimates of the
transformed sample means, with standard deviations of Si and ss, respectively.
Cyanobacterial biovolume (C) can then be expressed as the product of the relative
biovolume of cyanobacteria and total phytoplankton biovolume. After log-transforming, the
expression is as follows:
log(Cj) = log(fcCji) + log(pc [) + log(CWi)	(23)
where cyanobacterial biovolume in sample i is the sum of a log-transformed parameter kc, the
log-transformed cyanobacterial relative biovolume in the sample, and the log-transformed Chi o
concentration.
The final component of the model relates cyanobacteria biovolume to MC. Initial
exploration of the data indicated that MC increases at a rapid rate relative to cyanobacterial
biovolume at high levels of cyanobacteria. At low levels of cyanobacteria, however, microcystin
increases at a somewhat lower rate. To account for this change in rate, microcystin was
modeled with a piecewise linear model as follows:
41

-------
log (una) = gQg (Q))	(24)
where the response variable in this relationship is Hmcj, the estimated mean concentration of
microcystin in sample /'. The function g(.) is the piecewise linear function, which is characterized
by four parameters: the intercept, di, and slope, d2, of the first segment; the point along the
gradient at which the slope changes, cp; and the slope of the second segment, d3.
The distribution of observed MCs about the mean value was then modeled as a negative
binomial distribution as follows:
MCi~NB(nMC i, (p)	(25)
where MC, is the MC observed in sample i and NB(.) is a negative binomial distribution with
overdispersion parameter, cp. Because the negative binomial distribution specifies only
nonnegative integer outcomes, before fitting the model, the EPA multiplied microcystin
measurements by 10 and rounded to the nearest integer. Microcystin measurements below the
detection limit of 0.1 ng/L were set to zero (Yuan and Pollard 2017).
3.2.3.2 Results
A total of 2,352 observations of MC, cyanobacterial and phytoplankton biovolume, and
Chi a were available from the NLA data set for analysis. Those measurements were collected
from 1,116 different lakes spanning the conterminous U.S. An additional 112 samples of
laboratory replicates of phytoplankton and cyanobacterial biovolume measurements were
available to quantify measurement variability.
Three different relationships were estimated in the national model: (1) Chi a and
phytoplankton biovolume, (2) Chi a and cyanobacterial relative biovolume, and (3)
cyanobacterial biovolume and MC. (The relationship between phytoplankton biovolume,
cyanobacterial relative biovolume, and cyanobacterial biovolume required no statistical
estimation.) The observed relationship between Chi a and phytoplankton biovolume was
accurately represented as a line with a slope equal to 1 on log-log axes (left panel, Figure 21),
similar to the relationship estimated in the zooplankton model.
Cyanobacterial relative biovolume exhibited an increasing relationship with Chi a
(middle panel, Figure 21). The quadratic functional form allowed the model to represent the
steepening of the relationship at higher concentrations of Chi a. Mean MC increased with
cyanobacterial biovolume (right panel, Figure 21). The slope of the relationship increased at a
42

-------
cyanobacterial biovolume of 1.9 mm3/L, but the 90% credible intervals on the location of this
changepoint ranged from 0.5 to 5 mm3/L. At cyanobacterial biovolumes greater than the
changepoint, the slope of the mean relationship was statistically indistinguishable from 1,
whereas at cyanobacterial biovolumes less than the changepoint, the slope was 0.61, with 90%
credible intervals ranging from 0.51 to 0.69. Overall, the credible intervals about the
cyanobacteria-MC relationship were narrow compared to those estimated for the Chi o-
cyanobacterial relative biovolume relationship as shown.
_j
o>
O
2
o
o
o
0.001
0.01
0.1
1
10
o
o
o
"e
E
TJ
C
3
-Q
05
O
l
o
o
>
o
S2
5
o
c
O
o
c
o
I
o
d
o
d
o
o
0.1
1
10
100
1000
1
10
100
Chi a (ng/L)	Chi a (ng/L)	Cyano biovolume (mm3/L)
Figure 21. Modeled relationships for the microcystin model. Left panel: relationship between Chi a and
phytoplankton biovolume; open circles: observed measurements of Chi a and phytoplankton biovolume;
solid line: has a slope of 1. Middle panel: relationship between Chi a and cyanobacterial relative
biovolume; open circles: average cyanobacterial relative biovolume in ~20 samples at the indicated Chi a
concentration; solid line: estimated mean relationship; gray shading: 90% credible intervals about the
mean relationship; vertical axis: has been logit-transformed. Right panel: relationship between
cyanobacterial biovolume and MC; open circles: average MC in ~20 samples at the indicated
cyanobacterial biovolume; solid line: estimated mean relationship; gray shading: 90% credible intervals
about the mean relationship; small filled circles: Chi a bins in which MC in all samples was zero.
3.2.3.3 Chi a criteria derivation
Draft Chi a criteria to protect recreational uses and drinking water sources can be
derived from the estimated network of relationships by combining the model equations for total
phytoplankton biomass, cyanobacterial-relative biovolume, and microcystin and the uncertainty
inherent in each of the relationships (Figure 22). More specifically, based on a threshold
concentration for microcystin and an allowable exceedance frequency of that threshold,
Equation (25) can be used to compute the mean predicted MC that would be associated with
these values. Then, Equations (23) and (24) can be used to calculate the Chi a concentration
associated with this mean MC. This model is based on instantaneous measurements of Chi o,
cyanobacterial biovolume, and MC. To relate instantaneous Chi a concentrations to a seasonal
mean Chi a concentration, the EPA computed the variance of Chi a concentrations within lakes
43

-------
over the summer sampling season using repeat visits included in the NLA data set. Then, the
variance was used to estimate the probability of exceeding an instantaneous Chi a
concentration, based on the seasonal mean Chi a concentration.
Threshold concentrations for microcystin have been published, and those targeted
conditions can guide the use of the models to derive Chi a criteria. To protect sources of
drinking water, the EPA Health Advisory recommends a threshold concentration for microcystin
of 0.3 ju,g/L for preschool children less than 6 years old (US EPA 2015b). This threshold to protect
human health applies to finished drinking water; however, the EPA is aware that states or
authorized tribes apply water quality standards for protecting drinking water sources to either
the ambient source water before treatment or to the finished drinking water after treatment.
The ability of treatment technologies to remove microcystin is too variable (Westrick et al. 2010,
US EPA 2015c) for the EPA to set a national recommendation for a protective ambient source
water concentration that would yield a protective concentration after treatment. If a state or
authorized tribe applies the health advisory standard to finished drinking water, then they can
account for the expected treatment in their facilities and select a higher microcystin
concentration in the ambient source water that would result in the targeted microcystin
concentration in the finished drinking water. This will result in a concentration of Chi a in the
ambient source water that will protect human health from the effects of microcystin in the
finished drinking water. To protect recreational uses, the EPA recommends a threshold
concentration for microcystin of 8 ju,g/L to protect children (US EPA 2019).
44

-------
o
o
o
o
o
Oo
o
o
o
o
c o
C Q)
O "D O
9- A-
OO
Q_ c CM
CD
0.1	1	10	100	1000
Chi a (ng/L)
Figure 22. Example of derivation of Chi a criterion to protect recreational uses based on targeted MC of 8
Hg/L and exceedance probability of 1%. Top panel-open circles: observed values of microcystin and Chi a
for samples in which MC was greater than the detection limit; solid line: predicted MC that will be
exceeded 1% of the time for the indicated Chi a concentration; gray shading: 50% credible intervals about
mean relationship; solid vertical and horizontal line segments: candidate Chi a criterion based on targeted
MC. Bottom panel: proportion of samples for which microcystin was not detected in ~100 samples
centered at the indicated Chi a concentration.
After selecting the designated use of interest, calculating the corresponding Chi a
criterion requires two additional management decisions: selection of the allowable exceedance
probability of the threshold and selection of a credible interval of the model output. These
decisions are combined with a posterior simulation using the estimated distributions of the
model parameters to estimate Chi a criteria. The allowable exceedance probability can be
interpreted directly in terms of environmental outcomes as the probability of observing a
specified MC in a sample for a given seasonal mean Chi a concentration. For example, after
accounting for model uncertainty by selecting the 25th credible interval, MC in lakes with a
seasonal mean Chi a concentration of 22 ng/L would be expected to exceed a threshold of 8
Hg/L in 1% of samples (Table 4) (solid vertical line in Figure 22). The credible intervals express
the uncertainty in the model predictions of different exceedance probabilities. So, the shaded
area in Figure 22 shows the range over which at least 50% of the possible curves would be
45

-------
located that describe MCs that have a 1% probability of exceedance. Selection of lower credible
intervals yields more conservative criteria in terms of model uncertainty. An interactive tool
allowing the user to examine Chi a criteria associated with different combinations of microcystin
threshold, probability of exceedance, and the credible interval is available at https://chl-
microcystin-prod.app.cloud.gov.
Table 4. Illustrative Chi a criteria (ng/L) for different exceedance probabilities using the 25th credible
interval
Probability of
Microcystin threshold = 8 |ig/L
exceedance
to protect recreational uses
1%
22
5%
29
10%
35
3.2.4 Phosphorus-Chlorophyll a
ATP measurement is comprised of P contained within different compartments,
including P bound in phytoplankton, P bound to suspended sediment, and dissolved P (i.e.,
chemically dissolved P and P bound to particles small enough to pass through a filter) (Effler and
O'Donnell 2010). In many lakes, much of measured TP is associated with phytoplankton, and so,
differences in phytoplankton biomass among lakes can be associated with differences in both
Chi a and TP, yielding a strong correlation between the two (Lewis and Wurtsbaugh 2008). In
other lakes, high concentrations of suspended sediment can contribute to TP and affect
observed TP-Chl a relationships (Jones and Knowlton 2005). When TP-Chl a relationships are
being estimated, lakes with high concentrations of suspended sediment show low Chl:TP ratios
relative to the average pattern (Hoyer and Jones 1983, Jones and Knowlton 2005).
The EPA modeled the relationship between TP and Chi a by explicitly accounting for the
contributions of different compartments to observed TP, resulting in the positions of TP and Chi
a being reversed from the typical model formulations: The model explained variations in TP in
various compartments, rather than explaining variation in Chi a (Yuan and Jones 2020b).
46

-------
3.2.4.1 Statistical analysis
The EPA specified a model that estimates contributions to TP from different
compartments, where TP is modeled as the sum of contributions from dissolved P, P bound to
nonphytoplankton sediment, and P bound in phytoplankton (Figure 23).
p
1 diss

Chl
(26), ?
<27y
Turb
(31), (32)
Figure 23. Schematic representation of compartment model for TP. PdisS: dissolved P; Chi: Chlorophyll a;
Turb: total turbidity; Turbnp: turbidity attributed to nonphytoplankton sources. Shaded box for Turbnp: a
variable inferred by the model; numbers in parentheses: refer to equation numbers in the text. Equations
(28)(30) and equations (33)(35) describe the distributions of turbidity and TP measurements and are
not shown in the schematic.
Direct measurements of nonphytoplankton sediment were not collected during the NLA.
Instead, turbidity measurements were available that are associated with total suspended solids
and include contributions from both nonphytoplankton and phytoplankton components.
Because an estimate of nonphytoplankton sediment is needed to model TP, turbidity is modeled
as the sum of two components: (1) turbidity that is directly associated with phytoplankton
biomass, or autochthonous suspended sediment (Turbaut) and (2) turbidity associated with all
other sources, or nonphytoplankton turbidity (Turbnp). The second component of turbidity
includes turbidity associated with allochthonous sediment and sediment resuspended from the
lake basin (Hamilton and Mitchell 1996). The EPA modeled Turbaut as being directly proportional
to Chi a (Jones et al. 2008), a measure of algal biomass and, therefore, the components of
turbidity were expressed as follows:
E[Turb] = Turbnp + Turbaut = Turbnp + fChl	(26)
where E[Turb] indicates that the model applies to the expected value of turbidity (Turb). The
amount of Turbaut associated with each unit of Chi a is expected to vary with algal composition.
For example, small phytoplankton species would tend to scatter light differently than larger
species. Assuming that algal composition changes with trophic conditions (Godfrey 1982), the
47

-------
change in algal composition can be modeled by expressing the coefficient/as an unknown
function of Chi a. Also, assuming that f(Chl) can be modeled as a power function (f= bChlm), the
product of f(Chl) and Chi a can be written as follows as bChlk without any loss of generality:
E[Turb\ = Turbnp + f(Chl)Chl = Turbnp + bChlk	(27)
where the exponent, k, is equal to m+1.
Exploratory analysis indicated that concentrations of Turbnp varied with different lake
characteristics, but the predictor that accounted for the most variability was lake depth.
Therefore, 30 classes of lakes based on maximum depths were defined, and the value of Turbnp
within each of the classes was modeled as a log-normal distribution about a mean value specific
to that depth class as follows:
\og(Turbnp)~Normal(na>i, oa)	(28)
where jua /is the mean value of log(Turbnp) for depth class /', and oa is the standard deviation of
the distribution of individual measurements of Turbnp. The set of values for ^a , was then
assumed to be drawn from a single normal distribution as follows:
Ha>i~Normal(n, ctJ	(29)
where ^ and a^ are the mean and standard deviation of this distribution. The mean distribution
loosely constrains the possible values of iua,i, while allowing lakes with smaller amounts of data
to "borrow information" from lakes with larger amounts of data (Gelman and Hill 2007).
Finally, sampling variability for Turb was assumed to be log-normally distributed as
follows:
log (Turb) ~Normal(E[log(Turb)], aT)	(30)
where E[log(Turb)] is the expected value of log(Turb) expressed in Equation (27).
The EPA used results from the model for turbidity simultaneously to estimate
contributions to different components of TP. Recall that TP is modeled as being composed of
contributions from dissolved P (Pdiss), P bound to suspended sediment, and P bound to
phytoplankton. Based on this assumption, the following model equation can be written:
E[TP] = Pdiss + g^urb^ + g2Chl	(31)
48

-------
where the concentration of P bound to nonphytoplankton suspended sediment is modeled as
being directly proportional to Turbnp, and P bound to phytoplankton is modeled as being directly
proportional to Chi a. The coefficient gi quantifies the P content of Turbnp, while the coefficent
g2 expresses P concentration relative to Chi a concentration. P content is expected to vary with
the level of turbidity and the composition of the phytoplankton assemblage, so, similar to the
model for turbidity, the coefficients gi and g2 were allowed to vary as power functions of Turbnp
and Chi a, respectively. So, the final model equation can be written as follows:
E[TP] = Pdiss + djurb + d2Chln	(32)
Exploratory analysis indicated that dissolved P was associated with lake depth, so,
similar to Turbnp, different values of Pdiss were estimated for each of 30 lake depth classes as
follows:
log(P^[SS mn,i) ~ Nor7no,l(^n(iiss, 
-------
3.2.4.2 Results
Observations of turbidity were correlated with Chi a, and a distinct lower boundary in
the scatter of data was evident (Figure 24). The model relationship defining this lower boundary
can be computed by setting Turbnp to zero in Equation (27). Then, after log-transforming, the
equation can be written as log(Turb) = log(b) + klog(Chl). In other words, when Turbnp is
negligibly small, the relationship between Turbaut and Chi a is a straight line in the plot of
log(Chl) vs. log(Turb) (solid line in Figure 24). Deviations in sampled values above that line show
the contribution of Turbnptothe overall turbidity measurement. Mean values of fa and k
estimated from the model were 0.67 (0.62, 0.73) and 0.67 (0.65, 0.69) (90% credible intervals
shown in parentheses). Based on the functional form that was assumed for the relationship
between turbidity and Chi a, the contribution of phytoplankton to turbidity (i.e., Turbaut/Chl a)
was estimated as being proportional to Chi"033. That is, as Chi a increases, the amount of
turbidity associated with each unit of Chi a decreases, a trend that is consistent with a shift from
small-bodied, diatom-dominated assemblages to colonies of cyanobacteria cells (Scheffer et al.
1997).
o
I I 111	1	1 I I I I ll|	1I rTTTTTJ	1I I I I 1111	1I I I I 1111
0.1	1	10 100 1000
Chi a (ng/L)
Figure 24. Turbidity vs. Chi a. Solid line: the limiting relationship between Chi a and turbidity when
contribution of allochthonous sediment is negligible.
Estimates of Turbnp and mean dissolved P both exhibited decreasing relationships with
increasing depth (Figure 25). Turbnp decreased from approximately 1.4 nephelometric turbidity
units (NTU) in shallow lakes to nearly zero in deep lakes, while Pdiss varied from approximately
2.6 ng/L in shallow lakes to 1.6 ng/L in deep lakes. Both of these relationships are consistent with
a mechanism by which fine sediment from the lake bottom is likely to be collected in surface
50

-------
water samples in shallow lakes. In the case of Pdiss, the trend indicates that measurements of
dissolved and particulate components of TP are determined by filter size and P bound to
sediment fine enough to pass through the filter contributes to estimates of dissolved P.
CO
o
o
o
oou
TT-rr
1
-1	1	1IIn-|
10
Depth (m)
oooo o o o
I	1	1	1
CD
c\i
c\i
CM
c\i
CD
c\i
o q,
o o
CP
O o
O o
o o
t-^t
1
-1	1	1II I I I
10
Depth (m)
Figure 25. Relationship between Turbnp, Pdiss, and lake depth. Open circles: mean estimate of parameter
value in each of 30 lake depth classes.
The quantity of P bound to nonphytoplankton suspended sediment expressed by the
coefficient di exhibited substantial geographic variation (Figure 26). Coherent spatial patterns
could be discerned in the variation of di among different states, with relatively high levels of P
content in the upper midwest region of the country (e.g., Montana, North Dakota, and South
Dakota) as well as in parts of the western mountains. Comparatively lower levels of P content
were observed in the northeast region of the U.S. Mechanisms for these large-scale variations in
P content are likely related to the underlying geology of soils in each region (Olson and Hawkins
2013). Values of d2, the amount of P within phytoplankton, spanned a much narrower range
than estimated for di, only ranging from 1.6 to 4.5 per unit of Chi a. The relative difference in
regional variability in the coefficients indicates that spatial differences in the amount of P bound
to nonphytoplankton suspended sediment account for more of the variability in TP-Chl a
relationships than spatial differences in P within phytoplankton, and the amount of P residing in
phytoplankton is relatively constant.
51

-------
log(d1)
2	3	4	5
Figure 26. Ecoregion-specific values of loge(c/i), P bound to nonphytoplankton suspended sediment.
Limiting relationships that estimate the P content of phytoplankton biomass and Turbnp
can also be calculated (Figure 27). For phytoplankton biomass, the limiting relationship is
calculated by setting PdiSS and Turbnp in Equation (32) to zero, yielding the following log-
transformed relationship: log(TP) = log(d2) + n log (Chi). Different values of d2 were
estimated for each ecoregion, but the distribution of those values is characterized by an overall
mean value of 2.5 (2.0, 3.1), while the mean value of the parameter n was 0.87 (0.82, 0.92). The
straight line based on the two parameter values represents P associated with phytoplankton
biomass, as quantified by Chi o, and it tracks the lower limit of the observed data (solid line,
right panel, Figure 27). As a limiting relationship, one would expect that the majority of values of
TP would be greater than this line indicates, but variability associated with the value of d2 causes
some values of TP to fall below the limit.
For Turbnp, setting Pdiss and Chi a to zero yields the following relationship: log(TP) =
log(d1) + mlog (Turbnp). The mean value of the coefficient di was 31 (23, 40), and the value of
the exponent m was 0.35 (0.32, 0.40) (left panel, Figure 27). Overall, the RMS error for
predicting log(TP) was 0.48 for the model.
52

-------
0.001 0.1 1 10 100
Turbnp (NTU)
0.1	1	10 100 1000
Chi a (ng/L)
Figure 27, TP versus Turbnp and Chi a. Solid lines: the limiting relationship between the indicated variable
and TP; gray shaded areas: the 90% credible intervals about the mean relationship.
3.2.4.3 Phosphorus criteria
Two relationships between Chi a and TP that can be inferred from the TP model inform
the derivation of draft TP criteria. First, the limiting relationship between Chi a and TP estimated
from the model quantifies the amount of P that is bound to phytoplankton (Figure 27). This
relationship predicts TP concentration in samples in which suspended sediment and dissolved P
concentrations are very low and defines the minimum value of TP that is associated with a
targeted Chi a concentration. This limiting relationship can also be interpreted as the Chi a yield
of P (Yuan and Jones 2019) and could be used to predict the change in Chi a that would
potentially result from a change in the amount of biologically available P in the water column
(Reynolds and Maberly 2002).
A second relationship between TP and Chi a accounts for contributions from P bound to
nonphytoplankton sediment. If lake depth is specified, then the relationship estimated between
lake depth and nonphytoplankton sediment can be used to estimate an average contribution to
TP from these other compartments in the water column (Figure 25). The resulting relationship
then provides an estimate of the ambient TP concentration one would expect to observe as a
function of Chi a.
53

-------
1	10
Chi a (ng/L)
1000
Figure 28. Example of deriving TP criteria for a Chi a target of 10 ng/L for data from one ecoregion
(Southeastern Plains). Open circles: all data; filled circles: data from the ecoregion; solid line: limiting TP-
Chl a relationship from compartment model; dashed line: ambient TP-Chl a relationship taking into
account contributions from nonphytoplankton sediment for a 3-m deep lake; solid horizontal and vertical
line segments: Chi a target and associated TP criteria.
Table 5. Illustrative example of TP criteria corresponding to data shown in Figure 28. Example TP criteria
for illustrative Chi a targets. Ambient TP criteria calculated for a 3-m deep lake.
Chi a = 10 |ig/L	Chi a = 15 |ig/L
10th credible	25th credible	10th credible	25th credible
interval	interval	interval	interval
15	16	22	23
Limiting relationship
(TP Hg/L)
Ambient (TP ng/L)	23	25	30	32
Information from the two Chi a and TP relationships specifies a range of possible TP
criteria that can be associated with a desired concentration of Chi a (Figure 28). The prediction
of ambient TP that accounts for contributions from nonphytoplankton sediment provides an
estimate of the mean TP concentration that one would expect to observe for a given Chi a. As
such, this ambient TP concentration provides a candidate criterion. Note that contributions of
Pdiss are not included in predictions of ambient TP criteria. In many lakes Pdiss is composed of
more biologically available forms of P (e.g., soluble reactive P), and so, concentrations of Pdiss
should be near zero in lakes in which reductions in P loading would be expected to influence
phytoplankton abundance.
54

-------
The lower limiting relationship identifies the minimum possible TP concentration one
might expect to observe for a given Chi a. This limiting relationship between TP and Chi a can
also potentially be used to predict changes in Chi a from a change in loads of biologically
available P (Reynolds and Maberly 2002), information that can guide the development of waste
load allocation. Final uses of the range of values provided by these models depend on the
specific applications in each state and on the risk management decisions made by the state.
The interactive tool for computing different TP criteria associated with Chi a is available
at https://tp-tn-chl-prod.app.cloud.gov. This tool allows the user to specify the targeted Chi a
concentration and the lake depth of interest. Because the coefficients di and d2 vary among
ecoregions (Figure 26), users also can select a particular ecoregion for computing TP criteria.
Finally, users can select the confidence level, expressed as a lower credible interval, for
examining the effects of model uncertainty on the calculated criteria. Data selected for an
ecoregion are highlighted in the provided plots. The model then computes TP associated with
those conditions using a posterior simulation from the Bayesian model results. A lower credible
interval provides additional assurance that the calculated criterion is protective, given the data
and model uncertainty. For example, selecting the 25th credible interval implies that only 25%
of predicted TP concentrations at the selected Chi a concentration, given the data, were less
than candidate criterion value criteria. In statistical hypothesis testing, convention suggests that
p-values of 1% or 5% are statistically significant results. Those practices can also inform the
selection of the credible interval, but selection of the value of the lower credible interval as the
basis for the criteria is ultimately a water quality management decision.
3.2.5 Nitrogen-Chlorophyll a
Similar to the model for TP, each TN measurement is comprised of N contained within
three compartments: N bound in phytoplankton, dissolved inorganic N (i.e., nitrate, nitrite, and
ammonia), and dissolved organic N (DON). Unlike the TP model, exploratory analysis indicated
that the N content of inorganic suspended sediment was negligible (Yuan and Jones 2019).
3.2.5.1 Statistical analysis
Field measurements of the difference between TN and dissolved inorganic nitrogen (DIN
= NOx + ammonia) were modeled as follows:
E[TN - DIN] = ftChlkl + DON = fcCM"1 + f2D0C	(36)
55

-------
where variations in the observations of total N minus dissolved inorganic N (TN-DIN) are
attributed to two compartments: N bound in phytoplankton, modeled as fiChlkl and DON.
Exploratory analysis indicated that DON was closely associated with DOC, as they often originate
from the same watershed sources (Berman and Bronk 2003), so the concentration of DON was
modeled as being directly proportional to DOC.
As with the TP model, exploratory analysis indicated that the parameters fi and f2 varied
most strongly with geographic location. Because of those trends and to facilitate the use of this
model with local data sets, different values of fi, and f2 were specified for each Level III
ecoregion:
lg(/i,i) ~ Normal(nfl, afl)
lg(/2,i) ~ Normal(nf2, af2)	(37)
where the parameters Hfi and pif2 estimate the mean values of the distribution of/i and/2 while
Ofi and Of2 estimate the standard deviations.
The sampling distribution of TN-DIN was assumed to be log-normally distributed as
follows:
log(7W - DIN) ~Normal(E[log (TN - DIN)], aTN)	(38)
where oTn is the standard deviation of observed values of log(TN-DIN) about their expected
value.
3.2.5.2 Results
A total of 2466 samples collected from 1875 lakes were available for analysis. Values for
the coefficient, fi, quantifying phytoplankton N content ranged from 11 to 43 in different
ecoregions with an overall mean value of 18.3 (14.9, 22.3). The values estimated for/2 spanned
a greater range among ecoregions with a minimum value of 35 and a maximum value of 103.
The overall mean value of/2 was 64.9 (61.0, 68.9). The broad range in values of f2 indicates that
strong differences exist among different locations regarding the nature of the relationships
between DOC and DON. The mean value of the exponent, kl, was 0.90 (0.86, 0.94).
To visualize the variability in phytoplankton N among ecoregions, the concentration of N
bound in phytoplankton at the overall mean Chi a concentration of 9.3 ng/L is mapped (Figure
29). With the exception of one high value of 320 ng/L estimated for the Sand Hills, Nebraska
56

-------
ecoregion, N-content of phytoplankton exhibited only small variations among ecoregions. N
content ranged from 83 - 185 ng/L with a median value of 136 ng/L. Coherent spatial patterns
in the N-content of phytoplankton were not evident.
Phytoplankton N
T	I	I
100	150	200	250	300
Figure 29. Variation in the concentration of N bound in phytoplankton among Level III ecoregions at the
overall mean Chi a = 9.3 ng/L. Gray scale shows N concentrations in ng/L.
Estimated DON concentrations at the overall mean DOC concentration of 5.6 mg/L
ranged from 194 - 570 ng/L with a median concentration of 365 ng/L (Figure 30). Variations in
DON among ecoregions were substantially greater than observed for phytoplankton N. Spatial
patterns were also evident, with higher concentrations of DON in the upper Midwest regions of
the United States and lower concentrations in the mountains in the western and eastern regions
of the country.
I
57

-------
DON
 , ,
200	300	400	500
Figure 30. Variation in DON Level III ecoregions at an overall mean DOC = 5.6 mg/L. Gray scale shows N
concentrations in ng/L.
The EPA calculated limiting relationships that estimate the N content of phytoplankton
biomass with a procedure identical to that used for TP (Figure 31). In this case, the limiting
relationship was calculated by setting the contribution from DON in Equation (36) to zero,
yielding the following log-transformed relationship: log (TN  DIN) = log^) + k log (Chi).
The straight line based on those two parameter values represents N associated with
phytoplankton biomass, as quantified by Chi a, and it tracks the lower limit of the observed data
(solid line, left panel Figure 31).
Similarly, setting DIN and Chi a to zero in Equation (36) yields the following limiting
relationship for DON: log(7W) = log(/2) + log (DOC) (solid line, right panel Figure 31). The
mean value of f2 indicates that, on average, the concentration of DON was 0.065 times that of
DOC. Overall, the RMS prediction error for log(TN-DIN) was 0.37.
58

-------
o O
o
o
OcO
I
-	o
U)	o
=L	O
z
I-
oor
o
1
10
100
oo
o
o
I
^	o
O)	o
=L	O
z
H
oo
o
CO
Chi a (jj-9/L)	DOC (mg/L)
Figure 31. TN-DIN vs. Chi a and DOC. Solid lines: the limiting relationship between each variable and TN-
DIN; shaded area: the 95% credible intervals about this mean relationship.
3.2.5.3 Nitrogen criteria
As with TP, the model for TN-DIN provides two different predictions of TN-DIN
concentration, given the Chi a concentration. The prediction for the ambient concentration of
TN-DIN accounts for the increase in TN-DIN one would expect with increased Chi a, but also
includes contributions from DON (as estimated by DOC) and OSnp in the lake. Mean predictions
for TN-DIN can be computed for different values of Chi a that include average contributions
from other sources of N in the water column. The value of this ambient TN-DIN concentration
that is associated with a targeted Chi a concentration then provides a candidate criterion for TN-
DIN. The second prediction of TN-DIN can be estimated from the limiting relationship between
Chi a and TN-DIN (Figure 31). This relationship quantifies the amount of N that is bound in
phytoplankton, a quantity that is also referred to as the "Chi a yield of nitrogen" (Gowen et al.
1992). This limiting relationship can potentially be used to estimate the change in Chi a that
would result from a change in the amount of biologically available N in the water column
(Reynolds and Maberly 2002).
Criteria for N concentrations are commonly expressed in terms of TN rather than TN-
DIN. To convert a candidate criterion for TN-DIN to a criterion for TN, the availability of DIN for
phytoplankton uptake can be considered. More specifically, the components of DIN (NOx and
ammonia) are easily assimilated by phytoplankton and, when excess concentrations of DIN are
observed in a lake, it may indicate that factors other than N availability are limiting
phytoplankton growth. Therefore, controlling phytoplankton growth by reducing available N
59

-------
would first require that DIN concentrations are reduced to near zero and, when that occurs,
criteria expressed for TN-DIN would be the same as those for TN.
0.1	1	10	100	1000
Chla(ng/L)
Figure 32. Illustrative example of deriving TN criteria for a Chi a target of 10 ng/L for one ecoregion
(Southeastern Plains). Open circles: all data; filled circles: data from selected ecoregion; solid line: limiting
TN-DIN vs. Chi a relationship from compartment model; dashed line: mean ambient TN-DIN vs. Chi a
relationship taking into account mean DOC observed within the selected ecoregion: shaded area: 80%
credible intervals about mean relationships; horizontal and vertical solid line segments: Illustrative Chi a
target and associated TN criteria.
Table 6. Illustrative example of TN criteria corresponding to data shown in Figure 32.
Chi a = 10 |ig/L	Chi a = 15 |ig/L
10th credible	25th credible	10th credible	25th credible
interval	interval	interval	interval
Limiting relationship	
* . K	110	120	160	170
(TN ng/L)
Ambient (TN ng/L)	380	390	440	450
The same interactive tool for computing different TP criteria also provides TN criteria
associated with Chi a (https://tp-tn-chl-prod.app.cloud.gov). This tool allows the user to specify
the targeted Chi a concentration, DOC concentration, and an ecoregion of interest. Finally, users
can select the confidence level, expressed as a lower credible interval, for examining the effects
of model uncertainty on the calculated criteria. Data selected for an ecoregion are highlighted in
the provided plots.
60

-------
4 Characterization
4.1 Other Measures of Effect and Exposure
A variety of other measures of effect and exposure could be used for deriving nutrient
criteria associated with each of the pathways described in Figure 1 and Figure 2. In selecting the
responses for analysis, the EPA considered (1) available data, (2) the current state of scientific
understanding of each pathway, and (3) the degree to which a pathway and a response could be
applied broadly to most lakes. For many possible measures of effect and exposure, data
availability was a key consideration. For aquatic life, direct measurements of fish assemblage
composition and biomass were not collected during the NLA, and the lack of those data limited
the potential for considering several pathways such as evaluating alterations in fish assemblage
composition because of reduced visibility. Lake benthic communities also exhibit changes along
a eutrophication gradient (Vadeboncoeur et al. 2003), but none of those data were available.
For recreational and drinking water source uses, the effects of other cyanotoxins (e.g.,
cylindrospermopsin) might be important for certain lakes, but continental-scale data for those
other cyanotoxins were not available at the time of this analysis. In certain lakes, cyanobacterial
blooms have also been observed at depths below the surface layer (Jacquet et al. 2005), but
observations of phytoplankton at those depths were not available. Similarly, organic matter
generated by increased primary productivity can increase the concentrations of disinfection by-
products during the drinking water treatment process (Graham et al. 1998, Galapate et al.
2001), and chemicals produced during blooms of certain algal species can introduce unpleasant
taste and odors to drinking water (Graham et al. 2010). However, continental-scale data
pertaining to disinfection by-product precursors or taste and odor chemicals were not available.
Insufficient scientific understanding of a causal pathway also limited consideration of
certain measures of effect and exposure. For example, scientific consensus is currently lacking
on the precise level of cyanobacteria that is harmful to aquatic life. That information gap limited
the utility of using cyanobacterial abundance as a final response measurement, despite the fact
that increased cyanobacterial abundance occurs frequently with nutrient pollution (Dolman et
al. 2012). (Note, however, that cyanobacterial abundance measurements quantify a key step in
the model linking Chi a to microcystin.) Similarly, increased levels of cyanobacteria can cause
rashes on people who contact the water (Pilotto et al. 1997, Zhang et al. 2015, US EPA 2015b),
61

-------
potentially affecting the use of a lake for recreation. However, precise quantitative relationships
between the occurrence of rashes and cyanobacterial abundance are not currently available.
For certain measures of effect or exposure, data were available, but other factors
limited the degree to which the response could be applied. For example, Secchi depth data were
available in the NLA data set, and that measure of transparency could have informed an
assessment of the aesthetic appeal of different lakes for recreation. That is, increased nutrient
concentrations cause increases in the abundance of phytoplankton that reduce water clarity and
decrease the aesthetic appeal of a lake (Carvalho et al. 2011, Keeler et al. 2015). Aesthetic
considerations have been used by others as a basis of water quality criteria (Heiskary and Wilson
2008), but the aesthetic expectations for different lakes depend on geographic location
(Smeltzer and Heiskary 1990), and user perception survey data at the continental scale of this
analysis were not available. Similarly, reducing the frequency of phytoplankton blooms has been
cited as a motivation for controlling nutrient loads (Bachmann et al. 2003), but aesthetic
expectations regarding bloom frequency were not available at the national scale.
4.2 Incorporating State Data
State water quality managers are often interested in exploring relationships between
environmental factors and biological responses using locally collected monitoring data. In many
cases, leveraging knowledge from broader regional scales (e.g., national scale) can enhance local
understanding. This document describes draft recommended numeric nutrient criteria models
based on national data that link designated uses to Chi a, TN, and TP. The NLA data set provided
a comprehensive set of measurements collected from large numbers of sites with identical
protocols (US EPA 2011, Pollard et al. 2018), and the availability of consistent measures from
lakes spanning broad gradients facilitated the calculation of accurate national estimates of
relationships of interest. However, the number of samples is limited within the national data set
that is available to estimate relationships within any single state, and uncertainty in estimating
relationships specific to a single state is higher than that associated with the national models. In
contrast, monitoring conducted by state agencies can yield more intensive temporal sampling
over more sites, and hence, relationships estimated from those data can assist local
management decisions within that state. Data collected at the state level, however, can be
limited in the parameters that are measured, and the range of environmental conditions
sampled is limited by conditions occurring within the state boundaries.
62

-------
All the draft recommended criteria models described in this document are formulated
to facilitate consideration of state data. State-specific values for certain coefficients in each
model (e.g., Figure 33) have been estimated, and local state, monitoring data can be used to
refine the estimates of state-specific coefficients, while remaining consistent with national
trends. Appendices A, B, and C discuss three examples of case studies in which state monitoring
data have been combined with national data to refine draft recommended criteria. State
monitoring data sets are each unique, and the EPA is available to assist states in combining their
monitoring data with the national models.
4.3 Existing Nutrient-Chlorophyll a Models
Empirically estimated relationships between TP and Chi a concentrations have provided a
basis for lake water quality management for over four decades. This relationship was initially
identified in Connecticut and Japanese lakes (Deevey 1940, Sakamoto 1966), and subsequently
extended to a broad range of temperate lakes in the mid-1970s (Dillon and Rigler 1974, Jones and
Bachmann 1976, Carlson 1977). Those early analyses regressed Chi a on TP concentrations and
reported similar coefficients showing the ratio of Chl:TP increased with lake trophic state. Over
time, many studies have explored the veracity of that relationship and assessed sources of residual
variation, testing the limits of applicability to different regions and lake types (e.g., McCauley et al.
1989; Prairie et al. 1989; Jones and Knowlton 2005; Filstrup, Wagner, et al. 2014). Variations in the
relationship have been attributed to differences in lake depth (Pridmore et al. 1985), TN:TP ratio
(Smith 1982, Prairie et al. 1989, Molot and Dillon 1991), grazing by zooplankton and mussels
(Mazumder 1994, Mellina et al. 1995), landscape characteristics (Wagner et al. 2011), and light
limitation (Hoyer and Jones 1983, Knowlton and Jones 2000, Havens and Niirnberg 2004). Regional
studies have evaluated the relationship as influenced by edaphic and climatic factors in locations
such as Canada (Prepas and Trew 1983), Argentina, (Quiros 1990), the United Kingdom (Spears et
al. 2013), and Europe (Phillips et al. 2008). Recently, lake classifications have improved the
precision and accuracy of this relationship (Yuan and Pollard 2014).
As described in Sections 3.2.4 and 3.2.5, the EPA reformulated the nutrient-chlorophyll
models to account for variations in TP and TN, rather than in Chi a. The new models better
account for variability in measurements of TP and TN and are consistent with an understanding
of the components of TP and TN in the water column. The reformulated models cannot be
directly compared with earlier studies, including those cited previously. Estimates of N and P
63

-------
content of phytoplankton, however, are consistent with values reported elsewhere (Yuan and
Jones 2019).
4.4 Limitations and Assumptions
The draft recommended models for deriving numeric nutrient criteria are limited by the
nature of the data that underlie the analysis. First, nutrient data for each lake consisted of
samples collected at a single point, resulting in no information on within-lake spatial variability
in nutrient concentrations being included in the analyses. Nutrient concentrations within
particular lakes can vary considerably across different locations (Perkins and Underwood 2000),
resulting in criteria based on samples collected at the deepest point or midpoint of the reservoir
that might not be applicable to samples collected elsewhere. When deriving their criteria, states
may specify assessment methodologies to collect samples from different locations in the same
lake to address this issue and analyze those local data to account for spatial variability.
Similarly, nutrient and response data used in the current analysis were collected only in
the summer, so monitoring data assessed with respect to these draft recommended criteria
should also be limited to summer data. Nutrient concentrations in some lakes can vary
considerably between summer and winter (S0ndergaard et al. 2005), and states may specify
assessment protocols to ensure that only data collected in the summer are compared with
criterion concentrations.
As noted earlier, most of the draft statistical criteria models described here combine the
effects of spatial, temporal, and sampling variability and estimate a single value for each model
that is applicable to all lakes in the data set. The components of variability, however, might
differ across lakes and affect the resulting criteria. For example, spatial variability in complex,
dendritic reservoirs can be much greater than in simple, circular lakes (Gloss et al. 1980). In most
cases, local monitoring data can inform and potentially improve the parameter estimates both
for specific locations and for groups of lakes.
The uncertainty estimated for each modeled relationship is associated with the number
of samples used in the model, and consideration of sample size can affect the interpretation of
the resulting candidate criteria. For example, the number of NLA samples within a single Level III
ecoregion can be small. The hierarchical structure of the model does improve the precision of
model estimates in those ecoregions, but the precision of TP and TN criteria specific to
64

-------
ecoregions with small amounts of data could be further improved by including state monitoring
data. Additional national-scale data such as that from the 2017 NLA may also be incorporated as
they become available to improve model precision.
Draft recommended criteria based on the drinking water health advisories for
microcystin incorporate some conservative assumptions that affect the final values. The draft
recommended criteria are intended to reflect the ambient water quality conditions that protect
a drinking water use before treatment. They do not, however, account for the varying levels of
treatment a drinking water facility can implement to remove microcystin before generating
finished drinking water, the condition of the water to which the cyanotoxin health advisories
apply. As a precautionary step, a drinking water facility may implement treatment protocols that
minimize the breakage of cyanobacteria cells (Chow et al. 1999, Westrick et al. 2010) which, in
turn, would minimize the release of intracellular microcystin into the treated water. The EPA
based the draft recommended models on the total microcystin present in the NLA samples, both
dissolved in the water and within cyanobacterial cells, which necessitated the lysis of
cyanobacterial cells prior to microcystin quantification, a process that some drinking water
treatments for cyanotoxins are designed to limit. Criteria based on the draft national models
provide protective water quality conditions in the source water, but concentrations of
microcystin that slightly exceed health advisory values can be further reduced in the finished
drinking water through carefully engineered and operated source water treatment processes.
Draft recommended criteria derived using the models described here provide
concentrations that, when exceeded, are associated with a loss of support for designated uses,
but the draft models do not provide information regarding appropriate remediation actions.
Indeed, among lakes in which the criteria are exceeded, appropriate remediation actions will
likely differ. In some lakes, the magnitude of N loading from anthropogenic sources is small,
while P loading is large, and cyanobacteria supply N to the system via fixation (Schindler et al.
2008). In those lakes, reductions in P loading might be the appropriate water quality
management action. In other lakes, ample supplies of N from anthropogenic sources are
available, and management actions might need to focus on reducing both N and P loading
(Ferber et al. 2004). In some lakes, excess N in the form of inorganic nitrogen (NOx or ammonia)
is abundant, and the presence of high concentrations of DIN might provide insights into the
effects of different management interventions. For example, DIN is readily taken up by
65

-------
phytoplankton, so the presence of large concentrations of DIN might indicate that other factors
such as light availability limit phytoplankton growth. In those cases, initial reductions of N
loading to reduce NOx might be necessary before the effects of N control can be observed.
4.5	Deriving State-Specified Criteria
Criteria derived from the draft recommended national models vary with differences in
lake characteristics (e.g., depth and ecoregion), and specifying a single set of criteria applicable
to all lakes in a state might not account for those variations. Methods are already available for
deriving criteria that account for natural variations among water bodies that can be applied to
ensure that appropriate criteria are applied to different types of lakes. First, states can classify
water bodies and derive different criteria for each class of water body. The draft recommended
national models facilitate the classification of lakes by providing specific insights into the factors
that most affect the derivation of protective numeric nutrient criteria. Furthermore, the draft
national models can be used to compute candidate criteria for different lakes in a state to
provide information about the types of lakes for which criterion magnitudes are most similar.
For example, different draft recommended criteria for TP and TN are associated with different
Level III ecoregions; however, among the ecoregions within one state, the difference in criterion
magnitudes might be small enough to specify a single set of criteria applicable to multiple
ecoregions. Second, site-specific criteria can be specified for a small number of lakes with
characteristics that differ substantially from the rest of the lakes in a state. Here, too, the draft
recommended national models provide the means of deriving these criteria for individual lakes.
4.6	Duration and Frequency
The duration component of a water quality criterion is the length of time over which
discrete water samples are averaged to assess the condition of the water body. The frequency
component defines the number of times over a given time period that the specified magnitude
of the criterion can be exceeded while the water body is still assessed as being in compliance
with the criterion and maintaining designated uses. In conjunction with the magnitude of the
criterion, these additional components define a water quality criterion.
Specification of duration and frequency components of numeric nutrient criteria is
complicated by the fact that the ecological effects of elevated nutrient concentrations usually
arise from a sequence of events. For example, higher nutrient concentrations increase the
66

-------
abundance of phytoplankton. Over time, higher abundances of phytoplankton then increase the
amount of organic material in the deeper waters of a lake, and decomposition of the stored
organic material can reduce the concentrations of DO. In this case, the duration and frequency
of exceedance of a Chi a concentration in the mixed layer of the lake is related only indirectly to
the ecological effect of decreased DO and the ultimate reductions in the amount of habitat for
cool- and cold-water species. Contrast this example with the specification of duration and
frequency of toxic pollutants, for which the length of time and frequency of exposure to the
pollutant can be directly linked to effects on different organisms (e.g., mortality). A second
consideration arises from the variability of environmental measurements, for which estimates of
mean concentrations of Chi a, TN, and TP can only be estimated from a finite number of
samples. So, when specifying duration and frequency components of the draft recommended
numeric nutrient criteria, the EPA considered both the timescale of the ecological responses and
the statistical uncertainty in estimating mean values.
The draft recommended duration for draft Chi a criteria derived from the models
described in this document is a growing season (typically summer) geometric mean value,
consistent with the summary statistic used for Chi a in the stressor-response analyses. The
geometric mean was selected to account for the fact that Chi a measurements are frequently
log-normally distributed. The EPA used seasonal mean Chi a concentrations integrated over the
photic zone for analysis because timescales of ecological responses to increased nutrient
concentrations are long. For example, as described earlier, some of the increase in deepwater
oxygen demand arises from accrual of organic material over long time periods while other
oxygen demand arises from recently created organic matter that settles through the water
column. Mean Chi a concentration in the lake is associated with mechanisms acting at both
timescales, providing a measure of the average amount of organic material supplied by the
photic zone. Similarly, systematic changes in zooplankton composition would be expected to
occur at longer, seasonal timescales. For the microcystin model, the basic unit of analysis was an
individual sample, in which the model predicted the probability of different MCs in a sample,
given the sample's Chi a concentration. When estimating the relationship for computing criteria,
however, the EPA computed probabilities of different individual Chi a concentrations as a
function of seasonal mean Chi a concentration, again linking seasonal mean Chi a concentration
to the probability of deleterious effects.
67

-------
The unit of analysis for models relating Chi a to TN and TP concentrations was also the
individual sample, in that TN and TP concentrations measured within a water sample were
described as the sum of phytoplankton-bound N and P and other compartments in the sample
containing those nutrients. Because the models are expressed as simple sums of components,
each one remains applicable even if expressed in terms of seasonal averages. Hence, seasonal
geometric mean Chi a criteria can be converted to seasonal geometric mean TN and TP criteria
using the same model, and the draft recommended durations for TN and TP criteria are also
seasonal mean values.
The draft recommended frequency component of Chi a, TN, and TP criteria derived from
the national models is no exceedances, but the EPA recognizes that seasonal geometric mean
concentrations of Chi a, TN, and TP can vary among years about a long-term mean. As described
above, the timescale over which excess nutrients cause impacts to designated uses is long.
Conversely, occasional deviations of seasonal mean nutrient concentrations above a criterion
magnitude are unlikely to cause immediate, deleterious effects to uses. Draft recommended
nutrient criteria derived from these national models are intended to identify TN, TP, and Chi a
concentrations that, if maintained, on average will ensure protection of applicable designated
uses. Seasonal mean nutrient concentrations do vary, however. For example, a year with
particularly high precipitation might yield higher than average loads of TP to downstream lakes.
Similarly, a year with longer than average periods of sunshine might lead to higher rates of
accumulation of phytoplankton biomass and higher concentrations of Chi a. Hence, in lakes in
which long-term mean concentrations of Chi a, TN, and TP are below the criteria, occasional
seasonal mean concentrations might still exceed the criterion magnitude. Interannual variability
in seasonal mean concentrations can be addressed by characterizing this variability and
incorporating it into the expression of criteria or assessment methods. More specifically, states
may identify adjusted criterion magnitudes with allowable frequencies of exceedance based on
the observed interannual variability of nutrient concentrations. For example, seasonal mean
concentrations of TP would be expected to exceed a criterion magnitude that is equal to the
long-term mean in approximately 50% of the years, whereas less frequent exceedances of a
higher criterion magnitude would be expected. Appendix D provides examples of calculations
that identify different combinations of criterion magnitudes and frequencies.
68

-------
5 References
Arend, K. K., D. Beletsky, J. V. DePinto, S. A. Ludsin, J. J. Roberts, D. K. Rucinski, D. Scavia, D. J.
Schwab, and T. O. Hook. 2011. Seasonal and interannual effects of hypoxia on fish
habitat quality in central Lake Erie. Freshwater Biology 56:366-383.
Bachmann, R. W., M. V. Hoyer, and D. E. C. Jr. 2003. Predicting the frequencies of high
chlorophyll levels in Florida lakes from average chlorophyll or nutrient data. Lake and
Reservoir Management 19:229-241.
Baird, O. E., and C. C. Krueger. 2003. Behavioral thermoregulation of brook and rainbow trout:
Comparison of summer habitat use in an Adirondack River, New York. Transactions of
the American Fisheries Society 132:1194-1206.
Barnett, V., and A. O'Hagan. 1997. Setting Environmental Standards: The Statistical Approach to
Handling Uncertainty and Variation. Chapman and Hall/CRC, London, UK.
Bednarska, A., and P. Dawidowicz. 2007. Change in filter-screen morphology and depth
selection: Uncoupled responses of Daphnia to the presence of filamentous
cyanobacteria. Limnology and Oceanography 52:2358-2363.
Benndorf, Jii., W. Boing, J. Koop, and I. Neubauer. 2002. Top-down control of phytoplankton: the
role of time scale, lake depth and trophic state. Freshwater Biology 47:2282-2295.
Berman, T., and D. A. Bronk. 2003. Dissolved organic nitrogen: a dynamic participant in aquatic
ecosystems. Aquatic Microbial Ecology 31:279-305.
Burns, N. M. 1995. Using hypolimnetic dissolved oxygen depletion rates for monitoring lakes.
New Zealand Journal of Marine and Freshwater Research 29:1-11.
Cahill, K. L., J. M. Gunn, and M. N. Futter. 2005. Modelling ice cover, timing of spring
stratification, and end-of-season mixing depth in small Precambrian Shield lakes.
Canadian Journal of Fisheries and Aquatic Sciences 62:2134-2142.
69

-------
Carey, C. C., B. W. Ibelings, E. P. Hoffmann, D. P. Hamilton, and J. D. Brookes. 2012. Eco-
physiological adaptations that favour freshwater cyanobacteria in a changing climate.
Water Research 46:1394-1407.
Carlson, R. E. 1977. Atrophic state index for lakes. Limnology and Oceanography 22:361-369.
Carvalho, L., C. A. Miller (nee Ferguson), E. M. Scott, G. A. Codd, P. S. Davies, and A. N. Tyler.
2011. Cyanobacterial blooms: Statistical models describing risk factors for national-scale
lake assessment and lake management. Science of The Total Environment 409:5353-
5358.
Cheung, M. Y., S. Liang, and J. Lee. 2013. Toxin-producing cyanobacteria in freshwater: A review
of the problems, impact on drinking water safety, and efforts for protecting public
health. Journal of Microbiology 51:1-10.
Chorus, I. 2001. Cyanotoxins: occurrence, causes, consequences. Springer, Berlin, Germany.
Chow, C. W. K., M. Drikas, J. House, M. D. Burch, and R. M. A. Velzeboer. 1999. The impact of
conventional water treatment processes on cells of the cyanobacterium Microcystis
aeruginosa. Water Research 33:3253-3262.
Coker, G. A., C. B. Portt, and C. K. Minns. 2001. Morphological and ecological characteristics of
Canadian freshwater fishes. Fisheries and Oceans Canada Burlington, Ontario.
Colby, P. J., G. R. Spangler, D. A. Hurley, and A. M. McCombie. 1972. Effects of eutrophication on
salmonid communities in oligotrophic lakes. Journal of the Fisheries Research Board of
Canada 29:975-983.
Cole, J. J., S. R. Carpenter, J. KitchelI, M. L. Pace, C. T. Solomon, and B. Weidel. 2011. Strong
evidence for terrestrial support of zooplankton in small lakes based on stable isotopes of
carbon, nitrogen, and hydrogen. Proceedings of the National Academy of Sciences
108:1975-1980.
70

-------
Cornett, R. J. 1989. Predicting changes in hypolimnetic oxygen concentrations with phosphorus
retention, temperature, and morphometry. Limnology and Oceanography 34:1359-
1366.
Cornett, R. J., and F. H. Rigler. 1980. The areal hypolimnetic oxygen deficit: An empirical test of
the model: Areal hypolimnetic oxygen deficit. Limnology and Oceanography 25:672-
679.
Cornett, R. J., and F. H. Rigler. 1984. Dependence of hypolimnetic oxygen consumption on
ambient oxygen concentration: Fact or artifact? Water Resources Research 20:823-830.
Coutant, C. C. 1985. Striped bass, temperature, and dissolved oxygen: A speculative hypothesis
for environmental risk. Transactions of the American Fisheries Society 114:31-61.
Coutant, C. C., and D. S. Carroll. 1980. Temperatures occupied by ten ultrasonic-tagged striped
bass in freshwater lakes. Transactions of the American Fisheries Society 109:195-202.
Daly, C., M. Halbleib, J. I. Smith, W. P. Gibson, M. K. Doggett, G. H. Taylor, J. Curtis, and P. P.
Pasteris. 2008. Physiographically sensitive mapping of climatological temperature and
precipitation across the conterminous United States. International Journal of
Climatology 28:2031-2064.
De Robertis, A., C. H. Ryer, A. Veloza, and R. D. Brodeur. 2003. Differential effects of turbidity on
prey consumption of piscivorous and planktivorous fish. Canadian Journal of Fisheries
and Aquatic Sciences 60:1517-1526.
Deevey, E. S. 1940. Limnological studies in Connecticut; Part V, A contribution of regional
limnology. American Journal of Science 238:717-741.
Demers, E., and J. Kalff. 1993. A simple model for predicting the date of spring stratification in
temperate and subtropical lakes. Limnology and Oceanography 38:1077-1081.
71

-------
Demott, W., and D. Muller-Navarra. 1997. The importance of highly unsaturated fatty acids in
zooplankton nutrition: evidence from experiments with Daphnia, a cyanobacterium and
lipid emulsions. Freshwater Biology 38:649-664.
Diaz, R. J. 2001. Overview of hypoxia around the world. Journal of Environmental Quality
30:275-281.
Dillon, P. J., and F. H. Rigler. 1974. The phosphorus-chlorophyll relationship in lakes. Limnology
and Oceanography 19:767-773.
Dokulil, M. T., and K. Teubner. 2000. Cyanobacterial dominance in lakes. Hydrobiologia 438:1-
12.
Dolman, A. M., J. Riicker, F. R. Pick, J. Fastner, T. Rohrlack, U. Mischke, and C. Wiedner. 2012.
Cyanobacteria and cyanotoxins: The influence of nitrogen versus phosphorus. PLoS ONE
7:e38757.
Downing, J. A., S. B. Watson, and E. McCauley. 2001. Predicting cyanobacteria dominance in
lakes. Canadian Journal of Fisheries and Aquatic Sciences 58:1905-1908.
Duane, S., A. D. Kennedy, B. J. Pendleton, and D. Roweth. 1987. Hybrid monte carlo. Physics
Letters B 195:216-222.
Dumont, H. J., I. V. de Velde, and S. Dumont. 1975. The dry weight estimate of biomass in a
selection of Cladocera, Copepoda and Rotifera from the plankton, periphyton and
benthos of continental waters. Oecologia 19:75-97.
Eaton, J. G., and R. M. Scheller. 1996. Effects of climate warming on fish thermal habitat in
streams of the United States. Limnology and Oceanography 41:1109-1115.
Effler, S. W., and S. M. O'Donnell. 2010. A long-term record of epilimnetic phosphorus patterns
in recovering Onondaga Lake, New York. Fundamental and Applied Limnology / Archiv
fur Hydrobiologie 177:1-18.
72

-------
Elser, J. J. 1999. The pathway to noxious cyanobacteria blooms in lakes: the food web as the
final turn. Freshwater Biology 42:537-543.
Elton, C. S. 1927. Animal Ecology. University of Chicago Press.
Evans, D. O., K. H. Nicholls, Y. C. Allen, and M. J. McMurtry. 1996. Historical land use,
phosphorus loading, and loss offish habitat in Lake Simcoe, Canada. Canadian Journal of
Fisheries and Aquatic Sciences 53:194-218.
Ferber, L. R., S. N. Levine, A. Lini, and G. P. Livingston. 2004. Do cyanobacteria dominate in
eutrophic lakes because they fix atmospheric nitrogen? Freshwater Biology 49:690-708.
Ferguson, R. G. 1958. The preferred temperature of fish and their midsummer distribution in
temperate lakes and streams. Journal of the Fisheries Research Board of Canada
15:607-624.
Filstrup, C. T., H. Hillebrand, A. J. Heathcote, W. S. Harpole, and J. A. Downing. 2014a.
Cyanobacteria dominance influences resource use efficiency and community turnover in
phytoplankton and zooplankton communities. Ecology Letters 17:464-474.
Filstrup, C. T., T. Wagner, P. A. Soranno, E. H. Stanley, C. A. Stow, K. E. Webster, and J. A.
Downing. 2014b. Regional variability among nonlinear chlorophyllphosphorus
relationships in lakes. Limnology and Oceanography 59:1691-1703.
Fischer, W. J., I. Garthwaite, C. O. Miles, K. M. Ross, J. B. Aggen, A. R. Chamberlin, N. R. Towers,
and D. R. Dietrich. 2001. Congener-independent immunoassay for microcystins and
nodularins. Environmental Science & Technology 35:4849-4856.
Galapate, R. P., A. U. Baes, and M. Okada. 2001. Transformation of dissolved organic matter
during ozonation: Effects on trihalomethane formation potential. Water Research
35:2201-2206.
73

-------
Gelman, A. 2006. Prior distributions for variance parameters in hierarchical models (comment
on article by Browne and Draper). Bayesian Analysis 1:515-534.
Gelman, A., and J. Hill. 2007. Data analysis using regression and multilevel/hierarchical models.
Cambridge University Press, New York, NY.
del Giorgio, P. A., and J. M. Gasol. 1995. Biomass distribution in freshwater plankton
communities. The American Naturalist 146:135-152.
Gloss, S. P., L. M. Mayer, and D. E. Kidd. 1980. Advective control of nutrient dynamics in the
epilimnion of a large reservoir. Limnology and Oceanography 25:219-228.
Godfrey, P. J. 1982. The eutrophication of Cayuga Lake: a historical analysis of the
phytoplankton's response to phosphate detergents. Freshwater Biology 12:149-166.
Gorham, E., and F. M. Boyce. 1989. Influence of lake surface area and depth upon thermal
stratification and the depth of the summer thermocline. Journal of Great Lakes Research
15:233-245.
Gowen, R., P. Tett, and K. Jones. 1992. Predicting marine eutrophication: the yield of chlorophyll
from nitrogen in Scottish coastal waters. Marine Ecology Progress Series 85:153-161.
Graham, J. L., K. A. Loftin, M. T. Meyer, and A. C. Ziegler. 2010. Cyanotoxin mixtures and taste-
and-odor compounds in cyanobacterial blooms from the midwestern United States.
Environ. Sci. Technol. 44:7361-7368.
Graham, N. J. D., V. E. Wardlaw, R. Perry, and J.-Q. Jiang. 1998. The significance of algae as
trihalomethane precursors. Water Science and Technology 37:83-89.
Hamilton, D. P., and S. F. Mitchell. 1996. An empirical model for sediment resuspension in
shallow lakes. Hydrobiologia 317:209-220.
Haney, J. F. 1987. Field studies on zooplankton-cyanobacteria interactions. New Zealand Journal
of Marine and Freshwater Research 21:467-475.
74

-------
Hanson, P. C., D. L. Bade, S. R. Carpenter, and T. K. Kratz. 2003. Lake metabolism: Relationships
with dissolved organic carbon and phosphorus. Limnology and Oceanography 48:1112-
1119.
Havens, K. E., and G. K. Niirnberg. 2004. The phosphorus-chlorophyll relationship in lakes:
Potential influences of color and mixing regime. Lake and Reservoir Management
20:188-196.
Heathcote, A. J., C. T. Filstrup, D. Kendall, and J. A. Downing. 2016. Biomass pyramids in lake
plankton: influence of Cyanobacteria size and abundance. Inland Waters 6:250-257.
Heiskary, S., and B. Wilson. 2008. Minnesota's approach to lake nutrient criteria development.
Lake and Reservoir Management 24:282-297.
Hessen, D. O. 2008. Efficiency, energy and stoichiometry in pelagic food webs; Reciprocal roles
of food quality and food quantity. Freshwater Reviews 1:43-57.
Hessen, D. O., B. A. Faafeng, P. Brettum, and T. Andersen. 2006. Nutrient enrichment and
planktonic biomass ratios in lakes. Ecosystems 9:516-527.
Holmes, R., R. Norris, T. Smayda, and E. J. F. Wood. 1969. Collection, fixation, identification, and
enumeration of phytoplankton standing stock. Recommended procedures for measuring
the productivity of plankton standing stock and related oceanic properties. Washington
(DC): National Academy of Sciences: 17-46.
Hondzo, M., and H. G. Stefan. 1993. Lake water temperature simulation model. Journal of
Hydraulic Engineering 119:1251-1273.
Hoyer, M. V., and J. R. Jones. 1983. Factors affecting the relation between phosphorus and
chlorophyll a in Midwestern reservoirs. Canadian Journal of Fisheries and Aquatic
Sciences 40:192-199.
75

-------
Hutchinson, G. E. 1938. On the relation between the oxygen deficit and the productivity and
typology of lakes. Internationale Revue der gesamten Hydrobiologie und Hydrographie
36:336-355.
Jacobson, P. C., H. G. Stefan, and D. L. Pereira. 2010. Coldwater fish oxythermal habitat in
Minnesota lakes: influence of total phosphorus, July air temperature, and relative
depth. Canadian Journal of Fisheries and Aquatic Sciences 67:2002-2013.
Jacquet, S., J.-F. Briand, C. Leboulanger, C. Avois-Jacquet, L. Oberhaus, B. Tassin, B. Vingon-Leite,
G. Paolini, J.-C. Druart, O. Anneville, and J.-F. Humbert. 2005. The proliferation of the
toxic cyanobacterium Planktothrix rubescens following restoration of the largest natural
French lake (Lac du Bourget). Harmful Algae 4:651-672.
Jenny, J.-P., P. Francus, A. Normandeau, F. Lapointe, M.-E. Perga, A. Ojala, A. Schimmelmann,
and B. Zolitschka. 2016. Global spread of hypoxia in freshwater ecosystems during the
last three centuries is caused by rising local human pressure. Global Change Biology
22:1481-1489.
Jeppesen, E., J. P. Jensen, C. Jensen, B. Faafeng, D. O. Hessen, M. S0ndergaard, T. Lauridsen, P.
Brettum, and K. Christoffersen. 2003. The impact of nutrient state and lake depth on
top-down control in the pelagic zone of lakes: A study of 466 lakes from the temperate
zone to the arctic. Ecosystems 6:313-325.
Jones, F. E., and G. L. Harris. 1992. ITS-90 density of water formulation for volumetric standards
calibration. Journal of research of the National Institute of Standards and Technology
97:335-340.
Jones, J. R., and R. W. Bachmann. 1976. Prediction of phosphorus and chlorophyll Levels in lakes.
Journal (Water Pollution Control Federation) 48:2176-2182.
76

-------
Jones, J. R., and M. F. Knowlton. 2005. Chlorophyll response to nutrients and non-algal seston in
Missouri reservoirs and oxbow lakes. Lake and Reservoir Management 21:361-371.
Jones, J. R., M. F. Knowlton, D. V. Obrecht, and J. L. Graham. 2011. Temperature and oxygen in
Missouri reservoirs. Lake and Reservoir Management 27:173-182.
Jones, J. R., D. V. Obrecht, B. D. Perkins, M. F. Knowlton, A. P. Thorpe, S. Watanabe, and R. R.
Bacon. 2008. Nutrients, seston, and transparency of Missouri reservoirs and oxbow
lakes: An analysis of regional limnology. Lake and Reservoir Management 24:155-180.
Kasprzak, P., J. Padisak, R. Koschel, L. Krienitz, and F. Gervais. 2008. Chlorophyll a concentration
across a trophic gradient of lakes: An estimator of phytoplankton biomass? Limnologica
- Ecology and Management of Inland Waters 38:327-338.
Keeler, B. L., S. A. Wood, S. Polasky, C. Kling, C. T. Filstrup, and J. A. Downing. 2015. Recreational
demand for clean water: evidence from geotagged photographs by visitors to lakes.
Frontiers in Ecology and the Environment 13:76-81.
Klumb, R. A., K. L. Bunch, E. L. Mills, L. G. Rudstam, G. Brown, C. Knauf, R. Burton, and F.
Arrhenius. 2004. Establishment of a metalimnetic oxygen refuge for zooplankton in a
productive lake Ontario embayment. Ecological Applications 14:113-131.
Knowlton, M. F., M. V. Hoyer, and J. R. Jones. 1984. Sources of variability in phosphorus and
chlorophyll and their effects on use of lake survey data. Journal of the American Water
Resources Association 20:397-408.
Knowlton, M. F., and J. R. Jones. 2000. Non-algal seston, light, nutrients and chlorophyll in
Missouri reservoirs. Lake and Reservoir Management 16:322-332.
Kritzberg, E. S., J. J. Cole, M. L. Pace, W. Graneli, and D. L. Bade. 2004. Autochthonous versus
allochthonous carbon sources of bacteria: Results from whole-lake 13C addition
experiments. Limnology and Oceanography 49:588-596.
77

-------
Lawrence, S. G., D. F. Malley, W. J. Findlay, M. A. Maclver, and I. L. Delbaere. 1987. Method for
estimating dry weight of freshwater planktonic crustaceans from measures of length
and shape. Canadian Journal of Fisheries and Aquatic Sciences 44:s264-s274.
Lee, W. C., and E. P. Bergersen. 1996. Influence of thermal and oxygen stratification on lake
trout hooking mortality. North American Journal of Fisheries Management 16:175-181.
Leibold, M. A., J. M. Chase, J. B. Shurin, and A. L. Downing. 1997. Species turnover and the
regulation of trophic structure. Annual Review of Ecology and Systematics 28:467-494.
Lewis, W. M. 1983. A revised classification of lakes based on mixing. Canadian Journal of
Fisheries and Aquatic Sciences 40:1779-1787.
Lewis, W. M., and W. A. Wurtsbaugh. 2008. Control of lacustrine phytoplankton by nutrients:
Erosion of the phosphorus paradigm. International Review of Hydrobiology 93:446-465.
Lienesch, P. W., M. E. McDonald, A. E. Hershey, W. J. O'Brien, and N. D. Bettez. 2005. Effects of a
whole-lake, experimental fertilization on lake trout in a small oligotrophic arctic lake.
Hydrobiologia 548:51-66.
Loftin, K. A., M. Meyer, F. Rubio, L. Kamp, E. Humphries, and E. Whereat. 2008. Comparison of
two cell lysis procedures for recovery of microcystins in water samples from Silver Lake
in Dover, Delaware, with microcystin producing cyanobacterial accumulations. Open-File
Report, United States Geological Survey, Reston, VA.
Mackenzie-Grieve, J. L., and J. R. Post. 2006. Thermal habitat use by lake trout in two contrasting
Yukon Territory lakes. Transactions of the American Fisheries Society 135:727-738.
Marcus, M. D., W. A. Hubert, and S. H. Anderson. 1984. Habitat Suitability Index Models: Lake
Trout (Exclusive of the Great Lakes). U.S. Fish and Wildlife Service.
78

-------
Mazumder, A. 1994. Phosphorus-chlorophyll relationships under contrasting herbivory and
thermal stratification: predictions and patterns. Canadian Journal of Fisheries and
Aquatic Sciences 51:390-400.
McCauley, E. 1984. The estimation of the abundance and biomass of zooplankton in samples. A
manual on methods for the assessment of secondary productivity in fresh waters:228-
265.
McCauley, E., J. A. Downing, and S. Watson. 1989. Sigmoid relationships between nutrients and
chlorophyll among lakes. Canadian Journal of Fisheries and Aquatic Sciences 46:1171-
1175.
McMahon, T. E., J. W. Terrell, and P. C. Nelson. 1984. Habitat Suitability Information: Walleye.
Page 43. Division of Biological Services, Fish and Wildlife Service, U.S. Department of the
Interior, Washington, DC.
Mellina, E., J. B. Rasmussen, and E. L. Mills. 1995. Impact of zebra mussel (Dreissena
polymorpha) on phosphorus cycling and chlorophyll in lakes. Canadian Journal of
Fisheries and Aquatic Sciences 52:2553-2573.
Molot, L. A., and P. J. Dillon. 1991. Nitrogen/phosphorus ratios and the prediction of chlorophyll
in phosphorus-limited lakes in central Ontario. Canadian Journal of Fisheries and Aquatic
Sciences 48:140-145.
Molot, L. A., P. J. Dillon, B. J. Clark, and B. P. Neary. 1992. Predicting end-of-summer oxygen
profiles in stratified lakes. Canadian Journal of Fisheries and Aquatic Sciences 49:2363-
2372.
Muller, B., L. D. Bryant, A. Matzinger, and A. Wiiest. 2012. Hypolimnetic oxygen depletion in
eutrophic lakes. Environmental Science & Technology 46:9964-9971.
79

-------
Miiller-Navarra, D. C., M. T. Brett, A. M. Liston, and C. R. Goldman. 2000. A highly unsaturated
fatty acid predicts carbon transfer between primary producers and consumers. Nature
403:74-77.
Norton, S. B., D. J. Rodier, W. H. van der Schalie, W. P. Wood, M. W. Slimak, and J. H. Gentile.
1992. A framework for ecological risk assessment at the EPA. Environmental Toxicology
and Chemistry 11:1663-1672.
Olson, J. R., and C. P. Hawkins. 2013. Developing site-specific nutrient criteria from empirical
models. Freshwater Science 32:719-740.
Pace, M. L., J. J. Cole, S. R. Carpenter, J. F. Kitchell, J. R. Hodgson, M. C. Van de Bogert, D. L. Bade,
E. S. Kritzberg, and D. Bastviken. 2004. Whole-lake carbon-13 additions reveal terrestrial
support of aquatic food webs. Nature 427:240-243.
Paerl, H. W., and T. G. Otten. 2013. Harmful cyanobacterial blooms: Causes, consequences, and
controls. Microbial Ecology 65:995-1010.
Paerl, H. W., J. T. Scott, M. J. McCarthy, S. E. Newell, W. S. Gardner, K. E. Havens, D. K. Hoffman,
S. W. Wilhelm, and W. A. Wurtsbaugh. 2016. It takes two to tango: When and where
dual nutrient (N & P) reductions are needed to protect lakes and downstream
ecosystems. Environmental Science & Technology 50:10805-10813.
Paerl, H. W., and J. F. Ustach. 1982. Blue-green algal scums: An explanation for their occurrence
during freshwater blooms. Limnology and Oceanography 27:212-217.
Park, S., M. T. Brett, E. T. Oshel, and C. R. Goldman. 2003. Seston food quality and Daphnia
production efficiencies in an oligo-mesotrophic Subalpine Lake. Aquatic Ecology 37:123-
136.
80

-------
Perkins, R. G., and G. J. C. Underwood. 2000. Gradients of chlorophyll a and water chemistry
along an eutrophic reservoir with determination of the limiting nutrient by in situ
nutrient addition. Water Research 34:713-724.
Persson, J., M. T. Brett, T. Vrede, and J. L. Ravet. 2007. Food quantity and quality regulation of
trophic transfer between primary producers and a keystone grazer (Daphnia) in pelagic
freshwater food webs. Oikos 116:1152-1163.
Phillips, G., O.-P. Pietilainen, L. Carvalho, A. Solimini, A. Lyche Solheim, and A. Cardoso. 2008.
Chlorophyll-nutrient relationships of different lake types using a large European
dataset. Aquatic Ecology 42:213-226.
Phillips, G., N. Willby, and B. Moss. 2016. Submerged macrophyte decline in shallow lakes: What
have we learnt in the last forty years? Aquatic Botany.
Pilotto, L. S., R. M. Douglas, M. D. Burch, S. Cameron, M. Beers, G. J. Rouch, P. Robinson, M. Kirk,
C. T. Cowie, S. Hardiman, C. Moore, and R. G. Attewell. 1997. Health effects of exposure
to cyanobacteria (blue-green algae) during recreational water-related activities.
Australian and New Zealand Journal of Public Health 21:562-566.
Plumb, J. M., and P. J. Blanchfield. 2009. Performance of temperature and dissolved oxygen
criteria to predict habitat use by lake trout (Salvelinus namaycush). Canadian Journal of
Fisheries and Aquatic Sciences 66:2011-2023.
Pollard, A. I., S. E. Hampton, and D. M. Leech. 2018. The promise and potential of continental-
scale limnology using the U.S. Environmental Protection Agency's National Lakes
Assessment. Limnology and Oceanography Bulletin 27:36-41.
Prairie, Y. T., C. M. Duarte, and J. Kalff. 1989. Unifying nutrientchlorophyll relationships in
lakes. Canadian Journal of Fisheries and Aquatic Sciences 46:1176-1182.
81

-------
Prepas, E. E., and D. O. Trew. 1983. Evaluation of the phosphorus-chlorophyll relationship for
lakes of the Precambrian Shield in western Canada. Canadian Journal of Fisheries and
Aquatic Sciences 40:27-35.
Pridmore, R. D., W. N. Vant, and J. C. Rutherford. 1985. Chlorophyll-nutrient relationships in
North Island lakes (New Zealand). Hydrobiologia 121:181-189.
Qian, S. S., and R. J. Miltner. 2015. A continuous variable Bayesian networks model for water
quality modeling: A case study of setting nitrogen criterion for small rivers and streams
in Ohio, USA. Environmental Modelling & Software 69:14-22.
Quinlan, R., A. M. Paterson, J. P. Smol, M. S. V. Douglas, and B. J. Clark. 2005. Comparing
different methods of calculating volume-weighted hypolimnetic oxygen (VWHO) in
lakes. Aquatic Sciences 67:97-103.
Quiros, R. 1990. Factors related to variance of residuals in chlorophyll  total phosphorus
regressions in lakes and reservoirs of Argentina. Hydrobiologia 200-201:343-355.
R Core Team. 2017. R: A language and environment for statistical computing. R Foundation for
Statistical Computing, Vienna, Austria.
Reynolds, C. S., and S. C. Maberly. 2002. A simple method for approximating the supportive
capacities and metabolic constraints in lakes and reservoirs. Freshwater Biology
47:1183-1188.
Rognerud, S., and G. Kjellberg. 1984. Relationships between phytoplankton and zooplankton
biomass in large lakes. SIL Proceedings, 1922-2010 22:666-671.
Sakamoto, M. 1966. Primary production by phytoplankton community in some Japanese lakes
and its dependence on lake depth. Archivfur Hydrobiologie 62:1-28.
Sartory, D. P., and J. U. Grobbelaar. 1984. Extraction of chlorophyll a from freshwater
phytoplankton for spectrophotometric analysis. Hydrobiologia 114:177-187.
82

-------
Scavia, D., J. David Allan, K. K. Arend, S. Bartell, D. Beletsky, N. S. Bosch, S. B. Brandt, R. D.
Briland, I. Daloglu, J. V. DePinto, D. M. Dolan, M. A. Evans, T. M. Farmer, D. Goto, H. Han,
T. O. Hook, R. Knight, S. A. Ludsin, D. Mason, A. M. Michalak, R. Peter Richards, J. J.
Roberts, D. K. Rucinski, E. Rutherford, D. J. Schwab, T. M. Sesterhenn, H. Zhang, and Y.
Zhou. 2014. Assessing and addressing the re-eutrophication of Lake Erie: Central basin
hypoxia. Journal of Great Lakes Research 40:226-246.
Scheffer, M., and E. H. van Nes. 2007. Shallow lakes theory revisited: various alternative regimes
driven by climate, nutrients, depth and lake size. Hydrobiologia 584:455-466.
Scheffer, M., S. Rinaldi, A. Gragnani, L. R. Mur, and E. H. van Nes. 1997. On the dominance of
filamentous cyanobacteria in shallow, turbid lakes. Ecology 78:272-282.
Schindler, D. W., R. E. Hecky, D. L. Findlay, M. P. Stainton, B. R. Parker, M. J. Paterson, K. G.
Beaty, M. Lyng, and S. E. M. Kasian. 2008. Eutrophication of lakes cannot be controlled
by reducing nitrogen input: Results of a 37-year whole-ecosystem experiment.
Proceedings of the National Academy of Sciences 105:11254-11258.
Schindler, D. W., M. A. Turner, and R. H. Hesslein. 1985. Acidification and alkalinization of lakes
by experimental addition of nitrogen compounds. Biogeochemistry 1:117-133.
Smeltzer, E., and S. A. Heiskary. 1990. Analysis and applications of lake user survey data. Lake
and Reservoir Management 6:109-118.
Smith, V. H. 1982. The nitrogen and phosphorus dependence of algal biomass in lakes: an
empirical and theoretical analysis. Limnology and Oceanography 27:1101-1112.
Smith, V. H., S. B. Joye, and R. W. Howarth. 2006. Eutrophication of freshwater and marine
ecosystems. Limnology and Oceanography 51:351-355.
Smith, V. H., and D. W. Schindler. 2009. Eutrophication science: where do we go from here?
Trends in Ecology & Evolution 24:201-207.
83

-------
Snucins, E. J., and J. M. Gunn. 1995. Coping with a warm environment: Behavioral
thermoregulation by lake trout. Transactions of the American Fisheries Society 124:118-
123.
S0ndergaard, M., J. P. Jensen, and E. Jeppesen. 2005. Seasonal response of nutrients to reduced
phosphorus loading in 12 Danish lakes. Freshwater Biology 50:1605-1615.
Spears, B. M., L. Carvalho, B. Dudley, and L. May. 2013. Variation in chlorophyll a to total
phosphorus ratio across 94 UK and Irish lakes: Implications for lake management.
Journal of Environmental Management 115:287-294.
Stan Development Team. 2016. Stan Modeling Language Users Guide and Reference Manual,
Version 2.14.0.
Stefan, H. G., X. Fang, D. Wright, J. G. Eaton, and J. H. McCormick. 1995. Simulation of dissolved
oxygen profiles in a transparent, dimictic lake. Limnology and Oceanography 40:105-
118.
Stefan, H. G., M. Hondzo, X. Fang, J. G. Eaton, and J. H. McCormick. 1996. Simulated long-term
temperature and dissolved oxygen characteristics of lakes in the north-central United
States and associated fish habitat limits. Limnology and Oceanography 41:1124-1135.
Stemberger, R. S. 1995. The influence of mixing on rotifer assemblages of Michigan lakes.
Hydrobiologia 297:149-161.
Stewart, I., A. A. Seawright, and G. R. Shaw. 2008. Cyanobacterial poisoning in livestock, wild
mammals and birds - an overview. Pages 613-637 in H. K. Hudnell, editor.
Cyanobacterial Harmful Algal Blooms: State of the Science and Research Needs. Springer
New York.
Tessier, A. J., and J. Welser. 1991. Cladoceran assemblages, seasonal succession and the
importance of a hypolimnetic refuge. Freshwater Biology 25:85-93.
84

-------
US EPA. 1986. Ambient Water Quality Criteria for Dissolved Oxygen. Office of Water, US
Environmental Protection Agency, Washington, DC.
US EPA. 1998. Guidelines for Ecological Risk Assessment. Risk Assessment Forum, Washington,
DC.
US EPA. 2000a. Nutrient Criteria Technical Guidance Manual, Lakes and Reservoirs. Office of
Water, US Environmental Protection Agency, Washington, DC.
US EPA. 2000b. Nutrient Criteria Technical Guidance Manual, Rivers and Streams. Office of
Water, US Environmental Protection Agency, Washington, DC.
US EPA. 2001. Nutrient Criteria Technical Guidance Manual, Estuarine and Coastal Marine
Waters. Office of Water, US Environmental Protection Agency, Washington, DC.
US EPA. 2007. Survey of the Nation's Lakes: Field Operations Manual. Office of Water, U.S.
Environmental Protection Agency, Washington, DC.
US EPA. 2008. Nutrient Criteria Technical Guidance Manual, Wetlands. Office of Water, US
Environmental Protection Agency, Washington, DC.
US EPA. 2010a. Using stressor-response relationships to derive numeric nutrient criteria. Office
of Water, U.S. Environmental Protection Agency, Washington, DC.
US EPA. 2010b. National Lakes Assessment: A Collaborative Survey of the Nation's Lakes. Office
of Water and Office of Research and Development, Washington, DC.
US EPA. 2011. 2012 National Lakes Assessment. Field Operations Manual. Office of Water, US
Environmental Protection Agency, Washington, DC.
US EPA. 2012a. 2012 National Lakes Assessment. Laboratory Operations Manual. U.S.
Environmental Protection Agency, Washington, DC.
US EPA. 2012b. 2012 National Lakes Assessment Site Evaluation Guidelines. Office of Water,
Washington, DC.
85

-------
US EPA. 2014. Framework for Human Health Risk Assessment to Inform Decision Making. Risk
Assessment Forum, Washingon, DC.
US EPA. 2015a. Preventing Eutrophication: Scientific Support for Dual Nutrient Criteria. Office of
Water, Washington DC.
US EPA. 2015b. Drinking Water Health Advisory for the Cyanobacterial Microcystin Toxins. Office
of Water, U.S. Environmental Protection Agency, Washington, DC.
US EPA. 2015c. Recommendations for Public Water Systems to Manage Cyanotoxins in Drinking
Water. Office of Water, U.S. Environmental Protection Agency, Washington, DC.
US EPA. 2019. Recommended human health recreational ambient water quality criteria or
swimming advisories for microcystins and cylindrospermopsin. Office of Water,
Washington DC.
Vadeboncoeur, Y., E. Jeppesen, M. J. V. Zanden, H.-H. Schierup, K. Christoffersen, and D. M.
Lodge. 2003. From Greenland to green lakes: Cultural eutrophication and the loss of
benthic pathways in lakes. Limnology and Oceanography 48:1408-1418.
Vanderploeg, H. A., S. A. Ludsin, S. A. Ruberg, T. O. Hook, S. A. Pothoven, S. B. Brandt, G. A. Lang,
J. R. Liebig, and J. F. Cavaletto. 2009. Hypoxia affects spatial distributions and overlap of
pelagic fish, zooplankton, and phytoplankton in Lake Erie. Journal of Experimental
Marine Biology and Ecology 381, Supplement:S92-S107.
Vanni, M. J. 1987. Effects of nutrients and zooplankton size on the structure of a phytoplankton
community. Ecology 68:624-635.
Wagner, T., P. A. Soranno, K. E. Webster, and K. S. Cheruvelil. 2011. Landscape drivers of
regional variation in the relationship between total phosphorus and chlorophyll in lakes.
Freshwater Biology 56:1811-1824.
86

-------
Walker, W. W. 1979. Use of hypolimnetic oxygen depletion rate as a trophic state index for
lakes. Water Resources Research 15:1463-1470.
Westrick, J. A., D. C. Szlag, B. J. Southwell, and J. Sinclair. 2010. A review of cyanobacteria and
cyanotoxins removal/inactivation in drinking water treatment. Analytical and
Bioanalytical Chemistry 397:1705-1714.
Wetzel, R. G. 2001. Limnology, Third Edition: Lake and River Ecosystems. Academic Press, San
Diego.
Wood, S. N. 2006. Generalized additive models: An introduction with R. CRC Press, Boca Raton,
FL.
Yuan, L. L., and J. R. Jones. 2019. A Bayesian network model for estimating stoichiometric ratios
of lake seston components. Inland Waters 9:61-72.
Yuan, L. L., and J. R. Jones. 2020a. Modeling hypolimnetic dissolved oxygen depletion using
monitoring data. Canadian Journal of Fisheries and Aquatic Sciences.
Yuan, L. L., and J. R. Jones. 2020b. Rethinking phosphorus-chlorophyll relationships in lakes.
Limnology and Oceanography.
Yuan, L. L., and A. I. Pollard. 2014. Classifying lakes to improve precision of nutrient-chlorophyll
relationships. Freshwater Science 33:1184-1194.
Yuan, L. L., and A. I. Pollard. 2017. Using national-scale data to develop nutrient-microcystin
relationships that guide management decisions. Environmental Science & Technology
51:6972-6980.
Yuan, L. L., and A. I. Pollard. 2018. Changes in the relationship between zooplankton and
phytoplankton biomasses across a eutrophication gradient. Limnology and
Oceanography 63:2493-2507.
87

-------
Yuan, L. L., and A. I. Pollard. 2019. Combining national and state data improves predictions of
microcystin concentration. Harmful Algae 84:75-83.
Zhang, F., J. Lee, S. Liang, and C. K. Shum. 2015. Cyanobacteria blooms and non-alcoholic liver
disease: evidence from a county level ecological study in the United States.
Environmental Health 14:41.
88

-------
6 Appendix A: State Case Study: Chlorophyll a-Microcystin
This case study in Iowa describes chlorophyll a (Chi a) and microcystin data collected by
the Iowa Department of Natural Resources (IDNR) that are combined with national data to
estimate a stressor-response relationship for the state (Yuan and Pollard 2019).
6.1	Data
Chl a measurements in Iowa were collected as part of an ambient lake monitoring
program conducted by IDNR. Water samples were collected with an integrated water column
sampler above the thermocline, when present, to a maximum depth of 2 meters (m) at the
deepest point of each lake. Lake water samples were collected in the summer (May-
September). An aliquot of the water sample was analyzed for Chl a in the laboratory by non-
acidified fluorometry after filtering water samples through GF/C filters. In a separate IDNR
monitoring program, microcystin concentrations (MCs) are sampled regularly at swimming
beaches in Iowa during the summer. This sampling effort includes state park beaches and locally
managed beaches across the state. MC was quantified in composite water samples collected at
nine different locations on three transects spanning the swimming beach. On each transect,
samples were collected at depths of 0.15, 0.5, and 1.0 m. Chl a and MC samples were matched
by lake and sampling date for use in the analysis. To maximize the available data, MC and Chl a
measurements collected within 1 day of each other were included as matched samples.
6.2	Statistical Analysis
The structure of a statistical model that accommodates data collected at different
spatial scales must be defined to ensure that the available data appropriately inform model
estimates. Consider the case of a large national data set of approximately 1,000 samples and a
state data set of approximately 50 samples. If the two data sets were pooled, the national data
would dominate the state data simply because of the significantly larger sample size, and the
state data would exert a weak influence on the model. In any single state, however, only about
20 samples from the national data might be available, and we would expect the state data to
dominate estimates. Defining a hierarchical structure in the model helps ensure that each data
set exerts the appropriate influence on the model results (Gelman and Hill 2007).
89

-------
A second issue that arises from combining data sets is that different measurements are
often collected in the different data sets. This problem is addressed in the national models by
modeling a comprehensive network of relationships between different parameters to take
advantage of the many different measurements available in the National Lakes Assessment
(NLA) data (Qian and Miltner 2015). Then, state data sets in which only a subset of
measurements were collected could still be feasibly modeled by informing specific aspects of
the network.
State data from Iowa were included in the national model and inform estimates of
relationships in the same network. As mentioned earlier, however, only Chi a and MC
measurements were available in the Iowa state data set. To prevent overspecifying the model,
the EPA selected one of the relationships in the network that could be refined with data from
the state. The relationship between Chi a and the relative biovolume of cyanobacteria relied
most heavily on empirical calibration, so it was selected for refinement with state data. More
specifically, the national model was revised so that model coefficients specific to each state
were estimated (Equation (20)).
E[logit(pCii)] = fim + f2,mchli + f3ikli]chlf	(39)
where different values of each of the coefficients were estimated for each state in the United
States, k. The values of the coefficients for each state were constrained by normal distributions
defined by the parameters, /Uf and Of. For example, the set of state-specific coefficients for fi
were drawn from a single normal distribution as follows:
/1~WormaZ(Ju^i,oyi)	(40)
Identical expressions can be written for the set of f2 values and f3 values. These distributions
constrained the range of possible values so estimates of those parameters computed with
relatively small sample sizes within individual states can "borrow" information from estimates
computed from other states (Gelman and Hill 2007).
Iowa state data were included in the model by noting that the data should inform
estimates of the coefficients only in the state of Iowa. That is, estimates of fi,f2, and f3 from
Equation (39) in Iowa are based on both the Iowa state data set and NLA data collected in Iowa.
In other states, estimates of the coefficients are based only on NLA data. The influence of Iowa
state data on the national distributions of the coefficients (as characterized by /Uf and Of) is
90

-------
limited because the data affect only one element of the overall distributions of coefficients.
Within the state of Iowa, however, the coefficients can be fit to maximize the predictive
accuracy of the overall relationship linking Chi a to MC for both Iowa data and NLA data
collected in Iowa, while remaining consistent with the range of possible values observed across
all states.
One final difference in fitting the Iowa state data is that several sources of variability
modeled separately in the national model (e.g., Si and s2 in equations (21) and (22)) are
combined into one combined estimate of residual variability. This combination of error terms
reflects the data available from Iowa, in which no laboratory replicates or direct measurements
of cyanobacterial biovolume were available. Hence, one lumped source of variability was
estimated.
For comparison, a simple bivariate model was fit using only IDNR data, in which MC was
modeled as a quadratic function of Chi a.
6.3 Results
A total of 556 samples of Chi a were measured at 28 lakes in Iowa. In some lakes, MC
concentrations were sampled at different beaches, so 686 observations of MC were matched to
the Chi a measurements.
In the revised draft national model with state-specific relationships between Chi a and
the relative biovolume of cyanobacteria, coefficients varied substantially among states. Because
coefficient values for quadratic relationships are not easily interpreted, the predicted mean
cyanobacterial-relative biovolume at a Chi a concentration of 20 microgram per liter (ng/L) is
plotted to visualize the range of variation among states (Figure 33). For comparison, among all
the national data, mean cyanobacterial-relative biovolume was 0.18 at Chi a concentration of 20
Hg/L. Systematic changes in cyanobacterial-relative biovolume with latitude or longitude were
not evident, but some regional differences were observed. For example, cyanobacterial-relative
biovolume with a Chi a concentration of 20 ng/L in northeast states was generally lower than
elsewhere, whereas in midwest states, it was somewhat higher.
91

-------
PropCyano
-120	-100	-80
Longitude
Figure 33. Variation in the relationship between Chi a and cyanobacterial-relative biovolume among
states. PropCyano: predicted mean relative biovolume of cyanobacteria at an illustrative Chi a = 20 ng/L.
As described previously, the relationship between Chi a and cyanobacterial-relative
biovolume in Iowa was adjusted to maximize the accuracy of the predicted MC. Inclusion of
Iowa data reduced the magnitude of the slope of the relationship between Chi a and
cyanobacterial-relative biovolume, but increased the intercept (Figure 34). So, higher values of
cyanobacterial-relative biovolume were observed at Chi a concentrations less than about 10
Hg/L. At higher Chi a concentrations, inclusion of Iowa state data did not substantively change
the predicted cyanobacterial-relative biovolume. Overall, in Iowa, the estimated relationship
between cyanobacterial-relative biovolume and Chi a was statistically indistinguishable from a
constant value (Figure 34). The addition of the state data also narrowed the range of the
credible intervals, as would be expected.
92

-------
(D
lO
O
CD
1	10	100
Chi a (ng/L)
Figure 34. Comparison of Chi o/cyanobacterial-relative biovolume relationships in Iowa. Filled gray: 90%
credible intervals for estimate of relationship using only NLA data collected in Iowa; solid and dashed
lines: mean and 90% credible intervals for estimate of relationship using both Iowa state and NLA data.
The predicted mean relationship between Chi a and MC in Iowa from the state-national
model closely followed the observed data (left panel, Figure 35), exhibiting a slight increase in
slope as Chi a concentration increased. The 90% prediction intervals shown in the plot were
based on the mean values of repeated random draws of 15 samples from the predicted
distribution to replicate the plotted observed data. The intervals were broad and included most
of the estimated mean values. The curvature observed in the simple bivariate fit between Chi a
and MC using only Iowa data was opposite of that observed from the state-national model,
predicting that the rate of increase in MC was lower at high Chi a concentrations than at low Chi
a concentrations (right panel, Figure 35). The 90% prediction intervals of this fit also included
most of the observed mean values, but qualitatively, the simple bivariate model did not match
the observed data as closely as did the state-national model.
93

-------
1	10	100
Chi a (ug/L)
1	10	100
Chi a (ug/L)
Figure 35. Comparison of predicted relationship between Chi a and MC for the state-national model (left
panel) and a model using only Iowa state data (right panel). Open circles: average MC concentration
computed in ~15 samples at the indicated Chi a; solid lines: mean relationship; dashed lines: 5th to 95th
percentiles of distribution of means of 15 samples drawn from predicted distribution.
Three features inherent to the model combining state and national data are likely
responsible for the improved predictions of observations in the Iowa data set. First, the network
of relationships specified in the national model define a nonlinear function linking Chi a to MC
that yielded a curved mean response (left panel, Figure 35). When only Iowa data are available,
no information regarding the functional form of the relationship between Chi a and MC is
known. Hence, it is difficult for the model to identify the correct shape of the curve. Indeed, the
concavity of the mean relationship identified by the model using only Iowa data (right panel,
Figure 35) was opposite of that estimated in the combined state-national model. Second, the
network of relationships in the state-national model provided information regarding unobserved
variables and relationships that could be used in lieu of direct observations. In this example, the
relationships between Chi a and total phytoplankton biovolume and between cyanobacterial
biovolume and microcystin were supplied by the national model. The Iowa-only model lacked
the benefit of the additional information, and hence, for this model a direct relationship
between Chi a and MC had to be estimated that aggregated the different causal linkages.
Finally, the hierarchical structure of the national model placed constraints on the range of
possible values for parameters estimated within each state. These constraints limited model
parameters for the state data set to values that were generally consistent with national
parameters.
94

-------
6.4 Criteria Derivation
Derivation of a draft recommended Chi a criterion based on decisions such as allowable
exceedance rate, targeted MC, and model uncertainty follows an identical process as described
for the national model. The model based on both IDNR data and NLA data yields a slightly
different relationship from the model estimated from only the national data (Figure 36). Slightly
greater uncertainty accompanies the estimate of the mean relationship in the Iowa-NLA model
than the estimate in the NLA-only model (see Figure 22), and that uncertainty is reflected in a
broader range of possible Chi a criteria. In the example shown in Figure 36, to maintain a
maximum exceedance rate of 1% of MC of 8 ng/L, the Chi a criterion associated with the lower
25th credible interval was 14 ng/L.
O
O
o
o
o
o
=	O
O	Q)
'tl	0 ^
O	"D o
Q_ C= C\J
OO O
1	10	100
Chi a (ng/L)
Figure 36. MC and Chi a measurements in Iowa. Top panel-open circles: observed values of microcystin
and Chi a for samples in which MC was greater than the detection limit; solid line: predicted MC that will
be exceeded 1% of the time for the indicated Chi a concentration; gray shading: 50% credible intervals
about mean relationship; horizontal and vertical line segments: candidate Chi a criteria based on targeted
MC. Bottom panel: proportion of samples for which microcystin was not detected in ~100 samples
centered at the indicated Chi a concentration.
95

-------
7 Appendix B: State Case Study: Chlorophyll a-Hypoxia
This case study in Missouri describes national and state data that are combined to refine
estimates of the relationship between chlorophyll a (Chi a) and deepwater hypoxia. As
described in Section 3.2.2, mean concentrations of dissolved oxygen below the thermocline
(DOm) decrease with time during the period of summer stratification. The sampling design of the
National Lakes Assessment (NLA) allowed for one visit to most of the lakes, so estimating
temporal changes in deepwater DOm in the national model required a space-for-time
substitution. State monitoring data collected during multiple visits to a smaller number of lakes
provided an opportunity to directly estimate temporal changes in DOm and to compare the
relationship between eutrophication and the rate of oxygen depletion with estimates from NLA
data.
7.1	Data
The Missouri data considered in this case study were collected an average of 3-4 times
per year by the University of Missouri (MU) from 1989 to 2007 as part of a statewide monitoring
effort. Samples were collected near the dam for each reservoir (herein referred to as lakes for
simplicity), where vertical profiles for temperature and DO concentration were measured (YSI
model 51B or 550A meters). Composite water samples from a depth of approximately 0.25
meter (m) were transferred to high density polyethylene containers, placed in coolers on ice,
and transported to the MU Limnology Laboratory. There, a 250-milliliter aliquot was filtered
(Pall A/E) for determination of total chlorophyll a via fluorometry following pigment extraction
in heated ethanol (Knowlton et al. 1984, Sartory and Grobbelaar 1984). A total of 198
measurements of DOm were available for analysis, collected at 20 different lakes over 62 unique
lake-year combinations.
7.2	Statistical Analysis
The same model equations used in the national model were applied to data collected in
Missouri:
(41)
96

-------
where DO0 is the value of D0m at the start of spring stratification, volumetric oxygen demand
(VOD)k is the net imbalance in the volumetric oxygen budget for lake k corresponding to sample
/' expressed as milligrams per liter per day of DO (Burns 1995), t, is the date that sample /' is
collected, and t0j- is the date of the beginning of stratification for lake-year j. Observed values of
DOm were assumed to be normally distributed with a standard deviation of about the
expected value. Note that, like the national model, VOD is assumed to be constant for each lake,
but the date of the beginning of stratification varied by year and lake. The model equation
specifying the relationship between Chi a, dissolved organic carbon (DOC), and lake depth and
VOD was the same equation used in the draft recommended national model (see Equation (14)).
As with the national model, saturation DO concentrations at the minimum temperature in
Missouri were used to set the value of DO0.
The treatment of DO measurements less than 2 milligrams per liter (mg/L) in the
Missouri data differed from the approach used in the NLA. From 2 to 14 measurements of DOm
greater than 2 mg/L were available in the Missouri data set for each of the lake-years included in
the model, so data were available to directly estimate temporal changes in DOm. Because data
were available at each lake before DOm approached zero, measurements of DOm that were less
than 2 mg/L could be excluded without biasing the model results.
Two models were run to explore the effects of combining Missouri data with the
national model. In the first model, only Missouri data were used, and in the second model, both
Missouri and NLA data were used to estimate the parameter values.
7.3 Results
The range of values spanned by each of the covariates differed between the two data
sets. Missouri measurements were collected over a broader range of days than the NLA,
whereas lakes sampled by the NLA covered a broader range of Chi a concentrations (Figure 37).
Variations in DOC concentrations and depths below the thermocline were also narrower in the
Missouri data than in the NLA data. Those differences in the range of observations were
reflected in the strength of correlation between each covariate and DOm. For Missouri, sampling
day was most strongly correlated with DOm, whereas for the NLA, sampling day exhibited the
weakest correlation with DOm. Instead, in the NLA data, Chi a, DOC, and the depth below the
thermocline were all more strongly correlated with DOm.
97

-------
o
Q
O
Q
OO^Oo
I2e 8
rrTT~
0.1
1	'1 1 T
1	10
Chi (ug/L)
8 o
n O
Or
O qqO
0Oc5P.c
o 0 S *
 o?.f
O 00 &o  fc  o 
-IS?
X)0 O
& 8)
n	1	1	1	1	!	r
0	10 20 30 40 50 60
Depth below thermocline (m)
Figure 37. Observed DOm vs. Chi a, sampling day, DOC, arid depth below the thermocline. Open circles:
NLA data; filled circles: Missouri.
The first day of stratification for Missouri lakes was generally earlier than for most of the
dimictic lakes considered in the national model (Figure 38), a finding that is consistent with the
fact that Missouri is located at the southern end of the geographic distribution of dimictic lakes
(see Figure 7). Both the Missouri-only model and the NLA-only model yielded similar estimates
of the relationship between Chi a and VOD [d2 in Equation (14)) (Figure 39), and the estimate
based on the combined data sets improved further on the precision. Estimates of coefficients
characterizing the relationship between VOD and depth below the thermocline (d3) and DOC (d4)
were much more precise in the NLA-only data set than in the Missouri-only data set. Hence, the
estimate based on the combined data set mainly reflects the trends in the NLA data.
98

-------
o
-sf"
o
o
_qI
Ln	~
I	1	1
50	100 150
Day of year
o
-sf"
o
o
-ntttf
k n
i	1	1
50	100 150
Day of year
Figure 38. Estimated first day of stratification for Missouri lakes (left panel) and NLA lakes (right panel).
Intercept (d-i)
Combined
NLA
MO
1	1	1	1	1
-0.09 -0.08 -0.07 -0.06 -0.05
Coefficient value
Depth below thermocline (d3)
Combined -
NLA -
MO -
-O-
"i	1	1	1	1	r
-0.0005 0.0005 0.0015
Coefficient value
Chi (d2)
I	1	1
-0.03 -0.02 -0.01
Coefficient value
DOC (d4)
o
o
	r 	!	1	!	1	!	1	
0.0025 -0.08 -0.04 0.00 0.04
Coefficient value
Figure 39. Model coefficients estimated for models for Missouri data, NLA data, and combined data. Thick
line segment: 50% credible intervals; thin line segment: 90% credible intervals; vertical dashed line:
coefficient value of zero.
Qualitatively, the model accurately represented the decrease in DOm over time in
different lakes (Figure 40). The effects of differences in the timing of spring stratification was
manifested as differences in the vertical position of each line, and in some lakes, substantial
variation was observed across years.
99

-------
150
Day of year
200	250
Day of year
140 160 180
Day of year
T	1	1	T~
140	160	180	200
Day of year
"I	1	1	1	1	1	1	r
80 100 120 140 160 180 200 220
Day of year
"I	1	T
80 100 120 140 160 180 200
Day of year
Figure 40. Relationships between day of year and DOm for six Missouri lakes. Different line and symbol
colors in each panel correspond to data collected within different years with at least three samples. Open
gray circles: other samples collected at each lake.
7.4 Criteria Derivation
The utility of combining Missouri and NLA data to inform decision-making is evident
when one considers the predicted relationship between Chi a and DOm calculated using
parameter estimates from the Missouri data and from the combined Missouri-NLA data set
(Figure 41). In the example shown, the relationship is calculated based on illustrative values for
other covariates (depth below thermocline at 13 m, DOC at 3.5 mg/L, and time between spring
stratification and sampling at 120 days). Because use of both data sets improves the precision of
model parameters, the resulting mean relationship is also estimated with increased precision
and a targeted Chi a concentration can be identified with greater confidence. In this example,
the 50% credible interval for the targeted Chi a concentration corresponding to DOm = 0 extends
from 10 to 16 ng/L when the combined model is used. When using only Missouri data, the
interval expands to 8-22 ng/L.
100

-------
LO
CO
CNJ
O
2
3
4 5 6
7 8
20
10
Chi a (ng/L)
Figure 41. Relationship between Chi a and DOm in an illustrative lake with depth below thermocline at
13m, DOC at 3.5 mg/L, and 120 days after spring stratification. Solid line: mean relationship; gray shading:
50% credible intervals about mean relationship from combined Missouri-NLA model; dashed line: 50%
credible intervals about mean relationship from Missouri-only model; dotted line: DOm = 0 mg/L.
101

-------
8 Appendix C: State Case Study: Total Nitrogen-Chlorophyll a
This case study in Iowa examines how combining locally collected measurements of
total nitrogen (TN) and chlorophyll a (Chi a) with the national draft models can refine
predictions calculated from these local data sets.
8.1	Data
Data used for this case study were collected by the Iowa Department of Natural
Resources (IDNR) as part of their routine monitoring program. For each lake in the data set, TN,
NOx, Chi a, and dissolved organic carbon (DOC) values were measured. A total of 968
observations collected at 31 different lakes were available for analysis.
8.2	Statistical Analysis
The same model formulation provided in Equation (36) was applied to the IDNR data,
expressing TN-dissolved inorganic nitrogen (-DIN) as the sum of a phytoplankton compartment,
modeled as fiChlkl, and a dissolved organic nitrogen (DON) component, modeled as f2DOCk2; and
nitrogen (N) bound to organic sediment (equation is repeated below):
E[TN - DIN] = faChl11 + DON + 0Snp = fcCM" + f2DOC + 0Snp	(42)
DOC measurements were available only at a small proportion of Iowa lakes, so the EPA
simplified the national model to the following form for modeling Iowa data:
E[TN - DIN] = fcCM" + u	(43)
where u is a lake-specific constant representing the contributions of DON and nonphytoplankton
organic suspended sediment (OSnp) in each lake to observed values of TN-DIN. Recall also that,
in the national model, the coefficient fi varied across states. With the IDNR data set, multiple
samples were collected from each lake, so the model could be refined further to estimate a
value of fi for each lake as follows:
lg(/ij) ~ Normal{nfUA, afl)	(44)
where the index, j, refers to different lakes, and the mean value iXfijA is computed for data
collected in Iowa.
102

-------
To examine the effects of considering local state data in the context of the national
model, two models were fit. In the first model, only IDNR data were used to estimate the
coefficients. In the second model, relationships were fit to both the IDNR data and NLA data
simultaneously. The exponent k was modeled as being the same in both the IDNR and NLA data,
while the coefficients fi for each lake were estimated with IDNR data and NLA data collected
within Iowa, and the value of iXfijA was constrained by the national distribution among all the
states in the NLA data.
8.3 Results
Data collected during the NLA in Iowa and by IDNR spanned similar ranges of Chi o, TN-
DIN, and DOC (Figure 42). The limiting relationship between Chi a and TN-DIN estimated using
only IDNR data approximated the lower edge of the cloud of points (gray shading) but were
estimated with more uncertainty than when estimated using both IDNR and NLA data (solid
lines). The mean limiting relationships between Chi a and TN-DIN estimated with the two
models were statistically indistinguishable from one another.
Figure 42. Chi a vs. TN-DIN in Iowa. Open circles: data collected by Iowa DNR -filled circles: data collected
by NLA in Iowa; solid lines: 95% credible intervals for limiting relationships between Chi a and TN-DIN
estimated using both NLA and IDNR data; shaded gray area: 95% credible intervals for limiting
relationships estimated using only IDNR data.
The root mean square (RMS) prediction error of log(TN-DIN) measurements in the IDNR
data was the same for the models using only IDNR data (RMS = 0.27) and the combined Iowa -
NLA data (RMS = 0.27), indicating that imposing national constraints on the parameter values
did not reduce the accuracy of predictions at the scale of the local state data. Uncertainty about
103

-------
estimates of the relationship between TN-DIN and Chi a for individual lakes was very similar
(example shown in Figure 43), indicating that a sufficient number of samples was available for
each lake to estimate the relationship without the information provided by the national model.
_l
CO
=L
Z
~ O
i 
z
H
10	100
Chi a (ng/L)
Figure 43. Chi a vs. TN-DIN in Beeds Lake, Iowa. Open circles: observed data; gray shading: 90% credible
intervals for predicted relationship based on only IDNR data; solid lines: 90% credible intervals for
predicted relationship using both IDNR and NLA data.
8.4 Criteria Derivation
Because of the higher number of samples collected within each lake in the IDNR data
set, unique relationships between TN-DIN and Chi a for each lake could be calculated, and those
relationships, in turn, can be used to derive numeric nutrient criteria (Figure 44). Variations
across lakes in DON and OSnp and in the coefficients of the modeled relationship yield
differences in the estimated relationship between TN-DIN and Chi a. Then, resulting TN ambient
criterion differ as well. For an illustrative target Chi a concentration of 15 micrograms per liter
(Hg/L), the mean ambient TN criterion for the lake shown in the left panel of Figure 44 was 750
Hg/L, while the TN criterion for the lake in the right panel was 1260 ng/L.
104

-------
0.9
T	1	1I I I I |
3 4 5 7 9
-1	1	1ii i i i |
20 30 40 60 80
T	1	1 I I I |
4 5 7 9
-1	1	1ii i i i |
20 30 40 60 80
200
10
Chi a (ng/L)
100
10
Chi a (ng/L)
100
Figure 44. Lake-specific criteria derivation using combined Iowa-NLA model for two different lakes in
Iowa. Open circles: observed values of TN-DIN and Chi a in Iowa for each lake; gray shading: 50% credible
intervals about the mean relationship; solid line: mean relationship calculated using mean DOC
concentration in lake; horizontal and vertical line segments: TN criterion calculation for illustrative Chi a
target of 15 ng/L.
105

-------
9 Appendix D: Operational Numeric Nutrient Criteria
Operationally, chlorophyll a (Chi a), total nitrogen (TN), and total phosphorus (TP)
criteria can be specified to account for the effects of sampling and temporal variability on
observed mean concentrations (Barnett and O'Hagan 1997). In most cases, the condition of a
lake will be assessed by examining a small number of samples and the uncertainty in the
estimation of the true seasonal mean value from those data will be determined by the number
of samples, the temporal variability of nutrient concentrations in the lake, and the inherent
sampling variability of the measurement. By examining historical data from many different
lakes, sampling variability associated with TN and TP can be estimated and "operational" criteria
can be specified to account for this variability with adjusted criterion magnitudes and by
adopting a frequency component that allows for some excursions of the specified magnitude.
Ambient monitoring of nutrient concentrations provides the basis for determining
whether a lake complies with the specified numeric nutrient criteria. Because of logistical and
resource restrictions, the number of water quality samples available at different lakes can vary
from a single grab sample to weekly or monthly samples throughout the sampling season.
Statewide monitoring designs also vary in how often a lake is visited in different years. For
example, a typical rotating basin design might sample the same lake once every 5 years,
whereas other lakes might be sampled every year. Because of the differences in the frequency
of sample collection, a statistical analysis of available monitoring data might be necessary to
accurately assess compliance with the numeric nutrient criteria. This appendix describes a
statistical approach for deriving operational or realizable criteria magnitude, duration, and
frequency components.
This document provides tools to compute numeric nutrient criteria expressed as
seasonal mean values. Those criteria implicitly assumed that a large number of samples are
available for characterizing the condition of each lake and that the uncertainty in the
computation of the mean value is small (Barnett and O'Hagan 1997), a condition that is usually
not satisfied by routine monitoring data. Operational criteria incorporate statistical uncertainty
in estimating environmental conditions from a much smaller number of samples. The statistical
approach recommended here requires that one estimate the sampling and temporal variability
of nutrient concentrations within lakes for which criteria are specified.
106

-------
A variety of approaches are available that account for within-lake variability when
defining operational criteria, but they should all be designed to consider that nutrient
concentrations vary in space (e.g., at different points on a lake) and in time. Both sources of
variability account for a distribution of nutrient concentrations that will arise when a lake is
repeatedly sampled. For example, if a single sample of TP was collected from one lake every
year, over 10 years, the distribution of values might be as shown in Figure 45, in which observed
concentrations range from 30 to 80 micrograms per liter (|a,g/L). Given this example, the relevant
water quality management question is whether the lake complies with its specified numeric
nutrient criteria. Here, if the relevant criterion is 60 |a,g/L, a methodical approach for assessing
compliance can enhance the utility of the criterion. This section provides one example of an
approach for accounting for sampling variability and defining "operational" nutrient criteria.
"3" |		
CO 				
>
o
c
0
=3 CNJ 
cr
CD
1
LL
o I 				
I	T	T	I	I	I	1
20	40	80	160
TP WL)
Figure 45. Example distribution of 10 TP measurements. Note that the horizontal axis is log-scaled.
Estimates of sampling variability are needed to inform decisions on operational criteria,
and those estimates can be computed from historical data. For this example, the EPA analyzed
TP data extracted from the Storage and Retrieval Data Warehouse (STORET) that had been
collected in the summers from 1990 to 2011. From those data, lakes were identified in the U.S.
with at least 5 years of nutrient data, yielding 25,056 samples collected from 846 different lakes.
A statistical model was then used to estimate variance in nutrient measurements across
different samples collected in the same year and from the same lake (within-lake variability). A
model was fit to TP measurements that explicitly estimated intra-annual and interannual
variability as follows:
107

-------
log (TPi) = am + bm + rt	(45)
where TP, is measured in sample i at site j and in year k. So, observed TP in a sample is modeled
as being log-normally distributed about a mean value that is the sum of an overall site mean (ay)
and a random effect of year (bk). The random effect of year is assumed to be normally
distributed with a mean value of 0 and a standard deviation of syeor, and the intra-annual
variance (r,) is modeled as a normal distribution with a mean of 0 and a standard deviation of
Ssampie Intra-annual variance not only includes contributions from traditional sources of sampling
variability (e.g., measurement uncertainty), but also includes variability that could be attributed
to differences in TP concentrations among different locations in a lake and differences in TP
concentrations one might observe over the course of a single sampling season. Hence, intra-
annual variance was expected to differ among different lakes, so, the overall distribution of
different values of ssamp/e was modeled as a half-Cauchy distribution (Gelman 2006).
Fitting this model to the TP data collected from STORET yielded a mean estimate of 0.16
for intra-annual variability of log(TP). Among different lakes in the data set, this value ranged
from 0.10 to 0.27, so sampling variability varied substantially among the lakes in the data set.
Estimating intra-annual variability from local data collected in the lake of interest would help
ensure that the estimate correctly reflects variability in the lake.
Once intra-annual variability for the lake or lakes of interest has been estimated, this
information can be combined with the relevant criterion for that lake to estimate a distribution
of nutrient concentration values that would be observed if the lake complied with the criterion.
For example, if the standard deviation of the intra-annual variability of log(TP) in a particular
lake is estimated as 0.16 and the relevant TP criterion for the lake is 60 ng/L, we can infer the
characteristics of the cumulative distribution of TP values that would be observed at the lake if it
were exactly complying with its criterion (Figure 46). Then, based on this distribution,
operational criteria can be derived. For example, one might define an operational criterion that
corresponds with the 10th percentile of the distribution (TP = 37 ng/L) and assert that a single
TP observation below that value indicates the probability that the mean TP concentration in the
lake is greater than 60 ng/L is less than 10%. That is, a lake with an observation below that
threshold is likely in compliance with the criterion. Conversely, one might define a criterion at
the 90th percentile of the distribution (TP = 96 ng/L) and assert that a single TP observation that
exceeds that value indicates the probability that the mean TP concentration is lower than
108

-------
60 ng/L is less than 10%. That is, any lakes with an observation that exceeds that threshold is
likely to be out of compliance with the criterion. Different water quality management outcomes
(e.g., additional sampling) could be triggered at different threshold concentrations. Also,
different operational criteria can be developed depending on probabilities of error that are
acceptable to environmental managers.
-1	"l	1	1	1	1l |
30 40 50 60 80
100
TP (ng/L)
Figure 46. Example of defining an operational criterion magnitude. Solid line: the cumulative probability of
observing a single sample TP lower than or equal to the indicated value if the true annual mean was
exactly equal to the criterion (TP = 60 |J.g/L); dashed line: the cumulative probability for the average of
four samples; black arrows: operational criteria for one sample; gray arrows: operational criteria
associated with four samples.
This analysis also highlights the relative benefits of collecting additional samples from
each lake. More specifically, the standard error (s.e.) on the estimate of a summer mean
concentration is as follows:
s. e. =
ssample
\fW
(46)
where N is the number of samples collected and sSamPie is the sampling variability of the nutrient
concentration. Hence, additional samples increase the precision with which the annual average
nutrient concentration can be estimated. In Figure 46, the dashed line shows the cumulative
probability distribution of mean values computed using four samples. Because of the reduction
in the standard error, assessments for compliance can be made with much greater confidence.
The same 10% probabilities used above for single samples yield operational criteria of 47 and 76
Hg/L, when applied to the case of four measurements (gray arrows in Figure 46). Information
and procedures regarding the use of operational criteria in assessment might be described in a
state's assessment methodology to accompany criteria specified in the water quality standards.
109

-------