WATER POLLUTION CONTROL RESEARCH SERIES • 16070 DZV 02/71
ESTUARINE MODELING: AN ASSESSMENT
ENVIRONMENTAL PROTECTION AGENCY • WATER QUALITY OFFICE
-------
WATER POLLUTION CONTROL RESEARCH SERIES
The Water Pollution Control Research Series describes
the results and progress in the control and abatement
of pollution in our Nation's waters. They provide a
central source of information ,on the research , develop-
ment, and demonstration activities in the Water Quality
Office, Environmental Protection Agency, through inhouse
research and grants and contracts with Federal, State,
and local agencies, research institutions, and industrial
organizations.
Inquiries pertaining to Water Pollution Control Research
Reports should be directed to the Head, Project Reports
System, Office of Research and Development, Water Quality
Office, Environmental Protection Agency, Room 1108,
Washington, D. C. 20242.
For sate by the Superintendent of Documents, TJ.S. Government Printing Office
WMhlngton, B.C. 20402 - Price $4.60
Stock Number 5501-0129
-------
ERRATUM: 16070 DZV 02/71, S/N 5501-0129
THIS REPLACES THE TITLE PAGE (PG. i) IN THE ABOVE REPORT.
ESTUARINE MODELING: AN ASSESSMENT
Capabilities and Limitations for
Resource Management and Pollution Control
Edited by
George H. Ward, Jr. William H. Espey, Jr.
V
TRACOR, Inc.
6500 Tracor Lane
Austin, Texas 78721
for the
WATER QUALITY OFFICE
ENVIRONMENTAL PROTECTION AGENCY
Project 16070DZV
Contract 14-12-551
February, 1971
-------
ESTUARINE MODELING: AN ASSESSMENT
Capabilities and Limitations for
Resource Management and Pollution Control
by
TRAGOR, Inc.
6500 Tracer Lane
Austin, Texas 78721
for the
WATER QUALITY OFFICE
ENVIRONMENTAL PROTECTION AGENCY
Project 16070DZV
Contract 14-12-551
February, 1971
-------
EPA REVIEW NOTICE
This report has been reviewed by the Water
Quality Office, EPA, and approved for publi-
cation. Approval does not signify that the
contents necessarily reflect the views and
policies of the Environmental Protection
Agency, nor does mention of trade names or
commercial products constitute endorsement or
recommendation for use.
ii
-------
ABSTRACT
This report constitutes a technical review and critical appraisal of present techniques
of water quality modeling as applied to estuaries. Various aspects of estuarine modeling are
treated by a selection of scientists and engineers eminent in the field, and these essays are
supplemented by discussions from technical conferences held during the course of the report's
preparation. Topics discussed include one-, two-, and three-dimensional mathematical models for
estuarine hydrodynamics, water quality models of chemical and biological constitutents including
DO, BOD, nitrogen forms, phytoplankton and general coupled reactants, models of estuarine temp-
erature structure with special attention given the modeling of thermal discharges, and the
principles and applicability of physical models in estuarine analysis. Also included is a
review of solution techniques, viz. analog, digital and hybrid with a detailed discussion of
finite-difference methods, a brief survey of estuarine biota and present biological modeling
activities, and a collection of case studies reviewing several past estuarine modeling projects.
Conclusions about the existing state of the art of estuarine modeling and recommendations for
future research by EPA are advanced throughout the report and are summarized in the final
chapter.
This report was submitted in fulfillment of Project No. 16070DZV between the Environ-
mental Protection Agency and TRACOR, Inc., as part of the National Coastal Pollution Research
Program.
iii
-------
AUTHORS AND REVIEWERS
William E. Dobbins
Partner
Teetor-Dobblns, Consulting. Engineers
MacArthur Airport
Ronkonkoma, New York
John Eric Edinger
Associate Professor of Civil Engineering
Towne School of Civil and Mechanical
Engineering
University of Pennsylvania
Philadelphia, Pennsylvania
William H. Espey, Jr.
Director, Ocean Sciences and Water
Resources
TRACOR, Inc.
Austin, Texas
James A. Harder
Professor of Hydraulic Engineering
Department of Civil Engineering
University of California
Berkeley, California
Donald R. F. Harleman
Professor of Civil Engineering
Ralph M. Parsons Laboratory for Water
Resources and Hydrodynamics
Department of Civil Engineering
Massachusetts Institute of Technology
Cambridge, Massachusetts
Arthur T. Ippen
Institute Professor
Director, Ralph M. Parsons Laboratory for
Water Resources and Hydrodynamics
Department of Civil Engineering
Massachusetts Institute of Technology
Cambridge, Massachusetts
Jan J. Leendertse
The RAND Corporation
Santa Monica, California
Donald J. O'Connor
Professor of Civil Engineering
Department of Civil Engineering
Manhattan College
Bronx, New York
Gerald J. Paulik
Professor of Population Dynamics
Center for Quantitative Science in
Forestry, Fisheries and Wildlife
University of Washington
Seattle, Washington
Donald W. Pritchard
Professor of Oceanography
Director, Chesapeake Bay Institute
The Johns Hopkins University
Baltimore, Maryland
Maurice Rattray, Jr.
Professor of Oceanography
Chairman, Department of Oceanography
University of Washington
Seattle, Washington
Robert V. Thomann
Associate Professor of Civil Engineering
Department of Civil Engineering
Manhattan College
Bronx, New York
George H. Ward, Jr.
Ocean Sciences and Water Resources
TRACOR, Inc.
Austin, Texas
-------
WATER QUALITY OFFICE PARTICIPANTS
Donald J. Baumgartner
Chief, National Coastal Pollution
Research Program
Pacific Northwest Water Laboratory
Corvallis, Oregon
Richard J. Callaway
Chief, Physical Oceanography Branch
National Coastal Pollution Research
Program
Pacific Northwest Water Laboratory
Corvallis, Oregon
Kenneth D. Feigner
Systems Analysis and Economics Branch
Washington, D. C.
W. H. Fitch
Systems Analysis and Economics Branch
Washington, D. C.
Dan Fitzgerald
New England River Basins Office
Northeast Region
Needham Heights, Massachusetts
Thomas P. Gallagher
Southeast Water Laboratory
Athens, Georgia
Howard Harris
Southwest Region
Alameda, California
David R. Hopkins
South Central Region
Dallas, Texas
C. Robert Horn
Middle Atlantic Regional Office
Charlottesville, Virginia
Norbert A. Jaworski
Chesapeake Technical Support Laboratory
Middle Atlantic Region
Annapolis, Maryland
Malcolm F. Kallus
Federal Coordinator, Calveston Bay Project
South Central Region
Hous ton, Texas
Larry Olinger
Southeast Water Laboratory
Athens, Georgia
George D. Pence
Edison Water Quality Laboratory
Northeast Region
Edison, New Jersey
W. C. Shilling
Robert S. Kerr Water Research Center
South Central Region
Ada, Oklahoma
Roger D. Shull
Pollution Control Analysis Branch
Washington, D. C.
John Vlastelicia
Northwest Region
Portland, Oregon
T. A. Wastier
Estuarine and Oceanographic Programs Branch
Washington, D. C.
John Yearsley
Northwest Region
Por 11and, Oregon
VI
-------
ATTENDANCE AT CONFERENCES
24 June 1969
St. John's College
Annapolis, Maryland
D. J. Baumgartner
R. J. Callaway
W. E. Dobbins
W. H. Espey
K. D. Feigner
D. Fitzgerald
W. H. Fitch
T. P. Gallagher
J. A. Harder
D. R. Hopkins
C. R. Horn
R. J. Huston
N. A. Jaworski
M. F. Kallus
D. J. O'Connor
G. D. Pence
D. W. Pritchard
M. Rattray
W. C. Shilling
R. D. Shull
R. V. Thomann
J. Vlastelicia
G. H. Ward
T. A. Wastler
6 May 1970
Salishan Lodge
Gleneden Beach, Oregon
D. J. Baumgartner
R. J. Callaway
W. E. Dobbins
J. E. Edinger
W. H. Espey
K. D. Feigner
D. Fitzgerald
J. A. Harder
D. R. F. Harleman
H. Harris
D. R. Hopkins
C. R. Horn
N. A. Jaworski
M. F. Kallus
J. J. Leendertse
L. Olinger
G. J. Paulik
D. W. Pritchard
M. Rattray
W. C. Shilling
R. V. Thomann
J. Vlastelicia
G. H. Ward
T. A. Wastler
vii
-------
FOREWORD
The Water Quality Office of the Environmental Protection Agency has historically
supported research and development efforts for estuarine models, and is continuing to do so.
The purpose of this report is to provide a summary of presently available modeling methods used
in studying the causes of estuarine pollution problems and pollution abatement alternatives.
This current project has been undertaken by the National Coastal Pollution Research
Program as part of an overall mission to expand scientific understanding of the marine environ-
ment with respect to improvement and protection of marine water quality. This mission is
achieved by fundamental and applied studies on (a) physical transport and dispersion, (b bio-
logical, chemical and physical transformations and interactions, (c) water use damages,
(d) restoration of degraded marine waters, and (e) conservation of resources associated wir*
additions of heat and materials introduced from land, rivers, and the atmosphere.
The nature of estuarine water pollution problems demands consideration of solid wastes
which, through land use and ocean disposal, contribute materials to the estuarine system. Air-
borne materials can contribute to oceanic pollution—which can in turn influence or add to
existing estuarine concentrations. Estuaries act as a focal point for man's environmental
ills, screening, as they do, wastes discharged to inland watersheds and fed by small portions
of all the wastes accumulated in the oceans. With the establishment of the Environmental Pro-
tection Agency there now appears to be a framework for estuarine pollution control consistent
with the physical dimensions of the problem.
It was decided a report of this nature would be most beneficial if it were prepared
by a group of consultants and specialists with considerable experience with one or more modeling
topics. To make sure this assessment was directed to the problems pollution control agencies
were working on now and foresaw as impending problems, members of the Federal Water Quality
Administration regional office technical staffs were included.
This report presents a discussion of capabilities and limitations of both mathematical
and physical models, and of various solution techniques for the mathematical formulations.
Based on first-hand experience and a critical review of published reports by the authors and
participants, recommendations are made for research and development found necessary to over-
come deficiencies in present methods and to provide effective methods for handling more complex
problems in the future.
D. J. Baumgartner
R. J. Callaway
National Coastal Pollution Research Program
Pacific Northwest Water Laboratory
Water Quality Office
Environmental Protection Agency
Corvallis, Oregon
Vlll
-------
PREFACE
It Is appropriate to include here some prefatory remarks on what this report is and
what it is not, along with something about how it came to be. The purpose of this report is
to critically appraise the present state of the art of water quality modeling as applied to
estuaries. This purpose subsumes the question of how well do we understand estuarine processes,
and includes, further, the questions: how well can we predict these processes, how well need
we predict them for water quality regulation and management, and how do we best get from where
we are to where we need to be. The answers provided to these questions are not simple, as can
be judged from the bulk of this volume, nor in many cases are they complete or even clearly
defined. But this is a liability of attempting to obtain these answers, quickly enough to be
useful, from such a young, complex and vital field.
Although in a sense this report constitutes a cataloging of models, it is not intended
nor written to be a reference manual for the practicing engineer. Many of the arguments require
an examination of fundamentals, and return therefore to derivations and basic developments.
But assumed throughout is an at least superficial familiarity of the reader with estuarine
processes and their representation. Neither is this report a disjoint collection of technical
papers, but it is, rather, an organic sequence of essays on assigned topics written from a
consistent point of view.
The basic plan by which this report was prepared was to retain a group of consultants
with considerable experience and expertise in aspects of estuarine modeling, some of whom would
author the report and some of whom would provide a critical review. Preparation of the report
was coordinated insofar as the geographical distribution of the consultants permitted, and was
augmented by two general conferences attended by the consultants and those individuals in the
Water Quality Office of EPA with direct responsibility for the development or application of
models. These conferences were truly conferential, in that their purpose was to discuss in as
much technical depth as necessary the various aspects of the report and of water quality models
in general. In addition, throughout its preparation this report was open for discussions of
the contents and the inclusion of alternate viewpoints. The discussions, both formal and
informal, therefore constitute an integral and substantial part of the report.
The project was formally initiated in June 1969 and the first conference was held
within the month at St. John's College in Annapolis, Maryland. The purpose of the conference
was to agree on the content and scope of the report and to select the principal authors. The
Water Quality Office (then FWPCA) had submitted a suggested outline consisting of five chapters,
(i) mathematical models, (ii) analog and hybrid models, (ill) hydraulic models, (iv) use of
models in evaluating streamflow regulation, and (v) use of models in evaluating physiographic
modifications, which exemplifies the primary concerns of FWPCA at the time. At this conference
the federal programs in estuarine water quality modeling were reviewed, and the report outlined
in essentially its present form. (At the beginning of this work, the responsible agency was
the Federal Water Pollution Control Administration, which subsequently became the Federal
Water Quality Administration and is now the Water Quality Office of the Environmental Protection
Agency. The agency is therefore referred to variously as FWPCA, FWQA and EPA throughout the
report.)
ix
-------
After careful consideration of the report's intended content, scope and objectives,
the resources available for its preparation were apportioned in the following way: hydro-
dynamics 20%, water quality 207,,, temperature 67., physical models 10%, solution techniques 207,
case studies 10%, biological modeling 27,,, critical review 127,. While in fact the consultants
contributed more than their allotted time, this distribution does serve, however, to indicate
the relative weight given the various topics.
The authors of this report are men whose varied activities and scholarly pursuits
place a high demand on their time. With such contributors, printing deadlines are propositions
of naive optimism. (Out of similar undertakings has grown many an eloquent lament for the
editor who daily nurses an empty mailbox.) Although it was scheduled that the second conference
be convened in September and the report completed in November, or at least by the end of the
year, it was not until June of 1970 that the second conference met in Salishan Lodge at Gleneden
Beach, Oregon, nearly a year after the first conference. At this conference, the nearly com-
plete report was subjected to a thorough and wide-ranging discussion. Plans were made for
modifications and additions to the text and for its completion.
Editorial policy for the report was somewhat variable. The individual chapters were
only lightly edited. Although notations could be made uniform, duplication eliminated, and
the chapters modified to conform to a consistent style, we felt that in doing so, more would
be lost than gained, and that the overall value of the report would be greatest if the identity
and individuality of each author's work were preserved. The discussions, in contrast, required
heavy editing, as a precise rendering in print of spoken dialog is less a transcript than a
transmogrification. The contents of Chapter IX and the discussions appended to each chapter
have been freely rearranged to achieve a more useful organization. As a general rule, dis-
cussions which involve technical details or which are of specific pertinence to a particular
chapter are placed following that chapter, while discussions of a more general or recommenda-
tory nature are included in Chapter IX. In no instance, however, are statements separated
from the discussion context in which they occurred.
A number of acknowledgements are in order. Foremost are due D. J. Baumgartner and
R. J. Callaway of the National Coastal Pollution Research Program of EPA. Their interest in
the work and their resoluteness in achieving the project objective are in large measure respon-
sible for the quality of this report. We are grateful for their assistance and support
throughout the project.
The assistance and cooperation of our firm TRACOR, Inc., during the project is
gratefully acknowledged. Several of our associates here, especially R. J. Huston and
J. P. Buckner (now with the Coastal Bend Regional Planning Commission, Corpus Christ!), con-
tributed their assistance throughout the study. We wish to thank T. A. Wastler (Estuarine
and Oceanographic Programs Branch, EPA) for his interest and comments, and EPA scientists and
engineers, Norbert Jaworski, Howard Harris, George Pence, Tom Gallagher, Ken Feigner and
John Vlastelecia, for their individual assistance, particularly in the compilation of the case
studies (Chapter VII). E. G. Fruh (university of Texas) and B. J. Copeland (North Carolina
State University) read portions of Chapter VII and their reviews lead to many substantive improve-
ments. We are grateful to Colonel Frank P. Bender, Director of the Galveston Bay Project, for
permitting the use in Chapter VII of results from the Galveston Bay Project. Dr. Garbis H.
Keulegan of the Waterways Experiments Station has followed the preparation of the report since
its inception, although owing to other commitments he was unable to formally review the docu-
ment. We are grateful to him for his comments and for his great interest in the project.
-------
We appreciate receiving permission from the following publishers and authors to
reproduce copyrighted material: the American Society of Civil Engineers for several figures in
Chapters IV, V and VII which appeared in the Proceedings. McGraw-Hill Book Company for Figure
2.10 from Estuary and Coastline Hydrodynamics, the Water Pollution Control Federation for
Figures 3.11 - 3.17 from the Journal Water Pollution Control Federation, the American Fisheries
Society for the quotation of Prof. Hedgpeth in Chapter VIII from A Symposium on Estuarine
Fisheries. Donald Watts and 0. L. Loucks of the University of Wisconsin for Figure 8.1 from
Models for Describing Exchanges Within Ecosystems, the Controller of Her Britannic Majesty's
Stationery Office for the figures and tables in Chapter VII from Effects of Polluting Discharges
on the Thames Estuary, and the Controller of H. M. Stationery Office and the authors A. L. H.
Gameson and I. C. Hart for Figure 7.29 and Table 7.4 from "A Study of Pollution in the Thames
Estuary" in Chemistry and Industry. The data in Figure 7.29 and Table 7.4 are those of the
Greater London Council.
Finally, we wish to express our gratitude to the consultants, Drs. Dobbins, Edinger,
Harder, Harleman, Ippen, Leendertse, O'Connor, Paulik, Pritchard, Rattray and Thomann, who
retained their patience and generosity through the eighteen long months this report was in
preparation. This is their report. Its value and its merits are due to their integrity, their
industry, and their command of their fields.
G.H.W.
W.H.E.
Austin
November 1970
XI
-------
CONTENTS
Abstract 111
Authors and Reviewers v
Water Quality Office Participants vl
Attendance at Conferences Vll
FOREWORD D. J. Baumgartner and K. J. Callaway Vlll
PREFACE 1X
I. INTRODUCTION S. S. Hard and V. S. Bepey
1. Estuaries 1
2. Estuarine Models 2
3. Report Content and Organization 4
References 4
II. HYDRODYNAMIC MODELS
1. Three-Dimensional Models D. H. Pritahard 5
1.1 Introduction 5
1.1.1 The Coordinate System 5
1.1.2 The Concept of the Ensemble Average: S
Deterministic vs. Stochastic'Modes of Motion
1.1.3 The Notation 7
1.2 The Basic Equations 8
1.2.1 Conservation of Mass 8
1.2.2 Conservation of Momentum 10
1.2.3 Conservation of Dissolved Constituents 16
1.2.4 Some Possible Further Simplification of the 17
Three-Dimensional Conservation Equations for
Mass, Momentum and Salt
2. Two-Dinensional Models D. V. Pritahard 22
2.1 Introduction 22
2.1.1 The Coordinate System and the Notations to be Used 22
2.2 The Basic Equations for a Vertically Averaged 22
Two-Dimensional Model
2.2.1 Conservation of Mass 22
2.2.2 Conservation of Momentum 24
2.2.3 The Vertically Averaged Two-Dimensional 27
Salt Balance Equation
2.2.4 Summary of the Set of Vertically Averaged Dynamic and 28
Kinematic Equations Suitable to Serve as the Basis of
a Two-Dimensional Numerical Model of an Estuary
xii
-------
2.2.5 An Alternate Development of a Set of Vertically Averaged 30
Dynamic and Kinematic Equations in an Estuary
References for Sections 1 and 2 33
3. One-Dimensional Models D. P. F. Harleman 34
3.1 Introduction 34
3.2 Mathematical Models in Real Time 36
3.2.1 Mass Transfer Equations for a Non-Conservative Substance 37
3.2.2 Source and Sink Terms 39
3.2.3 Tidal Velocity 40
3.2.3.1 Continuity Equation for a Variable Area Estuary 42
3.2.3.2 Momentum Equation for a Variable Area Estuary 43
3.2.3.3 Solution Techniques for Tidal Hydraulic Problems 45
3.2.4 Longitudinal Dispersion in the Real Time Mass 53
Transfer Equation
3.2.4.1 Longitudinal Dispersion Coefficient in Regions 55
of Uniform Density
3.2.4.2 Longitudinal Dispersion in Salinity 57
Intrusion Regions
3.2.5 Discussion of the Real Time Mass Transfer Equation 61
3.2.6 Solution of the Real Time Mass Transfer Equation 63
in a Uniform Estuary
3.3 Mass Transfer Equations Using Non-Tidal Advective Velocities 67
3.3.1 Mass Transfer Equation: Time Average over a Tidal Period 67
3.3.2 Mass Transfer Equation: Slack Tide Approximation 67
3.3.3 Analytical Solutions of the Non-Tidal Advective 68
Mass Transfer Equations in a Uniform Estuary
3.3.4 Mathematical Models of Salinity Intrusion 70
3.3.4.1 Descriptive Studies 71
3.3.4.2 Predictive Studies 71
3.4 Comparison of Real Time and Non-Tidal Advective 76
Mathematical Models
3.4.1 Continuous Inflow of a Non-Conservative Pollutant 76
3.4.2 Salinity Intrusion 78
3.5 Summary and Conclusions 79
3.5.1 Mathematical Models in Real Time 79
3.5.2 Mass Transfer Equations Using Non-Tidal Advective 81
Velocities
3.5.3 Conclusions 81
References for Section 3 85
Discussion 90
WATER QUALITY MODELS: CHEMICAL, PHYSICAL AND
BIOLOGICAL CONSTITUENTS D. J. O'Connor and R. V. Thomann
1. Introduction 102
2. One-Dimensional Analysis 106
2.1 Solution Due to Instantaneous Release 107
2.2 Conservative Substances 109
2.3 Non-Conservative Substances 110
2.4 Consecutive Reactions (BOD-DO) 112
Kill
-------
2.5 Applications to Dissolved Oxygen Analysis 113
2.6 Multi-Stage Consecutive Reactions (Nitrification) 117
2.7 Data Requirements 119
3. Solution Techniques for One-Dimensional Steady State 123
3.1 Available Approaches 123
3.2 Continuous Solution Approach 125
3.3 Finite Section Approach 127
4. Two-Dimensional Steady State 131
4.1 Solution Approach: Finite Difference 131
5. One-Dimensional Time-Variable Models 135
5.1 Solution Approach 136
6. Verification 140
6.1 An Example: The East River (O'Connor 1966) 141
6.1.1 Description of the System 142
6.1.2 Field Data 144
6.1.3 Analysis of Data 145
6.1.4 Waste Water Discharges 147
6.1.5 BOD and DO Deficit Profiles 148
6.1.6 Discussion 151
7. Future Directions 152
7.1 Steady-State Feedback Modeling 152
7.2 Dynamic Phytoplankton Model 156
8. Discussion and Recommendations 163
References 167
Discussion I69
IV. ESTUARINE TEMPERATURE DISTRIBUTIONS J. E. Edinger
1. Natural Temperature Distributions 174
1.1 Similarity Between Salinity and Temperature Distributions 174
1.2 Influence of Surface Heat Exchange on Temperature Distributions 175
1.3 Additional Influences on Temperature Structure 180
2. Analytic Description of Estuarine Temperature Structure 182
2.1 The Vertically Mixed Case 182
2.2 Two-Layered Segmented Models 183
2.3 Continuous Vertical Temperature Structure 185
3. Modeling Temperature Distributions due to Large Heat Sources 187
3.1 The Temperature Excess 187
3.2 Initial Temperature Distributions 188
3.3 Intermediate Temperature Distributions 191
3.4 Large-Scale Temperature Distributions 194
3.5 Physical Hydraulic Modeling of Thermal Discharges 200
References 202
Discussion 207
XIV
-------
V. PHYSICAL HYDRAULIC MODELS D. K. F. Harleman
1. Introduction 215
2. Similitude of Momentum Transfer Processes in Tidal Motion 218
2.1 Physical Boundary Conditions for Tidal Models 218
2.2 Froude Scale Ratios for Tidal Motion 221
2.3 Reynolds Number Effects 225
3. Similitude of Mass Transfer Process 226
3.1 Similitude of Salinity Intrusion 226
3.1.1 Saline-Wedge Estuary 228
3.1.2 Mixed Estuary 230
3.2 Mass Transfer Similitude in Regions of Uniform Density 233
4. Model Verification 238
4.1 Tidal Verification 238
4.2 Salinity Verification 241
4.3 Sediment Transport Verification 241
5. Discussion of the Use of Physical Hydraulic Models in Estuarine 247
Water Quality Studies
6. Summary and Conclusions 251
References 252
Discussion 255
VI. SOLUTION TECHNIQUES
1. Analog and Hybrid Techniques j. A. Harder 264
1.1 Introduction 264
1.2 Analog Modeling 265
1.2.1 Analog Simulators and Analog Computers 265
1.2.2 Analog Simulators for Tidal Flows 265
1.2.2.1 Simulating Hydraulic Friction 269
1.2.3 An Analog Computer Element for Pollution Dispersion 270
1.3 Hybrid Computation 271
1.4 Economic Considerations in Analog Computing 272
1.5 System-Response Type Models 273
References for Section 1 276
2. Digital Techniques: Finite Differences J. J. Leendertee 277
2.1 Introduction 277
2.2 Considerations in Design of Computation Schemes 278
2.2.1 Conservation of Mass 278
2.2.2 Stability 279
2.2.3 Dissipative and Dispersive Effects 281
2.2.4 Effects of Noncentered Operations 284
2.2.5 Discontinuities 285
2.2.6 Comments on Higher Order Schemes 286
2.2.7 Use of Conservation Laws to Obtain Nonlinear Stability 287
2.3 Tidal Computations 290
2.3.1 One-Dimensional Computations 290
2.3.2 Two-Dimensional Tidal Computations 291
2.3.3 Numerical Simulation of Tidal Flow in Two Dimensions 295
XV
-------
2.4 Mass Transport Computations 2^6
2.4.1 One-Dimensional Water Quality Computations 296
2.4.2 Two-Dimensional Water Quality Computations 298
References for Section 2 301
Discussion
VII. CASE STUDIES 0. H. Hard and U. H. Espey
1. Introduction 31°
2. The Thames 311
2.1 Physiography and Hydrography 311
2.2 Polluting Loads 311
2.3 Modeling Techniques 317
2.3.1 Modeling of Temperature Distribution 323
2.3.2 Modeling of Dissolved Oxygen Distribution 325
2.4 Model Results 327
2.4.1 Calculated Profiles and Their Accuracy 327
2.4.2 Model Application 332
3. The DECS Model 341
3.1 The Delaware Estuary 341
3.1.1 Characteristics of the Estuary 342
3.1.2 Model Application and Results 542
3.2 The Potomac Estuary 348
3.3 Hillsborough Bay 35A
3.4 Summary 356
4. San Francisco Bay 361
4.1 Characteristics of the Bay-Delta System 361
4.1.1 Physiographic and Hydrographic Features 361
4.1.2 Water Quality 363
4.2 The Physical Model 366
4.2.1 Verification 368
4.2.2 Evaluation of the Effects of Solid Barriers on 371
Hydraulics and Salinities
4.2.3 Application of Dye Releases 377
4.3 Digital Computer Model 383
4.3.1 Theory and Computational Framework 383
4.3.2 Model Application and Verification 390
4.3.3 Comparison of Physical Model and Digital Computer Model 396
5. Galveston Bay 3"
5.1 Galveston Bay Project Mathematical Models 399
5.1.1 Hydrodynamic Model 401
5.1.1.1 Model Formulation and Application 401
5.1.1.2 Comparison of Mathematical Hydrodynamic 404
Model with Physical Model
5.1.2 Salinity Model 408
5.1.3 Thermal Discharge Model 411
XVI
-------
5.1.4 Application of the GBP Models to Ecological Parameters 420
5.2 Physical Models 424
References 434
VIII. BIOLOGICAL MODELING IN ESTUARIES: A NOTE C. J. Paulik
1. Characteristics of Animal Communities in Estuaries 438
2. Mathematical Models of Biological Communities 442
3. Use of Numerical Models to Conserve, Exploit and Exterminate 446
Animal Populations
4. Recommendations for Biological Modeling in Estuaries 447
4.1 Key Species Models 447
449
References
IX. DISCUSSIONS
1. Estuarine Hydrodynamic and Water Quality Models 451
2. Computational Aspects 4*"5
3. Data Collection and Verification 47°
4. Prepared Remarks 48°
4.1 Discussion G. a. Ward 48°
4.2 Discussion 7?. J. Callaway 482
4.3 Some Projections for the Immediate Future J. J. Leendertse 484
4.4 Discussion W. B. Espey 487
5. Review and Comments A. T. Ippen 489
X. CONCLUSIONS AND RECOMMENDATIONS
AQO
1. Conclusions *
1.1 Estuarine Hydrodynamic Models 492
1.2 Estuarine Water Quality Models 493
1.3 Utility of Physical Models 494
1.4 Solution Techniques 494
1.5 Applications 495
2. Recommendations for Future Research 496
2.1 Hydrodynamic Models 496
2.2 Water Quality Models 496
xvii
-------
CHAPTER I
INTRODUCTION
George H. Ward, Jr.
and
William H. Espey, Jr.
Growing concern with the preservation of our water resources, particularly within the
last decade, has precipitated the need for a balanced and astute management of these resources.
Even disregarding the political intricacies, water quality management is an exceedingly
difficult task. The difficulties arise from two essential facts. First, utilization of the
resources of a natural waterbody is subject, almost universally, to conflicts among potential
users, in that one use inherently delimits or excludes another. The problems are compounded
since, typically, each potential use is of some socio-economic importance. Much, if not
most, of the management task ultimately devolves to the resolution of these conflicts in a
satisfactory manner. Secondly, the processes--hydrographic, chemical, biological—which operate
in a waterbody are extremely complex and often poorly understood. Yet it is these processes
upon which management decisions must be based.
The application of water quality models has proved to be a powerful technique in
water resource management, as models can incorporate the complexity of the relevant processes
in the waterbody into a utilitarian form for management consideration. Through the use of
models both a diagnostic and predictive capability is provided, diagnostic in the sense of
permitting the identification and isolation of specific factors affecting the water quality,
and predictive, in permitting evaluation of future effects of proposed changes in the water-
body, its inputs, or its uses.
1. ESTUARIES
This report is concerned with water quality modeling as applied to a specific type of
waterbody, the estuary. (The general features of estuarine hydrography, viz., the defining
characteristics of estuaries, associated circulation patterns, physiography, etc., are well
known and may be assumed ab initio. See, e.g., Pritchard 1952, Cameron and Pritchard 1963,
Lauff 1967.) This apparently superficial circumscription perhaps deserves discussion.
Estuaries are, collectively, of singular importance. They are important biologically;
in contradistinction to the nutrient-rich water-poor land and the water-rich nutrient-poor
ocean, estuaries are typically both a nutrient-rich and water-rich environment, highly produc-
tive, constituting the prime habitat for a myriad of species and the nursery and spawning area
for many more. Besides their biological significance, however, estuaries are of specific
-------
importance to man, as is evidenced by their propensity for stimulating urban development.
Estuaries offer an accessibility to the sea, often with excellent shipping and harbor facilities,
while the feeding rivers provide a source of freshwater for domestic and industrial use, as well
as affording receiving water for municipal and industrial wastes.
Estuaries are similarly quite susceptible to pollution and degradative modifications
in quality. They often receive excessive volumes of wastes both from peripheral discharges and
from upstream discharges carried into the estuary by the inflows, and are often adversely
influenced by upstream freshwater consumption. In addition to affecting the estuary's own
delicate ecology, the effects of these practices may extend beyond the physiographic bounds of
the estuary itself.
Besides these very pragmatic reasons for a concern with estuaries, the estuary per
se poses modeling problems of extreme scope and complexity. Its hydrodynamics, for instance,
are determined by the interaction of the tides from the ocean and the influx of fresh water
from the rivers, modified by complicated boundary stresses due to the semi-enclosed physiography,
and further altered--or even dominated, as in fjords—by gravitational circulations arising from
density gradients. These factors, however, are sufficient to characterize only the transport
of a specific constituent in the estuary: the various processes and reactions to which that
constituent is subjected must be determined as well. These range from simple first-order
kinetics to reactions which depend nonlinearly upon additional variables and are possibly linked
with the biota.
It is evident therefore that estuaries, besides comprising a water resource of
critical importance, present modeling problems of intrinsic difficulty. This report attempts
to summarize these problems and assess present methods for their treatment. Although many of
the topics discussed can find application as well in coastal zones, In streams and rivers, or
in lakes, the discussions are inclined to those aspects fundamental to, or peculiar to,
estuarine processes.
2. ESTUARINE MODELS
The meaning of the term "model" as used in this report Is not precise, and, indeed,
four distinct meanings of the word as employed with regard to estuarine water quality may be
identified. It is appropriate here to indicate these usages in a somewhat pedestrian manner,
avoiding any judgments on their adequacy or correctness as such judgments can entail a lengthy
philosophical digression. (Recent considerations of the meaning and role of models In science
may be found In, e.g., Nagel, Suppes and Tarski 1962.)
In the strictest sense, a model is a formalism: the mathematical expression of the
relevant physical processes, e.g., the application of the principles of Newtonian fluid mechanics
with approximations appropriate to an estuary. A different but related usage is that a model
is a conceptual idealization or simplified representation of a physical process. The employment
of such a model is usually motivated by tractabillty and is Justified by the ability of the
model to reproduce observed results, even though It is ultimately unfaithful to the basic physics.
One example is the representation of the Reynolds flux of salinity
-------
the surface tension—even though the surface is known to be a transition region of high density
gradient.)
"Model" is also used to mean a simple (in some sense) analog of the real system. In
estuarine studies, the most familiar example is, of course, the physical or hydraulic model, a
highly distorted, small-scale reproduction of the estuary complete with freshwater inflows, a
sump of ocean salinity, and mechanical tide generators. This use of the term could as well
apply to a simulator consisting of electrical components, as described in Chapter VI.
Finally, a "model" may refer to a mathematical formulation (a "model" in either the
first or second sense above) together with a solution technique for the variables of concern.
This usage is implicit in the expressions "analog model" and "digital model," and may in some
cases even refer to a specific computer program, such as DECS-HI or the WRE Model (the "Orlob
Model"). This usage, which has some obvious faults, does not enjoy general approval by
scientists and engineers involved in estuarine water quality, but nevertheless has become more
frequent in recent years, particularly on the administrative or program level.
In the following pages all four uses of the term can be noted, sometimes in almost
confusing contiguity. Probably the only definition of "model" sufficiently general to encompass
these variations in meaning is the (superficial) statement that an estuarine model is a technique
or device for the representation of some property in an estuary.
The preponderance of present estuarine water quality models is concerned with the
concentration or distribution of a constituent, property or parameter in the estuary.
Accordingly, the processes represented in the model are generally classified as either transport
processes or reaction processes. Transport processes are basically hydrodynamic and include
advection, turbulent diffusion, and, when spatial averaging is involved, dispersion. Reaction
processes encompass the sources and sinks to which the parameter is subjected and may be
physical, chemical or biological, e.g., sedimentation and flocculation of organics, uptake of
oxygen in biochemical degradation of wastes, coliform mortality, algae loss by zooplankton
grazing, etc. The relative importance of transport processes and reaction processes obviously
varies from estuary to estuary, and depends upon the parameter modeled as well as the required
spatial and temporal refinement.
The question eventually arises in the development of an estuarine model as to
whether all of the relevant processes in the estuary have been identified or correctly re-
presented. Classically this is determined by model verification, the comparison of model
predictions with empirical measurements. As might be expected in a state-of-the-art assessment,
the topic of verification pervades this report. But a word of qualification (and, perhaps,
warning) is due, in that "verification" is often employed in a wider sense, and may mean the
use of empirical data for anything from calibration of the model to its actual validation. The
context must be considered, for a 'Verified" model may not be all it appears.
3. REPORT CONTENT AND ORGANIZATION
The basic arrangement of this report is to progress from model formulation to
solution to specific application. Such a segregation of topics, of course, cannot be strictly
observed and there is therefore considerable immixture. Estuarine hydrodynamics, which determine
the transport of substances in the estuary, are considered first and organized from the stand-
point of spatial dimensionality. Both advective and diffusive processes are discussed. Here,
also, the discussion of salinity distribution is undertaken, not only because of the value of
-------
salinity as a natural tracer, but because of its considerable influence on estuarine circulation.
Then the modeling of water quality parameters per se is treated, which involves the specifica-
tion of source and sink terms and their incorporation with the transport mechanisms into various
conservation relations. The discussion is divided into a treatment of general chemical and
biological parameters (Chapter III) and a treatment of temperature (Chapter IV). Temperature
is discussed separately for several reasons. It is a parameter which is currently receiving
particular attention because of its profound ecological importance. It is a parameter which,
like salinity, exhibits significant gradients in estuaries due to the mixing of river and
ocean water. Further, because of the dramatic dynamics of thermal discharges, e.g. the large
volumes of flow involved and the effects of buoyancy, their modeling requires special
consideration.
An appraisal of the use of physical hydraulic models in estuarine quality is
presented next (Chapter V) taking advantage thereby of the formulations of mass and momentum
transfer developed in the preceding chapters. Solution techniques are discussed in Chapter VI,
separated into analog and analog-hybrid techniques, which accomplish the integration of the
basic equations continuously (in one independent variable) using electronic components, and
numerical methods, in particular finite-difference techniques which involve integration of the
equations using finite increments in the independent variables and principally find application
with high-speed digital computers. In the latter topic, consideration is given to accuracy and
stability problems as well as numerical phenomena such as false diffusion and non-conservation
of mass or energy. Several case studies are presented in Chapter VII to provide exonples of
recent modeling investigations of estuarine water quality. The selection is not intended to be
exhaustive but rather to exemplify the present state of model application to estuarine problems.
Finally, a brief overview of the estuarine animal populace and present trends in biological
modeling in estuaries is given in Chapter VIII. At present, there is very little effort
toward incorporating or linking the higher life-forms in estuaries to water quality models,
but this is indubitably one of the ultimate directions of estuarine modeling.
In a strict sense, a complete state-of-the-art report on estuarine water quality
modeling would encompass a large part of existing knowledge in chemical and biochemical
kinetics, the dynamics of nonhomogeneous stratified fluids, boundary and transfer processes
at the air-sea interface, modeling of biological and ecological systems, turbulent transport,
and so on. For evident reasons, such an ambitious approach was not undertaken here. The
emphasis throughout this report is, rather, on the interaction of these phenomena and its
formulation. Indeed, the uniquely complex and difficult estuarine system can probably be best
characterized by that feature, interaction.
REFERENCES
Cameron, W. M. and D. W. Pritchard, 1963: Estuaries. Chapter 15, The Sea. _II, M. N. Hill (Ed.).
New York, Interscience Publishers.
Lauff, George H. (Ed.), 1967: Estuaries. American Association for the Advancement of Science,
Washington, D. C.
Nagel, Ernest, Patrick Suppes and Alfred Tarski (Ed.), 1962: Logic. Methodology and Philosophy
of Science. Stanford University Press.
Pritchard, D. W., 1952: Estuarine hydrography. Advances in Geophysics. I. New York, Academic
Press.
-------
CHAPTER II
HYDRODYNAMIC MODELS
1. THREE-DIMENSIONAL MODELS
D. W. Pritchard
1.1 INTRODUCTION
In this section the basic time-dependent equations in three spatial dimensions,
expressing the conservation of mass, of momentum, and of a dissolved or finely divided suspended
constituent of the waters in an estuary, are developed. In order to keep this section from
becoming overly long, the development will of necessity be abbreviated, with frequent use of
the phrase "it can be shown." However, some outline of the development of the basic equations
from first principles is desirable, in order to clearly indicate what terms and processes are
neglected in arriving at simplified forms of the equations which offer some possibility of
solution, either in terms of our present knowledge or of knowledge which might be gained in the
foreseeable future.
This development will include arguments for the neglect of certain terms, particularly
certain terms that arise in the averaging process. Again, these arguments must be brief, with-
out extensive observational proof of the statements of relative size of the pertinent terms
being considered.
1.1.1 The Coordinate System
Tensor notation will be employed in this section because of the economy it affords
in respect to the number of terms which must be written down in expressing the complex differ-
ential equations which develop from the first principles to be considered here. Thus, a
right-handed rectilinear coordinate system is envisioned with the three axes xt (i = 1,2,3)
oriented such that xl and x2 are in a horizontal plane, and x3 is the vertical axis
directed upwards (i.e., parallel to the force of gravity but positive in the direction of
-g ). The plane of x, = 0 is considered to coincide with mean sea level.
H.2 The Concept of the Ensemble Average: Deterministic vs. Stochastic Modes of Motion
The basic conservation principles will first be expressed in terms of the instan-
taneous field of motion and of the distribution of the pertinent properties of the water in
the estuary such as density, salinity and pressure. It has long been recognized that for the
real, turbulent conditions of a natural waterway these equations are intractable. Some type of
-------
averaging is required in order to separate what might be called the deterministic modes of
motion from that part of the motion which is by nature stochastic in character. The method of
averaging has considerable influence on the meaning of the individual terms of the resulting
equations.
The type of averaging to be Initially employed here is based on the following con-
ceptual argument. Assume that the instantaneous fields of motion and of the pertiment
properties could be observed, In detail, within an estuary over some time period. During this
time period the external processes, i.e., tidal variation at the mouth of the estuary, fresh-
water Inflow, radiation balance and meteorological conditions, are defined—that Is, are known
functions of time. Now assume that the instantaneous fields of motion and of the pertinent
properties in this same estuary were observed over another time period of the same length as
the first during which the time variations of the external processes were the same as in the
first period, and that the time histories of all the external processes prior to this second
period were the same as those prior to the first period, back as far in time as is required by
the Inherent dynamic "memory" of the system.
If such measurements were possible, it would be found that the fields of motion and
of the distribution of pertinent properties would differ in detail from one period to the
other, despite the equality of the time variations of the external processes. If k such
"identical" periods occurred, the Instantaneous field of motion and of the distribution of
pertinent properties for any one period would differ, in detail, from all the others.
A part of the instantaneous field of motion would be common to all such k "Ident-
ical" periods, and this part is called the deterministic mode of motion. A part of the
instantaneous field of motion, representing the difference between the instantaneous field of
motion for any one of the k "identical" periods and the deterministic mode, constitutes the
stochastic mode of motion. These "turbulent" departures in the instantaneous velocity field,
at any point in space and at any time relative to the start of each of the k "identical"
periods, are distributed among the k periods In accordance with some probability distribu-
tion- -that Is, they are "random" in character, though not necessarily "gaussian". A similar
deterministic mode and stochastic mode would exist for each of the Instantaneous distributions
of the pertinent properties such as density, salinity and pressure.
Note that each of these conceptual k "identical" periods starts at the same point
in the identical time histories of the external processes; for example, at the same tidal
epoch. If, at a given point in space, the time records of the velocity for each of the k
periods were lined up relative to this starting time, and the average at a given point in time
were taken across all k records, then the resulting average is called the ensemble average.
Such an ensemble average is time dependent, since it is taken at every point in time over all
the k periods. A similar operation would provide ensemble averages for the pertinent para-
meters such as density, salinity and pressure. If k is sufficiently large, the ensemble
averages may be taken as suitable approximations to the deterministic modes of motion and of
the distribution of the pertinent properties.
It has been more common practice to obtain averages of the instantaneous equations
by taking such averages over a time period AT small compared to the time scale at which one
wishes to resolve the field of motion or the distribution of properties. The problem with
this approach is that the stochastic or turbulent mode of motion encompasses a wide range of
time and spatial scales.
-------
The ensemble average has the advantage of smoothing out the stochastic mode of
motion at all time scales, while also retaining the variations in the deterministic mode at
all time scales. Furthermore, this is the average used in most theories of turbulence. The
fact that it is only conceptually attainable is not a disadvantage. The predicting equations
which can be developed will in any case provide an estimate of the most probable field of
motion and/or of properties in the estuary, and not the actual instantaneous velocities and
instantaneous values of the pertinent properties which would occur at a real instant of time.
1.1.3 The Notation
The instantaneous velocity at a point x^ and a time t during one of the k
"identical" periods is designated by u. . We then define
"i " ui + ui' (2a)
and
where Uj is the ensemble average of u^ over all k "identical" periods, and u^ is the
departure of the instantaneous velocity for any one of the k periods from the ensemble
average. Note that ( >k represents the operation of obtaining the ensemble average. Further
note that
-------
u = che ensemble average of Uj over all k.
s - Che ensemble average of 1 over all k .
p - Che ensemble average of "p over all k .
p - che ensemble average of p over all k .
£ - the ensemble average of ? over all k .
< ). » symbol for the ensemble averaging operation.
u,' » the "curbulent" deparcure of the instantaneous velocity from its ensemble
average.
s1 - the "turbulent" departure of the instantaneous salt concentration from its
ensemble average.
p' - the "curbulent" deparcure of the instantaneous density from its ensemble
average.
p' • the "turbulent" deparcure of the instantaneous pressure from its ensemble
average.
£' « the "turbulent" departure of the instantaneous concentration of a constituent
from its ensemble average.
1.2 THE BASIC EQUATIONS
1.2.1 Conservation of Mass
The basic statement of the principle of conservation of mass is simply that the mass
of a moving fluid particle, that is, a particle of fluid whose boundaries move with the
instantaneous velocity, remains constant. Thus
where M is the mass of the moving fluid particle. This statement leads to the differential
equation
As discussed in Sections 1.1.2 and 1.1.3 above, the ensemble average of this
equation is taken, after substituting
' ui
p " P + P
-------
resulting in
pl>-° (2-5)
Note that while the ensemble means of u.^' and p1 are zero, the ensemble mean of
the cross-product u.' p' is not necessarily zero. However, note that in this case
'>k = V,P'(k <"
where v i i designates the correlation coefficient between u^ and o' . The maximum value
the correlation coefficient can have is 1.0. The root-mean-square velocity deviation is on
the order of 10" l of the ensemble average, and the root-mean-square density deviation is at
most on the order of 10" 3 of the ensemble average density. Consequently
l' P'>k
< 10-
*
and the last term of Equation (2.5) can therefore be neglected, giving finally
(2.6)
(Note that this argument is essentially the same as the Boussinesq approximation.)
Making use of the equality
.
at
Equation (2.6) can also be written as
0 (2.7)
It is usually stated that, for an incompressible fluid, g£ = 0 and hence, since the
compressibility of water is very small, the equation of continuity applicable to an estuary
can be closely approximated by
-------
that is, the three-dimensional flow field must be non-divergent. This equation might be more
appropriately called the equation of volume continuity. However, more than simply incompressi-
bility must be invoked to obtain (2.8) from (2.7). The density of estuarine water is a function
of temperature, salinity and pressure. Thus
dp fdp\ dp . /Sp^ ds , /Sp^ dB
3t * Vdp/ »J dT vSe/ a~t
The statement that a fluid is incompressible means that (||) - 0. However, this
statement does not mean that (||) and (||) are zero. To obtain Equation (2.8) from
Equation (2.7) it is necessary to invoke a Boussinesq-type approximation; that is, that time
and spatial variations in density may be neglected except in terms in which they are multiplied
by gravity. In any case, Equation (2.8) is in fact a good approximation to the equation of
continuity for use in estuaries.
1.2.2 Conservation of Momentum
Newton's Second Law is the basic statement of the principle of conservation of
momentum; that is, the time rate of change of momentum of a moving elemental fluid particle is
equal to the sum of the forces acting on the particle. The differential equation expressing
this principle is called the equation of motion. For estuaries the instantaneous equation of
motion is given by
Ft <* V + $ \ V * - -2e ° P " + > « + * <2'9>
where e.. is the cyclic tensor, gt the acceleration of gravity, u the kinematic molecular
viscosity, and 0, the component of the angular velocity of the earth in the x. direction.
Now, substituting for the instantaneous value of each of the variables the sum of
the ensemble average value and a deviation term, that is
"i * ui + ui
P - P + P1
P • P + P1
then the x,, x~ and x., components of the ensemble-averaged equation of motion are
ut' u,')k (2.10)
u2' Uj'>k (2.11)
dtid
10
-------
where f - 2O, - 2n sin ? is the coriolls parameter. Here fi is the angular velocity of the
earth and k = p ) and from k = 0 .
Next, note that p1 does not appear in the above equations because the Boussinesq
approximation has been applied. The coriolis terms involving the ensemble average vertical
velocity u, have been neglected in the two horizontal component equations. Finally, all
acceleration terms, both local and inertial, are neglected compared to gravity in the vertical
component. Thus the vertical component equation (2.12) is reduced to the hydrostatic equation.
Note that by virtue of the equation of continuity, Equation (2.8), the field acceleration term
(the second term) in each of the horizontal components of the equation of motion can be written
in Equation (2.10), and
in Equation (2.11).
Consider now the pressure force term in Equation (2.10). We proceed by integrating
the hydrostatic equation, Equation (2.12), from the depth *3 , where the pressure equals p ,
to the water surface, x3 - TI , where the pressure equals pa , the atmospheric pressure.
That is,
f £r dx3' - -s f
J 9Xs J
Pa - P - -8 I
Hence
a + g a f pdx,' (2-13)
8 " *"
11
-------
Applying Leibnitz' Rule to the last term of (2.13) gives
*1 (2-14)
To the same approximation as the Bousslnesq approximation, we may replace (p)^ , the density at
the surface, by p , the density at depth z , in the second term on the right side of (2.14).
Hence
x.' (2.15)
x, l
Designating the Xj -component of the slope of the pressure surface at depth x3
relative to the slope of the water surface by i » tnen
The x, -component of the ensemble-averaged equation of motion can then be written
Likewise, the Xj-component of the ensemble- aver aged equation of motion becomes
where the pressure force term i ||j- has been replaced by
(2.18)
(2-19)
The relative isobaric slope terms i1>p>T, "nd ia.p.n are given by
i . I f^ SP d^ •
and <2'20>
2JP,1! P Jj. ^*3
The density depends on temperature, salinity and pressure. In the estuary, the
pressure dependence is unimportant. The density dependence on temperature and salinity is
12
-------
conveniently expressed in tables issued by the U. S. Navy Oceanographic Office (1966), in terms
of the oceanographic function sigma-t; i.e.
afc = 103 (p-1)
for p expressed in units of the c.g.s. system. Hence
10-3
Xi 3xj
and
90.
- ID'3 T-
Hence, the relative isobaric slope may be expressed as
xa' (2.21)
Xa
at
and
, IP"' ^ !t dx,' (2-22)
1+10- 3at J^ ax, "^
Given the spatial distribution of temperature and salinity, then the spatial distri-
bution of ac can be obtained from tabular compilations such as USNOO (1966) or from appropriate
numerical relationships which have been developed for use in high-speed computers. Equations
(2.21) and (2.22) can then be used to compute the relative isobaric slopes J-1>p>T1 and i2,p,ri
for use in the horizontal components of the equation of motion.
The slope terms i]_ ^ and ±2 ^^ arise because the vertical distance between any
two isobaric surfaces varies inversely as the mean density in the water layer between these
surfaces. These terms, which remain in modified form even in the vertically averaged, two-
dimensional equations and the sectionally averaged, one-dimensional equation, have been
neglected in numerical models which have been developed to date. In a strong, moderately mixed
coastal plane estuary, typical of those which occur along the Atlantic coast of the United
States, these terms are not negligible. Consider for example the James River estuary. The
mean density of the vertical water column at the mouth of the estuary is approximately
1.025 gm cm'3, while 30 miles up the estuary the mean density is approximately 1.000 gm cnr3.
The vertical distance between the water surface and the 10 decibar pressure surface (which
occurs at a depth of approximately 10 meters) is therefore about 25 cm (0.82 ft) greater at the
upstream location than at the mouth. The mean longitudinal slope of the 10 decibar surface
relative to the water surface over this 30 mile distance is therefore 5.5 x 10"«. Variations
in this slope with time over a tidal cycle are relatively small. In this same stretch of the
estuary, the longitudinal slope of the water surface g- varies with tidal phase from zero to
maximum values of about + 1.5 x 10"B .
These relative slopes of the pressure surfaces resulting from horizontal variations
in density, which are significant in estuaries having relatively strong horizontal salinity
13
-------
gradients, are the primary cause for the characteristic vertical shear in the horizontal current.
In such estuaries the seaward-directed flow during the ebb period typically decreases with depth,
while the flow directed up the estuary during the flood period typically increases with depth.
Consequently there is a net non-tidal circulation with seaward-directed flow in the upper layers
of the estuary and a flow directed up the estuary in the lower layers. The intensity of this
vertical shear in the current pattern is not necessarily related to the vertical gradient in
salinity, as has been pointed out by Hansen and Rattray (1966). That is, vertical mixing may be
strong enough in some estuaries to produce near vertical homogeneity in salinity, but the
horizontal gradients in salinity, and hence density, may still be sufficient to produce signifi-
cant vertical shear in the current. Such a vertical shear can be very important in producing
longitudinal dispersion of an introduced pollutant.
While, in this treatment, the vertical component of the equation of motion has been
reduced to the hydrostatic equation, the vertical velocity remains in certain terms in the
horizontal components of the equation; i.e., in the field acceleration terms and in the
Reynolds stress terms. Thus it is certainly not a priori evident that the terms u3 dUj/axj ,
and afc/a*3 arc negligibly small.
In order to put Equations (2.17) and (2.18) into a form which offers any hope of
solution, it is necessary to relate the Reynolds stress terms
(2.26)
dX,
. du2\ 3 / /8u2 8u3>\
I2v2,2 a3Ej; + 5xJ (V2,3 l5xj + al^y
14
-------
These equations have sometimes been further simplified by: (a) assuming
v = v = v, = v, , - VH , the horizontal "Austausch", or horizontal eddy viscosity;
1,1 1»2 i , 1 t-i*- "
(b) assuming v, , = v2 3 " v3 , the vertical eddy viscosity; and (c) neglecting Su^dx-j
and 3u3/8x2 in'comparison with 3u2/Sx3 . The horizontal components of the ensemble mean
equation of motion then become
a , ,u^ 2Yl
a^ IVH Va^ + axJ/J
(2.27)
au.,
and
(2.28)
d Xn I H dXnJ OXo v. J OXn*
There does not exist any clear direct observational or theoretical evidence for
further simplifying these equations, even though some indirect evidence will be used later in
this section in further modifying them. It should also be noted that the coefficients of eddy
viscosity are functions of position and, perhaps, time. They are dependent upon the intensity
of the turbulent velocity fluctuations, and hence probably on the velocity, the boundary
roughness and distance from the boundary, but the relationships are not known. The vertical
eddy viscosity must also be dependent upon the vertical stability, or more probably the
Richardson Number.
In any case, Equations (2.27) and (2.28), together with Equations (2.21) and (2.22)
and the equation of continuity (2.8), are the equations which one must start with in developing
a three-dimensional dynamic model of an estuary. If the goal is to develop a dynamic model
capable of computing the time dependent surface elevation n and the three components of the
time-dependent velocity ut , these equations, standing alone, in the form given, are insoluble,
since there are more unknowns than equations. However,-there are several possible procedures
for overcoming the difficulty.
In Section 2 of this chapter a two-dimensional vertical-averaged model is
developed. In that model the surface elevation TI appears as a dependent variable. Assuming
for the moment that such a two-dimensional model can be shown to provide satisfactory estimates
of the time and space variations in TI , even for estuaries with significant vertical and
horizontal variations in salinity (and hence density) and with vertical shear in the velocity
field, then such a two-dimensional model could be used to determine TI , as well as the vertical
mean values of the two horizontal velocity components. Thus, given TI as a function of time
and position, Equations (2.8), (2.27), and (2.28), together with the auxiliary equations (2.21)
and (2.22) for determining the relative slope terms, could theoretically be solved, by appro-
priate numerical means, to obtain the time-dependent three-dimensional velocity field.
In order to make such an approach possible, however, it will be necessary to develop
empirical and/or theoretical relationships for the coefficients of eddy viscosity in terms of
known or computable parameters.
15
-------
It should also be pointed out that the slope terms 1^ and i2 r depend upon
the spatial distribution of density and hence of salinity. Consequently the dynamic model
cannot be used to compute the field of motion, which then would be used in the determination
of the distribution of salinity. The equations of salt continuity, to be developed in the next
sub-section, must be solved simultaneously with the dynamic equations.
1.2.3 Conservation of Dissolved Constituents
The differential equation expressing the conservation of salt, in three-dimensions,
for instantaneous values of the salt concentration and the velocity, may be written
__
(2.29)
where D is the molecular diffusivity. Proceeding as in the cases for conservation of mass
and of momentum, the ensemble average salt balance equation is given by
(2.30)
Note that, using the equation of continuity Su^ax^ - 0 , the first term on the right side of
(2.30) may also be written
a , , as
dx. i i dx*
1. 1.
In order to make Equation (2.30) tractable, it is usual to replace the non-advective flux terms
by Fickian-type diffusion terms. That is,
where K. - K, for * = i and otherwise is zero. Then Equation (2.31) is written
f| .._|_(Ul8) ^Kf^l^l^^^ (2.32)
where the molecular diffusion term has been neglected compared to the eddy diffusion terms,
and it has been assumed that Kj - K2 - Kg , the horizontal eddy diffusivity.
As in the case of the coefficients of eddy viscosity, there does not exist an adequate
theoretical or empirical basis for relating the coefficients of eddy diffusivity to the ensemble
mean velocity distribution and to the distribution of density.
However, some limited evidence does exist which can be used to at least indicate
further simplifications of the three-dimensional conservation equations, and also a possible
format for relating the eddy coefficients to known or computable parameters. Before proceeding
with such a further treatment, however, consider the general form of the equation expressing
the balance for a dissolved constituent other than salinity, for which there may be internal
sources and sinks .
16
-------
In Equation (2.32), It has been assumed that salt is conserved within the water body;
that is, that there are no internal sources and sinks for salt. Sources and sinks may exist at
the boundaries, however. For example, the vertical diffusion term takes on a boundary value at
the water surface equal to the effective flux of salt across the boundary resulting from the
difference between evaporation and precipitation. In the case of a dissolved substance which
is not conserved within the water body, such as an oxidizable pollutant or a radioactive
material subject to decay, additional terms must appear in the equation.
Designating the ensemble-averaged concentration of a non-conservative material by C.
then the three-dimensional equation expressing the local time rate of change of concentration
is given by
H - - C-iO - k + » fr + 'g - *d <2'33)
where r is the mass rate of generation of substance per unit volume and rd is the mass
rate of decay of substance per unit volume. Again, expressing the non-advective flux k
in terms of an eddy diffusivity and incorporating the molecular diffusion term into this
turbulent diffusion term, Equation (2.33) takes the form
(2-34)
1.2.4 Some Possible Further Simplification of the Three-Dimensional Conservation Equations
for Mass. Momentum and Salt
Studies of one-dimensional and two-dimensional models of conservation of mass,
momentum and salt, and of one-dimensional and two-dimensional water quality models, suggest
the following:
(a) The greater the detail provided by the model with respect to variations
in velocity in time and space, the less significant are the horizontal
eddy diffusion terms, and probably the horizontal eddy viscous terms.
(b) In the vertically averaged two-dimensional equations of motion, the
bottom boundary frictional term is satisfactorily represented by a term
of the form
u(u" + v»)^
8
V
where C, is the Chezy coefficient.
(c) In a study of the James River estuary using a tidal mean, lateral-
averaged two-dimensional equation of motion, Pritchard (1956) showed that
the field acceleration term u"3 3^/8X3 is small compared to other terms
in the equation.
These observations suggest that the three-dimensional equations expressing conservation
of mass, momentum, and salt might be further simplified as follows:
17
-------
(a) In the Xj-component of the equation of motion, Equation (2.27), neglect
the field acceleration term u3 au^/axj and the horizontal eddy viscous
terms 3/BXj^ [2vH au^axj) and 3/3x2 [vjjCSUj/axj + auj/dx^ )
(b) In the x^-component of the equation of motion, Equation (2.28), neglect
the field acceleration term u3 3u2/3x3 , and the horizontal eddy viscous
terms 3/dXj [vH(au2/3x1 + au^axj) } and 3/8x2 {2vH 3u2/ax2}
(c) In the salt balance equation, Equation (2.32), neglect the horizontal
eddy diffusion terms 3/BXj^ {Kg ds/ax^ and 8/3x2 {Kg as/3x2]
The x,-component of the equation of motion then becomes
the xj-component of the equation of motion then becomes:
— . . _2 *"~ 0 a „ -2 _ _ -!! -. -i ' !"". (2-36)
and the salt balance equation becomes
as. a_ / )
at " ax^ v i*'
The auxiliary Equations (2.21) and (2.22) must still be used to compute the relative isobaric
slopes i, and 1- from the horizontal and vertical distribution of density (salinity)
Also, the equation of continuity in the form given in Equation (2.8), that is
au.
must be also solved simultaneously with the dynamic equations and the salt balance equation.
The solutions to these equations must satisfy the following vertical integral
conditions:
.T) ,au, du, au, i
I \ -^-^ + u, ^-^ + u_ T—=• + g ii „ „ - fu, r dx,
J , lat l 3x» 2 ax« I.P.TI *J J
"5 . t (2.38)
3u, i
axf + g ^.P.TI * ^ **3
t (2-39)
18
-------
"2 J-?
and
r-Tl
sir f. "i <«3 + sir f, »2 ** - -1? (2-40)
" "
ft r11 A r71 , a I"
^p J s dx- + jj- J (uj^s) dx.j + jj£— J
J_? " ""1 "-5 - -„
where -5 is the x3 value at the bottom, (E-P) is the difference between evaporation and
precipitation at the water surface, expressed in units of length per unit time, TW>I and TW> ,
are the x, and x- components of the wind stress at the surface, and the subscript h
implies a vertical average. That is,
Uj_ h - ^j- I"1 ^ dx3 (2.42)
1 r^ dx (2 44)
J s 3
The solutions to the three-dimensional equations must also satisfy the following
cross-section integral conditions, where in order to simplify the notation it is assumed that
the coordinate system is oriented such that . x.-^ is directed along the axis of the estuary,
and the cross section at any point x, has an area o and is in the plane perpendicular to
SI JJ
- !f
and
(E-P)b
(2.46)
For the coastal plane or drowned river valley estuaries characteristic of the Atlantic
Coast of the United States, which in their upper reaches extend into the fresh water tidal
river, and are bounded at their upper end by the position where the rise and fall of the tide
and the tidal oscillations in the current are reduced to zero, an additional integral condition
can be obtained. If Equations (2.45) and (2.46) are integrated from x.,^ = (XI)Q where the
integrated flow through the cross section equals the river flow QR to a cross section in
the estuary proper located at x.^ , then
(2.47)
J J J. ~R O i. W /_ \ *->
and
_-x (2^g)
19
-------
Equations (2.8), (2.35), (2.36) and (2.37), expressing the conservation of mass,
momentum and salt in three spatial dimensions, plus Equations (2.38), (2.39), (2.40), (2.41),
(2.45), (2.46), (2.47) and (2.48) expressing the integral boundary conditions, together with
the auxiliary Equations (2.21) and (2.22), represent a set of equations which, on the basis of
present knowledge, include the most important processes controlling the dynamics and kinematics
of an estuary. As such, these equations offer a suitable basis for a three-dimensional,
numerical model of an estuary.
The dependent variables to be computed from such a model are: (a) the surface
elevation of the estuary T) as a function of horizontal position (x^, x2) and of time t ;
(b) the three velocity components, u^, u2 and u3, as a function of horizontal position and
depth (x.p Xj, x3) and of time t ; and (c) the salinity s as a function of horizontal
position and depth (xj, Xg, x3) and of time t. The necessary inputs to the model are:
(a) the physical dimensions of the estuary, i.e., the horizontal boundaries and the distribution
of depth at, say, mean low water; (b) the temporal and spatial distribution of atmospheric
pressure and of surface wind stress over the period covered by the computation; (c) values of
the vertical coefficients of eddy viscosity Vj and eddy diffusivity K3 as a function of
position (x^, Xg, x^) and time during the period covered by the computations; (d) values of
the water density p as a function of position (x^, x2, x3) and time during the period
covered by the computations; (e) the fresh water inflow to the estuary as a function of time
during the period covered by the computations; (f) values of the difference, evaporation minus
precipitation, at the water surface of the estuary, as a function of horizontal position
(x^ x2) and time during the period covered by the computations; (g) the elevation of the water
surface r\ along the line marking the seaward boundary of the estuary, as a function of time
during the period covered by the computations; (h) values of the salinity s in the cross-
section marking the seaward boundary of the estuary, as a function of time during the period
covered by the computations; and (i) an initial set of values of the dependent variables at
all positions In the estuary.
The spatial and temporal distribution of density is required for computations of the
relative isobaric slopes ij and 12 from Equations (2.21) and (2.22). The density
of the estuarine water depends on temperature and salinity. In estuaries in which the relative
isobaric slopes contribute significantly to the dynamic balance, salinity variations are
generally much more important than temperature variations in determining density variations.
Consequently it should be adequate to assume either a constant water temperature throughout
the estuary, or a time independent spatial distribution of temperature characteristic of the
period covered by the computations. The initial distribution of salinity would then be used,
together with the assumed temperature, to determine the initial spatial distribution of density .
Equations (2.21) and (2.22) would then be used to compute initial values of the relative
isobaric slopes l1>p>T) and I2,p,r\ as a tv*0**-0* of position (xj, x2, x3) . After each time
step in the computation, the new distribution of salinity would be used with the assumed
distribution of temperature to determine the spatial distribution of density, from which new
values of ij and 12 would be found for use in the next time step.
The major remaining unknown input terms for the proposed model are values of the
vertical eddy viscosity v3 and the vertical eddy diffusivity IU. Some evidence is available
as to the relationship of K3 to the geometry of the estuary, the velocity shear and the
vertical stability. Presumably the eddy viscosity would depend on these parameters in a
similar manner. However, research is required to develop generally applicable relationships
such that VQ and K3 can be determined from known and/or computable parameters.
20
-------
Pritchard (1960) using the results of studies made in the James River (Pritchard
1954, Kent and Pritchard 1959) found that for tidal mean, lateral -aver aged data, the vertical
eddy diffusivity was related to the mean tidal velocity |u| , the water depth h , the
height H , period T , and wave length L of the surface wind waves, and the Richardson
Number R by the equation
(2.49)
where Z is the vertical distance below the water surface (i.e., Z - TI - x3 ) , h is the
vertical distance between the water surface and the bottom, and TI, ?, and B are non-
dimensional constants. For the James River estuary, the numerical value of these coefficients
were: TI - 8.59 x 10"3; 5 - 9.57 x 10~3; and p - 0.276. This equation is based on data
collected during one relatively short time interval in one segment of one estuary. It is
based on tidal mean data. The degree to which these results can be extended to other estuaries
and to conditions of time-varying currents within the tidal cycle is not known.
For an estuary which is vertically homogeneous and for calm wind conditions,
Equation (2.49) reduces to
K3(Z) - 9.59 x lO'3 |u| z2
-------
2. TWO-DIMENSIONAL MODELS
D. W. Pritchard
2.1 INTRODUCTION
The basic time-dependent equations expressing the conservation of mass, momentum, and
of a dissolved constituent of the waters of an estuary, In three spatial dimensions, as
developed in Section 1 of this chapter, will serve as the starting point for the treatment of
two-dimensional models of an estuary to be presented in this section. This section will include
discussions of the significance of terms which are neglected and of simplifying assumptions
which are made in order to develop a model composed of a tractable set of equations. Where
adequate descriptions are already available in scientific publications or technical reports,
existing models will be referred to, but not described in detail.
2.1.1 The Coordinate System and the notations to be Used
In Section 1 of this chapter tensor notation was employed in conjunction with a
right-handed rectilinear coordinate system, x^ , with the x^ and x2 axes lying In a
horizontal plane coincident with the plane of mean tide level, and the x^-axis directed
upwards, positive in the direction opposite to the acceleration of gravity. In this section
it will generally be more convenient to employ an x-y-z rectilinear coordinate system, and
a vector-scalar notation system. The x- and y-axes are considered to lie In the horizontal
plane coincident with the plane of mean tide level, and the z-axls is directed upwards. Thus
the x-axis corresponds to the x^-axis; the y-axls corresponds to the Xj-axls; and the
z-axls corresponds to the x^-axis. The ensemble average velocity components along the x-,
y-, and z-axls will be designated by u, v, and w, respectively. Thus u in the
notation to be employed in this section corresponds to u. In the notation used in Section 1.
Likewise, v corresponds to u2 , and w to u3 .
Much of title remaining notation used in Section 1 can be carried over directly to this
section, where confusion might result, special note is made in the text.
2.2 THE BASIC BQUATIONS FOR A VERTICALLY AVERAGED TWO-DIMENSIONAL MODEL
2.2.1 Conservation of Mass
We start with Equation (2.6), the ensemble mean equation of mass continuity In three
spatial dimensions as developed in Section 1. In the notation used In this section, that
equation is written
22
-------
ajE. . a (up) . a(vp) . a (wo) , 0 (2.52)
at ax ay az
Now consider that, at any point (x,y,z), the value of the density and of each of the velocity
components is equal to the sum of the vertical mean value of that variable taken over the water
depth plus a deviation term representing the difference between the value of that variable at
the point (x, y, z) and the vertical mean value of that variable. That is, we designate
P ' Ph + Ph'
u - uh+uh'
(2.53)
w - wh + wh
where
p dz
(2.54)
v dz - h- h • h - h-° <2-56)
Substituting from Equation (2.53) into Equation (2.52) and taking the integral of that equation
over the total depth h - (?-Hi) gives
'> + > " ° <2-57)
23
-------
Note that even though h and
, a(hvh> _ an (2.58)
ph at ax ay at
As discussed in Section 1, the first term in (2.58) may be neglected for a fluid such
as water which has a very small compressibility, and for which the Boussinesq approximation is
appropriate. Consequently, the equation of continuity of mass in two-dimensions is reduced to
.M (2.59)
ax
2.2.2 Conservation of Momentum
The three-dimensional ensemble-averaged equations of motion in the form given by
Equations (2.17) and (2.18) of Section 1 of this chapter serve as the starting point for the
development of a vertically averaged, two-dimensional model of the dynamics of an estuary. In
the notation u»ed in this section, these equations are (omitting the molecular viscous term)
ft
(2.60)
and
If - !SF + If + W - -• 0 -• S.P.I. - i ' fa - £ ^ - & k • k
(2.61)
Mow, we introduce into these equations
u = "h + "h"
v = vh + V
w - wh + V
24
-------
and integrate the equations from the bottom (z = -?) to the surface (z = TI) . The resulting
equations are
(2.62)
- E $7 ih< W + k>h} + Eb Tw,y - K^ TB,x
where T and T are the x- and y-components of the surface wind stress, while TB>X
and TB W'Xare the x- and y-components of the bottom stress. This last term in each equation
has usually been replaced by a term of the form
1 T , (2.64)
TB,x 8 ET^»
and
where C^ is the Chezy coefficient.
These forms of the equation of motion are similar to the equations used by Leendertse
(1967) and by Masch et a^. (1969), among others, in formulating a two-dimensional dynamic
model of an estuary, but with two rather significant differences.
Equations (2.62) and (2.63) each include a term which arises from the relative slope
of the pressure surfaces due to the horizontal variation in density (salinity) . The term
is the y-component of this same slope. The significance of this term in the drowned
river'valley estuaries of the Atlantic Coast was discussed for the case of the three-dimensional
equation in Section 1 of this chapter. At this tine we might note that in the James River
estuary, for example, the vertical-averaged value of the longitudinal component of the relative
isobaric slope has been observed to be about 2.7 x 10"6, and relatively independent of time
over a tidal cycle. This amounts to about 36% of the mean magnitude of the slope of the water
surface, averaged over a tidal cycle. That is, in this estuary the longitudinal slope of the
water surface dVax varies with tidal phase from zero to maximum values of about + 1.5 x 10 .
Each of these equations also contains terms of the type k*h and
(u^'v ' + . >h . The parts of these terms representing ensemble averages of the turbulent
velocity fluctuations, (v">h , n and h • "present the vertical averages of
the cross-products of the velocity departures at any depth from the vertical average velocity.
25
-------
If the velocity field has a vertical shear which is predominantly of one sign, then these terms
can be of significant size. In the typical drowned river valley estuary, the vertical shear
of the horizontal velocity is predominantly of one sign; i.e., if x is oriented along the
axis of the estuary toward the sea, then 3u/8z is predominantly positive.
It can be argued that terms of the type h can be replaced by the product of
a coefficient, having the dimensions of kinematic viscosity, times an appropriate component of
the deformation tensor, as was done with the ensemble average terms such as k . That is,
we might assume that
Su,,
" "
and
However, these coefficients certainly are not the same as the coefficients of eddy viscosity
which were employed in Section 1 of this chapter in the terms replacing the Reynolds stresses
in the three-dimensional equations.
For the present we have no basis for relating these terms to known or computable
parameters. Our best hope for the present is that, as they appear in the dynamic equations,
the terms containing the vertical mean cross-products of the velocity deviations can be
neglected. If this is so, then (2.62) and (2.63) become
(2.66)
T
Tw,
and
(2.67)
Tw,y
Equations (2.21) and (2.22) from Section 1 of this chapter must be used as auxiliary
equations to the above forms of the equation of motion in order that the relative isobaric
slopes, i^ and i , be computed from the distribution of density (salinity). Since
the temporaiPchanges in* thVdi attribution of salinity will alter the values of the relative
isobaric slope, the equation expressing conservation of salt in two-dimensions, to be developed
in the next sub-section, must be solved simultaneously with the equations of motion and the
two-dinensional form of the equation of continuity of mass.
26
-------
2.2.3 The Vertically Averaged Two-Dimensional Salt Balance Equation
Equation (2.30) given in Section 1 of this chapter, expressing the conservation of
salt in three dimensions, serves as a starting point for this development of the vertically
averaged salt balance equation. Omitting molecular diffusion, that equation is written in the
notation used in this section as
If ' ' & (U8) ' aV (vs> - £ (W8) - & >k - sy k - & k <2-68>
As was done in the treatment of continuity of mass and salt, we express the dependent variables
as the sum of a vertical average plus a difference term. For example
s = sh + sh'
where
sh = E J_5 s dx3 ' h
Integrating Equation (2.68) then gives
3 (hs,,) , -, --
-TF- - - & - 57 0«Vh> - & thk>h3
(2.69)
-aV^V-h' + k>h} + *h (E-P)
Again, this equation contains terms of the type k>h--Kx*iir
(2.70)
k>h--ViF
Then Equation (2.69) becomes
(2.71)
37 K* IF) + ^
27
-------
Coefficients such as K^* and K* , which arise from terms involving shear in the velocity
field and spatial gradients in the concentration distribution, are sometimes called "effective"
diffusion coefficients. It is to be expected, however, that these coefficients will differ
considerably from the eddy coefficients of diffusion which were introduced in the three-
dimensional salt balance equations. There is also less justification for neglecting the
diffusion terms in Equation (2.71) than for neglecting the horizontal viscous terms in the
two-dimensional momentum equations. It is reasonable to assume, however, that Kj^* - Ky* = KJJ*
the "effective" horizontal diffusivity, and Equation (2.71) would then be written
|y ChVfc) +
The vertically averaged equation expressing the conservation of a dissolved or finely
divided constituent, having sources and sinks in addition to those for salt, would then be
written
(2.73)
+ ^ K* ajr) + Ch + hrg>h - hrd>h
where r h is the vertically averaged mass rate of generation of the constituent per unit
volume (such as the discharge of a waste material) and rd h is the vertically averaged mass
rate of decay of the constituent per unit volume (such as oxidation of an organic waste).
2.2.4 Sumnary of the Set of Vertically Averaged Dynamic and Kinematic Equations Suitable to
Serve as the Basis of a Two-Dimensional Numerical Model of an Estuary
The vertically averaged equation of mass continuity,
<2-59>
together with the somewhat simplified vertically averaged equations of motion, given by
(2.66 j
and
28
-------
(2.67)
T B
Tw,y g
represent a set of three equations which could be solved simultaneously for the three dependent
variables uh, vh and TI , provided the necessary numerical values of the initial and boundary
conditions are available. However, we note that the vertical averages of the relative isobar ic
slopes, h and h , depend on the distribution of density and hence of salinity.
Thus the vertically averaged salt balance equation,
s (E-P) (2.72)
must be solved simultaneously with Equations (2.59), (2.66) and (2.67) to provide the horizontal
distribution of salinity. As discussed in Section 1 of this chapter, the horizontal distribu-
tion of salinity thus determined can be combined with an assumed distribution of temperature to
compute the horizontal distribution of density. The vertical mean values of the relative
isobaric slope terms Ux >h and p>T1>h can then be determined from the relationships
(obtained by taking the vertical average of'Equation 2.20 from Section 1 of this chapter)
and
The required inputs for numerical solution of Equations (2.59), (2.66), (2.67) and
(2.72) are: (a) the physical dimensions of the estuary; (b) the temporal and spatial distribu-
tion of atmospheric pressure and of surface wind stress over the period covered by the compu-
tations; (c) the fresh water inflow to the estuary as a function of time during the period
covered by the computations; (d) values of the difference, evaporation minus precipitation, at
the water surface of the estuary as a function of horizontal position and time during the
period covered by the computations; (e) an assumed horizontal distribution of vertical-averaged
water temperature in the estuary; (f) values of the Chezy resistance coefficient as a function
of position in the estuary; (g) values of the effective horizontal diffusivity as a function
of position in the estuary and of time over the period covered by the computations; (h) values
of the surface elevation TI and of the vertically averaged salinity sh along a line marking
the seaward boundary of the estuary, as a function of time during the period covered by the
computations; (i) an initial set of values of the dependent variables at all positions in the
estuary.
The dependent variables to be computed from the numerical solution of these equations
are: (a) the surface elevation T, ; (b) the vertically averaged horizontal components of the
velocity, uh and v. ; and (c) the vertically averaged salinity sh ; all as functions of
horizontal position (x, y) and of time over the period of the computations. Equations (2.74)
and (2.75) must be used as auxiliary equations of the computations to obtain values of the
29
-------
vertically averaged relative isobaric slopes h and h from the comPuted
distribution of salinity and the assumed distribution of temperature. The additional integral
condition, obtained from the Integration of the vertically averaged equation of continuity across
the estuary and then along the estuary from the landward boundary of the estuary (x - XQ) ,
where the time variation In water surface elevation is zero, must also be satisfied. That is,
<"*') (2-76)
where J dy indicates an integral across the width of the estuary b(x,t) ; and o(x,t) is
the cross-sectional area of the estuary perpendicular to the x-axis, which is considered here
to be oriented along the estuary.
2.2.5 An Alternate Development of a Set of Vertically Averaged Dynamic and Kinematic
Equations in an Estuary
Jan J. Leendertse (personal communication) has suggested an alternate approach to the
development of the vertically averaged equations which offers some computational advantages
over the development given above. Briefly, this approach Is based on the assumption that the
ensemble mean velocity components, u(x,y,z,t) and v(x,y,z,t) , and salinity s(x,y,
z,t) are related to their respective vertical-averaged values u^x.y.t) , vh(x,y,t) ,
and Sj^x.y.t) by the relationships
u - Uj, 11 + fu(z)]
v-vh(l + fv(z)} (2.77)
and
s - sh {1 + fg(z)}
where f (z) , f (z) and fs(z) are functions of the vertical coordinate z , and also perhaps
of the tine t , and describe all the vertical variations in u , v , and x , respectively.
Hote that this manner of relating the ensemble mean values of the variables to their respective
vertical- averaged values is more restrictive than the division given in Equation (2.53), since
Equation (2.77) requires that u, v and s are zero everywhere in the vertical when their
respective vertical mean values u^ vfa and sh are zero, a condition not generally true.
The vertical integrals of the functions fu(z), f^z) and fg(z) must all be zero.
However, cross-products involving these functions do not necessarily have vertical averages
equal to zero. The integral cross-products which appear in the vertically averaged equations
are
{1 + f..(z)l dz
-5 " "-5
11 {1 + fu2(z)} dz
30
(2.78)
-------
-5
and
I"*1 uv dz - u^ f1 {1 + fu(z)} {1 + fv(z)} dz
-? ' ~? (2.79)
us dz - u^ f {1 + fu(z)l {1 + fs5 ll+fs(z)1 dz
-? -5 (2.81)
Now designate
f U + fu2(z)]dz
1 !•*! ,.,^11 O 9>k\
= | [1 + f (z)*f (z) } dz v'-o'v
" -?
and
_ I T11 fi 4- f (,-\-f (z\'\ dz (2.85)
and a are called the distribution functions for the respective
ss-prctscne the^ HopefuUy, the spatial, and perhaps te.poral variations in
these parameters are sufficiently regular to be revealed from observations of the vertical
variation to the velocity and salinity in the estuary. The uncertainties in these parameters
substitute for the uncertainties in the effective horizontal coefficients of viscosity and
dlffuslvity. Hopefully, however, the variations In these in space and time would be more
readily related to known parameters than in the case of the effective coefficients.
In any case, the vertically averaged equation of continuity has the same form in
this development as that given earlier. That is, Equation (2.59) remains applicable. H>e
vertically averaged equations of motion become, however,
. _gh u _gh p^>h . £ gl
31
-------
and
ay 6 ay & v y.p.Vh Ph ay
(2.87)
where the vertical average of the ensemble mean cross-product terms Involving the turbulent
velocity deviations, such as k , have been neglected. There is considerably greater
justification for neglecting these terms in the development of Equations (2.86) and (2.87)
than for neglecting terms such as ("h'1^' + (u'v'^k^h as was ^0ne *n arr*ving at Equations
(2.66) and (2.67). Essentially, the distribution functions auu, a^ and auv take the
place of cross-products of the type ("h'Vh'^h ' Ihus Etluat*ons (2.86) and (2.87) are at least
conceptually more complete than Equations (2.66) and (2.67). The somewhat restrictive relation-
ship between u and u, , etc., as given in Equation (2.77) must, however, be kept in mind in
making this comparison.
The vertically averaged salt balance equation, using the relationships given in
Equation (2.77), then becomes
]
(2.88)
+ w ih *H w + -h
where the diffusion terms arise from the substitution
«u's'>k>h = -gHinr
and
The vertically averaged horizontal dif fusivity £g is thus probably much smaller than the
effective horizontal dif fusivity Kg* which appears in Equation (2.72), since the latter
coefficient arose from a substitution of the form given in Equation (2.70).
This alternate approach provides a model less dependent on our lack of knowledge of
the effective horizontal coefficients of viscosity and dif fusivity. However, in place of the
need for finding relationships between these coefficients and known parameters, it is necessary
to define the spatial and temporal variations in the distribution functions auu> aw. QUv'
aus, and avg. Certainly the observational requirements for establishing empirical relation-
ships for the distribution functions are less severe than in the case of the effective
horizontal coefficients of viscosity and dif fusivity. On the basis of available evidence, the
distribution functions will, for most estuaries, and for most times during the tidal cycle,
range between 0.5 and 1.5. The more nearly the estuary is vertically homogeneous in velocity
and salinity, the closer the distribution functions will be to unity.
32
-------
REFERENCES FOR SECTIONS 1 AND 2
Bowden, K. F., 1960: Circulation and mixing in the Mersey Estuary. Publication No. 51,
International Association of Scientific Hydrology, Commission of Surface Waters.
Hansen, Donald V., and Maurice Rattray, Jr., 1966: New dimensions in estuary classification.
Limn, and Oceanography, 11, No. 3 (July).
Kent, R. E., and D. W. Pritchard, 1959: A test of mixing length theories in a coastal plain
estuary. J. Mar. Res., 18, No. 1 (June).
Leendertse, Jan J., 1967: Aspects of a computational model for long-period water-wave
propagation. Rand Corporation Memorandum RM-5294-PR, May, 1967.
Masch, Frank D., 1969: A numerical model for the simulation of tidal hydrodynamics in shallow
irregular estuaries. Technical Report HYD 12-6901, Hydraulic Engineering Laboratory,
Department of Civil Engineering, The University of Texas at Austin.
Pritchard, D. W., 1954: A study of the salt balance in a coastal plain estuary. J. Mar. Res. ,
JL3, No. 1.
Pritchard, D. W., 1956: The dynamic structure of a coastal plain estuary. J. Mar. Res., 15,
No. 1.
Pritchard, D. W., 1960: The movement and mixing of contaminants in tidal estuaries.
Proceedings, 1st Conference on Waste Disposal in the Marine Environment, E. A. Pearson
(Ed.). Pergamon Press.
Sverdrup, H. U., Martin W. Johnson, and Richard H. Fleming, 1946: The Oceans, Their Physics,
Chemistry, and General Biology. New York, Prentice-Hall.
U. S. Naval Oceanographic Office, 1966: Handbook of Oceanographic Tables (Compiled by
E. L. Bialek). Special Publication SP-68, Washington, D. C. (Replaces H. 0. Publ.
No. 614.)
33
-------
3. ONE-DIMENSIONAL MODELS
Donald R. F. Harleman
3.1 INTRODUCTION
One-dimensional equations for the hydrodynamic and mass transfer processes in estuaries
have an obvious advantage of mathematical tractability in comparison with their multi-dimensional
counterparts. In addition, the one-dimensional equations utilize and predict information that
is related to available or accessible observational data. This is especially true in dealing
with the mass transfer equations for water quality parameters involving source and sink terms,
such as biochemical oxygen demand and dissolved oxygen.
An important objective of this section is the presentation, comparison and discussion
of the various one-dimensional hydrodynamic and mass transfer models which have been used in
estuaries. Thus it is necessary to consider the most general form of these equations in which
quantities such as concentration and velocity are functions of both longitudinal position and
real time. The general equations are employed as the basis for evaluation of simpler mathemat-
ical models which use non-tidal advective velocities.
The one-dimensional mass transfer equation is derived by a spatial integration of the
three-dimensional equation over the flow cross section. Thus any quantities such as velocity,
concentration or dissolved oxygen are spatial averages over a cross section corresponding to a
specified period of time averaging. The one-dimensional continuity and momentum equations are
derived by considering a material element occupying a cross section of the estuary. The one-
dimensional equations are well suited for estuaries displaying vertical and lateral homogeneity;
however, they have been successfully applied to estuaries with varying degrees of sectional
non-homogeneity. Applications to a fully stratified or saline-wedge estuary should be excluded.
The latter problems can be treated by the theory of two-layer stratified flows (Harleman 1960,
Keulegan 1966).
The various quantities, such as velocity, concentration, cross-sectional area, in the
one-dimensional mathematical models are assumed to be functions of a single spatial variable x
and time t . The longitudinal distance x is measured along the axis of the estuary and
cross-sectioaal areas are normal to the local direction of the axis as shown in Figure 2.1.
One-dimensional mathematical models encompass a wide range of time-averaging concepts.
This has been the subject of confusion in the literature due to the use of ambiguous terms for
time averaging. Adjectives such as "unsteady", "quasi-steady state" or "steady" are meaningless
unless precisely defined with respect to (i) the duration of the time period to be used in
averaging, and (ii) the quantity to which the adjective applies. Careful definition is particu-
larly important in estuarine problems because of the multitude of time-averaging concepts which
have been employed. As shown in Figure 2.2, the term "unsteady" when applied to the section-
mean velocity may imply as many as four possible periods of time averaging. These are defined
as follows:
34
-------
NATURAL BOUNDARY
REFERENCE DATUM
,fw MEAN WATER LEVEL
^^ —a MEAN LOW WATER
X
HYPOTHET 1 CAL CHANNEL
FOR TIDAL FLOW
BOTTOM OF SCHEMATIZED CHANNEL
Fig. 2.1 Longitudinal and transverse schematization of an estuary.
Fig. 2.2 Various periods of time-averaging the tidal velocity.
-------
(i) If it, is of the order of magnitude of a minute, the time average of the velocity
fluctuations associated with the fluid turbulence is obtained. This will be referred to as
the tidal velocity U . Its magnitude varies in real time between zero and the maximum tidal
velocity U , and the direction changes from flood (landward) to ebb (seaward) during each
tidal period!* The time of zero tidal velocity will be referred to as the time of slack tide;
high water slack if the velocity change is from flood to ebb and low water slack for the reverse.
The tidal velocity is always an unsteady quantity.
(ii) If At, is approximately six hours, the average magnitude of the tidal velocity during
either the flood or ebb portion of the tidal period is obtained. This will be referred to as
mean flood (or ebb) velocity. It is inherently unsteady because the magnitudes of successive
mean flood and ebb velocities change from day to day and week to week in accordance with the
spring-neap variations of the lunar month.
(iii) If At, is equal to a tidal period (12.4 hours for a semi-diurnal tide), the non-tidal
advective velocity is obtained. Whenever the time rate of change of estuary volume between the
section under consideration and the head of tide is a harmonic function of tidal period, the
non-tidal advective velocity will be equal to the velocity due to freshwater inflow, Uf
(Pritchard 1958). The non-tidal advective velocity may be steady or unsteady depending on the
time variation in the magnitude of the freshwater inflow.
(iv) If At. is equal to or greater than 25 hours the mean, non-tidal advective velocity for
the specified period of averaging is obtained.
Time-averaging terms may be applied to quantities other than velocity and care should
be taken to specify the intention. For example, consider the instantaneous or pulse injection
of a tracer into a steady, uniform flow. In this case the velocity is steady, however the con-
centration at a fixed point is a function of time and is therefore unsteady.
The various one-dimensional mathematical models employed In analysis of pollution
problems in estuaries can be classified according to the length of the time averaging associated
with the advective velocity. It is essential that good judgement be used in arriving at the
decision of the type of mathematical model to be used. An important function of a state-of-the-
art review is to provide a sound basis for such decisions. Mathematical models for the hydro-
dynamic and mass transfer processes for various time-averaging periods related to velocity are
presented in the following sections. Considerable attention will be given to the development
of models using real time in order to provide a basis for comparison with other models using
longer periods of time averaging.
3.2 MATHEMATICAL MODELS IN REAL TIME
If the events in an estuary are to be treated in real time, the tidal velocity and
concentration quantities In the mathematical models are those associated with measuring instru-
ments having integrating periods which are short compared to the tidal period. We consider an
estuary of arbitrary shape and derive the one-dimensional form of the mass transfer equation for
a non-conservative substance. The resulting differential equation describes the temporal and
spatial distribution of concentration of that substance in an estuary. In addition, the one-
dimensional continuity and momentum equations are derived. The latter pair of differential
equations may be used to determine the tidal velocity which appears in the advective term of the
mass transfer equation.
36
-------
3.2.1 Mass Transfer Equation for a Non-Conservative Substance
The one-dimensional mass transfer equation can be obtained by direct integration of
the three-dimensional equation over a cross section normal to the axis of the estuary. The
three-dimensional equation for turbulent flow, neglecting molecular diffusion is
(2.89)
where u, v, w are the time-averaged velocity components associated with turbulent flow, c
is the local concentration, ex, ey and ez are turbulent diffusivities , and rt is the time
rate of generation (or decay) of substance per unit volume of fluid. We define
u = U + u"
v = v"
w = w" (2-90)
where U - T [ " dA is the real time longitudinal tidal velocity averaged over the cross
A JA
section. The double-primed quantities u", v" and w" are spatial deviations of velocity from
the section-mean value due to the vertical and lateral velocity distribution. Note that
I «"
dA
A
and that v" and w" are not zero even though the section-mean velocities V and W are
zero in a one-dimensional flow. In a similar manner, we define
c = C + c" (2-91)
where C = - \ c dA is the concentration averaged over the cross section and c" is the
A JA
spatial deviation of the concentration due to the variation of concentration over the cross
section. Concentration is defined as the mass of substance per unit mass of solution (or weight
of substance per unit weight of solution) .
Equations (2.90) and (2.91) are introduced into the three-dimensional mass transfer
Equation (2.89), the products of sums are expanded and each term is integrated over the cross-
sectional area A. Inasmuch as the expanded equation contains at least thirteen terms a complete
development will not be presented. Rigorous derivations of this type have been given by Okubo
(1964) and Holley and Harleman (1965). After simplification, the one-dimensional equation may
be written as
(2.92)
The Integral term on the left side involving the cross product of the longitudinal velocity and
concentration deviations represents the mass transport associated with the non-uniform velocity
distribution. This is usually designated as the longitudinal dispersion term. The quantity e,
37
-------
Is the spatial-mean value of the turbulent dlffusivity. The term r^ represents the time rate
of Internal addition of mass per unit volume by generation of the substance within the fluid.
The term r is the time rate of addition of mass per unit volume externally across the lateral
boundaries (sides, bottom and water surface). The two latter quantities are usually known as
source and sink terms. As examples of internal terms, a positive r^ would result from oxygen
produced by suspended algae; and r. would be negative for a substance undergoing radioactive
decay or for the consumption of oxygen by suspended matter possessing BOD. In the case of mass
transfer across external boundaries, a positive value of rg could be due to the absorption
of atmospheric oxygen across the free surface of an oxygen deficient stream. A negative rg
would result from the adsorption of a tracer dye on the solid boundaries.
For steady, uniform flow, Taylor (1954) and Aris (1956) have shown that the advective
mass transport associated with the cross product of the spatial deviations u" and c" can be
represented as an analogous one-dimensional diffusive transport. To distinguish this process
from turbulent diffusion, which is associated with the temporal deviations u1 and c' , the
transport due to spatial deviations is called longitudinal dispersion. On this basis a coeffi-
cient of dispersion E is defined in terms of the concentration gradient as
u"c" dA - - AE (2.93)
The negative sign indicates mass transport in the direction of decreasing concentration. Taylor
(1954) has shown that E is more than two orders of magnitude larger than ?x ; therefore it
is convenient to add the two coefficients and refer to the sum as the longitudinal dispersion
coefficient Z^ , where
EL - E + ¥x (2.94)
Using Equations (2.93) and (2.94), Equation (2.92) becomes
l>+T + T <2-95)
Equation (2.95) is the most general form of the one-dimensional mass transfer for a non-
conservative substance. It is assumed to be valid for unsteady, non-uniform flows within the
restrictions of the one-dimensional approximation. The independent variables, tidal velocity
U , cross-sectional area A , and longitudinal dispersion coefficient E^ can be functions of
both x and t . In addition, the source and sink terms can be functions of x and t .
The dependent variable, concentration, refers to any substance of interest. In general, one
mass transfer equation of the form of Equation (2.95) must be written for each substance under
consideration. Typically, mass transfer equations are written for salinity, dye or radioactive
tracers, BOD and dissolved oxygen.
The following discussions refer specifically to the various terms in the mass trans-
fer equation.
38
-------
3.2.2 Source and Sink Terms
The nature of the internal and external source and sink terms in the mass transfer
Equation (2.95) depends on the substance under consideration. Their specification in quantita-
tive terms is probably one of the most difficult problems in water quality modeling; this is
treated in detail in Chapter III.
A substance is called conservative if the source and sink terms ^ and rg are
both zero. For example, salinity is usually considered to be a conservative substance in an
estuary, in that the source of salinity is specified as a maximum concentration at the ocean
end. The term r£ - 0 since no salinity is generated internally within an estuary and re = 0
if there is no influx of salinity along the external lateral boundaries.
The following discussion will serve to illustrate the proper dimensional form for
two typical source and sink terms. A common sink term of the internal type is a decay which
follows a so-called first-order reaction. This would be appropriate if the substance under con-
sideration were a radioactive tracer. The first-order reaction is stated as follows: the time
rate of decay (or generation) of a substance per unit volume |