MODELING WATER QUALITY AND BIOTA IN THE COLORADO RIVER BELOW AUSTIN, TEXAS
APPROVED BY DISSERTATION COMMITTEE:
MODELING WATER QUALITY AND BIOTA IN THE COLORADO RIVER BELOW AUSTIN, TEXAS
by
MICHAEL ERNEST BARBER, 8.5.C.E., M.S.C.E.
DISSERTATION
Presented to the Faculty of the Graduate School of
The University of Texas at Austin
in Partial Fulfillment
of the Requirements
for the Degree of
Doctor of Philosophy
THE UNIVERSITY OF TEXAS AT AUSTIN
December, 1991
To my family and friends
ACKNOWLEDGEMENTS
The author would like to express his appreciation to Dr. Neal E. Armstrong, Professor of Civil Engineering at The University of Texas, for the advice and guidance provided throughout the completion of this dissertation. The author would also like to thank his dissertation committee, Dr. Randall J. Charbeneau, Dr. Edward R. Holley, and Dr. George H. Ward, Jr., and other faculty such as Dr. John M. Sharp, Jr., for their review of this dissertation and their many helpful suggestions.
A very special note of thanks must go to my wife, Aleisa, and my parents, Joseph and Josephine, for their never ending support, understanding, and encouragement throughout my entire graduate program. The words "I Love You" don't seem nearly enough, yet I know of no others to say in their place.
I have had the opportunity to make a lot of great friends along the way. Friends like Wade and Rebecca Hathhorn, Chuck Downer, Su, Phil, Bao, Bruce and George helped make the journey a pleasant one and for that I thank them. And of course, a very special ’’skoal” goes out to my pseudo-brother Mark Adley who, through all the years and all the miles, has always been there willing to lend a sympathetic ear or share in the excitement. There will always be a special place in my heart for all of you.
This research was funded primarily through the gracious support of The City of Austin, Texas via contract number CIP #237435 G 87382. Computer facilities were provided by The Center for Research in Water Resources at The University of Texas at Austin.
This dissertation is dedicated in loving memory to my grandfather, Ernest John Canfield, who taught me, in a roundabout sort of way, to aim high at my goals otherwise I might miss and hit the hood. Papa never missed.
MODELING WATER QUALITY AND BIOTA IN THE COLORADO RIVER BELOW AUSTIN, TEXAS
Publication No.
Michael Ernest Barber, Ph.D. The University of Texas at Austin, 1991
Supervisor: Neal E. Armstrong
A model, capable of performing both steady state and dynamic analyses, has been developed and used to predict changes in potential population densities of submerged rooted aquatic vegetation (SRAV) and epiphytic algae as well as the concentration changes of other water quality parameters associated with eutrophication problems. The model has been applied to the Colorado River between Austin and Bastrop, Texas to evaluate the impacts on the downstream river ecosystem due to potential changes in waste water treatment levels and increased waste water discharges caused by population growth.
Nine steady state and twelve dynamic simulations were performed. Although phosphorus removal from Austin's waste water treatment plants reduced downstream phosphorus concentrations significantly, the SRAV growth potentials were not reduced proportionately. Similarly, the incremental difference in SRAV densities as a result of phosphorus and nitrogen removal suggested that
removal of both nutrients is also unnecessary.
The evidence does suggest a correlation between the level of waste water treatment and the growth potential of epiphytic algae. The largest reductions in growth potentials were due to phosphorus and nitrogen removal. However, less than half of the reductions were attributable solely to phosphorus removal.
The results of the analyses indicate that the constant rhetoric claiming phosphorus removal is need to reduce SRAV growth in the lower Colorado River is incorrect. Furthermore, reductions of epiphytic algae potentials due to the removal of phosphorus alone were relatively small for several scenarios. Consistently significant reductions in growth potentials were only obtained by phosphorus and nitrogen removal.
These results indicate that further scientific investigation is needed to determine the true “detrimental” effects of SRAV and epiphytic algae on the Colorado River before a cost/benefit analysis can be completed. Accordingly, the present conclusion would be not to expend the money and resources upgrading Austin's waste water treatment facilities solely on the basis of SRAV and epiphytic algae control.
Table of Contents
Page
ACKNOWLEDGEMENTS iv
ABSTRACT vi
LIST OF FIGURES xiii
LIST OF TABLES xvii
1.0 INTRODUCTION 1
1.1 Need for Study 1
1.2 Objectives 6
1.3 Scope 8
1.4 Uniqueness of Research 11
2.0 LITERATURE REVIEW 14
2.1 Introduction 14
2.2 Sources/Sinks of Nutrients 14
2.2.1 External Sources/Sinks 17
2.2.2 Internal Sources/Sinks 22
2.3 Nutrient Transport 31
2.3.1 Water Column Advection 31
2.3.2 Role of Sediment Transport 34
2.3.2.1 Sources of Sediment 35
2.3.2.2 Methods of Estimating Yield 38
2.3.3 Sediment Transport Theory 43
2.3.3.1 Suspended Sediment Equations 45
2.3.3.2 Total Sediment Load Equations 52
2.4 Summary of Existing Water Quality Models .. 60
2.5 Previous Water Quality Model of the Colorado River 63
Table of Contents (Continued)
2.5.1 Study Area Under Consideration 63
2.5.2 Norman Johns' Pre-WASP3 Model 64
2.5.3 Elizabeth Krousel's QUAL2E Model .... 66
3.0 MODEL DEVELOPMENT 88
3.1 Original EUTRO4 Program Development 88
3.1.1 Finite-Difference Solution of Mass Balance Equation 88
3.1.2 Applicability of One-Dimensional Model Formulation 94
3.2 Modifications to EUTRO4 Program 96
3.2.1 Final Nutrient-Oxygen-BOD- Phytoplankton Kinetic Equations 97
3.2.2 Macrophyte and Epiphytic Algae Kinetic Equations 112
3.3 Development and Calibration of Unsteady DAFLOW Model 121
3.3.1 Theory of Diffusion Analogy Model ... 122
3.3.2 Calibration of DAFLOW Program 127
3.3.3 Verification of DAFLOW Program 132
4.0 INPUT PARAMETER ESTIMATION 146
4.1 Need for Parameter Estimation 146
4.2 Water Column Advection for Hydrodynamic
Flow Simulations 146
4.2.1 Coefficients Required for Dynamic Flow Routing 147
4.2.2 Estimation of Ungaged Tributary Inflows 148
4.2.3 Estimation of Ground Water Inflows .. 149
4.3 Selection of Dispersion Coefficients 150
4.3.1 Advective Diffusion 151
4.3.2 Molecular Diffusion 151
4.3.3 Turbulent Diffusion 153
Table of Contents (Continued)
4.3.3.1 Longitudinal Dispersion in the Water Column 153
4.3.3.2 Pore Water Diffusion 157
4.3.4 Numerical Dispersion 159
4.4 Nutrient and Constituent Loadings 163
4.4.1 Initial Constituent Concentrations .. 163
4.4.2 Boundary Conditions 168
4.4.2.1 Under Steady State Flow Conditions 168
4.4.2.2 Under Unsteady Flow Conditions 169
4.4.3 External Sources 170
4.4.3.1 Historical Loadings 170
4.4.3.2 Projected Loadings 171
4.4.4 Internal Sources 172
4.5 Settling and Resuspension Velocities of Solids Systems 172
4.6 Coefficients and Constants for Epiphytic Algae, Macrophyte Kinetics, and Water Quality Variables 178
5.0 EXAMPLE PROBLEM 194
5.1 Description of Physical System 194
5.2 Input Requirements for Example Problem .... 195
5.3 Discussion of Results 198
5.3.1 Effect of Finite-Difference Solution Technique on Results 200
5.3.2 Effect of Time Step on Model Stability and Numerical Dispersion .. 201
6.0 SIMULATION RESULTS AND DISCUSSION 211
6.1 Calibration of Final EUTRO-M Model 211
6.1.1 Scenario Examined 212
Table of Contents (Continued)
6.1.2 Selection of Benthic Sediment Layer Properties 212
6.1.3 Comparison of Results 214
6.2 Model Verification 219
6.3 Discussion of Loading Scenarios 221
6.4 Steady State Flow Simulations for Alternative Loading Scenarios 223
6.4.1 Scenario 1 - Current Waste Water Treatment Plant Discharges 224
6.4.2 Scenario 2 - Projected Year 2000 WWTP Discharges 225
6.4.3 Scenario 3 - Permitted Year 2000 WWTP Discharges 227
6.4.4 Comparison of Steady State Scenarios 229
6.5 Steady State Sensitivity Analysis 230
6.5.1 Analysis of Model Coefficients 231
6.5.2 Effects of Ground Water Inflow 232
6.6 Hydrodynamic Flow Simulation 233
6.6.1 Simulation of Unsteady Flow 234
6.6.2 Calibration of Unsteady Water Quality Model 235
6.7 Evaluation of Unsteady Flow Simulations for Alternative Loading Scenarios 238
6.7.1 Scenario 1 - Current Waste Water Treatment Plant Discharges 239
6.7.2 Scenario 2 - Projected Year 2000 WWTP Discharges 240
6.7.3 Scenario 3 - Permitted Year 2000 WWTP Discharges 241
6.7.4 Scenario 4 - Permitted Year 2000 WWTP Discharges with a Minimum Instream Flow Requirement at Longhorn Dam 241
Table of Contents (Continued)
7.0 SUMMARY, CONCLUSIONS AND RECOMMENDATIONS 302
7.1 Summary of Analysis 302
7.2 Conclusions 303
7.3 Recommendations 306
APPENDICES
A. Listing of Modified WASPB Subroutine 311
B. Listing of Modified EUTRO Common Blocks ... 343
C. Example EUTRO-M Input File 351
D. Example EUTRO-M Output Files 363
E. Listing of Miscellaneous Support Programs . 394
REFERENCES 4 06
List of Figures
Number Title Page
2.2.1 Phosphorus Cycle in Rivers 67
2.2.2 Nitrogen Cycle in Rivers 68
2.2.3 Land Use versus Phosphorus Loading 69
2.2.4 Land Use versus Nitrogen Loading 70
2.2.5 Potential Nutrient Concentration Gradients in a Benthic Sediment Layer 71
2.3.1 Effect of SRAVs on Stream Velocity 72
2.3.2 Typical Suspended Sediment Concentration Profile 73
2.5.1 Location of Colorado River Basin within Texas 74
2.5.2 Generalized Map of Colorado River Study Area 75
2.5.3 Colorado River Between Austin and Bastrop . 76
2.5.4 Calibration Run for Pre-WASP3 Model 77
2.5.5 Predicted Potential Macrophyte Growth Using Pre-WASP3 Model 78
2.5.6 Predicting Potential Macrophyte and Epiphytic Algae Densities using QUAL2E Model 79
3.1.1 Notation Schematic for Finite-Difference Equation 134
3.2.1 Kinetics Linkages Between the Ecological Components 135
3.2.2 EUTRO-M State Variable Interactions 136
Number Title Page
3.2.3 Representation of Multi-Layer Water/Sediment System 137
3.2.4 Characteristic Phosphorus and Nitrogen Growth Coefficient Factors for Aquatic Vegetation 138
3.2.5 Characteristic Light and Temperature Growth Coefficient Factors for Aquatic Vegetation 139
3.3.1 Typical Time of Travel Hydrograph 140
3.3.2 Upstream Boundary Condition for DAFLOW Model Calibration Run 141
3.3.3 Original Statistical Calibration of DAFLOW Model 142
3.3.4 Modified Calibration Using Flows Greater Than 1,500 Cubic Feet Per Second 143
3.3.5 Upstream Boundary Condition for DAFLOW Model Verification Run 144
3.3.6 Verification of DAFLOW Model 145
4.2.1 Schematic of Point Source Inflows within the Study Area 180
4.3.1 Variation in Longitudinal Dispersion Coefficients Below Longhorn Dam 181
4.5.1 Relationship Between Stream Velocity, Particle Size, and Sediment Regimes 182
4.5.2 Particle Size Analysis of Suspended Sediment in Colorado River 183
4.5.3 Resuspension Rates as a Function of Bottom Shear Stress 184
5.1.1 Physical Layout of Example Problem 203
5.2.1 Schematic of Input Groups 204
Number Title Page
5.3.1 Effect of Finite-Difference Solution Procedure 205
5.3.2 Effect of Time Step on Numerical Dispersion 206
5.3.3 Concentration Oscillations as a Result of Model Instability 207
6.1.1 Calculated Versus Measured Concentrations and Densities for Steady State Calibration Simulation 243
6.2.1 Calculated Versus Measured Concentrations and Densities for Steady State Verification Simulation 250
6.4.1 Changes in Scenario 1 Concentrations Due to Loading Cases 256
6.4.2 Percent Changes in Maximum SRAV and Epiphytic Algae Growth Potentials for Scenario 1 262
6.4.3 Comparison of Predicted SRAV and Epiphytic Algae Growth Potentials for Scenario 1 .... 263
6.4.4 Changes in Scenario 2 Concentrations Due to Loading Cases 264
6.4.5 Percent Changes in Maximum SRAV and Epiphytic Algae Growth Potentials for Scenario 2 268
6.4.6 Comparison of Predicted SRAV and Epiphytic Algae Growth Potentials for Scenario 2 .... 269
6.4.7 Changes in Scenario 3 Concentrations Due to Loading Cases 270
6.4.8 Percent Changes in Maximum SRAV and Epiphytic Algae Growth Potentials for Scenario 3 275
Number Title Page
6.4.9 Comparison of Predicted SRAV and Epiphytic Algae Growth Potentials for Scenario 3 .... 276
6.4.10 Scenario Comparison of Predicted SRAV and Epiphytic Algae Growth Potentials for Case 1 277
6.4.11 Scenario Comparison of Predicted SRAV and Epiphytic Algae Growth Potentials for Case 2 278
6.4.12 Scenario Comparison of Predicted SRAV and Epiphytic Algae Growth Potentials for Case 3 279
6.6.1 Typical Daily Flow and Depth Variations in the Colorado River Below Longhorn Dam 280
6.6.2 Hourly July 1986 Discharges from Longhorn Dam 281
6.6.3 Calculated Versus Measured Gage Flows at Bastrop, Texas 283
6.6.4 Calculated Versus Measured Concentrations and Densities for Dynamic Calibration Simulation 285
6.6.5 Variations in Ammonia Concentrations with Time 291
6.7.1 Percent Changes in SRAV and Epiphytic Algae Growth Potentials Under Dynamic Flow Conditions 292
List of Tables
Number Title Page
2.2.1 Nutrient Concentrations in Domestic Waste Water 80
2.2.2 Nonpoint Source Loadings 81
2.2.3 Regional Nutrient Loadings 82
2.2.4 Critical Scour Velocities for SRAV 83
2.3.1 Relationship of Particle Size and Sediment Class 84
2.3.2 Reliability of Sediment Discharge Equations 85
2.4.1 Comparison of Nutrient Models 86
2.4.2 Comparison of Conventional Pollutant Models 87
4.2.1 Estimated Inflows from Ungaged Drainages .. 185
4.5.1 State Variables and Their Associated Solids Systems 186
4.6.1 Program Constants Used in EUTRO-M Steady State Analysis 187
4.6.2 Definition of Steady State Macrophyte and Epiphytic Algae Constants 190
5.3.1 Output Display Variables for EUTRO-M 208
6.3.1 Summary of Loading Scenarios for Steady State and Dynamic Flow Simulations 293
6.3.2 Walnut Creek WWTP Concentrations for Loading Cases 294
6.3.3 South Austin Regional WWTP Concentrations for Loading Cases 295
Number Title Page
6.3.4 Govalle WWTP Concentrations for Loading Cases 296
6.3.5 Bergstrom WWTP Concentrations for Loading Cases 297
6.4.1 Comparison of Average Steady State Aquatic Vegetation Growth Potentials 298
6.5.1 Summary of Steady State Sensitivity Study . 299
6.5.2 Summary of Steady State Sensitivity Analysis of Hydraulic Properties 300
6.7.1 Comparison of Percent Reduction in Dynamic Aquatic Vegetation Growth Potentials 301
CHAPTER 1 INTRODUCTION
1.1 NEED FOR STUDY
Despite continuing advances, a number of water quality problems still exist in our nation's rivers and streams. These problems are caused by excessive nutrient concentrations, high organic loads and lowered dissolved oxygen levels, high temperature, and hazardous or toxic materials. The problems associated with these types of pollution range from purely aesthetic nuisances such as bad odors to concentrations of substances that are lethal to aquatic life or even humans.
Nutrient concentrations are of concern mainly because of the role they play in the eutrophication process. Millions of dollars are spent each year trying to reduce or eliminate nutrient loadings from lakes and rivers. Of the various nutrients found in surface water, nitrogen and phosphorus have long been identified as the two nutrients involved in the growth of aquatic vegetation. Even though phosphorus has generally been considered to be the limiting nutrient in freshwater systems (Clark, Viessman, and Hammer, 1977; Stumm and Morgan, 1981), nitrogen has been found to be as or more limiting in Texas lakes, reservoirs, and rivers (Armstrong et al., 1987, 1989). Similar results have been found in other southern states. As a result, this research effort focuses on both nitrogen and phosphorus as they affect surface water ecosystems.
The growth of macrophytes, specifically submerged rooted aquatic vegetation (SRAV), have been of concern in river basins and have been thought to be an indication of eutrophication. Unlike phytoplankton and non-rooted macrophytes, SRAVs obtain most of their nutrients from the bottom sediments (Armstrong et al., 1987; Gunnison and Barko, 1989). This is an important distinction whose consequences must be fully understood in order to comprehend the fallacy of SRAV reduction by nutrient removal. Nutrients which adsorb to suspended sediment particles settling out of the water column provide one of several continuing sources of nourishment for the plants. In an effort to minimize ’’excessive” SRAV growth, a popular solution advocated by many engineers and scientists has been to reduce the point source nutrient loadings. This has prompted many environmental groups to call for tertiary treatment at new and existing wastewater treatment plants throughout the country.
However, studies indicate that a certain percent reduction in nitrogen and/or phosphorus loadings from a wastewater treatment facility does not necessarily translate into an equivalent percent reduction in SRAV growth (Armstrong et al., 1987; Johns, 1988). A number of studies have shown that macrophyte growth is also a function of nitrogen concentrations, light availability, temperature, substrate type, and stream currents (Thomann and Mueller, 1987; Wright and McDonnell, 1986; Barko and Smart, 1983). Consequently, it is important not only to understand the transport and fate of phosphorus, but also to comprehend the role of the other parameters in limiting plant growth. These complex interactions occurring within a river ecosystem make it difficult to predict the outcome of a proposed action prior to its implementation.
Another similar problem exists in many river basins due to the perceived problem involving epiphytic or attached algae. Epiphytic algae, sometimes referred to in the literature as periphyton, benthos, littoral or sessile algae (Butcher, 1940; Klasvic, 1974) are attached to substrate in the river such as rocks, bridge piers, or even other plants but obtain their nutrients from the water column much like phytoplankton and nonrooted macrophytes.
As an example of the potential problems associated with SRAVs and epiphytic algae, consider the long running debate between the City of Austin (Austin) and the Lower Colorado River Authority (LCRA) regarding discharges from Austin's wastewater treatment plants. LCRA contends that upgrading Austin's treatment facilities from secondary to tertiary plants would substantially reduce the amount of rooted macrophyte and epiphytic algae growth in the Colorado River between Austin and Bastrop. The estimated 1988 cost for renovating Austin's plants was approximately 126 million dollars (Johns, 1988) without any guarantee that such remediation measures would be successful. In fact, preliminary investigations indicated that tertiary treatment would not satisfactorily fulfill the intended goals. In this time of strained financial resources, it is imperative that money spent on controlling waste discharges be allocated prudently.
Advanced steady state methodologies that enable all of these interactions to be accurately modeled simultaneously exist in the model used by Armstrong et al. (1987) to model water quality and biota in the Colorado River below Austin. However, even this model did not take into account directly the effect of the benthic sediment layer on nutrient retention or phytoplankton kinetics. Furthermore, unsteady state models which combine the latest eutrophication kinetics and benthic layer capabilities into a model with macrophyte and epiphytic algae growth kinetics do not currently exist. Moreover, most existing methods attempt to solve only small portions of these problems using steady state assumptions rather than non-steady state assumptions, leaving the user to perform manual dynamic approximations.
Most of the newer models which report to be dynamic are, in fact, basically pseudo steady state. Pseudo steady state means that one or more of the processes included in the model are developed using the time derivative but other processes are handled using steady state assumptions. This type of model is prevalent in sediment transport models where unsteady typically means with respect to water discharge rather than sediment transport (Dawdy and Vanoni, 1986). Although this procedure may be valid for representing some parameters (such as sediment transport), caution should be exercised to insure that the averaging process does not adversely affect the results. As a hypothetical example, consider a case where the light extinction factor is averaged over the entire year. Time-averaging of this variable could lead to erroneous predictions if the growing season for macrophytes occurred primarily during the runoff season where high sediment concentrations might act to reduce light penetration.
There is a need to develop and combine dynamic nitrogen/phosphorus-sediment-SRAV/algae methodologies into a single model. The result would be a model capable of evaluating the effects of proposed developments, design guidelines, impoundments, runoff controls, or treatment plant expansions on phosphorus and nitrogen concentrations, macrophyte and algae growth rates, and sediment transport prior to their construction or implementation. The results of such studies could be used to generate reliable cost/benefit analyses.
In addition to solving questions currently being asked by agencies such as the City of Austin, the Lower Colorado River Authority, and the Texas Water Commission, the nitrogen/phosphorus-sediment-SRAV/algae methodologies would provide a vital link in bridging the gap between water resources and environmental engineering by combining state-of-the-art hydraulichydrologic concepts with unsteady water quality kinetics routines designed to answer both hydraulic and environmental concerns. The basic model could be applied to a number of similar watersheds and could also act as a base for further work in the areas of toxics and hazardous waste transport.
1.2 OBJECTIVES
The overall objective of this research is to develop a model capable of dynamically determining river bed and water column nitrogen and phosphorus concentrations in river basins where SRAVs and epiphytic algae are of concern. To accomplish this task, the following specific objectives were completed:
(1) Modified the original water quality model formulation from the Armstrong et al. (1989) study to be consistent with the latest version of the USEPA water quality model (EUTRO4). Included additional features such as a benthic sediment layer and phytoplankton as a state variable in the model. Converted the macrophyte, algae, and suspended sediment routines developed for the original model. Added options such as pore water transfer, ground water seepage, new macrophyte scour methodologies, and a refined light extinction equation.
(2) Performed calibration of the updated steady state EUTRO4 model by comparing the predicted results with those computed by the original model and the field data collected as part of the original City of Austin project.
(3) Verified the steady state water quality model using field measurements from data collected during similar flow and temperature conditions the previous year.
(4) Examined and evaluated the consequences of different nutrient loadings with special attention given to the growth rates of SRAVs and epiphytic algae. Performed detailed sensitivity analyses as part of this task.
(5) Developed programs to link the USGS DAFLOW dynamic flow model to the EPA EUTRO4 model and added a refined suspended sediment algorithm to the water quality model so that unsteady state conditions could be modeled with greater accuracy. Recalibrated the dynamic model. Used the calibrated model to evaluate different treatment and waste water discharge scenarios.
While somewhat general in nature, the methodologies and model derived as part of this research were calibrated and tested using data from the lower Colorado River drainage basin. The mainstem reach of the river was modeled from Longhorn Dam to the City of Bastrop in order to examine the effects of discharges from the City of Austin's wastewater treatment facilities on nitrogen and phosphorus concentrations. The results of this portion of the research can be used in the future to help develop a sound management strategy for nitrogen and phosphorus concentrations in the lower Colorado River basin.
1.3 SCOPE
Completion of the objectives outlined above was accomplished in two main phases. Phase I involved updating the previous pre-WASP3 steady state model to WASP 4, specifically EUTRO4, by incorporating a benthic sediment layer, and by adding suspended sediment as a system variable. The inclusion of factors such as a benthic sediment layer and phytoplankton as a state variable resulted in a significantly different model calibration. As a result, the pre-WASP3 and EUTRO4 nutrient concentrations, SRAV growth potentials, and algae densities were different.
Phase II of the study required the addition of a hydrodynamic segment to the model in order to accurately predict the long-term effects of nutrient loadings on the lower Colorado River. The model took full advantage of new program features by including options such as pore water exchange, ground water inflows, and updated nutrient loading and macrophyte scour techniques. A brief discussion of the proposed Phase II work effort is provided below. Details regarding the specific tasks are provided in Chapters 3 and 4.
The first step of Phase II was to replace the average values of the nutrient loadings used as initial conditions and boundary conditions by Johns (1988) with the dynamic loadings calculated using a more accurate method. The new method of estimating nutrient loading involves coupling the nutrient concentrations with the Soil Conservation Services' procedure for determining storm water runoff. The EUTRO4 model formulation also included pore water exchanges and ground water inflows.
Next, methodologies developed in the original study (Johns, 1988) for modeling the internal sources and sinks of nitrogen and phosphorus, caused primarily by the growth and senescence of macrophytes and epiphytic algae and the settling and resuspension of particulate nutrients, were combined with the scour equations developed by Vecchio (1989) and incorporated into the EUTRO4 model.
The next step was to include a suspended sediment routine in the water quality model. This was accomplished by using a number of sediment transport rating curves each based on the least squares relationship which best fit the data available from the USGS and other studies. The model allows the user to specify a different rating curve for each surface water segment. A one-dimensional hydrodynamic model was used to establish flow velocities. The amount of sedimentnutrient settling and/or resuspension was then estimated.
The concentrations of nitrogen and phosphorus in the river bed (both in the particulate and interstitial water phases) are also functions of the direction and quantity of groundwater movement into or out of the river. The research method must therefore be able to account for infiltration and exfiltration of groundwater. This was accomplished by examining the difference between the predicted and measured downstream gage flows and assuming the net difference was distributed uniformly along the length of the river bed.
Some of the nutrients which enter the river are used by the epiphytic algae and the phytoplankton. A portion of this aquatic plant life may settle out to become part of the sediment concentration; however, the remainder results in a net loss of nitrogen and phosphorus in the system. The kinetic equations in the model were modified to include both the transfer of nutrients to the sediment and the system loss due to the flushing effect.
The amount of sediment which accumulates in or was scoured from each river segment over a period of time and the appropriate partitioning coefficient were combined with the advective and diffusive components to calculate nutrient concentrations. These concentrations were combined with temperature and light extinction coefficients so that submerged rooted aquatic vegetation (SRAV) and epiphytic algae growth potentials could be calculated.
The completed model was calibrated and verified using data from the lower Colorado River. The calibrated model was then used to produce month long simulations of the river basin which indicated the response of the river to changes in nutrient loadings.
Field measurements or empirical estimates of state variables such as nitrogen and phosphorus concentrations, SRAV growth, epiphytic algae densities, sediment transport, and nonpoint source nitrogen and phosphorus loadings typically demonstrate a wide degree of scatter and uncertainty when compared to the values predicted by the governing kinetics equations. To account for this uncertainty, a thorough sensitivity analysis was conducted using the calibrated steady state model.
1.4 UNIQUENESS OF PROPOSED RESEARCH
The goal of this research was to use these widely accepted techniques to create a state-of-the-art dynamic water quality model. On a macro scale, the coupling of the hydrodynamic water quality model containing macrophyte and epiphytic algae components into a single model is an advancement of the current procedures. Several examples are needed here to explain the deficiencies of the existing technology. First, while a number of hydrodynamic water quality models have been developed in the past, the application of these models to riverine ecosystems concerned with eutrophication problems has been quite limited. Modeling of river systems is often simplified to steady state problems and the inclusion of rooted and attached macrophyte kinetics is not included.
Second, most sediment transport models have concentrated on the hydraulic aspects of the problem rather than the environmental problems. Consequently, total annual sediment yields into a reservoir are usually considered to be sufficient. The propensity of the suspended sediment to carry nutrients or other constituents is often ignored.
Finally, although some water quality models claim to handle sediment transport, their treatment of the process is usually quite simplified. Contaminant partitioning between sediment fractions is not allowed, light availability because of turbidity is averaged over the entire year, or the actual downstream transport of particles is ignored.
In addition to the novelty of the overall procedure mentioned above, there are several key elements of the research project which provided improvements to the previous solution procedures. These individual tasks can be described as follows:
(1) The proposed model includes parameters which incorporate rooted macrophyte and attached algae growth, scour, and senescence into a hydrodynamic transport model. The long-term effect of changing nutrient concentrations and loadings on submerged rooted aquatic vegetation was examined. Also, the effectiveness of using high flows to control macrophyte growth can be studied.
(2) A methodology for linking macrophyte growth to the nutrient concentrations in the benthic sediment layer and to the light extinction and scour coefficients in the water column was added so that better predictions of the aquatic vegetation growth potentials could be determined under unsteady hydraulic conditions.
(3) The model was modified to provide a refined procedure for linking the concentrations of
suspended sediment, phytoplankton and organic particulate matter with light extinction coefficients in order to determine the sensitivity of the system to light.
The completion of these tasks advanced the field of water quality modeling by eliminating the need for gross approximations in several critical areas. The above summary illustrates the uniqueness and need for improvement in each of the prescribed areas. A more thorough discussion of the consequences and importance of each of these developments is provided in the following chapters.
CHAPTER 2 LITERATURE REVIEW
2.1 INTRODUCTION
Since aquatic plant growth is a function of nutrient concentrations, light availability, and temperature, excessive aquatic plant growth is often considered to be an indicator of eutrophic conditions within a body of water. The photosynthetic and respiratory cycles of plants cause wide fluctuations in the dissolved oxygen concentrations which may be detrimental to the biota of the receiving water. The following literature review summarizes the findings of a number of studies concerning nutrient concentrations and light availability and describes several widely used methods for analytically modelling the natural systems.
2.2 SOURCES/SINKS OF NUTRIENTS
Nutrient loadings are a function of the density of population and livestock, methods or intensity of fertilization, type of cultivation (e.g., forests, grassland, cropland), the pedological characteristics of the soil, topography, climate and rainfall, and type of waste water treatment systems used in the area. While several elements are required for the growth of aquatic plants, the two most prevalent in eutrophication studies are phosphorus and nitrogen. Consequently, this literature review will also focus on these two elements when discussing the role of nutrients.
Total phosphorus is the measure of all the phosphorus present in both the dissolved and particulate phases. Within the dissolved and particulate phases are several forms of phosphorus. Metcalf and Eddy (1979) reported that the three most common forms of phosphorus present in surface waters are:
1. organic phosphorus 2. orthophosphates (H 2 PO 4 \ hpo 4 - 2 z po 4 - 3 ) 3. polyphosphates (Na 3 (PO 3 ) 6 , Na 5 P 3 O 10 , Na 4 P 2 0 7 )
However, not all of these phosphorus forms can be directly used by aquatic plants such as phytoplankton. Of these, only the dissolved orthophosphate ion is used in the photosynthesis process by the phytoplankton or macrophytes (Clark, Viessman, and Hammer, 1977). In the bed sediments, phosphate is supplied to the interstitial water through desorption of phosphate from ferrous iron (Ishikawa and Nishimura, 1989).
A schematic diagram of the phosphorus cycle is presented in Figure 2.2.1. As shown in the diagram, the other forms of phosphorus can undergo reactions which either convert the phosphorus to orthophosphate or to an insoluble precipitate.
The other key element is nitrogen. Like phosphorus, there are also several forms of nitrogen present in natural bodies of water. The most common forms are:
1. Ammonia (NH 3 ) 2. Nitrate (NO 3 ) 3. Nitrite (NO 2 ) 4. Organic Nitrogen
These four forms are referred to as total nitrogen (Thomann and Mueller, 1987). Ammonia, nitrate and nitrite are called total inorganic nitrogen and can be used by aquatic plants for growth. Total Kjeldahl nitrogen, TKN, is a combination of ammonia nitrogen and organic nitrogen and should not be confused with the term ’’total nitrogen”. Nitrite, in natural systems, does not occur in large concentrations for long periods of time. As a result, nitrite is not always treated as a separate state variable.
Although both ammonia and nitrate can be used by the aquatic plants, for physiological reasons, the ammonia phase is the form generally preferred by the aquatic vegetation (Ambrose, et al., 1988). Figure 2.2.2 illustrates a schematic of the nitrogen cycle. Notice that ammonia, in the presence of oxygen and nitrifying bacteria (Nitrosomonas and Nitrobacter), can be converted into nitrate. This process is called nitrification.
Denitrification can also occur in a natural system. This process involves the conversion of nitrate (NO 3 ) into N 2 and other nitrogen byproducts such as N 2 O and NO. Denitrification requires anaerobic conditions so that the organisms responsible can use N0 3 as an electron acceptor. Consequently, denitrification should occur
only in the sediment bed or under extremely stressed surface water conditions.
Because of their important role in the photosynthesis process, limiting nutrient loadings to a system may be used in some instances to control eutrophication. To properly estimate the impact of this action, it is first necessary to identify and quantify the potential sources and sinks of nutrients. The two broad categories of nutrient sources and sinks are called external and internal.
2.2.1 External Sources/Sinks
External refers to those sources/sinks of nutrients derived from outside the water body. The primary external sources of nutrients as described by Thomann and Mueller (1987) are:
1. Municipal wastes 2. Industrial wastes 3. Agricultural runoff 4. Urban and suburban runoff 5. Forest runoff 6. Atmospheric fallout.
These topics can be further defined as point and nonpoint sources. As the names imply, point sources are those which enter the system at a well defined point (such as a waste water discharge pipe) and nonpoint sources are those which enter the system over a diffuse area. In general, nitrogen and phosphorus loadings come from both point and nonpoint sources.
Point Sources
Municipal and industrial wastes discharges are the ones most frequently associated with the term point source although recently there has been a good deal of discussion concerning the possibility of treating some urban runoff as point sources. Point sources receive a lot of attention because they are relatively easy to identify and are usually the most cost effective remediation solution. Typical medium strength domestic waste water contains approximately 40 mg/1 of total nitrogen and 10 mg/1 of total phosphorus before treatment (Krenkel and Novotny, 1980). Table 2.2.1 reflects the significance of various types of treatment on final nutrient concentrations.
Nonpoint Sources
Nonpoint sources, derived primarily from urban, agricultural, and rural runoff, may be more difficult to quantify. This process is made more difficult by recent studies that indicate phosphorus settles out with the sediment under low flow conditions and is quickly resuspended during high flow conditions caused by storm runoff (Dorioz et al., 1988; Dong et al., 1983). On average, approximately 50 percent of the total phosphorus concentrations and an even greater percentage of the nitrogen loading originate from nonpoint sources (Novotny and Chesters, 1981). Of course, the percentages for any particular system may vary substantially from these averages. Much of the total phosphate from nonpoint sources enters the stream in particulate form. However, as discussed in Section 2.2, only dissolved phosphate can be readily used by the algae and the other aquatic plants. Phosphorus concentrations in solution are a function of adsorptiondesorption and precipitation-dissolution reactions.
The list of external sources contains three separate topics relating to runoff; agricultural, urban, and forest runoff. The reasoning behind these divisions is that each of the topics demonstrates different characteristics in regard to nutrient runoff. Novotny and Chesters (1981) report that nationwide studies found a distinct correlation between land use and nutrient concentrations in streams. As illustrated in Figure 2.2.3, phosphorus concentrations increase as the percentage of land used for agricultural and urban uses increases. Similar results can be seen in Figure 2.2.4 for nitrogen. Thomann and Mueller (1987) reported typical values for each component based on 1978 studies by the International Joint Commission and by Rast and Lee (1983). These values are presented in Table 2.2.2.
Using a regression model to relate runoff concentrations to agricultural and urban land use, Omernik (1977) developed nutrient loading information for three general regions. The data were obtained from a study performed as part of the National Eutrophication Survey conducted from 1972-1974 which involved 928 watersheds. The nitrogen and phosphorus concentrations determined by this study are summarized in Table 2.2.3. These results also indicate that nutrient concentrations increase with increases in agricultural and urban land uses.
A point must be made concerning the use of nutrient export coefficients. The amount of nutrients derived from surface runoff is highly correlated to rainfall frequency, duration, and intensity (Rast and Lee, 1983). Consequently, as pointed out by Rast and Lee (1983), the nutrient export coefficients reported in Tables 2.2.2 and 2.2.3 are ’’usually accurate within a factor of two” and thus should be replaced with measured loads whenever possible.
In a study of subsurface runoff from agricultural areas, Bottcher, Monke, and Huggins (1981) found that nitrates account for 75 percent of the nitrogen found in the runoff. Only 8 percent of the nitrogen was associated with the sediment. Soluble organic nitrogen and ammonia accounted for the remaining 17 percent. Conversely, the study found that approximately 70 percent of the phosphorus was associated with the sediments. The remaining 30 percent was divided nearly evenly between soluble inorganic and organic phosphorus phases.
The study also concluded that nutrient concentrations measured on or slightly after the first storm event following spring fertilization were abnormally high compared to other rainfall events. This was especially true for phosphorus adsorbed to sediment.
The process in which precipitation contributes to phosphorus loading (or any other pollutant loading) is often referred to as wet deposition or atmospheric fallout and should be considered as another potential external source of phosphorus. In undisturbed watersheds, aeolian sources of phosphorus (including both those caused by precipitation and dry fallout) can account for 50 percent or more of the phosphorus found in lakes (Galloway and Cowling, 1978). Galloway and Cowling (1978) also found that for developed basins, because of much higher concentrations due to man's activities, the percentage drops to approximately 3 percent (between 2 and 6 %). Novotny and Chesters (1981) site a Great Lakes study claiming that 26 percent of the phosphorus load is a result of precipitation.
In a compilation of several studies, Novotny and Chesters found that phosphorus concentrations in rain water typically ranged from 0.003 mg/L to 0.12 mg/L, although one study did report a value of 1.1 mg/L. The former values compare well to another study done on a lake in Ontario, Canada. The researchers here found the average concentrations of total phosphorus to be 0.026 mg/L with values ranging from 0.01 mg/L to 0.54 mg/L (Nicholls and Cox, 1978). Nicholls and Cox also found that the total atmospheric phosphorus load was 74.4 mg/m 2 ~yr. Of this, 39.0 mg/m 2 -yr was in the form of total dissolved phosphorus a portion of which, 20.6 mg/m 2 -yr, was in the dissolved orthophosphate form. The remaining 35.4 mg/m 2 -yr was in particulate form. The same study also found the total nitrogen loading to be 1600 mg/m 2 -yr.
Most all of the published studies deal with nutrient loadings to lake systems. However, there is no apparent reason why these values would be
inappropriate to use for river systems. The decision to include atmospheric fallout in a model should be based on the expected percent contribution to the waterway. If the contribution is small (as it usually is), it may not be necessary to add another source term due to the inherent uncertainty in the process.
2.2.2 Internal Sources/Sinks
Internal sources and sinks of nutrients are created by changes occurring within the system. Although internal nutrient sources can usually be traced back to some external supply, the retention time and physical nature of the instream sources allow them to be modeled as separate systems. The primary internal sources and sinks of nutrients are:
1. Sediment-nutrient fluxes 2. Phytoplankton growth 3. Macrophyte growth
Although numerous nutrient studies have been completed, investigations concerning the role of nutrients in the ecosystem of flowing waters are quite sparse. Most of the work done regarding sedimentnutrient relationships has involved lake samples. Still, the work completed to date appears to demonstrate that sediment-nutrient interaction in streams and rivers occurs in a complex two step process.
Studies indicate that the initial assimilation of soluble phosphorus is not due primarily to sediment adsorption. Two independent studies (Ball and Hooper, 1963; Kevern, Nelson and Wilham, 1966 unpublished) cited by Keup (1968) concluded that biological uptake of soluble phosphorus in river systems is due primarily to the periphyton (a bottom-covering film of microscopic plants). A similar study was also performed on a 65- mile reach of the South Platte River in Colorado. Keup describes the South Platte as a sandy-bottom stream with an autotrophic economy based on attached algae. In this study phosphorus was removed from the water in a logarithmic fashion with a decay constant of -0.011 per mile. However, the author was quick to point out that this value is quite variable from stream to stream.
Although the stream biota initially assimilate the nutrients from the water column, the periphyton do not play a significant role in the ultimate dispersal of the nutrients from the stream (McKee, et al., 1970). Bacterial decomposition of the organic material releases inorganic phosphorus which can either be taken up by other aquatic plants, or become bound to sediments, or precipitate out to the sediment layer, or changed into organic phosphorus by the bacteria. The assimilated nutrients bound to the sediments can then be transported downstream as part of the bed load (McKee, 1970; Keup, 1968). Some of the nutrients taken up by other aquatic plants may become part of the unmeasured floating debris which is transported downstream by advection.
Temporary bottom sediment storage can be significant during declining flow rates but long-term storage is limited due to degradation of most natural streams. In fact, both Dong, et al. (1983) and Harms, et al. (1978) found localized sharp phosphorus increases in the bottom sediments immediately downstream of point sources such as a secondary waste water treatment plant. Large increases in phosphorus concentrations at the beginning of a storm hydrograph are related to the scouring and resuspension of the previously deposited phosphorus bound sediments. However, even without scouring, some amount of nitrogen and phosphorus would be returned to the water column. The release of nutrients from the sediments is the result of a gradient in nutrient concentration between the water column and the interstitial water of the sediment (Thomann and Mueller, 1987). Figure 2.2.5 illustrates three hypothetical nutrient concentration gradients across a bed layer. In these examples, it can be seen that average nutrient concentrations over the entire bed layer are substantially different from the average concentration near the water/sediment interface.
The nutrient release process can be increased significantly if anaerobic conditions exist. When the overlying water column is anaerobic, the flux of phosphorus from the sediment increases significantly because of changes in the iron complexes at the sediment-water interface which result in increased diffusion. While the exact process is still unknown, it is theorized that phosphate, under oxidizing conditions and in the presence of iron, is precipitated on the sediment surface as insoluble ferric phosphate. When it is reduced (anaerobic conditions) soluble phosphorus is liberated back to the water column (McKee, et al, 1970).
This explanation has been quantitatively proven in a number of studies. Fillos and Swanson (1975), conducting laboratory experiments on river sediments, found that phosphorus release from the sediment increased from an average aerobic flux of 9.6 mg/m 2 -day to a maximum anaerobic flux of 96 mg/m 2 -day. Thomann and Mueller (1987) summarized the results of several studies which have shown dissolved phosphorus fluxes of 0.2 to 6.0 mg/m 2 -day under aerobic freshwater conditions which increase to 26 to 34 mg/m 2 -day under anaerobic conditions (excluding the Fillos and Swanson study).
Fillos and Swanson (1975) also reported that dissolved oxygen concentrations had no effect on the flux of ammonia from the sediment. Thomann and Mueller (1987) reported ammonia-nitrogen releases from the sediment in the range of 22 to 44 mg/m 2 -day.
The material composition of the benthic sediment layer can also make a significant difference in the amount of nutrients available for release. Barko and Smart (1986) found that total nutrient concentrations in the sediment were inversely proportional to the sand content for all nutrients. However, the interstitial water concentrations were essentially unrelated to sand content. Similarly, in an analysis of suspended sediment, Dong, et al. (1983) found that, on average, the 77 percent clay-sized particles in the sample contained 88 percent of the total phosphorus.
Another source or sink of nutrients is created by the growth of phytoplankton. Nutrients are depleted from the system during the photosynthesis process. The basic equation describing the photosynthesis process can be written as follows:
P0 4 + NH 3 + C0 2 + H 2 O sunlight > green plants
The photosynthetic process requires phosphorus in the dissolved orthophosphate form. However, as previously discussed, all forms of phosphorus can undergo reactions thus creating new forms. The general reactions are:
Organic P bacterial decomposition > P0 4 Polyphosphates hydrolysis in water > P0 4 P0 4 + multivalent excess coagulant > insoluble metal ions precipitates
Phytoplankton is the non-attached portion of the aquatic plant system which moves passively downstream with the water column (Thomann and Mueller, 1987). Their growth can be represented using the following mass balance equation:
V a (Ph) = (G Ph -Dp h ) Ph V-A v Ph Ph-Q in Ph in -Q out Ph out d u
(2.1)
where
Ph = phytoplankton chlorophyll (mg/L) G ph = growth rate of phytoplankton (1/day) D Ph = death rate of phytoplankton (1/day) Q = flow rate (into or out of segment) (L/day) A = area through which settling occurs (m 2 ) V Ph = phytoplankton settling rate (m/day).
Phytoplankton growth rate, G Ph , is a function of temperature, light and nutrient concentrations and can be expressed as:
= G(T)*G(L)*G(N)
(2.2)
where
G(T) = G„, (1.066) 1 - 20
T = temperature G nax = maximum growth rate (typically ranges from 1.5 to 2.5 with an average value = 1.8/day)
G(L) = J J F (I(i.z)) di dz t z
= .2.718 f r .«! _ -oi o "I K e H I J
(2.3)
where
« j e'** H 1 5
la « o = s
f = photo period (fraction of day light occurs) H = depth (m) I e = saturation light intensity (ranges from 100- 400 Ly/day with average value = 300 Ly/day) I a = total daily solar radiation/photo period.
Di Toro (1978) reported that the light extinction coefficient, K e , could be calculated using:
K e = 0.052 N + 0.174 D + 0.031 P
(2.4)
where
N = nonvolatile suspended solids (mg/L) D = detritus (dead organic material) (mg/L) P = chlorophyll a (/xg/L).
The nutrient growth rate function, G(N), is represented by a Michaelis-Menton type relationship while the death rate function, D Ph , is given by:
Dph Dpi (t) + D z
where
D pl = endogenous respiration rate (1/day) D z = death due to zooplankton grazing (1/day).
As previously stated, phytoplankton growth can be thought of as a potential source or sink of nutrients. Usually, because the phytoplankton are ’’flushed" from the system, a net reduction of water column nutrients is observed. However, in relatively quiescent waters, the phytoplankton might settle out of the water column and thus become a source of nutrients to the sediment layer.
The equations for macrophyte growth are presented in Chapter 3 of this dissertation and will not be repeated here. Like the growth of phytoplankton, macrophyte growth typically results in a net reduction of nutrients. However, for submerged rooted macrophytes, this reduction usually takes place in the sediment layer rather than the water column.
SRAV's are not ’’flushed” from the system unless the stream velocity can create enough drag force on the plant to scour them from the bottom or to break-off stems and leaves. Armstrong et al. (1987) reported SRAV densities in the lower Colorado River were visibly reduced in parts of some cross-sections when the average daily flow exceeded 43 cms (1500 cfs). Moreover, the report stated that average daily flows of 142 cms (5000 cfs) were capable of scouring all the vegetation from some river segments.
Koenig (1987), by examining aerial photographs, concluded that flows greater than 566 cms (20,000 cfs) on the lower Colorado River substantially reduced SRAV densities in most parts of the river to near zero.
In perhaps the most comprehensive work done to date concerning SRAV scour, Vecchio (1989) found that the critical velocity, V c , for Heteranthera dubia and Myriophyllum brasilience (the two most abundant species of submerged rooted aguatic vegetation on the lower Colorado River) could be calculated using the following eguation:
r 2 f c 1 °- 5 Vn = - c I PC d AJ
(2.5)
where
F c = critical drag force (N) C d = drag coefficient A = projected cross-sectional area (m 2 ) p = fluid density (kg/m 3 ) .
Using this equation the scour velocities shown in Table 2.2.4 were generated. The Vecchio study concluded that, based on these velocities, flows in excess of 1416 cms (50,000 cfs) would be needed to scour plants during senescence conditions (during cold weather). During the active growing season, however, flows between 71 and
1416 cms (2500-50,000 cfs) would be sufficient to cause scouring.
Based on this work, Vecchio (1989) developed the following equation for SRAV removal due to flow velocities greater than those indicated in Table 2.2.4:
(2.6)
= K(1-P)SRAV d t where
SRAV = vegetation concentration (gm/m 2 ) K = SRAV scour coefficient (1/hr) P = SRAV protection coefficient.
Furthermore, Vecchio found that the SRAV protection coefficient, P, could be represented by a Michaelis- Menten type relationship in the form of:
p _ SRAV K* * SRAV
where
K* = constant found experimentally by Vecchio to be 0.065 gm/m 2 .
It should also be noted that Vecchio's study reported that SRAV scour typically occurs as a result of breakage rather than root scour. However, it must also be noted that Vecchio's study did not include undisturbed SRAV samples. Therefore, the true effect of root scour may have been difficult to quantify.
2.3 NUTRIENT TRANSPORT
2.3.1 Water Column Advection
Once nutrients enter the stream channel they are transported downstream by advection of the water column subject to the internal reactions previously described. Consequently, the movement of water is an important parameter in determining downstream nutrient concentrations. The techniques used to represent stream flow depend upon flow classification which is a function of factors such as stream geometry, slope, and channel roughness. This section discusses the flow classifications, the different ways of mathematically representing open-channel flow, and the theoretical assumptions and limitations involved with each method.
There are several ways to classify the flow of water in an open-channel. The type ultimately used in the flow modeling process depends upon the assumptions deemed appropriate for the physical situation being modelled. Basically, flow classes can be described as follows (Chow, 1959):
A. Steady flow 1. Uniform flow 2. Varied flow a. Gradually varied flow b. Rapidly varied flow
B. Unsteady flow 1. Unsteady uniform flow (rare) 2. Unsteady varied flow a. Gradually varied unsteady flow b. Rapidly varied unsteady flow
Steady flow occurs when the velocity at every point in the flow is constant with respect to the time increment being considered (Chow et al., 1988). The flow is unsteady if the velocity varies with time. A uniform flow is one in which the mean velocity and depth remain constant, in magnitude and direction, throughout the entire fluid profile (no spatial variation in the direction of flow (Henderson, 1966)). Conversely, a nonuniform (varied) flow is one in which the velocity may vary from point to point (Henderson, 1966). Rapidly varied flow occurs when the depth changes abruptly over relatively short distances and typically refers to phenomenon such as flow beneath a sluice gate or a hydraulic jump or flow over a weir. Gradually varied flow is a flow in which the depth varies slowly along the length of the channel (Chow, 1959).
Despite the fact that steady uniform flow does not truly exist in natural river systems, the uniform flow assumption is used in many engineering applications. The two most well known solutions for solving steady uniform flow are Chezy's formula and Manning's equation. The Chezy formula can be written as:
v = C (RS)°- 5
(2.7)
where v is the flow velocity (L/T), R is the hydraulic radius (L), S is the slope of the energy line, and C is a flow resistance factor called Chezy's C.
Manning's equation can be written as:
v = (1.49/n) R 2/3 S l/2
(2.8)
where v is the flow velocity (ft/s), R is the hydraulic radius (ft), S is the slope of the energy line, and n is the coefficient of roughness (Manning's n). Note that the 1.49 constant in Equation 2.8 requires the use of English units rather than metric.
In some situations steady flow assumptions do not sufficiently describe the hydraulics of the problem. Processes that are sensitive to peak flow and duration such as macrophyte scour need a more detailed description of the flow environment. In these cases, solution of the basic partial differential equations for unsteady flow is required. These equations begin with the equations of continuity and momentum (motion). The level of sophistication desired dictates which terms can be neglected. For example, if the spatial variation in the velocity distribution with respect to depth and across the channel can be ignored, the flow process can be approximated as a one dimensional process. Such simplifications are the basis of the Saint-Venant equations (Chow, Maidment, and Mays, 1988). Derivations of these equations are presented in numerous texts and articles and consequently will not be repeated here. For further information on the subject the reader is directed to the following references (Chow, 1959; Henderson, 1966; Lai, 1986; Chow, et al. 1988; etc.).
Locations where bank and channel erosion substantially alters the channel geometry over time may require an alternative solution. Unsteady flows in open channels with moveable beds can be represented by three partial differential equations (Lai, 1988). The equations consist of (1) the continuity equation for sediment-laden flow, (2) the continuity equation for sediment, and (3) the equation of motion for sedimentladen flow. The simultaneous solution of these equations results in a coupled flow model.
Macrophyte beds can have an impact on the velocity profiles in a stream and should be taken into account when computing stream depths in densely vegetated areas. Figure 2.3.1 illustrates the affects of SRAV's on water column velocity and shear stress. The randomness of SRAV growth patterns in the lower Colorado River makes this computation unnecessary.
2.3.2 Role of Sediment Transport
Nutrients which adsorb to sediment may have different transport characteristics due to momentum transfer losses that occur between the fluid and the solid particles. Particulate phosphorus attached to sediment particles may settle out in slow moving reaches of the river or along the banks thereby increasing the bed concentrations and reducing the water column concentration of phosphorus. Yet, in other segments of the river, nutrient reservoirs may be scoured out resulting in increased water column concentrations.
Since one of the primary objectives of this dissertation is to develop a model capable of predicting nutrient balances, the fate and transport of nutrients would be incomplete if sediment transport processes were omitted. In addition, suspended sediment concentrations play an important role in determining the light extinction coefficient. In fact, Lung and Testerman (1989) found that turbid waters could reduce algal growth rates in the James Estuary by over eighty percent. Consequently, to help explain the sediment erosion and transport processes, the next subsections of this chapter identify the typical sources of sediment, present the different methods of estimating sediment yields, and discuss several of the more popular equations used to calculate sediment loads.
2.3.2.1 Sources of Sediment
Sediment loadings to streams and rivers originate from seven principal sources (Brown, 1949). These can be summarized as follows:
a. Sheet erosion b. Gullying c. Stream-channel erosion d. Mass movement of soil via landslides, slumps e. Flood erosion f. Erosion incident to cultural development g. Mining, industrial, and sewage waste discharges.
Brown concludes that sheet erosion is the dominant source of sediment in larger drainage areas which receive more than twenty inches of rainfall annually and in areas where land is used primarily for agricultural purposes. In forested and range areas receiving less than twenty inches of precipitation, gullying and
stream-channel erosion are the dominant sources.
Novotny and Chesters (1981) cite a more quantitative list of sources as presented in an article written by Brant et al. (1972) for Dow Chemical Company. Their list of sources includes:
a. Natural erosion 40 Tonnes/km 2 /yr b. Agricultural erosion 100-4000 Tonnes/km 2 /yr c. Urban erosion up to 50,000 Tonnes/km 2 /yr d. Highway erosion up to 50,000 Tonnes/km 2 /yr e. Stream bank erosion highly variable.
The authors chose not to estimate a value for stream bank erosion due to the highly variable nature of channel scour. In general, they state that stream bank erosion is a source of large grain size sediments and accounts for a smaller percentage of total erosion than the remaining sources. However, they are quick to point out that for some ephemeral streams in arid regions, the value can be quite high.
The conclusions drawn by Brant et al. (1972) regarding the importance and nature of stream bank erosion are not widely accepted. A U.S. Army Corps of Engineers (1978) study indicated that appreciable stream bank erosion occurs in approximately 25 percent of the streams in the United States (see citation in Dharamdial and Khanbilvardi, 1988). Other researchers have found that bank erosion contributes significantly to the suspended load (Odgaard, 1987; Great Lakes Commission, 1975; Harris, 1989). In fact, based on a mass balance approach, Harris (1989) found that approximately 90
percent of the sediment transported in the Colorado River was the result of stream bed/bank erosion. A similar study on the Sacramento River by the U.S. Army Corps of Engineers (1983) concluded that 6.2 of the 10.5 million tonnes (nearly 60 percent) of total sediment inflow was derived from bank erosion (1 metric ton = 1.1 US tons). Moreover, studies on two lowa rivers showed that silt and clay fractions made up about 75 percent of the material eroded from the stream bank (Odgaard, 1984). For the lower East Nishnabotna River in lowa this translates into approximately 245,000 and 73,000 tonnes/year of silt and clay, respectively (Odgaard, 1987). Odgaard estimated that this comprised about 30- 40 percent of the suspended load. Methods for estimating stream bank erosion are presented in the next section.
Mining of sand and gravel also contributes to the sediment load in the Colorado River. Although the Texas Water Commission (TWC) regulations require sand and gravel operations to construct ponds, ditches, berms or other structures to prevent wash water from entering the river and to file permits if the wash water is going to enter the river, the regulations are hard to enforce primarily because of jurisdictional debates. Mining operations on land are the responsibility of the TWC while in river operations are under the control of the Texas Parks and Wildlife Department. The difficulty occurs because so many operations are located near or on the banks of the river. In a 1988 telephone conversation with Rollin Macßae of the Texas Parks and Wildlife Department, Harris (1989) discovered that if
the berm surrounding the mining operation contains only one opening, the pit area is not considered as state waters. Moreover, Harris (1989) also discovered from Macßae that aerial photographs show steady plumes of muddy water entering the river from old mining sites even though no permits were filed with the TWC as of August 3, 1988.
2.3.2.2 Methods of Estimating Yield
Erosion is controlled by climate, soil properties, topography, vegetative cover, and human activities. Since these factors combine to cause erosion when surface runoff occurs, most creditable equations and methods of estimating sediment yield include them as variables. The most notable of these equations is the Universal Soil Loss Equation (USLE), although more recently a Modified Universal Soil Loss Equation (MUSLE) has become popular.
The original USLE was developed in 1965 by Wischmeier and Smith as a result of more than 40 years of field observations collected by the USDA (Novotny and Chesters, 1981). The MUSLE altered the original equation by replacing the rainfall energy factor with a runoff factor thereby eliminating the need for delivery ratios and improving the accuracy of sediment yield predictions (Williams and Berndt, 1977).
In its final form, the Modified Universal Soil Loss Equation can be written as (Williams, Nicks, and Arnold, 1985):
Y = 11.8 ( V « q ) 056 (K)(C)(PE)(LS)
(2.9)
where
Y = sediment yield from subbasin (tons/storm) V = surface runoff volume (m 3 ) K = soil-erodibility factor C = crop management factor PE = erosion control practice factor LS = slope length and steepness factor q p = peak flow rate for subbasin (m 3 /s).
The obvious criticism of Equation 2.9 is the use of all the factors. Despite the fact that most of the associated tables, charts, and equations used to estimate these values are derived from field data, short term sediment yield predictions can be off by a substantial amount (Ambrose, et al., 1988).
Harris (1990) computed suspended sediment loadings from the tributaries of the lower Colorado River using a computer program called HILOADSS. HILOADSS uses the Soil Conservation Service (SCS) method of estimating storm water runoff and concentration values to determine sediment loading. HILOADSS is a descendent of the HILOADS program developed by Miertschin (1986). In her analysis Harris found that only about seven percent of the total suspended sediment in the lower Colorado River was derived from the tributary areas.
More recently there have been attempts to develop empirical models. These models rely on model calibration against field data to verify the coefficients rather than on the extensive field data collection of the USDA. The Negev sediment generation model and the David & Beer model are two of the more well known models (Novotny and Chesters, 1981). The Negev model is based on the assumption that raindrop impact and overland flow account for the generation and transport. Novotny and Chesters (1981) report that this model appears to be most suited for analysis of fine grain soil particles such as silt and clay.
The David & Beer model assumes soil erosion can be classified into rill and interrill classes. Rainfall splash produces erosion in the interrill class whereas overland flow contributes to the soil movement along well defined rills or paths.
Harris (1990) identified the major sources of sediments in the lower Colorado River as those derived from stream bank erosion and bed (channel) scour. Stream bank erosion is influenced by a number of parameters such as shape geometry, flow velocity, tractive force (i.e., drag force exerted on banks by water), momentum of objects carried by water, waves, and climatic factors (Simons and Li, 1982). Bed scour is a function of the depth of flow, slope of energy gradient, density of sediment-water mixture, size of bed material, particle fall velocity, gradation of bed material, shape factor for the reach and cross section, subsurface flow through bank material, and piping of river bank (Simons and Li, 1982) .
Few sediment transport models are truly capable of handling bank scour. Dawdy and Vanoni (1986) claim that to properly model bank erosion as a dynamic process
requires a three-dimensional flow system model. Furthermore, the processes of bank erosion and erosion in bends would require different model formulations. Nevertheless, a few researchers have developed expressions for bank and channel erosion. Chang (1982) presented a mathematical model based on total stream power. Odgaard (1987) reported that several studies have attempted to relate bank erosion to either the channel width, the drainage area, or the ratio of the channel width to radius of curvature. Odgaard used the width/radius concept to develop an expression for the average rate of erosion in river bends. The equation can be written as:
r t 1 — = E M — 1 ♦ -Tr— F u r c L 2r Cj
(2.10)
in which
r _ e 3a n+l c « E " 8 2 k n+2 Fl)c ~ s.
Ur Fdo = 7- 7 — <p .?. -P* g D I p J
c 1 n r c *[, g ' F = 1 - exp -B —r 1- — d l t
R - K 2 b (m +1 ) 2 G c
where
v' = erosion rate (m/yr)
u = average mean velocity (m/yr) b = width of channel (m) e = erosion constant k = von Karman constant m = friction parameter D = particle diameter (m) g = gravity (m/s 2 ) d c = depth at center line (m) r c = radius of curvature at center line (m) £ = crossover angle © = Shields' parameter (=0.06) 4. = bend angle p = fluid density (kg/m 3 ) P s = sediment density (kg/m 3 ) o: = ratio of projected surface area to volume □□ = average center line velocity (m/s).
The volumetric rate of erosion can then be obtained by multiplying the erosion rate, v', by the areal factor and by the height of the cutbank.
In addition to bank erosion, sediment can also be scoured from the bottom of the river bed. The physical process of bed scour leads to a phenomena called armoring. As the stream bed degrades, smaller sediment particles are eroded and transported downstream more rapidly than larger grain sizes which results in the stream bed becoming coarser. The larger particles create a protective surface (armor) which inhibits further erosion (Dawdy and Vanoni, 1986). Although the bed armoring process is not well understood, several
attempts have been made to quantify the parameter. Karim and Kennedy (1982,1983) use a method involving several empirical coefficients and the fraction of the bed surface which is covered by immobile bed material.
2.3.3 Sediment Transport Theory
Up to this point the term sediment transport has been used to define the total sediment load in the channel. However, now it is necessary to distinguish between suspended sediment and bedload sediment. By definition, suspended sediment (sometimes incorrectly referred to as ’’washload” since the two are not identical) is the sediment that is always transported in suspension by the water. On the other hand, bedload is the material which is transported downstream predominantly in contact with the bottom of the stream. A problem occurs when one tries to determine the depth where the bedload layer ends and the suspended sediment layer begins. Presumably a fine line exists between the suspended sediment layer and the bedload, although, practically speaking, defining such a precise boundary becomes a very difficult task. Figure 2.3.2 illustrates a typical logarithmic suspended sediment concentration profile and the bedload layer concentration.
Since suspended sediment is being carried by the water, typical particle sizes are smaller than those found in the bedload material. Suspended sediment is usually made up of particles ranging from very fine clay to a coarse sand. Bedload is made up of particles
ranging from medium sand up to coarse gravel. A guideline for relating particle size to sediment class is given in Table 2.3.1. Because most adsorbed pollutants have been found to be transported in conjunction with suspended sediment rather than with bedload (Novotny and Chesters, 1981), this review will focus primarily on the role of suspended sediment.
The terminology found in the literature can be quite confusing because some authors subdivide suspended sediment into wash load and suspended sediment categories (Johnson, 1970). Using these definitions, wash load represents fine silt particles which do not usually settle out and suspended sediment consists of fine sand grains which can settle out. The primary reason for this appears to be the difference in the way runoff affects the timing of wash load and suspended sediment loadings. Wash load is the portion of the suspended sediment that is made up of grain sizes finer than those native to the stream bed as determined by the available upstream supply (Krenkel and Novotny, 1980). The smallest amount of overland flow can carry the entire wash load into the receiving stream. Therefore, the peak wash load concentrations usually occur before the peak of the storm runoff hydrographs.
It is probably an understatement to say that sediment transport is not an exact science. Hundreds of theoretical and empirical sediment transport equations can be found in the literature. Furthermore, the various assumptions made in one derivation can be slightly altered to form different sediment
concentration equations. Consequently, as stated previously in Chapters 1 and 2, this literature review does not attempt to present all of the different variations. Instead, an overview of the most popular and unique formulations that are applicable to the research goals is given.
2.3.3.1 Suspended Sediment Equations
Because the transport characteristics of suspended sediment and bedload are seemingly unique, the basic theoretical derivations and the resulting loading equations can initially be examined independently. Most theoretical equations for suspended sediment are based on the assumptions that solid particles are kept in suspension by turbulence and the effects of turbulence on the particles are analogous to the diffusiondispersion (also called diffusion-advection) process (Graf, 1971), (Nicholson, 1979). Armed with these general assumptions, the diffusion equation for incompressible, laminar flow can be written as:
= -U-V C + € M V 2 c at
(2.11 a
or, in tensor form as:
3C 11 30 X -|J . + €tn St l dxj 3X13X1
(2.11 b
where
i = x, y and z directions C = concentration (mg/L) U = velocity (m/s) t = time (s)
x = directional distance (m) = diffusion coefficient (m 2 /s).
For turbulent flow conditions, the concentration and velocity components used in Equation 2.11 b become:
C= C a +C' and U A = U ai + U' ai
where, C a and U ai represent the mean values of concentration and velocity at a given position, respectively, and C' and U' ai are their fluctuations.
Substituting these values into Equation 2.11 b and taking into account the time averaging properties associated with two of the derivative terms gives:
dC-a . । dCa _ 11- aC st a 1 ax i a 1 ax j ax iaxi
(2.12)
By defining a coefficient of turbulent diffusion as:
8 C = e t —— = -c- U' ai d X J
Equation 2.12 can be written as:
__ । । aC a 3 _ _ 3t a 1 3Xj SXj e-t dx j erq
(2.13)
Furthermore, assuming that molecular diffusion and turbulent diffusion are independent, we can add the two values together to form a diffusion tensor.
€ i, j +
Combining this with Equation 2.13 and dropping the time average distinctions leaves,
3C __ j. 0 + 3 (e ■ _ BC Bi 1 BXi + BXi Bxj
(2.14)
If diffusion is considered only along the principal axes, the three-dimensional form Equation 2.14 can be expanded into is:
= -u aC_ _ b aL _ u~— at ax ay z az
8 r 8 C s 3 ( 8 C 8 r 8 C 8x x 3x 8y 8z 8z
(2.15)
For the one-dimensional case where z is the vertical axis representing depth, Equation 2.15 becomes:
aC _ 11 aC 3 , aC x 17’ it’
(2-16)
A number of assumptions are necessary to relate Equation 2.16 to suspended sediment concentration. First, it must be assumed that steady state conditions exist (time derivative = 0). Second, the kinematic eddy diffusivity must be independent of depth. And third, the settling of sediment particles is assumed to be analogous to the flow velocity in the vertical direction of the advection-diffusion equation (a very large assumption in this author's opinion). With these assumptions in mind, Equation 2.16 can be simplified to:
Cw * e = 0 dz
(2.17)
where w = -U z = settling velocity of grains in fluid.
If additional constraints are placed on the kinematic eddy diffusivity by assuming its value is a constant, the solution to this eguation becomes:
C _ -w< z-a > /€ C a
(2.18)
where
Ca = initial Concentration at depth a (mg/L) a = depth at which C a was measured (m) z = depth (m) w = settling velocity (m/s) €= n = kinematic eddy diffusivity (m 2 /s). p
In natural streams, the turbulent intensity, and consequently the diffusivity (kinematic eddy diffusivity), will vary with the depth of flow. This fact led to the classical work of Prandtl and von Karman (Henderson, 1966). Prandtl and von Karman solved Equation 2.17 using a variable kinematic eddy viscosity and Prandtl's theory of ’’mixing length”. Their solution can be written as:
r T ru 1 C _ aIH - zj C a [z (H - a\
(2.19)
where
k = von Karman's Constant (0.4) V* = shear velocity H = total depth of flow.
Note: k = 0.4 assumes homogeneous fluid which may or may not apply to sediment-enriched flow.
The practical use of both Equations 2.18 and 2.19 is somewhat limited since both equations require that an initial concentration value at some depth be specified.
In 1955, Einstein and others attempted to modify the velocity distribution equation used in the Prandtlvon Karman formula by rationalizing that the sediment in suspension must influence the flow velocity (Graf, 1971). This allowed them to justify dividing the sediment-laden flow into a heavy-fluid zone and a lightfluid zone. The heavy-fluid zone is the layer close to the bed where the sediment concentration is high. Conversely, the light-fluid zone represents the layer where the sediment concentration does not significantly effect the fluid density.
Writing, combining, and simplifying the general diffusion equation presented as Equation 2.14 for both the rate of change of suspended matter and for water produces:
[e s * C (e; - e s )l — * (1 - OCw = 0 L J 3z
(2.20)
where the subscripts s and 1 represent the solid and liquid values of the diffusivity, respectively.
Letting the solid diffusivity equal the liquid diffusivity at equilibrium, Equation 2.20 becomes:
e s + (1 - C) C w = 0 8 Z
(2.21)
For the light-fluid zone this equation reduces to the same equation we developed in Equation 2.17:
3 C r - n e+C w = U az
(2.22)
However, for the heavy-fluid zone, the equation becomes:
r r r ' f 7 f C 1-Ca = 1-z/H B-(1-a/H) 1/2 _ B- (1 - z/H ) 1 /£L _,
(2.23)
where
B = an integration constant in the velocity distribution law <= 1. H = total depth.
In general, the two arbitrary constants in Equation 2.23 produce a slightly better fit for observed data (Graf, 1971). Even so, the added complexity of the equation does not always justify its use.
The previous equations give concentration values at a specified depth. To calculate the suspended load rate per unit of time, the following equation is used:
D g ss = J cu dz a
(2.24)
The equation used for predicting the suspended sediment concentration, C, in Equation 2.24 depends on the method selected as most appropriate. Similarly, the lower limit of integration (a) for some methods should be zero. Also, since each C refers to only one settling velocity, w, the spatial and temporal distribution of
each sediment grain size is required for the accurate determination of total sediment discharge.
Einstein (1964) attempted to quantify the lower limit to Equation 2.24 by claiming that "a” should be equal to a distance of about two grain diameter from the bed. He based his argument on the belief that eddies become too small to support the particles at distances less than two grain diameters from the bottom. The grain diameter used here is twice the D 6 5.
When sufficient field measurements are available, the most widely used method for estimating suspended sediment yield is the flow duration-sediment rating curve (Grenney and Heyse, 1985). The sediment rating curve is developed by plotting a log-log graph of sediment load against discharge. A regression line is then drawn through the points resulting in an equation in the form of:
S = a * Q b
(2.25)
where S is the sediment load (tons/yr); Q is the instantaneous discharge (cfs); and (a) and (b) are the constants representing intercept and slope.
There are two important facts regarding Equation 2.25 which should be discussed. First, as pointed out by Grenney and Heyse (1985), the regression between sediment load and discharge contains a spurious flow correlation which tends to produce a higher coefficient of determination, r 2, than actually exists between
sediment concentration and flowrate. Failure to understand this point could cause misinterpretation of the results. Second, the accuracy of the sediment concentration data used to construct the sediment rating curve should be closely examined. Much of the sediment concentration information currently available is derived from a single sample location taken near the middle of the river at a representative depth. Given the inherent spacial and temporal variability of suspended sediment concentrations, this too could lead to erroneous predictions.
2.3.3.2 Total Sediment Load Equations
Rather than attempting to divide sediment transport phenomena into separate suspended and bedload equations, many attempts to model sediment have combined the processes into one equation. These attempts have had varying degrees of success in part because ’’thus far, it has been impossible to develop any fully rational equation for total sediment load of a stream ...” (Morris and Wiggert, 1972). Still, in many instances, these total load equations can predict the transport of sediment reasonably well.
The most notable attempt at a total sediment equation came from the bedload work done by Einstein (Morris and Wiggert, 1972). He tried to expand upon his original work to predict the suspended sediment concentration at the boundary layer in order to establish an initial concentration value.
Unfortunately, his experiments produced rather poor results when compared to actual measurements (Morris and Wiggert, 1972).
R. A. Bagnold developed an expression for calculating total sediment load, q s , based on energy considerations (Yalin, 1977). In the original form the equation can be written as:
Rse b U s —_— = + e s (1 -e b ) — 'ToUrn lan w
(2.26)
where
U B = average transport velocity of the suspended solids (m/s) w = terminal velocity of the sand grains (m/s) e e = efficiency with respect to the suspended load = available energy.
Many of the parameters used in Equation 2.26 are unknown values which can only be determined experimentally. Bagnold conducted experiments and found that for sand grains with a diameter less than or equal to 0.50 mm, Equation 2.26 could be simplified to the following equation:
= [ 0. 17 ♦ 0.01 —] oUm w
(2.27)
where
U B = average flow velocity (m/s).
Another approach was attempted by Larsen who treated sediment load as a single phenomenon and then tried to develop an empirical relationship based on
model studies and dimensional analysis (Morris and Wiggert, 1972). The Larsen equation resulted in:
a . r v Tf f V 2 ] V Dl z , u 9s C v [Bg H J H <VsH V)
(2.28)
where
H = total depth (m) V = velocity (m/s) f = friction factor D = grain size of sediment (m) x, y, z, and C are found experimentally.
Albertson and Garde developed a similar equation based on the same logic as Larsen. Their equation is written as:
1.36 V 4 n s 9 s ~ 7 v (10) 5 ] D s/2 H
(2.29)
where all values are expressed in standard English units of feet, seconds, and pounds.
Both of these equations have been reported to produce reasonable results under the right set of conditions. However, a great deal of caution should be exercised in trying to apply these equations to any situation. This could be said about all sediment equations.
Ackers and White (1973) have proposed a different method for determining the sediment transport rate. They argue that many methods which rely on total shear stress have to separate the bed shear into two
components. The first component is the shear on the grains while the second is the nontransporting form loss. Since the transport rate is extremely sensitive to the shear on the grains, inaccuracies in the dividing procedure may lead to large errors in the estimate of sediment transport. Consequently, their "new approach" uses dimensional analysis and effective shear stress rather than total shear stress. First, they assume that the coarse sediment is transported primarily as a bed process and that the fine sediment is transported within the body of flow.
Next they introduce a term called the ’’mobility number” which is the ratio of the appropriate shear force on the unit area of the bed to the immersed weight of a layer of grains. They present this as:
-I 1-n F = VS V D(s- 1 ) 5> j oa |
(2.30)
where
d = mean depth (area/width) (m) D = sediment diameter (mm) g = gravity (m/s 2 ) s = mass density of sediment relative to fluid V = mean velocity (m/s) V* = shear velocity (m/s) ex = rough-turbulent coefficient (approximately 10) n = transition exponent = 0 for coarse sediments (sand diameter) = 1 for fine sediments.
For transitional sizes, n varies between 0 and 1. The actual value is determined by a dimensionless grain diameter term that Ackers and White (1973) present as:
, .. ] 1/3 D qr = D
(2.31 a
and,
n = 1 - 0.56 log(D gr )
(2.31 b
where the denominator in Equation 2.31 a is the kinematic viscosity of the fluid.
The useful work done in sediment transport in relationship to the stream power gives an expression for the efficiency of the transport process. Then, by combining the efficiency with the mobility number, one gets:
PF 1 n Ggr = C - 1
(2.32)
where the values for the coefficients for coarse and transitional (1 < D gr < 60) sediments are given by:
C = transport function coefficient C = 0.025 for coarse sediments Log C = 2.86 log(D gr ) - [log(D gr )] 2 - 3.53
A = value of F at initial motion A = 0.17 for coarse sediments A = 0.14 + [ 0.23/(D gr °’ 5 ) ]
in = transport function exponent m = 1.50 for coarse sediments
m = 1.34 + (9.66/D gr ) for transition range.
Notice that Ackers and White (1973) concede that fine sediments are transported as a function of total bed shear and thus, this method becomes inapplicable. Their definition of fine sediment is anything with a grain size diameter of less than 0.04 mm.
Finally, the value obtained for G gr in Equation 2.32 can be substituted into Equation 2.33 to determine the dimensional sediment flux value, X, in parts per million by weight of fluid flux.
r x d f M n b S r - s D V \ J
(2.33)
Arnold, Williams, and Nicks (1988) present a sediment transport procedure in their Simulator for Water Resources in Rural Basins (SWRRB) draft users manual which divides sediment routing into simultaneously occurring deposition and degradation components. Here, deposition is based on fall velocity and degradation is based on Bagnold's stream power concept.
Without rewriting the entire section from their manual, a complete derivation of their sediment transport equation would be impossible. The highlights of their solution and general definitions of variables are given below. For a more extensive review of the procedure, the reader is directed to the Arnold, Williams, and Nicks (1988) reference. Basically their approach is to develop equations for deposition (DEP) and degradation (DEG) and combine them with the sediment
yield derived from the MUSLE. These equations can be written in their final form as:
DEP = YDj (1-DR)
(2.34 a
DEG R = « sp y^- 5 DUwIPRf d w S w V c ) 1 - 5
(2.34 b
where
YD = sediment of particle size j entering channel (mm) DR = sediment delivery ratio PR = peak rate adjustment factor d w = flow depth (m) S w = surface water slope DU = duration of time increment (s) V c = channel velocity (m/s) a sp = parameter varying with maximum stream power.
In addition to river degradation, the bed material can also be scoured. This is taken into account using:
(2.35)
DEG b = K C * DEG R
where K and C are MUSLE factors for the stream channel.
The resulting equation for sediment reaching the basin outlet can be written as:
(2.36)
SED t = Y - DEP + DEG e + DEG b
There have been a number of computer models developed to solve the sediment transport equations presented above. In 1987 the Interagency Sedimentation Work Group found that some 48 sedimentation models had been developed in the United States alone (Fan, 1988). Most of these models require a great deal of professional judgement and are strongly data-dependent. Furthermore, almost all the models represent suspended and bedload as a lumped total load (Holly, 1988). Since this project is primarily concerned with the effects of suspended sediment, these types of models would be inappropriate.
A final note concerning the reliability of sediment transport models and equations should be made. In an assessment of the present state-of-the-art modelling techniques, Dawdy and Vanoni (1986) stated that ’’the choice of a model at this time is arbitrary...” and that the validity of the results depends more on the modeler than on the model. Their paper also presented a summary of predicted versus observed sediment concentrations performed by White, et al. in 1973. As indicated in Table 2.3.2, predictions made by using the ’’best" method only coincide with observed values 64 percent of the time. Furthermore, the acceptance range is 1/2 to 2 times the observed concentration. For an observed concentration of 100 mg/1 the acceptance range would be from 50 mg/1 to 200 mg/1.
In order to use the total load equations to predict suspended sediment load, a number of assumptions and hypotheses would have to be made. The amount of uncertainty demonstrated by the total load equations and models combined with the numerous suppositions would undoubtedly lead to criticism of these methods. For
these reasons, it is proposed that the widely used rating curve technique be used to estimate suspended sediment concentrations in the lower Colorado River.
2.4 SUMMARY OF EXISTING WATER QUALITY MODELS
The proliferation of water quality models has resulted in literally hundreds of different models. These models range in complexity from simple single element DO-BOD models to sophisticated three-dimensional programs containing advanced kinetic equations for multi-parameter analyses. Some are widely available through federal agencies such as the USEPA and others are proprietary in nature. The only common denominator for the models is that they are each based on the conservation of mass principle.
The lack of a uniform grading system often makes model comparison difficult. McCutcheon (1989) suggests that models should be classified into one of the following four complexity levels:
Level 1: Manual or graphical methods based on deterministic equations, used for crude screening over extensive areas.
Level 2: Simple computerized models based on equations which only approximate the basic processes involved, used for fine screening or crude planning.
Level 3: Computerized models of intermediate complexity incorporating some approximation of the basic processes, used as a fine planning model.
Level 4: Advanced mechanistic computerized model
requiring intensive data collection, used for detailed design and management of smaller areas.
However, while McCutcheon suggests four complexity levels, the WASP 4 manual identifies six levels of complexity for nutrient modeling (EUTRO4) and four complexity levels for toxicity modeling. Thomann and Mueller (1987), on the other hand, describe three levels of analysis for toxic studies.
Another way to classify models is by their division of computational elements. This method is much less subjective but reveals no information regarding the kinetic equations contained within the model. The classes have been defined by various authors (USEPA, 1985; McCutcheon, 1989) as:
Zero-dimensional: single element, completely mixed reactor, rarely used.
One-dimensional: lateral and vertical variations considered unimportant, series of elements in downstream direction, most common approach for stream and river analyses.
Two-dimensional: model describes variations in two directions, vertical direction usually neglected in stream and river applications.
Three-dimensional: all three gradient directions considered, used to describe complex mixing zones.
Table 2.4.1 shows a feature comparison of several of the most popular nutrient models. A complete review of each model is outside the scope of this dissertation. The time required to learn the features and limitations of every model would be enormous. An excellent review of computer models available through about 1982 is provided by McCutcheon (1989) for those who would like a more thorough discussion on the subject. However, it should be pointed out that most models undergo continuous revisions. As a result, many of the latest features would be excluded from McCutcheon's discussion. Table 2.4.2 contains additional information on the models' spatial domain, time domain, and state variable systems.
In selecting an appropriate model, one must consider both the intended application and the amount of data available for calibration. The extra time, data requirements, and costs generally associated with the more complex models may not be justifiable if there is insufficient field data and rough estimates are needed to fill in the gaps. The WASP 4 model (and subsequently the EUTRO4 model) was selected based on the widespread availability of the model, USEPA acceptance of the model, the dynamic and sediment layer capabilities, and the compartmentalized approach which made program modification more straightforward.
2.5 PREVIOUS WATER QUALITY MODELING STUDIES OF THE COLORADO RIVER
2.5.1 Study Area Under Consideration
As shown in Figure 2.5.1, the Colorado River
drainage basin runs through central Texas in a southeasterly direction from the Texas-New Mexico border to the Gulf of Mexico and has a total contributing drainage area of approximately 30,600 square miles (USGS, 1985). Although much of the flow is generated from drainages upstream of the study area, the study area itself covers only a small portion of the total drainage basin beginning immediately below Longhorn Dam in Austin, Texas and ending at the Town of Bastrop, Texas. This 58.4 mile (94 km) reach encompasses only about 970 square miles of drainage area.
Figure 2.5.2 illustrates a generalized representation of the Central Texas study region. As shown, the study area is located downstream from the City of Austin and the Highland Lakes chain of reservoirs. The locations of major tributary inflows and a more detailed map of the Lower Colorado River drainage basin are exhibited in Figure 2.5.3. While many of the tributaries represent large land areas, historic flow records and field data indicate that only two major sources of inflow exist during even short dry periods, those being Walnut Creek and Onion Creek. Furthermore, USGS gaging records indicate that during long dry periods, even the flows in these creeks may become insignificant (except for the waste water returns into Walnut Creek).
2.5.2 Norman Johns' Pre-WASP3 Model
In August of 1988, using a previous version of the USEPA's Water Quality Analysis Simulation Program (a prereleased version of WASP 3), Norman Johns completed a one-dimensional, steady state water quality model of the lower Colorado River from Longhorn Dam in Austin to Bastrop under the direction of Dr. Neal Armstrong of the University of Texas at Austin (Johns, 1988). The WASP 3 model was based on the conservation of mass principle and used finite different techniques to approximate the following mass balance equation:
aC, 1 a 2 E X AC<) i a(U x AC,) __J_ = 2 J 2 2—J_ ♦ s -U ( at A j X £ A ax j j
(2.37)
where
Cj = concentration of constituent j (mg/L) t = time (s) A = cross-sectional area (m 2 ) E x = longitudinal hydrodynamic dispersion coefficient (m 2 /s) U x = average water velocity (m/s) Sj = other sources, losses, or transformations of constituent j (mg/m 3 -s) Wj = point and nonpoint source loads of constituent j (mg/m 3 -s).
Two types of differencing schemes are used in the program. WASP 3 uses a backward-difference scheme in space and a forward-difference solution procedure in time to approximate the solution.
Using the methodology described above, Johns (1988) first calibrated the modified WASP 3 model and then
evaluated the consequences of three future scenarios. Figure 2.5.4 graphically illustrates the results of the calibration run. Johns reported that values for virtually all of the kinetics parameters were obtained from the literature. The values of the maximum plant growth rate and the light extinction factor for plants were adjusted until ’’good agreement was reached.”
Three different scenarios were examined in order to see if nutrient removal from the City of Austin's wastewater discharges would have a significant influence on potential macrophyte growth. Scenario 1 used projected 1990 discharges with no additional nutrient removal. Scenario 2 used 1990 discharges with a reduction in the total phosphorus concentration to 1.0 mg/L. Scenario 3 was the same as Scenario 2 except that nitrogen from the wastewater treatment plants was also removed. The consequences of these actions can be seen in Figure 2.5.5. Based on these results, the study concluded that nutrient removal at Austin's wastewater treatment facilities would not significantly reduce the potential densities of submerged rooted aquatic vegetation (macrophytes) and that light availability was the limiting factor governing SRAV growth.
2.5.3 Elizabeth Krousel's QUAL2E Model
Concurrent with the development of this work, Krousel (1991) modified a version of QUAL2E to investigate long-term quasi—steady state macrophyte and epiphytic algae growth rates using the Freundlich
isotherm to relate water column concentrations to sediment concentrations. Even though the boundary conditions and flows are considered constants in QUAL2E, Krousel was able to use the temporal variations in temperature and other factors to reasonably predict plant densities from April to July, 1986.
Krousel's conclusion was that light was the limiting growth factor effecting macrophyte or epiphytic algae growths rather than phosphorus or nitrogen. In fact, as illustrated in Figure 2.5.6, Krousel's model predicted virtually no change in SRAV concentrations and only minor changes in epiphytic algae densities as a result of various waste water loading scenarios.
k = height of vegetation = undeflected height of vegetation k n - n = minimum deflected height of vegetation U = vertical velocity profile Y’ = effective channel bottom T = total shear = soil boundary shear T v = shear transferred to vegetation T e = effective shear *T O = bottom shear if no vegetation were present
L - Average stream velocity H- Average stream depth Op- Sediment discharged in bed-load layer B-B 1 -Hy pot he t ica I plane separating bed-load layer from, suspend e c - I c a d layer x-Distance in the longitudinal direction y- Depth from the water surface C-Suspended-sediment concentration 0 ।-Sediment concentration in bed-load layer Op-Sediment concentration at bottom of suspend e d - I o a d layer a'-D e p t h of bed-load layer
a = Hicher removal are attained when coagulating chemicals are usee b = COD lnflu .„t - [BOD u (removed)/0.9) c = Nmfxu.nx - 0-12 (excess biological sludge), lb; v - 0 026 (excess biological sludge), lb. r influent w. w , . . . • _ . d = Depends on resin used, molecular state, anc efficiency cesirec
All values in kg/ha-yr
A S .=n,to to. to model mW. Ac me numoer of user-sped Bed cousmuenis that may ...dude BOD. DO. and
Figure 2.2.1 Phosphorus Cycle in Rivers (After Baca and Arnett, 1976 in USEPA, 1985)
Figure 2.2.2 Nitrogen Cycle in Rivers (After Baca and Arnett, 1976 in USEPA, 1985)
Figure 2.2.3 Land Use versus Phosphorus Loading (After Novotny and Chesters, 1981)
Figure 2.2.4 Land Use versus Nitrogen Loading (After Novotny and Chesters, 1981)
Figure 2.2.5 Potential Nutrient Concentration Gradients in a Benthic Sediment Layer
Figure 2.3.1 Effect of SRAVs on Stream Velocity (After Thompson and Roberson, 1976)
Figure 2.3.2 Typical Suspended Sediment Concentration Profile
Figure 2.5.1 Location of Colorado River Basin within Texas (After Leifeste and Lansford, 1968)
Figure 2.5.2 Generalized Map of Colorado River Study Area (After Harris, 1989)
Figure 2.5.3 Colorado River Between Austin and Bastrop (After Harris, 1989)
Figure 2.5.4 Calibration Run for Pre-WASP3 Model (After Johns, 1988)
Figure 2.5.5 Predicted Potential Macrophyte Growth Using Pre-WASP3 Model (After Johns, 1988)
Fiaure 2.5.6 Predicting Potential Macrophyte and Epiphytic Algae Densities Using QUAL2E Model (After Krousel, 1991)
Process BOD COD ss N P TDS Sedimentation, % removal 10-30 - 50-90 - - - Flotation* , % removal 10-50 - 70-95 - - - Activated sludge, mg/1 <2 5 b <20 c c - Aerated lagoons, mg/1 <50 - >50 — - - Anaerobic ponds, mg/1 >100 — <100 — — — Carbon adsorption, mg/1 <2 <10 <1 - - - Ammonia stripping, % removal Denitrification and >95 nitrification, mg/1 <10 — — <5 — Chemical precipitation, mg/1 — — <10 — <1 d Ion exchange, mg/1 d d
Table 2.2.1 Nutrient Concentrations in Domestic Waste Water (After Krenkel and Novotny, 1980)
Type Total Phosphorus Total Nitrogen Approx. mean Approx. range Approx. mean Approx. range Forest, natural 0.4 (0.01-0.9) 3.0 (1.3-10.2 ) Atmospheric ra infall 0.2 (0.08-1.0) 8.0 dry fallout 0.8 16.0 Urban 1.0 (0.1-10) 5.0 (1-20) Agricultural 0.5 (0.1-5) 5.0 (0.5-50)
Table 2.2.2 Nonpoint Source Loadings (After Thomann and Mueller, 1987)
Nitrogen Regional mooets Nationwide mooe. East Central Vves) *n agncuitura! * uroar. Average A' 67‘„ limns 67‘„ Average A limns Average A' 675« limns Average A' limns 0 0.57 25 0.90 50 1.45 75 2.3) JOO 3.69 (0.35-0 92) (0.56-1.46) (0.89-2.34) (1 43-3.74) (2.28-5.97) 0 51 (0.33-0 771 0.87 (0.57-1.32) ) 49 (0.98-2.26) 2.55 (1.68-3.88) 4.37 (2.87-6.64 0.57 0.8? 1.22 1.78 2.61 (0.37-0 881 (0.54-1.28) (0 79-1.88) (1.15-2.74) (1.69-4.02) 0.63 0.95 1.42 2.13 3.19 (0.34-1.I6> (0.52-1 75) (0.77-2.61) (1.16-3.92) (1.73-5.87) 5merr. Hom Omernn (1977) Phosphorus Regional mooeis Nationwide mooe East Central Wes: % agricultural -r •« urcan Average p con: 67*« limns 67'ft Average p limns Average p 67‘« limns Average p 675, limns C 0.020 25 0.033 50 0.052 75 0.053 IOC 0.133 (0.010-0.042) (0.0)6-0.067) (0.025-0.106) (0.041-0.170) (0.065-0.272) 0.0)5 (0.008-0.0261 0.026 (0.0I4-0.O48) 0.045 (0.024-0.0S3) 0.078 (0.O42-0.144) C.136 (0.073-0.25!) 0.018 0.032 0.058 0.103 0.185 (0.008-0.037) (0.015-0.066) (0.028-0.120) (0.049-0.213) (0.089-0.382) 0.028 0.048 0.083 0.14] 0.240 (0.013-0.060) (0.022-0.103) (0.038-0.178) (0.065-0.303) (0.111-0.516) Mouret: Hom Omemik (19T7).
Table 2.2.3 Regional Nutrient Loadings (mg/L) (After Thomann and Mueller, 1987)
Plant Species Cold Weather Warm Weather Heteranthera dubie 2.13 m/s 1. 58 m/s Mvrioohvllum brasilience 2.80 m/s 1. 32 m/s
Table 2.2.4 Critical Scour Velocities for SRAV (after Vecchio, 1989)
Size Standard Sieve Class Millimeters : Inches: 4000-2000 160-80 Very large boulders 2000-1000 80-40 Large boulders 1000-500 40-20 Medium boulders 500-250 20-10 Small boulders 250-130 10-5 Large cobbles 130-64 5-2.5 Small cobbles 64-32 2.5-1.3 Very coarse gravel 32-16 1.3-0.6 Coarse gravel 16-8 0.6-0.3 Medium gravel 8-4 0.3-0.16 5 Fine gravel 4-2 0.16-0.08 10 Very fine gravel Microns: 2000-1000 18 Very coarse sand 1000-500 35 Coarse sand 500-250 60 Medium sand 250-125 120 Fine sand 125-62 230 Very fine sand 62-31 Coarse silt 31-16 Medium silt 16-8 Fine silt 8-4 Very fine silt 4-2 Coarse clay 2-1 Medium clay 1-0.5 Fine clay 0.5-0.24 Very fine clay
Table 2.3.1 Relationship of Particle Size and Sediment Class (After Brown, 1949)
Author of Sediment Discharge Relation Percent of Concentration Ratios Calculated/Observed in Range 1/2 to 2 Ackers and White [1973] 64 Engelund and Hansen [1967] 58 Rottner [1959] 53 Laursen [1958] 44 Einstein total load [1950] 44 Bishop et al [1973] 42 Einstein bed load [1950] 41 Graf [1971] 40 Bishop et al [1965] 39 Toffaleti [1968] 37 Einstein-Brown [1950] 35
Table 2.3.2 Reliability of Sediment Discharge Equations (After Dawdy and Vanoni, 1986)
Kode 1 ( Au Ihur ) Muir lent $ Modeled hut r1 ent 1 ones inorganic Mlirogei. i ui nt C A P SI OKIwd. Dhl ed. Or 94 n 1 c Part k Or 9411 lc Sed 1 - aenti Al 9*t 2 00- pl ankion Other Organ1 mh 3 mo ? HU j lota) Avail AQUA 1 V Cl -QUAl -KI LI I Ah Cl I AMI * KS Cl I ANl K bl A O0SAG3 (AM IS IL CO IIMLOMC-) mspf L AM CO Mil Artwork QUAL -11 KLC( IV- 11 SSAM IV MASK wo RS B Itrman J I 1 I 1 1 1 I I I J X X 1 X I X X X I X X X X X 1 X X X X X X J X X X X X X X 1 X X X X X X X X I X X X X X X I X X X X X X X X I X X X X X A X X X A X X X I X X X X X X X X > X X X X X X X X X X* I I p X I X* X X 1 X X X 1 X 1 1 X X X X X X X X X X X X X X X X X X X X X X X X I X X 1 X X X X X X X X I X X X X X X X X X X X X X X I X X X X X X X X X X X X I t I X X I X X X X X X X X X X X X X X Jorgensen L e lauan A ybo b» X x x X X X X X X X X X X X h X X X X X X X X X X X X I X X X X X Spec » 1 y flu
Table 2.4.1 Comparison of Nutrient Models (After USEPA, 1985)
Spatial domain Time domain State variable systems Model i; stream fl stream * - - s - —' L. — r AL'TO-QUAL AUTO-QD Bauer and Benneti DOSAG-1 DOSAG-IJI DOSC1 EXPLORE-! GENQUAL G475 HSPF LTM MIT-DNM MTJ-DNM Overton and Meadows PIONEER QUAL-1 Oual-ii RECEJV (SWMM) RECEIV-I] RJBAM RIVSCI SNSIM SMM SSAM-IV WASP WASP/SU1SAN W1RQAS (Veiz) WRECEV wqam c L. x ex <z> & • z ■ z • • X WQMM/Chowan WQRRS • • • • • • • • • -
Table 2.4.2 Comparison of Conventional Pollutant Models (After McCutcheon, 1989)
CHAPTER 3 MODEL DEVELOPMENT
3.1 ORIGINAL EUTRO4 PROGRAM DEVELOPMENT
EUTRO4 is a combination of the USEPA's Water Quality Analysis Simulation Program —4 (WASP 4 and their WASPB subroutine. WASP 4 contains a FORTRAN instruction set that simulates the movement and interaction of pollutants within the ecosystem (Ambrose, et al., 1988). The WASPB sub-model contains FORTRAN code to simulate eutrophication kinetics. The model is a simplified version of the Potomac Eutrophication Model developed originally by Thomann and Fitzpatrick (1982). Like most water quality models, EUTRO4 is based on the conservation of mass principle. The heart of the model uses a finite-difference approximation to compute changes in the mass balance equation. The solution procedure is discussed in the next section.
3.1.1 Finite-Difference Solution of Mass Balance Equation
EUTRO4 is based on a general mass balance equation written for an infinitesimally small control volume of fluid. Accounting for changes in concentration due to water column advection, pore water transfer, and dispersion, as well as physical, chemical and biological transformations (kinetics), the mass balance equation can be written as:
(H I j - TT = “17 (U » C) ’ “ C) ~77 (U - C)
+ -A_ (E AL) + A_ ( E AL) + rF aC 3 X 3x Sy y 3 y 8 Z Z d Z
4 + Sg + Sx
(3.1)
where
C = concentration (mg/L) t = time (days) U = water column velocity (m/s) E = diffusion coefficients (m 2 /day) S L = direct and diffuse loading rate (g/m 3 -day) S B = boundary loading rate (g/m 3 -day) S K = total kinetic transformation rate (g/m 3 -day) x,y,z = denote longitudinal, transverse, and vertical directions respectively.
The EUTRO4 model implements a finite-difference form of Equation 3.1 by expanding the control volume into larger segments and by specifying proper transport, loading, and transformation parameters (Ambrose et al., 1988). Simplifying the three-dimensional form of the difference equation by assuming vertical and lateral homogeneity and integrating over the depth and width (z and y respectively) leads to the following expression:
C) = — (-U x AC+ E x A + A (S L + S B ) *A S K at ax a* ♦
(3.2)
where A is the cross-sectional area (m 2 ) .
Taylor Series expansions of concentration can be performed if the derivatives of the concentration are single-valued, finite, and continuous functions of x (Smith, 1985). Assuming that these assumptions are valid, the resulting expansions are:
C(x.ix) - C(x) • AX SUX2. , t a s C(x I 4 2 ax* 6 ax*
(3.3)
C(x-Ax) = C(x) -Ax -^ c -- x ' ♦ - + ax* 6 ax*
(3.4)
Subtracting Equation 3.4 from Equation 3.3, neglecting higher order terms with exponents greater than two, and rearranging the result, leaves the following central-difference equation:
aC (x) _ C(x■*-Ax) - 0(x- Ax ) ax 2 Ax
(3.5)
Inserting the central-difference approximation shown in Equation 3.5 into the advection term of Equation 3.2 (the first term on the right hand side of the equation) leads to the following expression:
3CQ C) _ ax 2 ax
(3.6)
Likewise, substituting the central-difference approximation into the dispersion term (the second term on the right hand side of the equation) yields.
a [e x A |£] dx ” 2
(3.7)
The central-difference expressions for the derivatives on the right side of the above equation can be written as:
3CIx + A x ) _ C(x +2A x ) - C (x) 2 Ax aC(x-Ax) = C(x) - C(x-2Ax) ax 2 Ax
which can be inserted into Equation 3.7 to yield:
i x sxl E(x+Ax) A(x+Ax) £ ■ x Ctx ) S J = 2 ax 2 ax
E(x -Ax ) A(x -Ax ) 2 ax 2 Ax
(3.8)
A general mass balance equation can be created by substituting Equations 3.6 and 3.8 into Equation 3.2 for the appropriate terms. To convert the resulting mass balance equation into an expression which reflects the segmented nature of EUTRO4, a notation for describing the locations within the segments must be developed. This notation can be partially explained by examining the diagram shown in Figure 3.1.1. As illustrated in the figure, the center of the segment, x o , can be replaced by the subscript ”j”. Similarly, the downstream and upstream interfaces between segments can be represented by ”j,j+l" and respectively, while the center nodes of the downstream and upstream segments are designated as ”j+l” and "j-1”. The new mass balance equation for segment j can be written:
at 2 Ax 2 Ax
(E (E A) i-! ,K - I I / J - ' J 2 &x 2 Ax ’
+ A j (S L + Sg + S K )j
(3.9)
Rearranging the terms in Equation 3.9 and defining several new variables leads to:
a(V . C,) H = + Qj-i.j c j-uj + - Cj]
Rj-i-J [Cj ' + V j * S B + S K )j
(3.10)
where
V = 2 Ax A = segment volume (m°)
R = (E A)/(2 △x) = dispersive flow (n 3 /day).
Notice that Equation 3.10 contains several references to the concentration at the segment interfaces (the use of subscripts such as j, j+l denote an interface). EUTRO4 only computes average concentrations within the segment. Therefore, estimates of interfacial concentrations must be made. This is accomplished using the following:
Cj,j+l “ v £j + i * d Cj
(3.11)
v + Cj-l
(3.12)
where
V = advection factor (between 0 and 1).
Neglecting program stability the advection factor can be assigned any value between zero and one. If the advection factor is set to zero, the backward-difference equation is used for the advection term. If the advection factor is set to one-half, the centraldifference equation is used for the advection term. The effect of this parameter on model stability is discussed in a subsequent chapter. For the time being, however, it should be pointed out that the WASP 4 Users Manual recommends a value of between 0.0 and 0.4 for this parameter (Ambrose et al., 1988).
Equation 3.10 is a simplified form of the mass balance equation used by EUTRO4; however, it can easily be expanded into the multi-dimensional form employed by the model. For the situation where "i” segments adjoin segment ”j” and segment interfaces are expressed as ”ij”, the general mass balance equation becomes:
=- V Q,J C,j . Z R >J (C > - C J> •S Vj Slj c T j j x-
(3.13)
• ZVj S BJ .S Vj Skj b K
where
q a = flow, positive when leaving segment j, and negative when entering j (m 3 /day).
This equation is used to calculate the mass derivative for each stream and sediment segment. Given the initial concentrations and segment volumes at some initial time, EUTRO4 computes new constituent masses at the next time step using the one-step Euler method.
This equation can be expressed as:
(Vj = (Vj Cj) t . M 3 t
(3.14)
Finally, the new constituent masses calculated using Equation 3.14 are divided by the new segment volumes to determine the concentrations. This procedure is then repeated for each simulation time step.
3.1.2 Applicability of One-Dimensional Model Formulation
Strictly speaking this application of EUTRO4 is formulated as a two-dimensional (width averaged) model since it contains segmented layers in the 2-direction (vertical). From a more pragmatic perspective, however, the water column is modeled as a single layer forcing the model to behave like a one-dimensional formulation. Therefore, discussion of EUTRO4 limitations will treat the program as if it were a one-dimensional model.
One-dimensional approximations of river systems are generally believed to be sufficient for water quality modeling purposes (Renaldi, et al., 1979). However, severe errors may occur if the limitations and consequences of one-dimensional representations are not carefully considered (Ward and Fischer, 1971). To determine the applicability of a one-dimensional model, one must calculate the transverse mixing length (McCutcheon, 1989). The transverse mixing length can be defined as the distance downstream from a source required for the source to mix with the receiving waters
such that the resultant stream concentration is essentially uniform throughout the entire cross-section.
The best way to determine the transverse mixing length, X, in a river is to conduct dye studies over the range of flows being considered in the project. However, in the absence of field measurements, Holley and Jirka (1986) present the following predictive equation first developed by Fischer:
x = Xd V 52 - x d V B 2 e y ' <x y H(g R S f ) 0 5
(3.15)
where
X = distance from source to point where discharge has been well mixed laterally (m) H = average stream depth (m) V = average stream velocity (m/s) B = stream width (m) g = gravity (9.81 m/s 2 ) R = hydraulic radius S f = friction slope (m/m) X d = dimensionless longitudinal distance = dimensionless mixing coefficient.
By making several assumptions primarily concerning the hydraulics of the problem, a rough estimate for the transverse mixing distance can be computed using Equation 3.15. For a discharge point located on the bank of the stream (like the effluent pipe at the South Austin Regional WWTP) the distance required for ninety percent mixing in a flow of 45.34 cms (1600 cfs) is.
Y _ Q. 3*o. 3Q3 *(84. 7) 2 _ 24.3 km = 15 miles 0.5*1.766*[9.81*1.6953*0.0000554] °- 5
This distance indicates that the mixing zones of the major point source loads (Walnut Creek, Onion Creek, South Austin Regional Waste Water Treatment Plant, etc.) are well within the study area. Therefore, it will be assumed that a one-dimensional problem formulation will produce satisfactory results. However, this mixing length should be considered when designing monitoring network, sampling water quality parameters, and evaluating the results of point source samples vs segment averaged concentrations generated by EUTRO4 or other water quality models. Moreover, exact criteria concerning the applicability of a one-dimensional models should be explored further if the area of interest is close to any significant point source.
3.2 MODIFICATIONS TO EUTRO4 PROGRAM
Johns (1988) modified the original eutrophication portion of the pre-WASP3 computer program to include equations to predict the potential growth rate four species of aquatic plants and the steady state suspended sediment concentration. As part of this research, these equations were further modified to be compatible with the current WASP4VI4 code and WASPB subroutine. To distinguish the updated model from EPA's original eutrophication model, the name EUTRO-M will be used to designate the new macrophyte capabilities.
3.2.1 Final Nutrient-Oxygen-BOD-Phytoplankton Kinetic Equations
The final three terms of the mass balance equation (Equation 3.1) account for the biological, chemical, and physical processes occurring in rivers. Collectively, these processes are often referred to as the kinetics term. The nutrient, oxygen, BOD, phytoplankton kinetic equations contained in EUTRO-M represent the complex interactions between systems which occur in natural systems. Figure 3.2.1 illustrates these basic linkages between the different ecological components of the Colorado River. While Figure 3.2.1 is practical for illustrative purposes, the schematic linkage of system parameters shown in Figure 3.2.2 may be more useful in understanding the modeling approach used to represent the interactions between the state variables.
The kinetic equations and descriptions presented below are for water column computations of the nine state water quality constituents. For most parameters, the equations are essentially the same for the sediment layers except that losses due to settling in the water column are generally considered as sources to the benthic sediment layer.
To avoid confusion and unnecessary duplication, the following definitions are used repeatedly throughout the development of the kinetic equations. Namely,
C x = ammonia concentration (mg/L) q — nitrate concentration (mg/L) C 3 = orthophosphate concentration (mg/L)
C 4 = phytoplankton concentration (mg/L) C 5 = carbonaceous biochemical oxygen demand (mg/L) C 6 = dissolved oxygen concentration (mg/L) C 7 = organic nitrogen concentration (mg/L) C 8 = organic phosphorus concentration (mg/L) C 9 = total suspended sediment concentration (mg/L) C IO = epiphytic algae potential (gm/m 2 ) C XI = macrophyte species no. 1 potential (gm/m 2 ) C l 2 = macrophyte species no. 2 potential (gm/m 2 ) C l 3 = macrophyte species no. 3 potential (gm/m 2 ) H = depth of segment (m) T = temperature (°C) t = time (days).
Other variables which are used in the kinetic equations are defined as needed.
It should also be point out that the kinetic equations presented subsequently are general in nature, meaning that they apply to both the surface water segments and the benthic sediment layer segments. The appropriate terms in the equations are set to zero whenever the rules governing the kinetics of the situation dictate. For example, the kinetic equation for ammonia indicates that the concentration is a function of epiphytic algae growth, SRAV growth, nitrification, mineralization, and phytoplankton growth. However, in surface water segments the concentrations of the SRAVs are zero (rooted vegetation assumed to grow only in sediment layer) so these terms drop out of the equation. Likewise, since epiphytic algae and phytoplankton grow only in the surface water segments, when computing the ammonia concentration derivative in the sediment layers, the algae sinks drop out.
The difficulty in this procedure stems from the fact that several water column parameters are required to calculate the light extinction coefficient during the computation of macrophyte growth. Consequently several new program variables were added. Program modifications were made which stored water column parameters such as segment depth and velocity as well as phytoplankton, suspended sediment, and BOD concentration until they were needed.
1. Ammonia (NH3I
Four types of nitrogen are included in the model: (1) ammonia, (2) nitrate, (3) phytoplankton nitrogen, and (4) organic nitrogen. The linkages and transformations of these species were presented in Figure 3.2.2. The first state variable is ammonia (NH 3 ) . As shown in Equation 3.16, the main source of ammonia is from the mineralization of organic nitrogen ( see system 7). The terms on the right hand side of the equation account for mineralization of organic nitrogen, phytoplankton growth, nitrification, phytoplankton death, and epiphytic algae and macrophyte growth, respectively.
= *k 71 e?Z 20 X prc C 7 - G pl anc 8 t
, .T-20 I Cl • Dpi Cc ano (l-fon) k l 2 612 [Knit * C 6 J P
UI 10 C lO BMON gj Uln On BIIIN £ 2 H ' H
UI 12 C l 2 BII2N e 2 UI 1S Cis BU3N 62 H ’ H
(3.16)
where
a nc = nitrogen to carbon ratio (mg N/mg C) Dpi = phytoplankton loss rate (1/day) G pl = phytoplankton growth rate (1/day) k l 2 = nitrification rate (1/day) k 7l = organic nitrogen mineralization rate (1/day) BIION = nitrogen/total cell mass fraction for epiphytic algae BIIIN = nitrogen/total cell mass fraction for Heteranthera dubia BII2N = nitrogen/total cell mass fraction for Myriophyllurnbrasilience BII3N = nitrogen/total cell mass fraction for Potamogeton pectinatus UI 10 = growth rate of epiphytic algae (1/day) UI 11 = growth rate of Heteranthera— dubia (1/day) UI 12 = growth rate of Myriophyllum brasilience (1/day) UI 13 = growth rate of Potamogeten pectinatus (1/day) fon = fraction of dead phytoplankton returned to organic nitrogen pool K t = half saturation constant for oxygen limitation of nitrification (mg 0 2 /L) X = phytoplankton mineralization coefficient = ammonia preference factor for phytoplankton 3 = ammonia preference factor for epiphytic algae and macrophytes, respectively @ l2 = temperature coefficient e = temperature coefficient.
Explicitly, the ammonia preference factor for phytoplankton can be represented by:
= — —-—— - + * CiHKhN * C 2) (Gj ♦ 63) tK NN - Cz)
(3.17)
where
- Michaelis ammonia limitation factor (mg N/L).
Similarly, the ammonia preference factor for the macrophyte and epiphytic algae species is given by:
p _ PREF Ci [PREF Ci (1-PREF) C 2 ]
(3.18)
where
PREF = ammonia preference factor (between 0 and 1).
The effects of phytoplankton mineralization are computed using:
y - Xprc K nPC *C 4
(3.19)
where
= half saturation constant for phytoplankton limitation of phosphorus recycle (mg C/L)
2. Nitrate (NO3)
The second state variable modeled by EUTRO-M is nitrate nitrogen. Nitrate nitrogen is created by a twostep oxidation process whereby ammonia is transformed first to nitrite and then to nitrate. This process, called nitrification, is a complex phenomenon depending on dissolve oxygen, pH and flow conditions (Ambrose et al., 1988). Nitrate sinks can be attributed to growth
of aquatic vegetation and denitrification. An equation based on a simplification of the processes is included in the EUTRO-M model. The equation used is:
■ -f'' - *ki2 p —p — C’l - (1-Pnh 3 ) Snc at l-Knit + u 6 J
, , T-20 - K2O 620 -p ? — 3 +
Ulio Cjo BIION (1-Pi) Ulu C tl BIIIN (l-g 2 ) H ' H
(3.20)
UII2 Cl 2 BII2N (l-g 2 ) Ulis Cis BII3N (l-g 2 H ’ H
where
K 2O = denitrification rate (1/day) K NO3 = Michaelis denitrification constant (mg 0 2 /L) e 2O = temperature coefficient.
3. Orthophosphate (0P04)
Three types of phosphorus are modeled: (1) orthophosphate (inorganic), (2) phytoplankton phosphorus, and (3) organic phosphorus. Orthophosphate is the type used directly by the aquatic vegetation. It is a function of the mineralization rate, vegetation growth, and settling. The equation used in EUTRO-M to model orthophosphate kinetics is:
3Cg , T-2O xr i* r' c D 3 p —7— “ + kgs ©B3 " a pc " m O L
UI 10 C lO BIIOP Uln Ch BIIIP H ' H
UI 12 c l 2 BII2P UI 13 Cis BU3P H ' H
(3.21)
where
a pc = phosphorus to carbon ratio (mg P/mg C) f D3 = fraction of dissolved inorganic phosphorus in the water column k B3 = dissolved organic phosphorus mineralization rate (1/day) V SS = inorganic sediment settling velocity (m/day) BIIOP = phosphorus/total cell mass fraction for epiphytic algae BIIIP = phosphorus/total cell mass fraction for Heteranthera dubia BII2P = phosphorus/total cell mass fraction for Myriophyllum brasilience BII3P = phosphorus/total cell mass fraction for Potamogeton pectinatus e 8 3 = temperature coefficient.
4. Phytoplankton
Phytoplankton is modeled as the fourth state variable in EUTRO-M. Phytoplankton assumes a critical role in the eutrophication process because it interacts with most of the other systems. For internal computational purposes, the model uses phytoplankton carbon as a measure of algal biomass. However, for calibration or external purposes, chlorophyll a values are calculated using a carbon to chlorophyll ratio.
Similarly, the effect of phytoplankton on the nutrient species is modeled using the following two equations:
31C 4 Sp t ) f V s 4) = - Upl '—J ape
(3.22)
3( C 4 a n c ) f ~ n 1 a = b pl " -Upl ' ~pj d nc
(3.23)
where
V S4 = settling velocity of phytoplankton (m/day).
Notice that both the growth and death rates of phytoplankton in Equations 3.22 and 3.23 are each represented by single terms. In fact, phytoplankton growth and death rates are functions of many complex interactions. These were discussed in some detail in Chapter 2 and will not be repeated here.
5. Carbonaceous Biochemical Oxygen Demand (CBOD)
The oxidation of carbonaceous material is modeled as the fifth state variable in the EUTRO-M program. In addition to man-made sources, the primary sources of CBOD are those contributed by detrital phytoplankton carbon, epiphytic algal carbon, and SRAV carbon. The equation used in EUTRO-M is:
K r , n T-20 r TT - - d oc k D 7“ oX KB Q £ + U g
V ss (1-F D 5) 5 32 , Knq, H Cs ■T H k2c e2c K NOa •C 6 Cs
(RESP Io *MORT lo ) C lO BIIOC (RESPu'MORTu) On BIIIC H * H
(RESP 12 *M0RT 12 ) C l 2 BII2C (RESP 12 *MORT 1S ) C IS BII3C * ~ H * H
(3.24)
where
a oc = oxygen to carbon ratio (mg 0 2 /mg C) f DS = fraction dissolved CBOD = half saturation constant for oxygen limitation (mg 0 2 /L) k D = deoxygenation rate (1/day) k ID = non-predatory phytoplankton death rate (1/day) k 2D = denitrification rate (1/day) V S3 = organic matter settling velocity (m/day) BIIOC = carbon/total cell mass fraction for epiphytic algae BIIIC = carbon/total cell mass fraction for Heteranthera dubia BII2C = carbon/total cell mass fraction for Myriophyllum brasilience BII3C = carbon/total cell mass fraction for Potamogeton pectinatus MORT lo = epiphytic algae mortality (g/m 2 -day) MORT 1X = Heteranthera dubia mortality (g/m 2 -day) MORT 12 = Myriophyllum brasilience mortality (g/m 2 -day) MORT 13 = Potamogeton pectinatus mortality (g/m 2 -day) RESP 1O = epiphytic algae respiration (g/m 2 -day) RESP 1X = Heteranthera dubia respiration (g/m 2 -day) RESP 12 = Myriophyllum brasilience respiration (g/m 2 -day) RESP 13 = Potamogeton pectinatus respiration (g/m 2 -day) = temperature coefficient © 2D = temperature coefficient.
6. Dissolved Oxygen (DO)
Dissolved oxygen is an important component of any water quality model. It plays a critical role in the oxidation processes. Again, dissolved oxygen is a function of many complex interactions. In EUTRO-M dissolved oxygen is modeled by the following equation:
- l oT-zo ri b , -ka ©a (Cs CgJ kc ©D 1/ x r at + J
64 , U-20 r 32 T r SOD "TT k l 2 e l 2 . Cj .j Ci - 12- k IB ei R C 4 - —
r f 32 48 a nc n J 2.667 UI 10 C lO BIIOC + IyY + 14 11-r N H 3 J | + h
2.667 Uln Ch BIUC 2.667 UI 12 C l 2 BU2C + H + H
2.667 UI 1S C IS BII3C 1.87 RESP 10 C lO BIIOC H ’ H
1.87 RESPn Cu BIUC 1.87 RESP 12 C l 2 BII2C H ' H
1.87 RESP ls C IS BII3C H
(3.25)
where
CS = saturation concentration for oxygen (mg/L) ka = reaeration rate (1/day) kIR = algal respiration rate (1/day) SOD = sediment oxygen demand (g/m 2 -day).
7. Organic Nitrogen (ON)
Organic nitrogen is modeled as the final species of nitrogen. The primary sources of organic nitrogen in the kinetic eguation below come from the death of aquatic vegetation (phytoplankton and SRAV's). Organic nitrogen must undergo mineralization to ammonia before the aquatic vegetation may use it as a nutrient.
~~~T~ “ * Dpi D 4 a nc fort - k 7l ©71 2 ° D 7 o X
V sS (RESP Io+ MORT lo ) C lO BIION h 4 n
(RESPn+MORTu) C n BIIIN (RESP 12 *M0RT 12 ) C l 2 BII2N + H + H
(RESP ls +MORT 13 ) C lg BII3N H
(3.26)
8. Organic Phosphorus (OP)
Organic phosphorus is modeled as the eighth state variable in EUTRO-M. Like organic nitrogen, the primary sources are functions of the death and respiration of phytoplankton, epiphytic algae, and macrophytes. Organic phosphorus must undergo mineralization to orthophosphate before it can be utilized as a nutrient for the aguatic plants. This process is shown as a sink in the following eguation:
- + Dpi D 4 apc “
V (I'Fpg) (RESP io*MORT iq) Oto BIIOP ' H Ce * H
(RESPn+MORTu) Cu Bill? (RESP 12 *M0RT 12 ) C l 2 BII2P H * H
(RESP 1S +MORT IS ) C IS BII3P H
(3.27)
9. Total Suspended Sediment (TSS)
There are a number of reasons for modeling suspended sediment in a river system. Two of the main reasons in the lower Colorado River basin are (1) the shading effects of suspended sediment and (2) the adsorption of nutrients to sediments. Shading by the suspended sediment reduces the amount of light that can penetrate the water column. This influences the growth rate of macrophytes, epiphytic algae, and, to some extent, the phytoplankton. Nutrient adsorption is important because phosphorus (and to a lesser extent nitrogen) has a strong tendency to sorb to sediments and therefore its fate is related to the transport of the sediment.
Suspended sediment is not included as a state variable in the EUTRO4 eutrophication kinetics provided by the USEPA. However, because of its importance in determining the overall availability of light, the WASPB subroutine was modified such that total suspended sediment is now included as a state variable (System 9).
In his pre-WASP3 study, Johns (1988) calculated total suspended solids concentrations using a formula based on water column velocity and settling and resuspension characteristics of the suspended particles. This equation was modified and included in EUTRO-M so that the change in sediment concentrations for steady flow problems is calculated using the following derivative equation:
aCg _ Vel * Resup TSS * VS at " 1000 * H ” H
(3.28)
where
Vel = flow velocity (m/s) Resup = resuspension rate (mg/day)/(m 2 )/(m/s) H = average water depth (m) TSS = total suspended sediment (mg/L) VS = settling velocity (m/s).
This is not the only means of computing suspended sediment discharge. Numerous theoretical and empirical sediment transport eguations based on dimensional analysis, laboratory flume studies, models, and field measurements can be found in the literature. It would not be practical or prudent to attempt, here or in the literature review, to discuss all of the methods. Based on the limitations in the availability of field data, this simple method seems appropriate for the current steady state investigations.
Although sufficient for steady state analysis, Equation 3.28 is not appropriate for dynamic flow simulations. The rapid changes in flow rates do not permit the suspended sediment concentrations predicted by Equation 3.28 to reach steady state. Therefore, the calculated concentrations would usually be incorrect.
where
C = sediment concentration (mg/L) Q = flow rate (cms) a = regression coefficient b = regression coefficient e = Monte Carlo Uncertainty coefficient (mg/L).
At first Equation 3.29 may seem like a rather simplistic relationship. The model could have an option that would allow the user to select from several of the suspended sediment transport equations discussed in Chapter 2. However, the data requirements and accuracy of these equations do not currently justify the added confusion and effort.
The regression coefficients "a” and "b” are specified for each model segment. These coefficients were determined by using Equation 3.28 to predict total suspended solids concentrations over a wide range of steady state discharges and then performing linear regression on the results for each stream segment. Field measurements of sediment concentration versus flowrate were used to calibrate the parameters used in the equation.
This method was chosen over the suspended sediment relationship developed by Harris (1990). Harris found that under “normal” flow conditions suspended sediment concentrations on the lower Colorado River could by approximated by the following equation:
Ml(x) = 123 [1 - exp(-0.024X)]
(3.30)
where M1(X) = suspended sediment concentration (mg/L) X = distance downstream of Longhorn Dam (km) 123 = estimated low flow carrying capacity of river (mg/L) 0.024 = net resuspension coefficient (1/km).
The application of this equation is limited due to the fixed value of the carrying capacity. As pointed out by Celik and Rodi (1991), the carrying capacity of a stream is a function of stream velocity and shear stress, and therefore, should not be assumed constant over a wide range of unsteady discharges. (3.29) Given the constraints of the suspended sediment relationship, the alternative proposed here is to use a widely accepted sediment rating curve technique. Plots of suspended sediment concentration versus stream discharge often reveal a logarithmic relationship. The sediment transport function currently being used is: Ideally a linked sediment/water transport model could be used to calculate unsteady suspended sediment concentrations. However, most all of the available sediment transport models are based on lumped total sediment load equations (Holly, 1988). To use these types of models to predict suspended sediment concentrations, an empirical relationship would need to be developed. Given the amount of variability which exists naturally, this relationship would be rather subjective.
C = a*Q b + e
3.2.2 Macrophyte and Epiphytic Algae Kinetics Equations
There are literally hundreds of different genera of algae and macrophytes in the world today. Bauman (1978) found ninety-three genera of algae alone in his survey of streams and ponds around the Austin, Texas area. Seagle (1977) found five types of epiphytic algae (Characium, Closterium, Palmellococcus, Scenedesmus, and Stigeoclonium) at sampling locations on the Comal and Guadalupe Rivers. Bauman (1978) reported finding all but Palmellococcus on one or more occasions on Onion, Shoal, Waller, and Williamson Creeks. For the purpose of this investigation, all genera of attached algae have been lumped into a single category called epiphytic algae.
The species of SRAV being studied can be a very important factor in developing the mathematical expressions for macrophyte growth. For example, studies performed by Chambers and Kalff (1987) indicated species specific effect of light and sediment composition. In their investigation, growth of erect plants such as Potamogeton pectinatus depended more on the availability of light while the growth of canopy producing species such as Myriophyllum spicatum was influenced more by
tamogeto
sediment composition. The plant species considered to be the most important on the lower Colorado River are epiphytic algae, Heteranthera dubia, Myriophyllum brasilience, and Potamogeton pectinatus. All four plant species were modeled using the same basic eguation with different coefficients.
The equation used to predict SRAV and epiphytic algae growth potentials in the EUTRO-M model can be written as:
=(G-R-S-M)«C a at
(3.31)
where
C A = plant density per unit area (gm/m 2 ) G = plant growth rate (1/day) R = plant respiration rate (1/day) S = loss of plants due to scour (1/day) M = loss of plants due to mortality (1/day).
The interdependency of Equation 3.31 to the other state variables can be seen by examining the four terms inside the parenthesis on the right hand side of the equation. Each of the variables are functions of parameters which are themselves functions of still other factors. For example, the plant growth rate, G, can be represented by:
G = Gmax * f(T) * f(N) * f(P) * f(L)
(3.32)
where Gmax = maximum growth rate (1/day) f(T) = temperature factor f(N) = nitrogen factor
f(P) = phosphorus factor f(L) = light factor.
It must be pointed out that Equation 3.32 contains most of the important parameters influencing macrophyte growth which can be reasonably quantified in a large scale field environment. As such, the macrophyte kinetic equations used in EUTRO-M are typical of those used to model other aquatic plants. It should not be assumed however, that these are the only parameters which affect SRAV growth. For example, as a result of laboratory analysis, Barko and Smart (1986) found that macrophyte growth decreased substantially as sediment organic matter increased and when the sand content of the sediment was more than seventy-five percent. Moreover, Duarte and Canfield (1990) found that the amount of canopy cover significantly influenced the growth of macrophytes.
The temperature factor, f(T), is a two-tailed function whose value ranges from zero at low temperatures to nearly one over an optimal temperature range and then back to zero for extremely high water temperatures. The function can be calculated using the following equation:
mx - T) f( T ) = A T~ * A T“ 1 +Ki > _l| 1+K 4 [e V2<T max ' T> -1]
(3.33)
where
1 T Y 1 "f T T i Ln Ki(l-K 2 ) * cpt 1 * min j
1 , k s (l-k 4 ) Y 2 = “ “TT k 4 (i-k s f p hiax 1 opt 2 J
T optl = lower limit of optimum temperature (°C) T opt2 = upper limit of optimum temperature (°C) K x = multiplier near lower limit of growth Tmin K 2 = controls rate of asymptotic approach of f(T) K 3 = controls rate of asymptotic approach of f(T) K 4 = multiplier near upper limit of growth Tmax Yi = coefficient for lower side of curve (1/°C) Ys = coefficient for upper side of curve (1/°C).
As in this investigation, K 2 and K 3 are often assumed to have the value of 0.99.
The nitrogen factor, f(N), is formulated as a Michaelis-Menton relationship involving both ammonia and nitrate nitrogen. The equation used in EUTRO-M to determine the macrophyte/epiphytic algae nitrogen factor is:
r f N s _ ft * * Cl-g) * Cn3 * Cni + (1-p) * Cns + Kn)
(3.34)
where
K n = Michaelis-Menton constant for nitrogen (mg/L) C NI = ammonia nitrogen concentration (mg/L) C N3 = nitrate nitrogen concentration (mg/L) = ammonia preference factor.
Similarly, the phosphorus factor, f(P), is also written as a Michaelis-Menton relationship. The equation is:
C p i f(P) = ~ (Cpi + KpJ
(3.35)
where
C PI = available phosphorus concentration (mg/L) K P = Michaelis-Menton coefficient (mg/L).
The light factor, f(L), averages the incident light over the water column depth and over the course of a 24- hr day. Its value also varies between 0 and 1. The equation used to represent this factor is given by:
PHOTO , [ *lq J ex H * r K L +lo e>< p(
(3.36)
where
PHOTO = photoperiod (fraction of day) H = average water depth (m) K l = Michaelis-Menton light constant (Ly/day) I o = average surface light intensity (Ly/day) = total light extinction coefficient (1/m)
Each individual component of the light extinction coefficient is accounted for separately in the model using a modified form of the formula developed by DiToro (1978). In his study, DiToro found that organic detritus and algal chlorophyll contribute to the light extinction coefficient mainly via absorption while suspended solids contribute by both absorption and scattering. The equation used in EUTRO-M accounts for the self-shading characteristics of the aquatic vegetation and the light scattering properties of pure water. The overall extinction coefficient used in the model is calculated using the following equation:
K. = K„ + K es *TSS + K a , trI *CBOD + K pbyto * PHYTO
+ K shade *SßAV
(3.37)
where
K e = total light extinction coefficient (1/m) K w = pure water light extinction coefficient (1/m) K ss = coefficient for suspended sediment (m 2 /mg) TSS = total suspended sediment (mg/L) K detri = coefficient for organic material (m 2 /mg) CBOD = carbonaceous biochemical oxygen demand (mg/L) K phyto = coefficient for phytoplankton (m 2 /mg) PHYTO = phytoplankton concentration (mg/L) K shade = macrophyte self-shading coefficient (m/gm) SRAV = macrophyte density (gm/m 2 ).
Respiration can be thought of as the inverse of photosynthesis. It is the process whereby organic matter and oxygen are converted back into inorganic carbon with the release of chemical energy (Jahnke, 1990). Using several of the factors previously defined, the plant respiration rate, R, can be given by:
R = f(T)*[ R min + {K E *f (N) *f (P) *f (L) ) ]
(3.38)
where
R = respiration rate of plants (1/day) R min = minimum respiration rate (1/day) K R = maximum respiration constant (1/day)
At high water velocities, macrophytes can lose stems and leaves or the entire plant may be removed which results in a decrease in plant density. This process is called plant scour, S, and was included in the model using the following equation developed by
Vecchio (1989) :
S = K (1 - SRAV/ (K* + SRAV) )
(3.39)
where
K = coefficient for plant scour (1/day) K* = Michaelis-Menton protection factor coefficient (gm/m 2 ) SRAV = submerged rooted aguatic vegetation density (gm/m 2 ).
Vecchio's results may be biased toward conservative estimates of macrophyte scour by the very nature of the plants sampled for laboratory analysis. Plants growing in the field may already be devoid of tender leaves and branches which would easily break off under lower flow velocities. Wright and McDonnell (1986 b determined estimates of advective losses by placing screens upstream and downstream of macrophyte beds and measuring the guantities of biomass scoured from the river segment directly. Their investigation did not indicate a threshold velocity. The resulting kinetics equation therefore permitted macrophyte scour even at low velocities. Nevertheless, since Vecchio sampled and tested plants actually found in the river basin, his methodology has been incorporated into EUTRO-M.
As alluded to above, the scour term is only implemented when the water column velocity exceeds the threshold water velocity for a particular plant species. By performing experiments on macrophyte samples taken from the lower Colorado River, threshold scour velocities for Heteranthera dubia and Myriophyllum brasilience were determined in the 1989 Vecchio study
and were presented in Chapter 2 of this dissertation. Since no value was determined for Potamogeton pectinatus, an average of the threshold velocities of the other two SRAV species was used. Vecchio (1989) also determined that the values for K and K‘ should be 9600/day and 0.065 gr/m 2 , respectively.
Klasvik (1974) found that stream velocities greater than 1.0 m/s severely retarded or prevented the growth of attached algae on artificial surfaces placed on the stream bed. Therefore, this value was used for the epiphytic algae threshold velocity.
The plant mortality rate, M, was assumed to be a function of temperature. The equation can be formulated as:
M = Mmax * f'(T)
(3.40)
where
Mmax = maximum mortality reached in optimal temperature range (1/day) f'(T) = identical to temperature factor for temperatures less than T opt2 . Above T opt2 value remains at 1.
The final EUTRO-M model is capable of handling the various exchanges between a multi-layered water/sediment system. An example of this type of system is shown in Figure 3.2.3. The model assumes that epiphytic algae grow only in the water column segments (obtain nutrients from the water column) and SRAV's grown only in the upper bed segments (obtain nutrients from upper sediment layer). The growth kinetics routines are only called
for the appropriate segments thereby increasing computational efficiency.
The EUTRO-M program was compiled, debugged, and run using Lahey's extended memory F77L-EM/32 FORTRAN compiler 1 and an IBM PS/2 Model P7O PC 2 . A copy of the updated eutrophication kinetics subroutine (WASPB) is presented in Appendix A. The program requires three distinct ’’include” files for proper execution. The ’’include" files contain most of the common blocks and variable declaration statements used in the program. Copies of the final modified include files are presented in Appendix B.
The temperature, nitrogen, phosphorus and light growth coefficients presented in Equation 3.32 and explained in the subsequent text are complex functions involving many different variables and program constants. While the precise values of these coefficients cannot be determined without benefit of a simulation, it is possible to demonstrate the general shapes and characteristics of the growth coefficients. By assuming typical values for all but one important parameter in each of the coefficient equations, the identifiable trends can be plotted. These trends can be seen in Figures 3.2.4 and 3.2.5.
Figure 3.2.4 illustrates two interesting
phenomenon. The phosphorus curves show how little of the nutrient is needed before the related growth coefficient approaches its maximum. The nitrogen (ammonia) curves are relatively flat indicating the plants ability to substitute nitrate for ammonia. The nitrogen growth coefficient goes to zero only when the nitrate concentration is removed from the system.
The temperature and light coefficients are shown in Figure 3.2.5. The temperature coefficient curve is almost flat near the optimum temperature range but the tail drop off steeply at high and low temperatures. The light coefficient curve illustrates how water depth and sediment concentration influence the magnitude of the growth coefficient.
1 F77L-EM/32 is a trademark of Lahey Computer Systems, Inc.
2 PS/2 Model P7O is a trademark of IBM
3.3 DEVELOPMENT AND CALIBRATION OF UNSTEADY DAFLOW MODEL
To quantify the effects of unsteady flow on the water quality environment, it is necessary to determine the flow rate (and subsequently the water column velocity) in each model element. As currently written, WASP 4 itself does not actually route the stream discharge. Instead, the program allows the stream flow to be specified an a series of time dependent functions. The program uses the step like function to change the flows in all the downstream segments at once. However, as an option, flows for each segment can be read into the program from a data file.
Unlike the steady state case, wave propagation
causes the flow rate to change as the hydrograph travels downstream. Since gage flows are not available at each segment, an unsteady flow routing model had to be used to estimate these values. For this investigation, the USGS model known as DAFLOW (diffusion analogy) was selected (Jobson, 1989). A number of other suitable models exist, including EPA's DYNHYD and the National Weather Service's DWOPER and FLDWAV models, making the decision on which model to use somewhat arbitrary. Ultimately, factors such as the author's familiarity with DAFLOW, the model's stability, and the relatively simple data requirements of the model, led to the selection of the DAFLOW model.
3.3.1 Theory of Diffusion Analogy Model
The following section highlights the basic principals of the unsteady flow modeling technique used by the diffusion analogy flow model (DAFLOW). A comprehensive discussion concerning numerical solutions of the Saint-Venant equations is beyond the scope of this report. For more information on the subject the reader is referred to the following references (Chow et al., 1988; Strelkoff, 1970; Saint-Venant, 1871; Lai, 1986) .
Jobson (1989) uses the differential equations developed by Barre de Saint-Venant in 1871 as the nucleus of the one-dimensional unsteady flow model. Assuming no lateral inflow, the conservative forms of the Saint-Venant continuity and momentum equations can be written respectively as:
aQ + a A = n ax at
(3.41 a
1 3 U + U 3 U + 3 z + q _ _ o g St g ax 3x
(3.41 b
where
A = area of flow (m 2 ) g = gravity (m/s 2 ) Q = volumetric rate of flow (m 3 /s) S f = friction slope (m/m) S o = bed slope (m/m) U = velocity (m/s) x = longitudinal distance along channel (m) z = depth of flow (m).
No analytical solution to these equations exist except for cases where the channel geometry is uniform and the nonlinear properties of the equations can be neglected or linearized (Strelkoff, 1970; Jobson, 1989). As a result, numerical solutions are required. The complexity of these solutions depends on the number of terms needed in the momentum equation to accurately model the system.
The five terms in the momentum equation from left to right represent the (1) local acceleration, (2) convective acceleration, (3) pressure force, (4) friction force, and (5) gravity force. There are three general classes of solutions to the momentum equation (Chow et al., 1988). A solution which includes all five components is called a dynamic wave model. Neglecting the local and convective acceleration terms and using the remaining three forces results in a diffusion wave
model. Models which include only the gravity and friction force terms are called kinematic wave models.
For most river systems not influenced by tidal action or severe backwater effects, the first and second terms of the momentum equation can be neglected without impacting the accuracy of the solution (Jobson, 1990). Therefore, the diffusion analogy model developed by Jobson (1989) solves the diffusion wave form of the momentum equation. However, since the model makes several simplifying assumptions, the designation diffusion analogy rather than diffusion wave model is used.
As an example of the relative magnitudes of the terms in Equation 3.41 b the following example is presented. If the discharge in the lower Colorado River is increased from 45.34 cms (1600 cfs) to 90.68 cms (3200 cfs) then the water column velocity in Segment 38 (segment with the highest velocity) increases from 0.62 to 0.88 m/s. Assuming that this change takes an hour, the value of the local acceleration term is:
= —i—- | — = 7.25 E-06 *O. g at g At 9. 81 3600 J
Integrating the remaining part of Equation 3.41 b with respect to "X” yields:
* Z|? - hid? - E k |! = 0
(3.42)
where
U = water column velocity (m/s) g = gravity (m/s 2 ) Z = depth of flow (m) E b = datum to bed elevation (m) h L = head loss (m).
Inserting information from the segments with the highest and lowest velocities, the magnitudes of the remaining terms can be found to be:
U 2 /(2g) = ( 0.88 2 -0.31 2 ) /[ 2 ( 9.81) ] = 0.035 m Z = 3.35 -1.40 = 1.95 m E b = 128 - 122 = 6.0 m h L = 4.015 m (by addition)
These numbers clearly indicate the local acceleration and convective acceleration terms can be neglected without introducing significant error to the problem. Therefore, the simplifications and assumptions required to use a diffusion wave model on the lower Colorado River are valid.
Two of the most notable simplifications concern approximations for the cross-sectional area and the discharge. DAFLOW assumes that these parameters can be estimated using the following equations:
A = A t + A o
(3.43 a
Q = Q DF — s ax
(3.43 b
where
A = area of flow (m 2 ) A o = average cross-sectional area at zero flow (m 2 ) Ai = hydraulic geometry coefficient for area A 2 = hydraulic geometry exponent for area DF = wave dispersion coefficient Q = volumetric rate of flow (m 3 /s) Q s = normal discharge defined as the steady-state discharge that corresponds to an area A (m 3 /s) x = longitudinal distance along channel (m).
Substituting Equation 3.43 b into Equation 3.41 a and simplifying the resulting equation yields:
+ C JA _ df = Q at ax
(3.44)
where
r - Qs 1 ~ Qs2 A 2 - Al aA
DF = S 2 S c W
Inserting Equation 3.43 a into Equation 3.44 and simplifying the derivative terms yields:
. c - DF = 0 at ax a x 2
(3.45)
DAFLOW solves for the space-time variation in discharge represented by Equation 3.45 using a threestep, finite-difference approach. Specific information regarding this process may be found in the program users manual (Jobson, 1989). The main program and the associated support programs are written in Fortran 77 and can easily be run on an IBM-PC or compatible.
Only a couple of minor modifications were made to the DAFLOW code for this project. These changes focused on the output section of the model so that the information generated by DAFLOW could be read into EUTRO-M without many manual revisions.
3.3.2 Calibration of DAFLOW Program
The first step in determining unsteady flow velocities was to calibrate the DAFLOW model. Calibration of the program was performed using hourly gage flow data from the USGS gaging stations in Austin and Bastrop, Texas obtained by Armstrong et al. (1987) for during a previous study of the Colorado River. Initially, estimates of the input parameters were obtained from the wave celerity (travel time) and from typical experimental values given by Jobson (1989). This information was input into the wave celerity model provided by the USGS (CEL.F77) to estimate the hydraulic geometry coefficients (A o , A IZ A 2, W x and W 2) and the wave dispersion coefficient (DF). The time of travel and associated flow rate can be found by determining the time it takes a single peak moving past the USGS gaging station at US Highway 183 to reach the Bastrop gage. This concept is illustrated in Figure 3.3.1. As indicated, it takes 30 hours for an initial flow of
approximately 2900 cfs to travel the 56.9 river miles (91.5 km) between gages.
Flows were input into the model at user specified nodal points. As previously stated, Armstrong et al. (1987) had obtained flow readings from the USGS for the gaging stations located near the upstream end of the study area and at Bastrop, Texas for the 1986 water year. For the calibration phase, the period from March 1, 1986 to March 20, 1986 was used. As shown in Figure 3.3.2, power releases by the LCRA make the discharges quite irregular at the upstream end of the river. Hourly inflow data were not available for the Boggy Creek, Walnut Creek, and Onion Creek tributaries. Instead, the average daily flow from each of these streams was input as tributary inflows. Similarly, discharges from the Govalle, Walnut Creek, and South Austin Region WWTPs were estimated from the annual average flows reported by each of these treatment facilities (Armstrong et al., 1987). These flows were then routed downstream using the DAFLOW model.
Initially, calibration of the model was completed by adjusting the hydraulic geometry coefficients (A x and A 2) and the wave dispersion coefficient (DF) until the mean error and root-mean-square error between the calculated and measured stream flows at Bastrop, Texas were minimized. Based strictly on these criteria, the mean error and the root mean square error were determined to be 2.7 cms and 7.1 cms (95 cfs and 250 cfs), respectively. The statistics were calculated using a slightly modified version of the FLWPLT model
supplied by the USGS. The screen graphics generated by the original FLWPLT model were not compatible with PC FORTRAN compiler and were therefore removed from the program.
In examining the results of the initial calibration phase an interesting and important trend was identified. The computed arrival times of the hydrograph peaks for flows greater than 42.5 cms (1,500 cfs) were consistently shorter than the measured data at Bastrop indicated. Conversely, the computed hydrograph peaks for flows less than 42.5 cms (1,500 cfs) arrived after the measured flows. This phenomenon can be seen clearly in Figure 3.3.3.
The hourly observed gage flows at Bastrop for the verification period (May 1986) and the simulation period (July 1986) indicated a series of daily peak flows greater than 42.5 cms (1,500 cfs). In fact, the average daily flow for the simulation period was slightly over 45.3 cms (1,600 cfs). Using the parameters found in the original calibration phase caused the computed hydrograph peaks in the verification and simulation runs to arrive 3 to 4 hours before the measured hydrograph. With this in mind, the model parameters were adjusted with the primary goal of matching the timing and magnitude of the flows greater than 42.5 cms (1,500 cf s) .
The results of the final calibration run can be seen in Figure 3.3.4. In general, the higher flows predicted by DAFLOW match those measured at the USGS gage in Bastrop very well. However, the timing adjustment required to match the larger hydrograph peaks substantially worsened the timing of the lower flows increasing the mean error to 4.4 cms (155 cfs) and the RMS error to 8.6 cms (303 cfs). All of the calculated low flows arrive at Bastrop before the measured values. Because it was more important to match the simulation flows, the parameters found in the final calibration phase were used throughout the flow modeling. This points out a limitation of the DAFLOW model: it must be calibrated with flows similar to those being used during the simulation run.
There are other discrepancies between the calculated and measured flows which point out a couple of other interesting characteristics of the problem. For example, the model seems to slightly under predict the lag of some of the larger peaks. This can be explained in a couple of different ways. First of all, the accuracy of the gage flows might be questioned. The USGS classifies gage flow records as (1) excellent (within 5 percent of actual flows), (2) good (within 10 percent), (3) fair (within 15 percent), or (4) poor (greater than 15 percent error). A discrepancy at the upstream gaging station of 10-15 percent, the accuracy limit established by the USGS (1986) for the US 183 gaging station, could impact the predicted values at Bastrop.
A second possibility steins from the fact that DAFLOW is unable to account for the changes in flow due to bank storage. As pointed out by Sharp and Larkin
(1988), under high flow conditions the Colorado River can supply water to the surrounding alluvium. The continuity equation embedded in DAFLOW requires that the entire flow be routed downstream rather than delayed through exfiltration into the streambank.
Another problem evident in Figure 3.3.4 occurs during low flow periods. During these periods the model consistently under predicts the discharges measured at Bastrop. Inflows from ungaged tributaries and ground water infiltration were not estimated for the hydrodynamic calibration phase of the project. Ground water studies of the area indicate that substantial quantities of water may enter the river from the alluvium and the Carrizo-Wilcox aquifer (Brune and Duffin, 1983; Follett, 1970; Woodward, 1989; Sharp and Larkin, 1989). Before the DAFLOW modeling results were used in the water quality modeling phase, conservative allowances were made to estimate each of these quantities. The procedures used for allocating the inflows are discussed in subsequent sections of this report.
It must be pointed out that although the magnitude of the total inflow can be determined by examining the flows deficits at the Bastrop gage, the inflow locations are much less certain. The USGS operated gaging stations on a couple the tributaries from the late 1970's and early 1980's. Although currently discontinued, the historic flows indicate that the stream flow in some of the ungaged tributaries can be extremely variable.
3.3.3 Verification of DAFLOW Program
To avoid the potential problems found during the calibration phase, the DAFLOW model was verified using upstream boundary conditions more similar to those which existed during the simulation time frame. A period from May 1, 1986 through May 10, 1986 was used for program verification. Figure 3.3.5 illustrates the cyclical pattern of the summer releases by the LCRA with daily peak flows varying between 85 and 113 cms (3,000 and 4,000 cfs). The exceptions to this pattern occurred due to storm events on May 1 and May 9-10. These events are also illustrated in the figure.
The simulated results of the measured versus calculated flows at Bastrop are generally much better for this simulation than in the calibration run. The timing and magnitudes of the flows between 56.7 cms and 85.0 cms (2,000 and 3,000 cfs) are extremely close reducing the mean error to 4 cms (140) cfs. However, the root mean square error increased to 16.9 cms (600 cfs). The cause of most of these errors is the May 9-10 storm event. Because average daily flows were used at the Boggy Creek, Walnut Creek, and Onion Creek gages, the timing and magnitude of the flood hydrograph is incorrect. As seen in Figure 3.3.6, several of the calculated values between hour 230 and hour 240 are hundreds or thousands of cfs too low. Furthermore, the USGS gage at Onion Creek was inoperative during that period and the discharge estimated by the USGS was classified even by themselves as poor.
Verification of the calibrated DAFLOW model using an independent data set determined that the parameters selected during calibration were representative of the river system. Since the model produced acceptable results under both sets of conditions, it was assumed that the model was valid and that DAFLOW could be used to predict the surface water flows at each surface water segment during the actual study period. The application of this model to the July 1986 data is discussed in Chapter 6.
Figure 3.1.1 Notation Schematic for Finite-Difference Equation (After Ambrose, et al., 1988)
Figure 3.2.1 Kinetics Linkages Between the Ecological Components (After Johns, 1988)
Figure 3.2.2 EUTRO-M State Variable Interactions (After Ambrose, et al., 1988)
Figure 3.2.3 Representation of Multi-Layer Water/Sediment System
Figure 3.2.4 Characteristic Phosphorus and Nitrogen Growth Coefficient Factors for Aquatic Vegetation
Figure 3.2.5 Characteristic Light and Temperature Growth Coefficient Factors for Aquatic Vegetation
Figure 3.3.1 Typical Time of Travel Hydrograph
Figure 3.3.2 Upstream Boundary Condition for DAFLOW Model Calibration Run
Figure 3.3.3 Original Statistical Calibration of DAFLOW Model
Figure 3.3.4 Modified Calibration Using Flows Greater Than 1,500 Cubic Feet per Second
Figure 3.3.5 Upstream Boundary Condition for DAFLOW Model Verification Run
Figure 3.3.6 Verification of DAFLOW Model
CHAPTER 4 INPUT PARAMETER ESTIMATION
4.1 NEED FOR PARAMETER ESTIMATION
Although the modified EUTRO-M model formulation and the other support models can be applied to a number of water quality modeling situations, much of the input data used in the EUTRO-M simulation is site specific. Whenever possible, field measurements on the lower Colorado River have been made to justify the values used in the model. This chapter describes the input parameters used by the various models that were not specifically measured during the field investigations and the methods used for determining the appropriate values.
4.2 WATER COLUMN ADVECTION FOR HYDRODYNAMIC FLOW SIMULATIONS
Operation of the Highland Lakes chain of reservoirs significantly effects the downstream flow rates which makes forecasting future short duration hydrologic conditions (such as daily) flows extremely difficult. For this reason, the unsteady portion of this study was performed using historic flow patterns measured at the USGS gaging station downstream of Longhorn Dam from July 1 to July 31, 1986 with the understanding that these identical conditions will never be exactly duplicated.
4.2.1 Coefficients Required for Dynamic Flow Routin
As stated in Section 3.3, simplicity of input was one of the primary reasons for selecting the DAFLOW model. Very few input parameters were required to use the model. The development of these coefficients was discussed in a previous section. After calibration, the final values of the coefficients for flows in cubic feet per second were:
AO = 0.0 Al = 0.6 A 2 = 1.0 W 1 = 40.80 W 2 = 0.26 DF = 10800.0
A couple of notes should be made about the selection of these parameters. First, the area of dead pools, AO, does not significantly influence the prediction of the water hydrograph. Furthermore, since AO can only be determined using the results of a comprehensive dye study (which was not available during the calibration, verification or simulation study periods), the value was arbitrarily set to zero. The selection of a zero value is not meant to insinuate that there are no pools in the Colorado River. There are, of course, numerous pools within the study reach. In fact, future studies may wish to investigate the significance of such pools. The selection was made simply because the hydrodynamic routing is insensitive to this parameter.
Similarly, as pointed out by Jobson (1990), the DAFLOW model is extremely insensitive to the values selected for W 1 and W 2. Consequently, these values were not changed during the model calibration. The significance of these two width parameters is only evident when modeling temperature or a highly volatile substance with the other models available from the USGS as part of the DAFLOW model package.
4.2.2 Estimation of Unpaged Tributary Inflow
Figure 4.2.1 is a schematic of the study area showing the possible point sources of inflow to the Colorado River drainage basin. Obviously, the difference between the USGS gage flows at Bastrop and at the US 183 bridge represents the net amount of inflow that occurs within the basin once the attenuations of any flood hydrographs have been taken into account by the model. Some of the inflows come from measured or gaged sources and are therefore easy to quantify (within the accuracy of the measurement technique). However, the location and quantity of the other inflows are very difficult to predict. The process is further complicated by the fact that groundwater inflows (or outflows) are occurring along the main stem of the river (Woodward, 1989).
Despite these drawbacks, it was necessary to devise some method for assigning inflows. For the ungaged surface water tributary inflows this was accomplished by calculating the average yield per square mile of the
three gaged drainages over the same time period. Applying the arithmetic average of these yield factors to the other basins allowed for a rough estimate of the other tributary inflows. Table 4.2.1 illustrates the yield factors, tributary areas, and estimated flow rates.
4.2.3 Estimation of Ground Water Inflows
Even with the addition of the ungaged tributary inflows, there remained a discrepancy between the predicted flows and the measured USGS gage flows at Bastrop, Texas. It was assumed these differences were the direct result of ground water inflow. The amounts of these flows are very dependent on season. In March, for example, tens to hundreds of cubic feet per second were unaccounted for in the simulation. However, by July most of the differences had disappeared. While the magnitude of this fluctuation was not known, the general result was expected. Hibbs (1990) discussed the problems of modeling ground water in this area due to seasonal pumping for irrigation and evapotranspiration losses from the alluvium.
As alluded to earlier, while the volume of this flow has well definable values, the locations of the sources/sinks are not easily identifiable. The problem is further complicated because the lower Colorado River flows over a few different geologic formations within the study area. From Austin to just below Webberville the river flows over the relatively impermeable Navarro- Taylor formation. In this reach, ground water interaction is primarily a function of the alluvium. After the river enters Bastrop County, the underlying aquifer changes to the Carrizo-Wilcox. The ground water heads in this aquifer tend to force water into the alluvium surrounding the river (Hibbs, 1990). Therefore, all other things being equal, the ground water interaction in this area would be expected to be somewhat greater than in the upstream reach.
In spite of (or because of) these general geologic observations, definitive conclusions concerning the locations of the ground water flows are not currently available. Consequently, the interactions were distributed uniformly over the river reach between the two main stem gages.
4.3 SELECTION OF DISPERSION COEFFICIENT
There were several different types of dispersive action which had to be considered while developing the proper input conditions including (1) advective diffusion, (2) molecular diffusion, (3) turbulent diffusion, and (4) numerical dispersion. The literature regarding dispersion can be quite confusing. Many authors use different names for describing the same physical processes or the same name for different processes. For example, the term hydrodynamic dispersion is often used in groundwater to define the macroscopic outcome of a combined advective and molecular diffusion term (Bear, 1972). However, Johns
(1988) used hydrodynamic dispersion to represent turbulent and molecular dispersion in a surface water application. While these definitions may seem different at first, further investigation reveals that Bear included both the laminar flow regime and the turbulent flow regime in his definition of convective diffusion. The significance of each type of diffusion and the magnitudes used in the EUTRO-M simulations are discussed in the following sections.
4.3.1 Advective Diffusion
Advective or convective diffusion occurs as the result of variations in the velocity field. For surface water problems this variation can be caused by the noslip condition imposed at stream boundaries. Similarly, variations in the velocity distribution within each pore and changes in permeability may lead to convective diffusion in a ground water segment. Advective diffusion will not be considered as a separate dispersion component in the EUTRO-M program. For diffusion in the water column this will be included as part of the longitudinal dispersion coefficient. For groundwater applications it is included in the pore water diffusion term.
4.3.2 Molecular Diffusion
Even when a volume of fluid is considered to be "at rest", the molecules within the fluid element continue to move randomly with respect to each other. Consequently, a localized high concentration of pollutant will tend to eventually spread out with time. This process is called ordinary molecular diffusion (Li, 1972). The term ordinary is used to distinguish this process from diffusion that occurs as a result of strong temperature or pressure gradients.
The effect of molecular diffusion on the overall dispersion in ground water is more significant at low flow velocities (Bear, 1972). For water quality studies in rivers involving water column transport, molecular diffusion is usually small when compared to turbulent diffusion (Holley and Jirka, 1986; Li, 1972; Rinaldi, et al., 1979). Furthermore, for the time step used in this model, the molecular diffusions in the water column are several orders of magnitude smaller than the numerical dispersion values (see Section 4.3.4). With this in mind, molecular diffusion in the water column was neglected.
Bear (1972) summarizes the work performed by several researchers concerning the relationship between convective dispersion and molecular diffusion in porous media such as the upper and lower sediment layers. These researchers found that the overall coefficient of hydrodynamic dispersion is a function of the Peclet number of molecular diffusion, Pe = Vd/D d = (q/n)d/D d , where d is the mean grain size of the soil layer, n is porosity, and D d is the coefficient of molecular diffusion. Moreover, the conclusion drawn from these studies is that molecular diffusion must be considered for Peclet numbers of diffusion less than five.
4.3.3 Turbulent Diffusion
Turbulent (or eddy) diffusion can be an important component of surface water diffusion because fluid flows in nature are almost always turbulent (Fischer, et al., 1979; Li, 1972). Turbulent diffusion is caused by differing velocities within the numerous eddies that make up the main flow. As stated previously, the surface water component of turbulent diffusion will be incorporated into a longitudinal dispersion coefficient and the groundwater component will be part of the pore water dispersion.
4.3.3.1 Longitudinal Dispersion in the Water Column
For an essentially one-dimensional problem, water column dispersion is commonly referred to as longitudinal dispersion. Longitudinal dispersion includes the combined effects of differential advection associated with the transverse velocity distribution and transverse mixing (Holley and Jirka, 1986). Ambrose et al. (1988) report that longitudinal dispersion may play an important role in the dilution of peak concentrations caused by unsteady loads.
The dispersive exchange between two water column segments, i and j, at time t is given by:
(4.1)
aM - E j ; (t) A; ; °' 1 i _ 1 J 1 J f r p । T <d L L c j t J where
M ± = mass of pollutant in segment i (g) C = total concentration of constituent (mg/L) E i:) (t) = time dependent dispersion coefficient function for segments i-j (m 2 /d) Aij = interfacial area between segments i-j (m 2 ) L cij = characteristic mixing length (m) .
The interfacial area, dispersion coefficient, and characteristic mixing length are specified by the user. The characteristic mixing length represents the length over which mixing occurs and is defined for surface water segments as the distance between adjoining segment midpoints. For interactions between the benthic sediment layer and the water column this distance is typically defined as the thickness of the sediment layer.
Both analytical and empirical attempts have been made to find generalized expressions for the dispersion coefficient (Holley and Jirka, 1986). Johns (1988) used two equations cited by the USEPA (1983) to find the average dispersion coefficients for the lower Colorado River during low flow conditions (100 cfs). Equation 4.2 was developed by Fischer (1975) and Equation 4.3 by Liu (1977). Holley and Jirka (1986) report Fischer's equation was based on measured velocity distributions and transverse mixing studies and Liu's equation on fifteen measured values of the dispersion coefficient (some of the experiments were the same as those used by Fischer). The equations are:
V 2 p 2 E = 0.011 H U*
(4.2)
r e V 2 B s E = ~ U* A
(4.3)
with
g = 0.18 (U^/V) 372
where
V = average velocity (m/s) U* = shear velocity (m/s) = (gHS) 1/2 B = river width (m) H = average depth (m) A = cross-sectional area (m 2 ) g = gravity (m 2 /s) S = channel slope (m/m).
Note: shear velocity assumed equal to (gHS) 1/2 where g is gravity (m/s 2 ), and S is channel slope (m/m).
In a 1987 study of the Colorado River below Longhorn Dam, the Texas Water Commission (1987) used the following equation to estimate longitudinal dispersion:
E = 18.53 n V H°' 833
(4.4)
where
n = Manning's roughness coefficient V = average velocity (m/s) H = average depth (m) .
Using Equation 4.4 with a Manning's n value of 0.035 (the value recommended by the TWC) produces much smaller dispersion coefficients than those predicted by the
equations proposed by Fischer or Liu (on the order of 1000 times smaller). While not explicitly stated by the TWC, their equation for longitudinal dispersion is based on Harleman's adaptation of G. I. Taylor's and J. W. Elder's work (Holley, 1990). Taylor originally developed a dispersion equation for flow in a pipe. Elder modified this equation based on an infinitely wide plane and on the assumption that the vertical velocity profile could be approximated by the von Karman's logarithmic velocity formula (Fischer et al., 1979). By solving Manning's equation for the friction slope, substituting this into the expression for shear velocity, and replacing the shear velocity term in Elder's equation with the resultant, Harleman derived Equation 4.4.
Elder's equation does not apply to flows in natural channels because it bases dispersion only on the vertical logarithmic distribution of velocity instead of taking into account the contribution of the transverse velocity (Holley, 1990; Fischer et al., 1979; USEPA, 1985). It was shown by Fischer (1967) that the transverse profile is at least 100 times more important in producing longitudinal dispersion than the vertical profile for natural channels with width to depth ratios greater than about ten. For width to depth ratios greater than about fifty, it has been found that Fischer's equation tends to over predict dispersion (Holley, 1990).
For the reasons cited above, the longitudinal dispersion coefficients used in this study were based on the average of the Fischer and Liu equations. The values were calculated using a spreadsheet and then transferred to the EUTRO-M input file. The typical magnitudes of the dispersion coefficients are shown in Figure 4.3.1. The figure also illustrates how the coefficients vary with flow rate and distance below Longhorn Dam. Note that the dispersion coefficients vary with distance because changes in the channel characteristics cause different flow velocities, depths, widths and shear velocities.
It should be pointed out that these average dispersion coefficients were used even in the dynamic simulations. In reality, the rapid changes in stream discharge would cause the dispersion coefficients to change each time step. However, the data requirements for such an analysis would be enormous. This was assumed to be especially true in light of the variability exhibited by the dispersion equations and the lack of any measured field data.
4.3.3.2 Pore Water Diffusion
Pore water diffusion, the exchange of dissolved pollutant between a benthic sediment layer and the overlying water column, can significantly influence benthic pollutant concentrations (Ambrose, et al., 1988). Whether this type of diffusion acts as a source or sink of pollutant to the sediment layer depends on the gradient of the dissolved concentration. If the concentration is higher in the water column, then the
material will be diffused to the sediment layer thus resulting in a sink for the water column but a source for the sediment layer.
EUTRO-M uses the following equation to calculate diffusive exchanges:
aM 1 E i ;( t ) A j ; n 1 ; = —J- —J J- (f Cj /nj . f DI Ci / ni ) wrere Ecij /n ij
(4.5)
Mi = mass of pollutant in segment i (g) n i3 = average porosity at interface i-j (L w /L) f Di ,f D j = dissolved fraction of chemical in i and j Eij(t) = time dependent dispersion coefficient function for segments i-j (m 2 /d) A i;J = interfacial area between segments i-j (m 2 ) L cij = characteristic mixing length (m) .
As reported by Thomann and Mueller (1987), a summary of available data by Di Toro et al. (1981) indicated that the resistance to exchange across the sediment/water interface occurs primarily in the sediment layer because the vertical dispersion in the water column is several orders of magnitude greater than that in the sediment. Consequently, the dispersion coefficients used to account for the diffusive exchange are affected by fluid velocity and pore configurations in the sediment (Bouwer, 1978). In fact, studies of porous media problems indicate that longitudinal dispersion coefficients increase exponentially with fluid velocity and with particle size (Harleman et al., 1963; Bruch, 1970).
Harleman et al. (1963) found longitudinal dispersion coefficients in porous media ranging from approximately 0.017 to 0.130 m 2 /day for sands with an average particle size (D 50 of 0.45 to 1.4 mm and fluid velocities from 0.008 to 0.12 cm/s. Similarly, Shamir and Harleman (1967) reported transverse dispersion coefficients between 0.0009 and 0.017 m 2 /day for sands with an average particle size of 0.48 to 1.67 mm for a range of fluid velocities. For sediment/water interfaces, Di Toro et al. (1981) found the dispersion coefficient to be on the magnitude of molecular diffusion (10~ 4 m 2 /day) . Because the work summarized by Di Toro et al. dealt directly with sediment/water interfaces, it was assumed 10' 4 m 2 /day would be a more appropriate estimate of the dispersion coefficient. The sensitivity analysis examines the effect of this estimate.
4.3.4 Numerical Dispersion
Use of finite difference techniques to approximate derivatives produces artificial mixing called numerical dispersion (Ambrose et al., 1988). The problem is created because higher order terms in the Taylor Series expansion are generally neglected when developing the difference equation. Thus the term local truncation error is sometimes used to describe numerical dispersion (Smith, 1985). As pointed out by Holley and Jirka (1986), obtaining a numerical scheme that does not significantly contribute to the apparent system dispersion is not a trivial matter. Selection of
appropriate spacial and time grids and solution technique can effect the magnitude of the dispersion error.
This application of EUTRO-M used a backward difference approximation to represent the change in concentration with respect to distance, x, and a forward difference approximation to represent the change in concentration with respect to time, t. The total numerical dispersion, E nun , created by this solution scheme is given by (Ambrose, et al., 1988):
Enun = (L - V At)
(4.6)
where
V = average stream velocity (m/s) L = length of segment (m) t = time step (s).
Looking at Equation 4.6, it should be evident that by selecting a time step equal to L/V (or Volume/Flow if multiplied by the area) the numerical dispersion goes to zero. Unfortunately, such a selection causes instability in the EUTRO-M model. To maintain stability at a segment, the advected, dispersed and transformed mass must be less than the resident mass. This constraint translates into a restriction on the maximum time step which can be given by the following equation:
Volume; Ai max = Min < — — — ZE Qij *ZE R ij * 2Z ( SjK Volumej/Cj) i i k
(4.7)
where
Q i:j = advective flow between segments i-j (m 3 /day) R i:) = dispersive flow between segments i-j (m 3 /day) Cj = concentration of constituent (g/m 3 ) S jk = kinetic transformation (g/m 3 -day).
It is important to realize that both sets of criteria must be met for every segment in the system. In general, large time steps reduce numerical dispersion and small segments produce more restrictive limits on time step criteria as long as the length component of Equation 4.6 is greater than the velocity component. A time step of 0.001 days satisfied both these conditions for the steady state simulations.
Using Equation 4.6, a time step of 0.001 days, and a flow of 45.3 cms (1600 cfs), the numerical dispersion values in the water column segments ranged from a low of approximately 52 m 2 /s to a high of 244 m 2 /s. These values are several orders of magnitude larger than typical molecular diffusion coefficients. Molecular diffusion coefficients have values in the range of 0.0001 to 0.015 cm 2 /s. As discussed earlier, this is one reason that molecular diffusion in the water column segments was neglected.
It should also be pointed out that by changing an input parameter called the advection factor, ADFAC, the solution technique used in EUTRO-M for the distance derivative term can be changed from a backward difference approximation (ADFAC=O.O) to a central
difference approximation (ADFAC=O.S) or something in between. However, since forward and central difference approximations are unstable in EUTRO-M, an ADFAC value greater than 0.4 is not recommended (Ambrose et al., 1988). Dresnack and Dobbins (1968) explain that this instability is caused by the fact that under pure advection the concentration C i/n+l could not physically be dependent upon the downstream concentration C i+lzn , which occurred in the previous time step. Moreover, since Bella and Grenney (1970) found that errors in the advection equation caused by finite difference solutions are minimized using a backward differencing scheme, an ADFAC of zero was selected. The effect of this decision is discussed in Chapter 5.
In a previous study of the Colorado River, Johns (1988) claimed that the numerical dispersion was approximately equal to longitudinal dispersion. Therefore, by ignoring the numerical dispersion and by setting the longitudinal dispersion to zero, the correct solution was obtained. McCutcheon (1989) reports that this technique is a common approach to the problem. However, Holley and Jirka (1986) argue that this method should be avoided whenever possible.
For unsteady flow simulations the practice employed by Johns is not a practical solution. As was previously shown, longitudinal dispersion varies depending on the magnitude of the flow while numerical dispersion is a function of flow and time step. Since both dispersion terms changed throughout the simulation, there was no guarantee that the terms would remain equal. Therefore an attempt was made in these simulations to include longitudinal dispersion terms and select time steps that would minimize numerical dispersion while maintaining program stability.
4.4 NUTRIENT AND CONSTITUENT LOADINGS
As previously discussed, phosphorus, macrophyte, and sediment concentrations in a river system are influenced by many complex factors. One of the most important factors to quantify is the magnitude of the nutrient loadings. As discussed in Chapter 2, nutrient sources can be derived from a number of external and internal sources. The ones used in the Colorado River study are presented below.
4.4.1 Initial Constituent Concentrations
The EUTRO-M model requires initial conditions at each model segment be specified for each system being modeled. The final output of the model is not extremely sensitive to these initial conditions. For steady state analysis, the closer the initial conditions are to the final concentrations, the faster steady state conditions will be achieved. For unsteady analysis, the initial state variable concentrations will be flushed rather quickly and will therefore have little effect on the final concentrations. The notable exceptions are the macrophyte and epiphytic algae systems. Their growth rates change quite slowly compared to the other systems. Therefore, one would like to begin the simulation with
initial conditions close to the final values. If desired, EUTRO-M will calculate these initial steady state macrophyte and epiphytic algae potentials based on the initial nutrient profiles using a subroutine called SECANT. This option was invoked during these analysis.
As a first approximation, the initial steady state conditions for the surface water nutrient concentrations in EUTRO-M were those concentrations calculated in Johns' (1988) Scenario 1 simulation. However, to shorten the required simulation time and to facilitate closer approximation of initial macrophyte densities, the concentrations were adjusted during the simulations so that the initial values were close to the steady state values.
The initial conditions specified for the phytoplankton system were based on the chlorophyll A samples measured during a previous study of the Colorado River (Armstrong, et al., 1987). The coefficients representing phytoplankton were adjusted within acceptable ranges until the predicted results matched the field data. The concentrations for model segments located between the sampling stations were estimated using linear interpolation.
Each dynamic simulation scenario used the steady state concentrations described above as initial conditions. However, because of dilution errors near the point sources, several of the values were adjusted manually.
For the upper sediment layer, the initial nutrient conditions for the calibration run were estimated from measured field data. For the alternative loading scenarios, the nitrogen and phosphorus concentrations were estimated using the Freundlich isotherm. It is commonly assumed that the correlation between water column concentrations and bed sediments can be determined using the following relationship:
kl * ( C water )
(4.8)
where
C water = concentration of nutrient in water (mg/L) = concentration of nutrient in sediment, measured on a total bed mass basis (mg/kg) XI = empirical constant ({mg/kg}/{mg/L}) X 2 = empirical constant.
Based on field measurements of water and sediment concentrations on the Lower Colorado River, Johns (1988) found the values of these coefficients to be:
XI X 2 Ammonia 152. 1.000 Nitrate 152. 1.000 Phosphorus 104. 0.822
Unfortunately there is a flaw in Johns' original hypothesis which allowed him to relate the sediment concentrations of nitrate to the water column concentration via the use of a Freundlich isotherm. The nitrate coefficients were based on field data which
measured only Total Kjeldahl Nitrogen in the sediment layer. In developing the coefficients for nitrate, Johns neglected to account for the kinetic transformation which occurs as nitrate moves from an aerobic water column to an anaerobic sediment layer.
As discussed by Brezonik (1977), there are three possible sources of nitrate to the sediments: (1) diffusion of nitrate from the water column, (2) seepage of nitrate bearing groundwater through the sediment, and (3) nitrification of ammonia in the oxidized sediment layer. However, once nitrate enters the reduced sediment layer, it is quickly converted to either molecular, ammonia or organic nitrogen. The reasons for this have been discussed in the literature published by various researchers.
First, since Nitrosomonas and Nitrobacter are obligate aerobes, nitrification of ammonia to nitrate will only occur in aerobic sediments (typically the top few millimeters). In fact, studies indicate that nitrification will not occur when dissolved oxygen levels are less than about 0.3 mg/L (Keeney, 1973).
Second, Kamp-Nielsen and Andersen (1977) reported that particulate, organic nitrogen accounts for the majority of nitrogen present in sediments. In his analyses of 13 sediment samples, Keeney (1973) found that organic nitrogen made up 82.5 percent of the total nitrogen. And, while interstitial water in sediments usually consist of much higher soluble ammonia and organic nitrogen concentrations than the overlying water column, Konrad et al. (1970) reported in their study of sediments, that the nitrate concentration were usually very low.
Finally, although ground water seepage accounts for a large portion of the flow in high runoff months and could therefore transport a lot of nitrate, very little net inflow is occurring during July. A number of researchers have found that denitrification, the conversion of nitrate to molecular nitrogen, occurs very rapidly near the water-sediment interface. Edwards and Rolley (1965) found that for river sediments, nitrate could be consumed at 120-1440 mg N/m 2 -day although this rate is thought to be highly dependant on pH. Using lake sediments, Chen et al. (1972) found that up to 90 % of the added nitrate disappeared from the system within 2 hours.
As a result of these findings, the nitrate and ammonia coefficients were recalculated. Using a conservative estimate (based on Keeney's 1973 study), 85 percent of the total nitrogen in the sediment was assumed to be organic nitrogen. The remainder was assumed to be either ammonia or nitrate. Since macrophytes can use both forms of nitrogen the exact percentages were not considered to be important. However, to be consistent with the literature, only a small amount of nitrate was assumed to be present in the sediment layer due to the thin aerobic portion at the water/sediment interface. Therefore, 99 percent of the remaining nitrogen was assumed to be ammonia.
Inserting these values into the Freundlich isotherm along with the water column concentrations, new coefficients for ammonia and nitrate were determined. These coefficients were used in subsequent simulations to determine the initial sediment concentrations.
4.4.2 Boundary Conditions
In EUTRO-M, every interface where water enters or leaves the system requires a boundary condition for each model constituent. This includes the headwater segment, tributary inflow segments, point source loading segments, ground water inflow segments and the downstream boundary condition. EUTRO-M allows these conditions to be specified as a function of time if the user so desires. The following two subsections describe the boundary conditions used in the model and how these values were determined.
4.4.2.1 Under Steady State Flow Conditions
Water column concentrations were measured at selected stations during a previous study of the Colorado River (Armstrong et al., 1987). The calibration phase of the model used the conditions measured during a July 8-13, 1986 sampling trip. Wherever appropriate, these values were used as boundary conditions. However, these samples were not actually taken during steady state flow conditions so interpretation of the values was sometimes necessary.
All of the tributary drainages within the study area were treated as a type of point source load. Each point of inflow was treated as a boundary condition and assigned a discharge and a concentration for each water quality parameter. In the case of Walnut Creek, this required averaging the concentrations due to the waste water treatment plant discharge and the natural flow of the creek. This was accomplished using a weighted average based on percentage of flow.
Boundary conditions at the ground water inflow interfaces for nitrate, ammonia, and phosphorus were estimated from typical values reported in ground water wells in Travis and Bastrop Counties (Follett, 1970; Baker et al., 1986). Except for phosphorus, these values were found to be highly variable. Therefore, further research in the area of alluvial water quality might be necessary in the future.
4.4.2.2 Under Unsteady Flow Conditions
Under typical unsteady flow conditions, nutrient concentrations and other water quality parameters would be expected to vary quite significantly as the flow at the upstream boundary increased from 50 cfs to 3500 cfs (a typical range for the July 1986 flows). However, in this case, the flows originate from Town Lake and are the result of scheduled releases rather than storm events. As a result, the upstream boundary conditions were assumed to be constants. This assumptions appears to be justified since Miertschin and Armstrong (1986)
reported adequate results when they modeled Town Lake as a continually stirred tank reactors (CSTR) thus indicating that the concentrations in the lake are relatively uniform. Furthermore, the July 1986 precipitation records for the Austin area indicated the only rainfall event of the month came on July 7 when the weather station in Austin recorded 0.45 inches of precipitation. The lack of a substantial number of rainfall events reduces the likelihood of highly variable runoff loads occurring from the surrounding drainage basins.
The boundary conditions at the ground water inflow interfaces were assumed to be constant. Therefore, the boundary conditions used in the steady state analysis were specified for the hydrodynamic simulations. This may not be an accurate assumption given that ground water inflows in this area are highly variable and often tied directly to precipitation events and surface water discharge. However, the lack of field measurements made this type of hypothesis necessary.
4.4.3 External Sources
4.4.3.1 Historical Loadings
Point source loads within the lower Colorado River study area consisted primarily of Austin's waste water treatment plant discharges. During the July 1986 field investigation period, these discharges came from the Govalle, Walnut Creek, Hornsby Bend, and South Austin
Regional waste water treatment facilities. Historic records of these flows and the associated concentrations were obtained from the Texas Water Commission and used as model inputs during the calibration phase (Texas Water Commission, 1987 b).
The original studies estimated nonpoint loadings using methodologies adapted from the National Eutrophication Survey done by the USEPA in the early 1970's (Armstrong et al., 1987). However, more accurate methods of approximating loads are now available including the Rational method, the Mockus (SCS) method, the unbiased stratified ratio estimator, and the Espey variable flow-variable concentration method (Armstrong, 1989 b). Mean average values for the loadings generated using the SCS method have been calculated and calibrated in previous studies. These mean values were used to determine nonpoint source tributary loadings to the system.
4.4.3.2 Projected Loadings
During the evaluation of the alternative scenarios (See Section 6.3), projected loadings from Austin's waste water treatment facilities were used as model input. These flows and concentrations were calculated in a previous study (Armstrong et al., 1987) using factors such as population estimates made by the Texas Department of Water Resources, historic loading per capita information, permitted plant capacities and anticipated plant operations.
4.4.4 Internal Sources
Methodologies developed in the original study (Johns, 1988) for modeling the internal sources and sinks of nitrogen and phosphorus, caused primarily by the growth and senescence of macrophytes and epiphytic algae and the settling and resuspension of particulate nutrients, were combined with the scour equations developed by Vecchio (1989) and incorporated into the WASP 4 model. The internal sources are modeled after the external sources because they are greatly influenced by the external sources and the transport of nutrients within the system. Actually, the origin of most internal sources can be traced back to external sources. However, the long time delay of their return to the system often makes it simpler to treat them as internal sources.
4.5 SETTLING AND RESUSPENSION VELOCITIES OF SOLIDS SYSTEMS
EUTRO-M permits up to three solids transport fields which allow for the vertical movement of the particulate matter and the adsorbed constituents. These fields can represent any type of division between solids. For example, if the appropriate information were available, one could model three specific grain sizes of inorganic sediment. The solids fields were modeled in this study as organic, phytoplankton, and inorganic particles.
Whether net sediment scour or settling actually occurs is a function of the particle size distribution, cohesiveness, consolidation of surficial benthic deposits, and the water column velocity (shear stress) (Graf, 1971; Ambrose et al., 1988). These criteria require a knowledge of the stream velocities prior to the final water quality simulation. Therefore, a spread sheet was set-up to calculate stream velocities for a number of discharges. Then, using Figure 4.5.1, an initial determination was made concerning whether erosion, sedimentation, or transportation of inorganic particles was occurring. However, as pointed out by Ambrose, et al. (1988), site specific calibration of these criteria is always necessary.
Waste water discharges can contribute organics to the river system. Settling and scour criteria for organic particles are slightly different than those for inorganics. Studies have shown that scour of organic particles occurs when the stream velocity reaches about 0.3 to 0.5 m/s (1.0 to 1.5 ft/sec) and deposition occurs when the stream velocity is less than approximately 0.2 to 0.3 m/s (0.6 to 1.0 ft/sec) (Imhoff and Fair, 1956; Krenkel and Novotny, 1980). These same criteria were used to determine the transport status of phytoplankton even though phytoplankton has a much smaller specific gravity.
The criteria described in the previous paragraphs indicate that erosion of inorganic particles with an average grain size diameter of 0.004 mm will not occur until the mean stream velocity exceeds 1 m/s and that sedimentation of these particles will not occur at any velocity. Erosion would exist only in a couple of
segments when the discharge exceeds 115 cms (4100 cfs). Furthermore, for the overall average segment velocity to exceed 1 m/s, the discharge must be greater than 220 cms (7750 cfs). The result of this being that neither net settling nor net resuspension of particles would occur given the range of flows used in this study. This is contradictory to the conditions which actually exist on the lower Colorado River. Therefore, allowances had to be made in order to reflect the results of previous field investigations.
For organic particles and phytoplankton, the above criteria require a discharge of 35.4 cms (1250 cfs) to produce velocities sufficient for net scour. For net settling to occur, the stream discharge must be less than 14 cms (495 cfs). For the discharges used in this study, only net scour and transport of particles can occur using this criteria.
Once the status of the particles has been determined, settling and resuspension velocities have to be designated as functions of time for each of the particle classes. Remember, like all systems within EUTRO-M, sediment fluxes are based on mass. Therefore, if net resuspension is occurring, a combination of the resuspension velocity multiplied by the sediment concentration must be greater than the settling velocity multiplied by the water column sediment concentration. If, on the other hand, the stream velocity permits net settling over a certain time interval, then a settling velocity which produces a greater mass transfer to the sediment layer is input into the model.
The settling rate for inorganic particles was determined using Stokes' Law which can be written as
8.64 g 2 v s = iBp P P Pw) d P
(4.9)
where
V 6 = settling velocity (m/day) g = gravity (981 cm/s 2 ) d p = particle diameter (mm) p = absolute viscosity of water (0.01 g/cm 3 -s) p p = density of particle (g/cm 3 ) p w = density of water (g/cm 3 ).
Particle density was assumed to be 2.65 g/cm 3 . The mean particle diameter was determined from a 1983 USGS analysis of suspended sediment performed on the Colorado River at the Columbus, Texas gaging station. Figure 4.5.2 illustrates the results of their study. As shown, the mean particle diameter is 0.004 mm. This diameter would classify the average particle as a medium to fine silt.
Settling velocities for organic particles and phytoplankton were taken from data obtained by Thomann and Mueller (1987). In general, these particles are less dense than the inorganic particles and thus have much smaller rates of settling. The values used for each of the settling velocities are:
Flow Field Solids Type Settling Velocity 3 Organic 0.20 m/day
4 Phytoplankton 0.15 m/day 5 Inorganic 1.25 m/day
The velocity given for inorganic settling is based on the assumption that Stokes' Law is valid for the range of particle sizes found in the lower Colorado River. As pointed out by Hawley (1982), settling velocities based on Stokes' Law for particles with diameters less than about 0.1 millimeters are generally very low compared to observed values. Therefore, in the future, further analysis of local conditions may be needed to verify the rate.
A review of the literature revealed that a number of theoretical constructs exist for determining the rates of bottom scour and stream bank erosion (Chang, 1982; Jain and Park, 1989; Odgaard, 1987). For example, Krenkel and Novotny (1980) state that the mass transfer of material scoured (resuspended) from the sediment layer to the water column can be estimated by knowing the shear stress acting on the sediment layer. The method assumes the bed shear stress must exceed a critical shear stress for the incipient motion of the soil particles. Once the critical shear has been exceeded, the resuspension rate varies linearly with the bottom shear stress. Figure 4.5.3 illustrates this relationship. Using this figure yields a resuspension velocity of 6.8 E-06 m/day.
Most, if not all, of the scour models are not capable of accurately predicting the complex combination of channel degradation, stream bank erosion, and
armoring without the benefit of site calibration. Moreover, even with such data, many equations provide only a rough approximation (Schnoor, et al., 1987). In the absence of specific information, Thomann and Mueller (1987) suggest a method proposed by O'Connor be used to predict resuspension velocities. This method asserts that the resuspension velocity can be found using the following equation:
V s Wi(x) - exp f- Vu = 7 : ."V ■ pp(l-n 2 )10 6 1 - exp
(4.10)
where
V £ = settling velocity (m/day) V u = resuspension velocity (m/day) m x = solids concentration (mg/L) n 2 = porosity H = depth (m) U = velocity (m/day) X = distance downstream from origin {point 0} (m) p p = density of particle (g/cm 3 ).
Using Equation 4.10 with average values for stream depth and velocity produces an estimate of 3.9 E-05 m/day for the resuspension velocity.
A final check on the resuspension velocities was performed using measured field data from the Colorado River. The resuspension velocities were adjusted until the model concentrations matched the sample data. This produced scour rates ranging from approximately 6.4 E-06 m/day to 2.7 E-05 m/day. It is important to remember the model predicted steady state concentration while the
measured values represent field conditions which were far for steady state. Nevertheless, the range of values appeared to be consistent with the theoretical values which Schnoor et al. (1987) reported to be between 2.74 E-6 and 8.20 E-5 m/day. Because of this, an average value of 1.67 E-05 m/day was used in the EUTRO-M model.
The settling or scour of material contributes to the mass transfer of each water quality constituent across the sediment-water interface. Each state variable in EUTRO-M can be transported with only one of the solids fields. The solids field and partitioning coefficient used in the model for each variable are presented in Table 4.5.1.
4.6 COEFFICIENTS AND CONSTANTS FOR EPIPHYTIC ALGAE, MACROPHYTE KINETICS, AND WATER QUALITY VARIABLES
The equations used to predict macrophyte and epiphytic algae growth contain a number of parameters which must be determined by the user prior to execution of the program. Some of these constants were computed from measured values, others were taken from the literature, and still others were adjusted during model calibration. These values are entered in Data Group H of the model's main input data file. Table 4.6.1 contains the values of all the program constants used in these analysis. Additional columns in Table 4.6.1 indicate which constants were adjusted during calibration and, where suitable information was available in the literature, the typical ranges of the constants. It is important to point out that in many
cases researchers have found situations where the measured or calibrated values are several multiples higher or lower than the specified "typical" ranges.
Many of the constants shown in Table 4.6.1 are those which were included in the original USEPA model. Table 4.6.2 defines only the new coefficients added to the EUTRO-M program. The reader is directed to the WASP 4 manual (Ambrose, et al., 1988) for a complete description of the other program constants. It must be pointed out that many of the constants have default values which are set within the program if no value is specified by the user. Therefore, just because a constant is not presented in Table 4.6.1, one should not assume a value of zero.
Figure 4.2.1 Schematic of Point Source Inflows within the Study Area
Figure 4.3.1 Variation in Longitudinal Dispersion Coefficients Below Longhorn Dam
Figure 4.5.1 Relationship Between Stream Velocity, Particle Size, and Sediment Regimes
Figure 4.5.2 Particle Size Analysis of Suspended Sediment in Colorado River
Figure 4.5.3 Resuspension Rates as a Function of Bottom Shear Stress
Tributary Drainage Area fmi 2 ) Gaged Flows 1 (cfs) Yield Factor (cf s/mi 2 ) Estimated Flow (cf s) 1. Boggy Cr. 13.1 0.42 0.0321 n.a. 2. Walnut Cr. 51.3 6.43 0.1253 n.a. 3. Onion Cr. 321.0 21.00 0.0654 n.a. sum 0.2228 average 0.0743 4. Carson Cr. 17.1 n.a. 0.0743 1.27 5. Gilleland Cr 74.0 n.a. 0.0743 5.50 6. Dry Cr. 53.0 n.a. 0.0743 3.94 7. Wilbarger Cr 183.0 n.a. 0.0743 13.59 8. Big Sandy Cr 114.0 n.a. 0.0743 8.47 9. Piney Cr. 38.3 n.a. 0.0743 2.84 n.a. = not applicable 1 Average Daily Flows for March 1-20, 1986
Table 4.2.1 Estimated Inflows from Ungaged Drainages
State Variable 59 Associated lids System Dissolved Fraction Particulate Fraction 1 . Ammonia Nitrogen (WC) 3 1.0 0.0 Ammonia Nitrogen (S) 3 0.5 0.5 2 Nitrate Nitrogen 1.0 0.0 J Inorganic Phosphorus 5 0.5 0.5 4 Phytoplankton Carbon 4 0.0 1.0 c Carbonaceous BOD 3 0.3 0.7 6 Dissolved Oxygen 5 1.0 0.0 7 . Organic Nitrogen (WC) 3 0.9 0.1 Organic Nitrogen (S) 3 0.4 0.6 E Organic Phosphorus - 0.2 0.8 WC - water column S - sediment layer
Table 4.5.1 State Variables and Their Associated Solids Systems
Constant Value Adjusted During Calibration Typical Range BI10C 0.400 No o i o U1 BI11C 0.400 No 0.4-0.5 BI12C 0.400 No 0.4-0.5 BI13C 0.400 No 0.4-0.5 BI10N 0.150 No 0.07-0.10 BI11N 0.100 No 0.07-0.10 BI12N 0.100 No 0.07-0.10 BI13N 0.100 No 0.07-0.10 BI10P 0.010 No 0.01-0.012 BI11P 0.010 No 0.01-0.012 BI12P 0.010 No 0.01-0.012 BI13P 0.010 No 0.01-0.012 CCHL 30. No 20-37 FON 0.850 Yes 0-1 FOP 0.850 Yes 0-1 GAMLO 0.820 No GB10 0.000 No GB11 0.000 No GB12 0.000 No GB13 0.000 No GCF1 0.052 No 0.046-0.057 1/m GCF2 0.031 No 0.03-0.04 1/m GCF3 0.174 No 0.170-0.182 1/m GCF10 0.200 Yes GCF11 0.030 Yes GCF12 0.055 Yes GCF13 0.040 Yes INCLT 492. No 200-750 ly/d IS1 350. No 173-518 ly/d K1C 1.680 Yes 1.0-2.5 1/d KID 0.175 Yes 0.003-0.17 1/d K1RC 0.125 No 0.05-0.20 1/d K1RT 1.045 No KIT 1.050 Yes K1013C 0.020 Yes 0.02-0.024 1/d Table 4.6.1 Program Constants Used in EUTRO-M Steady State Analysis Constant Value (Continued) Adjusted During Calibration Typical Range K1013T 1.020 Yes 1-1.03 K1320C 0.030 Yes 0.025-0.20 1/d K1320T 1.030 Yes 1-1.08 K140C 0.100 Yes 0.02-0.10 1/d K140T 1.080 No 1.02-1.09 K58C 0.220 Yes 0.02-0.70 1/d K58T 1.080 No 1.02-1.09 KBOD 0.400 No KDC 0.200 No 0.1-0.5 1/d KDSC 4.230 Yes KDST 1.080 No KDT 1.050 No 1.02-1.09 KMNG1 0.125 No KMPG1 0.001 No KMPHYT 0.000 No KNIT 2.000 No KNO3 0.100 No KONDC 0.0004 No 0 .0004-0.10 1/d KONDT 1.080 No 1.02-1.09 KOPDC 0.0004 No KOPDT 1.080 No KPZDC 0.020 No KPZDT 1.080 No LGHTSW 1.000 No MMLHT 30. No MMN10 0.020 Yes 0.01-0.4 mg/L MMN11 5.100 Yes 3-10 mg/L MMN12 5.100 Yes 3-10 mg/L MMN13 5.100 Yes 3-10 mg/L MMP10 0.005 Yes 0.001-0.08 mg/L MMP11 1.335 Yes 0.3-2.4 mg/L MMP12 1.335 Yes 0.3-2.4 mg/L MMP13 1.335 Yes 0.3-2.4 mg/L MRTAL 0.050 Yes 0.0-0.8 1/d MRTMC 0.100 Yes NCRB 0.250 No 0.17-0.43 OCRB 2.670 No Table 4.6.1 Program Constants Used in EUTRO-M Steady State Analysis (Continued) Constant Value Adjusted During Calibration Typical Range PCRB 0.025 No 0 . 024-0.05 PHIMAX 720. No PHOTO 0.577 No PREFI 0.100 No 0-1 PREF2 0.880 No 0-1 RESPK 0.050 Yes 0.05-0.15 1/d RMIN 0.015 Yes TK1 0.250 No TK2 0.990 No 0.98-1.0 TK3 0.990 No 0.98-1.0 TK4 0.250 No TMAX 34. No 32-40 °C TMIN 10. No 0-10 °C TOPT1 25. No TOPT2 30. No UMX10 0.330 Yes 0.2-1.5 1/d UMX11 0.770 Yes 0.2-1.5 1/d UMX12 0.770 Yes 0.2-1.5 1/d UMX13 0.770 Yes 0.2-1.5 1/d VS 0.305 No 0.05-2.0 m/d VSC10 9600. No VSC11 9600. No VSC12 9600 . No VSC13 9600. No VTHR10 1.000 No VTHR11 1.580 No VTHR12 1.320 No VTHR13 1.450 No VPFC10 0.065 No VPFC11 0.065 No VPFC12 0.065 No VPFC13 0.065 No XKC 0.017 Yes
Table 4.6.1 Program Constants Used in EUTRO-M Steady State Analysis
Constant Value Description BI10C 0.400 Carbon/Total Cell Mass Fraction for System 10 BI11C 0.400 Carbon/Total Cell Mass Fraction for System 11 BI12C 0.400 Carbon/Total Cell Mass Fraction for System 12 BI13C 0.400 Carbon/Total Cell Mass Fraction for System 13 BI10N 0.150 Nitrogen/Total Cell Mass Fraction for System 10 BI11N 0.100 Nitrogen/Total Cell Mass Fraction for System 11 BI12N 0.100 Nitrogen/Total Cell Mass Fraction for System 12 BI13N 0.100 Nitrogen/Total Cell Mass Fraction for System 13 BI10P 0.010 Phosphorus/Total Cell Mass Fraction for System 10 BI11P 0.010 Phosphorus/Total Cell Mass Fraction for System 11 BI12P 0.010 Phosphorus/Total Cell Mass Fraction for System 12 BI13P 0.010 Phosphorus/Total Cell Mass Fraction for System 13 EPSIL1 0.002 Recompute Steady State Plant Growth at this Time from End of Run (days) GAMLO 0.820 Pure Water (Base) Light Extinction Coefficient (1/m) GB1O 0.000 Grazing Rate Loss for System 10 (1/d) GB11 0.000 Grazing Rate Loss for System 11 (1/d) GB12 0.000 Grazing Rate Loss for System 12 (1/d) GB13 0.000 Grazing Rate Loss for System 13 (1/d) GCF1 0.052 Light Extinction Coefficient for Suspended Sediment (L/mg-m) GCF2 0.031 Light Extinction Coefficient for Phytoplankton (L/mg-m) GCF3 0.174 Light Extinction Coefficient for CBOD (L/mg-m) Table 4 .6.2 Definition of Steady State Macrophyte and Epiphytic Algae Constants (Continued) Constant Value Description GCF10 0.200 Self-Shading Light Extinction Coefficient for System 10 (m 2 /gm-m) GCF11 0.030 Self-Shading Light Extinction Coefficient for System 11 (m 2 /gm-m) GCF12 0.055 Self-Shading Light Extinction Coefficient for System 12 (m 2 /gm-m) GCF13 0.055 Self-Shading Light Extinction Coefficient for System 13 (m 2 /gm-m) INCLT 492. Surface Light Intensity (ly/d) MMLHT 30. Michaelis-Menten Light Constant for Systems 10-13 (ly/d) MMN10 0.020 Michaelis-Menten Inorganic Nitrogen Constant for System 10 (mg-N/L) MMN11 5.100 Michaelis-Menten Inorganic Nitrogen Constant for System 11 (mg-N/L) MMN12 5.100 Michaelis-Menten Inorganic Nitrogen Constant for System 12 (mg-N/L) MMN13 5.100 Michaelis-Menten Inorganic Nitrogen Constant for System 13 (mg-N/L) MMP10 0.005 Michaelis-Menten Orthophosphorus Constant for System 10 (mg-P/L) MMP11 1.335 Michaelis-Menten Orthophosphorus Constant for System 11 (mg-P/L) MMP12 1.335 Michaelis-Menten Orthophosphorus Constant for System 12 (mg-P/L) MMP13 1.335 Michaelis-Menten Orthophosphorus Constant for System 13 (mg-P/L) MRTAL 0.050 Maximum Mortality Rate for Epiphytic Algae (1/d) MRTMC 0.100 Maximum Mortality Rate for Rooted Macrophytes (1/d) PHOTO 0.577 Photo Period for Plant Growth (fraction of day) PREFI 0.100 Epiphytic Algae Ammonia Preference Factor (0.-1.;l=prefer NH3, 0=N03) PREF2 0.880 Macrophyte Ammonia Preference Factor (0.-1.;l=prefer NH3, 0=N03) RESPK 0.050 Respiration Increase Under Optimal Conditions (1/d) Table 4.6.2 Definition of Steady State Macrophyte and Epiphytic Algae Constants (Continued) Constant Value Description RMIN 0.015 Minimum Base Respiration Rate (1/d) SECSW 1.000 Steady State Macrophyte Flag SEDMET 0.000 Method of Sediment Analysis TK1 0.250 Growth Rate Coefficient near TMIN TK2 0.990 Growth Rate Coefficient near TOPT1 TK3 0.990 Growth Rate Coefficient near TOPT2 TK4 0.250 Growth Rate Coefficient near TMAX TMAX 34. Maximum Temperature for Aquatic Plant Growth (°C) TMIN 10. Minimum Temperature for Aquatic Plant Growth (°C) T0PT1 25. Lower Optimal Temperature for Plant Growth (°C) TOPT2 30. Upper Optimal Temperature for Plant Growth (°C) UMX10 0.330 Maximum Growth Rate at 20 Degrees C for System 10 (1/day) UMX11 0.770 Maximum Growth Rate at 20 Degrees C for System 11 (1/day) UMX12 0.770 Maximum Growth Rate at 20 Degrees C for System 12 (1/day) UMX13 0.770 Maximum Growth Rate at 20 Degrees C for System 13 (1/day) VS 0.305 Setting Velocity of Suspended Sediment (m/d) VSC10 9600. Scour Coefficient for System 10 (sec/m-d) VSC11 9600. Scour Coefficient for System 11 (sec/m-d) VSC12 9600. Scour Coefficient for System 12 (sec/m-d) VSC13 9600. Scour Coefficient for System 13 (sec/m-d) VTHR10 1.000 Threshold Velocity for Scouring System 10 (m/s) VTHR11 1.580 Threshold Velocity for Scouring System 11 (m/s) VTHR12 1.320 Threshold Velocity for Scouring System 12 (m/s) Table 4.6.2 Definition of Steady State Macrophyte and Epiphytic Algae Constants (Continued) Constant Value Description VTHR13 1.450 Threshold Velocity for Scouring System 13 (m/s) VPFC10 0.065 Michaelis-Menten Protection Factor Coefficient for System 10 (gm/m 2 ) VPFC11 0.065 Michaelis-Menten Protection Factor Coefficient for System 11 (gm/m 2 ) VPFC12 0.065 Michaelis-Menten Protection Factor Coefficient for System 12 (gm/m 2 ) VPFC13 0.065 Michaelis-Menten Protection Factor Coefficient for System 13 (gm/m 2 )
Table 4.6.2 Definition of Steady State Macrophyte and Epiphytic Algae Constants
CHAPTER 5 EXAMPLE PROBLEM
5.1 DESCRIPTION OF PHYSICAL SYSTEM
Because of the size of the actual data set used to model the Colorado River study area (6900 lines, 300kb), it would not be practical to include a detailed description of each input variable nor would it be practical to include a complete listing of the input file. Nevertheless, the complexity of the EUTRO-M water quality model requires that some information be provided so that others may be able to use the model in the future. This goal is accomplished through the use of a hypothetical example problem.
A simple problem was created that would best illustrate the capabilities and options used in the actual simulations. The problem consists of a main flow, a tributary flow, a ground water inflow, longitudinal water column dispersion, boundary conditions, program constants, and initial conditions. Figure 5.1.1 illustrates the layout of the physical system and the numbering scheme used to describe each of the segments. Segment numbers 1-4 are referred to as upper surface water segments, segment numbers 5-8 as the upper sediment layer, and segment numbers 9-12 as the lower sediment layer. The input requirements for such a system are discussed in the next section.
5.2 INPUT REQUIREMENTS FOR EXAMPLE PROBLEM
The input requirements for EUTRO-M are separated into ten lettered groups (A through I). Figure 5.2.1 illustrates a general schematic of the information groups required. Appendix C contains a complete listing of the input data set for the example problem described above. Each group is briefly described below in relation to the example problem. For specific information, the reader is referred to the WASP 4 User Manual (Ambrose, et al., 1988).
Group A consists of data set identification and system control parameters. The first two lines of the input file contain title cards which can be used to describe the data set after which system control data such as the type of simulation, number of segments, number of state variables, print time, and time step must be specified by the user.
Data Group B contains dispersion coefficients for longitudinal water column dispersion and for molecular diffusion. The example problem illustrates each of these cases. The user must first specify the type of exchange, the number of different coefficients which will be used, and the number of segment that will use a given dispersion coefficient. Then, the surface area, mixing length, and the two segments between which the diffusion is occurring is entered. And finally, the time dependant diffusion coefficient is specified.
Group C contains the segment volume, segment type, and channel geometry coefficients for each segment. The geometry coefficients are set-up so that the average depth and water velocity in each segment can be determined by knowing the total discharge passing through the segment.
The flow information is contained in Group D. A main surface water discharge of 25 cms flows into segment number 1 via an upstream boundary condition before proceeding downstream through segment numbers 2 and 3 and finally out of the system through segment number 4. A tributary inflow of 2 cms is assumed to enter segment number 3 and also exit the system through segment number 4.
In addition to the surface water flows, a ground water flow enters the system through segments 9, 10, 11 and 12. It is important to note how this is handled in EUTRO-M. First, the flow is routed across the boundary condition into the lower sediment segment, then from the lower sediment segment to the upper sediment layer, and finally from the upper sediment layer to the surface water segment. Once the flow enters the surface water segment, another flow path is required to route the water downstream. EUTRO-M does not route the flow automatically. Without the companion flow, the surface water segments would continually increase in volume with time. Also, one cannot simply route the groundwater flow downstream. The flow path has to be divided into two flows because interstitial waters transferred as ground water only transport the dissolved fractions of the pollutants.
The next block of data, Group E, contains boundary conditions. Boundary conditions in EUTRO-M are required anywhere that a flow crosses an external boundary such as the upstream headwater segment, the outflow segment, inflow from a point source, or groundwater inflows. The boundary conditions may vary with time if so desired, however, for the example problem, all of the concentrations are held constant.
In some modeling situations, it may be desirable to specify additional constituent waste loading from point or nonpoint sources without any additional flow. Data Group F handles these waste loads. This option was not used in the EUTRO-M model study so each of the loads have been specified as zero.
Environmental parameters are specified in Group G. These parameters allow the modeler to input segment specific data parameters such as salinity, temperature, nutrient flux multipliers, time function pointers, sediment oxygen demands, and suspended sediment properties. This section was modified in EUTRO-M so that additional properties such as sediment resuspension coefficients could be read into the model.
Group H permits the user to specify a number of program constants. This section was modified so that additional constants regarding the new state variables could be specified. EUTRO-M contains numerous default values for most of these constants. Users should refer to the users manual to make sure these default values are appropriate for the particular ecological system
being modeled.
Group I contains time functions for various environmental and kinetic equations. Seasonal or monthly variations in certain parameters requires that the model be flexible enough to allow for changes. Data Group I contains multipliers or other time dependant parameters which permit these changes.
The final data group, Group J, allows the user to specify initial conditions. Initial concentrations for each of the thirteen state variables are required for every segment. The user specifies the segment, the concentration, and the dissolved fraction for each entry. As shown in the example input file, the information is input with each line containing up to three segments.
5.3 DISCUSSION OF RESULTS
Four separate output files are created during the simulation. Each has the same root name as the input file but different extensions. The extensions are ” .DMP”, ’’.OUT”, ”.TRN", and ’’.MSB”. The ”.DMP" file contains forty-seven display variables at each of the print times specified in Data Group A. A list of these variables is presented in Table 5.3.1. EUTRO-M writes information to this file using a sequential unformatted file type. Two advantages are derived from this type of file. By writing the data in binary form, the I/O is faster and the resulting file is considerably smaller
than that produced using standard ASCII text. In order to process this information, the display program provided by the USEPA was slightly modified in order to read a sequential file type.
The ".OUT” file contains an printout of the input file in tabular format which makes checking the input parameters easier along with any simulation error messages which may have been generated during the simulation. The ”.TRN" file contains transport associated information including factors such as the segment volume, discharge, dispersion coefficient and numerical dispersion coefficient. The ".MSB" file contains mass balance information at each print time for whatever state variable is specified in Data Group A. This includes information such as advection, dispersion, and segment loading into and out of the system. An option flag in EUTRO-M also allows for the creation of a restart ”.RST" file. However, this option was not used in current analysis.
Since the input values used in the example problem are not representative of any actual system, the details of the output are not specifically discussed in this section. Instead, examples of the ’’.OUT”, ”.TRN" and ’’.MSB" files have been included in Appendix F. While the data specified in the example problem is purely hypothetical and should not be mistaken as representative of the Colorado River system, the output contained in these three files is helpful in explaining the type of output one should expect from EUTRO-M. The output would also be helpful in debugging if the model
is ever converted to a different computer system.
5.3.1 Effect of Finite-Difference Solution Technique on Results
As discussed in Section 4.3.4, the advection factor variable in EUTRO-M was set so that the model used a backward differencing scheme to solve the system derivatives. This was accomplished by setting the variable ”ADFAC” equal to zero. Although this produces the most stable solution, a nonzero advection factor could be used to reduce the numerical dispersion caused by a particular velocity, time step, and segment length combination. The consequences of this action can be seen in Figure 5.3.1.
The example input file was run three times adjusting only the advection factor. Concentrations of ammonia and orthophosphate were plotted against time and against segment number. As illustrated in Figure 5.3.1, altering the advection factor causes the program to calculate different steady state concentrations of both constituents in segment number one. However, by segment number three, the three simulations predict the nearly the same concentrations. Given this small amount of discrepancy, it was decided that model stability was a more important consideration and, therefore, the backward differencing technique was considered appropriate.
5.3.2 Effect of Time Step on Model Stability and Numerical Dispersion
The time step selected by the user may have a direct effect on the simulation results. In addition to the length of time required to perform a simulation, the chief area of concern is the trade-off between numerical dispersion and program stability. Therefore, the selection of a time step should not be an arbitrary decision. Instead, a trial and error approach should be used so that an appropriate time step may be specified.
In terms of both simulation time and numerical dispersion one would like to specify a large time step. Figure 5.3.2 illustrates how numerical dispersion decreases in a linear fashion as a function of the time step for each of the four surface water segments used in the example problem. Unfortunately, there is an upper limit to the size of the time step due to program instability. Figure 5.3.3 demonstrates what happens to ammonia and orthophosphate concentrations when the time step is increased gradually from 0.001 days to 0.025 days. For time steps of 0.001, 0.010 and 0.015 days, the simulations results are nearly identical and approach a steady state concentration in about a half a day. However, when the time step is increase to 0.020 days, the simulation results begin the oscillate. These oscillations propagate with time and eventually lead to numeric overflow errors.
This example illustrates the processes involved in determining the proper time step. As illustrated in Figure 5.3.2 there is a definite trade-off between the time step and model stability. In general, a modeler would like to select a time step as large as possible to reduce simulation time and reduce numerical dispersion. However, one needs to be cognizant of the fact that too large a time step will result in unstable solutions and ultimately overflow errors. This should also point out that the selection of a 0.001 day time step in the Colorado River analysis was not simply an arbitrary number pulled out of thin air. The use of larger time steps caused model instability and the use of smaller numbers increases simulation time.
Figure 5.1.1 Physical Layout of Example Problem
Figure 5.2.1 Schematic of Input Groups
Figure 5.3.1 Effect of Finite-Difference Solution Procedure
Figure 5.3.2 Effect of Time Step on Numerical Dispersion
Figure 5.3.3 Concentration Oscillations as a Result of Model Instability
Number Variable Definition 1 SEG DEPTH Depth in segment (m) 2 WIND Wind velocity over segment (m/s) 3 WATER VEL Water velocity in segment (m/s) 4 I TOT Incoming solar radiation (Ly/d) 5 SEG TEMP Temperature within segment (°C) 6 SEG TYPE Segment type (1, 2, 3 or 4) 7 PHYT Phytoplankton biomass carbon (mg/L) 8 RESP Phytoplankton respiration rate (1/d) 9 DEATH Phytoplankton death rate (1/d) 10 LIMIT Nutrient limitation indicator 11 TCHLAX Chlorophyll a concentration (mg/L) 12 XEMP1 Nitrogen limitation factor for phytoplankton 13 XEMP2 Phosphorus limitation factor for phytoplankton 14 GP1 Light/nutrient limited phytoplankton growth rate constant (1/day) 15 RLIGHT Light limitation factor for phytoplankton 16 RNUTR Nutrient limitation factor for phytoplankton 17 PNH3G1 Ammonia preference factor for phytoplankton Table 5.3.1 Output Display Variables for EUTRO-M (Continued) Number Variable Definition 18 NH3 Ammonia (mg/L) 19 NO3 Nitrate (mg/L) 20 ON Organic nitrogen (mg/L) 21 TIN Total inorganic nitrogen (mg/L) 22 TOT N Total nitrogen (mg/L) 23 TON Total organic nitrogen (mg/L) 24 CN Total inorganic nitrogen (mg/L) 25 OP Organic phosphorus (mg/L) 26 OPO4 Orthophosphate (mg/L) 27 TIP Total inorganic phosphorus (mg/L) 28 TOP Total organic phosphorus (mg/L) 29 RATIO Inorganic nitrogen to phosphorus ratio (mg/mg) 30 DO Dissolved oxygen (mg/L) 31 CBOD Carbonaceous biochemical oxygen demand (mg/L) 32 BOD5 5-day biochemical oxygen demand (mg/L) 33 UBOD Ultimate 30-day BOD (mg/L) 34 DOMIN Minimum diurnal dissolved oxygen (mg/L) 35 DOMAX Minimum diurnal dissolved oxygen Table 5.3.1 Output Display Variables for EUTRO-M Number (Continued) Variable Definition 36 CS Dissolved oxygen saturation (mg/L) 37 KDC Specific carbonaceous BOD deoxygenation rate (1/d) 38 DELO2 Diurnal dissolved oxygen variation (mg/L) 39 KA Reaeration rate (1/d) 40 SOD1D Sediment oxygen demand (gm/m 2 -d) 41 D3 Dummy print variable 42 D4 Dummy print variable 43 TSS Total suspended sediment (mg/L) 44 ALGAE Algae growth density (gm/m 2 ) 45 TMAC1 Macrophyte species 1 growth density (gm/m 2 ) 46 TMAC2 Macrophyte species 2 growth density (gm/m 2 ) 47 TMAC3 Macrophyte species 3 growth density (gm/m 2 )
Table 5.3.1 Output Display Variables for EUTRO-M