ytDst%
I s
IJ
PRO^
Mercury Inputs and Cycling in
Devil s Lake, Wisconsin:
A Pilot Study for Conducting a Total Maximum
Daily Load (TMDL) Analysis for an
Atmospherically-Derived Pollutant
-------
Final Report
March2006
This page left blank.
-------
Final Report
March2006
Mercury Inputs and Cycling in Devils Lake,
Wisconsin:
A Pilot Study for Conducting a Total Maximum Daily Load
Analysis for an Atmospherically-Derived Pollutant
March, 2006
Watershed Branch (4503-T)
Office of Wetlands, Oceans and Watersheds
U.S. Environmental Protection Agency
1200 Pennsylvania Avenue, NW
Washington, DC 20460
Document posted at http://www.epa.gov/owow/tmdl/techsupp.html
-------
Final Report
March2006
Acknowledgments
This document was adapted from a report prepared for the U.S. Environmental Protection Agency (EPA)
under Contracts 68-C-01-022 and 68-W-03-028.
The Devil's Lake, Wisconsin, mercury TMDL pilot study was a collaborative effort between EPA's Office
of Water and Office of Air and Radiation; Wisconsin's Department of Natural Resources (WDNR), and
EPA Region 5. Ruth Chemerys, Dwight Atkinson, and Randy Waite, the EPA co-leads, would like to
thank Curt Pollman (Aqua Lux Lucis, Inc., formerly of Tetra Tech, Inc.), Reed Harris (Reed Harris
Environmental Ltd., formerly of Tetra Tech, Inc.), and Tom Myers (ICF International) for their excellent
work in developing this report and conducting the supporting analyses. We also thank Martin Burkholder
(WDNR), Todd Nettesheim (Great Lakes National Program Office), and Erin Newman (EPA Region 5) for
assistance in coordinating this project, and many others who contributed to the effort. Cover photos (top
two) are courtesy of Randy Waite, U.S. EPA; bottom photo courtesy of Carlton Nash, U.S. EPA.
Tetra Tech, Inc., and Systems Application International, Inc./ICF Consulting (SAI) coordinated in
conducting the pilot project analyses. The final report was compiled by Tetra Tech, with key components
on the air modeling prepared by SAL SAI conducted simulations of the atmospheric transport,
transformation, and deposition of mercury using the Regional Modeling System for Aerosols and
Deposition (REMSAD) model (Myers et al., 2006). Tetra Tech conducted simulations of the fate and
transport of the mercury within the aquatic environment using the Dynamic Mercury Cycling Model (D-
MCM; Harris et al., 2005). Initial results from the pilot project were subjected to external peer review in
2004, and this report reflects updates in response to peer review comments. Technical support documents
with additional details on the water and air modeling are available upon request. Data and other
information in this report reflect conditions at the time the pilot study was initiated, although information
on Wisconsin fish consumption advisories since the pilot began has been updated.
-------
Final Report
March2006
TABLE OF CONTENTS
1 INTRODUCTION 1
2 PROBLEM STATEMENT 4
2.1 Overview Of The Problem 4
2.1.1 Nature of the Devils Lake Mercury Problem and Selection as the
Pilot Study Site 4
2.2 Description Of The TMDL Process 5
2.2.1 Relationship of Pilot to TMDL Elements 7
2.3 Watershed Description 8
2.4 Existing Conditions and Summary of Monitoring Data 10
2.5 Identification of Pollutants Being Addressed and Why 11
3 NUMERIC TARGETS AND WATERSHED INDICATORS 13
3.1 Numeric Water Quality Standards 13
3.2 Fish Consumption Guidelines 14
3.3 Numeric Target for the Devils Lake Pilot Study 14
4 AIR AND WATERSHED MODELING: SOURCE ANALYSIS FOR
LOADINGS, MERCURY MASS BALANCE, LINKAGE OF STRESSORS
TO WATER QUALITY ENDPOINTS 16
4.1 Technical Approach 16
4.2 Atmospheric Transport and Deposition Modeling 16
4.2.1 Overview of the Air Deposition Model, REMSAD 16
4.2.2 Modeling Domain 21
4.2.3 Emission Inputs and Background Assumptions 21
4.2.4 Meteorological Inputs for Revised (2001) Modeling 23
4.2.5 Analysis of Year-to-Year Variation 25
4.2.6 Results and Model Performance 27
4.2.7 Source Attribution—Tagging 29
4.3 Aquatic Mercury Cycling Modeling 31
-------
Final Report
March2006
4.3.1 Model Description 31
4.3.2 Modeling Approach 34
4.3.3 Linkage of Mercury Loads to Fish Tissue Concentrations 47
5 UNCERTAINTY ANALYSIS 51
5.1 Atmospheric Modeling Uncertainty Analyses 51
5.1.1 Variations due to Meteorological Variability 52
5.1.2 Variations Due to REMSAD Model Uncertainty for Monte Carlo
Analysis 54
5.2 Aquatic Mercury Cycling Model Uncertainty Analyses 55
5.2.1 "Minimum-Maximum" Sensitivity Analysis Results on Effects of
Parameter Uncertainty on Endpoint Prediction 56
5.2.2 Monte Carlo Analysis of Effects of Parameter Uncertainty on
Endpoint Prediction 60
5.2.3 Effects of Annual Variablity in Atmospheric Deposition Rates on
Endpoint Predictions 64
6 DISCUSSION AND CONCLUSIONS 70
6.1 Atmospheric Modeling 70
6.2 Aquatic Modeling 71
6.3 Load-Target Relationships 72
6.3.1 Linearity of Response 72
6.3.2 Factors Affecting Uptake of Mercury Inputs 73
6.4 Uncertainty in Model Predictions 75
6.4.1 Uncertainty in REMSAD Predictions 75
6.4.2 Uncertainty in D-MCM Predictions 76
7 REFERENCES 78
//
-------
Final Report
March2006
LIST OF FIGURES
Figure 1. Site location map for Devils Lake, Wisconsin 9
Figure 2. Map of watershed basin boundaries for Devils Lake, Wisconsin. From
Krohelski and Batten, 1995 10
Figure 3. Flow of Control through Key Subroutines in a REMSAD Mercury Simulation. 19
Figure 4. Observed Annual Rainfall Totals (cm) for 2001 24
Figure 5. MM5-Derived Annual Rainfall Totals (cm) for 2001 24
Figure 6. Simulated and Observed Annual Total Wet Deposition at MDN Monitors
(gm/km2) 27
Figure 7. Simulated and Observed State Averages of Annual Total Wet Deposition at
MDN Monitors (gm/km2) 28
Figure 8. Simulated Contributions to Annual Total Deposition (Wet and Dry) of Mercury
at Devil's Lake, Based on the Revised 2001 Simulation. (Wet deposition
contributes 5.9 g/km2; dry deposition contributes 4.6 g/km2.) Simulated
concentrations and deposition are affected by uncertainties in model input data
and model formulation. See uncertainties discussion in Section 5 31
Figure 9. Contributions to Non-Background Mercury Deposition at the Devil's Lake
Site, Based on the Original 1998 Modeling. Simulated concentrations and
deposition are affected by uncertainties in model input data and model
formulation. See uncertainties discussion in Section 5 31
Figure 10. Mercury cycling in D-MCM 32
Figure 11. Estimated mean monthly surface water temperatures for Devils Lake (Data
source: Lathrop 1999a, Foster 1999a) 38
Figure 12. Estimated mean monthly precipitation for Devils Lake. (Source: Krohelski
and Batten 1995) 39
Figure 13. Estimated inflows and outflows for Devils Lake (derived from Krohelski and
Batten 1995) 39
Figure 14. Estimated loading of Hg(II) directly to lake surface, originating from wet and
dry deposition 40
Figure 15. Estimated inflow Hg(II) loads to Devils Lake. Loading fluxes have been
normalized to the lake surface area 42
Figure 16. D-MCM model calibration results for epilimnetic concentrations of Hg(II). .44
in
-------
Final Report
March2006
Figure 17. D-MCM model calibration results for epilimnetic concentrations of
methylmercury 45
Figure 18. D-MCM model calibration results for hypolimnetic concentrations of Hg(II).45
Figure 19. D-MCM model calibration results for hypolimnetic concentrations of
methylmercury 46
Figure 20. D-MCM model calibration results for walleye methylmercury concentrations
as a function of fish weight. Included are observations for May, 1986, April,
1987, and April 1993 46
Figure 21. Predicted Hg concentrations in age 5 walleye as a function of different long
term constant annual rates of wet and dry Hg(II) deposition to Devils Lake.
Two curves are shown: one shows the predicted response with no changing
inputs from the lake watershed; the other shows the predicted response
assuming that watershed loads of methyl and total mercury are reduced an
equivalent fractional amount as the total atmospheric deposition flux of Hg(II).48
Figure 22. Relative response in age 5 walleye methylmercury concentrations predicted
by D-MCM as a function of time. Plot compares two scenarios: the red curve
shows the response assuming a 3 50
Figure 23. Box and whisker plot of variations in monthly wet and dry deposition of Hg at
Devils Lake, Wisconsin. Variations in deposition based on CART analysis of
the effects of meteorological variability on REMSAD predictions of Hg wet
and dry deposition for 1996, and extending the predictions for monthly Hg
fluxes over the period 1990 to 2000 using ambient meteorology. Boxes show
median and the range of central 50% of the simulated values. Whiskers define
the range of values that lie within 1.5 times the interquartile range. Asterisks
define values that lie beyond 1.5 times but less than 3 times the interquartile
range. Open circles lie beyond 3 times the interquartile range (SPSS, 1998). 53
Figure 24. Results of Minimum/Maximum Sensitivity Analysis 59
Figure 25. Time-dependent results for Monte Carlo simulation of the 50% load reduction
scenario. Results show as a function of time the mean value, minimum,
maximum, and range encompassed by ± 2a for age 5 walleye. The 2a
distributions were computed by assuming that the simulation results were log-
normally distributed (n = 100). Load is reduced as a step function at l = 0
years. Note: lower limit defined by -2a from the mean is nearly identical to
and obscured by the minimum 63
Figure 26. Steady state (t = 150 years) results for the Monte Carlo simulations showing
the predicted changes in age 5 walleye tissue concentrations of methylmercury
as a function of load reduction. Load reductions simulated are 10, 25, 50, 75,
IV
-------
Final Report
March2006
and 95%. Plot shows the range of predicted walleye mercury concentrations
encompassed by 95% of the simulations in conjunction with expected dose-
response relationship from the calibrated model (see Figure 26). Plot also
shows the EPA tissue criterion of 0.3 ppm for reference 64
Figure 27. Frequency distribution of total mercury atmospheric deposition fluxes used to
simulate effects of year-to-year variations in atmospheric deposition rates on
age 5 walleye mercury concentrations. Number of years simulated = 500 65
Figure 28. Results from D-MCM simulation of effects of annual variability in
atmospheric deposition on age 5 walleye methylmercury concentrations. Plot
shows the results for both watershed scenarios (no response and instantaneous
response to changes in atmospheric inputs) 66
Figure 29. Results from D-MCM simulation of effects of annual variability in
atmospheric deposition on age 5 walleye methylmercury concentrations shown
in Figure 33 expressed as the percent cumulative distribution 67
Figure 30. Comparison of the magnitude of dynamic changes in atmospheric deposition
and age 5 walleye response. Walleye simulation assumes instantaneous
watershed response. Both parameters are shown as the response relative to the
mean over the 500-year simulation period 68
Figure 31. Coefficient of variation for methylmercury concentrations in walleye in
owing to year-to-year variations in atmospheric deposition as a function of age
class of fish 69
-------
Final Report
March2006
LIST OF TABLES
Table 1. General Characteristics of Devils Lake, Wisconsin 10
Table 2. Mercury water quality criteria for surface waters in Wisconsin. Derived from
WAC NR 105. All concentrations expressed as ng/L 13
Table 3. Comparison of Existing Conditions to Water Quality Endpoints 15
Table 4. Major Differences between Original 1998 REMSAD Simulation and the Revised
2001 REMSAD Simulation 17
Table 5. Mercury Chemical Mechanism in REMSAD, Version 8 20
Table 6. Summary of 2001 Mercury Emissions for Revised REMSAD Modeling
{emissions in tons/year) 21
Table 7. Average Speciation of Mercury Emissions 22
Table 8. Total Annual U.S. Pollutant Emissions (106 ton/year) in PM inventory: 2001
Base Year 22
Table 9. Key Classification Parameters for Wet and Dry Deposition {not necessarily in
order of importance) 25
Table 10. Deposition Estimates (g/km2) for Devil's Lake 26
Table 11. Statistical Evaluation of Simulated Total Annual Wet Deposition of Mercury
against Observed Data from MDN Sites for the Year 2001 28
Table 12. Compartments and mercury forms in D-MCM 33
Table 13. Summary of data inputs by major category and type 35
Table 14. Statistical properties of 10-year CART simulation of wet and dry deposition of
Hg at Devils Lake, WI. CART analysis based on initial REMSAD simulation
of mercury deposition in 1998. Table shows mean values by month for wet,
dry, and total deposition, as well as ln-transformed (geometric) means.
Untransformed fluxes expressed as ng/m2 53
Table 15. Statistical properties of 10-year CART simulation of wet and dry deposition of
Hg at Devils Lake, WI. CART analysis based on REMSAD simulation of
mercury deposition in 1998. Table shows standard deviation values by month
for arithmetic and ln-transformed wet, dry, and total deposition. Units for
untransformed fluxes are ng/m2 54
Table 16. Uncertainty in long-term estimates for wet, dry and total deposition of Hg at
Devils Lake, WI. Wet and dry estimates derived from CART analysis {n =
10), using standard error of the mean to estimate uncertainty. Untransformed
vi
-------
Final Report
March2006
fluxes in |j,g/m2-yr. Total deposition simulated based on uncertainty in average
wet and dry deposition, and assuming that the error in modeling with
REMSAD was 45.3% for each flux component 55
vii
-------
Final Report
March2006
LIST OF ABBREVIATIONS USED IN TEXT
CART
Classification and Regression Tree
CB-IV
Carbon-Bond IV
CMAQ
Community Multiscale Air Quality model
CWA
Clean Water Act
DOC
Dissolved organic carbon
FSL
NOAA Forecast Systems Laboratory
GIS
Geographic Information System, computer mapping system
GMT
Greenwich Mean Time (5 hours earlier than USA east coast time)
Hg
The chemical notation for the element mercury, derived from the Greek
Hydrargyrum
Hg(0)
Elemental mercury, the silvery metal liquid at room temperature
Hg(II)
divalent mercury, a form of inorganic mercury; mercuric ion, RGM
Hg(p)
Atmospheric particulate-associated mercury
Hgt, Hgtot
Total mercury, i.e. lab analysis of all forms of mercury
mb
millibars of pressure in the atmosphere
MCM
Mercury Cycling Model
MDN
Mercury Deposition Network, a sub-network of the National
Atmospheric Deposition Program
MeHg
Methylmercury
METAALICUS
Mercury Experiment To Assess Atmospheric Loading in Canada
and the United States
MSRTC
Mercury Study Report to Congress (USEPA 1997)
MWC
Municipal Waste Combustors, (aka., Municipal Solid Waste
Incinerator)
MWI
Medical Waste Incineration
NADP
National Atmospheric Deposition Program
ng/L
nanograms per liter, unit of measure of concentration, ppt
NCEP
NOAA National Center for Environmental Prediction
NGM
Nested Grid Model
NOAA
National Oceanic and Atmospheric Administration
NOAA - ARL
National Oceanic and Atmospheric Administration - Air Resources
T oK
NPDES
Li3D
National Pollutant Discharge Elimination System
NTI
National Toxics Inventory
NWS
National Weather Service
OAQPS
EPA Office of Air Quality Planning and Standards
PPb
parts per billion, a unit of measure of concentration,
ppm
parts per million, a unit of measure of concentration
ppt
parts per trillion, a unit of measure of concentration
RAMS
Regional Atmospheric Modeling System
RELMAP
Regional Lagrangian Model of Air Pollution
REMAP
Regional - Environmental Monitoring Assessment Program
REMSAD
Regional Modeling System for Aerosols and Deposition
RUC
Rapid Update Cycle
RGM
Reactive Gaseous Mercury
TGM
Total Gaseous Mercury
TMDL
Total Maximum Daily Load
UAM-V
Urban Airshed Model
USEPA
U.S. Environmental Protection Agency
USFWS
U.S. Fish and Wildlife Service
USGS
U.S. Geological Survey
WDH
Wisconsin Department of Health
VIII
-------
Final Report
March2006
WDNR Wisconsin Department of Natural Resources
WI Wisconsin
l-ig/L micrograms per liter, a unit of measure of concentration, ppb
l-icq/L microequivalents per liter, a unit of measure of concentration
CB micro-Carbon Bond
-------
Final Report
March2006
This page left blank.
-------
Final Report
March2006
1 INTRODUCTION
In 1999 EPA began a pilot study to investigate the relationship between air emissions of mercury and water
quality impacts. Two sites were selected for analysis - Devils Lake in Wisconsin, and a portion of the
Florida Everglades. The study was conducted as a voluntary and cooperative effort with the States of
Wisconsin and Florida.
The underlying goal of the pilot study was to examine modeling and other technical approaches for taking
atmospheric sources into account when developing total maximum daily loads (TMDLs) for mercury,
especially how to integrate the results from atmospheric modeling with aquatic cycling modeling.
Although the pilot study is not a TMDL, it does include as a key second objective an analysis of what
uncertainties should be considered when developing a mercury TMDL. TMDLs are a product of The Clean
Water Act (CWA) Section 303(d), which requires States to: (1) identify and list waterbodies where State
water quality standards are not being met; and (2) establish "Total Maximum Daily Loads" for these
waters. TMDLs specify the amount of a pollutant that may be present in a water body and still allow it to
meet State water quality standards. TMDLs allocate pollutant loads among pollution sources, and include a
margin of safety that accounts for the uncertainty in the relationship between loadings of the pollutant and
the response of the system. The pilot study for the Florida Everglades was completed in 2002 (FDEP,
2002). This study presents the results from the pilot TMDL study for Devils Lake.
Mercury represents a good candidate pollutant for assessing how to conduct a TMDL analysis for an
impaired water body where the pollutant of interest is largely or wholly atmospheric in origin. Many of the
water bodies in the United States with elevated mercury concentrations of concern in biota are otherwise
pristine seepage lakes that have no known sources of mercury within their watersheds and derive the
majority of their inputs of mercury from atmospheric deposition directly to their surfaces. Moreover, there
is growing evidence establishing the link between atmospheric inputs of mercury and the aquatic cycling of
mercury. For example, wet deposition of mercury in southwestern Sweden, apparently in response to
declines in European emissions, has decreased by perhaps 50% since the late 1980's (Munthe et al., 2001).
At the same time, both mercury accumulation rates in lake sediments and fish concentrations show
evidence of recent declines (Bindler et al., 1991 and Johansson, 2001). In Canada, early results from the
Mercury Experiment To Assess Atmospheric Loading in Canada and the United States (METAALICUS)
indicate that the aquatic cycling of mercury responds very rapidly to new inputs of mercury (R. Harris,
pers. comm.). In Florida, declines in local emissions of mercury approximating 90% have been paralleled
by declines in mercury in fish within portions of the Florida Everglades approximating 80% (Pollman and
Porcella, 2003).
Mercury is regarded as a pollutant that operates at global, regional, and local scales (Fitzgerald et al.,
1998). Mercury occurs in primarily two oxidation states in the environment - elemental mercury (Hg°) and
Hg2+ - and its global cycling reflects the fundamental difference in the physical and chemical
characteristics of these two species. Elemental mercury vaporizes easily and and is generally believed tp
react slowly in the atmosphere, leading to large transport distances once it is released and an atmospheric
1
-------
Final Report
March2006
residence time approximating one year (Hall, 1995; Schroeder and Munthe, 1998; Lamborg et al., 2002)1.
Hg2+, on the other hand, is quite reactive and rapidly deposited once it is either released into or formed in
the atmosphere. Once it enters aquatic ecosystems, Hg2+ can become methylated given the appropriate
biogeochemical conditions. It is this form of mercury that ultimately poses the greatest concern to human
and ecological health. The primary route of methylmercury exposure to humans, which bioconcentrates in
aquatic food chains by factors in excess of a million, is through consumption of fish. Methylmercury that
is consumed in fish can cause harm to the brain and the developing fetus because it readily crosses the
placenta and blood brain and is neurotoxic. Developing fetal nervous systems are particularly sensitive to
methylmercury exposure, which can cause mental retardation and cerebral palsy (Chan et al., 2002). It also
has been linked to reduced fecundity (Arakawa et al., 2004) and a reduction in fine motor skills (Barbone et
al., 2004).
Emissions of mercury into the atmosphere derive from both natural and anthropogenic sources. Natural
emissions include volcanic releases and diffusion from ore bodies (Nriagu, 1979), while anthropogenic
emissions are related to sources such as stationary combustion processes and non-ferrous metal production
(Pacyna and Pacyna, 2002). Anthropogenic sources are believed to dominate natural sources by
approximately a ratio of 2:1. For example, of an estimated total annual global emission flux of 6,386
Mg/yr, perhaps 2,127 Mg/yr comes from direct anthropogenic emissions, 1,064 Mg/yr derives from natural
terrestrial and oceanic sources, an additional 2,130 and 1,065 Mg/yr derives from re-emissions of
previously deposited anthropogenic and naturally emitted mercury respectively (Seigneur et al., 2003).
One critical area of scientific uncertainty is how to incorporate the varying spatial scales of the atmospheric
mercury cycle into an understanding of mercury bioaccumulation in aquatic systems. Whether atmospheric
inputs of mercury to a specific waterbody derive from global scale sources, local sources, or regional
sources intermediate in proximity is in part related to the prevailing meteorology influencing the receptor
water body. In addition, because the geographic scales of mercury transport differ greatly according to the
chemical species under consideration, the chemical speciation of the emissions from a source greatly
influences whether the scale of influence of its emissions on atmospheric deposition fluxes will be
predominantly local, or more globally distributed. Explicitly tracking source emissions from a given area
("source tagging") using the Regional Modeling System for Aerosols and Deposition (REMSAD) model
was explored in this pilot study as one method for resolving the source attribution question.
This pilot study also explicitly examined the inherent uncertainty in predicted fish mercury concentrations
owing to uncertainty in the inputs to the model used to simulate the aquatic cycling of mercury in Devils
Lake - the Dynamic Mercury Cycling Model (D-MCM). The effects of parameter uncertainty for both
"forcing function" type variables (e.g., hydrologic variables and atmospheric loadings) as well as model
1 Hedgecock and Pirrone (2004) modeled the residence time of Hg(0) in the marine boundary layer using recently
measured gaseous Hg(0) oxidation reactions. Modeled residence times were far lower than the currently accepted
value of ca. 1 yr. For example, the residence time between the equator and 60° N was ca. 10 days during typical
summer conditions. The oxidation reactions considered by Hedgecock and Pirrone depend heavily upon bromine
production and loss in the atmosphere which takes place over marine waters and near coastal environments and do not
generally apply to continental air masses.
2
-------
Final Report
March2006
rate constants were explored using Monte Carlo analysis. These results, in conjunction with parameter
sensitivity analyses, are used to identify which uncertainties need to be considered when conducting a
mercury TMDL. In addition, the effects of natural year-to-year variations in atmospheric loadings on
variations in fish tissue mercury concentrations were examined. Inherent in this latter analysis is the
recognition that, even under a regime of constant mercury emissions, meteorologic variations from one
year to the next will result in varying atmospheric loads to Devils' Lake. For this pilot study, classification
and regression tree analysis (CART), was used in conjunction with REMSAD output for a single
simulation year to statistically identify the most important meteorological parameters and conditions
associated with certain ranges of daily wet and dry mercury deposition and then ultimately infer year-to-
year variations in mercury deposition due to changing meteorological conditions across a ten-year period.
3
-------
Final Report
March2006
2 PROBLEM STATEMENT
2.1 Overview Of The Problem
High concentrations of mercury in piscivorous fish in Wisconsin lakes were identified as early as 1983
(Wiener, 1983). By the mid-1980's, it became clear that mercury levels were typically greatest in
softwater, acidic lakes found in northcentral Wisconsin (Wiener, 1987). By 1987, the State of Wisconsin
posted health advisories for individual lakes recommending that human consumption be limited or avoided
or avoided based on observed concentrations exceeding 0.5 and 1.0 parts per million (wet weight) in edible
muscle offish (WDNR and WDH, 1987).
In 2001, the State of Wisconsin adopted a statewide advisory for most inland waters regarding fish
consumption and mercury exposure (WDNR, 2003).2 This action was predicated in part on two findings.
First was the publication of the National Research Council's report, "Toxicological Effects of
Methylmercury" (2002), which concluded that the standard the U.S. Environmental Protection Agency's
(EPA) reference dose3 is appropriate to determine whether to issue consumption advice for fish. The use of
this new reference dose required that consumption advice be issued when fish exceeded 0.05 ppm mercury.
While mercury concentrations differ by species, nearly all of the 1,200 water bodies tested to date within
Wisconsin contain some species that tend to have fillet mercury concentrations in excess of 0.05 ppm. As a
result, most - but not all - inland waters carry the same statewide advice that varies by species, ranging
from "unlimited" to "do not eat", recommending how much fish people can safely eat while reducing their
risk of exposure to mercury. The advice differs for different groups of people (WDNR, 2003). For waters
where higher concentrations of mercury are found in fish, more stringent advice or exceptions to the
general statewide advice are provided.
2.1.1 Nature of the Devils Lake Mercury Problem and Selection as the Pilot
Study Site
In 1998, Devils Lake was listed on the State of Wisconsin 303(d) list of impaired waters because of a fish
consumption advisory relating to high concentrations of mercury for walleye between 10 and 30 inches (25
to 76 cm).4 The current advisory for Devils Lake (based on the statewide advisory) is primarily based on
2 Before 2001, Wisconsin issued fish consumption advisories only for specific waters. In 2001, the state began issuing statewide
consumption advice for all waters, with exceptions or special advice for waters with higher fish mercury concentrations. For further
information on Wisconsin fish consumption advisories, see http://dnr.wi.gov/fish/consumption/generaladvice.html.
3 A reference dose is defined as "mi estimate of a daily exposure to the human population (including sensitive subpopulations) that is
likely to be without a risk of adverse effects when experienced over a lifetime."
4 Prior to 2000, Devil's Lake did have specific fish consumption advice for walleye ranging from 13 meals per year to "do not eat."
Under the approach adopted in 2001, Devil's Lake falls under the general statewide advice, rather than the special advice, as fish
mercury concentrations do not exceed the criteria for special advice. The state's 2008 303(d) listing methodology indicates that
waters are included on the 303(d) list when added to the special advisory, and de-listed when the special advisory no longer applies
(see http://dnr.wi.gov/org/water/wm/wqs/303d/2008/2008methodology.htmV
4
-------
Final Report
March2006
limiting exposure to high concentrations of mercury via ingestion of fish to developing fetuses and children
under the age of 15, as well as women of childbearing years.
Although WDNR (2003) cited 84 waterbodies with more restrictive advisory limits than Devils Lake, there
are a number of factors considered in deciding to use Devils Lake for this pilot exercise. First, in order to
conduct the pilot TMDL, it was desirable to use a site with existing data, instead of having to collect new
data. The Devils Lake site was selected after an initial data gathering effort indicated that a comparatively
large amount of data relating to mercury monitoring (relative to other lakes with posted consumption
advisories) was available for the lake. These data include relatively good mercury emissions, deposition,
water column chemistry and biota data for Devils Lake. The Wisconsin emissions inventory includes data
on mercury emissions reported by individual facilities for 1993 and 1996. In addition, mercury
concentrations and fluxes in bulk deposition was monitored between 1995 and 1999. Another site that
monitored ambient air concentrations of both total gaseous and aerosol mercury for the period 12/94 -
12/96 is located 40 miles from the Lake. In addition, the Mercury Deposition Network (MDN) maintains a
suite of six stations throughout Wisconsin to monitor wet deposition of mercury, including a site at Devils
Lake that began monitoring in January 2001. Water column mercury chemistry and fish tissue data using
appropriate ultra-trace metal sampling and analytical techniques have been collected by the WDNR. These
data include water column mercury chemistry data for 1994-1996 (whole water total and methyl Hg;
particulate total and methyl Hg), and data on mercury levels in walleye, blue gills, shiner, and Daphnia.
Second, Devils Lake has no known sources of mercury originating within the watershed. Atmospheric
deposition is considered the only pathway through which mercury enters the lake and its watershed, and
there are a number of facilities near Devils Lake that emit mercury, including a coal-fired power plant in
Portage (22.5 kilometers [km] to the east/northeast); boilers, incinerators, and power plants in Madison (32
km SSE); a chlor-alkali plant and paper mills approximately 110 km to the north; and coal-fired boilers
145-161 km to the west/southwest.
Third, Devils Lake is a reasonably well-characterized aquatic ecosystem. A hydrologic budget for the lake
- a critical prerequisite for developing a mercury mass balance - was developed by Krohelski and Batten
(1995). A two-year water intensive lake quality study was conducted in 1986 and 1987 (basic limnological
measurements including pH and alkalinity, plus fish, macrophytes, and plankton data). General
limnological sampling also was conducted for the period 1988-1998.
2.2 Description of the TMDL Process
Water quality standards are established to protect the designated uses of Wisconsin's waters. When States,
Tribes or local communities identify problems in meeting water quality standards, a TMDL can be a
framework for addressing those problems. The purpose of this pilot study is to explore the utility of using
atmospheric and mercury-cycling models within the TMDL framework and to provide the stakeholders
with technical information that may be used to develop a water quality plan to address mercury issues in
water bodies within Wisconsin for which atmospheric deposition of mercury is considered the source cause
of impairment.
Section 303(d) of the Clean Water Act (CWA) requires states to identify the waters for which the effluent
limitations required under the National Pollutant Discharge Elimination System (NPDES) or any other
enforceable limits are not stringent enough to meet any water quality standard adopted for such waters.
5
-------
Final Report
March2006
The states must also prioritize these impaired water bodies for TMDL development, taking into account the
severity of the pollution and the beneficial uses of the waters.
A TMDL represents the maximum amount of a given pollutant that a particular waterbody can assimilate
without exceeding surface water standards. The TMDL can be expressed as the total mass of pollutant that
can enter the water body within a unit of time. For this pilot TMDL, it is the total mass of atmospherically
and water-borne mercury that enters Devils Lake. In most cases, the TMDL determines the allowable mass
per day of a pollutant and apportions it among the various pollution sources in the watershed as waste load
(i.e., point source discharge) and load (i.e., non-point source) allocations. The TMDL also accounts for
natural background sources (e.g., atmospheric deposition derived from global sources) and provides a
margin of safety.
Although this document is not a TMDL determination per se, the elements of TMDLs as described in
USEPA guidance are described here in order to put this pilot effort in context. Specifically, the
components of the pilot project are compared with the TMDL elements, and the report describes where
additional work may be needed to fully develop a TMDL using this pilot project results. Some of these
elements are required under the Clean Water Act, while others are elements recommended in USEPA
guidance. This information can be summarized as follows:
1. Identify waterbody, pollutant of concern, pollutant sources, and priority ranking. The TMDL
should identify all pollutant sources, both point and nonpoint sources, as well as the the magnitude
and location of sources. The TMDL should also describe any important assumptions made in
developing the TMDL.
2. Describe water quality standards and numeric water quality target: The TMDL must include
a description of the applicable water quality standards, including designated uses to be protected,
numeric and/or narrative criterion, and antidegradation policy. The TMDL must also include a
numeric target used to determine whether the water quality standard is attained. Where the
pollutant or target is different from the water quality standard, the TMDL must describe the
linkage(s) between the target and water quality standards.
3. Identify reductions for point and nonpoint sources of pollution. The TMDL must identify
the portions of the loading capacity allocated to point sources (wasteload allocation) and the
portion of the loading capacity allocated to nonpoint sources (load allocation).
4. Describe the linkage between water quality endpoints and pollutants of concern. The TMDL
must identify the loading capacity of the waterbody for the pollutant of concern. The TMDL
should also describe the method used to determine the relationship between the numeric targets
and the pollutant(s) of concern, usually a water quality model.
5. Identify margin of safety, seasonal variations, and critical conditions. The TMDL must
include a margin of safety that accounts for lack of knowledge regarding the relationship between
the allocations and water quality standards. The TMDL must also account for seasonal variations.
6. Provide implementation recommendations for pollutant reduction actions and a monitoring
plan. EPA does not require implementation plans as part of TMDLs, although some states do
require such plans. However, EPA encourages States to provide an implementation plan. A
6
-------
Final Report
March2006
monitoring plan is also recommended in order to track TMDL effectiveness, especially where
management actions will be phased in over time and to assess the achievement and validity of the
pollutant reduction goals.
7. Include an appropriate level of public involvement in the TMDL process. This is usually met
by publishing public notice of the TMDL, circulating the TMDL for public comment, and holding
public meetings in local communities. Public involvement must be documented in the state's
TMDL submittal to USEPA.
The required elements of a TMDL are loading capacity, wasteload allocation, load allocation, margin of
safety, and seasonal variation (see http://www.epa.gov/owow/tmdl/guidance/final52002.html').
2.2.1 Relationship of Pilot to TMDL Elements
This pilot analysis addresses of five of the seven elements of the TMDL process, described below. Because
this is a pilot TMDL study, not all of the elements of a TMDL have been completely addressed in this
report.
Identify waterbody, pollutant of concern, pollutant sources, and priority ranking. Devil's Lake
characteristics and conditions are described in Sections 2.3, 2.4, and 2.5. The primary source of mercury to
Devils Lake and its watershed is atmospheric deposition. Source attribution based on REMSAD modeling
of emissions, transport, and deposition indicates that background sources (i.e., sources outside the model
domain of southern Canada, the 48 conterminous states of the US, and northern Mexico) were the dominant
contributor to overall mercury deposition to Devils Lake at over 90%. Of the non-background fraction,
nearly 2/3 is estimated to derive from Canadian and other US sources, with only less than 3% of the total
flux originating from within Wisconsin. Of the within-state components, the primary source was coal-fired
utilities (1.2% of the total). Illinois and Minnesota sources contributed ca. 2 and 1% respectively. It is
important to note that these findings apply to the relatively isolated (from significant mercury emission
sources) Devils Lake location and just to the 1998 meteorological year modeled using the specific
inventory available to the study.
Describe water quality standards and numeric water quality target. Wisconsin's water quality
standards, fish consumption guidelines, and numeric targets for the pilot are discussed in Section 3. The
present quantitative endpoints for mercury in Wisconsin waters derive from its propensity for
bioaccumulation in aquatic food webs, presenting chronic health risks to humans and wildlife that eat large
amounts of fish. Wisconsin water quality criteria for protection of wildlife and human health are presented
in Table 2. In addition, WDNR has established guidelines for issuing mercury health advisories for fish
consumption to protect human health. For Devils Lake, the fish consumption advisory level is more
restrictive than the ambient water quality criteria and, for the purposes of this pilot study, the fish
consumption advisory level was selected for the target endpoint. It should be noted that this target may or
may not be the same target that would be used for an actual TMDL.
Identify reductions for point and nonpoint sources of pollution. For the purposes of this pilot exercise,
the pollutant reduction goal is the reduction of atmospheric mercury deposition necessary to achieve
mercury deposition needed to achieve a concentration in Devils Lake fish of less than 0.05 ppm for age 5
walleye. Walleye were selected as the target species because they are a top predator within the aquatic
food chain of Devils Lake, and the initial fish consumption advisory for the lake had also been issued for
7
-------
Final Report
March2006
walleye. Also considered were the reductions necessary to reach EPA fish tissue methylmercury criterion
of 0.3 ppm. Within the limitations of the assumptions of the air and water models, the pollutant reduction
goal for atmospheric deposition from all sources has been estimated.
Describe the linkage between water quality endpoints and pollutants of concern. This is described by
the modeling presented in Chapter 4 of this analysis, which suggests a strong direct relationship between
loads of atmospherically deposited mercury to Devils Lake and its watershed (Section 4.3.3) and mercury
concentrations in walleye. In addition, modeling suggests that, depending on the actual exchange dynamics
between the sediments and the water column, fish tissue concentrations of mercury may respond relatively
rapidly to changes in load (Section 4.3.3).
Identify margin of safety, seasonal variations, and critical conditions. A margin of safety could be
incorporated into the analysis as regards a human health reference dose for mercury, and in the atmospheric
and aquatic modeling. An analysis of various components of uncertainty in both REMSAD and D-MCM
predictions and the effects of model and parameter uncertainty on endpoint predictions are discussed
Section 5.
2.3 Watershed Description
Devils Lake is in the Driftless Area of southwestern Wisconsin, approximately 65 km northwest of
Madison (Figure 1). The lake is part of the Devils Lake State Park, which is one of the most heavily used
state parks in Wisconsin (Krohelski and Batten, 1995). The lake has a surface area of approximately 148
hectares (ha), and a drainage area of 839 ha (see Figure 2). The lake receives nearly half (47%) of its
hydrologic income from precipitation directly on the lake surface. The remaining 53% of the hydrologic
inputs derive from surface runoff from an intermittent stream that flows into the lake from the southwest.
The lake has no surface outflow, and water exits the lake via evaporation and seepage (37 and 63%
respectively). The lake is bounded on the eastern and western shores by steep talus slopes of boulder-sized
Precambrian quartzite, and on the northern and southwestern shores by gently sloping beach areas
underlain by glacial, fluvial, and lacustrine sediment (Krohelski and Batten, 1995).
8
-------
Final Report
March2006
y
~P 7"
Devil's
Lake,
WI
Figure 1. Site location map for Devils Lake, Wisconsin.
Because of its depth. Devils Lake stably stratifies each summer, and anoxia typically develops in the
hypolimnion in the late summer. Phosphorus is released from the sediments during anoxia, which is then
introduced into the euphotic zone during thermal destratification. These pulsed releases of phosphorus
contribute to algal blooms during late summer and early fall and, as a result, WDNR begun withdrawing
hypolimnetic waters in 2002 as a means to mitigate this problem. Phosphorus-enriched waters from the
hypolimnion are removed from early September until mid-October (the onset of turnover) as a means of
phosphorus removal from the sediments (the source of phosphorus enrichment in the hypolimnion during
anoxia) and ultimately reduce internal releases of phosphorus. During 2002 and 2003, 616,000 and
213,000 m' of water were withdrawn from the lake. In late 2003, WDNR began replacing the water
withdrawn from the lake with water from Babbling Brook, a nearby intermittent stream. If successful, the
onset of anoxia should be either delayed or eliminated after a time frame of treatment approximating 15
years or so (R. Lathrop, pers. comm.); this in turn may also positively affect the cycling of mercury in
Devils Lake by precluding the release of iron-bound mercury that likely occurs during anoxia (see Sections
4.3.3 and Appendix II, Section 7.1).
Devils Lake has an average depth of 9.3 in. and a maximum depth of ca. 15.2 m. The average (outflow-
based) hydraulic residence time of the lake is 8.1 years.
9
-------
Final Report
March2006
Figure 2. Map of watershed basin boundaries for Devils Lake, Wisconsin. From Krohelski and
Batten, 1995.
2.4 Existing Conditions and Summary of Monitoring Data
General limnological characteristics including a summary of relevant water chemistry for Devils Lake are
included in Table 1. The lake is considered eutrophic, experiencing blue-green algal blooms during late
summer and early fall (Krohelski and Batten, 1995). Total mercury concentrations average less than 0.6
ng/L, which is less than 1/2 the current WDNR standard for total Hg in surface water to protect wildlife.
Concentrations of mercury in walleye, a top predatory fish, average greater than 0.5 ppm, which is an order
of magnitude greater than the current WDNR advisory limit of 0.05 ppm.
Table 1. General Characteristics of Devils Lake, Wisconsin
Parameter
Value
Surface Area
1.47 km2
Area of Drainage basin
8.4 km2
10
-------
Final Report
March2006
Maximum Depth
15 m
Mean Depth
9.3 m
Productivity
Eutrophic
Flow pattern
Seepage outflow
Thermal stratification in summer
Yes
Anoxia in hypolimnion in summer
Yes
Dissolved organic carbon
Assumed low
Surface water pH
-7.7
Surface water sulfate
-190 |a,eq/L
Sedimentation rate in profundal zone:
- 1 mm/yr
Inorganic Hg(II) concentrations in surface waters:
Hg(II): 0.47 ng/L
unfiltered (mean from
1993-96 dataset, n= 43)
Methylmercury concentrations in surface waters
MeHg: 0.036 ng/L
unfiltered (mean from
1993-96 dataset, n= 8)
Top predator fish
Walleye
Mercury advisory
Yes
2.5 Identification of Pollutants Being Addressed and Why
Devils Lake falls under Wisconsin's general statewide fish consumption advisory to protect against
exposure to high concentrations of mercury, particularly in women of childbearing years, and children
under the age of 15. As a result, this pilot study addresses the trace element mercury in its various
environmental forms. The behavior of mercury in aquatic environments is highly complex with each of the
several chemical forms behaving differently. Methylmercury, formed from inorganic mercury by sulfate-
reducing bacteria in the sediments is the most biologically active form. Once formed, methylmercury is
readily taken up and retained by organisms and tends to increase in concentration in higher trophic levels
(i.e., biomagnifies). Although at ambient levels methylmercury does not appear to significantly affect
plants, invertebrates, or fishes, when biomagnified, as in the Devils Lake by over 1-million fold, it poses
the risk of chronic neurotoxicity to birds and mammals, including humans.
11
-------
Final Report
March2006
Devils Lake is essentially a seepage lake with a very small watershed that provides some surface runoff.
No known sources of mercury exist within the watershed, and the only inputs of mercury are from
atmospheric deposition. 72% of the mercury flux into Devils Lake itself derives from atmospheric
deposition directly to the lake surface; the remaining 28% derives from surface runoff (Appendix II).
12
-------
Final Report
March2006
3 NUMERIC TARGETS AND WATERSHED
INDICATORS
TMDLs are developed to meet applicable water quality standards. These may include numeric water
quality standards, narrative standards for the support of designated uses, and other associated indicators of
support of beneficial uses. This section summarizes the existing water quality standards and fish
consumption advisory in Wisconsin, and identifies a target for the purposes of the pilot study. Note that
this is a pilot study, and the selection of a target for this study is not intended to presuppose what the target
may be in an actual TMDL.
3.1 Numeric Water Quality Standards
Wisconsin currently has in place a set of numerical water quality standards for mercury based on the type
of exposure (e.g., chronic or acute) to aquatic life, or exposure to different types of other species (i.e.,
mammals and birds) as well as humans (Table 2). The values for acute and chronic exposure for aquatic
organisms are for total recoverable mercury in the water column and are 0.83 and 0.44 |.ig/L respectively.
The most restrictive water quality criterion for Devils Lake is that for the protection of wildlife, established
at 1.3 ng/L.
Table 2. Mercury water quality criteria for surface waters in Wisconsin. Derived from WAC NR
105. All concentrations expressed as (ig/L.
Toxicity Category Cold Water Warm water sportfish, Limited aquatic life
warm water forage
and limited forage fish
Acute* 0.83 0.83 0.83
Chronic* 0.44 0.44 0.44
Wildlife 0.0013 (all surface waters)
Human threshold (public water 0.0015 0.0015
supply)
Human threshold (non-public 0.0015 0.0015 336
water supply)
Total recoverable mercury
At the same time, mercury concentrations in aquatic biota are a function of both mercury loads and the
biogeochemistry of the aquatic ecosystem itself. There is thus often little relationship between ambient
water column concentrations and fish tissue concentrations observed across different lakes and streams.
Water column mercury levels may be below the ambient water quality criterion and yet still reach high
levels in fish. For example, in Devils Lake, ambient concentrations of mercury currently average 0.51 ng/L
13
-------
Final Report
March2006
(or below the most restrictive water quality standard in Wisconsin of 1.3 ng/L), while concentrations in
walleye average 0.69 ppm.5
The level of methylmercury in fish is now thought to be the more appropriate endpoint in TMDLs for the
protection of human health, rather than the ambient water column criterion in place in many states.
Methylmercury is the form of mercury which accumulates in fish and wildlife and the form of mercury
most toxic to humans and wildlife. As a result, EPA issued a recommended fish consumption advisory of
0.3 ppm to protect human health. Because consumption of contaminated fish and shellfish is the primary
route of human exposure to methylmercury, EPA has expressed this water quality criterion as a fish and
shellfish tissue value, rather than a water column value. EPA is developing guidance on how states would
implement the new criterion. Wisconsin has not yet adopted a methylmercury water quality criterion.
3.2 Fish Consumption Guidelines
Although Wisconsin has not adopted a water quality criterion for methylmercury, Wisconsin does have a
fish consumption advisory in place for due to high levels of mercury in fish. In 1998, Devil's Lake was
listed on the State of Wisconsin's 303(d) list of impaired waters because of a fish consumption advisory
relating to high concentrations of mercury in walleye between 10 and 30 inches. Wisconsin's fish
consumption advisory, 0.05 ppm, is based primarily on limiting exposure to high concentrations of mercury
via ingestion of fish to developing fetuses and children under the age of 15, as well as women of
childbearing years. Advice varies from unlimited consumption to "do not eat", and meal-frequency advice
is based on keeping ingestion below health protection values.
In 2001, Wisconsin adopted a statewide advisory for most inland waters regarding fish consumption and
mercury exposure. Nearly all of the 1,200 water bodies tested within Wisconsin contain some species with
fillet mercury concentrations in excess of 0.05 ppm. Statewide advice for how many meals of fish people
can safely eat to reduce risk of exposure to mercury is based on statewide levels of mercury found in the
different species of fish. There are separate recommendations for sensitive populations including pregnant
women and young children and for all other individuals. Wisconsin also has special advisories for specific
waters where fish mercury levels are higher than in waters under the general statewide advisory. See
http://dnr.wi.gov/fish/consumption/generaladvice/html for further information on Wisconsin's fish
consumption advisories.
3.3 Numeric Target for the Devils Lake Pilot Study
As stated earlier, Wisconsin has in place ambient water quality criterion for mercury, but has not
established a fish-tissued based criterion for methylmercury. For purposes of the pilot study, the EPA
recommended methylmercury fish tissue criterion of 0.3 ppm was used as the target. This is only for
purposes of illustration and does not indicate what the target would be in an actual TMDL. For the
purposes of the pilot study, the endpoint predictions for different mercury loading rates are also compared
to the Wisconsin fish consumption advisory limit of 0.05 ppm.
5 Mercury concentrations can be strongly associated with age of length of fish from a particular waterbody. Such a relationship should
be considered when examining changes in mercury concentrations. Some mercury TMDLs have used a standard-length fish to take
into account changes in mercury concentrations with age and length. For a summary of mercury concentrations by species in
Wisconsin, see http://dnr.wi.gov/fish/consumption/1990-2005MercurvSummarv.pdf.
14
-------
Final Report
March2006
Note that a state may also use a fish consumption advisory as the basis for determining whether a
waterbody is impaired and developing a TMDL. According to an EPA policy memo issued in 2002, states
may use a fish consumption advisory to demonstrate the non-attainment of water quality standards if it
meets all of the following conditions:
• The advisory must be based on fish tissue data;
• The data are from the specific waterbody in question; and
• The risk assessment parameters of the advisory or classification are cumulatively equal to or less
protective than those in the state water quality standards.
If the advisory satisfies all of these conditions, the waterbody must be listed as impaired under Section
303(d) and a TMDL developed for that waterbody.6
Table 3 illustrates how the conditions in Devil's Lake compare to various mercury-related water quality
endpoints.
Table 3. Comparison of Conditions to Water Quality Endpoints.
Existing Value
Parameter
(Mean and range)
Water Quality Endpoint
Comments
Total mercury in
0.51 (0.1-0.8)
1.3 ng/L
This value is designed to protect
ambient water
Wisconsin's wildlife from adverse
(ng/L)
effects from ingesting surface water
or aquatic organisms taken from
surface water of the state.
Methylmercury in
0.69 (0.28- 1.7)*
WDNR statewide fish
This value is basis for WDNR's
edible sport fish
consumption advisory:
general statewide fish consumption
tissue
0.05 ppm"
advisory.
(PPm)
USEPA criterion:
0.3 ppm
EPA guidance; not an exceedance of
water quality standards unless a
State adopts it.
Data courtesy of Wisconsin DNR. Samples collected between September 1983 and April 1993; n = 46.
"Wisconsin issues special advisories for waters where fish mercury levels are higher than in waters subject to the general statewide
advisory (see http://dnr.wi.gov/fish/consumption/FishAdv08MercurvList.pdf.')
6 Wisconsin's 2008 criteria for listing a water as impaired are, in general, lppm for gamefish and 0.21 for panfish species (which are
the criteria for special advisories). Previously, waters were included on Wisconsin's 303(d) list where there were mercury fish
consumption advisories for those specific waters. Currently, most waters in the state are subject to the general statewide advisory.
After 2002, waters were included on the 303(d) list when they are added to the fish consumption advisory, and, if the special advisory
no longer applies, the waters are no longer listed. See http://dnr.wi.gov/org/water/wm/wqs/303d/2008/2008methodology.htm.
15
-------
Final Report
March2006
4 AIR AND WATERSHED MODELING: SOURCE
ANALYSIS FOR LOADINGS, MERCURY MASS
BALANCE, LINKAGE OF STRESSORS TO
WATER QUALITY ENDPOINTS
This section summarizes the methods and results of the atmospheric modeling (Myers el al., 2006) and
aquatic mercury cycling (Harris et al., 2005) modeling efforts used in this study. Full details of the
modeling efforts can be found in the Technical Reports provided as Appendix I (atmospheric modeling)
and Appendix II (aquatic mercury cycling modeling) of this report.
There are four primary forms of mercury present in the environment: three of these are methylmercury
{MeHg}, divalent mercury salts and compounds (Hg(II)} and elemental mercury (Hg(0)}. Hg(II) is
explicitly defined here as all divalent inorganic mercury (other than particulate associated divalent
inorganic mercury). The fourth form is Hg(p) representing divalent inorganic mercury (Hg(II)} that is
associated with suspended particulate matter in water or air. REMSAD explicitly treats thee elemental,
divalent gaseous, and divalent particulate forms of mercury only, while MCM includes all four primary
forms. Provisions have been made in the MCM for some of the particulate Hg(II) on non-living solids in
water to exchange slowly, while the remainder is assumed to exchange rapidly enough to assume
instantaneous partitioning.
4.1 Technical Approach
Two modeling components were employed for this pilot study. The first modeling approach was to
simulate the atmospheric transport, transformation, and deposition of mercury using the Regional Modeling
System for Aerosols and Deposition (REMSAD) model (Myers et al., 2006; Appendix I). REMSAD is a
three-dimensional grid model that was developed by Systems Applications International (SAI) for EPA
under separate contract (SAI, 2002). The second model, the Mercury Cycling Model (MCM), was used to
simulate the fate and transport of the mercury within the aquatic environment (Harris et al., 2005;
Appendix II). The MCM provides estimates of the rates of mercury methylation and bioaccumulation
through the food chain to a top-level predatory fish (e.g., walleye). Detailed descriptions of these models
and uncertainties are provided in the referenced technical support documents for this analysis.
4.2 Atmospheric Transport and Deposition Modeling
4.2.1 Overview of the Air Deposition Model, REMSAD
The Regional Modeling System for Aerosols and Deposition (REMSAD) was used to provide estimates of
wet and dry mercury deposition over Devil's Lake, Wisconsin, for this analysis. REMSAD was developed
by SAI with EPA support beginning in 1994 and has undergone continuous development since then (SAI,
2005). REMSAD is a three-dimensional grid model designed to simulate the atmospheric transport,
transformation, and deposition of criteria and toxic pollutants. To date, REMSAD has been configured to
address mercury, cadmium, dioxin, atrazine, and polycyclic organic matter (POM), and lead. In addition, it
can be used to provide deposition analysis of nitrogen and sulfur species and to examine particulate matter
(PM) air quality issues. The mercury modeling capabilities of REMSAD, as well as the algorithms
16
-------
Final Report
March2006
pertaining to the other pollutants, have undergone external peer review (Seigneur et al., 1999)
independently of and prior to initiation of the TMDL Pilot Project. The comments and suggestions raised in
that peer review pertaining to deposition modeling of mercury have been incorporated into the version of
REMSAD used in this application, Version 8. The REMSAD model is publicly available.
In this current document, which is a revision of a report submitted to EPA in July 2003, information is
provided for use in evaluating approaches for developing Total Maximum Daily Loads (TMDLs) involving
mercury from air deposition.
To help the reader in interpreting the results presented here, a brief timeline of developments that led to
this report is provided.
1. October 1999 - EPA sponsored peer review of the REMSAD modeling system is completed.
2. October 2000 - As an aid in conducting source attribution analyses, SAI develops and implements
mercury tagging in the REMSAD system.
3. July 2001 - Based on the above peer review, SAI implements updates to the REMSAD modeling
system.
4. October 2002 - SAI completes 48-state mercury modeling (tagging simulations of mercury emissions
from each of the lower 48 states) using the 1998 annual simulation as the basis.
5. July 2003 - SAI submits report for the WI TMDL Pilot Project, using the 1998 annual simulation as the
basis.
6. March 2004 - External peer review of the WI TMDL Pilot Project report is conducted.
7. January 2006 - Additional information is submitted in response to the peer review of the TMDL Pilot
Project report, using the 2001 annual simulation as the basis.
References to the REMSAD peer review mean the review of the REMSAD modeling system conducted in
1999 under a separate effort. References to the review of the pilot project refer to the review of the WI
TMDL Pilot Project report conducted in 2004. In Section 4 and the remainder of this integrated report,
references to the modeling will be to the revised modeling, conducted after and in response to the March
2004 review of the pilot project, unless otherwise noted. Technical support documents on both the original
modeling and the revised modeling are available on request. Major differences between the original 1998
simulation and the revised 2001 simulation are presented in Table 4.
Table 4. Major Differences between Original 1998 REMSAD Simulation and the Revised 2001
REMSAD Simulation
Feature
1998 Simulation
2001 Simulation
Emissions
Meteorology
1996 NTI w/ updates for utility testing and
certain other source updates.
1998 meteorological files based on RUC
forecast/analysis files from NCAR
2001 CAMR inventory
2001 meteorological files based on MM5
modeling by EPA OAQPS (as used for
CAIR and CAMR)
Boundary
concentrations
Static estimates based on literature and
consultation with experts
Spatially and temporally varying estimates
based on outputs from GEOS-CHEM
global mercury modeling, as developed
for CAMR modeling
17
-------
Final Report
March2006
Feature
1998 Simulation
2001 Simulation
Gas phase
chemical
mechanism
microCB
Full CB-V
Mercury
chemistry
Based on Lin and Pehkonnen (1999) As in 1998, but with updates based on
recent literature (In particular, gas phase
reaction of elemental mercury with OH
radical added.)
Tagging setup
6 emisisons categories in Wisconsin plus Several individual states,
several individual states
A flow diagram of the REMSAD sub-routines is provided in Figure 3. REMSAD was based upon and
derived from the variable-grid Urban Airshed Model (UAM-V) modeling platform, also developed by SAI
(SAI, 1999). One key area in which the original REMSAD model differed from UAM-V was the manner in
which it calculated estimates of atmospheric oxidants such as ozone. Because many pollutants of interest,
including mercury, can react with oxidants in such ways that affect their transport and deposition, it is
necessary to have knowledge of oxidant concentrations. UAM-V uses the Carbon Bond mechanism 5 as the
foundation for the atmospheric chemistry calculations (an updated version of the well-documented Carbon-
Bond 4 (CB-IV) mechanism (Gery et al., 1989)). CB-IV is used by many of today's airshed models
including the EPA's Community Multiscale Air Quality (CMAQ) model (Byun et al., 1999). In response to
the review of the pilot project, REMSAD version 8 was modified to include the full CB-V mechanism. For
its atmospheric chemistry calculations, the revised modeling used CB-V rather than the micro-Carbon
Bond (|iCB) mechanism (SAI, 2005) that was used in the original modeling.
18
-------
Final Report
March2006
Figure 3. Flow of Control through Key Subroutines in a REMSAD Mercury Simulation.
19
-------
Final Report
March2006
One of the primary suggestions from the initial 1999 REMSAD peer review was that the treatment of
mercury chemistry in the model should be expanded to reflect the significant findings reported in recently
published literature. Prior to conducting the original 1998 modeling for the WI TMDL pilot project, a
mechanism was incorporated into REMSAD based primarily on Lin and Pehkonen (1999). Since that time,
further updates to the mechanism have been made based on recently published information. The mercury
reactions and associated reaction rates currently implemented in REMSAD Version 8, the version used in
the revised 2001 modeling, are presented in Table 5. The species HGO represents elemental mercury, HG2
represents divalent gas mercury, and HGP represents divalent particulate mercury.
Table 5. Mercury Chemical Mechanism in REMSAD, Version 8.
Reaction
Rate (unit)
For HGO
Gas phase
HGO + 03 -> !/2 HGP + !/2 HG2
3.0e-20 (cm3molecule"1s"1)
hgo + no3->hgp
4.0e-15 (cm3molecule"1s"1Xcurrently disabled)
HGO + H202 -> HG2
8.5e-19 (cm3molecule"1s"1)
HGO + OH —> !/2 HGP + '/2 HG2
7.7e-14 (cm3molecule"1s"1) (lower bound of range in Pal &
Ariya, 2004; also within range in Sommar et al., 2001)
Aqueous phase
HGO + 03 -> HG2
4.7e+7(M"1s"1)
HGO + OH —> HG2
2.0e+9 (M-'s"1)
HGO + Claq -> HG2
(See eq. 8 in Lin and Pehkonen, 1999)
HgS03 -> HGO
T e(31971T - 12595)/T s"1 (Van Loon et al., 2000)
For HG2
Aqueous phase
HG2 + H02 -> HGO
1.7e+4(M"1s"1)
HG2 + S032" <-» HgS03
5.e+12 (M"1)
HgSOj + S032" ^ Hg(S03)22"
2.5e+ll (M"1)
Hg2 + OH" <-» Hg(OH)+
4.27e+10 (M"1)
Hg2 + 2 OH" <-» Hg(OH)2
1,74e+22 (M"1)
Hg2 + OH" + CI" ^ HgOHCl
1.78e+18(M"2)
Hg2 + CI" <-» HgCl+
2.0e+7 (M"1)
Hg2 + 2 CI" ^ HgCl2
l.e+14 (M"2)
20
-------
Final Report
March2006
Hg2 H
H3Cr<->HgCl3-
l.e+15 (M"3)
Hg2 H
H4Cl-^HgCl42-
3.98e+15 (M"4)
Source: Lin and Pehkonen, 1999, except as noted.
The treatment of mercury in REMSAD also reflects the belief that some divalent mercury is adsorbed to
soot particles (e.g., see Seigneur et al., 1998). REMSAD uses a parameterized formula to simulate this
phenomenon whereby a certain portion of dissolved divalent mercury in aqueous phase is assumed to be
adsorbed to soot particles (as indicated by the modeled concentration of primary elemental carbon). This is
further discussed in Section 2 of the report.
4.2.2 Modeling Domain
Because the regional and long-range transport aspects of mercury deposition to Devil's Lake was of interest
to the project team, the air quality modeling domain includes the entire coterminous U.S. and extends
approximately 400 miles north into Canada, approximately 300 miles south into Mexico, and roughly 200
miles east and west over the Atlantic and Pacific Oceans. In the Devil's Lake modeling application,
REMSAD was operated with a combination of 36x36 km coarse grids over most of the domain with 12x12
km fine grids over selected areas of interest such as Wisconsin.
4.2.3 Emission Inputs and Background Assumptions
The emissions inventory that was used in this modeling was the 2001 inventory utilized by EPA in the
Clean Air Mercury Rule (CAMR) modeling. Further information on the mercury emissions inventory can
be found in "Emissions Inventory and Emissions Processing for the Clean Air Mercury Rule (CAMR)"
(EPA, 2005). Mercury emission data for the year 2001 for area and point sources were included in the
CAMR inventory. U.S. emissions therefore represent the same calendar year (2001) as the meteorological
data described in Section 3A of Appendix I. Table 6 summarizes the U.S., Canada, and, separately, WI
mercury emissions by key source categories. The associated average speciation of the emissions among
elemental, divalent gaseous, and divalent particle-bound forms is presented in Table 7.
Table 6. Summary of 2001 Mercury Emissions for Revised REMSAD Modeling (emissions in
tons/year).
Location
Category Area
Point
Wisconsin
Coal Utilities 0.00
1.13
Wisconsin
Coal Industrial Boilers 0.00
0.05
Wisconsin
Chlor-alkali Plants 0.00
0.54
Wisconsin
Medical Waste 0.00
8.4E-06
Wisconsin
Municipal Waste Incinerators 0.00
0.07
Wisconsin
All other WI sources 0.10
0.30
All other states
All 6.43
101.7
All U.S.
All 6.53
103.79
21
-------
Final Report
March2006
Location Category Area Point
Canada All 2.29 5.78
Table 7. Average Speciation of Mercury Emissions.
Location
Category
HG0
(% Total)
HG2
(% Total)
HGP
(% Total)
Wisconsin
Coal Utilities*
68
32
1
Wisconsin
Coal Industrial Boilers
52
29
19
Wisconsin
Chlor-alkali Plants
95
5
0
Wisconsin
Medical Waste
5
75
20
Wisconsin
Municipal Waste Incinerators
22
58
20
Wisconsin
All other Wl sources
65
20
15
All other states All
60
31
9
All U.S.
All
60
31
9
Canada
All
55
35
10
* Average speciation profiles were not used in modeling coal utility emissions; rather, unique speciation profiles from the CAMR
mercury emission inventory (EPA 2005c) based upon coal type, boiler type, and air pollution controls were used.
Estimates of emissions for key pollutants that play a direct or indirect role in the atmospheric reactions
involving mercury, e.g., ozone precursors, sulfur dioxide, and PM (collectively referred to here as the "PM
inventory," were the same ones used to support the Agency's recent modeling for the Clean Air Interstate
Rule (CAIR) and reflect the 2001 calendar year (EPA, 2005b). The emissions data for the 2001 base-year
(for mobile, area, point-, and biogenic sources) are described in detail by E. H. Pechan (2000). The annual
total emissions of criteria pollutants are summarized in Table 8.
Table 8. Total Annual U.S. Pollutant Emissions (106 ton/year) in PM inventory: 2001 Base Year.
Species
Point
Area
Mobile
Biogenic
Total
NOx
7.89
1.61
12.13
1.05
22.68
VOC
1.63
7.74
7.49
41.84
58.71
S02
13.87
1.29
0.71
15.86
PM2.5
1.30
3.22
0.47
4.99
PMC
0.36
8.43
0.07
8.87
22
-------
Final Report
March2006
NH3 0.09 3.32 0.28 3.69
CO 4.42 11.76 83.95 5.82 105.96
* Total emissions of on-road and non-road mobile sources
Because the air mass entering the US already contains mercury from numerous natural and man-made
sources from around the world, it is necessary to approximate the presence of this mercury in the modeling
exercise. For the revised (2001) simulation the results of a global model simulation of mercury were
available to set the boundary concentrations. The global model GEOS-CHEM (Yantosca, 2004) was run by
Harvard University and the outputs processed as boundary concentration files for the EPA OAQPS CAMR
simulations (EPA, 2005c). The GEOS-CHEM-based boundary files was chosen for the revised modeling.
The boundary concentrations derived from the GEOS-CHEM model output vary horizontally along each edge
and vertically. The boundary concentrations have been averaged in time so that they vary monthly. The
average concentrations of HG0 range from about 1 X 10"7 ppmto 1.4 X 10"7 ppm. The average concentrations
of HG0 in February are somewhat higher with values close to 2 X 10"7 ppm near the surface. The HG2
concentrations exhibit a relatively wide range of concentration depending on which boundary is considered.
The North and West boundaries have average boundary concentrations ranging from about 0.8 X 10~8 to 1.5 X
10~8 ppm in July. The average February HG2 concentrations are about half the July values. The average HGP
boundary concentrations in the revised modeling are typically less than 2 X 10"10 ppm.
4.2.4 Meteorological Inputs for Revised (2001) Modeling
For the revised modeling using REMSAD, the project used existing databases at 36-km resolution for
meteorology that had been developed by EPA for use in their evaluation of emissions rules. The
meteorological data files were prepared by EPA from processed MM5 outputs. Data files for the calendar
year 2001 were available. The last several days of December 2000 were also available for a model "spin-
up" period. These 2001 meteorological files were used in the EPA's evaluation of the Clean Air Interstate
Rule (CAIR) and the Clean Air Mercuiy Rule (CAMR) (EPA, 2005b; EPA, 2005c).
Because rainfall is a key parameter in modeling mercury deposition, Figure 4 and Figure 5 were assembled
to provide insight into how closely the MM5-based data files matched observations from the National Acid
Deposition Program (NADP) network. Figure 4 shows the total annual rainfall derived from NADP
monitors. Figure 5 shows the total annual rainfall as represented in the REMSAD input files. As can be
seen, the pattern of projected precipitation from the MM5 derived files compares favorably with
measurements. The observed precipitation plot was obtained from the National Atmospheric Deposition
Program (NADP) web site (http://nadp.sws.uius.edu) (NADP, 2005).
23
-------
Final Report
March2006
Total precipitation, 2001
1« 28
Sites not pictured:
AK01 35 cm
AK03 24 cm
H199 374 cm
VI01 105 cm
i y
••r J
71
t.95
r A i
, TUO 100
82 73*
^ 99 » -
69 $6 • gj
• " n •
" „• fc B . £.
*5)
W, «
--V"M
3 97
3 -»
«
83 30
74^75
»• *» .
• " *
io® 119 t«'102 •
• i -Vm 9*
36 *118 *99
?T
*90
it »Hf « & . ;
\L ^ t7 • 88 Precipitation
J34 3<
* 88 101
(cm)
SB
100,
__jH
V13I
" 119 *
132
National Atmospheric Deposition Program/National Trends Network
http://nadp.sws.uiuc.edu
£ 20
20-40
40-60
60-80
80-100
100-120
120-140
140 -160
160-180
180-200
] >200
Figure 4. Observed Annual Rainfall Totals (cm) for 2001.
LEVEL 1 RAIN (cm) + MAXIMUM = 544.4 cm (145,5)
Time: 0 Apr 11. 2010-0 Apr 10. 2020 - MINIMUM = 0.4 cm (30,32)
Km
-2?30. -2028. - 1296. -576. 144. 864. 15B4. 2304.
v % *b «b %% % % % % •*
^ O O -o 0 O Q O O Q $0
° *0 so % % % \ % % % °
O O O O o -Q Q Q G
Total Precipitation for 2001 MMb REM SAD input file (cm)
Figure 5. MM5-Derived Annual Rainfall Totals (cm) for 2001.
-------
Final Report
March2006
4.2.5 Analysis of Year-to-Year Variation
An analysis of year-to-year variation due to meteorological conditions was conducted based on the original
(1998) REMSAD modeling. Meteorological data for this REMSAD application were obtained from Rapid
Update Cycle (RUC) analysis and forecast system operated by the National Oceanic and Atmospheric
Administration (NOAA) for the calendar year 1998. (See http://maps.fsl.noaa.gov/, maintained by NOAA
Forecast Systems Laboratory (FSL), RUC development group for more information.) This database was
chosen after discussions among the research and applications community. At the time the project was
begun in 1999 the RUC was unique among models using data assimilation in its ability to cover the entire
modeling domain for an entire year with the required data and resolution. Further discussion of the RUC
meteorology can be found in Section 3B of this report. In the following discussion of year-to-year
variability, the reader should bear in mind that the analysis was based on the 1998 REMSAD modeling
using the RUC meteorological files rather than the revised 2001 modeling using MM5 files.
Regardless of how well REMSAD replicates wet and dry mercury deposition for a given year, a key
question relevant to the TMDL Pilot application is, "How much would one expect deposition to vary at
Devil's Lake from year to year given a typical variation in meteorological parameters affecting
deposition?" The issue of the representativeness of the original baseline year, 1998, was examined through
the use of the Classification and Regression Tree (CART) statistical software (Breiman, et al., 1984;
Steinberg, et al., 1997) for the 10-year period of 1992 to 2001. This statistical classification method is used
to provide an estimate of the variations in deposition (due to meteorological variations) over that 10-year
period.
CART was used to statistically identify the most important meteorological parameters and conditions
associated with certain ranges of daily wet and, separately, dry mercury deposition from the 1998
REMSAD outputs. Key surface and upper-air meteorological parameters are listed in Table 9. From this
analysis, a set of approximately 30 bins with wet deposition values and another set of approximately 30
bins with dry deposition values were defined. Each of these bins can be characterized by the associated
meteorological parameters from the 1998 calendar year RUC database. A second frequency-weighted
dataset of these same meteorological parameters was then constructed for each year from 1992 to 2001
using data collected from National Weather Service (NWS) stations near Devil's Lake (Madison, WI NWS
for surface meteorological variables and Green Bay, WI NWS for upper-air variables.) The average
deposition value (from the 1998 REMSAD analysis) associated with each CART bin was used to estimate
deposition for each day and construct pseudo deposition totals for each year in the ten-year time period.
The CART-based deposition values were supplied as inputs to the water modeling component of the
Devil's Lake TMDL Pilot Project.
Table 9. Key Classification Parameters for Wet and Dry Deposition (not necessarily in order of
importance).
Wet Deposition
Dry Deposition
Total daily rainfall
Maximum surface temperature
Average temperature at 850 mb
Average surface wind direction
Average moisture at 700 mb
Wind speed for the morning sounding at 850
mb
25
-------
Final Report
March2006
Wind direction at 850 mb on the previous! Temperature lapse rate between the surface
afternoon and 900 mb for the morning sounding
Total daily rainfall
The estimated total annual wet and dry deposition based on the CART analysis is presented in Table 10. In
reviewing the results of this analysis, discussed in more detail in Section 6 of the report, the reader is reminded
that the CART-based analysis examines the basic question, "How much variation in year-to-year deposition of
mercury would one expect due solely to changes in meteorological conditions?" In answering this question,
emissions were assumed to be constant throughout the 1992 to 2001 time period. While it is believed that the
mercury inventory used in the original 1998 modeling is representative of actual emissions in the 1996-1999
timeframe, the emissions estimates may or may not be applicable to years considerably earlier or later than
this.
The CART analysis provides deposition estimates for a ten-year period from 1992 to 2001. Table 10 presents
the results from the REMSAD simulated wet and dry mercury deposition for 1998 to Devil's Lake together
with the CART-based estimates of deposition for each year from 1992 to 2001. The REMSAD results for
1998 indicate wet deposition of 7.1 g/km2 and dry of 12.0 g/km2, while the CART analysis shows, on average,
a nearly equal contribution of wet and dry deposition to Devil's Lake for the ten-year time period, 9.3 g/km2
and 10.3 g/km2 for wet and dry deposition, respectively. Since the 1998 simulation results differ somewhat
from the results of the revised modeling for 2001, the reader should use the range values in Table 10 as an
indication of the potential variation in deposition over the 10-year period, but should consider the values
found in the revised modeling (see the next section) for estimating the magnitude of deposition.
Table 10. Deposition Estimates (g/km2) for Devil's Lake.
Year
Wet Deposition
Dry Deposition
REMSAD
1998
7.1
12.0
CART-based Analysis
1992
5.1
9.2
1993
11.1
10.0
1994
10.8
11.2
1995
9.7
11.1
1996
6.4
9.3
1997
8.5
9.8
1998
11.1
10.4
1999
11.1
10.7
2000
8.5
10.4
2001
9.3
10.6
10-year CART Average*
1992-
9.3
10.3
2001
26
-------
Final Report
March2006
These ten-year averages were supplied as input to the water modeling
portion of the project.
4.2.6 Results and Model Performance
The REMSAD simulation results provide estimates of the total wet and dry deposition of mercury for 2001.
The REMSAD results for 2001 indicate wet deposition of 5.9 g/km2 and dry of 4.6 g/km2.
In order to evaluate the performance of REMSAD, wet deposition measurements of mercury from the
Mercury Deposition Network (MDN) were gathered from the 53 stations with data covering 2001. Figure
6 is a scatter plot showing the REMSAD modeled wet deposition values versus the MDN measured values
at these locations. On the basis of this comparison, the REMSAD wet deposition performance is quite good
with a correlation (R2) of 0.66 and a normalized bias of 23%. Figure 7 shows the results of a similar
analysis of MDN and REMSAD data in which measurements from states with more than one MDN station
are combined and plotted with data from states with a single MDN station in order to examine state-and
regional-specific performance. This plot indicates that REMSAD performs well in modeling wet deposition
of mercury in all parts of the U.S. that had MDN data in 2001. A statistical summary of performance is
provided in Table 11. (For a discussion of performance for the original 1998 simulation, see Section 5B of
Appendix I.)
£
£
3
"O
4)
+*
3
£
'(/>
Observed (gm/km2)
Figure 6. Simulated and Observed Annual Total Wet Deposition at MDN Monitors (gm/km2).
27
-------
Final Report
March2006
Observed (gm/km2)
Figure 7. Simulated and Observed State Averages of Annual Total Wet Deposition at MDN
Monitors (gm/km2).
Table 11. Statistical Evaluation of Simulated Total Annual Wet Deposition of Mercury against
Observed Data from MDN Sites for the Year 2001.
All sites
State
Averages
Number of data pairs
53
24
R2
0.66
0.76
Normalized % error
35
31
Normalized % bias
23
22
Average observed value (gm/km2)
9.3
8.5
Average simulated value (gm/km2)
11.1
10.3
Because there is no comparable network of dry mercury deposition monitors, it is not possible to
adequately evaluate REMSAD performance for dry deposition by comparisons to observations. It should be
noted, however, that the general approach to dry deposition in REMSAD is based upon the long-established
methodology applied in the Regional Acid Deposition Model (RADM) (Wesley, 1989). Under this
scheme, the dry deposition flux of a pollutant is proportional to its concentration in the lowest model layer
and its deposition velocity. Deposition velocities differ for gaseous and particulate species and depend in
part on such parameters as atmospheric stability, air density and viscosity, the molecular weight of the
depositing species, and the type of surface in question (e.g., water, forest, etc.). A more thorough
discussion of dry deposition in REMSAD and similar models can be found in the REMSAD User's Manual
28
-------
Final Report
March2006
(SAI, 2005) and in Wesley (1989). (The REMSAD User's Manual is publicly available; see
www. remsad. com.)
It should also be noted that recently published research suggests that many bodies of water are super-
saturated with elemental mercury (Schroeder and Munthe, 1998). The presence of mercury in such
quantities would result in no net deposition of dry-deposited mercury. Due to such evidence, it was
assumed, like other modelers have done (Pai et al., 1997), that there is no net dry deposition of elemental
mercury.
4.2.7 Source Attribution—Tagging
As discussed in the last section, evidence from the sensitivity analysis suggests that background
concentrations of mercury can play an important role in deposition at locations relatively distant from or
upwind from local emission sources of mercury. In order to further investigate this premise and to provide
the project team with insight into the sources most responsible for deposition to Devil's Lake, as well as to
other areas of Wisconsin affected by mercury deposition, a "tagging" scheme was developed and
implemented into REMSAD as part of the TMDL Pilot Project. Note that this new feature, tagging, had
not been incorporated into REMSAD prior to the October 1999 REMSAD peer review mentioned earlier,
but was used in the 1998 modeling results that underwent peer review as part of the WI TMDL pilot
project.
Traditionally, a "zeroing-out" methodology has been employed to investigate the relative contributions of
specified sources to areas with air quality and/or deposition impairments. Under the zero-out methodology,
if one wanted to learn the relative contribution of Canada's mercury emissions to deposition at a certain
location, for example, a set of two model runs would be required. With the first, the inventory for the
entire domain would be used. A second run would use the same inventory except that Canada's emissions
would be removed. The difference in deposition between the two runs would reflect Canada's contribution.
Because the project team desired analysis of the relative contributions to mercury deposition from each of
six different source categories within Wisconsin as well as from each of several nearby states in the
Midwest region, Canada, and from background, a method that reduced the computational and time
requirements compared to the zeroing-out technique was desired.
With tagging, emissions from specified sources, source categories, or geographic regions are labeled as
unique pollutants in the model's internal accounting while maintaining the appropriate chemistry and
physics treatment for each species. For example, an inventory that has been prepared with tags could have a
unique identifier for elemental mercury emissions from Canada. Those emissions would be subjected to the
same algorithms affecting transport, transformation, and deposition as any other elemental mercury
emissions in the inventory. However, from an accounting point of view, Canada emissions would be
considered to be essentially a different pollutant under tagging, facilitating analysis of its relative
contribution to deposition to any grid square in the domain. Tagging can thus facilitate developing
attribution information from numerous sources simultaneously in a greatly reduced number of model runs,
resulting in savings of time and resources. A run with 50 tags, for example one for each of the lower 48-
states plus one each for Canada and background, would require approximately a month's time to model an
entire year with REMSAD on a current generation high-end PC. Running 50 separate simulations would
take a substantially greater amount of time.
29
-------
Final Report
March2006
During the development of the tagging procedure, individual tags were defined such that the sum of the
emissions of all tags was equal to the total mercury emissions. Including the boundary and initial
concentration tag, the sum of all tag contributions (to concentration or deposition) theoretically would be
the same as the result of modeling all mercury emissions and boundary and initial conditions combined. A
comparison of results of the sum of the tags and simulations including all emissions showed discrepancies
in the mercury totals that were typically only 1 - 2%. This consistency serves as a source of verification
that the tagging method is functioning correctly.
Figure 8 shows simulated contributions to deposition at Devil's Lake from Wisconsin, Canada, and other
individual states in the region for 2001. As suggested with the sensitivity analysis, background was found
to be the dominant contributor to overall mercury deposition to Devil's Lake at 76%. (At other locations in
Wisconsin that are closer to sources, the importance of background was found to be nearer 60%.) Several
states contribute to the remainder of the deposition at Devil's Lake including Wisconsin and Illinois and to
a lesser extent Minnesota, Iowa, and Michigan. States that are not immediate neighbors of Wisconsin
("Other states") contribute 6%, collectively, to the deposition of mercury at Devil's Lake. Emissions from
Canada contribute another 1%.
Figure 9, based on the original 1998 modeling, presents estimates of contributions to deposition from
several source categories within Wisconsin. (Estimates for individual source categories are not available
from the 2001 simulation.) Of the non-background components at Devil's Lake, Wisconsin's utilities
contributed more to mercury deposition at that location than any other Wisconsin source—15% of the non-
background fraction. The estimate of the contribution of Wisconsin emissions to deposition at Devil's Lake
is 35% of the non-background component based on the 1998 simulation. The revised modeling estimates a
somewhat lower contribution (about 25%, disregarding re-emission), but it was expected that the relative
contributions among source categories would remain about the same. It is important to note that these
findings apply to the relatively isolated (from significant mercury emission sources) Devil's Lake location
and just to the 1998 meteorological year modeled using the specific inventory available to the project. As
noted in the sensitivity analysis discussion, results can be affected by uncertainties in these and other
factors. Additional findings from the tagging analysis are discussed in Section 7 of Appendix I.
30
-------
Final Report
March2006
Total deposition =10.5 g/km2
~ Wisconsin
¦ Minnesota
~ Illinois
~ Iowa
¦ Michigan
~ Other states
¦ Canada
~ Background
¦ Re-emission
Figure 8. Simulated Contributions to Annual Total Deposition (Wet and Dry) of Mercury at Devil's
Lake, Based on the Revised 2001 Simulation. {Wet deposition contributes 5.9 g/km'; dry deposition
contributes -1.6 g/km2.) Simulated concentrations and deposition are affected by uncertainties in
model input data and model formulation. See uncertainties discussion in Section 5.
~ Wl coal utilities
¦Wl Industrial boiler
~Wl Chlor-alkali plant
~Wl Medical waste
¦Wl Municipal waste incinerator
~ Other WI
¦ Other states and Canada
Figure 9. Contributions to Non-Background Mercury Deposition at the Devil's Lake Site, Based on
the Original 1998 Modeling. Simulated concentrations and deposition are affected by uncertainties
in model input data and model formulation. See uncertainties discussion in Section 5.
4.3 Aquatic Mercury Cycling Modeling
4.3.1 Model Description
The Dynamic Mercury Cycling Model (D-MCM) is a Windows 95/98/NT™ based simulation model for
personal computers (EPRI, 2002). It predicts the cycling and fate of the major forms of mercury in lakes,
including methylmercury, inorganic Hg(II), and elemental mercury. D-MCM is a time-dependent
mechanistic model, designed to consider the most important physical, chemical and biological factors
affecting fish mercury concentrations in lakes, including mercury loading rates. It can be used to develop
31
-------
Final Report
March2006
and test hypotheses, scope field studies, improve the understanding of cause/effect relationships, examine
responses to changes in mercury loading, and to help design and evaluate mitigation options.
The development of D-MCM has been funded by EPRI and the Wisconsin Department of Natural
Resources. It is an extension of previous mercury cycling models developed by Tetra Tech, including the
original Macintosh-based MCM models developed during the EPRI-sponsored Mercury in Temperate
Lakes Project in Wisconsin (Hudson et al., 1994). and the subsequent steady state Regional Mercury
Cycling Model (Tetra Tech Inc.. 1996).
Using a mass balance approach. D-MCM predicts time-dependent concentrations for three forms of
mercury in water and sediments (dissolved and particulate phases) and a simplified food web (see Figure
10 and Table 12).
Figure 10. Mercury cycling in D-MCM.
Mercury Forms
There are three primary mercury forms in the D-MCM: methylmercury, Hg(II) and elemental mercury.
Hg(II) is defined here as all mercury which is neither methylmercury nor elemental mercury. Provisions
have also been made for some of the particulate Hg(II) on non-living solids to exchange slowly if desired,
while the remainder is assumed to exchange rapidly enough to assume instantaneous equilibrium
partitioning.
Model Compartments
w
Mercury Cycling in D-MCM "%J?"
Wet and dry
Deposition
ttU. \ Volatilization
32
-------
Final Report
March2006
Model compartments include the water column, sediments and a food web that includes three fish
populations. Mercury concentrations in the atmosphere are input as boundary conditions to calculate fluxes
across the air/water interface (gaseous, wet deposition, dry deposition). Similarly, watershed/upstream
loadings of inorganic Hg(II) and methylmercury are input directly as time-series data, not modeled. The
user provides inputs for flow rates (surface and groundwater) and associated mercury concentrations, which
are combined to determine the watershed mercury loads.
Table 12. Compartments and mercury forms in D-MCM.
Compartment
Mercury Form
MeHg
Hg(H)
Elemental Hg
Solid HgS
Water Column (abiotic)
Dissolved
•
•
•
Non-living suspended particles*
•
•
•
Sediments
Sediment Porewater
•
•
•
Sediment Solids*
•
•
•
Food Web
Phytoplankton
•
•
Zooplankton
•
•
Benthos
•
•
Non-Piscivore Fish Cohorts (up to 20)
•
Omnivore Fish Cohorts (up to 20)
•
Piscivore Fish Cohorts (up to 20)
•
Includes slowly and rapidly exchanging components for Hg(II)
33
-------
Final Report
March2006
The food web consists of six trophic levels (phytoplankton, zooplankton, benthos, piscivore fish, omnivore
fish and non-piscivore fish). Specific fish species can be selected. Fish mercury concentrations tend to
increase with age, and are therefore followed in each year class. Methylmercury fluxes for individual fish
were coupled to fish bioenergetics equations from Hewett and Johnson (1992) as described by Harris and
Bodaly (1998). Fluxes were then scaled up to represent year classes and entire populations.
D-MCM Processes
A detailed description of model processes and equations in D-MCM is provided in the user's guide (EPRI,
2002). An overview of the major processes involved in the mercury cycle in lakes are shown in Figure 10.
These processes include surface inflows and outflows, vertical groundwater flow, instantaneous
methylmercury partitioning between abiotic solids and dissolved complexes, instantaneous and slower
adsorption/desorption kinetics for Hg(II) on abiotic solids, particulate settling, resuspension and burial,
atmospheric deposition, air/water gaseous exchange, in-situ transformations (methylation, demethylation,
MeHg photodegradation, Hg(II) photoreduction), mercury kinetics in plankton, and methylmercury fluxes
in fish populations (uptake via food and water, excretion, egestion, mortality, fishing).
The predictive capability of D-MCM is evolving but is currently constrained by a number of deficiencies in
the state-of-the-art understanding of the processes that govern the biogeochemical cycling of mercury.
These knowledge gaps include:
• The true rates and governing factors for methylation and demethylation, including whether all Hg(II)
present is available for methylation (i.e., is methylation driven largely by newly deposited Hg, or is the
entire Hg(II) pool of both new and previously deposited available). Moreover, the precise locus of
where these two processes occur is not well understood.
• Rates of inorganic Hg(II) reduction, which are likely linked to DOC concentrations in the water
column and light penetration.
• Site factors affecting methylmercury uptake at the base of the food web.
Furthermore, the relationship between mercury loading and fish mercury concentrations is not adequately
known when considering loading changes that might occur due to reductions in atmospheric inputs.
4.3.2 Modeling Approach
Overall, the aquatic modeling approach was designed to estimate the magnitude and timing of the response
of fish mercury concentrations to changes in atmospheric Hg(II) deposition, while considering the model
sensitivity to different inputs, uncertainty regarding true Hg(II) loading rates and other model inputs, and
year-to-year variations in mercury deposition in the absence of a long-term record of mercury deposition at
the site.
Aquatic modeling for the Devils Lake pilot mercury TMDL included the following key components:
• Selection of a lake identified by the State of Wisconsin 303(d) process as impaired due to
excessively high fish mercury concentrations and where the source of mercury is principally
atmospheric deposition
34
-------
Final Report
March2006
• Data assembly
• Calibration of D-MCM using estimates of long-term average annual conditions for Devils Lake.
The critical end point was mercury in walleye, but the calibration examined total and mercury
concentrations in each compartment for which data were available.
• Development of a long-term steady- state dose-response curve relating predicted long-term
average fish mercury concentrations to different levels of long-term continuous atmospheric
Hg(II) deposition. For example, if atmospheric deposition decreased to 50% of current levels and
was maintained at the lower value, one question that immediately arises is at what levels would
fish mercury concentrations ultimately stabilize? Model runs were carried out for several mercury
deposition rates to develop the steady-state response curve.
• Assessment of the predicted timing of the response of fish mercury concentrations to different
loadings of inorganic Hg(II).
• Sensitivity analysis of D-MCM predictions to various model input parameters, including
atmospheric deposition rates of mercury.
• Uncertainty analysis using a Monte Carlo approach to consider the uncertainty of predicted fish
mercury concentrations associated with the combined uncertainty associated with model inputs
related to mercury loads, site conditions and model constants.
• Assessment of the effects of year-to-year variations of atmospheric Hg(II) deposition.
Further information on the approach to model calibration, sensitivity analysis, uncertainty analysis and
assessment of year-to-year variations is provided below. Data assembly, model calibration and modeling
results are presented in later sections of this document.
Input Data
Sources of data used to parameterize and calibrate D-MCM to Devils Lake are summarized in
Table 13.
Table 13. Summary of data inputs by major category and type.
Data Type
Parameter Estimate and Source
Hydrologic Data
Precipitation
Measured (Krohelski and Batten, 1995). Period of record
is 1980 to 1991.
Surface water elevations
Measured (Krohelski and Batten, 1995). Period of record
is 1980 to 1991.
Surface flow
Derived from precipitation data and runoff coefficients.
Simulated values using calibrated hydrologic model
(Krohelski and Batten, 1995).
Physical Data
Bathymetric data
Supplied by Wisconsin DNR (Lathrop 1999c)
35
-------
Final Report
March2006
Surface water temperature
Observations from Wisconsin DNR 1986-97 (Lathrop
1999a, Foster 1999a)
Mercury Loadings
Wet Hg(II) deposition
Used monthly mean estimates from modeling by Myers et
al. (2003).
Dry Hg(II) deposition
Used monthly mean estimates from modeling by Myers
et. (2003).
Stream surface water concentrations - Hg(II)
Estimated from limited data: Herrin et al. 1998. n=2
Stream surface water concentrations - MeHg (unfiltered)
Estimated from limited data: Herrin et al. 1998. n=2
Surface Water Chemistry
DOC
No data
PH
Observations from Wisconsin DNR 1986-87 (Lathrop
1999a)
SO42-
Wisconsin DNR (Lathrop 1999b)
Chloride
Wisconsin DNR (Lathrop 1999b)
Dissolved oxygen
Observations from Wisconsin DNR 1986-97 (Foster
1999a)
Hg Concentrations (abiotic)
Water column water Hgtotal and MeHg (filtered and
unfiltered, epilimnion and hypolimnion)
Seasonal data from 1993-1996 (Wisconsin DNR)
Elemental Hg (DGM)
No data
Sediment Hg
No data
Sediment porewater chemistry
No data
Food web and associated mercury
Fish growth (walleye) and Hg concentrations
No site data; Calibrated based on regional growth rates
for northern and southern Wisconsin (Source Wisconsin
DNR)
Fish diets
No data
Fish biomasses
No data
Shiner Hg concentrations
Observations from Gorski(1999)
Bluegill Hg concentrations
Observations from Herrin etal. (1998)
Walleye Hg concentrations
Observations from Wisconsin DNR
Particle Dynamics
Sediment accumulation rates
Assigned ~ 1 mm/year on the basis of limited information
Sediment decomposition rates
Calibrated to result in non-depositional littoral zone and
~1 mm per year sedimentation in profundal zone.
Additional information describing inputs used in simulations is provided below:
Sediment Layer Thickness
36
-------
Final Report
March2006
The dynamics of a contaminant such as mercury in aquatic ecosystems, (e.g., its rate of response to
perturbations in loadings) is governed by its overall residence time. That residence time characteristically
is dictated by the size of the reservoir of labile mercury available in the surficial sediments to support
methylation as well as the rate at which that reservoir is buried. The size of the sediment pool in turn
reflects to what depth below the sediment-water interface mixing occurs such that particle reworking
(largely due to benthic infauna; see Boudreau, 1998) continually reintroduces previously buried Hg(II) into
the zone of active methylation. This depth of active mixing (which assumes all particles within the zone
are equally and continually well-mixed) can significantly affect predicted response dynamics when site
conditions or mercury loading regimes change. For example, the time predicted by D-MCM for fish
mercury concentrations to respond to changes in atmospheric mercury deposition rates can vary by decades
depending on the depth of the active sediment layer assigned (EPRI 2003). Thinner sediment layers
accelerate the predicted response in the lake system, particularly if methylation is assumed to occur
primarily in sediments.
Diffusion of mercury dissolved in the sediment porewater is calculated as a function of the difference for a
given mercury species between the overlying lakewater and the porewater:
JHgi = V ' I ^ ^i ] lakewater I ^ ^i ] porewater
whereis the areal flux of mercury species /' across the sediment-water interface (|ag/m2-yr), u is a mass
transfer coefficient (m/yr) and concentration is in |-ig/m3 (ng/L). u is related to the effective diffusion
coefficient, D^, used to describe mass transfer of solutes across the sediment-water interface (e.g., Lerman,
1979; Berner, 1981; Chapra and Reckow, 1983) by:
IK„
u = -
c/2
where x is the thickness of the mixed sediment layer. Typical values for u used in D-MCM range from
0.01 to 0.03 m/d coupled with a sediment layer thickness of 3 cm. To test the sensitivity of model
predictions of response dynamics following reductions in mercury loading, scenarios were also tested using
a model calibration with a 3 mm sediment layer.
Water Temperature
Monthly water temperatures for deterministic simulations were estimated using field measurements from
1986-1997 and are shown in Figure 11. Surface temperatures during the coldest months averaged between
1 to 2°C, while during the warmest months, temperatures approximated 24°C
37
-------
Final Report
March2006
n n
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
Month
Figure 11. Estimated mean monthly surface water temperatures for Devils Lake (Data source:
Lathrop 1999a, Foster 1999a).
Hydrology
Monthly precipitation values were derived from Krohelski and Batten (1995). Estimated mean monthly
precipitation values are shown in Figure 12. It is important to know that precipitation estimates are not
necessarily indicative of mercury loading to the lake surface. During periods of ice cover, mercury in
precipitation accumulated but was not loaded into the lake until snow melted, as described in subsequent
sections.
38
-------
Final Report
March2006
0.14
0.12
0.1
0.08
0.06
o
0)
o. 0.04
0.02
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
Month
Figure 12. Estimated mean monthly precipitation for Devils Lake. (Source: Krohelski and Batten
1995)
Inflows and outflows to Devils Lake were estimated on the basis of precipitation, using the hydrologic
budget of Krohelski and Batten (1995) to derive empirical relationships between precipitation and surface
inflow (Appendix II). Mean monthly inflows used in simulations are shown in Figure 13.
7000
6000
5000
5 4000
CO
E.
3 3000
2000
1000
-Inflow
Outflow
50 100 150 200 250
Day of Year (1 = January 1)
300
350
400
Figure 13. Estimated inflows and outflows for Devils Lake (derived from Krohelski and Batten 1995).
39
-------
Final Report
March2006
External Mercury Loadings
Mercury Loads Directly to the Lake Surface
Mean monthly wet and dry Hg(II) loading rates to Devils Lake were estimated using the following
approach:
• Monthly wet and dry Hg(II) deposition rates to Devils Lake were estimated using the results from the
original REMSAD simulations (Myers et al., 2003 )7.
• During months when the lake surface was covered with ice, atmospheric inputs of Hg to the lake
suraface were, in essence, suspended. Instead, atmospheric deposition during periods of ice cover was
assumed to accumulate in the snow/ice pack and enter the lake during the month ice breakup occurred
on the lake surface. Mercury accumulating on the lake surface in snow was reduced by 40% due to the
potential forHg(II) photoreduction in snow (M. Amyot, personal communication).
• Hg(II) loadings from wet and dry deposition were assumed to be equally available for methylation and
other reactions as they entered the lake water column. The resulting Hg(II) loading estimates directly
to the lake surface used in the simulations with current atmospheric Hg(II) deposition are shown in
Figure 14.
50
45
40
? 35
3 25
J 20 ,
^ 15 ¦
0 "I 1 1 1— 1 1 1 1 1 1 1
0 30 60 90 120 150 180 210 240 270 300 330 360
Day of year (1 = January 1)
Load originating from wet Hg(ll) deposition Load originating from dry Hg(ll) deposition
Figure 14. Estimated loading of Hg(II) directly to lake surface, originating from wet and dry
deposition.
7 As noted in Section 4.2.1, the REMSAD simulations were revised in response to peer review comments. Because of
timing constraints, and because the goal of this study is to demonstrate the approach to conducting a TMDL for an
atmospherically driven pollutant rather than conduct an actual TMDL, the aquatic modeling analysis based on the
original REMSAD calculations has not been revised to reflect the current simulated deposition results.
Ice cover ^ T , , • ,,
i/ Load during month
Ice
with ice breakup
cover
h-.
40
-------
Final Report
March2006
Methylmercury concentrations in precipitation resulted in annual wet deposition fluxes on the order of 0.1
(.ig/nr-vr. The assumed mean annual concentration used for methylmercury concentrations in precipitation
was 0.16 ng/L. This flux is consistent in magnitude with low methylmercury fluxes reported in wet
deposition by Lamborg el al. (1995) on the order of 0.06 (.ig/nr-vr.* Sensitivity analyses demonstrated that
the model is insensitive to variations in concentrations at these comparatively low levels when the other
atmospheric loadings are input at current levels. Dry deposition rates for methylmercury were assumed to
be negligible. Whether photodemethylation of accumulating MeHg in the snowpack occurs is unknown
but, based on the low ambient inputs of MeHg in wet deposition compared to the predicted in situ rates of
production, the effect of any errors associated with this flux on the MeHg budget for the lake are low.
Mercury inflows to Devils Lake
Inorganic Hg(II) and methylmercury loads to Devils Lake via inflows were estimated using the hydrology
described above in combination with assigned mercury concentration in the inflows. Measurements of
mercury concentrations in inflows were very limited: two observations for each of methylmercury and total
mercury from Messenger Creek on March 11 and March 20 in 1995 (Herrin el al., 1997). The average
concentrations of Hg(II) and methylmercury were 7.3 and 0.13 ng/L unfiltered respectively. These
concentrations were assumed to apply year round in lieu of any additional data. Estimated runoff of Hg(II)
to Devils Lake is shown in Figure 15 expressed as a flux normalized to the lake surface area.
sThese estimates for MeHg loadings from wet deposition are also supported by more recent measurements conducted at
four MDN sites in Wisconsin between January 1997 and May 2002 (B. Brunette, personal communication). VWM
concentrations across the period of record ranged from 0.085 to 0.108 ng/L, while annualized deposition fluxes ranged
from 0.06 to 0.08 |ig/m2.
41
-------
Final Report
March2006
14
12
-o 10
.><
CM
^ 8
O)
^ 6
O)
^ 4
2
0
30 60 90 120 150 180 210 240 270 300 330 360
Day of year (1 = January 1)
Figure 15. Estimated inflow Hg(II) loads to Devils Lake. Loading fluxes have been normalized to the
lake surface area.
Approach to Model Calibration
Calibration of D-MCM was conducted sequentially by optimizing the fit of the predicted values to the
observed data, one compartment at a time (i.e., surface water concentrations, sediment concentrations, and
biota concentrations). The first step was to calibrate the model with respect to observed epilimnetic
concentrations for the various mercury species. This was followed by calibrating with respect to observed
hypolimnetic concentrations of mercury. The next step was to calibrate the predicted sediment
concentrations. The final step was to calibrate the model with respect to the observed biota concentrations.
The process was also iterative - as a new compartment was calibrated, the calibration with respect to the
previous compartments was rechecked and recalibrated if necessary.
The calibration of D-MCM was initially conducted on the basis of estimated long-term average conditions.
Although Devils Lake was selected because of the comparative richness of mercury data available for
model calibration relative to other lakes that might require a mercury TMDL, there were nonetheless
several key gaps in data that made it to difficult to test the model calibration against observations. These
gaps included a lack of both sediment mass accumulation rates and sediment mercury concentrations, no
DOC concentrations, and no information on fish diets. The combined effect of limited site data and gaps in
the state of knowledge regarding mercury cycling was that alternative calibrations were possible to arrive at
the same end result (fish mercury concentration). Nonetheless, as the state of knowledge regarding
mercury cycling improves, the predictive strength of D-MCM should also improve through model
calibrations to a variety of sites. This would reduce the overall uncertainty associated with model
predictions for TMDLs at sites with the types and amount of data typically available.
Input loadings of mercury from atmospheric deposition (wet and dry) were based on the original simulated
wet and dry deposition fluxes predicted by REMSAD (Myers et al., 2003). No measurements of in situ
mercury fluxes were available. Stream and runoff loads were based on the long-term hydrologic budget
42
-------
Final Report
March2006
developed for the lake by Krohelski and Batten (1995) and limited measurements of stream mercury
concentrations (Herrin el al. 1998).
A key assumption that was dictated by the availability of data was that the extant datasets reasonably
represented current mercury fluxes and concentrations at Devils Lake. Secondly, it was assumed these
conditions were reasonably stable. Simulations were then run for 100 years with annual deposition patterns
and site conditions repeating year over year, often with monthly frequencies for inputs. The resulting
mercury concentrations and fluxes were then examined once the system had effectively stabilized (i.e.,
concentrations were not changing year to year). These results were reported on a weekly basis for the 101st
year of the simulation, so the the seasonality of the predictions could be examined.
Observed surface water concentrations of total and methylmercury were generally low. The model
calibration reasonably reflected the magnitude of both mercury species for concentrations in the epilimnion
(Figure 16 and Figure 17). Clear seasonal trends for Hg(II) are not apparent in the field data or model
results. There is a spike in predicted surface water Hg(II) concentrations in the spring resulting from the
mercury load derived from snowmelt on the lake surface. As discussed previous, the simulations assumed
that mercury deposition in winter months accumulated during ice cover and was then released to the lake
during the month that the spring thaw and ice-out occurred, having first reduced the accumulated Hg(II)
mass by 40% because of photo-reduction and re-emission from the snow/ice pack. Note that it was not
expected that the simulations of long-term average surface water mercury concentrations would match
precisely the observations on specific dates sampled from given years.
Seasonal trends were apparent for observed Hg(II) and methylmercury concentrations in the hypolimnion
in Devils Lake (Figure 18 and Figure 19, respectively). In both cases, concentrations started to increase at
some point during stratification, and continued to do so until overturn occurred in the fall. The model was
calibrated to attempt to generate comparable trends. Mechanisms that could be invoked to generate
hypolimnetic increases in Hg(II) or methylmercury concentrations included:
• Decomposition and release of Hg(II) and methylmercury settling into the hypolimnion from the
epilimnion on particles.
• Effects of anoxia:
o Triggering the dissolution of iron or manganese oxides at the sediment interface or settling
from the epilimnion, resulting in the release of Hg(II) and methylmercury originally adsorbed
to the oxidized (insoluble) forms of these two metals. Sulfides also will stably develop in the
hypolimnion during anoxia due to increased sulfate reduction. These sulfides could then
compete favorably to remove mercury from settling particles or from particles at the
sediment-water interface. In either case the result would be a tendency to build up unfiltered
hypolimnetic mercury concentrations (Hg(II) or methylmercury).
o Increasing methylation rates in hypolimnetic waters. This could occur as a result of increased
sulfate reduction or changes to the bioavailable fraction of Hg(II) available to methylate in
hypolimnetic waters. This mechanism would increase methylmercury concentrations in the
hypolimnion but not Hg(II) concentrations.
43
-------
Final Report
March2006
Herrin et al. (1998) have demonstrated that higher concentrations Hg(II) and methylmercury develop in the
hypolimnion of Devils Lake during summer stratification. These increases appear linked to the onset of
anoxia rather than stratification alone. For example, hypolimnetic concentrations of Hg(II) and
methylmercury showed no consistent tendency to increase following the onset of stratification until anoxia
occurred. In addition, increases in methylmercury concentrations were also greater when anoxia started
sooner. If these increases were driven primarily by the decomposition of particles from the epilimnion and
the subsequent release of mercury, one would not expect the correspondence of these trends with periods of
anoxia.
There was insufficient field information to resolve which of these potential mechanisms governed the
hypolimnetic buildup of Hg(II) and methylmercury in Devils Lake. D-MCM thus was calibrated such that
a combination of factors produced increasing mercury concentrations in the hypolimnion. These factors
included:
• A decrease in the ability of settling particles to compete for Hg(II) and methylmercury in the
hypolimnion, relative to the epilimnion. Since D-MCM is not yet capable of fully handling the
effects of sulfides, iron or manganese on Hg complexation, the desired effect was achieved by
reducing the mercury binding strengths of particles in the hypolimnion.
• Invoking some decomposition of settling particles at the sediment water interface.
• Invoking methylation in hypolimnetic waters, but only after anoxia occurred. Methylation in
epilimnetic sediments continued as before, but methylation in hypolimnetic sediments stopped in
the model when methylation began in the overlying hypolimnetic waters.
As a result of the above approach, the D-MCM calibration reasonably matched observations for
hypolimnetic methylmercury concentrations.
1.00
0.75
0.50
0.25
0.00
0.2 0.4 0.6
Time (yrs)
~ Observed 1990
¦ Observed 199'
A Observed 199!:
~ Observed 1996
D-MCM
Figure 16. D-MCM model calibration results for epilimnetic concentrations of Hg(II).
44
-------
Final Report
March2006
Figure 17. D-MCM model calibration results for epilimnetic concentrations of methylmercury.
0.2 0.4 0.6
Time (yrs)
~ Observed 1990
¦ Observed 199'
A Observed 199!:
« Observed 1996
—— D-MCM
Figure 18. D-MCM model calibration results for hypolimnetic concentrations of Hg(II).
45
-------
Final Report
March2006
Time (yrs)
Figure 19. D-MCM model calibration results for hypolimnetic concentrations of methylmercury.
No sediment mercury concentrations (total or methylmercury) were available for comparison with model
predictions. Measurements of mercury concentrations in walleye, bluegill and mimic shiner obtained from
the Wisconsin DNR (Gorski, 1999; Lathrop, 1999d) were used for calibrating biota concentrations. Figure
20 shows the comparison between predicted and observed concentrations of methylmercury in walleye as a
function of fish size.
Weight (kg)
Figure 20. D-MCM model calibration results for walleye methylmercury concentrations as a function
of fish weight. Included are observations for May, 1986, April, 1987, and April 1993.
46
-------
Final Report
March2006
4.3.3 Linkage of Mercury Loads to Fish Tissue Concentrations
A fundamental question to examine in this pilot TMDL study was the expected relationship between
atmospheric Hg(II) deposition and long term fish mercury concentrations. Once the model was calibrated
to the current atmospheric Hg(II) deposition estimate of 19.6 (.ig/nr-vr.. a series of simulations were
conducted to define the relationship predicted by D-MCM between various long-term atmospheric loading
rates and fish tissue Hg concentrations. Long-term simulations were conducted at atmospheric inputs
reduced by 10, 15, 25, 50 and 75% of the current total deposition estimate of 19.6 (.ig/nr-vr. Assumptions
also had to be made regarding the response of mercury loadings coming from surface runoff in response to
changing atmospheric inputs. Although Devils Lake is generally considered a seepage lake, it typically
receives over 50% of its hydraulic income from surface runoff. Because surface runoff also contains fairly
high concentrations of total mercury, it constitutes approximately 28% of the total mercury loading to the
lake. D-MCM currently is not capable of simulating mercury dynamics within the watershed soils.
Nonetheless the uncertainty in the magnitude of the response of Devils' Lake can be bounded by making
two key assumptions: (1) changing atmospheric inputs of mercury have no effect on watershed or surface
runoff loadings; and (2) changing atmospheric inputs of mercury have an equivalent (percentage) and
instantaneous response. For the second scenario, Hg(II) and methylmercury concentrations in inflows were
adjusted in proportion to Hg(II) deposition. This assumes that the methylmercury in surface runoff
originates from in situ production, and not from atmospheric inputs. Because there was no basis for
changing atmospheric loadings of methylmercury, that input was left unchanged. Predicted fish mercury
concentrations were compared after each simulation had run 150 years, producing essentially steady state
conditions. Annual cycles of site conditions and mercury deposition were repeated throughout the
simulation period.
Figure 21 shows the predicted long-term relationship between atmospheric mercury deposition and
mercury concentrations in age 5 walleye in Devils Lake for both watershed response scenarios. Both
scenarios produced a linear relationship with non-zero intercepts - in other words, even in the absence of
any atmospheric loadings of Hg(II), there is still some mercury being supplied to the lake, resulting in
methylation and incorporation into the food chain. For the no watershed response scenario, the large
intercept (0.2 ppm) reflects the constant input of Hg(II) and methylmercury from the watershed for all load
reduction scenarios. For the instantaneous watershed response, the small intercept (ca. 0.03 ppm) reflects a
residual source of mercury - principally exchange of ambient gaseous methylmercury across the air-
lakewater interface and the flux of methylmercury via atmospheric deposition assumed in MCM.
47
-------
Final Report
March2006
Load Reduction %
Figure 21. Predicted Hg concentrations in age 5 walleye as a function of different long term constant
annual rates of wet and dry Hg(II) deposition to Devils Lake. Two curves are shown: one shows the
predicted response with no changing inputs from the lake watershed; the other shows the predicted
response assuming that watershed loads of methyl and total mercury are reduced an equivalent
fractional amount as the total atmospheric deposition flux of Hg(II).
The predicted linear response between changes in atmospheric loadings and lake response needs to be
further verified by field data. The linearity of the D-MCM predictions reflects several important
assumptions, including:
1. It was assumed that two limiting factors govern methylation and demethylation rates: the supply
of available mercury and the rate of activity of the methylating and demethylating microbes. It
was also assumed that methylation occurred in sediments and that it was porewater Hg(II) or a
fraction of it being methylated. It is possible that cinnabar formation or other geochemical
constraints could buffer current levels of porewater Hg(II) at relatively constant concentrations
over a range of Hg(II) loading conditions. In this case, methylation rates in sediments might not
respond in a linear manner to a change in Hg(II) loading. Work also is needed to clarify the
location of methylation and demethylation in Devils Lake and whether the concentrations of
mercuiy available to drive each reaction respond linearly to different Hg(II) loadings to the
system. Although the limited time series data available for Devils Lake indicates that methylation
occurs in the hypolimnion (Appendix II, Section 7.1), it is unclear whether the locus of
methylation is in the sediments or the water column.
48
-------
Final Report
March2006
2. It was also assumed that microbial methylation and demethylation rates were limited by their
respective mercury substrates and organic matter decomposition rates. It is possible that at some
point these reactions could be limited by other factors such as the availability of carbon or a
micronutrient in short supply. For example, it was assumed that no effects of changing sulfate
loads that may accompany reductions in mercury emissions. Hrabik and Watras (2002) recently
attributed part of the observed declines in fish tissue mercury in Little Rock Lake in
northernWisconsin to reductions in atmospheric loadings of sulfate. In addition, in an effort to
mitigate eutrophication in Devils Lake and attendant problems such as anoxia, WDNR began in
2002 withdrawing hypolimnetic waters to remove phosphorus from the lake and ultimately the
sediments to prevent anoxia and reduce sediment releases of phosphorus. By delaying or
precluding the onset of anoxia, this mitigative action also likely will have a beneficial effect on
mercury cycling in Devils Lake by reducing methylation, either due to a reduction in sulfate-
reducing bacteria activity, or through a reduction in release of iron-bound Hg(II) when anoxia
occurs.
3. For the watershed response loading response scenario, it was assumed that Hg(II) and MeHg
concentrations in inflows would respond linearly to changes in atmospheric deposition. In fact, if
other factors emerge which suggest that the concentrations of Hg(II) or MeHg in Devils Lake do
not respond linearly to changes in atmospheric Hg(II) deposition, then the assumed linear
relationship for inflowing mercury would not be appropriate either. This would further contribute
to a non-linear response.
A second fundamental question addressed by the pilot mercury TMDL was: How rapidly will fish mercury
concentrations change following reductions in mercury loading? To address this question, simulations
were initially conducted for 100 years to ensure stable initial conditions; atmospheric loadings were then
reduced instantaneously as a step function and the simulations were continued for an additional 150 years.
The analyses included all the loading scenarios used to develop the long-term dose-response curve shown
in Figure 21. These analyses all assumed that the active exchange depth of the surficial sediments is 3 cm.
The relative response time required to approach steady state9 was also examined assuming that the
exchange depth is only 3 mm - i.e., one-tenth of what was assumed in all the other simulations. This
second simulation was run because there is increasing evidence from stable isotope studies such as
METAALICUS that aquatic ecosystems react very rapidly to new inputs of mercury (R. Harris, personal
communication).
Figure 22 shows the results of the time response analysis for both sediment exchange depth simulations for
a 50% load reduction. The analysis assumes no watershed response. [It should be noted that all load
reduction scenarios for a given sediment exchange depth regime yield the same relative response curve.] In
addition, the response curve is the same whether loadings from the watershed are reduced instantaneously
in concert with atmospheric loadings, or kept constant. This is because, under both watershed response
9 Relative response time is defined as here as the time required to reach a certain fraction or percentage of the predicted
total change in concentration that occurs between t =0 years and t = 150 years.
49
-------
Final Report
March2006
scenarios, the temporal response of the lake is governed solely by the hydraulic residence time of the lake
and the turnover rate of the surficial sediments.
The response time to reach 95% of the steady state concentration (xo.95) of mercury in age 5 walleye
following a load reduction is approximately 52.5 years. Not surprisingly, reducing the assumed exchange
depth to 3 mm has a profound effect on t0 95 by reducing it to approximately 9.6 years. This reduction in
response time largely results from a decrease in the size of the methylmercury pool available in the
sediments (i.e., all other factors being equal, reducing the size of the pool by 90%)). Approximately 2/3 of
the methylmercury production flux in Devils Lake is due to methylation occurring in the hypolimnion, with
10% and 25% originating from methylation in the surficial sediments and the watershed, respectively.
Recycling of methylmercury (whether originally produced in situ or derived from external inputs),
however, is dominated by diffusion of methylmercury across the sediment-water interface into the
hypolimnion. As organic matter containing methylmercury accretes in the sediments due to settling,
mineralization releases this methylmercury into the porewater, establishing a concentration gradient that
drives the upward diffusion of methylmercury from the sediments into the hypolimnion. Reducing the
sediment exchange depth reduces the size of this pool and thus decreases the residence time of
methylmercury in the lake. These results also illustrate that further research needs to be conducted on the
sediment exchange dynamics of mercury to produce more accurate predictions of the likely response time
to a given load reduction.
Figure 22. Relative response in age 5 walleye methylmercury concentrations predicted by D-MCM
as a function of time. Plot compares two scenarios: the red curve shows the response assuming a 3
mm sediment exchange depth; the blue curve shows the response assuming a 3 cm sediment exchange
depth.
50
-------
Final Report
March2006
5 UNCERTAINTY ANALYSIS
One of the critical elements in conducting a TMDL is addressing the degree of uncertainty inherent in
predicting a load required to produce an acceptable response in the numeric target of interest. In general,
sources of uncertainty lie with delineating the numeric target itself, uncertainty or error in the structure of
the predictive models because they are not an exact representation of the system being simulated,
uncertainty in the model parameters and forcing functions used in the models, and uncertainty in field
measurements. Often confounding the uncertainty due to measurement error is variability due to natural
processes (e.g., temporal or seasonal variations) and spatial heterogeneity. Characterizing both the effects
of uncertainty and natural variability on the requisite target loads is critical in developing a margin of safety
(MOS) that is appropriately protective of the water body. This section presents several different types of
analyses that were conducted as part of this pilot study to examine both the effects of uncertainty and
natural variability on model predictions.
5.1 Atmospheric Modeling Uncertainty Analyses
As part of the mercury atmospheric deposition analysis for the Devil's Lake TMDL Pilot, the importance of
several key modeling parameters and assumptions were examined. Most notably, these include alternative
speciation of mercury emissions and alternative background concentrations. These sensitivity simulations
were based on the original 1998 REMSAD simulations and data files.
The overall emission speciation profile for chlor-alkali facilities in the 1998 modeling was 97/3/0 % for
elemental, divalent gaseous, and particle-bound mercury, respectively. As there is only one chlor-alkali
facility in WI and its overall emissions, as reported by the facility, are among the highest in the state, it was
chosen as the subject of additional analysis into the sensitivity of modeling results to mercury emissions
speciation. In this sensitivity analysis, an alternative speciation profile of 70% elemental, 30% divalent, and
0% particle-bound was chosen, the same profile used in the EPA Mercury Report to Congress (EPA, 1997).
(Note that the 70/30/0 profile predates speciated testing performed at the Olin facility in Georgia by EPA
which formed the basis for the current profile of 97/3/0.) The impact of this change in speciation was
evaluated for Devil's Lake and areas surrounding the chlor-alkali facility (located approximately 100
kilometers north of Devil's Lake) by re-running REMSAD for the month of July 1998. It was found that such
a change could be expected to increase deposition near that particular source from approximately 4 to over 20
g/km2. However, due to the low level of emissions from the chlor-alkali plant (essentially considered to be
ground level or fugitive) the deposition of divalent mercury emissions drops off rapidly with distance from the
source and deposition beyond about 50 km from the source is very small. Hence, deposition at Devil's Lake
(which is also in a direction that is predominately up-wind from the source) was minimally affected by this
change in speciation.
The importance of considering background concentrations of mercury was highlighted by an annual
simulation in which the background was set to zero. Overall REMSAD performance as measured by
comparisons of modeled wet deposition values with measurements from the MDN stations indicated a
significant decrease in R2 from approximately 0.71 to approximately 0.41. (A similar result is also found in
the revised 2001 modeling. See section 5 A.) Background mercury was thus found to play an important role
in explaining mercury deposition at the largely remote locations (i.e., not near major mercury emission
51
-------
Final Report
March2006
sources) where MDN monitors are typically placed. In more urban settings, or at locations near major
sources, one would expect the relative importance of background to be lessened and local source influence
to increase.
5.1.1 Variations due to Meteorological Variability
D-MCM for Devils Lake was calibrated using simulated atmospheric deposition fluxes for both wet- and
dry-deposited Hg provided by REMSAD. Several different sources of uncertainty are inherent in these
predictions. First, there is uncertainty underlying the estimates of average annual wet and dry deposition
arising from (but not limited to):
1. Analytical and field collection errors for sites used to evaluate REMSAD and, in particular, poor
characterization of dry deposition fluxes of Hg,
2. Errors in modeling emission source and deposition relationships, including uncertainty in emission
estimates and/or their speciation, and uncertainty in simulating meteorological conditions.
Superimposed on the uncertainty in the true total Hg deposition flux are those variations likely to occur
even when sources (local and long-range) are constant and measurement errors are absent. These
variations are due to changes in both meteorology and rainfall patterns. Thus, when establishing target
loading rates for a TMDL, it is also important to consider the effects of these long-term fluctuations on the
target end-points.
Myers et al. (2006) used Classification and Regression Tree (CART) analysis to estimate the variability in
Hg deposition (wet and dry) owing to meteorological variations alone for a ten-year period. CART is an
objective procedure for reducing daily meteorology data for a specific site to a suite of meteorological
clusters or event types that produce a characteristic amount of Hg wet and dry deposition. Daily
meteorology and Hg deposition fluxes output from REMSAD for Devils Lake for 1998 were used by
Myers et al. (2005) to develop a reduced set of event types. By ascribing a given day to an event type
based on the ambient meteorological data, CART provides the basis for estimating Hg deposition in other
years for which REMSAD results were not available.
Figure 23 shows a box-and-whisker plot of the simulated dry and wet deposition fluxes by month across the
10-year simulation period. Table 14 and Table 15 provide the statistics (average and standard deviation,
respectively) for the simulated wet, dry and total deposition fluxes. In addition, the geometric mean and
the standard deviation of the ln-transformed fluxes are provided. Visual inspection of the frequency
distributions of the data indicate that the data are more appropriately represented by a ln-normal rather than
a simple normal distribution. As a result, subsequent calculations were based on the ln-transformed results.
To assess the effect of annual fluctuations in wet and dry deposition of mercury, a 500-year wet- and dry
deposition data set was synthesized using the statistical properties of the ln-transformed wet and dry
deposition results from the CART analysis summarized in Table 14 and Table 15. Wet and dry deposition
fluxes were computed independently for 500 years by randomly producing a monthly flux based on the
mean and standard deviation using a random number generator. The monthly fluxes were then summed for
each year to produce a 500-year series of total deposition fluxes. The geometric and arithmetic means of
the total deposition fluxes produced in the simulation were 18.6 and 19.9 (.ig/nr-vr. respectively. The
propagated relative standard deviation for the annual fluxes (expressed arithmetically) was 45.1%.
52
-------
Final Report
March2006
CM
E
O)
c
o
-4—»
'
-------
Final Report
March2006
5
0.625
-0.712
1.008
0.000
1.633
0.456
6
1.644
0.331
1.129
0.109
2.772
0.964
7
2.745
0.879
1.258
0.214
4.003
1.319
8
1.622
0.406
1.183
0.146
2.805
0.990
9
0.679
-0.711
0.967
-0.043
1.646
0.451
10
0.232
-1.623
0.765
-0.274
0.997
-0.011
11
0.283
-1.606
0.556
-0.598
0.840
-0.185
12
0.162
-2.098
0.604
-0.511
0.766
-0.287
Table 15. Statistical properties of 10-year CART simulation of wet and dry deposition of Hg at Devils
Lake, WI. CART analysis based on REMSAD simulation of mercury deposition in 1998. Table
shows standard deviation values by month for arithmetic and ln-transformed wet, dry, and total
deposition. Units for untransformed fluxes are ng/m2.
Standard Deviation
Month
Wet
Ln(Wet)
Dry
Ln(Dry)
Total
Ln(Total)
1
0.192
0.962
0.122
0.172
0.247
0.268
2
0.205
0.658
0.098
0.153
0.261
0.295
3
0.198
0.937
0.090
0.124
0.240
0.225
4
0.309
0.763
0.133
0.173
0.396
0.301
5
0.531
0.723
0.136
0.133
0.523
0.256
6
0.889
0.650
0.185
0.161
0.956
0.355
7
1.351
0.568
0.226
0.187
1.495
0.401
8
0.620
0.438
0.276
0.214
0.843
0.308
9
0.560
0.873
0.135
0.140
0.563
0.318
10
0.142
0.606
0.089
0.116
0.131
0.133
11
0.182
1.058
0.094
0.157
0.126
0.153
12
0.106
0.873
0.075
0.119
0.157
0.213
5.1.2 Variations Due to REMSAD Model Uncertainty for Monte Carlo
Analysis
In the original REMSAD simulation analysis, model performance error was evaluated by comparing
simulated annual wet deposition fluxes against fluxes measured at MDN sites and calculating the root mean
54
-------
Final Report
March2006
square error (RMSE). Errors in speciation were not considered because the MDN network does not
speciate different forms of wet-deposited Hg. The analysis of errors also was restricted to annual fluxes
because the MDN data for the simulation year (1998) has many data gaps for individual months such that,
for some months, some monitors may only have a week of data available.
For all available MDN sites (33 sites) across the entire US, the RMSE is 3.9 |-ig/m2 compared to an average
observed value for wet deposition is 10.0 |-ig/m2. Bias is only about 1%. For Wisconsin and Minnesota
sites only (n = 8), the RMSE is 4.0 |-ig/m2 and the average observation is 8.8 |-ig/m2. which is equivalent to
an error of 45.3 %. This estimate of error also includes a bias low in the model predictions of 2.7 |-ig/m2.
The error associated with the original REMSAD predictions on annual wet deposition was propagated by
assuming a coefficient of variation of 45.3%. Because no data sets yet exist to adequately evaluate the
performance of REMSAD with respect to its ability to predict dry deposition, the uncertainty (as reflected
in the coefficient of variation) in REMSAD predictions of dry deposition was assumed equivalent to that
for wet deposition.
The goal of the D-MCM Monte Carlo experiment in part is to test the effect of the overall uncertainty in the
estimate for the long-term total Hg deposition flux on predicted concentrations of mercury in walleye.
From a long-term perspective, the uncertainty in the estimate for the mean value for total deposition is of
interest. As a result, the average annual fluxes derived from the CART analysis for wet and dry deposition
were varied according to the standard error of the respective mean value. Superimposed on that variability
was the assumed error associated with REMSAD modeling, viz., 45.3% for wet and dry deposition
individually. Table 16 gives the resultant statistical properties of the simulated distribution for both
untransformed values, and ln-transformed values.
Table 16. Uncertainty in long-term estimates for wet, dry and total deposition of Hg at Devils Lake,
WI. Wet and dry estimates derived from CART analysis (n = 10), using standard error of the mean
to estimate uncertainty. Untransformed fluxes in |j.g/m2-yr. Total deposition simulated based on
uncertainty in average wet and dry deposition, and assuming that the error in modeling with
REMSAD was 45.3% for each flux component.
Metric Wet Deposition Dry Deposition Total Deposition
Average 9.261 10.258 19.053
Mean Standard Error 0.670 0.219
Standard Deviation 6.060
5.2 Aquatic Mercury Cycling Model Uncertainty Analyses
Uncertainty in model predictions can be identified using a variety of techniques. One approach is to
conduct sensitivity analysis in some form to examine and identify which variables in the D-MCM exert the
greatest effect on predicted fish mercury concentrations for Devils Lake. For this pilot study, two types of
sensitivity analyses were conducted to address this question. The first approach was to conduct a
traditional type of analysis where each parameter is varied by the same relative amount, without regard to
whether this value is actually likely to occur or is appropriate for the system of interest. The second
approach considered the range of actual or (if such information was not available) likely values a parameter
55
-------
Final Report
March2006
can assume. This latter approach is essentially a "minimum-maximum" analysis that examines sensitivity
in the context of likely or actual parameter distributional ranges. It thus defines the bounds of uncertainty
in model response related to a given variable. Details of the approaches and a presentation of the input
parameter values used in the analyses are presented in Appendix II.
5.2.1 "Minimum-Maximum" Sensitivity Analysis Results on Effects of
Parameter Uncertainty on Endpoint Prediction
It is important to understand which inputs in the model have the greatest influence on predicted fish
mercury concentrations. One traditional approach is to conduct a sensitivity analysis by varying selected
input parameters a given amount, e.g. + 50%. This type of sensitivity analysis does not account for the
actual variability/uncertainty associated with each input. Some inputs may vary proportionately more than
others, thus the analysis can either overestimate or underestimate the sensitivity of a given model
application to the uncertainty in parameter values. It also does not account for the fact that some inputs co-
vary, such as DOC concentration and light extinction.
An alternative approach considers the range of likely values an input parameter can assume. These limits
are used to help define the bounds of uncertain model response related to a given variable. The objective of
this study was to implement this approach with the D-MCM calibration to Devils Lake, and identify those
parameters whose bounds of uncertainty influence predicted Hg concentrations in walleye at the pilot
TMDL site to the greatest degree.
The approach used was a "minimum-maximum" analysis that considers sensitivity in the context of likely
parameter distributional ranges. To conduct the analysis, minimum and maximum expected values were
assigned to a subset of model inputs expected to have the potential to significantly affect predicted fish
mercury concentrations. For each input, D-MCM was run using the high and low limits for the input, while
all other parameters are held constant at their nominal values. The sensitivity index for the input (SI)
(Hoffman and Gardner, 1983) was calculated as:
SI = 1 - {low output result for the input/high output result for the input}
Note that the closer the SI is to 1 (i.e., the larger the difference between the high and low results), the more
sensitive the model is to the range in parameter uncertainty. This simple technique reflects both classical
sensitivity (the partial derivative of output Y to parameter X), and the range of variability of the input
parameter.
The results are shown in Figure 24. In some cases single input parameters were altered, while in other
cases multiple inputs were varied simultaneously because of their fundamental interdependency. There are
three key points to highlight in terms of the interpretation of the results of this sensitivity analysis:
1. This analysis provides an indication of how sensitive model predictions are to uncertainties associated
with various input parameters. It does not necessarily indicate which input parameters have the most
impact on predictions. For example, the concentration of Hg(II) in wet deposition can have a large
impact on predicted fish Hg levels, particularly in systems such as Devils Lake where the primary
input of Hg is atmospheric deposition. In other words, the model predictions are influenced strongly
by changes in the magnitude of this parameter. Since, however, the concentration of Hg(II) in wet
deposition is well-quantified and does not appear to vary substantially from year-to-year based on the
56
-------
Final Report
March2006
CART analysis (Section 4.2.5), the results indicate that the D-MCM model predictions are
comparatively insensitive to the inherent uncertainty in this parameter. Thus, the approach used in this
study provides insight into where future efforts to reduce input parameter uncertainties will produce
the greatest benefit towards reducing uncertainty in model predictions.
2. The sensitivity index does not necessarily address model parameter interdependencies. In reality,
groups of dependent parameters may move in concert with each other.
3. Some reaction rates and constants are tuned to obtain reasonable "net rates" for two reactions
comprising both forward and backward reaction paths -methylation/demethylation and photochemical
oxidation/reduction in surface waters are two such examples. Different combinations of these inputs
can result in the same net result - e.g., a combination of two very fast forward and backward reaction
rate constants can be selected to produce precisely the same rate of net production as two very slow
forward and backward rate constants. Until the understanding of the true unidirectional or gross (i.e.,
forward and backward, rather than net) reaction rates is improved, the model sensitivity to such inputs
is subject to the rates assumed in the model.
The parameters examined in the sensitivity analysis can be divided into the following categories:
• Inputs related to mercury loading
• Rate constants involving the production and removal of methylmercury;
• Inputs related to particle budgets and associated Hg fluxes;
• Inputs related to fish methylmercury uptake; and
• Inputs related to water quality
Parameter sensitivities for these categories are discussed below:
Inputs related to mercury loading
The model sensitivity to uncertainty associated with inputs related to atmospheric mercury loading was
moderate to low (precipitation rates, Hg concentrations in precipitation, dry Hg(II) deposition rates). This
seems unexpected at first, but it must be remembered that this analysis considers the sensitivity of predicted
fish mercury concentrations to the inherent uncertainty surrounding long-term mercury loading rates. This
is different than the issue of whether predicted fish mercury concentrations are sensitive to the actual rate of
mercury loading. Although Devils Lake receives most of its mercury supply from atmospheric Hg
deposition, the predicted range of wet and dry deposition fluxes over a 10-year period was low (see Section
4.2.5).
The model sensitivity to uncertainty associated with inputs related to watershed mercury loading was
moderate (inflow rate, Hg concentrations in inflows). In this case, uncertainty surrounding inflow Hg(II)
concentrations was high but, since the inflowing Hg(II) supplied less than 1/3 of the overall Hg(II) supply
to the lake, the effect of uncertainty in this flux was mitigated to a more modest level.
57
-------
Final Report
March2006
Rate Constants Ultimately Affecting Methylmercury Production and Destruction
Overall, model predictions were most sensitive to uncertainty associated with in situ methylmercury
production. Although results are shown for only for methylation, the model would be similarly affected by
uncertainty surrounding true biological demethylation rates. This is because: (1) in situ methylation is
predicted to be a major source of methylmercury for fish in Devils Lake; and (2) there is considerable
uncertainty regarding true rates. This result is not surprising, and is consistent with previous assessments
of research and development needs to better understand factors affecting fish mercury concentrations.
Uncertainty in rates of photochemical reduction of Hg(II) also emerged as important in the analysis. This
was because Hg(II) photoreduction was calibrated to be a major sink for Hg(II) in Devils Lake. Variations
in Hg(II) photoreduction rates affected concentrations of Hg(II) predicted to be available for methylation.
It should be noted that it was assumed there was no photo-oxidation of elemental Hg back to Hg(II), and
the calibrated Hg(II) photo-reduction rate is really a better measure of net reduction/oxidation than gross
Hg(II) reduction.
Uncertainty associated with true rates of methylmercury photodegradation also significantly affected
predicted fish mercury concentrations. Overall, uncertainty surrounding values for inputs affecting the in-
situ production or destruction of methylmercury had prominent effects on predicted fish mercury
concentrations. These model inputs warrant further attention.
Inputs related to particle budgets and associated Hg fluxes
Somewhat surprisingly, uncertainty regarding the amount of settling solids in the water column, and thus
the rate of sedimentation, did not markedly alter predicted fish mercury concentrations in Devils Lake.
This is because (1) sedimentation was not an important mechanism to directly bury methylmercury in
simulations (although it was important for total mercury), and increasing particle fluxes had nearly
offsetting predicted effects on methylation. More settling solids simultaneously resulted in more
decomposition and lower Hg(II) concentrations available to methylate. These trends tended to offset each
other in the model equation for methylation. Whether these model outcomes in fact approximate what
really transpires in terms of particle budget effects on methylation is not known.
Fish Mercury Uptake
Uncertainties associated with model inputs affecting fish mercury exposure via the diet moderately to
strongly affected predicted fish mercury concentrations. Fish typically obtain 90% or more of their
methylmercury via food, with gill uptake of methylmercury being a secondary factor. Benthos and
plankton are often both predicted to contribute significantly to the dietary methylmercury uptake by
predatory fish. This may occur directly by eating plankton and benthos, or indirectly by eating fish whose
mercury burdens originated from benthos or plankton. For the Devil's Lake sensitivity analysis, there were
no site measurements of benthic methylmercury concentrations or the appropriate value of partitioning to
apply for MeHg in benthos relative to sediment solids. This situation resulted in the predicted fish mercury
concentrations being sensitive to uncertainty associated with benthic methylmercury concentrations and the
partitioning of methylmercury into benthos in the model (bioconcentration factor or BCF.
58
-------
Final Report
March2006
Water Quality
The model was very sensitive to uncertainty associated with DOC concentrations in surface waters in the
lake. This is because (1) predicted fish mercury concentrations are sensitive to DOC in D-MCM, and (2)
there were no site data available for DOC for Devils Lake - i.e., the uncertainty regarding the true DOC
value was high. This situation could be remedied with some direct measurements of DOC in surface
waters, the hypolimnion, and porewaters in the lake.
Overall, the predicted concentrations of mercury in age 5 walleye in Devils Lake are most sensitive to
uncertainties associated with inputs related to the in situ production and destruction of methylmercury,
partitioning of methylmercury into the dietary items eaten by fish (e.g. benthos) and DOC concentrations
(due to a lack of data).
O)
T3
o
In-situ methylation constants
DOC in surface waters
BCF for MeHg in benthos
Hg(ll) photoreduction constant
MeHg photodegradation constant
Hg(ll) inflow concentration
BCF for MeHg in zooplankton
Hgll in precipitation
pH (surface and porewaters)
Flowrate
Fish growth rate (Wmax) using age basis
Dry Hgll deposition
Settling solids in water column
0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80
Sensitivity Index (Higher value = more sensitive)
Figure 24. Results of Minimum/Maximum Sensitivity Analysis.
59
-------
Final Report
March2006
5.2.2 Monte Carlo Analysis of Effects of Parameter Uncertainty on
Endpoint Prediction
The capability to conduct probabilistic "Monte Carlo" simulations has been added to D-MCM. The Monte
Carlo approach allows a model to be run in a manner that systematically attempts to estimate the
uncertainty associated with model predictions, based on considering simultaneously the uncertainty in the
various model inputs or parameters. In other words, the objective of the Monte Carlo analysis is to
examine the cumulative effects of the uncertainty in estimating the correct or true value for each of a suite
of input parameters or rate coefficients. This is fundamentally different to a deterministic simulation, for
which inputs and outputs have no associated uncertainty and sensitivity analyses where the effect of
variations in a parameter is assessed on a single-parameter basis.
The probabilistic version of D-MCM currently includes three types of probability distributions: normal,
lognormal and rectangular. Algorithms from Knuth (1990) were used for the normal and lognormal
distributions. Model input forms allow the user to enter distribution types and values for mean, minimum
and maximum values for input variables. Specification of minima and maxima for variables truncates
distributions so that unwanted extreme or negative values are avoided.
The Monte Carlo approach was used in this study to attempt to address uncertainty in predicted fish
mercury concentrations specifically associated with:
• Uncertainty in model input estimates of existing site conditions
• Uncertainty in input estimates for model constants (e.g. for rate reactions and Hg partitioning)
• Uncertainty in model input estimates for atmospheric mercury deposition and watershed mercury
runoff (see 5.1.2)
Uncertainty associated with estimates of existing site conditions
Datasets describing the physical, chemical and biological site conditions for Devils Lake were relatively
limited. Thus uncertainty existed for many model inputs when assigning values as estimates of current site
conditions. In cases where site-specific datasets existed spanning several years, distributions could be
developed directly. This was done for epilimnetic and hypolimnetc pH, water temperature and dissolved
oxgen. The analyses of uncertainty for all of these site conditions were conducted in a similar manner.
Monthly averages were obtained for each variable. The mean annual value was estimated as the mean of
the monthly values:
where //, is the estimate of the mean for the i'h month
60
-------
Final Report
March2006
The standard error of the mean was employed to characterize the uncertainty in the monthly mean values.
The error was then propagated to obtain an annual estimate, as follows:
where ct, is the estimate of the variance in the mean for the ith month.
Linear interpolation was used to obtain values for both mean and the standard error of the mean where
insufficient data existed (in order to compute annual estimates). Otherwise, an approach was taken to
assign normal distributions about an estimated mean value, with the standard deviation assumed to be 25%
of the mean value. Thus when the Monte Carlo simulation was carried out, approximately 95% of the
selected values for an input should lie within two standard deviations of the mean, ranging from 50% to
150% of the mean. Minimum and maximum values of 50% and 150% of the mean were also assigned,
meaning that all selected input values were within this range. This measure prevented unrealistic outlier
values to be selected for inputs. The inherent uncertainty in most parameters used as input to D-MCM for
Devils Lake has not been characterized, forcing the use of the latter approach. Clearly uncertainty for all
inputs is not the same, and under ideal circumstances, the inherent uncertainty would be quantified.
However, since this a pilot exercise demonstrating the application of the approach rather than a fully-
developed TMDL, making assumptions about parameter uncertainty when no information is available gives
us the opportunity to test this method in a TMDL-like framework.
Uncertainty associated with model constants
Most true Hg reaction rates have not been quantified (e.g., methylation, photoreduction), and model
constants have significant uncertainty associated with them outside the context of the Devils Lake
calibration. For the Devils Lake calibration, however, reaction rates were tuned to D-MCM to generate
mercury concentrations consistent with field observations. For example, if the methylation rate constant
was increased by 10 fold in the model calibration, fish mercury concentrations would have clearly been
outside the realm of observations. To allow for uncertainty associated with model constants, the analysis
used a similar approach as for site conditions with limited or no data - viz., assign normal distributions
about an estimated mean value, with the standard deviation assumed to be 25% of the mean calibrated
value. Minimum and maximum allowable values of 50% and 150% of the mean were also assigned,
meaning that all selected input values were within this range.
Uncertainty associated with watershed mercury loads to Devils Lake
Watershed mercury loads are the product of water flows and mercury concentrations in those flows. The
water budget for Devils Lake was derived from measured changes in lake stage, precipitation, and other
hydrologic parameters in conjunction with hydrologic modeling conducted by Krohelski and Batten (1995).
Uncertainty associated with inflows to Devils Lake was estimated by variations in rainfall measured
between 1980 and 1991 and empirical relationships derived between precipitation and surface water inflow
(see Appendix II). Data on mercury concentrations in inflows to Devils Lake were based on two
61
-------
Final Report
March2006
observations for both methyl and total mercury from Messenger Creek on March 11 and March 20 in 1995
(Herrin el al., 1998). To allow for uncertainty associated with inflowing concentrations of Hg(II) and
methylmercury, the same approach was used as for model constants and some site conditions - viz., assign
normal distributions about an estimated mean value, with the standard deviation assumed to be 25% of the
mean calibrated value. Minimum and maximum values of 50% and 150% of the mean were also assigned,
meaning that all selected input values were within this range.
Inter dependencies between model inputs
One of the key concerns surrounding the use of the Monte Carlo method for environmental simulations is
that unrealistic combinations of inputs can emerge if all inputs are assumed to behave independently, when
in fact some inputs are often related to others. To help reduce the potential for this situation to occur, D-
MCM has been designed such that some inputs can be made to depend on other inputs for probabilistic
simulations.
For this study, the following inputs were made to depend on other inputs:
• Water inflows, outflows and lake elevations were functions of precipitation.
• Light extinction was a function of dissolved organic carbon
• Hypolimnetic DOC was a function of epilimnetic DOC
• Hypolimnetic suspended solids was a function of epilimnetic suspended solids
The equations used to represent the above dependencies are documented in Appendix II. Monte Carlo
simulations were run for the base case calibration of current conditions, and for load reduction scenarios
involving a 3 cm thick sediment layer.
The Monte Carlo simulations were conducted for load reduction scenarios of 10, 25, 50, 75, and 95%. 100
simulations were conducted for each load reduction scenario. Watershed response for this simulation
experiment was considered to be instantaneous and proportional to the atmospheric load reduction. Time-
dependent results of the Monte Carlo simulations involving 50% load reductions are shown in Figure 25.
This simulation was based on the model calibration with a 3 cm sediment layer, and sediment layer
thickness was not a probabilistic input. The graph shows a considerable range represented by plus and
minus two standard deviations about the mean probabilistic result (i.e., encompassing approximately 95%
of the predicted distribution of values for age 5 walleye). The predicted fish mercury concentrations closely
approximated log-normal distributions. Similar trends emerged for the Monte Carlo simulations using other
load reductions.
These results suggest that the approach of assigning many inputs normal distributions with values ranging
from 50% to 150% of mean values (arithmetic range) and a standard deviation of 25% of the mean value,
resulted in significant uncertainty in predicted mercury concentrations in fish.
62
-------
Final Report
March2006
1.6
1.4
^ 1.2
O)
j*:
1.0
g 0.8
0
^ 0.6
^ 0.4
0.2
0.0
0 20 40 60 80 100
Time (yrs)
Figure 25. Time-dependent results for Monte Carlo simulation of the 50% load reduction scenario.
Results show as a function of time the mean value, minimum, maximum, and range encompassed by
± 2ct for age 5 walleye. The 2ct distributions were computed by assuming that the simulation results
were log-normally distributed (n = 100). Load is reduced as a step function at t = 0 years. Note:
lower limit defined by -2ct from the mean is nearly identical to and obscured by the minimum.
This uncertainty is also illustrated by the results from the Monte Carlo simulations for the relationship
between mercury load reduction and age 5 walleye tissue concentrations of methylmercury under quasi
steady state conditions (I = 150 years; Figure 26). Figure 26 shows the range of predicted walleye mercury
concentrations encompassed by 95% of the simulations (i.e., the plot excludes the lowest and highest 2.5%
of the simulated results) in conjunction with expected dose-response relationship from the calibrated model
(see Figure 21). [Note that the calibrated model dose-response for a 50% reduction in atmospheric loadings
differs somewhat from the response shown in Figure 25 because the latter response curve is based on the
geometric rather than the arithmetic mean from the Monte Carlo simulations.] Using the EPA tissue
criterion of 0.3 ppm as an example target, the results of this simulation experiment indicate that mercury
loads to Devils Lake would have to be reduced by nearly 79% to ensure with 95% confidence that the
target criterion would be met under steady state conditions, compared to an estimate of a requisite 50%
reduction indicated by the dose-response curve derived assuming no uncertainty. Note that these results of
course apply within the limitations of the assumptions imposed on the model analysis (e.g., instantaneous
watershed response and parameter uncertainty estimates).
63
-------
Final Report
March2006
20 40 60 80 100
Hg Load Reduction (%)
Figure 26. Steady state (t = 150 years) results for the Monte Carlo simulations showing the predicted
changes in age 5 walleye tissue concentrations of methylmercury as a function of load reduction.
Load reductions simulated are 10, 25, 50, 75, and 95%. Plot shows the range of predicted walleye
mercury concentrations encompassed by 95% of the simulations in conjunction with expected dose-
response relationship from the calibrated model (see Figure 21). Plot also shows the EPA tissue
criterion of 0.3 ppm for reference.
5.2.3 Effects of Annual Variablity in Atmospheric Deposition Rates on
Endpoint Predictions
Year-to-year variations in atmospheric Hg(II) deposition are expected to influence biota mercury
concentrations over the same time scale, particularly for species that are shorter-lived and occupy lower
trophic level niches. The extent that likely year-to-year variations in Hg(II) deposition to Devils Lake
might influence walleye mercury concentrations was evaluated by first synthesizing a set of 500 annual
Hg(II) deposition rates based on variations in meteorology and rainfall obtained from CART analysis (see
Section 4.2.5). This resulted in a synthesized data set for total deposition that varies each year about a
geometric mean value of 18.6 (.ig/nr-vr (arithmetic mean of 19.9 (.ig/nr-vr: Figure 27). Wet
methylmercury deposition also was varied proportionally in the same manner.
As discussed previously (Section 4.3.3), how mercury export from the watershed might respond to changes
in atmospheric loadings is unknown. Because of this uncertainty, two simulations were run: the first
assumed an instantaneous and proportional response to surface runoff loadings of mercury; the second
assumed that surface runoff loadings remained constant (see Section 4.3.3). Each simulation was run with
an initial 100-year "warm-up" period with Hg deposition rates held constant at current levels to ensure
model stability before imposing any changes in atmospheric deposition. The simulations then were
continued for an additional 500 years with year-to-year variations in Hg loadings based on the variations in
64
-------
Final Report
March2006
wet and dry Hg(II) deposition rates predicted by the CART analysis. Each annual input flux obtained from
the CART analysis was held constant as a step function for the entire year.
200
150
| 100
o
50
0
5 15 25 35 45 55 65 75 85 95 105 115 125135
Hg Deposition (|ag/m2-yr)
Figure 27. Frequency distribution of total mercury atmospheric deposition fluxes used to simulate
effects of year-to-year variations in atmospheric deposition rates on age 5 walleye mercury
concentrations. Number of years simulated = 500.
The results of the 500-year simulations are shown in Figure 28. The greatest variability was predicted for
the instantaneous watershed response scenario, with methylmercury concentrations in age 5 walleye
ranging from a low of 0.48 to a high of 0.72 ppm wet muscle. The variations for the no watershed response
were more muted, with a range of 0.51 to 0.66 ppm. Both watershed scenarios exhibited a mean
concentration of 0.56 ppm, with a coefficient of variation of 4.4 and 6.9% for the no watershed and
instantaneous watershed response scenarios respectively. Cumulative frequency distributions are shown in
Figure 29, which illustrates, for example, that 95% of the simulated age 5 walleye mercury concentrations
were below 0.61 and 0.64 ppm (i.e., 8.6 and 13.9%) for the no response and instantaneous watershed
response scenarios, respectively.
Compared to the variance in the deposition rates, the walleye response for both watershed scenarios was
considerably damped (Figure 30), and reflects: (1) buffering against short-term changes by the sediment
pool of mercury which, under the currently assumed sediment exchange depth of 3 cm, greatly exceeds the
amount of mercury entering the system on an areal basis any given year; and (2) the fact that fish tissue
concentrations reflect an integrated response over the 5 years of varying exposure, rather than simply
reflecting the current year deposition rate. For younger age class fish that integrate their mercury response
65
-------
Final Report
March2006
over short time intervals, the relative magnitude of their response to variations in atmospheric deposition
fluxes is greater (Figure 31). The water column, which contains a total mercury areal burden of ca. 5
|j.g/nr compared to annual average deposition flux of 19.9 |ag/m2; essentially provides no buffering on an
annual basis given that the mercury turnover rate is ca. 4 times per year.
0.45
0 100 200 300 400 500
Time (years)
t 1 1 i |—i i i i 1—i 1 1 i |—i i i i 1—i r
Instantaneous Watershed
Response
Average
Figure 28. Results from D-MCM simulation of effects of annual variability in atmospheric
deposition on age 5 walleye methylmercury concentrations. Plot shows the results for both
watershed scenarios (no response and instantaneous response to changes in atmospheric inputs).
66
-------
Final Report
March2006
CD
0.75
0.70
0.65
0.60
0.55
0.50
0.45
i r
Average
J L
.01 .1 1
"i—i—i—r
"i 1—i—i—i 1 r
Instantaneous Watershed
Response
No Watershed
Response
J I I L
_L
J I I L
_L
_L
5 10 2030 50 7080 9095
% Cumulative Distribution
99 99.999.99
Figure 29. Results from D-MCM simulation of effects of annual variability in atmospheric deposition
on age 5 walleye methylmercury concentrations shown in Figure 28 expressed as the percent
cumulative distribution.
67
-------
Final Report
March2006
c
o
w
o
Q.
0
Q
c
0
1
c
0
o
c
o
o
0
>
ro
0
200 300
Time (years)
500
Figure 30. Comparison of the magnitude of dynamic changes in atmospheric deposition and age 5
walleye response. Walleye simulation assumes instantaneous watershed response. Both parameters
are shown as the response relative to the mean over the 500-year simulation period.
68
-------
Final Report
March2006
10
8
6
4
0
0
-i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 r
Instantaneous Watershed
Response
No Watershed
Response
_i I i i i I i i i I i i i I i i i_
4 6
Age
8
10
Figure 31. Coefficient of variation for methylmercury concentrations in walleye in owing to year-to-
year variations in atmospheric deposition as a function of age class of fish.
69
-------
Final Report
March2006
6 DISCUSSION AND CONCLUSIONS
The primary objective of this pilot study was to explore how a TMDL or critical loading analysis could be
conducted when the pollutant of interest (in this case, mercury) is largely atmospheric in origin, and enters
aquatic ecosystems through atmospheric deposition either directly to the lake surface or indirectly to the
watershed. Similar to the parallel pilot study that was conducted in Florida (FDEP, 2002), this pilot study
was structured upon an integrated modeling approach linking an atmospheric emissions and deposition
model with an aquatic mercury cycling model that, in concert, were used to predict the relationship
between mercury emissions, deposition and, ultimately, concentrations in top predator fish. The
atmospheric modeling was designed to examine the spatial and temporal relationships between emission
sources and mercury deposition, while the aquatic cycling modeling was designed to examine the
relationship between atmospheric loadings of mercury to the receptor ecosystem and uptake of mercury by
top predator fish, and what biogeochemical processes influence that uptake. The inherent uncertainty in
these predicted relationships in turn was examined by analyzing the effects of both natural variability and
error in parameter estimation.
6.1 Atmospheric Modeling
The REMSAD modeling system was used in conjunction with statistical analysis to provide estimates of
wet and dry atmospheric mercury deposition inputs to water modeling efforts as part of the Devil's Lake
Mercury TMDL Pilot. In response to comments from reviewers, revised air modeling was conducted using
updated databases. The revised modeling provided estimates of a typical annual load of mercury to Devil's
Lake from wet and dry deposition as 5.9 g/km2 and 4.6 g/km2, respectively. (The original modeling had
given estimates of wet and dry deposition of 9.3 g/km2 and 10.3 g/km2, respectively.)
The algorithms used by REMSAD to describe atmospheric transport, transformation, and deposition of
mercury have undergone external peer review and the model's capabilities enhanced to reflect comments
raised in that review. The resulting performance of REMSAD is quite good on the basis of comparisons of
modeled wet deposition of mercury to measured wet deposition values. Performance evaluation of dry
deposition estimates is hindered by the lack of measurements for dry deposition.
Sensitivity analyses indicate that modeled deposition estimates can be greatly impacted by changes to
mercury emission speciation and to changes in background concentration assumptions, indicating that these
are areas that should be targeted for further study and analysis. In the revised modeling, new estimates of
boundary concentrations were obtained from the output of global mercury modeling.
A particularly useful component of REMSAD was its built-in capability to individually track emissions,
transport, and deposition of mercury derived from different sources throughout the model domain. This
"source tagging" enabled REMSAD to quantitatively track the effects of different sources on given
receptor sites. Such an approach is critical to the underlying requirement for TMDL's to ascribe or allocate
acceptable loads to sources contributing to the impairment of the water body. On the basis of existing
emission and meteorological inputs, mercury falling on Devil's Lake is largely dominated by background
70
-------
Final Report
March2006
concentrations10. In other parts of the modeling domain that are nearer sources of mercury, the role of
background concentrations is less significant.
6.2 Aquatic Modeling
Devils Lake was selected for analysis in this pilot study for three reasons: (1) the state of Wisconsin had
previously issued a fish consumption advisory because of excessively high concentrations of mercury in
walleye; (2) the lake derives all of its mercury inputs from atmospheric deposition; and (3) relative to other
lakes that might require a mercury TMDL, Devils Lake was comparatively rich with data for calibrating an
aquatic mercury cycling model. Modeling of the aquatic cycling of mercury in Devils Lake was conducted
using the Dynamic Mercury Cycling Model (D-MCM). Calibration of D-MCM was based on estimated
long-term average conditions and assumed the lake was essentially at steady state with respect to an annual
total loading of mercury of 19.6 g/km2 obtained from the original REMSAD simulations.
There were several key gaps in data that made it to difficult to test the model calibration against
observations. These gaps included limited data on total and methylmercury concentrations on surface
runoff discharged to the lake, a lack of both sediment mass accumulation rates and sediment mercury
concentrations, no data on DOC concentrations, and no information on fish diets. The combined effect of
limited site data and gaps in the state of knowledge regarding mercury cycling was that alternative
calibrations were possible to arrive at the same end result (fish mercury concentration) and that the aquatic
mercury cycling model was underconstrained11.
Internal production of MeHg is predicted to be the key source of MeHg to Devils Lake. Methylation in
Devils Lake appears to be closely linked to thermal stratification and the onset of anoxia in the
hypolimnion. Measured concentrations of both Hg(II) and MeHg both tend to increase in the hypolimnion
once anoxia occurs, and the seasonal timing of MeHg increases in the hypolimnion was linked to the timing
of anoxia. The temporal relationship between stratification, anoxia, and increases in MeHg suggest that
iron cycling (i.e., the cycle of oxidation and reduction of iron that results in iron precipitation [during oxic
conditions] and iron dissolution [during anoxic conditions] may be important in governing the availability
of Hg(II) for methylation, although other explanations may be possible (see Appendix II, Section 7.1).
This, in turn, suggests that methylation occurs in the hypolimnion rather than in just the sediments (as is
assumed during other parts of the year) in response to releases of sedimentary Hg(II) into the hypolimnion
10 Background concentration here is used to refer to effects of the boundary concentrations used for the model
simulations. Boundary concentrations represent contributions from all natural and anthropogenic emissions that are
external to the modeling domain and any recirculation of U.S. and Canada emissions that may occur. Since the tagging
simulations were not global in scope, it cannot be quantified at this time how much of background may be
recirculation.
11 In other words, insufficient data on key components of the mercury mass balance were available to ensure that the
model calibration is robust. For example, if the model for a given calibration was overpredicting total Hg in the water
column, this problem could be resolved simply by increasing particulate mercury removal via settling and burial.
Because there were no data to bound sediment accumulation rates of mercury, this pathway was essentially
unconstrained in the model, and there was little basis for choosing one value relative to another.
71
-------
Final Report
March2006
coupled with the ability of sulfate reducing bacteria to methylate mercury once anoxia occurs. By invoking
methylation in hypolimnetic waters (but only after anoxia occurred), the D-MCM calibration reasonably
matched observations for hypolimnetic methylmercury concentrations. For methylmercury, in situ
methylation was calibrated to supply -75% of the overall methylmercury supply to the lake, most of which
was calibrated to occur via methylation in the hypolimnion (65% versus 10% in sediments). Whether
significant rates of methylation can occur in the water column is controversial, and further research is
needed to establish the appropriate roles of the Fe/Mn redox cycle, settling particle dynamics and
decomposition rates on methylmercury dynamics in Devils Lake.
6.3 Load-Target Relationships
6.3.1 Linearity of Response
The D-MCM model currently predicts a linear relationship between atmospheric mercury deposition and
mercury concentrations in walleye. Because Devils Lake derives its water almost equally from two sources
- rainfall directly to the lake surface and surface runoff - a significant fraction (25-30%) of the mercury
inputs for the lake comes from surface runoff. Although this latter source likely is still ultimately derived
from atmospheric deposition, the time scale of changes in mercury fluxes from the watershed relative to
changes in atmospheric deposition is unknown. Depending upon whether mercury export from the
watershed eventually fully responds to changes in atmospheric inputs, the slope of the atmospheric input-
walleye response deviates from equivalency (i.e., 1:1; see Figure 21) with of course the greatest deviation
(less effective response) predicted for the scenario where watershed mercury export fluxes show no effect
of changing atmospheric inputs. The relatively small intercept of the instantaneous watershed response
scenario in Figure 21 reflects small inputs of methylmercury via gaseous atmospheric methylmercury
exchange into the lake that were assumed to persist in all scenarios. There is little information available to
quantify this flux or how it might change with mercury emissions to the atmosphere.
Based on the linear relationship between atmospheric mercury inputs and walleye mercury concentrations
predicted by D-MCM, a reduction of about 96 to 97% of current total annual mercury atmospheric
deposition rates would be needed for the mercury concentrations in a 5-year old walleye to be reduced to
the WDNR's present fish consumption advisory action level of 0.05 ppm. To meet the EPA tissue criterion
of 0.3 ppm, mercury loads (assuming an instantaneous watershed response) would have to be reduced by
approximately 50%.
Recent MDN data suggest that methylmercury concentrations in wet deposition collected in the Upper
Midwest may be significantly higher than previously believed (B. Brunette, personal communication), so
the effect of higher inputs of wet deposited methylmercury on mercury loads to Devils Lake should be
further examined. This is all the more important given that the current D-MCM model estimates indicate
that in excess of 95% of the mercury loadings to Devils Lake have to be eliminated to meet a concentration
of 0.05 ppm in age 5 walleye, without considering any effects of uncertainty.
The linearity of the D-MCM predictions to changes in atmospheric deposition rates reflects several
important assumptions, including:
• Methylation and demethylation were assumed to respond to the concentrations of available Hg(II)
or methylmercury, respectively, to act upon (1st order reactions). It should also be noted that
72
-------
Final Report
March2006
changes to Hg deposition could be accompanied by changes to sulfate loading. This could affect
methylation rates but was not part of the analysis.
• For the watershed response loading response scenario, it was assumed that Hg(II) and MeHg
concentrations in inflows would respond linearly to changes in atmospheric deposition. In fact, if
other factors emerge which suggest that the concentrations of Hg(II) or MeHg in Devils Lake do
not respond linearly to changes in atmospheric Hg(II) deposition, then the assumed linear
relationship for inflowing mercury would not be appropriate either. This would further contribute
to a non-linear response.
It is possible that cinnabar formation or other geochemical constraints could buffer current levels of
porewater Hg(II) at relatively constant concentrations over a range of Hg(II) loading conditions. It is also
possible that the fraction of dissolved Hg(II) that is available to methylate does not change linearly as
mercury loads increase. This could occur if some of the binding sites associated with dissolved phases or
particles reach or approach saturation as the levels of Hg(II) increase in the system. For these cases,
methylation rates in sediments likely would not respond in a linear manner to a change in Hg(II) loading.
Work also is needed to clarify where methylation and demethylation occur in Devils Lake and whether the
concentrations of mercury available for the reactions change in response in a linear way to different Hg(II)
loading to the system.
6.3.2 Factors Affecting Uptake of Mercury Inputs
Since most of the mercury entering Devils Lake is inorganic Hg(II), while most of the mercury in fish is
methylmercury, there are a number of steps that are required before changes in atmospheric Hg deposition
translate into relatively stable (steady-state) fish mercury concentrations. These include:
• Changes and stabilization in the watershed export of Hg(II) and methylmercury; and
• Changing in situ methylation rates as a result of new concentrations of Hg(II) that can be
methylated.
• After the methylation rate changes and reaches a new steady state level, it takes time for these
changes to cascade through the lake system, ultimately up the food web to top predators.
The rate at which mercury concentrations in Devils Lake are predicted to decline is largely a function of
three variables: (1) the hydraulic residence time of the lake; (2) the rate at which watershed fluxes of
mercury respond to changes in atmospheric inputs; and (3) the turnover rate of the pool of Hg(II) driving
methylation. Only the hydraulic residence time of Devils Lake is reasonably understood. Currently there
is almost no knowledge of the rate of watershed response, nor is there a clear understanding of the depth of
active exchange for Hg(II) (or MeHg) in surficial sediments. Since this latter pool is believed to supply
much of the Hg(II) supporting methylation, its magnitude is critical in governing how quickly
concentrations of methylmercury in the water column of the lake will react to changing atmospheric loads.
Likewise, the turnover time for sedimentary mercury is also affected by the net rates of sedimentation. The
response time predicted by D-MCM for age-5 walleye to reach 95% of the steady state concentration (t0.95)
of mercury following a load reduction is approximately 52.5 years given the current model calibration with
assumes a sediment exchange depth of 3 cm. Reducing the assumed exchange depth to 3 mm also reduced
T095 to approximately 9.6 years.
73
-------
Final Report
March2006
Other factors such as the true rates for key reactions such as methylation, demethylation, and Hg(II)
photoreduction followed by evasion can influence the in-lake response to changes in atmospheric loading.
Considerable uncertainty remains regarding true Hg reaction rates, and in the case of Devils Lake, data are
limited to quantify sedimentation rates. Likewise, a verified predictive method for reliably estimating
evasional fluxes of elemental mercury derived from the photoreduction of Hg(II) is not available.
Furthermore, the location of in situ methylation, depth of active sediments assumed for the model, and
structure of the food web supplied MeHg to top predators can also affect predicted response dynamics.
Uncertainty also remains regarding each of these factors for Devils Lake, and further research on these
model-sensitive parameters is required to more adequately define the relationship between loading and the
timing of biotic response.
The potential effects of watershed response time have not been included in this component of the analysis,
and the timing of the predicted response of fish mercury concentrations to changes in atmospheric
deposition in Devils Lake should be regarded as quite uncertain until mercury cycling in aquatic and
terrestrial systems is better clarified. Reducing mercury loads may prove to be only one component of a
multi-pronged strategy for reducing the effects of mercury loadings on fish tissue concentrations in Devils
Lake. Methylmercury concentrations were observed to build up in the hypolimnion of Devils Lake,
particularly following the onset of anoxia (Herrin el ol., 1998). There are several plausible links between
anoxia and methylmercury concentrations, including (1) under anoxia, solubilization of colloidal iron and
manganese in the surficial sediments could result in a pulsed release of Hg(II) into the uppermost
porewaters and the water column, thereby increasing methylation rates; (2) the same iron-manganese cycle
of colloidal solubility and dissolution tied to anoxia can also result in pulsed releases of sorbed
methylmercury; (3) under anoxia, sulfides could build up and better compete with settling solids and solids
at the sediment interface to hold Hg(II) and methylmercury in the dissolved phase; and (4) anoxic
conditions are required for biologically-mediated methylation to occur. In 2002, WDNR began withdrawal
of phosphorus-enriched waters from the hypolimnion as a means of controlling phosphorus releases in the
lake and ultimately to delay or eliminate the onset of anoxia in the hypolimnion and the associated iron
cycle dynamics of sorption and release of phosphorus.12 This could have the added benefit of reducing
methylmercury concentrations in the hypolimnion and the subsequent exposure of fish when this
methylmercury is mixed into the water column when fall mixing occurs. In addition, the possible role of
12 Withdrawal of hypo limnetic waters constitutes an additional sink for total and methylmercury, while replacement of
these waters with waters derived from Babbling Brook constitute an additional source. Assuming that (1) the rates of
withdrawal average around 613,000 m3/yr (the volume withdrawn in 2002); (2) there are no changes in the current
dynamics of total and methylmercury in the hypolimnion; and (3) assuming that total Hg concentrations in Babbling
Brook average 7 ng/L or less (equivalent to concentrations measured in surface runoff to Devils Lake), the net flux of
total mercury to the lake will increase by ca. 14% or less during the years the restoration project is in place. For
methymercury, however, withdrawal and replacement would result in a net loss of 0.06 |_ig/m2-yr, which is essentially
equivalent to the current loading from atmospheric deposition directly to the lake surface.
74
-------
Final Report
March2006
anthropogenic sulfate inputs to lakes in the Upper Midwest as an exacerbating influence on the
accumulation of methylmercury in fish is gaining more attention. For example, Hrabik and Watras (2002)
have attributed observed declines in mercury concentrations in fish tissue in Little Rock Lake in northern
Wisconsin to, in part, a decline in atmospheric inputs of sulfate.
6.4 Uncertainty in Model Predictions
One of the objectives of this pilot study was to examine the degree of uncertainty inherent in predicting a
load required to produce an acceptable response in the numeric target of interest. This uncertainty can lay
with several factors, including uncertainty or error in the mathematical structure of the predictive models
(i.e., the algorithms used in the model to simulate a process), uncertainty in the model parameters and
forcing functions used in the models, and uncertainty in field measurements. Often confounding the
uncertainty due to measurement error is variability due to natural processes (e.g., temporal or seasonal
variations) and spatial heterogeneity.
6.4.1 Uncertainty in REMSAD Predictions
The algorithms used by REMSAD to describe atmospheric transport, transformation, and deposition of
mercury have undergone external peer review and the model's capabilities enhanced to reflect comments
raised in that review. The resulting performance of REMSAD for wet deposition is quite good on the basis
of comparisons of modeled wet deposition of mercury to measured wet deposition values. The bias
between observed wet deposition fluxes and those predicted by REMSAD across the United States,
however, is sensitive to the boundary conditions assumed for reactive gaseous mercury in the atmosphere.
This assumption directly affects the amount of mercury deposited at a particular location such as Devils
Lake that is predicted to originate from outside the model domain. The concentration of RGM used for the
North and West boundaries of the model domain ranged from about 0.8 x 10~8 to 1.5 x 10~8 ppm in July, and
should be validated by further field measurements. Performance evaluation of dry deposition estimates is
hindered by the lack of measurements for dry deposition.
Sensitivity analysis with REMSAD indicates that mercury deposition from local sources at a given location
can be greatly influenced by the assumed speciation of low-altitude fugitive source emissions (such as a
chlor-alkali plant). These predictions are less sensitive for sources with large stacks from which emissions
would travel farther. In addition, the relative impact of this assumption on deposition fluxes decreases
rapidly with distance from the emission source; for example, REMSAD indicates that changes in deposition
due to emission speciation changes would be restricted to within about 50 km of the source. Since Devils
Lake lies approximately 100 km south of its major local source, the effects of emission speciation on
predicted deposition fluxes in this particular case would be barely noticeable. In order for a change in
speciation to show an effect at Devils Lake, the source would have to be very close, have a large emissions
rate, or be in a predominantly upwind location from the lake.
The above discussion touches on two factors (boundary concentrations and speciation of emissions) that
can strongly influence the REMSAD simulation results. Other factors can also influence the simulation
results, but typically these factors are not subject to adjustment in the current application. For instance, the
magnitude of emissions affects the simulation results, but these emissions represent the best estimate
available for use in the simulation and are not subject to variation in the current study. Likewise,
meteorological inputs represent data developed through established and tested methodologies and hence are
75
-------
Final Report
March2006
considered outside the realm of parameters that might be modified for this simulation. The REMSAD
model itself is formulated to be a tool adequate for the current application, subject to time constraints and
available data. Nevertheless, model formulation, particularly the air chemistry of mercury, is subject to
revision based on the availability of new or updated information. Depending on the nature of these
revisions, effects on simulation results could be large or small.
Unlike other types of uncertainty, the tagging methodology has an inherent means of estimating
uncertainty, at least uncertainty arising from strictly numerical effects. The REMSAD simulations included
a simulation of the overall mercury concentrations and deposition in addition to tags which are summed to
encompass all sources of mercury. Therefore, in theory, the sum of the simulation results for all tags
should be equal to the results of the simulation of the overall mercury concentrations. In practice, exact
equality cannot be achieved because the behavior of certain algorithms in the model are affected by the
magnitude and gradients of concentrations in the simulation. The simulation results show that for the
results on the 12-km grid, the sum of the tags differs from the overall mercury result by 2% or less.
Therefore, it is deduced that uncertainties in the tagging results of numerical origin are 2% or less.
Uncertainty in the REMSAD simulated wet deposition values was estimated through statistical comparison
with the MDN wet deposition observations. Uncertainty in the simulated dry deposition was assumed to be
comparable to the uncertainty in wet deposition since dry deposition observation data were not available.
The coefficient of variation utilized for the purpose of supplying deposition estimates to the D-MCM was
45.3%. This coefficient represents the uncertainty in the simulation results for 1998. In addition, since the
simulation results are used to drive D-MCM simulations for many years, an estimate of variability in the
results due to year-to-year variation in meteorology was also made. Based on a statistical analysis of 10
years of meteorological data, the relative standard deviation for this source of variation was set at 45.3%.
Further detail on these topics is presented in Section 5.1
6.4.2 Uncertainty in D-MCM Predictions
The effects of parameter uncertainty on the aquatic modeling were considered using several different
methods, including classical single parameter sensitivity analysis and analyzing uncertainty in toto using a
Monte Carlo module that had been recently developed and tested (Harris et ol., 2005). Previously,
parameter uncertainty analysis with D-MCM has been restricted to single parameter analyses. Although
very useful in terms of identifying for which parameters uncertainty critically influences model
performance, the single parameter uncertainty analysis provided only limited insight regarding the overall
uncertainty in D-MCM model predictions.
The D-MCM sensitivity analysis for Devils Lake indicates that the magnitude of predicted fish mercury
concentrations is most sensitive to uncertainties associated with methylmercury production, Hg(II)
photoreduction rates and DOC concentrations in surface waters, and methylmercury partitioning in the
lower food web, including benthos. The sensitivity to methylmercury production is not surprising, but the
sensitivity to DOC, as well as the absence of some other factors warrants comment. These findings are
based not on the model sensitivity to the values for a given parameter, but rather on the model sensitivity to
the uncertainty associated with these values. For example, predicted fish mercury concentrations are
indeed sensitive to rates of Hg(II) loading to the lake. However, based on the results from the REMSAD
CART analysis, the uncertainty assigned to atmospheric Hg deposition was low in the minimum-maximum
sensitivity analysis, and this factor did not emerge as exerting a large degree of uncertainty in predicted fish
76
-------
Final Report
March2006
mercury concentrations. Conversely, the model is sensitive to DOC concentrations in surface waters
because this is a model input that can have an important effect on predictions and for which no site
measurements were available. Thus uncertainty associated with DOC resulted in significant uncertainty in
predicted fish mercury concentrations. Uncertainty surrounding the response of watershed Hg export to
changes in atmospheric Hg deposition could strongly impact predictions of fish mercury levels following
load reductions to some lakes, but for Devils Lake most of the Hg(II) supply was from the atmosphere,
reducing the influence of this factor.
Monte Carlo analysis considered the effects of uncertainty in the following factors on predicted
concentrations of mercury in age five walleye:
• watershed loadings of mercury (due to hydrologic uncertainty, uncertainty in measured inflow
concentrations of total and methylmercury), and uncertainty in atmospheric inputs;
• model rate constants (e.g., methylationand demethylation rate constants); and
• existing site conditions (e.g., water column chemistry and fish growth rates)
In addition, several key interparameter relationships were developed so that parameters that intrinsically
co-vary (e.g., light extinction and the presence of dissolved organic carbon; lake hydrology and stage as a
function of precipitation) were not allowed to vary independently (and thus unrealistically) during the
course of the Monte Carlo analysis. In order to focus on the effects of uncertainty on long-term predictions
of mercury concentrations, the standard error of the mean rather than the sample variance alone was used to
characterize uncertainty in a given parameter. For parameters with insufficient data to characterize the
uncertainty in its mean value, it was assumed that the coefficient of variation was 25%, and that the range
in parameter values did not extend beyond 50 to 150% of the mean.
The Monte Carlo simulations were conducted for load reduction scenarios of 10, 25, 50, 75, and 95%. 100
simulations were conducted for each load reduction scenario. Watershed response for this simulation
experiment was considered to be instantaneous and proportional to the atmospheric load reduction. Using
the EPA tissue criterion of 0.3 ppm as an example target, the results of the Monte Carlo indicate that
mercury loads to Devils Lake would have to be reduced by nearly 79% to ensure with 95% confidence that
the target criterion would be met under steady state conditions, compared to an estimate of a requisite 50%
reduction indicated by the dose-response curve derived assuming no uncertainty.
As discussed above, for this pilot study, a significant degree of uncertainty is associated with model
predictions of the response of fish mercury concentrations to changes in atmospheric deposition, both in
terms of the magnitude and timing of the response. This is evident through the Monte Carlo analysis and
the simulations using the alternative calibration with a 3 mm sediment layer. The uncertainty emerges from
both gaps in the state of the science regarding mercury cycling and by limited data. As a result, a number
of substantively different calibrations of the D-MCM to Devils Lake are possible, leaving the model
underconstrained at this juncture. However, the tools used in this pilot study can be used to assess where
further data gathering efforts need to be undertaken in the future to develop a mercury TMDL. Such data
gathering in turn may assist in reducing the levels of uncertainty in the predicted relationship between
mercury loadings and target endpoint response.
77
-------
Final Report
March2006
7 REFERENCES
Arakawa, C., J. Yoshinaga, K. Nakai, H. Satoi, and K. Okamura. 2004. Effects of methylmercury
exposure on human reproduction. RMZ - Materials and Geoenvironment 51(1): 339-342.
Benjamin, S. G., D. Devenyi, S. S. Weygandt, K. J. Brundage, J. M. Brown, G. A. Grell, D. Kim, B. E.
Schwartz, T. G. Smirnova, T. L. Smith, and G. S. Manikin. 2003. An Hourly Assimilation/Forecast Cycle:
The RUC. (Submitted for publication).
Barbone, F., F. Valent, F. Pisa, M. Horvat, F. Daris, V. Fajon, D. Gibicar, and M. Logar. 2004. Risks and
benefits of fish consumption in an Italian population moderately exposed to methylmercury. RMZ -
Materials and Geoenvironment 51(1): 344-345.
Berner, R.A. 1981. Early diagenesis: A theoretical approach. Princeton University Press, Princeton, NJ.
Boudreau, B.P. 1998. Mean mixed depth of sediments: The wherefore and why. Limnol. Oceanogr. 43:
524-526.
Breiman, Leo, Jerome Friedman, Richard Olshen, and Charles Stone. Classification and Regression Trees.
Pacific Grove: Wadsworth, 1984.
Byun, D.W., Ching, J.K.S. (Eds). "Science Algorithms of the EPA Models-3 Community Multiscale Air
Quality (CMAQ) Modeling System." EPA-600/R-99/030. U.S. Environmental Protection Agency,
Research Triangle Park, NC. 1999.
Chan, H.M., A.M. Schumacher, A. Ferran, C. Loupelle, J. Holloway, and S. Weech. 2003. Impacts of
merecury on freshwater fish-eating wildlife and humans. Human and Ecological Risk Assessment 9:867-
883.
Chapra, S.C. and K.H. Reckhow. 1983. Engineering approaches for lake management. Volume 2:
Mechanistic modeling. Ann Arbor Science. Boston.
Douglas, S., B. Hudischewskyj, S. Beckmann, and T. Myers. 2005. Memorandum, Re: Comparison of
MM5- and RUC-Based Meteorological Input Fields for REMSAD Mercury Modeling (Revised), 14
December 2005.
E.H. Pechan & Assoc. 2000. Procedures for Developing Base Year and Future Year Mass and Modeling
Inventories for the Heavy-Duty Engine and Vehicle Standards and Highway Diesel Fuel (HDD)
Rulemaking. E.H. Pechan & Associates, Inc., Durham, North Carolina.
EPA. 1997. Mercury Study Report to Congress; Volume III: Fate and Transport of Mercury in the
Environment. EPA OAQPS/ORD, EPA-452/R-97-005,
(http://www.epa.gOv/ttn/atw/l 12nmerc/volume3 .pdf).
EPA. 2005. Emissions Inventory and Emissions Processing for the Clean Air Mercury Rule (CAMR), EPA
OAQPS.
78
-------
Final Report
March2006
EPA. 2005b. Technical Support Document for the Final Clean Air Interstate Rule: Air Quality Modeling,
EPA OAQPS (http://www.epa.gov/cair/pdfs/finaltech02.pdf).
EPA. 2005c. Technical Support Document for the Final Clean Air Mercury Rule: Air Quality Modeling,
EPA OAQPS (http://www.epa.gOv/ttn/atw/utility/aqm_oar-2002-0056-6130.pdf).
EPRI (Electronic Power Research Institute). 2002. Dynamic Mercury Cycling Model for Windows
98/NT/2000/XPTM: A Model for Mercury Cycling in Lakes. D-MCM Version 2.0. December 2002.
Electric Power Research Institute, Palo Alto, California.
. 2003. Factors Affecting the Predicted Response of Fish Mercury Concentrations to Changes in
Mercury Loading. Technical Report 1005521. Final report September 2003. Electronic Power Research
Institute, Palo Alto, California
Fitzgerald, W.F., D.R. Engstrom, R.P. Mason, and E.A.. Nater. 1998. The case for atmospheric
contamination of mercury in remote areas. Environ. Sci. Technol. 32: 1-7.
Foster, T.A. 1999. Devils Lake dissolved oxygen profiles.xls. Unpublished data. Wisconsin Department of
Natural Resources.
Friedlander, S. K. 1977. Smoke, Dust, and Haze: Fundamentals of Aerosol Behavior. John Wiley & Sons,
New York.
Gery, M. W., G. Z. Whitten, J. P Killus, and M. C. Dodge. 1989. "A Photochemical Kinetics Mechanism
for Urban and Regional Computer Modeling." J. Geophys. Res. 94:12,925-12,956.
Gorski, P. 1999. Shiner.xls. Unpublished data. Wisconsin Department of Natural Resources.
Gray, H. A., M. P. Ligocki, G. E. Moore, C. A. Emery, R. C. Kessler, J. P. Cohen, C. C. Chang, S. I.
Balestrini, S. G. Douglas, R. R. Schulhof, J. P. Killus, and C. S. Burton. 1991. Deterministic Modeling in
the Navajo Generating Station Visibility Study. Systems Applications International, San Rafael, California
(SYSAPP-91/045).
Hales, J. M., and S. L. Sutter. 1973. "Solubility of Sulfur Dioxide in Water at Low Concentrations." Atmos.
Environ. 7: 997-1001.
Hall, B. 1995. The gas phase oxidation of elemental mercury by ozone. Water Air Soil Pollution 80: 301-
315.
Harris, R.C. and R.A. Bodaly. 1998. Temperature, growth and dietary effects on fish mercury dynamics in
two Ontario Lakes. Biogeochemistry 40: 175-187.
Harris, R., C. Pollman, and D. Hutchinson. 2005. Wisconsin Pilot Mercury Total Maximum Daily Load
(TMDL) Study: Application of the Dynamic Mercury Cycling Model (D-MCM) to Devils Lake,
Wisconsin. Final Report Submitted to United States Environmental Protection Agency, Office of
Wetlands, Oceans and Watersheds, Washington, DC. Tetra Tech, Inc., Lafayette, CA.
Hedgecock, I.M. and N. Pirrone. 204. Chasing Quicksilver: Modeling the atmospheric lifetime of Hg°(g) in
the marine boundary layer at various latitudes. Environ. Sci. Technol. 38: 69-76.
79
-------
Final Report
March2006
Herrin, R.T., R.C. Lathrop, P.R. Gorski, and A.W. Andren. 1997. Hypolimnetic methylmercury and its
uptake by plankton during fall destratification: a key entry point of mercury into the food chain of Devil's
Lake, Wisconsin. Wisconsin Department of Natural Resources, Bureau of Science Services, unpublished
report.
Herrin, R.T., R.C. Lathrop, P.R. Gorski, and A.W. Andren. 1998. Hypolimnetic methylmercury and its
uptake by plankton during fall destratification: A key entry point of mercury into lake food chains? Limnol.
Oceanogr. 43(7) 1476-1486.
Hewett, S.W., and B.L. Johnson. 1992. Fish Bioenergetics Model 2. Published by the University of
Wisconsin Sea Grant Institute (WIS-SG-91-250).
Hoffman, F.O. and Gardner, R.H. 1983. Evaluation of Uncertainties in Environmental Radiological
Assessment and Environmental Dose Assessment Models, In: Radiological Assessment: A Textbook on
Environmental Dose Assessment, J.E. Till and H.R. Meyer (eds.), NUREG/CR-3332, ORNL-5968, U.S.
Nuclear Regulatory Commission, Washington D.C.
Hrabik, T.R. and C.J. Watras. 2002. Recent declines in mercury concentration in a freshwater fishery:
Isolating the effects of de-acidification and decreased atmospheric mercury deposition in Little Rock Lake.
Sci. Total Env. 297: 229-237.
Hudson R.J.M., S.A. Gherini, C.J. Watras, and D.B. Porcella. 1994. Modeling the Biogeochemical Cycle of
Mercury in Lakes: The Mercury Cycling Model (MCM) and Its Application to the MTL Study Lakes. In
Mercury Pollution - Integration and Synthesis, ed. C.J. Watras and J.W. Huckabee. CRC Press Inc., Lewis
Publishers.
ICF/SAI. 2002. Memorandum, Subject: Summary of Results of Mercury Sensitivity Simulations, 2 July
2002.
Keeler, G.J., F.J. Marsik, K.I Al-Wali, and J.T. Dvonch. 2001. Modeled Deposition of Speciated Mercury
to the SFWMD Water Conservation Area 3A: 22 June 1995 to 21 June 1996. Project Description and
Results. The University of Michigan Air Quality Laboratory, Ann Arbor, Michigan 48109.
Kessler, E. 1969. "On the Distribution and Continuity of Water Substance in Atmospheric Circulation".
Meteor. Monogr., No. 32, 84 pp.
Knuth, D.E. 1990. The Art of Computer Programming, Vol. 2., 2nd ed. Addison-Wesley.
Krohelski, J.T. and W.G. Batten. 1995. Simulation of stage and the hydrologic budget of Devils Lake,
Sauk County, Wisconsin. Undated US Geological Survey report.
Lamborg, C.H., W.F. Fitzgerald, G.M. Vandal and K.R. Rolfhus. 1995. Atmospheric Mercury in Northern
Wisconsin: Sources and Species. Water, Air and Soil Pollution 80:189-195.
Lamborg, C.H., W.F. Fitzgerald, A.W.H. Damman, J.M. Benoit, P.H. Balcom, and D.R. Engstrom.
Modern and historic atmospheric mercury fluxes in both hemispheres: Global and regional mercury cycling
implications. Global Biogeochemical Cycles 16: 51-1 - 51-11.
80
-------
Final Report
March2006
Lathrop, R.C. 1999a. Devils Lake 1986 & 1987, pH & Alkalinity.xls. Unpublished Data: Wisconsin
Department of Natural Resources.
. 1999b. Personal communication (email) sent to R. Harris July 27, 1999.
. 1999c. Devils Lake Water Levels 1980-1999.tif. Unpublished data.
. 1999d. Wisconsin Department of Natural Resources - Mercury Samples - Mercury in Fish from
Devils Lake.
Lathrop, R. 2004. Personal communication. Wisconsin Department of Natural Resources.
Lathrop, R.C., K.C. Noonan, P.M. Guenther, T.L. Brasino, and P.W. Rasmussen. 1989. Mercury Levels in
Walleyes from Wisconsin Lakes of Different Water and Sediment Chemistry Characteristics. Technical
Bulletin No. 163. Wisconsin Department of Natural Resources.
Lerman, A. 1979. Geochemical processes: Water and sediment environments. Wiley-Interscience. New
York.
Lin, C. A. and S.O. Pehkonen, "The Chemistry of Atmospheric Mercury: A Review," Atmos. Environ, 33,
2067-2079, 1999.
Louis, J. F. 1979. "Parametric Model of Vertical Eddy Fluxes in the Atmosphere." Boundary Layer
Meteorology, 17:187-202.
Munthe, J., and A. Palma. 2003. Chapter 10: Atmospheric cycling of mercury and persistent organic
pollutants. Overview of Subproject MEPOP.
http://www.gsf.de/eurotrac/publications/et2 fin rep/fr 10 mepop.pdf.
Myers, T. 2004. Memorandum to Dwight Atkinson, Subject: Recent revisions to REMSAD mercury
chemistry routine, 20 January 2004
Myers, T. C., Y. H. Wei, G. Glass, J. Mangahas, and B. Hudischewskyj. 2002. Application of the
REMSAD Modeling System to Estimate the Deposition of Nitrogen and Mercury (Draft). ICF
Consulting/SAI, San Francisco, CA.
Myers, T.C.. Y.H. Hua Wei, B. Hudischewskyj, and S. Douglas. 2003. Application of the REMSAD
Modeling System to Estimate the Deposition of Mercury in Support of the Wisconsin TMDL Pilot (Draft).
Draft report submitted to USEPA OW, Washington, DC. Submitted by Systems Applications International,
Inc./ICF Consulting, San Rafael, CA.
Myers, T.C.. Y.H. Hua Wei, B. Hudischewskyj, and S. Douglas. 2005. Application of the REMSAD
Modeling System to Estimate the Deposition of Mercury in Support of the Wisconsin TMDL Pilot
(Revised). Final report submitted to USEPA OW, Washington, DC. Submitted by Systems Applications
International, Inc./ICF Consulting, San Rafael, CA.
Pacyna, E.G. and J.F. Pacyna. 2002. Global emission of mercury from anthropogenic sources inn 1995.
Water Air Soil Poll. 137: 149-165.
81
-------
Final Report
March2006
Pai, P., P. Karamchandani, and C. Seigneur. 1997. "Simulation of the Regional Atmospheric Transport and
Fate of Mercury Using a Comprehensive Eulerian Model." Atmos. Env., 31:2717-2732.
Pai, P., P. Karamchandani, and C. Seigneur. 1999. "Sensitivity of Simulated Atmospheric Mercury
Concentrations and Deposition to Model Input Parameters." JGR, 104:13855-13868.
Pollman, C.D. and D.B. Porcella. 2003. Assessment of trends in mercury-related data sets. J. Phys. IV
France 107: 1083-1086.
Ryaboshapko, A., R. Bullock, R. Ebinghaus, I. Ilyin, K. Lohman, J. Munthe, G. Petersen, C. Seigneur, and
I. Wangberg. 2002. "Comparison of mercury chemistry models," Atmos. Env., 36:3881-3898.
SAI. 2005. User's Guide to the Regional Modeling System for Aerosols and Deposition (REMSAD),
Version 7. ICF Consulting/SAI, San Francisco, CA. (Available atwww.remsad.com)
SAI. 2002. User's Guide to the Regional Modeling System for Aerosols and Deposition (REMSAD),
Version 7. ICF Consulting/SAI, San Francisco, CA.
SAI. 1999. User's Guide to the Variable-Grid Urban Airshed Model (UAM-V). Systems Applications
International, San Rafael, California..
Schroeder, W.H. and J. Munthe. 1998. Atmospheric mercury: An overview. Atmos. Environ. 32: 809-
822.
Schliiter, K. 2000. "Review: evaporation of mercury from soils. An integration and synthesis of current
knowledge." Env. Geol. 39:249-271.
Scire, J. S. 1991. A Review of the UAM-V Dry Deposition Algorithm and Recommendations for Dry
Deposition Modeling in the LMOS Study Region. Sigma Research Corporation, Westford, Massachusetts
(Document A195-100).
Scott, B. C. 1978. "Parameterization of Sulfate Removal by Precipitation." J. Appl. Meteor., 17, 1375-
1389.
Seigneur, C., H. Abeck, G. Chia, M. Reinhard, N. Bloom, E. Prestbo, and P. Saxena. 1998. "Mercury
adsorption to elemental carbon (soot) particles and armospheric particulate matter." Atmos. Env., 32:2649-
2657.
Seigneur, C., G. Hidy, I. Tombach, J. Vimont, and P. Amar. 1999. Scientific Peer-Review of the
Regulatory Modeling System for Aerosols and Deposition (REMSAD). The KEVRIC Company, Inc.
Durham, NC.
Seigneur, C., K Vijayaraghavan, K. Lohman, P. Karamchandani, C. Scott, and L. Levin. 2003. Simulation
of fate and transport of mercury in north America. J. Phys. IV France 107: 1209-1212.
Shia, R., C. Seigneur, P. Pai, M. Ko, and N. D. Sze. 1999. "Global simulation of atmospheric mercury
concentrations and deposition fluxes." J. Geophys. Res., 99 :23,747 - 23,760.
82
-------
Final Report
March2006
Slinn, W.G.N., L. Hasse, B. B. Hicks, A. W. Hogan, D. Lai, P. S. Liss, K. O. Munnich, G. A. Sehmel, and
O. Vittori. 1978. "Some Aspects of the Transfer of Atmospheric Trace Constituents Past the Air-Sea
Interface." Atmos. Environ., 14:1013-1016.
Sofiev, M.A., and M.V. Galperin. 2000. "Long-Range Transport Modeling of the Particle-Carried and
Persistent Toxic Pollutants." Proceedings of the Millennium NATO/CCMS ITM, May 15-19, 2000,
Boulder, CO, AMS, pp. 148-155.
Steinberg, D. and P. Colla. 1997. CART—Classification and /Regression Trees. San Diego, CA: Salford
Systems.
Syrakov, D. 1998. "Long-Term Calculation of Bulgarian Impact in Mercury Pollution over SE Europe."
Bulgarian contribution to 1998 Annual Report for 1998, NIMH, EMEP/MSC-E, Sofia-Moscow.
Tetra Tech, Inc. 1996. Regional Mercury Cycling Model: A Model for Mercury Cycling in Lakes (R-MCM
Version 1.0b Beta), Draft User's Guide and Technical Reference. Prepared for the Electric Power Research
Institute and Wisconsin Department of Natural Resources. December 1996.
Vijayaraghavan, K., K. Lohman, P. Karamchandani, and C. Seigneur. 2002. Modeling Deposition of
Atmospheric Mercury in Wisconsin. AER, Inc., San Ramon, California.
Van Loon, L. L., E. A. Mader, and S. L. Scott. 2001. "Sulfite Stabilization and Reduction of the Aqueous
Mercuric Ion: Kinetic Determination of Sequential Formation Constants." J. Phys. Chem. A 105:3190-3195
Wesley, M. L. and B. M. Lesht. 1988. Comparison of the RADM Dry Deposition Module with Site-
Specific Routines for Inferring Dry Deposition. US EPA, Research Triangle Park, NC (EPA/600/54-
88/027).
Wesley, M. L. 1989. "Parameterization of Surface Resistances to Gaseous Dry Deposition in Regional-
Scale Numerical Models." Atmos. Environ., 23:1293-1304.
Wiener, J.G. 1983. Comparative analyses of fish populations in naturally acidic and circumneutral lakes in
northcentral Wisconsin. FWS/OBS-80/40.16 US Fish and Wildlife Service, Eastern Energy Land Use
Team.
Wiener, J.G. 1987. Metal contamination of fish in low-pH lakes and potential implications for piscivorous
wildlife. Trans. N. Am. Wildlife Resources Conf. 52: 645-657.
WDNR and WDH (Wisconsin Department of Natural Resources and Wisconsin Division of Health). 2003.
Health advisory for people who eat sport fish from Wisconsin waters. Madison, WI. 7 pp.
Yantosca, B. 2004. "GEOS-CHEMv7-01-02 User's Guide," Atmospheric Chemistry Modeling Group,
Harvard University.
83
------- |