United States Environmental Protection Agency Office of Water (4305) EPA-820-R-12-015 August 2012 &EPA AQUATOX (RELEASE 3.1) MODELING ENVIRONMENTAL FATE AND ECOLOGICAL EFFECTS IN AQUATIC ECOSYSTEMS VOLUME 2: TECHNICAL DOCUMENTATION ------- {This Page Left Blank, Back of Cover} ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION AQUATOX (RELEASE 3.1) MODELING ENVIRONMENTAL FATE AND ECOLOGICAL EFFECTS IN AQUATIC ECOSYSTEMS VOLUME 2: TECHNICAL DOCUMENTATION Richard A. Park1 and Jonathan S. Clough2 AUGUST 2012 U.S. ENVIRONMENTAL PROTECTION AGENCY OFFICE OF WATER OFFICE OF SCIENCE AND TECHNOLOGY WASHINGTON DC 20460 Modeling, Diamondhead MS 39525 2Warren Pinnacle Consulting, Inc., Warren VT 05674 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION DISCLAIMER This document describes the scientific and technical background of the aquatic ecosystem model AQUATOX, Release 3.1. Anticipated users of this document include persons who are interested in using the model, including but not limited to researchers and regulators. The model described in this document is not required, and the document does not change any legal requirements or impose legally binding requirements on EPA, states, tribes or the regulated community. This document has been approved for publication by the Office of Science and Technology, Office of Water, U.S. Environmental Protection Agency. Mention of trade names, commercial products or organizations does not imply endorsement or recommendation for use. ACKNOWLEDGMENTS This model has been developed and documented by Richard A. Park of Eco Modeling and by Jonathan S. Clough of Warren Pinnacle Consulting, Inc. under subcontract to Eco Modeling. The work was funded with Federal funds from the U.S. Environmental Protection Agency, Office of Science and Technology under contract numbers 68-C-01-0037 to AQUA TERRA Consultants, Anthony Donigian, Work Assignment Leader; and EP-C-12-006 to Horsley Witten Group, Inc., Nigel Pickering, Work Assignment Leader. Integration of Interspecies Correlation Estimation (Web-ICE) was made possible due to the work of US. EPA Office of Research and Development Gulf Breeze, the University of Missouri-Columbia, and the US Geological Survey. The assistance, advice, and comments of the EPA work assignment manager, Marjorie Coombs Wellman of the Standards and Health Protection Division, Office of Science and Technology have been of great value in developing this model and preparing this report. Further technical and financial support from David A. Mauriello, Rufus Morison, and Donald Rodier of the Office of Pollution Prevention and Toxics is gratefully acknowledged. Marietta Echeverria, Office of Pesticide Program, contributed to the integrity of the model through her careful analysis and comparison with EXAMS. Amy Polaczyk and Marco Propato, Warren Pinnacle Consulting, made valuable contributions to the model formulation and testing. Release 2 of this model underwent independent peer review by Donald DeAngelis, Robert Pastorok, and Frieda Taub; and Release 3 underwent peer review by Marty Matlock, Damian Preziosi, and Frieda Taub. Their diligence is greatly appreciated. ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION TABLE OF CONTENTS DISCLAIMER ii TABLE OF CONTENTS iii PREFACE viii 1. INTRODUCTION 1 1.1 Overview 1 1.2 Background 4 1.3 The Multi-Segment Version 6 1.4 The Estuary Model 6 1.5 The PFA Model 7 1.6 AQUATOX Release 3.1 Overview 7 1.7 Comparison with Other Models 9 1.8 Intended Application of AQUATOX 10 2. SIMULATION MODELING 12 2.1 Temporal and Spatial Resolution and Numerical Stability 12 2.2 Results Reporting 15 2.3 Input Data 16 2.4 Sensitivity Analysis 17 2.5 Uncertainty Analysis 19 2.6 Calibration and Validation 24 3. PHYSICAL CHARACTERISTICS 41 3.1 Morphometry 41 Volume 41 Bathymetric Approximations 44 Dynamic Mean Depth 47 Habitat Disaggregation 47 3.2 Velocity 48 3.3 Washout 49 3.4 Stratification and Mixing 50 Modeling Reservoirs and Stratification Options 54 3.5 Temperature 55 3.6 Light 56 Hourly Light 58 3.7 Wind 59 3.8 Multi-Segment Model 60 Stratification and the Multi-Segment Model 62 State Variable Movement in the Multi-Segment Model 62 4. BIOTA 64 4.1 Algae 65 iii ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION Light Limitation 69 Adaptive Light 74 Nutrient Limitation 75 Current Limitation 77 Adjustment for Suboptimal Temperature 78 Algal Respiration 80 Photorespiration 81 Algal Mortality 82 Sinking 83 Washout and Sloughing 84 Detrital Accumulation in Periphyton 89 Chlorophyll a 89 Phytoplankton and Zooplankton Residence Time 90 Periphyton-Phytoplankton Link 91 4.2 Macrophytes 92 4.3 Animals 96 Consumption, Defecation, Predation, and Fishing 98 Respiration 102 Excretion 105 Nonpredatory Mortality 106 Suspended Sediment Effects 107 Gamete Loss and Recruitment 123 Washout and Drift 125 Vertical Migration 126 Migration Across Segments 127 Anadromous Migration Model 128 Promotion and Emergence 129 4.4 Aquatic Dependent Vertebrates 130 4.5 Steinhaus Similarity Index 130 4.6 Biological Metrics 131 Invertebrate Biotic Indices 135 5. REMINERALIZATION 137 5.1 Detritus 137 Detrital Formation 141 Colonization 142 Decomposition 144 Sedimentation 147 5.2 Nitrogen 150 Assimilation 152 Nitrification and Denitrification 153 lonization of Ammonia 155 Ammonia Toxicity 157 5.3 Phosphorus 158 5.4 Nutrient Mass Balance 160 Variable Stoichiometry 160 Nutrient Loading Variables 161 iv ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION Nutrient Output Variables 161 Mass Balance of Nutrients 162 5.5 Dissolved Oxygen 169 Diel Oxygen 173 Lethal Effects due to Low Oxygen 174 Non-Lethal Effects due to Low Oxygen 180 5.6 Inorganic Carbon 181 5.7 Modeling Dynamic pH 184 5.8 Modeling Calcium Carbonate Precipitation and Effects 187 6. INORGANIC SEDIMENTS 189 6.1 Sand Silt Clay Model 189 Deposition and Scour of Silt and Clay 191 Scour, Deposition and Transport of Sand 194 Suspended Inorganic Sediments in Standing Water 196 6.2 Multi-Layer Sediment Model 197 Suspended Inorganic Sediments 199 Inorganics in the Sediment Bed 199 Detritus in the Sediment Bed 201 Pore Waters in the Sediment Bed 201 Dissolved Organic Matter within Pore Waters 202 Diffusion within Pore Waters 203 Sediment Interactions 204 7. SEDIMENT DIAGENESIS 207 7.1 Sediment Fluxes 209 7.2 POC 212 7.3 PON 214 7.4 POP 214 7.5 Ammonia 214 7.6 Nitrate 216 7.7 Orthophosphate 217 7.8 Methane 218 7.9Sulfide 220 7.10 Biogenic Silica 221 7.11 Dissolved Silica 222 8. TOXIC ORGANIC CHEMICALS 224 8.1 lonization 231 8.2 Hydrolysis 232 8.3 Photolysis 234 8.4Microbial Degradation 236 8.5 Volatilization 237 8.6 Partition Coefficients 240 Detritus 240 Algae 243 Macrophytes 244 v ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION Invertebrates 245 Fish 245 8.7 Nonequilibrium Kinetics 246 Sorption and Desorption to Detritus 247 Bioconcentration in Macrophytes and Algae 248 Macrophytes 248 Algae 249 Bioaccumulation in Animals 252 Gill Sorption 252 Dietary Uptake 255 Elimination 256 Bioaccumulation Factor 260 Linkages to Detrital Compartments 261 8.8 Alternative Uptake Model: Entering BCFs, Kl, and K2 261 8.9 Half-Life Calculation, DT50 and DT95 262 8.10 Chemical Sorption to Sediments 263 8.11 Chemicals in Pore Waters 265 8.12 Mass Balance Capabilities and Testing 267 8.13 Perfluoroalkylated Surfactants Submodel 269 Sorption 269 Biotransformation and Other Fate Processes 269 Bioaccumulation 269 Gill Uptake 270 Dietary Assimilation 271 Depuration 272 Bioconcentration Factors 273 9. ECOTOXICOLOGY 275 9.1 Lethal Toxicity of Compounds 275 Interspecies Correlation Estimates (ICE) 275 Internal Calculations 277 9.2 Sublethal Toxicity 280 9.3 External Toxicity 283 10. ESTUARINE SUBMODEL 286 10.1 Estuarine Stratification 286 10.2 Tidal Amplitude 287 10.3 Water Balance 288 10.4 Estuarine Exchange 289 10.5 Salinity Effects 290 Mortality and Gamete Loss 290 Other Biotic Processes 290 Sinking 291 Sorption 293 Volatilization 293 Reaeration 293 Migration 295 vi ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION 10.6 Nutrient Inputs to Lower Layer 295 11. REFERENCES 296 APPENDIX A. GLOSSARY OF TERMS 315 APPENDIX B. USER SUPPLIED PARAMETERS AND DATA 318 vn ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION PREFACE The Clean Water Act- formally the Federal Water Pollution Control Act Amendments of 1972 (Public Law 92-50), and subsequent amendments in 1977, 1979, 1980, 1981, 1983, and 1987- calls for the identification, control, and prevention of pollution of the nation's waters. Data submitted by the States to the U.S. Environmental Protection Agency's WATERS (Watershed Assessment, Tracking & Environmental Results) database (http://www.epa.gov/waters/) indicate that a very high percentage of the Nations waters continue to be impaired. As of early 2009, of the waters that have been assessed, 44% of rivers and streams, 59% of lakes, reservoirs and ponds, and 35% of estuaries were impaired for one or more of their designated uses. The five most commonly reported causes of impairment in rivers and streams were: pathogens, sediment, nutrients, habitat alteration and organic enrichment/dissolved oxygen depletion. In lakes and reservoirs the five most common causes were mercury, nutrients, organic enrichment/dissolved oxygen depletion, metals, and turbidity. In estuaries the five most common causes were pathogens, mercury, organic enrichment/oxygen depletion, pesticides and toxic organics. Many waters are impaired for multiple uses, by multiple causes, from multiple sources. New approaches and tools, including appropriate technical guidance documents, are needed to facilitate ecosystem analyses of watersheds as required by the Clean Water Act. In particular, there is a pressing need for refinement and release of an ecological risk methodology that addresses the direct, indirect, and synergistic effects of nutrients, metals, toxic organic chemicals, and non-chemical stressors on aquatic ecosystems, including streams, rivers, lakes, and estuaries. The ecosystem model AQUATOX is one of the few general ecological risk models that represents the combined environmental fate and effects of toxic chemicals. The model also represents conventional pollutants, such as nutrients and sediments, and considers several trophic levels, including attached and planktonic algae, submerged aquatic vegetation, several types of invertebrates, and several types of fish. It has been implemented for experimental tanks, ponds and pond enclosures, streams, small rivers, linked river segments, lakes, reservoirs, linked reservoir segments, and estuaries. Vlll ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 1 1. INTRODUCTION 1.1 Overview The AQUATOX model is a general ecological risk assessment model that represents the combined environmental fate and effects of conventional pollutants, such as nutrients and sediments, and toxic chemicals in aquatic ecosystems. It considers several trophic levels, including attached and planktonic algae and submerged aquatic vegetation, invertebrates, and forage, bottom-feeding, and game fish; it also represents associated organic toxicants (Figure 1). It can be implemented as a simple model (indeed, it has been used to simulate an abiotic flask) or as a truly complex food-web model. Often it is desirable to model a food web rather than a food chain, for example to examine the possibility of less tolerant organisms being replaced by more tolerant organisms as environmental perturbations occur. "Food web models provide a means for validation because they mechanistically describe the bioaccumulation process and ascribe causality to observed relationships between biota and sediment or water" (Connolly and Glaser 1998). The best way to accurately assess bioaccumulation is to use more complex models, but only if the data needs of the models can be met and there is sufficient time (Pelka 1998). It has been implemented for experimental tanks, ponds and pond enclosures, streams, small rivers, linked river segments, lakes, reservoirs, linked reservoir segments, and estuaries. It is intended to be used to evaluate the likelihood of past, present, and future adverse effects from various stressors including potentially toxic organic chemicals, nutrients, organic wastes, sediments, and temperature. The stressors may be considered individually or together. The fate portion of the model, which is applicable especially to organic toxicants, includes: partitioning among organisms, suspended and sedimented detritus, suspended and sedimented inorganic sediments, and water; volatilization; hydrolysis; photolysis; ionization; and microbial degradation. The effects portion of the model includes: sublethal and lethal toxicity to the various organisms modeled; and indirect effects such as release of grazing and predation pressure, increase in detritus and recycling of nutrients from killed organisms, dissolved oxygen sag due to increased decomposition, and loss of food base for animals. AQUATOX represents the aquatic ecosystem by simulating the changing concentrations (in mg/L or g/m3) of organisms, nutrients, chemicals, and sediments in a unit volume of water (Figure 1). As such, it differs from population models, which represent the changes in numbers of individuals. As O'Neill et al. (1986) stated, ecosystem models and population models are complementary; one cannot take the place of the other. Population models excel at modeling individual species at risk and modeling fishing pressure and other age/size-specific aspects; but recycling of nutrients, the combined fate and effects of toxic chemicals, and other interdependences in the aquatic ecosystem are important aspects that AQUATOX represents and that cannot be addressed by a population model. ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 1 Figure 1. Conceptual model of ecosystem represented by AQUATOX Suspended and bedded sedim ents Settling, scour Plants Phytoplankton Attached algae M acrophytes o *= O. Any ecosystem model consists of multiple components requiring input data. These are the abiotic and biotic state variables or compartments being simulated (Figure 2). In AQUATOX the biotic state variables may represent trophic levels, guilds, and/or species. The model can represent a food web with both detrital- and algal-based trophic linkages. Closely related are driving variables, such as temperature, light, and nutrient loadings, which force the system to behave in certain ways. In AQUATOX state variables and driving variables are treated similarly in the code. This provides flexibility because external loadings of state variables, such as phytoplankton carried into a reach from upstream, may function as driving variables; and driving variables, such as temperature, could be treated as dynamic state variables in a future implementation. Constant, dynamic, and multiplicative loadings can be specified for atmospheric, point- and nonpoint sources. Loadings of pollutants can be turned off at the click of a button to obtain a control simulation for comparison with the perturbed simulation. ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 1 Figure 2. State variables in AQUATOX as implemented for Cahaba River, Alabama. Forage Fish shiner, bluegill Bottom Fish stoneroller Piscivore bass Zoobenthos grazers mayfly, riffle beetle Zoobenthos susp. feeders caddisfly Predatory Invertebrate crayfish, stonefly Zoobenthos molluscs snail, mussel Corbicula Penphyton diatom, green Phytoplankton diatom, green Macrophyte moss Periphyton blue-green Nitrate & Nitrite Carbon Dioxide Labile Diss. Detritus Refractory Susp. Detritus Labile Susp. Detritus Refractory Diss. Detritus Labile Sed. Detritus Buried Refrac Sed. Detritus Refractory Sed. Detritus Sediments sand, silt, clay (or TSS) The model is written in object-oriented Pascal using the Delphi programming system for Windows. An object is a unit of computer code that can be duplicated; its characteristics and methods also can be inherited by higher-level objects. For example, the organism object, including variables such as the LC50 (lethal concentration of a toxicant) and process functions such as respiration, is inherited by the plant object; that is enhanced by plant-specific variables and functions and is duplicated for four kinds of algae; and the plant object is inherited and modified slightly for macrophytes and moss. This modularity forms the basis for the remarkable flexibility of the model, including the ability to add and delete given state variables interactively. AQUATOX utilizes differential equations to represent changing values of state variables, normally with a reporting time step of one day. These equations require starting values or initial conditions for the beginning of the simulation. If the first day of a simulation is changed, then the initial conditions may need to be changed. A simulation can begin with any date and may be for any length of time from a few days, corresponding to a microcosm experiment, to decades, corresponding to an extreme event followed by long-term recovery. The process equations contain another class of input variables: the parameters or coefficients that allow the user to specify key process characteristics. For example, the maximum consumption rate is a critical parameter characterizing various consumers. AQUATOX is a ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 1 mechanistic model with many parameters; however, default values are available so that the analyst only has to be concerned with those parameters necessary for a specific risk analysis, such as characterization of a new chemical. In the pages that follow, differential equations for the state variables will be followed by process equations and parameter definitions. Finally, the system being modeled is characterized by site constants, such as mean and maximum depths. At present one can model lakes, reservoirs, streams, small rivers, estuaries, and ponds- and even enclosures and tanks. The "Generalized Parameter" screen is used for all these site types, although some, such as the hypolimnion and estuary entries, obviously are not applicable to all. The temperature and light constants are used for simple forcing functions, blurring the distinctions between site constants and driving variables. Table 1. Model Overview Summary (also see Section 2.1) Category: Reporting Time Step Differentiation Output Averaging Conceptual Approach Horizontal Spatial Resolution Vertical Spatial Resolution Sediment Bed Boundary Conditions Ecological Complexity Chemical Complexity Mass Balance Tracking Summary: Daily or Hourly Variable time-step Runge Kutta (with fixed time-step option) Variable Kinetic; biomass model Point model, or ID and 2D with linked segments Vertically stratified water column when relevant Multiple sediment bed options Inflows and outflows of all state variables (dissolved oxygen, nutrients, biota, detritus, and toxic organics) Variable user can model representative groups or individual species Zero to 20 organic chemicals For nutrients and chemicals Notes: time-step over which equations are solved smaller step sizes than the reporting time- step may be utilized to reduce relative error editable by user no longer a fugacity option for chemicals; individual organisms are not modeled modeled units can be a lake, river, reservoir, stream segment, estuary, or enclosure user-specified or model calculated dates of stratification active layer only, multi-layer sediments, sediment diagenesis submodels water inflow, point sources, nonpoint sources, direct precipitation, separate tributary inputs can model abiotic conditions or single macrophyte species in a water tank up to dozens of plant and animal species in a complex river or reservoir system biotransformation to daughter products may be modeled 1.2 Background AQUATOX 3.1 is an update to AQUATOX Release 3 (see section 1.6 below for a list of updates in Release 3.1). AQUATOX Release 3 was the result of an effort to combine all of the various versions of AQUATOX into a single consolidated version. Models that were combined to produce Release 3 included: AQUATOX Multi-Segment version AQUATOX Estuarine Version ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 1 . AQUATOX PFA Model (Perfloroalkylated Surfactants) Each of these versions is discussed in a separate section below. AQUATOX is the latest in a long series of models, starting with the aquatic ecosystem model CLEAN (Park et al., 1974) and subsequently improved in consultation with numerous researchers at various European hydrobiological laboratories, resulting in the CLEANER series (Park et al., 1975, 1979, 1980; Park, 1978; Scavia and Park, 1976) and LAKETRACE (Collins and Park, 1989). The MACROPHYTE model, developed for the U.S. Army Corps of Engineers (Collins et al., 1985), provided additional capability for representing submersed aquatic vegetation. Another series started with the toxic fate model PEST, developed to complement CLEANER (Park et al., 1980, 1982), and continued with the TOXTRACE model (Park, 1984) and the spreadsheet equilibrium fugacity PART model. AQUATOX combined algorithms from these models with ecotoxicological constructs; and additional code was written as required for a truly integrative fate and effects model (Park, 1990, 1993). The model was then restructured and linked to Microsoft Windows interfaces to provide greater flexibility, capacity for additional compartments, and user friendliness (Park et al., 1995). The current version has been improved with the addition of constructs for sublethal effects and uncertainty analysis, making it a powerful tool for probabilistic risk assessment. This technical documentation is intended to provide verification of individual constructs or mathematical and programming formulations used within AQUATOX. The scientific basis of the constructs reflects empirical and theoretical support; precedence in the open literature and in widely used models is noted. Units are given to confirm the dimensional analysis. The mathematical formulations have been programmed and graphed in spreadsheets and the results have been evaluated in terms of behavior consistent with our understanding of ecosystem response; many of those graphs are given in the following documentation. The variable names in the documentation correspond to those used in the program so that the mathematical formulations and code can be compared, and the computer code has been checked for consistency with those formulations. Much of this has been done as part of the continuing process of internal review. Releases 2 and 3 of the AQUATOX model and documentation have undergone successful peer reviews by external panels convened by the U.S. Environmental Protection Agency. Release 3 has also been described in the peer-reviewed literature (Park et al. 2008). Release 3 has significant additional capabilities compared to Release 2.2: Link to WEB-ICE (Interspecies Correlation Estimates) database and graphics Sediment diagenesis based on the Di Toro model Optional hourly time step with diel oxygen, light, and photosynthesis; Low oxygen effects Toxicity due to ammonia Suspended and bedded sediment effects on organisms; % embeddedness Calcium carbonate precipitation and removal of phosphorus Adaptive light limitation for plants ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 1 Linked periphyton and phytoplankton compartments Conversions for many units in input screens User-specified seasonally varying thermocline depth User-specified reaeration constant in addition to alternative estimation procedures Improved CBOD to organic matter estimation Estuarine reaeration incorporating salinity Sensitivity analysis with tornado diagrams Correlation of variables in uncertainty analysis Sediment oxygen demand (SOD) output Enhanced graphics including log plot, duration and exceedance graphs, and threshold analysis Option to export all graphs to Microsoft Word Output of statistics for all graphed model results Output of trophic state indices and ecosystem bioenergetics such as gross primary productivity and community respiration Integrated user's manual and context-sensitive help files 1.3 The Multi-Segment Version The AQUATOX Multi-Segment version was developed and applied for the EPA Office of Water in support of the Modeling Study of PCB Contamination in the Housatonic River. Capabilities introduced with this version include the linkage of individual AQUATOX segments into a single simulation. Segments can be linked together in a manner that allows feedback into the upstream segment or a one-way "cascade" linkage can be created. More information about the physical characteristics of linked segments may be found in Section 3.8 of this document. Additionally, a sediment submodel was added to the AQUATOX model to enable tracing the passage of toxicants within a multi-layered sediment bed. Specifications for this multi-layer sediment model may be found in section 6.2 of this document. 1.4 The Estuarine Submodel The Risk Assessment Division (RAD), EPA Office of Pollution Prevention and Toxics, is responsible for assessing the human health and ecological risks of new and existing chemicals that are regulated under the Toxic Substances Control Act (TSCA). RAD has partially funded AQUATOX from its initial conceptualization. Many of the industrial chemicals regulated under TSCA are discharged into estuarine environments. Therefore, AQUATOX's capabilities were enhanced by adding salinity and other components (including shore birds) that would be needed to simulate an estuarine environment. The estuarine version of AQUATOX is intended to be an exploratory model for evaluating the possible fate ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 1 and effects of toxic chemicals and other pollutants in estuarine ecosystems. The model is not intended to represent detailed, spatially varying site-specific conditions, but rather to be used in representing the potential behavior of chemicals under average conditions. Therefore, it is best used as a screening-level model applicable to data-poor evaluations in estuarine ecosystems. Complete documentation for the AQUATOX estuarine submodel may be found in Chapter 10 of this document. 1.5 The PFA Submodel The bioaccumulation and effects of a group of chemicals known as perfluorinated surfactants has been of recent interest. There are two major types of perfluorinated surfactants: perfluoro- alkanesulfonates and perfluorocarboxylates. The perfluorinated compounds of interest as bioaccumulators are the perfluorinated acids (PFAs). Perfluoroctane sulfonate (PFOS) belongs to the sulfonate group and perfluorooctanoic acid (PFOA) belongs to the carboxylate group. These persistent chemicals have been found in humans, fish, birds, and marine and terrestrial mammals throughout the world. PFOS has an especially high bioconcentration factor in fish. The principal focus was on PFOS because of its prevalence and the availability of data. Because both chemical classes contain high- and low-chain homologs, AQUATOX will be useful in estimating the fate and effects of a wide range of molecular weight components where actual data are not available for every homolog. Complete documentation for the AQUATOX Perfloroalkylated Surfactants model may be found in Section 8.13 of this document. 1.6 AQUATOX Release 3.1 Overview Additional capabilities are available in Release 3.1 as compared to Release 3. Some highlights follow: Addition of sediment-diagenesis "steady-state" mode to significantly increase model speed; Modification of denitrification code in order to simplify calibration and to achieve alignment with other models; Enabled importation of equilibrium CO2 concentrations to enable linkage to CO2SYS and similar models; New CBOD to organic matter conversion relying on percent-refractory detritus input; Input and output BOD is clarified to be "carbonaceous" BOD. Floating plants refinements Added floating option for plants other than cyanobacteria (formerly known as "blue-green algae) Converted the averaging depth for floating plants to the top three meters to more closely correspond to monitoring data Floating plants now explicitly move from the hypolimnion to the epilimnion when 7 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 1 a system is stratified. Modifications to PFA (perfluoroalkylated surfactants) model to increase flexibility: Uptake rates (Kls) and elimination rates (K2s) are visible and editable for animals and plants New interface to estimate animal Kls and K2s as a function of chain length Improved gill-uptake equation for invertebrates. Bioaccumulation and toxicity modeling improvements: Optional alternative elimination-rate estimation for animals based on Barber (2003); Updated ICE (toxicity regressions), based on new EPA models released in February 2010 and improved AQUATOX ICE interface; Addition of output of Kl, K2, and BCF estimates. Improved sensitivity and uncertainty analyses "Output to CSV" option for uncertainty runs so that complete results for every iteration may be examined; Option for non-random sampling for "statistical sensitivity analyses"; A "reverse tornado" diagram (effects diagram) that shows the effects of each parameter change on the overall simulation; Nominal range sensitivity analysis has been added for linked segment applications. Database Improvements AQUATOX database search functions dramatically improved. "Scientific Name" field added to Animal and Plant databases. Interface and Data Input Improvements Software and software installer is 64-bit OS compatible; Added an option in the "Setup" screen to trigger nitrogen fixation based on the N to P ratio. Addition of output variables to clarify whether photosynthesis is sub-optimal due to high-light or low-light conditions. Time-varying evaporation option in the "Site" screen with linkage from the "Water Volume" screen Grid mode within a study so that all animal, plant, and chemical parameters in a study can be examined, edited, and then exported to Excel Added capability to input time-series loads of organisms based on fish stocking Updated HSPF WDM file linkage to be more generally applicable (does not require use of WinHSPF). Enabled hourly loadings for the following variables: all nutrients, CO2, Oxygen, Inorganic suspended sediments (sand/silt/clay), TSS, Light, Organic Matter "Graph Setup" window now enabled for linked-mode graphics. Other minor interface improvements. ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 1 Documentation for each of these enhancements may be found in this technical documentation volume or in the User's Manual. 1.7 Comparison with Other Models The following comparison is taken from Park et al. (2008): The model is perhaps the most comprehensive aquatic ecosystem simulation model available, as can be seen by comparison with other representative dynamic models being used for risk assessment (Table 2). All the models, with the exception of QSim and CASM, are public domain. The closest to AQUATOX in terms of scope is the family of CATS models developed by Traas and others (Traas et al., 1996; Traas et al., 1998; Traas et al., 2001); these ecotoxicology models have simple representations of growth and are not as suitable as AQUATOX for detailed analyses of eutrophication effects. CASM (DeAngelis et al., 1989; Bartell et al., 1999) is similar to CATS, with simplified growth terms, but it lacks a toxicant fate component. QUAL2K (Chapra et al., 2007) and WASP (Di Toro et al., 1983; Wool et al., 2004) are water quality models that share many functions with AQUATOX, including benthic algae (Martin et al., 2006); WASP also models fate of toxicants. The hydraulic and water quality models EFDC (Tetra Tech Inc., 2002) and HEM3D (Park et al., 1995a) are often combined; EFDC has also been used to provide the flow field for linked segments in AQUATOX, resulting in a similar representation. AQUATOX, QUAL2K, WASP, and EFDC include the sediment diagenesis model for remineralization (Di Toro, 2001). WASP and the bioaccumulation model QEAFdChn (Quantitative Environmental Analysis, 2001) have been combined in the Green Bay Mass Balance (GBMB) study (U.S. Environmental Protection Agency, 1989), which Koelmans et al. (2001) considered to be more accurate for portraying bioaccumulation than AQUATOX. However, GBMB does not include an ecotoxicology component. BASS (Barber, 2001) is a very detailed bioaccumulation and ecotoxicology model; it provides better resolution than AQUATOX in modeling single species, but so far it has only been applied to fish and does not include ecosystem dynamics. The German model QSim (Schol et al., 1999; Schol et al., 2002; see also Rode et al., 2007) has detailed ecosystem functions and has been applied in studying impacts of both eutrophication and hydraulics on river ecosystems. Similar to AQUATOX, it has been used to analyze relationships between plankton and mussels and impacts of oxygen depletion. Further comparison of models can be found in a book by Pastorok et al. (2002). ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 1 Table 2. Comparison of AQUATOX with other representative dynamic models used for risk assessment (Park et al. 2008). State variables and processes. Nutrients Sediment diagenesis Detritus Dissolved sicyger. DC effects on biota PH Ntti toxitity Sand'silt/day Sediment effects Hydraulics Heat budget salinity Phytopiaritton Periphyton Macrophytes ZoaplanJtton Zoobenthos Fish Bacteria Pathogens Organic toxicant late Organic toxicants in Sediments Stratified sediments Phyttiplanktan Periphytosn Macropbytes Zooplafikton Zoobenthos Fish Birds oc other animals Eoatoxicity Linked segments AQUATOX X X X X X X X. X X X X X X X X X X, X X X X X X X X X X X CATS CASM Qual2K WASP7 EFDC-HEM3D QEAFdChn XXIX X XX X XXIX X XXX X I X X X XX X X X X X X X X XXX X X 3C X X X X X X X I X X X X XX X X X X X X X X X X I X X X XXX X BASS QSim X J. X X X X X X X X X X X X X X X X X 1.8 Intended Application of AQUATOX AQUATOX is intended to be used at any one of several levels of application. Like any model, it is best used as one of several tools in a weight-of-evidence approach. The level of required precision, rigor, data requirements and user effort depend upon the goals of the modeling exercise and the potential consequences of the model results. Perhaps its most widespread use is as a screening-level model requiring few changes to default studies and parameters. In fact, it was originally developed as an evaluative model to assess the fate and effects of pesticides and industrial organic chemicals in representative or "canonical" environments; these include ponds and pond enclosures, experimental streams, and a representative estuary. It is especially useful in taking the place of expensive, labor-intensive mesocosm tests. It has been calibrated and validated with data from pond enclosures, experimental streams, and a polluted harbor. In one early application, AQUATOX was driven with predicted pesticide runoff into a farm pond adjacent to a corn field using the field model PRZM. Also, with little effort the model can provide insights into the potential impacts of invasive species and the possible effects of control measures, such as pesticide application, on the aquatic ecosystem. In recent years AQUATOX has been applied as part of the process of developing water quality 10 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 1 targets for nutrients, and comparing model-derived values with regional criteria developed empirically. This application has involved setting up the model and calibrating with available data for rivers and reservoirs receiving nutrients from wastewater treatment plants, agricultural runoff, and background "natural" loadings. It has been our experience that this entails a substantial level of effort, especially if the system is spatially heterogeneous, which then requires application of linked segments. A certain amount of site-specific biotic, water quality and flow data is required, as well as pollutant loading data, for calibration. However, once the model is set up and calibrated for a site, it is relatively easy to represent a series of loading scenarios and determine threshold nutrient levels for deleterious impacts such as nuisance algal blooms and anoxia. This process is facilitated by the fact that the model has been calibrated across nutrient, turbidity, and discharge gradients, resulting in robust parameter sets that span these conditions. This is important because the intent of setting water quality targets is to model ecological communities under changing conditions as a result of environmental management decisions; this would give better assurance that the sometimes costly nutrient reduction actions would render the desired environmental result. The most intensive, time-consuming application of AQUATOX is in environmental remediation projects, such as SUPERFUND. Because of the likely litigation and the potential for costly remediation, this level of application requires site-specific calibration and validation using quality-assured data collected specifically for the model. In dynamic systems, linkage to an equally well calibrated and validated hydrodynamic model is essential to represent, for example, burial and exhumation of contaminated sediments. Several of the more powerful features of the model, such as the linked segments and IPX layered-sediment submodel, were developed for this type of application. Unfortunately, the one remediation application performed by the model developers cannot be published because of continuing litigation. 11 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 2. SIMULATION MODELING 2.1 Temporal and Spatial Resolution and Numerical Stability AQUATOX Release 3 is designed to be a general, realistic model of the fate and effects of pollutants in aquatic ecosystems. In order to be fast, easy to use, and verifiable, it was originally designed with the simplest spatial and temporal resolutions consistent with this objective. Releases 1X?.. ... ... , . r J Model is run with a daily or hourly Simulation Modeling: Simplifying Assumptions Each modeled segment is well- maximum time-step. Results are trapezoidally integrated may still be run as a non-dimensional point model. However, unlike previous versions of AQUATOX, in Release 3 spatial segments may be linked together to form a two- or three- dimensional model if a more complicated spatial resolution is desired. The model generally represents average daily conditions for a well-mixed aquatic system. Each segment in a multi-dimensional run is also assumed to be well-mixed in each time-step. AQUATOX also represents one-dimensional vertical epilimnetic and hypolimnetic conditions for those systems that exhibit stratification on a seasonal basis. Multi-segment systems also can be set up with vertical stratification. Furthermore, the effects of run, riffle, and pool environments can be represented for streams. Results may be plotted in the AQUATOX output screen with the capability to import observed data to examine against model predictions. While the model is generally run with a daily maximum time-step, the temporal resolution of the model can also be reduced to an hourly maximum time-step. This capability was added so that AQUATOX can represent diel oxygen. See sections 3.6 and 5.5 for more information on how this choice of hourly time-step affects AQUATOX equations. The reporting step can be as long as several years or as short as one hour; results are integrated to obtain the desired reporting time period. According to Ford and Thornton (1979), a one-dimensional model is appropriate for reservoirs that are between 0.5 and 10 km in length; if larger, then a two-dimensional model disaggregated along the long axis is indicated. The one-dimensional assumption is also appropriate for many lakes (Stefan and Fang, 1994). Similarly, one can consider a single reach or stretch of river at a time. Usually the reporting time step is one day, but numerical instability is avoided by allowing the step size of the integration to vary to achieve a predetermined accuracy in the solution. (This is a numerical approach, and the step size is not directly related to the temporal scale of the ecosystem simulation.) AQUATOX uses a very efficient fourth- and fifth-order Runge-Kutta integration routine with adaptive step size to solve the differential equations (Press et al., 1986, 1992). The routine uses the fifth-order solution to determine the error associated with the fourth- order solution; it decreases the step size (often to 15 minutes or less) when rapid changes occur and increases the step size when there are slow changes, such as in winter. However, the step size is constrained to a maximum of one day (or one hour in hourly simulations) so that short- term pollutant loadings are always detected. The reporting step, on the other hand, can be as 12 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 long as several years or as short as one hour; results are integrated to obtain the desired reporting time period. As an alternative, the user may specify an exact step size that is used thoughout the simulation. This is similar to the way that many models solve the differential equations. The disadvantage is that the accuracy of the solution may not be maintained. However, it is useful under some circumstances and is discussed more fully later in this section. The temporal and spatial resolution is in keeping with the generality and realism of the model (see Park and Collins, 1982). Careful consideration has been given to the hierarchical nature of the system. Hierarchy theory tells us that models should have resolutions appropriate to the objectives; phenomena with temporal and spatial scales that are significantly longer than those of interest should be treated as constants, and phenomena with much smaller temporal and spatial scales should be treated as steady-state properties or parameters (Figure 3; also see O'Neill et al., 1986). AQUATOX uses a longer time step than dynamic hydrologic models that are concerned with representing short-term phenomena such as storm hydrographs, and it uses a shorter time step than fate models that may be concerned only with long-term patterns such as bioaccumulation in large fish. Figure 3. Position of ecosystem models such as AQUATOX in the spatial-temporal hierarchy of models. Rule-based habitat models succession, urbanization, sea-level rise Ecosystem models Population models High-resolution process models flood hydrograph diurnal pH Changing the permissible relative error (the difference between the fourth- and fifth-order solutions) of the simulation can affect the results. The model allows the user to set the relative error, usually between 0.005 and 0.01. Comparison of output shows that up to a point a smaller error can yield a marked improvement in the simulation, although execution time is longer. For example, simulations of two pulsed doses of chlorpyrifos in a pond exhibit a spread in the first pulse of about 0.6 ug/L dissolved toxicant between the simulation with 0.001 relative error and the simulation with 0.05 relative error (Figure 4); this is probably due in part to differences in the timing of the reporting step. However, if we examine the dissolved oxygen levels, which combine the effects of photosynthesis, decomposition, and reaeration, we find that there are pronounced differences over the entire simulation period. The simulations with 0.001 and 0.01 relative error give almost exactly the same results, suggesting that the more efficient 0.01 relative error should be used; the simulation with 0.05 relative error exhibits instability in the oxygen 13 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 simulation; and the simulation with 0.1 error gives quite different values for dissolved oxygen (Figure 5). The observed mean daily maximum dissolved oxygen for that period was 9.2 mg/L (U.S. Environmental Protection Agency, 1988), which corresponds most closely with the results of simulation with 0.001 and 0.01 relative error. Figure 4. Pond with chlorpyrifos in dissolved phase. Figure 5. Same as Figure 4 with Dissolved Oxygen. 0) 06/19/88 06/30/88 07/12/88 06/24/88 07/06/88 07/18/88 0.001 -0.01 0.05 12 r 11 10 o O 06/19/88 06/30/88 07/12/88 06/24/88 07/06/88 07/18/88 0.001 0.01 0.05 0.1 A common use of AQUATOX is to determine the impact of a perturbation in a perturbed simulation when compared with a control simulation. For example, the model is often run with and without a potentially toxic organic chemical, and a percent difference graph is plotted showing how the two simulations differ. Of particular interest is whether there are likely to be significant differences in state variables or other environmental indices at very low concentrations of the chemical. Because a simulation with the toxicant may require a decrease in step size to capture the dynamics of the fate of the toxicant as opposed to a simulation without the toxicant, there may be a mismatch in the step sizes of the two simulations, and the simulations may differ solely on the basis of the difference in numerical resolution. Although decreasing the relative error may decrease the mismatch, there may still be a difference that prevents determination of the "no effects" level of the chemical. In the example that follows, toxicity of PFOS has been turned off by setting all LC50 and EC50 parameter values to 0. In Figure 6 the default variable step size option has been used with a very small relative error. In Figure 7 the simulation is the same except the constant step size option has been used, and it is readily apparent that there is no difference between the perturbed and control runs. 14 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 Figure 6. Percent differences in fish biomass between perturbed simulation with PFOS (toxicity turned off) and control simulation without PFOS, using variable step size (relative error = 0.0001). The differences are purely artifacts of the numerical method. 50.0 40.0 30.0 20.0 LU £ 10.0 LU aL t °-° J -10.0 -20.0 -30.0 -40.0 -50.0 Conasauga River, GA (Difference) "-N / ^~J Stoneroller Spotted Bass Catfish Shiner Bluegill 3/3/2008 6/1/2008 8/30/2008 11/28/2008 2/26/2009 Figure 7. There is no difference in fish biomass between perturbed simulation with PFOS (toxicity turned off) and control simulation without PFOS using a constant step size (0.1 day). 5.00E+01 4.00E+01 3.00E+01 2.00E+01 LU g 1.00E+01 oL LL O.OOE+00 LL ° -1.00E+01 5* -2.00E+01 -3.00E+01 -4.00E+01 -5.00E+01 Conasauga River, GA (Difference) Stoneroller Spotted Bass Catfish Shiner Bluegill 3/3/2008 6/1/2008 8/30/2008 11/28/2008 2/26/2009 2.2 Results Reporting The AQUATOX results reporting time step may be set to any desired frequency, from a fraction of an hour to multiple years. The Runge-Kutta differential equations solver produces a series of results of variable frequency; this frequency may be either greater than or less than the reporting time-step. To standardize AQUATOX output, the user has two options, the trapezoidal 15 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 integration of results (default) or the output of "instantaneous" concentrations. Using either of these options, AQUATOX will produce output with time-stamps that match the reporting time- step precisely. When instantaneous concentrations are requested (in the model's setup screen) AQUATOX returns output precisely at the requested reporting time-step through linear interpolation of the nearest Runge-Kutta results that occur before and after the relevant reporting time-step. When results are trapezoidally integrated, AQUATOX calculates results by summing all of the trapezoids that can be produced by linear interpolation between Runge-Kutta results and dividing by the results-reporting step-size to get an average result over the reporting step. In Figure 8, for example, the areas of the four shaded trapezoids are summed together and this sum is divided by the results reporting step to achieve an average result over that reporting step. When trapezoidal integration is selected, AQUATOX output is time-stamped at the end of the interval over which the integration is taking place. For example, if a user selects a 366.25 day time-step, the results at the end of the first year will be reflective of all time-steps calculated within that year. Figure 8. An example of trapezoidal integration. Runge Kutta Step Size (Variable) Results- Reporting Step Results may be plotted in the AQUATOX output screen including the capability to import observed data to examine against model predictions. 2.3 Input Data AQUATOX accepts several forms of input data, a partial list of which follows: Point-estimate parameters describing animals, plants, chemicals, sites, and remineralization. Default values for these parameters are generally available from included databases (called "libraries"). The full list of these parameters, their units, and 16 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 their manner of reference in the interface, this document, and the source code may be found in Appendix B of this document. Time series (or constant values) for nutrient-inflow, organic matter-inflow, and gas- inflow loadings. Time series for inorganic sediments in water, water volume variables, and the pH, light, and temperature climates. Time series of chemical inflow loadings and initial conditions. A feeding preference matrix must be specified to describe the food web in the simulation. Additional parameters may be required depending on which submodels are included (e.g. additional sediment diagenesis parameters.) Nearly all point-estimate parameters may be represented by distributions when the model is run in uncertainty mode (see section 2.5). For more discussion of AQUATOX data requirements please see the "Data Requirements" section in the AQUATOX Users Manual (or in the context sensitive help files included with the model software). Furthermore, a Technical Note on Data Requirements is available. For time-series loadings, when a value is input for every day of a simulation, AQUATOX will read the relevant value on each day. If missing values are encountered by the model, a linear interpolation will be performed between the surrounding dates. If the AQUATOX simulation time includes dates before or after the input time-series the model assumes an annual cycle and tries to calculate the appropriate input value accordingly. Please see the "Important Note about Dynamic Loadings" in the AQUATOX Users Manual (integrated help-file) for a complete description of this process. 2.4 Sensitivity Analysis "Sensitivity" refers to the variation in output of a mathematical model with respect to changes in the values of the model inputs (Saltelli 2001). It provides a ranking of the model input assumptions with respect to their relative contribution to model output variability or uncertainty (U.S. Environmental Protection Agency 1997). Simplifying Assumptions: Parameters are treated as independent Feeding preference matrices are not included Sensitivity is compared for the last step of the simulation Caution 10% change is appropriate, a large change can exceed reasonable values and give misleading results AQUATOX includes a built-in nominal range sensitivity analysis (Frey and Patil 2001), which may be used to examine the sensitivity of multiple model outputs to multiple model parameters. The user first selects which model parameters to vary and which output variables to track. The model iteratively steps through each of the parameters and varies them by a given percent in the positive and negative direction and saves model results in an Excel file. A sensitivity statistic may then be calculated such that when a 10% change in the parameter results in a 10% change in the model result, the sensitivity is calculated as 100%. 17 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 Sensitivity = Result Pos-Result Basetine + ResultNeg-ResultBaselme 100 2- Result, 'aseline PctChanged where: Sensitivity J\eSUltScenario PctChanged normalized sensitivity statistic (%); averaged AQUATOX result for a given endpoint given a positive change in the input parameter, a negative change in the input parameter or no change in the input parameter (baseline) percent that the input parameter is modified in the positive and negative directions. Sensitivity is computed for the last time step of the simulation, so one usually sets the reporting time step to encompass a year or the entire period of the simulation. For each output variable tracked, model parameters may be sorted on the average sensitivity (for the positive and negative tests) and plotted on a bar chart. The end result is referred to as a "Tornado Diagram." Tornado diagrams may automatically be produced within the AQUATOX output window (Figure 9). When interpreting a tornado diagram, the vertical line at the middle of the diagram represents the deterministic model result. Red lines represent model results when the given parameter is reduced by the user-input percentage while blue lines represent a positive change in the parameter. An "effects diagram" that illustrates the effects of a single parameter change on all tracked outputs can also be created. See the User's Manual (or context-sensitive help) for more information on how to create and interpret these types of output. When sensitivity analysis is run on a multi-segment model, the user must choose either parameters that are relevant to all segments (e.g. animal or plant parameters) or individual segments (e.g. state-variable initial conditions, or the segment's physical characteristics). The segment for which the parameter is relevant may be selected in the sensitivity analysis setup window (see the User's Manual for more information) Any number of global or segment- specific parameters may be selected for a single sensitivity-analysis; output files will be written for each segment in the simulation. 18 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 Figure 9. An example tornado diagram showing calculated sensitivity statistic. Sensitivity of Chironomid (g/m2 dry) to a 20% change in the 15 most sensitive (tested) parameters 220% - Site: Ave. Epilimnetic Temperature (deg C) 163% - Chironomid: Maximum Temperature (deg. C) 137% - Peri, Green: Max Photosynthetic Rate (1/d) - 1 1 5% - Temp: Multiply Loading by I ' =^i ^^^ ^^^^^ i I 1 1 ^^^~ ^^^ 1.5 2 2.5 3 3.5 4 4.5 5 Chironomid (g/m2 dry) 2.5 Uncertainty Analysis There are numerous sources of uncertainty and variation in natural systems. These include: site characteristics such as water depth, which may vary seasonally and from site to site; environmental loadings such as water flow, temperature, and light, which may have a stochastic component; and critical biotic parameters such as maximum photosynthetic and consumption rates, which vary among experiments and representative organisms. In addition, there are sources of uncertainty and variation with regard to pollutants, including: pollutant loadings from runoff, point sources, and atmospheric deposition, which may vary stochastically from day to day and year to year; physico-chemical characteristics such as octanol-water partition coefficients and Henry Law constants that cannot be measured easily; chemodynamic parameters such as microbial degradation, photolysis, and hydrolysis rates, which may be subject to both measurement errors and indeterminate environmental controls. Increasingly, environmental analysts and decision makers are requiring probabilistic modeling approaches so that they can consider the implications of uncertainty in the analyses. AQUATOX Uncertainty Analysis: Strengths Use of Latin hypercube sampling is more efficient than brute-force Monte Carlo analysis Nearly all variables and parameters may be represented as distributions Variables can be correlated Simplifying Assumptions: Feeding preference matrices are not included Modeled correlations cannot be perfect (e.g. 1.0) due to limitations of the Iman & Conover method 19 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 provides this capability by allowing the user to specify the types of distributions and key statistics for almost all input variables. Depending on the specific variable and the amount of available information, any one of several distributions may be most appropriate. A lognormal distribution is the default for environmental and pollutant loadings. In the uncertainty analysis, the distributions for constant loadings are sampled daily, providing day-to-day variation within the limits of the distribution, reflecting the stochastic nature of such loadings. A useful tool in testing scenarios is the multiplicative loading factor, which can be applied to all loads. Distributions for dynamic loadings may employ multiplicative factors that are sampled once each iteration (Figure 10). Normally the multiplicative factor for a loading is set to 1, but, as seen in the example, under extreme conditions the loading may be ten times as great. In this way the user could represent unexpected conditions such as pesticides being applied inadvertently just before each large storm of the season. Loadings usually exhibit a lognormal distribution, and that is suggested in these applications, unless there is information to the contrary. Figure 11 exhibits the result of such a loading distribution. Figure 10. Distribution screen for point-source loading of toxicant in water. Distribution Type: i Triangular Uniform C Normal '' Lognormal 0.673 5.83 Distribution Parameters: Mean [i Stcl. Deviation fo.6 ' Probability r Cumulative Distribution In an Uncertainty Run: <' Use Above Distribution Use Point Estimate Choice of distribution: A sequence of increasingly informative distributions should be considered for most parameters. If only two values are known and nothing more can be assumed, the two values may be used as minimum and maximum values for a uniform distribution (Figure 12); this is often used for parameters where only two values are known. If minimal information is available but there is reason to accept a particular value as most likely, perhaps based on calibration, then a triangular distribution may be most suitable (Figure 13). Note that the minimum and maximum values for the distribution are constraints that have zero probability of occurrence. If additional data are available indicating both a central tendency and spread of response, such as parameters for well-studied processes, then a normal distribution 20 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 may be most appropriate (Figure 14). The result of applying such a distribution in a simulation of Onondaga Lake, New York, is shown in Figure 15, where simulated benthic feeding affects decomposition and subsequently the predicted hypolimnetic anoxia. Most distributions are truncated at zero because negative values would have no meaning (Log Kow is one exception). Figure 11. Sensitivity of bass (g/m2) to variations in loadings of dieldrin in Coralville Lake, Iowa. Largemouth Ba2 (g/m2 4/17/2009 2:16:20 PM 0.7 Mean Minimum Maximum Mean - StDev Mean + StDev Deterministic 10/24/1969 12/23/1969 2/21/1970 4/22/1970 6/21/1970 8/20/1970 Figure 12. Uniform distribution for Henry's Law constant for esfenvalerate. 0.00 6.1 E-8 1.53E-6 3E-6 Figure 13. Triangular distribution for maximum consumption rate for bass. 0.04 0.00 0.015 0.0425 0.07 21 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 Figure 14. Normal distribution for maximum consumption rate for the detritivorous invertebrate Tubifex. (Distribution Information 0.04 0.03 0.02 0.01 0.25 048.3 V Probability r Cumulative Distribution Distribution Typo: r Triangular r uniform ( Normal C Lotjnormal In an Uncertainty Run: ff Use Above Distribution r Use Point Estimate Distribution Parameters: Mean ]0.25 Stii. Deviation 10.1 V OK X Cancel Figure 15. Sensitivity of hypolimnetic oxygen in Lake Onondaga to variations in maximum consumption rates of detritivores. 01/01/89 09/24/89 06/17/90 05/14/89 02/04/90 10/28/90 Minimum Mean Maximum - - Deterministic Efficient sampling from the distributions is obtained with the Latin hypercube method (McKay et al., 1979; Palisade Corporation, 1991). Depending on how many iterations are chosen for the analysis, each cumulative distribution is subdivided into that many equal segments. Then a 22 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 uniform random value is chosen within each segment and used in one of the subsequent simulation runs. For example, the distribution shown in Figure 14 can be sampled as shown in Figure 16. This method is particularly advantageous because all regions of the distribution, including the tails, are sampled. A non-random seed can be used for the random number generator, causing the same sequence of numbers to be picked in successive applications; this is useful if you want to be able to duplicate the results exactly. The default is twenty iterations, meaning that twenty simulations will be performed with sampled input values; this should be considered the minimum number to provide any reliability. The optimal number can be determined experimentally by noting the number required to obtain convergence of mean response values for key state variables; in other words, at what point do additional iterations not result in significant changes in the results? As many variables may be represented by distributions as desired. Correlations may be imposed using the method of Iman and Conover (1982). By varying one parameter at a time the sensitivity of the model to individual parameters can be determined in a more rigorous way than nominal range sensitivity offers. This is done for key parameters in the following documentation. Figure 16. Latin hypercube sampling of a cumulative distribution with a mean of 25 and standard deviation of 8 divided into 5 intervals. 4 0.8 0.6 0.4 0,2 t 0 - An alternate way of presenting uncertainty is by means of a biomass risk graph, which plots the probability that biomass will be reduced by a given percentage by the end of the simulation (Mauriello and Park 2002). In practice, AQUATOX compares the end value with the initial condition for each state variable, expressing the result as a percent decline: Decline = \ 1 - EndVal StartVal \-100 (1) where: Decline EndVal percent decline in biomass for a given state variable (%); value at the end of the simulation for a given state variable (units depend on state variable); 23 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 StartVal initial condition for given state variable. The results from each iteration are sorted and plotted in a cumulative distribution so that the probability that a particular percent decline will be exceeded can be evaluated (Figure 17). Note that there are ten points in this example, one for each iteration as the consecutive segments of the distribution are sampled. Figure 17. Risk to bass from dieldrin in Coralville Reservoir, Iowa. Biomass Risk Graph 4/17/2009 2:17:42 PM r Largemouth Ba2 -40 -30 -20 -10 0 Percent Decline at Simulation End 10 Uncertainty analysis can also be used to perform statistical sensitivity analysis, which is much more powerful than the screening-level nominal range sensitivity analysis. Parameters are tested one at a time using the most appropriate distribution of observed parameter values. The time- varying and mean coefficient of variation can be calculated in an exported Excel file using the mean and standard deviation results for a particular endpoint. Examples will be published in a separate report. 2.6 Calibration and Validation Rykiel (1996) defines calibration as "the estimation and adjustment of model parameters and constants to improve the agreement between model output and a data set" while "validation is a demonstration that a model within its domain of applicability possesses a satisfactory range of accuracy consistent with the intended application of the model." A related process is verification, which is "a subjective assessment of the behavior of the model" (I0rgensen 1986). The terms are used in those ways in our applications of AQUATOX. Endpoints for comparison of model results and data should utilize available data for various ecosystem components, preferably covering nutrients, dissolved oxygen, and different trophic 24 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 levels, and toxic organics if they are being modeled. Although AQUATOX models a complete food web, often the only biotic data available are chlorophyll a values. The model converts biomass predictions to chlorophyll a values to facilitate comparison. Likewise, Secchi depth is computed from the overall extinction coefficient for comparison with observed data. Verification should consider process rates to confirm that the results were obtained for the correct reasons (Wlosinski and Collins 1985). Rate information that can be assessed for reasonableness and compared with observations includes sediment oxygen demand (SOD), the fluxes of phosphorus, nitrogen, and dissolved oxygen, and all biotic process rates. These can be presented in tabular and graphical form in AQUATOX. There are several measures of model performance that can be used for both calibrations and validations (Bartell et al. 1992, Schnoor 1996). The primary difficulty is in comparing general model behavior over long periods to observed data from a few points in time with poorly defined sample variability. Recognizing that evaluation is limited by the quantity and quality of data, stringent measures of goodness of fit are often inappropriate; therefore, we follow a weight-of- evidence approach with a sequence of increasingly rigorous tests to evaluate performance and build confidence in the model results: Reasonable behavior as demonstrated by time plots of key variablesis the model behavior reasonable based on general experience? Are the end conditions similar to the initial conditions? This is highly subjective, but when observed data are lacking or are sparse and restricted to short time periods it provides a limited reality check (Figure 18, Figure 19). Visual inspections of data points compared to model plotsdo the observations and predictions exhibit a reasonable concordance of values (Figure 20, Figure 21)? Visual inspection can also take into consideration if there is concordance given a slight shift in time. Do model curves fall within the error bands of observed data (Figure 22)? Alternatively, if there are limited replicates, how do the model curves compare with the spread of observed data? Do point observations fall within predicted model bounds obtained through uncertainty analysis? This has the limitation of being dependent on the precision of the model; the greater the model uncertainty, the greater the possibility of the data being encompassed by the error bounds (Figure 23). Regression of paired data and model resultsdoes the model produce results that are free of systematic bias? What is the correlation (R2)? See Figure 24, which corresponds to the results shown in Figure 20. Overlap between data and model distributions based on relative bias (rB) in combination with the ratio of variances (F)how much overlap is there (Figure 25)? Relative bias is a robust measure of how well central tendencies of predicted and observed results correspond; a value of 0 indicates that the means are the same (Bartell et al. 1992). The F test is the ratio of the variance of the model and the variance of the data. A value of 1 indicates that the variances are the same. 25 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 Do the observed and predicted values differ significantly based on their cumulative distributions (Figure 26)? The Kolmogorov-Smirnov statistic, a non-parametric test, can be used; however, the two datasets should represent the same time periods (for example, one should not compare predicted values over a year with observed values taken only during spring and summer). Figure 18. Predicted biomass patterns for animals in a hypothetical farm pond in Missouri. FARM POND, ESFENVAL (CONTROL) Run on 01-25-09 10:48 AM 0.88 0.80 0.72 0.64 0.56 0.48 0.40 0.32 0.24 0.16 0.08 - - Daphnia(mg/Ldry) Mayfly (Baetis (g/m2dry) - Gastropod (g/m2 dry) Shiner (g/m2dry) Largemouth Bas (g/m2dry) Largemouth Ba2 (g/m2 dry) 5/16/1994 8/14/1994 11/12/1994 2/10/1995 Figure 19. Sediment oxygen demand predicted for Lake Onondaga, using Di Toro sediment diagenesis option; this is an example of using rates for a reality check. ONONDAGA LAKE, NY (CONTROL) Run on 06-20-08 2:47 PM (Hypolimnion Segment) | SOD(gO2/m2d)| 0.3 0.0 1/12/1989 5/12/1989 9/9/1989 1/7/1990 5/7/1990 9/4/1990 1/2/1991 26 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 Figure 20. Comparison of predicted and observed (Oliver and Niemi 1988) PCB congener bioaccumulation factors in Lake Ontario lake trout. 11 10- 9 - 2 8 ? 7H 6 - 5 - 4 Lake Trout Prod Obs = 0.97 +/-1.03 * Observed Predicted 6 7 Log KOW 8 Figure 21. Predicted biomass and observed numbers of chironomid larvae in a Duluth, Minnesota, pond dosed with 6 ug/L chlorpyrifos. 0.40 0.36 0.32 0.28 ^0.24 -D IN 0.20 ^0.16 0.12 0.08 0.04 0.00 CHLORPYRIFOS 6 ug/L (PERTURBED) Run on 01-25-09 3:09 PM [ I 6/27/1986 7/27/1986 8/26/1986 rti_- -j i i n j i 1000 Obs. Chironomids (no./sample) -900 600 27 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 Figure 22. Predicted and observed benthic chlorophyll a in Cahaba River, Alabama; bars indicate one standard deviation in observed data. 1/1/01 7/20/01 2/5/02 8/24/02 Figure 23. Visual comparison of the envelope of model uncertainty, using two standard deviations for each of the nutrient loading distributions, with the observed data for chlorophyll a in Lake Onondaga, NY. 120 100 -Min Chloroph (ug/L) Mean Chloroph (ug/L) -Max Chloroph (ug/L) -Det Chloroph (ug/L) -ObsChl 28 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 Figure 24. Regression shows that the correlation between predicted and observed (Oliver and Niemi 1988) PCB congener bioaccumulation factors in Lake Ontario trout may be very good, but the slope indicates that there is systematic bias in the relationship. See Figure 20 for another presentation of these same results. 10 00 o 9 - 8 H o 7 4) 7 6 -\ LAKE ONTARIO TROUT ^=0.915 6789 Obs Log BAF 10 29 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 Figure 25. Relative bias and F test to compare means and variances of observed data and predicted results with AQUATOX. The isopleths correspond to the probability that the distributions of predicted and observed, as defined by the combination of the rB and F statistics, are similar. The isopleths assume normal distributions. Statistical Comparison of Means and Variances rB rB = 0.242, f = 0.400 predicted and observed distributionsare similar j-j *J pred obs Figure 26 Comparison of predicted and observed chlorophyll a in Lake Onondaga, New York (U.S. Environmental Protection Agency 2000). The Kolmogorov-Smirnov p statistic = 0.319, indicates that the distributions are not significantly different. 100 i 90 i ^ 80 1 u 70 i -£ 60: I 3 o 40 J 30 i 20 i 10 ^ 0 : 20 40 60 80 100 Chlorophyll a (ug/L) Observed Predicted 120 30 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 Data are often too sparse for adequate calibration at a given site. However, AQUATOX can be calibrated simultaneously across sites using an expanded state variable list representative of a range of conditions and using the same parameter set. In this way the observed biotic data can be pooled and the resulting state variable and parameter sets, being applicable to diverse sites, are assured to be robust. This is an approach that we have used on the Cahaba River, Alabama (Park et al. 2002); on three dissimilar rivers in Minnesota (Park et al. 2005); and on 13 diverse reaches on the Lower Boise River Idaho (CH2M HILL et al. 2008). The Minnesota rivers application is discussed below. Time series of driving variables for the Minnesota rivers were obtained from several sources with varying degrees of resolution and reliability. Results of watershed simulations with HSPF (Hydrologic Simulation Program- Fortran, a watershed loading model) were linked to AQUATOX, providing boundary conditions (site constants and drivers) for the Blue Earth and Crow Wing Rivers (Donigian et al. 2005). HSPF was not run for the Rum River; however, a U.S. Geological Survey (USGS) gage is located at the sample site and both daily discharge and sporadic water quality data were available from the USGS Web pages (search on "National Water Information System"). AQUATOX interpolates between points, and this feature was used to compute daily time series of nutrient concentrations from USGS National Water Information System (NWIS) observed data. Total suspended solids (TSS) are critical because the daily light climate for algae is affected. Therefore, we derived a significant relationship by regressing TSS against In-scaled discharge and used that to generate a daily time series for the Rum River (Figure 27). 31 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 Figure 27. TSS at Rum River: a) linear regression against daily flow at gage; b) resulting simulated daily time series (line), and observed values (symbols). a) 25 20 10 b) 140 120 - 100 - 60 H 40 - 20 - 400000 800000 1200000 Daily MeanFlow (m3/d) 1600000 2000000 0 01/99 05/99 08/99 12/99 04/00 08/00 12/00 After calibration we evaluated the efficacy of generating daily time series for TN using a regression of TN on discharge. The relationship is statistically significant and yielded a more realistic time series than the interpolation with sparse data that we had used (Figure 28). 32 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 However, calculation of the different limitations on photosynthesis indicates that N is not limiting in the Rum River (Figure 29), so we kept the simpler approach and did not repeat the calibration (see section 4.1 for an explanation of the reduction factor as an expression of nutrient limitation) . TP did not exhibit a statistically significant trend with discharge (R2 = 0.124) so the simple interpolation was also kept. Figure 28. TN at Rum River site: a) In-linear regression against daily flow at gage; b) interpolated TN observations (red) and time series (black) estimated from discharge regression. ,1! b) 1.8 1.6 1-1 ~ 1.2 | & * \ 0.8 £ 0.6 0.4 0.3 2.50 0.00 R2 = 0.6861 500000 1000000 1500000 Discharge (tni/d) 2000000 -EstTN -ObsTM 33 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 Figure 29. Predicted nutrient limitations for the dominant algal group in the Rum River. Note that N is not limiting. 0.9 0.8 0.7 0.6 0.5 u £o.5 0.4 0.3 0.2 0.1 Rum R. 18 MN (CONTROL) Run on 01-27-09 7:35 AM Peri High-Nut NJ.IM (frac) Peri High-Nut PO4J.IM (frac) - Peri High-Nut CO2_LIM (frac) 3/21/1999 7/19/1999 11/16/1999 3/15/2000 7/13/2000 11/10/2000 In almost all cases parameter values were chosen from ranges reported in the literature (for example, Le Cren and Lowe-McConnell 1980, Collins and Wlosinski 1983, Home and Goldman 1994, Jorgensen et al. 2000, Wetzel 2001). However, because these often are broad ranges and the model is very sensitive to some parameters, iterative calibration was necessary for a subset of parameters in AQUATOX. Conversely, some parameters have well established values and default values were used with confidence. A few parameters such as extinction coefficients and critical force for sloughing of periphyton are poorly defined or are unique to the AQUATOX formulations and were treated as "free" parameters subject to broad calibration. For example, some periphyton species are able to migrate vertically through the periphyton mat, and others have open growth forms; therefore, they could be assigned extinction coefficient values without regard to the physics of light transmission through biomass fixed in space. As noted earlier, sensitivity analysis can help determine how much attention needs to be paid to individual parameters. Sensitivity analysis of five diverse studies has shown that the model is sensitive to optimal temperature (TOpf) for algae and fish, maximum photosynthesis (PMax) for algae, % lost in periphytic sloughing, and log octanol-water partition coefficient (KOW). It is advisable to perform sensitivity analysis when the initial calibration is complete in order to identify parameters and driving variables requiring additional attention. Although not used in this application, if modeling a toxic chemical, there are several published sources (for example, Lyman et al. 1982, Verscheuren 1983, Schwarzenbach et al. 1993), and there are a couple excellent online references, including the US EPA ECOTOX site and the USDA ARS Pesticide Properties Database, which can be found with an Internet search engine. Calibration of AQUATOX for the Minnesota rivers used observed chlorophyll a as the primary target for obtaining best fits. Because there were only five to eight sestonic chlorophyll a observations in each of the two target years and only one benthic chlorophyll a observation at 34 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 each location, calibration adequacy was evaluated subjectively, based on generally expected behavior (e.g. blooms occurring during summer) and approximate concordance with observed values (in terms of both magnitude and timing), as determined through graphical comparisons of model output and data (Figure 30). The central tendencies are similar for predicted and observed distributions for all three sites, as shown by the relative bias (Figure 31). Despite the fluctuations in predicted chlorophyll a, the predicted and observed variances are similar for the Crow Wing River and Rum River simulations. Predicted periphyton sloughing events played a major role in determining the timing of chlorophyll a peaks in both simulations. The variance in predicted values is too high in the Blue Earth River simulation, where summer peak concentrations in 1999 appear to be overestimated by a factor of about two. The reason for this is not known, but may be related to inherent uncertainties in the simulated flow and TSS values, the sparseness of water chemistry sampling data, and/or limitations of model algorithms. Given the wide range in degree of enrichment among these three rivers, and the fact that the model was calibrated against all three data sets using a single set of parameters, a two-fold error during one period of the Blue Earth River simulation seems to be acceptable. The combined probability that the Blue Earth River predictions and observations have the same distribution, based on both central tendency and dispersion, is greater than 0.8. For the purpose of this analysis, we judged the calibration to be adequate for the three rivers. 35 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 Figure 30. Observed (symbols) and calibrated AQUATOX simulations (lines) of chlorophyll a in three Minnesota rivers: a) Blue Earth at mile 54, b) Rum at mile 18, c) Crow Wing at mile 72. Note the order- of-magnitude range in scale among the figures. a) 400 350 300 250 950 01/99 05/99 08/99 12/99 04/00 08/00 12/00 b) 01/99 05/99 08/99 12/99 04/00 08/00 12/00 c) 30 25 - 20 - 110 - CJ 5 - 0 01/99 05/99 08/99 12/99 04/00 08/00 12/00 36 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 Figure 31. Overlap between model and data distributions based on relative bias and ratio of variances, F; 1 = Blue Earth River, 2 = Crow Wing River, 3 = Rum River. Isopleths indicate the probability that the predicted and observed distributions are the same, assuming normality. -202 Relative Bias The calibrated algal model was also applied to three dissimilar sites on the Lower Boise River, Idaho, without modification from the Minnesota calibration. This provided additional verification of the generality of the parameter set. The three sites cover a broad range of nutrient and turbidity conditions over 90 km. Eckert is a low-nutrient, clear-water site upstream of Boise; Middleton receives wastewater treatment effluent and is a nutrient-enriched, clear-water site; and Parma is a nutrient-enriched, turbid site impacted by irrigation return flow from agricultural areas. Although the model overestimated periphyton at the Eckert site, the fit of the initial application (Figure 32) provided an excellent basis for further river-specific calibration. 37 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 Figure 32. Predicted (line) and observed (symbols) benthic chlorophyll a (a) at Eckert Road, (b) near Middleton, (c) near Parma, Lower Boise River, Idaho, using Minnesota parameter set. b) c) 2'JO 1/1/98 250 11, yG 36296 36796 9/27/00 5/16/99 9/27/00 38 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 As a limited validation, the calibrated model was applied to a site on the Cahaba River south of Birmingham, Alabama, with modifications to only two parameters, critical force for periphyton scouring and optimal temperature for algae. The Crow Wing and Rum Rivers have cobbles and boulders and are more sensitive to higher current velocities than the bedrock outcrops in the Cahaba River. Not only is the bedrock stable, it also provides abundant crevices and lee sides that are protected refuges for periphyton. For these reasons greater water velocity is expected to be required to initiate periphyton scour in the Cahaba River than in the Crow Wing and Rum Rivers, thus the critical force (Fcrit) for scour of periphyton was more than doubled in the Cahaba River simulation. Also, between Minnesota and Alabama one would expect different local ecotypes in resident algal species, with differing adaptations to temperature. Based on professional judgment, the optimum temperature values (Topt) for green algae and cyanobacteria were therefore increased by 5°C to 31°C and 32°C respectively. The resulting fit to observed data (Figure 22) was good. Furthermore, the fish and zoobenthos fits were acceptable (Figure 33, see also Figure 70). Note that the bluegill are predicted to exhibit ammonia toxicity in 2001, an observation made possible by viewing biotic process rates. (Within rates graphs, animal mortality rates may be broken down into their various constituents, see (112)). Figure 33. Predicted and observed fish in Cahaba River, Alabama; predicted shiner mean biomass = 0.6 g/m2 compared to observed 0.5 g/m2. Cahaba River AL (CONTROL) Run on 01-26-09 12:07 PM -Shiner (g/m2 dry) - Bluegill (g/m2 dry) - Stoneroller (g/m2 dry) - Smallmouth Bas (g/m2 dry] Smallmouth Ba2 (g/m2 dry) Obs stonerollers (g/m2dry) Obs shiners (g/m2 dry) Obs bluegill (g/m2 dry) Obs bass (g/m2 dry) 8/26/2000 2/24/2001 8/25/2001 2/23/2002 8/24/2002 In another validation, published PCB data from New Bedford Harbor, Massachusetts, were used to verify the generality of the estuarine ecosystem bioaccumulation model. The observed concentrations of total PCBs in the water and bottom sediments in the Massachusetts site were set as constant values in a simulation of Galveston Bay, Texas. The predicted PCB concentrations in the various biotic compartments at the end of the simulation were then compared to the observed means and standard deviations in New Bedford Harbor (Figure 34). Considering that the sites and some of the species were different, the concordance in values provides a validation of the model for assessing bioaccumulation of chemicals in a "canonical" or representative estuarine environment. 39 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 2 Figure 34. Predicted and observed concentrations of PCBs in selected animals based on ecosystem calibration for Galveston Bay, Texas and exposure data (Connolly 1991) for New Bedford Harbor, Massachusetts. 10 0.1 CO u Q_ >- T3 O -Q | 0.01 eS^ £ 4^ x? x>* .><> ^ *v »o vCr vO x^ v* <>y <>y ^ >^ A third example of a validation is shown in Figure 21, which provides a visual comparison of predicted biomass and observed numbers per sample of chironomid larvae with dosing by an insecticide. No calibration was performed for either the fate or toxicity of the chemical. 40 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 3 3. PHYSICAL CHARACTERISTICS 3.1 Morphometry Volume Volume is a state variable and can be computed in several ways depending on availability of data and the site dynamics. It is important for computing the dilution or concentration of pollutants, nutrients, and organisms; it may be constant, but usually it is time varying. In the model, ponds, lakes, and reservoirs are treated differently than streams, especially with respect to computing volumes. The change in volume of ponds, lakes, and reservoirs is computed as: Morphometry: Simplifying Assumptions Base flow equation assumes a rectangular channel Site shapes are represented by idealized geometrical approximations Mean Depth may be held constant or user varying depth may be imported dVolume dt = Inflow - Discharge - Evap (2) where: dVolume/dt Inflow Discharge Evap derivative for volume of water (m3/d), inflow of water into waterbody (m3/d), discharge of water from waterbody (m3/d), and evaporation (m3/d), see (3). AQUATOX cannot successfully run if the volume of water in a site falls to zero. To avoid this condition, if the site's water volume falls below a minimum value (which is defined as a fraction of the initial condition using the parameter "Minimum Volume Frac." from the "Site" screen), all differentiation of state variables is suspended (except for the water volume derivative) until the water volume again moves above the minimum value. Differentiation of all state variables then resumes. A time series of evaporation may be entered in the "Site" screen in units of cubic meters per day. Otherwise, evaporation is converted from an annual value for the site to a daily value using the simple relationship: Evap = MeanEvap 365 0.02 54-Area (3) where: Evap MeanEvap 365 0.0254 Area mean daily evaporation (m /d) mean annual evaporation (in/yr), days per year (d/yr), conversion from inches to meters (m/in), and area of the waterbody (m2). 41 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 3 The user is given several options for computing volume including keeping the volume constant; making the volume a dynamic function of inflow, discharge, and evaporation; using a time series of known values; and, for flowing waters, computing volume as a function of the Manning's equation. Depending on the method, inflow and discharge are varied, as indicated in Table 3. As shown in equation (2), an evaporation term is present in each of these volume calculation options. In order to keep the volume constant, given a known inflow loading, evaporation must be subtracted from discharge. This will reduce the quantity of state variables that wash out of the system. In the dynamic formulation, evaporation is part of the differential equation, but neither inflow nor discharge is a function of evaporation as they are both entered by the user. When setting the volume of a water body to a known value, evaporation must again be subtracted from discharge for the volume solution to be correct. Finally, when using the Manning's volume equation, given a known discharge loading, the effects of evaporation must be added to the inflow loading so that the proper Manning's volume is achieved. (This could increase the amount of inflow loadings of toxicants and sediments to the system, although not significantly.) Table 3. Computation of Volume, Inflow, and Discharge Method Constant Dynamic Known values Manning Inflow InflowLoad InflowLoad InflowLoad ManningVol - State/dt + Discharge + Evap Discharge InflowLoad - Evap DischargeLoad InflowLoad - Evap + (State - Known Vals)/dt DischargeLoad The variables are defined as: InflowLoad DischargeLoad State KnownVals dt ManningVol user-supplied inflow loading (m /d); user-supplied discharge loading (m3/d); computed state variable value for volume (m3); time series of known values of volume (m3); incremental time in simulation (d); and volume of stream reach (m3), see (4). Figure 35 illustrates time-varying volumes and inflow loadings specified by the user and discharge computed by the model for a run-of-the-river reservoir. Note that significant drops in volume occur with operational releases, usually in the spring, for flood control purposes. 42 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 3 Figure 35. Volume, inflow, and discharge for a 4-year period in Coralville Reservoir, Iowa. 6.0E+07 O.OE+00 2.5E+08 2.0E+08 1.5E+08 1 .OE+08 5.0E+07 E 13 o Oct-74 Oct-75 Nov-76 Dec-77 Apr-75 May-76 Jun-77 Jul-78 O.OE+00 Inflow Discharge- Volume The time-varying volume of water in a stream channel is computed as: where: 7 CLength Width ManningVol = Y CLength Width dynamic mean depth (m), see (5); length of reach (m); and width of channel (m). (4) In streams the depth of water and flow rate are key variables in computing the transport, scour, and deposition of sediments. Time-varying water depth is a function of the flow rate, channel roughness, slope, and channel width using Manning's equation (Gregory, 1973), which is rearranged to yield: Y = Q Manning -3/5 (5) where: Q Manning Slope Width flow rate (m /s); Manning's roughness coefficient (s/m1/3); slope of channel (m/m); and channel width (m). The Manning's roughness coefficient is an important parameter representing frictional loss, but it is not subject to direct measurement. The user can choose among the following stream types: 43 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 3 . concrete channel (with a default Manning's coefficient of 0.020); dredged channel, such as ditches and channelized streams (default coefficient of 0.030); and natural channel (default coefficient of 0.040). These generalities are based on Chow's (1959) tabulated values as given by Hoggan (1989). The user may also enter a value for the coefficient. In the absence of inflow data, the flow rate is computed from the initial mean water depth, assuming a rectangular channel and using a rearrangement of Manning's equation: IDepth5/3 JSlope Width QBase = (6) Manning where: Idepth = mean depth as given in site record (m). QBase = base flow (m3/s); and The dynamic flow rate is calculated from the inflow loading by converting from m3/d to m3/s: 86400 where: Q = flow rate (m3/s); and Inflow = water discharged into channel from upstream (m3/d). Bathymetric Approximations The depth distribution of a water body is important because it determines the areas and volumes subject to mixing and light penetration. The shapes of ponds, lakes, reservoirs, and streams are represented in the model by idealized geometrical approximations, following the topological treatment of Junge (1966; see also Straskraba and Gnauck, 1985). The shape parameter P (Junge, 1966) characterizes the site, with a shape that is indicated by the ratio of mean to maximum depth. : ZMax Where: ZMean = mean depth (m); ZMax = maximum depth (m); and P = characterizing parameter for shape (unitless); P is constrained between -1.0 and 1.0 Shallow constructed ponds and ditches may be approximated by an ellipsoid where Z/ZMax = 0.6 and P = 0.6. Reservoirs and rivers generally are extreme elliptic sinusoids with values ofP 44 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 3 constrained to -1.0. Lakes may be either elliptic sinusoids, with P between 0.0 and -1.0, or elliptic hyperboloids with P between 0.0 and 1.0. Not all water bodies fit the elliptic shapes, but the model generally is not sensitive to the deviations. Based on these relationships, fractions of volumes and areas can be determined for any given depth (Junge, 1966). The AreaFrac function returns the fraction of surface area that is at depth Z given Zmax and P, which defines the morphometry of the water body. For example, if the water body were an inverted cone, when horizontal slices were made through the cone looking down from the top one could see both the surface area and the water/sediment boundary where the slice was made. This would look like a circle within a circle, or a donut (Figure 36). AreaFrac calculates the fraction that is the donut (not the donut hole). To get the donut hole, 1 - AreaFrac is used. AreaFrac = (1-P)- ^ + P (-^ / (9) ZMax ZMax 6.0 3.0-(1.0-P)-(^f-2.0-P-(^f VolFrac = ZMax ZMax ZMax / j Q\ 3.0 + P where: AreaFrac = fraction of area of site above given depth (unitless); VolFrac = fraction of volume of site above given depth (unitless); and Z = depth of interest (m). For example, the fraction of the volume that is epilimnion can be computed by setting depth Zto the mixing depth. Furthermore, by setting Z to the depth of the euphotic zone, where primary production exceeds respiration, the fraction of the area available for colonization by macrophytes and periphyton can be computed: T , ZEuphotic (ZEuphotic\ _. FracLit = (l-P) + P-\ (11) ZMax \ ZMax ) A relatively deep, flat-bottomed basin would have a small littoral area and a large sublittoral area (Figure 36). Figure 36. 45 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 3 If the site is an artificial enclosure then the available area is increased accordingly: Area + EnclWallArea FracLittoral = FracLit Area otherwise (12) where: FracLittoral ZEuphotic Area EnclWallArea FracLittoral = FracLit = fraction of site area that is within the euphotic zone (unitless); = depth of the euphotic zone, is assumed to be 1% of surface light and calculated as 4.605/Extinct (m) see (40); = site area (m2); and = area of experimental enclosure's walls (m2). Figure 37. Area as a function of depth RESERVOIR (P = -0.6) 1 3 5 7 9 11 13 15 17 19 21 23 25 2 4 6 8 10 12 14 16 18 20 22 24 DEPTH(m) Figure 38. Volume as a function of depth RESERVOIR (P = -0.6) 1 3 5 7 9 11 13 15 17 19 21 23 25 2 4 6 8 10 12 14 16 18 20 22 24 DEPTH (m) If a user wishes to model a simpler system, the bathymetric approximations may be bypassed in favor of a more rudimentary set of assumptions via an option in the "site data" screen. When the user chooses not to "use bathymetry" the system is assumed to have vertical walls; the system is assumed to have a constant area as a function of depth; the system's depth may be calculated at any time as water volume divided by surface area. This option may be useful when linking data from other models to AQUATOX as the horizontal spatial domain of AQUATOX remains unchanged over time. However, a system will not undergo dynamic stratification based on water temperature unless the more complex bathymetric approximations are utilized ((8) to (11)). 46 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 3 Dynamic Mean Depth AQUATOX normally uses an assumption of unchanging mean depth (i.e., mean over the site area). However, under some circumstances, and especially in the case of streams and reservoirs, the depth of the system can change considerably over time, which could result in a significantly different light climate for algae. For this reason, an option to import mean depth in meters has been added. A daily time-series of mean depth values may be imported into the software (using an interface found within the "site" screen by pressing the "Show Mean Depth Panel" button.) A time-series of mean depth values can be estimated given known water volumes or can be imported from a linked water hydrology model. The user-input dynamic mean depth affects the following portions of AQUATOX: Light climate, see (43); Calculation of biotic volumes for sloughing calculations, see (74); Calculation of vertical dispersion for stratification calculations, Thick in equation (18); Calculation of sedimentation for plants & detritus, Thick in (165); Oxygen reaeration, see (190); . Toxicant photolysis and volatilization, Thick in (320) and (331). Habitat Disaggregation Riverine environments are seldom homogeneous. Organisms often exhibit definite preferences for habitats. Therefore, when modeling streams or rivers, animal and plant habitats are broken down into three categories: "riffle," "run," and "pool." The combination of these three habitat categories make up 100% of the available habitat within a riverine simulation. The preferred percentage of each organism that resides within these three habitat types can be set within the animal or plant data. Within the site data, the percentage of the river that is composed of each of these three habitat categories also can be set. It should be noted that the habitat percentages are considered constant over time, and thus would not capture significant changes in channel morphology and habitat distribution due to major flooding events. These habitats affect the simulations in two ways: as limitations on photosynthesis and consumption and as weighting factors for water velocity (see 3.2 Velocity). Each animal and plant is exposed to a weighted average water velocity depending on its location within the three habitats. This weighted velocity affects all velocity-mediated processes including entrainment of invertebrates and fish, breakage of macrophytes and scour of periphyton. The reaeration of the system also is affected by the habitat-weighted velocities. Limitations on photosynthesis and consumption are calculated depending on a species' preferences for habitats and the available habitats within the water body. If the species preference for a particular habitat is equal to zero then the portion of the water body that contains that particular habitat limits the amount of consumption or photosynthesis accordingly. 47 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 3 where: HabitatLimit species = Preference habitat = Per cent habitat = s > Per cent habitat \ 100 (13) fraction of site available to organism (unitless), used to limit ingestion, see (91), and photosynthesis, see (35), (85); preference of animal or plant for the habitat in question (percentage); and percentage of site composed of the habitat in question (percentage). It is important to note that the initial condition for an animal that is entered in g/m is an indication of the total mass of the animal over the total surface area of the river. Because of this, density data for various benthic organisms, which is generally collected in a specific habitat type, cannot be used as input to AQUATOX until these values have been converted to represent the entire surface area. This is especially true in modeling habitats; for example, an animal could have a high density within riffles, but riffles might only constitute a small portion of the entire system. 3.2 Velocity If the user has site-specific velocity data, this may be entered on the "site data" screen in units of cm/s. Otherwise, velocity is calculated as a simple function of flow and cross-sectional area: T, , AvgFlow Velocity = s XSecArea 86400 -100 (14) where Velocity AvgFlow XSecArea 86400 100 velocity (cm/s), average flow over the reach (m3/d), cross sectional area (m2), s/d, and cm/m. where: Inflow Discharge AvgFlow = Inflow + Discharge 2 (15) flow into the reach (m /d); flow out of the reach (m3/d). It is assumed that this is the velocity for the run of the stream (user entered velocities are also assumed to pertain to the run of the screen). No distinction is made in terms of vertical differences in velocity in the stream. Following the approach and values used in the DSAMMt model (Caupp et al. 1995), the riffle velocity is obtained by using a conversion factor that is 48 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 3 dependent on the discharge. Unlike the DSAMMt model, pools also are modeled, so a conversion factor is used to obtain the pool velocity as well (Table 4). Table 4. Factors relating velocities to those of the average reach. Flows (Q = discharge) Q < 2.59e5 m3/d 2.59e5 m3/d < Q < 5.18e5 m3/d 5. 18e5 m3/d < Q < 7.77e5 m3/d Q > 7.77e5 m3/d Run Velocity 1.0 1.0 1.0 1.0 Riffle Velocity 1.6 1.3 1.1 1.0 Pool Velocity 0.36 0.46 0.56 0.66 Figure 39. Predicted velocities in an Ohio stream according to habitat. 400 350 CD 00 . Riffle Pool 3.3 Washout Transport out of the system, or washout, is an important loss term for nutrients, floating organisms, and dissolved toxicants in reservoirs and streams. Although it is considered separately for several state variables, the process is a general function of discharge: where: Washout State Washout= _ Discharge Volume State loss due to being carried downstream (g/m -d), and concentration of dissolved or floating state variable (g/m3). (16) 49 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 3 3.4 Stratification and Mixing Thermal stratification is handled in the simplest form consistent with the goals of forecasting the effects of nutrients and toxicants. Lakes and reservoirs are considered in the model to have two vertical zones: epilimnion and hypolimnion (Figure 40); the metalimnion zone that separates these is ignored. Instead, the thermocline, or plane of maximum temperature change, is taken as the separator; this is also known as the mixing depth (Hanna, 1990). Dividing the lake into two vertical zones follows the treatment of Imboden (1973), Park et al. (1974), and Straskraba and Gnauck (1983). The onset of stratification is considered to occur when the mean water temperature exceeds 4 deg. and the difference in temperature between the epilimnion and hypolimnion exceeds 3 deg.. Overturn occurs when the temperature of the epilimnion is less than 3 deg., usually in the fall. Winter stratification is not modeled, unless manually input. For simplicity, the thermocline is generally assumed to occur at a constant depth. Alternatively, a user-specified time-varying thermocline depth may be specified, see the section on modeling reservoirs below. Figure 40. Thermal stratification in a lake; terms defined in text Stratification: Simplifying Assumptions Two vertical zones modeled; metalimnion is ignored Flowing waters are assumed not to stratify Stratification occurs when vertical temperature difference exceeds three degrees Winter stratification is not modeled Thermocline occurs at constant depth except when user enters time series Wind action is implicit in vertical dispersion calculations Epilimnion Thermocline VertDispersion There are numerous empirical models relating thermocline depth to lake characteristics. AQUATOX uses an equation by Hanna (1990), based on the maximum effective length (or fetch). The dataset includes 167 mostly temperate lakes with maximum effective lengths of 172 to 108,000 m and ranging in altitude from 10 to 1897 m. The equation has a coefficient of determination r2 = 0.850, meaning that 85 percent of the sum of squares is explained by the regression. Its curvilinear nature is shown in Figure 41, and it is computed as (Hanna, 1990): \og(MaxZMix) = 0.336- \og(Length) - 0.245 (17) where: 50 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 3 MaxZMix = maximum mixing depth under stratified conditions (thermocline depth) for lake (m); and Length = maximum effective length for wave setup (m, converted from user- supplied km). Figure 41. Mixing depth as a function of fetch MAXIMUM MIXING DEPTH 100 11500 22900 34300 5800 17200 28600 40000 LENGTH (m) Wind action is implicit in this formulation. Wind has been modeled explicitly by Baca and Arnett (1976, quoted by Bowie et al., 1985), but their approach requires calibration to individual sites, and it is not used here. Vertical dispersion for bulk mixing is modeled as a function of the time-varying hypolimnetic and epilimnetic temperatures, following the treatment of Thomann and Mueller (1987, p. 203; see also Chapra and Reckhow, 1983, p. 152; Figure 42): VertDispersion = Thick I HypVolume rr,t-l rr,t+l -L hypo ~ J- h lypo \ ThermoclArea Deltat T'epi - (18) rypo where: VertDispersion = Thick = vertical dispersion coefficient (m2/d); distance between the centroid of the epilimnion and the centroid of the hypolimnion, effectively the mean depth (m); 51 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 3 HypVolume ThermoclArea Deltat ji t-i ji t+i 1 hypo ' * hypo T~T t T~T t 1 epi > * hypo volume of the hypolimnion (m ); area of the thermocline (m2); time step (d); temperature of hypolimnion one time step before and one time step after present time (deg. C); and temperature of epilimnion and hypolimnion at present time (deg.C). Stratification can break down temporarily as a result of high throughflow. This is represented in the model by making the vertical dispersion coefficient between the layers a function of discharge for sites with retention times of less than or equal to 180 days (Figure 43), rather than temperature differences as in equation 11, based on observations by Straskraba (1973) for a Czech reservoir: VertDispersion = 1.3 7 J Q4 Retention -2.269 (19) and: Retention = Volume TotDischarge (20) where: Retention = Volume = TotDischarge = retention time (d); volume of site (m3); and total discharge (m3/d). Figure 42. Vertical dispersion as a function of temperature differences 25 iu O o O (A 0.01 0 -I 1 1 1 1 1 1 1 1 1 1 1 h1- 0.001 12/30 02/28 04/29 06/28 08/27 10/26 12/25 DAY Epilimnion Temp. -*- Hypolimnion Temp. Vert. Dispersion (sq m/d) 4 degrees 52 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 3 Figure 43. Vertical dispersion as a function of retention time 100 Q w tr LJJ D_ 10 . OH LU 0.1 VERTICAL DISPERSION 180 162 144 126 108 90 72 54 36 171 153 135 117 99 81 63 45 27 RETENTION TIME (d) The bulk vertical mixing coefficient is computed using site characteristics and the time-varying vertical dispersion (Thomann and Mueller, 1987): BulkMixCoeff = VertDispersion ThermoclArea Thick (21) where: BulkMixCoeff = bulk vertical mixing coefficient (m3/d), ThermoclArea = area of thermocline (m2). Turbulent diffusion of biota and other material between epilimnion and hypolimnion is computed separately for each segment for each time step while there is stratification: T hPi'ff _ BulkMixCoeff L UtOLJlJJ e^ ^ ' ( (^OTICcompartment, hypo ~ ^OWC compartment, epi/ Volumeepi _ BulkMixCoeff _/ Hi ULJIJJ -faypO \ \^ OTIC compartment, epi ~ ^ OTIC compartment, hypo/ (22) (23) where: TurbDiff Volume Cone turbulent diffusion for a given zone (g/m3-d); volume of given segment (m3); and concentration of given compartment in given zone (g/m3). 53 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 3 The effects of stratification, mixing due to high throughflow, and overturn are well illustrated by the pattern of dissolved oxygen levels in the hypolimnion of Lake Nockamixon, a eutrophic reservoir in Pennsylvania (Figure 44). Figure 44. Stratification and mixing in Lake Nockamixon, Pennsylvania as shown by hypolimnetic dissolved oxygen 14 10 S, 8 o> "o 4 00 CO onset of stratification , overturn jf '' high throughflow 01/01182 03/07/82 05/11182 07/15/82 09/18/82 11 /22/S2 Modeling Reservoirs and Stratification Options Stratification assumptions and equations based on lake characteristics may not be appropriate for modeling reservoirs. Moreover, a lake may have a unique morphometry or chemical composition that renders inappropriate the equations presented above. For this reason, a "stratification options" screen is available (through the "site" screen or "water-volume" screen) that allows a user to specify the following characteristics of a stratified system: a constant or time-varying thermocline depth; options as to how to route inflow and outflow water; and the timing of stratification. Water volumes for each segment are calculated as a function of the overall system volume and the thermocline depth (see (10)). Because of this, if a time-varying thermocline depth is specified, water from one segment must usually be transferred into the other segment, along with the state variables within that water. In this manner, specifying a time-varying thermocline depth has the potential to promote mixing between layers. Alternatively, using the linked-mode model, two stratified segments may be specified with water volumes that are calculated independently from the thermocline depth; see section 3.8 for more details about stratification in linked-mode. By default, AQUATOX routes inflow and outflow to and from both segments as weighted by volume. For example, if the hypolimnion has twice as much volume as the epilimnion, twice as 54 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 3 much inflow water will be routed to the hypolimnion as to the epilimnion (and twice as much outflow water will be routed from the hypolimnion). The user has the option to route all inflow and outflow waters to and from either segment. In this case, all of the nutrients, chemicals, and other loadings within the inflow water will be routed directly to the specified segment and will not be transferred to the other segment except through turbulent diffusion or overturn. Atmospheric and point-source loadings are assumed to be routed to the epilimnion in all cases (unless a linked-mode model is used in which case more flexibility is present). Additionally, if a user has information about the timing of stratification, this may be specified on the "stratification-options" entry screen. This can be used to specify winter stratification, for example, or precise periods of stratification for each year modeled. If only one year of stratification dates are entered and multiple years are modeled, all years are assumed to stratify and overturn on the dates specified in the user input (regardless of the year specified). 3.5 Temperature Temperature is an important controlling factor in the model. Virtually all processes are temperature-dependent. They include stratification; biotic processes such as decomposition, photosynthesis, consumption, respiration, reproduction, and mortality; and chemical fate processes such as microbial degradation, volatilization, hydrolysis, and bioaccumulation. On the other hand, temperature rarely fluctuates rapidly in aquatic systems. Default water temperature loadings for the epilimnion and hypolimnion are represented through a simple sine approximation for seasonal variations (Ward, 1963) based on user-supplied observed means and ranges (Figure 45): r + r ** / i n TempRcmge Temperature = TempMean + (-1.0 2 (24) (sin(0.0174533 (0.987 (Day + PhaseShift) - 30))))] where: Temperature = average daily water temperature (deg. C); TempMean = mean annual temperature (deg. C); TempRange = annual temperature range (deg. C), Day = day of year (d); and PhaseShift = time lag in heating (= 90 d). Observed temperature loadings should be entered if responses to short-term variations are of interest. This is especially important if the timing of the onset of stratification is critical, because stratification is a function of the difference in hypolimnetic and epilimnetic temperatures (see Figure 42). It also is important in streams subject to releases from reservoirs and other point- source temperature impacts. 55 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 3 3.6 Light Light is important as the controlling factor for photosynthesis and photolysis. The default incident light function formulated for AQUATOX is a variation on the temperature equation, but without the lag term: Light: Simplifying Assumptions Ice cover is assumed when the average water temperature drops below 3 degrees centigrade. Photoperiod is approximated by Julian date (day of year) Average daily light is the program default, although hourly light may be simulated Solar = LightMean + Ll8htRan8e . 01 74533 . Day . L 76) . FracLi Light Fracught =1.0- 0.98(Canopy) where: Solar = LightMean = LightRange = Day Frac Light Canopy = average daily incident light intensity (ly/d); mean annual light intensity (ly/d); annual range in light intensity (ly/d); day of year (d, adjusted for hemisphere); fraction of site that is shaded; and user input fraction of site that is tree shaded. (25) The derived values are given as average light intensity in Langleys per day (Ly/d = 10 r\ kcal/m -d). An observed time-series of light also can be supplied by the user; this is especially important if the effects of daily weather conditions are of interest. If the average water temperature drops below 3 deg.C, the model assumes the presence of ice cover and decreases transmitted light to 15% of incident radiation. (This has changed from 33% in Release 2.2.) This reduction, due to the reflectivity and transmissivity of ice and snow, is an average of widely varying values summarized by Wetzel (2001). For estuaries, average water temperature must fall below -1.8 deg.C before the model assumes ice cover due to the influence of salinity. Shade can be an important limitation to light, especially in riparian systems. A user input "fraction of site subject to shade from a canopy" parameter can be entered either as a constant or as a time-series within the "Site" input screen. This parameter can be left as zero for no shading effects on light. Transmission of light through a riparian (stream-side) canopy is a combination of diffuse and direct transmission (Canham et al. 1990). The average of four forest types from closed hemlock to open spruce (and cypress) forests is 2% of incident radiation (Canham et al. 1990). Detailed studies in a Midwestern mixed deciduous forest confirm this value for the summer months, although transmission increased to 40% in winter (Oliphant et al. 2006). A value of 2% transmission for a closed canopy is used in AQUATOX. If the density of canopy 56 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 3 varies during the year, then a time-series should be provided, keeping in mind that the 2% transmission will still apply to the fraction of canopy that is indicated. Photoperiod is an integral part of the photosynthesis formulation. It is approximated using the Julian date following the approach of Stewart (1975) (Figure 46): 12 + A- cos (380 Photoperiod = 365 248) 24 (26) where: Photoperiod = A Day fraction of the day with daylight (unitless); converted from hours by dividing by 24; hours of daylight minus 12 (d); and day of year (d, converted to radians). A is the difference between the number of hours of daylight at the summer solstice at a given latitude and the vernal equinox, and is given by a linear regression developed by Groden (1977): where: Latitude Sign A = 0.1414 Latitude - Sign 2.413 latitude (deg., decimal), negative in southern hemisphere; and 1.0 in northern hemisphere, -1.0 in southern hemisphere. (27) Figure 45. Annual Temperature Figure 46. Photoperiod as a Function of Date TEMPERATURE IN A MIDWESTERN POND 77 153 229 305 39 115 191 267 343 JULIAN DAY : 0.65 0.55 Q "5 0.45 tz I 0.4 ro LL 0.35 1 53 105 157 209 261 313 365 27 79 131 183 235 287 339 Julian Date Latitude 40 N Latitude 40 S 57 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 3 Hourly Light When the model is run with an hourly time-step, solar radiation is calculated as variable during the course of each day. The following equation is used to distribute the average daily incident light intensity over the portion of the day with daylight hours. Solar, Solar,,, Photoperiod sin brae Day Passed 1-Photoperiod 2 Photoperiod (28) Frac Light where: Solar houriy Photoperiod FracDayPassed Frac Light solar radiation at the given time-step (ly/d); average daily incident light intensity (ly/d), see (25); fraction of the day with daylight (unitless); see (26); fraction of the day that has passed (unitless) fraction of site that is un-shaded, (frac., 1.0-user input shade); A user may enter a constant or time-series shade variable in the site window ("Fraction of Site that is Shaded"). When this input is utilized then the Frac Light variable is calculated. Figure 47: Average light per day is distributed during daylight hours in a semi-sinusoidal pattern based on photoperiod. 1200 1000 - 800 - 600 - 400 - 200 - 0 0 10 20 30 Hours 50 58 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 3 3.7 Wind Wind is an important driving variable because it determines the stability of blue-green algal blooms, affects reaeration or oxygen exchange, and controls volatilization of some organic chemicals. Wind also can affect the depth of stratification for estuaries. Wind is usually measured at meteorological stations at a height of 10 m and is expressed as m/s. If site data are not available, default variable wind speeds are represented through a Fourier series of sine and cosine terms; the mean and twelve additional harmonics seem to effectively capture the variation (Figure 48): Wind: Simplifying Assumptions If site data are not available a Fourier series is used to represent wind loadings Win d = 365 sinCoeff . Sm where: Wind = wind speed; amplitude of the Fourier series (m/s); CosCoeffo = cosine coefficient for the 0-order harmonic, which is the mean wind speed (default = 3 m/s); CosCoeffn = cosine coefficient for the n* -order harmonic; Day = day of year (d); SinCoeffn = sine coefficient for the nth-order harmonic; Freqn = selected frequency for the nth- order harmonic. This default loading is based on an annual cycle of data taken from the Buffalo, NY airport. Therefore, it has a 365-day repeat, representative of seasonal variations in wind. Frequencies were selected to ensure that the standard deviation of the Fourier series and the data were closely matched. The frequency of wind-speeds of less than three meters per second were also precisely matched to observed data as well as the periodicity of wind-events. The Fourier approach is quite useful because the mean can be specified by the user and the variability will be imposed by the function. If ice cover is predicted, wind is set to 0. A user also may input a site-specific time series, which may be important where the timing of a cyanobacteria bloom or reaeration is of interest. 59 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 3 Figure 48. Default wind loadings for Onondaga Lake with mean = 4.17 m/s. ONONDAGA LAKE, NY (CONTROL) Run on 04-24-08 9:07 PM (Epilimnion Segment) Wind (m/s) 1/12/1989 5/12/1989 9/9/1989 1/7/1990 5/7/1990 9/4/1990 1/2/1991 3.8 Multi-Segment Model AQUATOX Release 3 includes the capability to link AQUATOX segments together, tracking the flow of water and the passage of state variables from segment to segment. Some general guidelines for using this model follow: Multi-Segment Model: Simplifying Assumptions All linked segments have an identical set of state variables Each segment is well mixed Linkages between segments may be unidirectional or bidirectional Dynamic stratification does not apply; stratified pairs of segments must be specified by the user . All linked segments must have an identical set of state variables. (State variables that do not occur in one segment may be set to zero there.) Parameters pertaining to animal, plant, and chemical state variables (i.e. "underlying data") are considered global to the entire linked system. If the user changes one of these parameters in one segment, this parameter changes within all segments. On the other hand, "site" parameters, initial conditions, and boundary conditions are unique to each segment. State variables can pass from segment to segment through active upstream and downstream migration, passive drift, diffusion, and bedload. Mass balance of all state variables is maintained throughout a multi-segment simulation. There are two types of linkages that may be specified between individual segments, "cascade links" and "feedback links." A cascade link is unidirectional; there is no potential for water or 60 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 3 state variable flow back upstream. Segments that are linked together by cascade linkages are solved separately from one another moving from upstream to downstream. This is particularly useful when modeling faster flowing rivers and streams. A feedback link allows for water or state variables to flow in both directions. For bookkeeping purposes, water flows are required to be unidirectional (i.e. entered water flows over a feedback link must not be negative). However, two feedback links may be specified simultaneously (in opposite directions) to allow for bidirectional water flows. Feedback links may also be subject to diffusion; a diffusion coefficient, characteristic length, and cross section must be entered for diffusion to be calculated, see (32). Segments that are linked together by feedback links are solved simultaneously. There may only be one contiguous set of segments linked together by feedback linkages within a simulation (i.e. the model will not solve a "feedback" set of segments followed by downstream cascade segments followed by more feedback segments below that.) Figure 49 gives an example of a simulation in which cascade segments and feedback segments are both included. In this case, AQUATOX solves the simulation from the top down, solving each segment 1-4, 6, and 6b individually before moving on to solve the feedback segments simultaneously. Finally, segments 11-14 are be solved individually using the results from the simultaneous segment run. Figure 49: An example of feedback and cascade segments linked together. \ x ) Feedback Seg. Cascade Seg. Feedback Link Cascade Link 61 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 3 Stratification and the Multi-Segment Model Dynamic stratification as described in section 3.4 does not apply to the multi-segment model. Instead, a user may specify two linked segments as a stratified pair. In this case, the segments must be linked together with a feedback linkage. A "stratification" screen within each segment's main interface allows a user to specify whether a segment is part of a stratified pair and, if so, whether it is the epilimnion or the hypolimnion segment. When two segments are set up as stratified together, the thermocline area is defined by the user- entered cross section between. Annual cycles of stratification and overturn may be specified using the time varying water flows and dispersion coefficients. As was the case in the dynamic stratification model, fish automatically migrate to the epilimnion in the case of hypoxia in the lower segment. Sinking phytoplankton and suspended detritus in the epilimnion segment fall into the designated hypolimnion segment. The light climate of the bottom segment is limited to that light which penetrates the segment defined as the epilimnion. When the linked system has enough specified throughflow between the epilimnion and hypolimnion segment, it is considered to be "well mixed." This is defined as when the average daily water flow between segments is greater than 30% of the total water volume in both segments. In this case, fish are assumed to have an equal preference to both segments and they migrate to equality in a biomass basis. (This allows fish to return to the hypolimnion if it had earlier been vacated due to anoxia.) Another implication of a well-mixed stratified system (in linked mode) is that a weighted average of light climate is used when calculating plant productivity. The calculation of LightLimit for plants (38) is based on a thickness-weighted average of algal biomass and sediment throughout the entire thickness of the system. This prevents unreasonable model results due to the light climate in a very thin epilimnion, for example. Because the system is well-mixed, suspended algae should instead be subject to the light climate throughout the water column. State Variable Movement in the Multi-Segment Model To maintain mass balance, all state variables that are subject to washout or passive drift are also added to any downstream linked segments. The calculation for this process is as follows: Washin= y Washoutupanam VolumeUpstream FracWash,,^ / J -TT J \ S upstream link ' ^" Dowstream Segment In the case of toxicants that are absorbed to or contained within a drifting state variable, the following equation is used: VolumeUpstream FracWas^^ ' ' T/ / upstream links ' OlUmt Dow ,stream Segment 62 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 3 where: Washin = inflow load from upstream segment (unit/Ldownstream'd); Washout upstream = washout from upstream segment (unit/LUpStream'd), see (16); Volume segment = volume of given segment (m3); FracWashThiSLink = fraction of upstream segment's outflow that goes to this particular downstream segment (unitless); WashinToxCamer = inflow load of toxicant sorbed to a carrier from an upstream Segment (Hg/Ldownstream'd); Washout earner = washout of toxicant carrier from upstream (mg/Lupsireanrd); PPBcarrier = concentration of toxicant in carrier upstream (|j,g/kg), see (310); le-6 = units conversion (kg/mg) This Washin term is added to all derivatives for state variables that are suspended in the water column and subject to drift or "washout." Dissolved state variables are subject to diffusion across feedback links. r»vcr DiffCoeff Area i \ Diffusion = JJ JJ (Conc0 -Cone ) (32) CharLength where: DiffusionThisSeg = gain of state variable due to diffusive transport over the feedback link between two segments, (unit/d); DiffCoeff = dispersion coefficient of feedback link, (m2 /d); Area = surface area of the feedback link (m ); CharLength = characteristic mixing length of the feedback link, (m); = concentration of state variable in the relevant segment, (unit/m3); 63 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 4 4. BIOTA The biota consists of two main groups, plants and animals; each is represented by a set of process-level equations. In turn, plants are differentiated into algae and macrophytes, represented by slight variations in the differential equations. Algae may be either phytoplankton or periphyton. Phytoplankton are subject to sinking and washout, while periphyton are subject to substrate limitation and scour by currents. Bryophytes and freely-floating macrophytes are modeled as special classes of macrophytes, limited by nutrients in the water column. These differences are treated at the process level in the equations (Table 5). All are subject to habitat availability, but to differing degrees. Table 5. Significant Differentiating Processes for Plants Biota: Simplifying Assumptions Biomass is simulated but not numbers of individual organisms Responses are simulated as averages for the entire group Plant Type Phytoplankton Periphyton Benthic Macrophytes Rooted-Floating Macrophytes Free-Floating Macrophytes Bryophytes Nutrient Lim. J J J a Current Lim. J Sinking J Washout J J Sloughing J Breakage J J J a Habitat J J J J J a Animals are subdivided into invertebrates and fish; the invertebrates may be pelagic invertebrates, benthic insects or other benthic invertebrates. These groups are represented by different parameter values and by variations in the equations. Insects are subject to emergence and therefore are lost from the system, but benthic invertebrates are not. Fish may be represented by both juveniles and adults, which are connected by promotion. One fish species can be designated as multi-year with up to 15 age classes connected by promotion. Differences are shown in Table 6. In addition, a bioaccumulative endpoint such as bald eagle, dolphin, or mink that feeds on aquatic compartments can be simulated; it is defined by feeding preferences, biomagnification factor, and clearance rate. 64 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 4 Table 6. Significant Differentiating Processes for Animals Animal Type Pelagic Invert. Benthic Invert. Benthic Insect Fish Washout J Drift J J Entrainment J J a Emergence J Promotion/ Recruitment a Multi-year a 4.1 Algae o 9 -expressed as g/m for phytoplankton, but as g/m for periphyton- The change in algal biomass- is a function of the loading (especially phytoplankton from upstream), photosynthesis, respiration, excretion or photorespiration, nonpredatory mortality, grazing or predatory mortality, sloughing, and washout. As noted above, phytoplankton also are subject to sinking. If the system is stratified, turbulent diffusion also affects the biomass of phytoplankton. Plants: Simplifying Assumptions Photosynthesis is modeled as a maximum observed rate multiplied by reduction factors. The reduction factors are assumed to be independent of one another. Intracellular storage of nutrients is not modeled; constant stoichiometry within species is assumed For each individual nutrient, saturation kinetics is assumed Algae exhibit a nonlinear, adaptive response to temperature changes Low temperatures are assumed not to affect algal mortality The ratio between biovolume and biomass is assumed to be constant for a given growth form Constant chlorophyll a to biomass ratios are assumed within algae groups Phytoplankton-speciflc Phytoplankton other than cyanobacteria are assumed to be mixed throughout the well-mixed layer unless specified as "surface floating." In the event of ice cover, all phytoplankton will occur in the top 2 m Sinking of phytoplankton is modeled as a function of physiological state Phytoplankton are subject to downstream drift as a simple function of discharge To model phytoplankton (and zooplankton) residence time, an implicit assumption may be made that upstream reaches included in the "Total River Length " have identical environmental conditions as the reach being modeled Cyanobacteria-specific By default cyanobacteria are specified as "surface floating" in which case they are assumed to be located in the top 0.1 m unless limited by lack of nutrients or sufficient wind occurs in which case they are located within the top 3 m. This default assumption (that cyanobacteria float) can be changed by the user. The averaging depth for "surface floating" plants is three meters to more closely correspond to monitoring data. Cyanobacteria are not severely limited by nitrogen due to facultative nitrogen fixation (if N less than 1A KN) Periphyton-specific Periphyton are limited by slow currents that do not replenish nutrients and carry away senescent biomass Periphyton are assumed to adapt to the ambient conditions of a particular channel Periphyton are defined as including associated detritus; non-living biomass is modeled implicitly Macrophyte-speciflc Macrophytes occupy the littoral zone Rooted macrophytes are not limited by nutrients but are assumed to take up necessary nutrients from bottom sediments 65 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 4 dBiomass Phyto dt = Loading + Photosynthesis - Respiration - Excretion -Mortality - Predation ± Sinking ± Floating - Washout + Washin ± TurbDiff + Diffusion Slough Seg (33) dBiomass ~dt Pen = Loading + Photosynthesis - Respiration - Excretion -Mortality - Predation + Sedperi - Slough (34) where: dBiomass/dt Loading Photosynthesis Respiration Excretion Mortality Predation Washout Washin Sinking Floating TurbDiff Diffusionseg Slough Sedpen = change in biomass of phytoplankton and periphyton with respect to T r\ time (g/m -d and g/m -d); = boundary-condition loading of algal group (g/m3-d and g/m2-d); T r\ = rate of photosynthesis (g/m -d and g/m -d), see (35); T r\ = respiratory loss (g/m -d and g/m -d), see (63); = excretion or photorespiration (g/m3-d and g/m2-d), see (64); T r\ = nonpredatory mortality (g/m -d and g/m -d), see (66); = herbivory (g/m3-d and g/m2-d), see (99); = loss due to being carried downstream (g/m3-d), see (129); = loadings from upstream segments (linked segment version only, g/m3-d), see (30); = loss or gain due to sinking between layers and sedimentation to bottom (g/m3-d), see (69); = loss from the hypolimnion or gain to the epilimnion due to the floatation of "surface-floating" phytoplankton. 100% of "surface- floating" phytoplankton that arrive in the hypolimnion through loadings or water flows are set to immediately float. turbulent diffusion (g/m3-d), see (22) and (23); = gain or loss due to diffusive transport over the feedback link between two segments, (g/m3-d), see (32); = Scour loss of Periphyton or addition to linked Phytoplankton, see (75); and = Sedimentation of Phytoplankton to Periphyton, see (83). Figure 50 and Figure 51 are examples of the predicted changes in biomass and the processes that contribute to these changes in a eutrophic lake. Note that photosynthesis and predation dominate the diatom rates, with respiration much less important during the growing season. 66 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 4 Figure 50. Predicted algal biomass in Lake Onondaga, New York ONONDAGA LAKE, NY (PERTURBED) Run on 04-23-08 2:59 PM (Epilimnion Segment) 1/12/1989 5/12/1989 9/9/1989 1/7/1990 5/7/1990 9/4/1990 1/2/1991 Cyclotella nan (mg/L dry) Greens (mg/L dry) - Phyt, Blue-Gre (mg/L dry) Cryptomonad (mg/L dry) Figure 51. Predicted process rates for diatoms in Lake Onondaga, New York ONONDAGA LAKE, NY (PERTURBED) Run on 04-23-08 2:59 PM (Epilimnion Segment) 3/11/1989 9/9/1989 3/10/1990 9/8/1990 Cyclotella nan Photosyn (Percent) Cyclotella nan Respir (Percent) Cyclotella nan Excret (Percent) - Cyclotella nan Other Mort (Percent) Cyclotella nan Predation (Percent) Cyclotella nan Washout (Percent) Cyclotella nan Sediment (Percent) Cyclotella nan TurbDiff (Percent) Cyclotella nan SinkToHypo (Percent) 67 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 4 Photosynthesis is modeled as a maximum observed rate multiplied by reduction factors for the effects of toxicants, habitat, and suboptimal light, temperature, current, and nutrients: Photosynthesis = PMax PProdLimit Biomass HabitatLimit SaltEffect (35) The limitation of primary production in phytoplankton is: PProdLimit = LtLimit NutrLimit TCorr FracPhoto (36) Periphyton have an additional limitation based on available substrate, which includes the littoral bottom and the available surfaces of macrophytes. The macrophyte surface area conversion is based on the observation of 24 m2 periphyton/m2 bottom (Wetzel, 1996) and assumes that the observation was made with 200 g/m macrophytes. PProdLimit = LtLimit NutrLimit VLimit TCorr FracPhoto (37) ( FracLittoral + SurfAreaConv BiomassMacroPhytes) where: Pmax = maximum photosynthetic rate (1/d); LtLimit = light limitation (unitless), see (38); NutrLimit = nutrient limitation (unitless), see (55); Vlimit = current limitation for periphyton (unitless), see (56); TCorr = limitation due to suboptimal temperature (unitless), see (59); HabitatLimit = in streams, habitat limitation based on plant habitat preferences (unitless), see (13). SaltEffect = effect of salinity on photosynthesis (unitless); FracPhoto = reduction factor for effect of toxicant on photosynthesis (unitless), see (421); FracLittoral = fraction of area that is within euphotic zone (unitless) see (11); SurfAreaConv = surface area conversion for periphyton growing on macrophytes (0.12 m2/g); BiomassMacm = total biomass of macrophytes in system (g/m ); and Biomassperi = biomass of periphytic algae (g/m ). Under optimal conditions, a reduction factor has a value of 1; otherwise, it has a fractional value. Use of a multiplicative construct implies that the factors are independent. Several authors (for example, Collins, 1980; Straskraba and Gnauck, 1983) have shown that there are interactions among the factors. However, we feel the data are insufficient to generalize to all algae; therefore, the simpler multiplicative construct is used, as in many other models (Chen and Orlob, 1975; Lehman et al., 1975; J0rgensen, 1976; Di Toro et al., 1977; Kremer and Nixon, 1978; Park et al., 1985; Ambrose et al., 1991). Default parameter values for the various processes are taken primarily from compilations (for example, J0rgensen, 1979; Collins and Wlosinski, 1983; Bowie et al., 1985); they may be modified as needed. 68 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION CHAPTER 4 Light Limitation Because it is required for photosynthesis, light is a very important limiting variable. It is especially important in controlling competition among plants with differing light requirements. Similar to many other models (for example, Di Toro et al., 1971; Park et al., 1974, 1975, 1979, 1980; Lehman et al., 1975; Canale et al., 1975, 1976; Thomann et al., 1975, 1979; Scavia et al., 1976; Bierman et al., 1980; O'Connor et al., 1981), AQUATOX uses the Steele (1962) formulation for light limitation. Light is specified as average daily radiation. The average radiation is multiplied by the photoperiod, or the fraction of the day with sunlight, based on a simplification of Steele's (1962) equation proposed by Di Toro et al. (1971). The equation is slightly different when the model is run with a daily versus an hourly time-step: e Photoperiod (LtAtDepthn ., -LtAtTopn ., ) PeriphytExt LtLimit Dailv = 0.85 ^-^ ^^ Extinct (Depth Bottom - DepthTop ) ,,,. ., = g (LtAtDepthHourly - LtAtTopHourly)- PeriphytExt LiLimii Hourl I-*") Extinct ( Depth Bottom-DepthTop) where: LtLimitTimeStep = light limitation (unitless); e = the base of natural logarithms (2.71828, unitless); Photoperiod = fraction of day with daylight (unitless), see (26); Extinct = total light extinction (1/m), see (40), (41); DepthBottom = maximum depth or depth of bottom of layer if stratified (m); if periphyton or macrophyte then limited to euphotic depth; Depthrop = depth of top of layer (m); LtAtTop = limitation of algal growth due to light, (unitless) see (44), (45); LtAtDepth = limitation due to insufficient light, (unitless), see (43); PeriphytExt = extinction due to periphyton; only affects periphyton and macrophytes (unitless), see (42). Because the equation overestimates by 15 percent the cumulative effect of light limitation over a 24-hour day, a correction factor of 0.85 is applied to the daily formulation (Kremer and Nixon, 1978). When AQUATOX is run with an hourly time-step, the correction factor of 0.85 is not relevant, nor the inclusion of photoperiod. Light limitation does not apply to free-floating macrophytes as these are assumed to be located at the surface of the water. Even when the model is run with an hourly time-step, two algal equations utilize the daily light limit equation (38) as most appropriate. First, when calculating algal mortality, the stress factor for suboptimal light and nutrients (68) is expecting the input of daily light limitation (i.e. the plants do not all die each night). Secondly, when calculating the sloughing of benthic algae (75) 69 ------- AQUATOX (RELEASE 3.1) TECHNICAL DOCUMENTATION _ CHAPTER 4 the calculation of suboptimal light is calibrated to daily light limitation, not the instantaneous absence or presence of light (i.e. sloughing is not more likely to occur when it is dark). Extinction of light is based on several additive terms: the baseline extinction coefficient for water (which may include suspended sediment if it is not modeled explicitly), the so-called "self- shading" of plants, attenuation due to suspended parti culate organic matter (POM) and inorganic sediment, and attenuation due to dissolved organic matter (DOM): Extinct = WaterExtinction + PhytoExtinction + ECoeffDOM DOM + ECoejfPOM-l |