PB81-218349
Water Temperature Dynamics in Experimental
Field Channels: Analysis and Modeling
Minnesota Univ., Minneapolis
Prepared for
Environmental Research Lab
Duluth, MN
Feb 81
-------
Page Intentionally Blank
-------
EPA 600/3-81-008
February 1981
PB31-213349
WATER TEMPERATURE DYNAMICS
IN EXPERIMENTAL FIELD CHANNELS: ANALYSIS AND MODELING
by
Heinz G. Stefan,
John Gulliver,
Michael G. Hahn,
and
Alec Y. Fu
University of Minnesota
St. Anthony Falls Hydraulic Laboratory
Mississippi River at 3rd Ave. S. E.
Minneapolis, Minnesota 55414
Grant Nos. R 80368601 and R 80473601
Project Officer
Kenneth E. F. Hokanson
Monticello Ecological Research Station
Box 500, Monticello, Minnesota 55362
ENVIRONMENTAL RESEARCH LABORATORY - DULUTH
OFFICE OF RESEARCH AND DEVELOPMENT
U. S. ENVIRONMENTAL PROTECTION AGENCY
DULUTH, MINNESOTA 55362
-------
NOTICE
THIS DOCUMENT HAS BEEN REPRODUCED
FROM THE BEST COPY FURNISHED US BY
THE SPONSORING AGENCY. ALTHOUGH IT
IS RECOGNIZED THAT CERTAIN PORTIONS
ARE ILLEGIBLE, IT IS BEING RELEASED
IN THE INTEREST OF MAKING AVAILABLE
AS MUCH INFORMATION AS POSSIBLE.
-------
TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing)
1 REPORT NO
EPA-600/3-81-008
3 REfiJfJgNT S ACCESSION NO
21834 9
4,TITLE AND SUBTITLE
Water Temperature Dynamics in Experimental Field
Channels: Analysis and Modeling
5 REPORT DATE
FEBRUARY 1981 ISSUING DATE.
6 PERFORMING ORGANIZATION CODE
7 AUTHOR(S)
Heinz G. Stefan, John Gulliver, Alex Y. Fu, and
Michael G. Hahn
8 PERFORMING ORGANIZATION REPORT NO
9 PERFORMING ORGANIZATION NAME AND ADDRESS
University of Minnesota
St. Anthony FalTs Hydraulic Laboratory-
Mississippi-River _at-3rd Ave._.S.E.
Minneapolis, Minnesota55414
10 PROGRAM ELEMENT NO
11 CONTRACT/GRANT NO
"R8036360-1 and" R8047 3601-
12 SPONSORING AGENCY NAME AND ADDRESS
U.S. Environmental Protection Agency
Environmental Research Laboratory-Duluth
6201 Congdon Boulevard
Duluth, Minnesota 55804
13 TYPE OF REPORT AND PERIOD COVERED
14 SPONSORING AGENCY CODE
EPA-600/03
15 SUPPLEMENTARY NOTES
16 ABSTRACT
This study is on water temperature dynamics in the shallow field channels of the
USEPA Monticello Ecological Research Station (MERS). The hydraulic and temperature
environment in the MERS channels was measured and simulated to provide some back-
ground for several biological studies at the Research Station, including one of the
effects of artificially high water temperatures on fish and invertebrate populations.
An analysis of the temperature measurement problem, the channel temperature regime,
tnicrohabitat conditions, and temporal and spatial water temperature dynamics was
essential for the investigation. The results of this study are also applicable
to the study of temperature dynamics in thermally polluted shallow streams, to the
design of canals for cooling water disposal, and to the study of natural water
temperature regimes of small streams.
A computer simulation program, TSOIL, to predict the unsteady water temperature
distribution in the three-layered system was also developed.
17
KEY WORDS AND DOCUMENT ANALYSIS _
DESCRIPTORS
b IDENTIFIERS/OPEN ENDED TERMS
c COSATI Field/Croup
18 DISTRIBUTION STATEMENT
RELEASE TO PUBLIC
19 SECURITY CLASS jThts Report)
UNCLASSIFIED
21 NO OF PAGES
20 SECURITY CLASS (Thispage)
UNCLASSIFIED
22 PRICE
EPA Form 2220-1 (Rev 4-77) PREVIOUS EDITION is OBSOLETE ,'
-------
DISCLAIMER
This report has been reviewed by the Environmental Research Laboratory-
Duluth, U. S. Environmental Protection Agency, and approved for publication.
Approval does not signify that the contents necessarily reflect the views
and policies of the U. S. Environmental Protection Agency, nor does mention
of trade names or commercial products constitute endorsement or recommenda-
tion for use.
11
-------
PREFACE
Shortly after the Monticello Ecological Research Station (MERS), a field
laboratory of the USEPA Environmental Research Laboratory - Duluth, went into
operation in August 1973 a research program on the effects of elevated water
temperatures, resulting from waste heat discharges on fish and other aquatic
organisms, was initiated by the USEPA staff. The purpose of these studies was
(a) to validate water quality criteria data produced in the laboratory under
semi-natural field (mesoscale) conditions and (b) to identify significant
response by aquatic organisms under field monitoring conditions. The field
channels in which the experimental studies were conducted were supplied with
water, heated artificially by waste heat, from Northern States Power Company's
Monticello Generating Plant.
To provide assistance in these ecological studies, it was necessary to
document, analyze, and predict water temperature characteristics in the field
channels. The study described herein fulfilled that purpose. The study of
water temperature characteristics in the 520 m (1700 ft) long field channels
required investigation of several components such as surface heat exchange
processes , longitudinal dispersion, stratification, and heat transfer to the
cnannel bed. A summary of these studies is provided in this report. Additional
details on the substudies can be found in the following memoranda and theses.
1. Physical Characteristics of the Experimental Field Channels
at the USEPA Ecological Research Station in Monticello,
iMinnesota, by Michael G. Hahn, John S. Gulliver, and H- Stefan,
Memorandum No. 156. St.-Anthony Falls Hydraulic Laboratory,
University of Minnesota, April, 1978.
2. Operational Water Temperature Characteristics in Channel No. 1
of the USEPA Monticello Ecological Research Station, by Michael
Hahn, John Gulliver, and H. Stefan, Memorandum No. 151,
St. Antnony Falls Hydraulic Laboratory, University of Minnesota,
January, 1978.
3. User's Manual for Operational Water Temperature Statistics
Computer Programs WTEMP1 and WTEMP2, by Alec Y. Fu and H. Stefan,
Memorandum No. 162, St. Anthony Falls Hydraulic Laboratory,
University of Minnesota, July, 1979.
-------
4. Water Temperature Data Processing for the Experimental Field
Channels at the USEPA Ecological Research Station in Monticello,
Minnesota, by Alec Y. Fu and H. Stefan, Memorandum No. 167,
St. Anthony Falls Hydraulic Laboratory, University of Minnesota,
April, 1979.
5. Soil Thermal Conductivity and Temperature Prediction in the
Bed of the Experimental Field Channels at the USEPA Ecological
Research Station in Monticello, Minnesota, by J. S. Gulliver
and H. Stefan, Memorandum No. 165, St. Anthony Falls Hydraulic
Laboratory, University of Minnesota, January, 1980.
6. - Pore Water Temperatures and Heat Transfer in a Riffle (Rock)
Section of the Experimental Field Channels at the USEPA
Ecological Research Station at Monticello, Minnesota, by H. Stefan
and J. S. Gulliver, Memorandum No. 166, St. Anthony Falls Hydraulic
Laboratory, University of Minnesota, April, 1980.
7. Analysis of Surface Heat Exchange and Longitudinal Dispersion
in a Narrow Open Field Channel with Application to Water Tempera-
ture Prediction, by John S. Gulliver, M.S. thesis, University
of Minnesota, August, 1977.
8. Experimental Studies of Vertical Mixing in Temperature Stratified
Waters, by Michael G. Hahn, M.S. thesis, University of Minnesota,
March, 1978.
9. Stochastic Water Temperature Characteristics in an Outdoor
Experimental Channel, by Alexander C. Demetracopoulos, Plan B
paper for M.S. degree, University of Minnesota, November, 1976.
The overall geometrical, hydraulic, and soil thermal characteristics of
the Monticello field channels have been described in Item No. 1 above. A
first summary of operational water temperature characteristics in the field
channels was given in Item 2. The surface heat exchange and the longitudinal
dispersion in the artificially heated channels subjected to real (unsteady)
weather conditions were described in Item 7. The water temperature stratifi-
"catrrorr~rn'the"-llpoo±su--of--the--£teld-channel3--w-a3--described--and- analyzed- in- - - -
Item &. - A-computer progrant was developed for the statistical analysis of 3-
"hour water temperatures recorded at several stations in-four-of the-eight
field channels. A user's manual was prepared and made available to the MERS
staff (Item 3). Water temperature data processing was described in Item 4.
Water temperature statistics were computed for the four channels equipped with
water temperature instrumentation for the periods from December 1975 to
September 1977. The tabular computer printout covers many pages and copies
IV
-------
are available for purchase. Earlier statistical analysis of some partial
records led to an unpublished, but fairly extensive, progress report by H.
Stefan, C. Gutschik, A. Demetracopoulos and J. Gulliver entitled "Water
Temperature Data Retrieval and Statistical Analysis for the Experimental
Channels at the USEPA Monticello Field Station," September, 1976, 53 pp.
A spectral analysis of some water temperature records is given in Item 9.
Finally, the heating and cooling of the channel bed by the water above has
been analyzed and a computer program for soil temperature prediction was
described in Item 5. Pore water temperatures in_ .the.riffles.are_described_ in
Item 6.
The work performed for the USEPA Monticello Ecological Research Station
(MERS) under this grant had two basically different purposes. One was to
provide background material and general information in support of the biological
work conducted by the MERS field station staff. The second purpose was to
investigate surface heat transfer, longitudinal dispersion, and water tempera-
ture predictions in streams under field conditions. The second purpose is
basic and includes verification of the applicability of relationships developed
for larger water bodies to smaller streams.
This report has been divided into sections which address the two above
purposes separately. Section 4 - Hydraulic and Thermal Characteristics of the -
MERS Field Channels provide general background information, Section 5, Longi-
tudinal Dispersion, Section 6 - Heat Budget, and Section 7 - Water Temperature
Prediction, give results of basic experimental and analytical studies.
-------
fc ABSTRACT
This report summarizes some of the morphological, hydraulic, thermal
and meteorological characteristics"of the experimental field channels at the
USEPA Monticello Ecological Research Station. It contains an overview of
measurements and parameters which characterized the physical operation condi-
tions of the channels from December 1975 through December 1977.
Reported herein are measured values of hydraulic channel roughness,
permeability and porosity of rock sections, and thermal diffusivity of bed
materials. Recorded water temperatures have been statistically analyzed to
give local diel, longitudinal,and vertical water temperature variations. A
data bank of 3-hour water temperature data was established. Mean, standard
deviation, skewness, and water temperature values at 5 per cent and 95 per cent
probability of occurrence have been calculated on sliding, weekly samples of
these 3-hour data.
Various aspects of the water temperature regime of the experimental
field channels of the USEPA Monticello Ecological Research Station (MERS)
were studied theoretically. Longitudinal dispersion and heat transfer
relationships which were included in a dynamic water temperature model had
to be established. Water temperatures at various channel locations and
Tieteorological parameters were continuously recorded. A convective water
surface heat transfer relationship (wind function) was determined from the
recorded longitudinal temperature-proflie and meteorological parameters during
quasi-steady state periods. Progressive'-heat fronts were-used with tracer
theory to determine longitudinal dispersion coefficients for the field "
channels.
The wind function and the longitudinal dispersion coefficient were
incorporated in an implicit finite difference computer model, MNSTREM, for
the highly dynamic water temperature prediction in the very shallow field
channels. The water temperature model results were verified against
3-hour water temperature measurements over four time periods of up to one
month.
VI
-------
Development of a model for diel water temperature variations at a
one to three hour time scale goes beyond conventional stream water tempera-
ture modeling. The value and necessity of 3-hour to 6-hour weather data in
making accurate predictions has been demonstrated. A standard error of
0.24 C between predictions and measurements could be achieved.
VII
-------
TABLE OF CONTENTS
Preface m
Abstract vi
List of Figures- .~. x
List of Tables T. xm
-Abbreviations arid-Symbols -sn-n. '.~~l 7 —.-. .%-7.-»-» .-.:. .^»-*777. xv
Acknowledgements xviii
1. INTRODUCTION 1
2. CONCLUSIONS 5
3. RECOMMENDATIONS 9
4. HYDRAULIC AND THERMAL CHARACTERISTICS OF THE MERS FIELD
CHANNELS 10
4.1. Channel Geometry 10
4.2. Channel Hydraulics 14
Flow Rates and Residence Times 14
Roughness and Water Surface Slope 14
Flow Velocities 17
4.3. Soil and Rock Characteristics 21
Grain Size Distribution 21
Permeability of Riffle Rock Sections 23
Soil Thermal Diffusivity and Soil Temperature in
Pool Sections 23
Rock/Water Thermal Diffusivity and Temperatures
in Riffle Sections 28
4.4. Observed Operational Water Temperature Regime 31
Water Temperature Instrumentation 31
Water Temperature Data Processing 31
Water Temperature Data Analysis and Results 37
5. A S-TTO*- OF-LONGITODINAL DISPERSION IN THE MERS-FIELD CHANNELS.. 57
":- "".5.1. Background '..." — .- '.-..".";...."... .\ .. .v;-; _.:.'. rr."T7;~T '"5 Tv
---- 5.2. -~Methbd-of Longitudinal Dispersion-Coefficient - - .'.
Determination from Transient Water Temperatures 60
5.3. Formulation of Dimensionless Longitudinal Dispersion
Numbers 64
5.4. Results of Temperature Routing Tests 66
5.5. Effects of Variable Cross Sectional Area on Longitudinal
Dispersion 67
5.6. Comparison with Other Longitudinal Dispersion Measure-
ments 69
5.7. Effect of Roughness (Mannings 'n') on Longitudinal
Dispersion in the MERS Channels 70
viii
-------
6. HEAT TRANSFER ACROSS THE WATER SURFACE IN THE MERS FIELD
CHANNELS 71
6.1. Background 71
6.2. Net Shortwave (Solar) Radiation 72
6.3. Net Longwave Radiation 75
6.4. Net Evaporative and Convective Heat Transfer 76
6.5. Determination of Wind Related Heat Transfer Coefficient
from Measured Steady-State Longitudinal Temperature
Profiles 78
General Procedure 78
Errors in Computed Heat Transfer Coefficients 79
Determination of Wind Function Coefficients 81
7. DYNAMIC WATER TEMPERATURE PREDICTION IN OPEN CHANNELS 91
7.1. General Objective 91
7.2. Formulation of One-Dimensional Finite Difference Equa-
tions for Water Temperature Prediction 91
7.3. Description of Numerical Finite Difference Computer
Program for Water Temperature Prediction 96
7.4. Examples of Water Temperature Predictions in the MERS
Field Channels 97
References 115
Appendices
A. Temperature Front Data for Determination of Longitudinal
Dispersion Coefficient 119
B. Data and Computations of Bulk Surface Heat Transfer Coefficient
and Wind Function Parameters from Steady-State Longitudinal
Water Temperature Profiles 126
C. User's Guide to the Minnesota Stream Water Temperature
Prediction Model (MNSTREM) 133
D. Program Listing for the Minnesota Stream Water Temperature
Model MNSTREM "... 137
S. Sample Input Data for MNSTREM 148
F. Partial Sample Output from MNSTREM 153
G. MERS Weather Station 157
H. Sample Output from Program WTEMP1 Water Temperature Statistics. 162
-------
LIST OF FIGURES
Number Page
1.1 Aerial photograph of the Monticello Station experimental channels,
*" Northern States Power Company Monticello Power plant, and
Mississippi River 2
4.1 MERS channels showing the transition from pool to riffle to pool. H
4.2 Channel cross sections 12
4.3 Cross sections of riffle 5, channel 1 on April 6, 1977 13
4.4 Velocity profiles in Channel 1, riffle 5, looking downstream.
Measured October 30, 1976 18
4.5 Velocity profiles in Channel 4, riffle 11, looking downstream,
with bottom free of weed. Measured May 9, 1977 19
4.6 Velocity profile in Channel 5, Pool 14, looking downstream.
Measured October 30, 1976 20
4.7 Size distribution for material in Channel 1 at MERS 22
4.8 Soil profile below pool and location of thermistor probes 25
4.9 Schematic cross section of one channel pair drawn to approximate
scale 25
4.10 Measured soil temperature profiles in pool 6 during heat wave
test: warming period 26
4.11 Comparison of measured and predicted soil temperatures in Pool 12
...... _ during- heat- wave test , V 27
4.12 Vertical temperature profiles in Riffle 17 during period of
rock temperature study 30
4.13 Thermal probe placement in MERS experimental channels 32
4.14 Seasonal water temperature plots of Channel 1 39
4.15 Seasonal water temperature plots of Channel 5 41
-------
Number Page
4.16 Diurnal water temperature plots of Channel 1 44
4.17 Diurnal water temperature plots of Channel 5 48
4.18 Vertical temperature stratification in Channel 1, pool 6.
Measured August 11, 1976 54
4.19 Per cent days with stratification AT > 1.0°C and AT > 2.0°C
in Channel 1, pool 6 and pool 12 55
4.20 Per cent hours with stratification AT > 1.0°C—and AT >-l.O^C ~
in Channel 1, pool 6 and pool 12 - 56
5.1 Schematic streamlines at transition from riffle to pool in
Monticello experimental channels without and with transverse
winds 61
5.2 Schematic outflow concentration curve for ;jet mixing in experi-
mental pool 61
5.3 Temperature fronts in MERS channel 1 on November 17, 1976 63
5.4 Computed temperature fronts (normalized) which occur at location
17 in Channel 1 with and without heat transfer 65
6.1 Comparison of Eq. 6-46 and 6-47 with computed values of the
wind function, F 84
w
6.2 Comparison of Eq. 6-48 with computed values of the wind function
minus natural convection term 85
6.3 Comparison of Eq. 6-49 with computed values of the wind function,
F , minus natural convection term 36
w
6.4 Comparison of Eq. 6-50 with computed values of the wind function,
F , minus natural convection term 87
w
7.1 Control volume in time-distance coordinates for formulation of
finite differences 53
> ,__ .Gejieral__flo_w__chart _for_ MNSTEH,^ _tHe_ Minnesota Stre.am Water
.".Temperature "Prediction Model ......................... .....
7.3 Run A. Hourly water temperatures at three stations in channel 1
from January 31 through February 4, 1976 ........................
7.4 Run B. Three" hour water temperatures at three stations at
-------
Page
Run C. Three hour water temperatures at Station 3 (upstream)
and Station 17 (downstream) in Channel 8 106
7.5b Run C. Three hour water temperatures at Station 3 (upstream)
and Station 17 (downstream) in Channel 8 107
7.5c Run C. Three hour water temperatures at Station 3 (upstream)
and Station 17 (downstream) in Channel 8 108
7.6a Run D. Three hour water temperatures at Station 3 (upstream)
and Station 17 (downstream) in Channel 1 110
7.6b Run D. Three hour water temperatures at Station 3 (upstream)
and Station 17 (downstream) in Channel 1 Ill
B-l Steady state water temperature profile and "best fit" curve
for November 19, 1976 130
G-l Sample strip charts for wind velocity, wind direction, and
solar radiation 158
G-2 Sample hygrothermograph strip chart for degrees Celcius air
temperature and percent relative humidity 159
G-3 Weather station with hygrothermograph and shelter, pyrano-
meter, wind vane, and anemometer 160
G-4 One-year cycle of air temperatures. February 1976
through January 1977 161
-------
LIST OF TABLES
Number Page
4.1 Water Surface Slopes, Cross Sectional Areas, Wetted Perimeters,
Manning's Roughness Coefficients, Shear Velocities, and Shear
Stresses on the Bottom of Channel 1 on April 18, 197-7 ........... ^5
4.2 Water Surface Slopes, Cross Sectional Areas, Wetted Perimeters,
Manning's Roughness Coefficients, Shear Velocities, and Shear
Stresses on the Bottom of Channel 1 on June 15, 1977 ............ 16
4.3 Summary of Water Temperature Information Stored on Magnetic Tape. 35
4.4 Monthly Mean and Standard Deviation of Daily Minimum and Daily
Maximum Longitudinal (Spatial) Temperature Difference, and of
Diurnal (Temporal) Temperature Difference at the Outflow for
Channel 1, December 1975 - November 1976 ........................ 53
4.5 Frequency Analysis of Temperature Stratification in Channel 1
(January 1 - Decemoer 31, 1976) ................................. 53
5.1 Longitudinal Dispersion in MERS-Channels from Temperature
Front Data [[[ 68
5.2 Average Pool and Riffle.D from Temperature Front Routing Tests. 68
L
5.1 Wind Function Formulas Determined by Various Investigators ...... 88
7.1 Standard Error of Water Temperature Predictions when Weather
Parameters are Averaged Over 1, 3, 6, and 12 Hr Periods ........ 112
7.2 Standard Error of Water Temperature Predictions with Instan-
taneous Readings of Weather Data Taken at 3 , 6, and 1 Hr
Increments ....... ...... ........................................ 114
A-l Channel 1 Cross-Sectional Areas and Surface Widths for
Temperature Front Experiments .................................. 120
A-2 Residence Times to Each Station in Temperature Front
Experiments [[[ 120
A-3 Measured 0 Values for 11/17/76, 11/18/76, 11/22/76, 12/6/76,
-------
Number Page
B-l Data and Computational Results for Determination of Bulk
Surface Heat Transfer Coefficient and Wind Function
Parameters from Steady State Longitudinal Water Temperature
Profiles 127
B-2 Maximum and Minimum Values of K , F , and F -5.52{A6 )
S W W v
in sample computation due to maximum possible errors in data 132
xiv
-------
ABBREVIATIONS AND;SYMBOLS
A = cross sectional area (L )
A = mean cross sectional area over entire reach
B = surface width (L)
B = mean channel width over entire reach
C = concentration (ML )
C = fraction cloud cover (-)
c
C = cloudiness ratio
r _X o -1
c = specific heat of water (EM C )
P
a,b,c = coefficients
D = thermal diffusivity (L T )
D = 2D /U
L -T
D* = D BQ ~
L _L _]_
D* = D_ A B Q A
2-1
D = longitudinal dispersion coefficient (L T )
Ij
e = actual air pressure (mb)
a
e = saturation vapor pressure at dewpoint (mb)
a2
e = saturated vapor pressure (mb)
sw
f = bottom friction factor
h = mean channel depth (L)
n = mean channel depth over entire reach
-2 -1
H = convective heat loss(EL T )
_9 -i
H, = evaporative heat_loss ^EL ~__T! )
H = net long wave radiation (EL T )
-2 -1
H = net short wave solar radiation (EL T )
H = incoming solar radiation (EL~ T~ )
H = reflected solar radiation (EL~ T"1)
S L
i = distance node
3 = time node
k = Darcy permeability coefficient (LT~ )
EC = 2m/U = 2 K p"1 c ~1 h'1 U~l
XV
-------
K = bulk surface heat transfer coefficient
S - -1
L = latent heat of vaporization (EM )
ra = K p"1 c ~1 h'1
s * p
M = mass of tracer injection (M)
n = Mannings roughness coefficient
NS = 3T/3t
P = wetted perimeter (L)
p = atmospheric pressure (mb)
Q = flow rate (L3 T*1)
r = reflectivity or albedo of water surface
R = hydraulic radius (L)
RH = relative humidity (-)
S = rate of surface heat input (EL T )
S = slope of energy grade line (-)
SS = Q B"1 h"1 3T/3x
t = time (T)
T = water temperature ( C)
o o
T = air temperature ( C or K)
a
T = dew point temperature (°C)
T = equilibrium temperature (°C)
E
TQ = T(x=0)
T = surface water temcerature ( C or K)
3 " _!
U = mean channel flow velocity (LT )
U = mean channel flow velocity over entire reach
u^ = shear velocity (LT~ )
(Wftu) = wind speed function
z
W = wind velocity at elevation z above ground (L
z
x = distance along channel (L)
a = solar angle (degrees)
A8 = virtual temperature difference (°C)
- porosity of bed material (-)
o -4 -2
e = emissivity of atmosphere (E K L )
-2 o -4
£ = emissivity of water (EL K )
xvi
-------
Q = dimensionless water temperature = [T(t)-T(0)/ T(°o)-T(0)]
9 = dimensionless temperature = [T(x)-T ]/(T -T )
_ E O E
p = water density (ML )
a = Stefan-Boltzman constant
T = bed shear stress (FL )
o
xvi i
-------
ACKNOWLEDGEMENTS
The very understanding cooperation and planning by the Project Officer,
Dr. Kenneth E. Hokanson is gratefully acknowledged. The planning and per-
formance of many of the.field, experiments, required some^involvement and. help
from the permanent staff of the field station (MERS), expecially Mr. Thomas P.
Henry and Mr. Jack Arthur, Acting Chief of MERS. The cooperation of
Mr. Charles F. Kleiner in providing and maintaining some of the required
instrumentation deserved particular recognition and gratitude.
Funding for this study was provided through grants from the U. S.
Environmental Protection Agency, Environmental Research Laboratory - Duluth,
Minnesota 55362.
Reviews of a draft copy of the report by Professor John A. Hoopes,
University of Wisconsin, Madison, Wisconsin, and Professor Frank H. Verhoff,
West Virginia University, Morgantown, West Virginia, are gratefully acknow-
ledged.
XV 111
-------
SECTION 1
INTRODUCTION
The Monticello Ecological Research Station (MERS) is a Branch of the U.S.
Environmental Protection Agency's Environmental Research Laboratory in Duluth,
Minnesota. The MERS is located approximately 40 miles northwest of Minneapolis,
near Monticello, Minnesota, on 34 acres adjacent to a Nuclear Power Generating
Plant owned and operated by Northern States Power Company. An aerial photo-
graph of the Monticello Field Station is given in Fig. 1.1. The MERS has eight
soil bottom experimental open channels of approximately 520 m (1700 ft) length
each. The Mississippi River serves as the main source of water for the channels
which operate in a once-through mode. Heat exchangers using waste heat from
the power plant can be used to artificially heat the water in any of the chan-
nels. The channels are in an elevated area above the Mississippi River and
not normally subjected to flooding.
Flow rates in each of the channels as well as water stages can be con-
trolled. The discharge from the channels is returned to the Mississippi River
without any provision for treatment.
Heated cooling water effluents from power facilities create artificially
high temperatures in natural waters. These thermal additions also intensify
the dynamic character of natural water temperatures. Analysis of the higher,
more variable water temperatures is needed to aid in determining effects upon
the ecology of natural water's. Such studies are "conducted at the MERS.
Artificially heated water in "a river, lake, pond, estuary, or coastal"
area is cooled primarily through heat transfer with the air. Such heat trans-
fer has been thoroughly studied for large water surface areas but the relation-
ships developed may not apply well to shallow and narrow open channels. The
quantitative effect of wind velocity, which enhances heat transfer, is different
for narrow channels and wide surface areas.
-------
Fig. 1.1. Aerial photograph of the Monticello Station
experimental channels- (center); Right - Northern
States Power Company Monticello Power Plant.
Bottom - Mississippi River.
-------
Water temperatures in small open channels are highly dynamic due to the
time dependent character of the controlling processes: advective transport,
wind and flow induced mixing, and heat exchange with the air and the channel
bed. In the MERS field channels the unsteady character of water temperature
is especially acute since neither the mixing characteristics nor the effect
of advective transport are negligible in the channels. The channels are ex-
posed to a wide range of air temperatures, humidities, wind velocities, and
solar radiation. The corresponding heat transfer with the environment is
-signif icarvtt— a differ.ence_in_ temperature between the upper and lower end as
high as 15°C has been recorded and diurnal variations above 5 C have been
observed at the downstream end. The experimental channels provide an excellent
medium for the study of unsteady water temperatures in open channel flow and
their relationship to meteorological parameters.
The energy transport equation applied to an open channel of constant
cross section takes the following form:
h
p
D is a dispersion coefficient in the direction of flow (x-direction) , S is
L
a source or sink term which includes heat transfer with the surrounding
environment, t is time, U is mean channel velocity, h is mean depth, p is the
density of water, and c is the heat capacity of water. The potential im-
portance of the dispersion mechanism and of heat exchange with the air is
evident in Eq. 1-1.
Continuous water temperature records and meteorological data were
gathered^ and usedMio^jiietermine longitudinal dispersion and adapt known air-
wa'tef " heat transfer rela-ti-cnships to narrow open channels with well defined
non-uniform cross sections (crenelataons) . The effect of wind upon longitudinal
dispersion in the cnannels was studied. This information was needed to develop
a dynamic model of channel water temperatures. It was intended that the model
should be able to predict the diurnal water temperature variations -on a time-
scale of 1 to 3 hours, as well as the longitudinal gradients. The longitu-
dinal gradients were induced by adding artificially heated water to the
channels and oy the diurnal heating and cooling cycle at the channel surface.
-------
A time scale of 3 hours was found adequate to simulate the observed diurnal
variations.
The model simulation had to pay attention to several factors not usually
considered in stream temperature modeling: (a) the particularly complex
geometry of the channels caused by the alternation of riffles and pools,(b)
the presence of weeds and intermittent stratification and (c) the unusually
short time scale required. This led to a water temperature model which
is more dynamic than any other model known to the authors.
This study of the physical environment in the MERS field channels was
intended to provide some background for several biological studies at the
Research Station, including one on the effects of artificially high water
temperatures on fish and invertebrate populations. An analysis of the
temperature measurement problem, channel temperature regime, microhabitat
conditions, and temporal and spatial dynamics was essential for such investi-
gations. The planning of future channel experiments may be facilitated by
the results of this study. In a broader sense, this study also contains
material applicable to the study of temperature dynamics in thermally polluted
streams, to the design of irrigation canals for cooling water disposal, and to
the natural water temperature regime of small streams.
-------
SECTION 2
CONCLUSIONS
1. The field channels at the Monticello Ecological Research Station represent
a meso-scale ecosystem whose temperature regime and hydraulic charac-
teristics have been measured. The water temperature regime has also been
simulated in a dynamic numerical model using channel morphological and
hydraulic parameters as an input and air/water heat exchange as a water
temperature forcing function.
2. Development and operation of instrumentation for continuous monitoring
of flow rates and water temperatures in the channels from 1975 to 1977
was in the hands of MERS staff, in particular Charles F. Kleiner. The
investigators participated in the selection of the temperature sampling
sites shown in Fig. 4.13. The investigators also collected intermittently
data on flow velocity profiles, water temperature stratification, surface
water surface slope and flow cross sections. These data have been used to
describe the temperature regime of the channels, derive a number of
relevant physical parameters, and verify techniques for the analysis
of temperature dynamics as described below.
3. The nydraulic characteristics include:
Mannings roughness: n = 0.06 in the clear channels
n = 0.3 with weed growth
Mean flow velocities: U = 0.1 m/s in riffles
U = 0.02 m/s in pools
depending on flow rate and downstream channel control. Velocity profiles
are strongly affected by growth of filamentous algae in the rock sections
and macrophytes in the pools.
Darcy permeability coefficient of rock k = 0.18 m/s
Rock ,porosity e = 0.4-2.
-------
4. The soil profile below pool bends is composed of four layers: a highly
organic surface layer, water saturated sandy loam, water saturated
bentonite clay layer, and a dry sandy loam, in order, from surface
downward. The depth of the highly organic surface layer is variable and
increases with time. The thermal diffusivity of each layer was determined
by fitting soil temperature prediction to temperature measurements.
D = 0.0017 cm /s in the organic surface layer
D = 0.013 cm /sTn'fche wet sandy. loam_.helow.-
D.3_-5-.0013^ cm./s in the dry sandy loam below J:he_benton,ite- seal
The thermal diffusivity of the riffle was estimated to be D = 0.0031 cm /s
in the rock/water system without natural convection.
5, Vertical stratification of the pool sections was studied in the artificially
heated channel 1 and was found to occur most frequently in August. The
temperature gradient was found to occur primarily in the lowermost 20 cm
of the pool below a 60 cm isothermal layer; the temperature difference
was up to 5°C strong. Analysis of some 1975/76 data in pools 6 and 12
showed that a temperature differential of more than 1 C developed during
33 per cent of all hours in August and less frequently during other
months. The more downstream pool showed a stronger tendency toward strati-
fication, presumably because the cumulative effect of atmospheric heating
of the water in it is stronger.
5. A dynamic soil temperature model for the pool sections was developed,
primarily to assist in invertebrate studies. The model solves the
unsteady heat conduction equation and predicts the soil temperature
profile for a specified time variation of the water temperature at the
_ _gpil/wateJL interface^.. It was_verified (Fig. 4.11} that spil_temp_er_atures _.
__,„ jnay..be accurately, predicted with this "model.
'7. Three-hour "water temperature measurements in four of "eight fieia channels
have been recorded, processed and statistically analyzed for a period
from December 1975 through September 1977. A computer program and
User's Manual were provided to EPA.
8, Longitudinal dispersion in the MERS field channels was studied separately
using temperature fronts. The longitudinal dispersion coefficient was
found to be best approximated by
-------
DL = 7.47 * f ,
where Q is flow rate and B is mean surface width. The
effects of pool/riffle interaction appeared to be minor and longitudinal
dispersion was similar to that in a channel of uniform cross section.
In the range of very low mean flow velocities investigated (0.04 m /s
discharge and U = 0.02 m/s yi pool to 0.1 m/s in riffle), the typical
wind and vegetation effects could not be distinguished separately. A
typical value of D was 0.1 m~/s, with more detail provided in Tables
LJ
5.1 and 5.2.
Evaporative and convective surface heat exchange were found to be
dependent on wind speeds and natural convection potential. A wind speed
function
Wftn = 0.0096 (A6 ) + 0.0083 W2
was found to fit the experimental field data well. A9 represents a
virtual temperature differential defined by Eq. 6-30 and W is the wind
speed at 2 meters above the ground in m/s. The windspeed function is
used to calculate evaporative surface heat transfer in accordance with
Eq. 6-22 and convective surface heat transfer according to Eq. 6-26.
10. The water temperature regime was found to be highly dynamic for several
reasons:
(a) the shallow mean depth of the channels (< 1 m)
(b) the residence time of several hours in the channels
(c) the open cycle mode of operation whereby water from the Mississippi
River is fed (heated or unheated through the MERS channels in a
once-through mode)
Since the nydraulic conditions on thermal responses in the Mississippi
River and in the MERS field channels are very different, channel water
temperatures are usually in transient conditions with ooth strong
longitudinal gradients (several degrees C) and significant diurnal
variations (several degrees C) at the downstream end. Taole 4.4 provides
specific figures. Upstream conditions were more stable because the
inlet water came from the Mississipi River, a deeper water oody than
the channels.
-------
11. A finite difference, implicit computer model MNSTREM for the simulation
of the very dynamic channel water temperature regime was developed. It
was shown capable of simulating 1 hr, 3 hr or 6 hr water temperatures
over periods up to one month with standard errors of 0.2 to 0.3°C
between measurements and predictions. Weather data, specifically solar
radiation, air temperature, dewpoint and wind velocity are required as
input at least every six hours. Three hour input data are standard for
the model. If intervals of more than six hours are used, the accuracy
of the prediction suffers (Table 7.1). MNSTREM has unusually high time
resolution meaning that it can predict rapid water temperature changes
in very shallow water.
-------
SECTION 3
RECOMMENDATIONS
Water temperature variations "in time ancT~space"are a characteristic of
"any aquatic environment-and especially the MERS field channels.. Standard
methods for the characterization of water temperature dynamics and incorpora-
tion of these characterizations into water temperature criteria for pollution
control are needed. The MERS field channels lend themselves to studies of
dynamics very well. They should be used for continued studies of fate
and effect of various pollutants in an ecosystem at the meso-scale.
The focus of this study was on water temperature as a water quality
parameter. For the model formulation, a dispersion parameter and surface
exchange relationship adequate for water temperature modeling were derived.
Studies of pollutant (mass) transportiin the MERS channels require more
refined and detailed information on the flow and dispersion in a bed of
macrophytes since uptake and precipitation are added factors (not present
in the temperature model).
If further ecosystem studies in the MERS channels are conducted,
detailed investigations of the hydrodynamics and transport dynamics of a
pool and of a riffle, including the effects of vegetation must be conducted.
For example, tne effect of vegetation on longitudinal dispersion needs
-investigation.- -Slow. af_f ec.ts_must fee ._ss.r_t.ed._pu.t.Jr.om. the. f ield_ observations
to derive-the crue biological and chemical kinetics of the s_ystem.
The 3-hour water temperature~data base derived from continuous recordings
snould be analyzed further by time series (spectral, correlation) analysis.
In addition, weather data and flow data snould be digitized.
Continuously recording flow metering and channel stage recording
devices at each end of a channel should be installed in the MERS channels.
Flow rates were recorded manually once a day and stages only intermittently.
-------
SECTION 4
HYDRAULIC AND THERMAL CHARACTERISTICS OF THE MERS FIELD CHANNELS
4.1. CHANNEL GEOMETRY
Each of the eight experimental channels at the MERS is approximately
520 m (1700 ft) long, with alternating "pool" and "riffle" sections of
approximately 30.5 m(100 ft) length each and one 90 degree bend from an
east-west to north-south alignment. There are nine pools and eight riffles
in each channel. The approximate geometry of the pool sections is 3.5 m
(12 ft) surface width and 0.8 m (3 ft) maximum depth. The pool cross section
is approximately parabolic. The riffle sections are contracted areas to
increase flow velocity. They are trapezoidal in shape and formed by 38 mm
(1.5 inch) gravel placed in the channel. The riffle width at the water
surface is approximately 2.5 m (8 ft). Riffles are very shallow, e.g. 0.3 m
(12 inches). As the water stage can be controlled by a bulkhead at the end
of each channel, the exact water depths and surface widths can also be changed.
The channels are arranged in pairs, separated by a high berm. There is
a gravel road between channel pairs. Figure 4.1 shows one pair of channels.
Locations in the channels are numbered, with the number 1 assigned to
the inlet, numoer 2 to the first pool, number 3 to the first riffle, etc.r
up to the outlet which is numbered 19. In addition, each channel was given a
number from 1 to 3. Figure 4.2 gives measured average riffle and pool cross
.sections in_Channel 1 .(northernmost channel). Cross sections, are not uniform
in all cases. One example of a non-uniform- riffle cross section is given in
Fig. 4.3.
10
-------
Fig. 4.1. MERS channels showing the transition from pool
(foreground) to riffle to pool. The channel on
the left is flowing full, while the channel..on
the right is partially drained, exposing the gravel
sides and bottom of the riffle section.
11
-------
F
305
I*-30 5
185.9
208
"T
193
|< 62.*
62
.5 —*|
Average riffle cross section for Channel 1 on April 6, 1977.
457 —»j
457
216 4
Average riffle cross section for Channel 1 on June 15, 1977«
341 4
Average pool cross section for Channel 1 calculated from
measurments made from July 13, 1977 through July 19, 1977.
Fig. 4.2. Channel cross sections, (dimensions in centimeters)
12
-------
I-1
LJ
Upstreom
Section
Middle
Section
Downstreom
Section
'00m
Scolei '
05m
Fig. 4.3. Cross sections of riffle 5, channel 1 on April A, 1977.
-------
4.2. CHANNEL HYDRAULICS
Flow Rates < and Residence Times
Flow rates in each channel are controlled and metered in the control room
of the office building. V-notch weirs are installed at the head (inlet) of
each channel and are available for installation at the end of each channel.
Maximum design flow for a channel is on the order of 37.9 i/s (600 gpm) .
Actual flow rates have been as high as 39.1 £/s (620 gpm). In its first year
of operation the flow rate in Channel 1 was approximately 18.9 H/s (300 gpm),
and" the total hydraulic residence time of the water in the channel on the
order of 12 hours. In the second year, beginning in August 1976, the flow rate
was approximately 28.4 2,/s (450 gpm) to 37.9 H/s (600 gpm) , with residence
time on the order of 4 hours. A lower stage was used with the higher flow.
Roughness and Water Surface Slope
The hydraulic roughness was determined in Channel 1 on three different
occasions. Water surface slopes were measured in Channel 1, enabling calcula-
tion of roughness coefficients, wall shear stresses on the wetted perimeter,
and shear velocities. Two sets of values are presented in Tables 4.1 and 4.2
where the station numbers refer to pool and riffle locations. (Counting in a
downstream direction, the odd-numbered stations are riffles and the even-
numbered stations are pools.) Computations were made for three different flow
rates and for clear and weed choked conditions in the channel.
On April 18, 1977 the channel was relatively free from weeds and obstruc-
tions. On June 15, 1977 heavy macrophyte growth in both the riffles and pools
had raised the water surface level substantially from the April level.
Manning's roughness coefficient 'n' is defined by the Manning equation
(see e.g. Olson, 1973 or-Barnes, 19"67 or Limerinos, 1970).
,4-1,
where Q = flow rate in m /s,
A = cross sectional area of the channel in m ,
R = hydraulic radius (m) of the channel cross section, defined as
the cross section area A divided by the wetted perimeter P, and
SE = slope of'-energy grade line
14 - -
-------
TABLE 4.1. WATER SURFACE SLOPES (S), CROSS SECTIONAL AREAS (A), WETTED
PERIMETERS (P), MANNING'S ROUGHNESS COEFFICIENTS (n), SHEAR
VELOCITIES (Uj, and SHEAR STRESSES ON THE BOTTOM OF
CHANNEL 1, T , ON APRIL 18, 1977. (Flow rate = 0.0320 m /s)*
Section
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
S
0.000537
0
0.000263
0.000104
0.000208
0.000116
0.000436
0.000053
0.000312
0
0.000573
0
0.000260
0.000104
0.000210
A P n
(m ) (m)
.316 2.13 0.064
.372 2.04 0.061
.372 2.16 0.052
.307 2.04 0.057
.297 1.92 0.047
.279 1.S9 0.058
.297 1.92 0.043
.214 1.71 0.024
U T
° 2
(m/s) (N/iri )
0.028 0.78
0.022 0.48
.
0.019 0.36
0.025 0.63
0.022 0.48
0.029 0.84
0.020 0.40
0.016 0.26
* After Hahn et al, 1978b.
15
-------
The transitions between pools and riffles appeared to produce an insigni-
ficant loss compared to the roughness of the riffle sections.
It is noteworthy that in the clear channels (without growth) the hydraulic
headless occurs predominantly in the riffle sections. With heavy weed growth
the contributions of pools and riffles may, however, be equal.
The average values of Manning's roughness coefficient were n = 0.056 on
April 18, 1977, and n = 0.29 on June 15, 1977. The associated measured head-
losses along the channel were 0.097 m (3.8 in.) on April 18, 1977, and 0.21 m
(8.3 in.) on June 15, 1977. The measurements excluded the first and the last
pool.
It should be noted that the roughness coefficients calculated for the
channel with heavy weed growth account not only for the additional resistance
to flow caused by the weeds, but also implicitly include the effects of the
weeds in altering the channel cross section and length of flow. The overall
pool and riffle cross sections with no area reduction were used in computing
the 'n' values. In parts of the weedy sections the flow may have been appor-
tioned between several subsections, thus producing a divided flow with several
cross sections and hydraulic radii.
Flow Velocities
Mean flow velocities in riffles and pools vary, of course, with flow rate
and stage. Mean values applicable to conditions observed on April 18, 1977
with a flow of .033 m /s (508 gpm) were:
Riffles: .108 m/s (.335 ft/s)
Pools: .019 m/s (.062 ft/s)
Velocity distributions as a function of depth were measured on several occasions
with the aid of a micro currentmeter. Sample results are shown in. Figs. 4.4
through 4.6. The effect of the attached filamentous algae growth on the
velocity distribution is very apparent.
17
-------
178cm
B
algae
growjh
0 10 20 30
0 -
E 10
8 '5
20
25L
algae
growth
Velocity (cm/sec)
0 10 20 30
0 -
Location A
10
15
20
25
10 20 30
0 -
Location B
10
15
20
25
growth
Location C
.Fig. 4.4. Velocity profiles in Channel 1, riffle 5, looking downstream.
Measured October 30, 1976.
-------
10 20
- 5
e
f '0
-------
320cm
10
o
Velocities
too small
to measure
algae-
growl h
7"
Velocity (cm/sec)
0 25 50
0 -
20
£
u
a.
-------
4.3. SOIL AND ROCK CHARACTERISTICS
Grain Size Distribution
Size distributions of the pool and riffle bed materials and of the sedi-
ment deposited within the first riffle section of Channel 1 were measured.
The rock size distribution for the riffles was made on a random 325 N (73 Ib)
sample of the gravel collected from the sides of all riffles in Channel 1.
Two in. (50.8 mm), 1.5 in. (38.1 mm), 1.0 in. (25.4 mm), and 0.75 in. (19.0 mm)
sieves were used. The results are plotted in Fig. 4.7.
A visual inspection of the rock material in the first (farthest upstream)
riffle of Channel 1 indicated that it contained a considerable amount of fine
sediment in the form of a mud or slurry. No such sediment was observed in the
last (farthest downstream) riffle of Channel 1. Permeability measurements
made in the two riffles indicated a lower permeability of the first riffle
compared to the last, as would be expected if there were sediment deposition
upstream. The origin of the sediment accumulation is presumably the Mississippi
River water. The first pool of each channel acts as a settling basin for
suspended sediment carried by the river. The sediment removal in the first
pool is, however, not complete and some sediment is apparently deposited in
the rock of the first riffle.
A sample of the sediment in the first riffle was obtained by digging a
hole into the gravel cross section and collecting samples of the gravel and
the slurry from througnout the vertical cross section. The sediment adhering
to the rocks was later collected by washing the rocks.
The mechanical analysis of the sediment was performed according to the
standard method (AASHO Designation: T88-57) described in the Soils Manual
for Design of Asphalt Pavement Structures published oy the Asphalt Institute.
The resulting sediment size distribution is shown in Fig. 4.7.
The size distribution data for the pool bottoms in Channel 1 obtained in
1975 and available at MERS were also plotted in Fig. 4.7. The bed material
was previously sampled in each pool and a size distribution was determined
for each individual sample. The per cents finer than in each category were
averaged for the nine pools and the resulting size distribution curve was
plotted. Also, the per cent finer than in each category was determined for
21
-------
U S Stondord Sieve Numbers
25 18 10
O Sediment in grovel of first riffle
O Upper bound for soil in pool bottoms
A Average pool soil size distribution
O Lower bound for soil in pool 'bottoms
O Gravel from riffles
I i I I 11 ll
I I I I
...I
02 05 .1 2 512
Sediment Size in Millimeters
10 20
50
Fig. 4.7. Size distributions for material in Channel 1 in MBRS.
-------
each pool, and the highest and lowest per cents among all the pools in each
size category were plotted as upper and lower bounds for the sediment sizes.
Permeability of Riffle Rock Sections
The hydraulic (Darcy) permeability of the rock was determined by separate
experiments for the downstream part of the last riffle in Channel 1 and for
both the upstream and the downstream part of the first riffle in Channel 1.
In the case of the last riffle, the water level in the pool immediately—
downstream was lowered so that all flow was only- through the-rocks- and-piezo--
metric elevations were measured to determine the water table profile along the
riffle section as a function of location and of time. Experiments were con-
tinued for several hours. For the downstream riffle a quasi-steady state
condition was reached. The Darcy permeability coefficient of the rock was
found to be 18 cm/s. With a channel flow on the order of 0.0315 m /s (500 gpm),
the mean flow velocity through the pores of the rock was calculated to be
approximately 0.02 cm/sec with a rock porosity on the order of 0.42.
Porosity of the rock was found by measuring void space in a container
of 19.6 liters total volume. Average porosity from several samples was 0.44.
Soil Thermal Diffusivity and Soil Temperature in Pool Sections
A study of soil thermal diffusivity and temperatures below the bed of the
pools was made for two reasons. One reason was to determine if the heat
exchange between the pool water and the bed had an important influence on the
water temperature regime in the channels. The influence was found to be
negligible. The second reason was that aquatic biologists studying benthic
organisms desired to know something about the temperature regime in the mud.
The soil temperature regime was studied oy heat transfer analysis (unsteady
and one-dimensional)., - - ' .- " •
- The soil - temperature profile as a function of depth and time was predicted-
for known water temperature in the pool above the soil/water interface. A
numerical finite difference model of the unsteady, one-dimensional heat con-
duction equation was developed (Gulliver and Stefan, 1980). Measurements of
soil temperature profiles were obtained from thermistor arrays placed in the
bed of pools 6 and 12 in Channel 1.
23
-------
Analysis of soil composition and of measured soil temperature profiles
led to the conclusion that a four-layered system had to be used: the first
layer of the pool bed was a highly organic layer of variable depth (0-15 cm).
It appeared that this layer had developed since the channels went into
operation and that the organic material was the residue of several growth
seasons. The second layer (below) was a water saturated sandy loam of 38 cm
(15 in.) thickness. This second layer was backfill material which rested
on a 5 cm (2 in.) bentonite layer which had been installed to prevent the
channels from leaking. The fourth layer below was a dry sandy loam. Figures
4.8 and 4.9 give cross sections.
Several heat wave experiments were performed. The soil surface was
suddenly exposed to heated water, after a cold period, and the penetration of
the warm front into the soil with time was recorded by the thermistor arrays.
An example of a record is given in Pig. 4.10. Analysis of several such fronts
made it possible to refine initial estimates of soil thermal diffusivities.
Although there were three thermal diffusivities and one layer depth to
be determined, the heat wave experiment encompassed two unsteady and two
quasi-steady periods and the problems of fitting four parameters could be
solved. The surface layer thermal diffusivity and surface layer depth
were highly interacting, but the use of unsteady periods helped resolve the
problem to a reasonable degree of accuracy. The following soil properties
were determined by fitting predictions of the numerical finite difference
model to soil temperature measurements:
Highly organic (surface) layer depth = 12.5 cm (5 in.)
Thermal diffusivities:
Highly organic surface layer D = 1.7 x 10 cm /sec
Water saturated sandy loam" • D "= 1.3 x 10~ cm /sec"
and bentonite layer
Dry sandy loam D = 1.3 x 10 cm /sec
(lowest layer)
The above thermal diffusivities were used to predict soil temperatures
in pool 12 during other heat wave tests. A set of predicted soil temperatures
is compared with thermistor soil temperature measurements in Fig. 4.11. A
24
-------
Woter
^Bottom of Pool
\\\\\\\\\\\\\\\\\\
Highly organic surface layer
of variable depth " 2.5 cm
\
/
water
sandy
f
i
saturated
loam
38 cm
O- /ft
0
n —
*
-5'm I0*cm f '
T >y 1
23cm
JT
i
36 cm
5cm _be_ntomte_ _ _ _ ^V^
dry sandy loam to end of domain
Thermisfor Probes
Fig. 4.8. Soil profile below pool and location of thermistor probes.
Channel 2
Channel 1
Dry sandy loam
Sentonife layer
Fig. 4.9. Schematic cross section of one channel pair
drawn to approximate scale.
25
-------
Temperature (°C)
Fig. 4.10. Measured soil temperature profiles in
... , . p^ol_fi_duxing_nea±_wave_test:. .warming.,- ...
_„_ - -„, . period. - - ,
26
-------
Measured Soil Temperature
O 6 cm depth
24 j 25 | 26
December 1977
I 2 I 3
Jonuory 1978
Fig. 4.11. Con,parison of measured and predicted (solid lines) soil temperatures in Pool 12 during
heat wave test. Notation for depth of predicted soil temperatures: (1) soil surface,
(2) 2.5 cm, (3) 6 cm, (4) 16 cm, (5) 28 cm, (6) 41 cm, (7) 60 cm, (8) 1 m.
-------
surface organic layer of 6.25 cm (2.5 in.) was assumed. The standard error of
the predicted soil temperatures in the example is 0.48 C.
More information on the soil temperature investigations conducted at the
MERS is given in Gulliver and Stefan, 1980.
Rock/Water Thermal Diffusivity and Temperatures in Riffle Sections
Water temperatures in the voids of the rocks forming the riffle section
were studied for the same reasons as the soil temperature profile in the pool.
It,may be recalled that the rocks are of 30 mm average size and form a system
of high hydraulic permeability.
The transfer of heat from the water flowing over the riffle into the
water and rock below is potentially by the following processes.
(a) Horizontal flow of water through the pores of the rock. The
hydraulic permeability was determined previously to be 18 cm/sec
—4
and the effective pore flow velocity on the order of 0.4 * 10
to 2.5 * 10 cm/sec, assuming a porosity of 0.44.
(b) Vertical conduction through the water and the rock mass.
Estimates of thermal diffusivity of the water and the rock
as found in the literature are 5.9 x 10~ cm /sec for dry
gravel and 1.4 x 10 cm /sec for water.
(c) Vertical forced convection through the voids of the rock
system. This could also be looked at as a downward hydro-
dynamic diffusion process.
(d) Vertical natural convection through the voids of the rock
system, when unstable water temperature stratification
-- occurs-
" " "The penetration of a heat wave from the bed surface into the rock of
riffle 17 was studied experimentally from December 19, 1977, through
January 2, 1978 (Stefan and Gulliver, 1980).
A vertical array of six thermocouples at 2.5, 7.6, 15, 25, 38, and 51 cm
below the rock surface was installed in the center of riffle section 17 in
Channel 1. The location is approximately 15 m from the upstream end of the
riffle. The lowermost thermocouple is in the sandy loam below the rock.
28
-------
Temperatures at the six sensor locations were read once daily during each work-
ing day following an artificially imposed surface water temperature rise. The
rock temperature profiles indicated by the sensors are shown in Fig. 4.12.
The temperature response in the upper 15 cm (6 in.) appears to be with a
time lag on the order of only a few hours. This rapid response as well as the
shape of the temperature distributions with depth suggests that a very effec-
tive vertical mixing mechanism in the upper 15 cm-of _the rock exis-ta. This. —
mechanism could be a hydrodynamic diffusion or forced convection process induced
by the bed shear stress of the free surface flow "in the riffle section. ~ -
Below 15 cm depth the temperature response is delayed but still quite
rapid. It is noteworthy that between December 20 and 22 and over a 48 hour
period the temperatures in the rock rose by the same amount as the surface
water temperature, namely 2 C per 24 hours. This observation would suggest
that the rate of downward heat flow during that time period was constant and
that the sandy loam (below the rock) was the main heat sink.
The observed rapid changes in rock temperature profiles cannot be explained
by horizontal flow through the rock. The horizontal advance of a temperature
front based on uniform horizontal flow and measured permeability would be
about 5 m/day in terms of real pore flow velocity.
Density stratification of the water in the rock iiass is associated with
vertical temperature gradients. Stratification may be responsible for
the steeper temperature gradient below 15 cm than above. It reflects a
smaller vertical thermal diffusivity than above 15 cm depth.
During the cooling phase the vertical temperature profiles are more
--un-i-feem than-during, the heating. This_reflects _the effect of natural vertical
---eoRvecfcipn ,induced,bit__d,e_n.5j-.ty_.instability. _ -
An analysis of"the data described" by Stefan and Gulliver -(1980), gave
-3 2
a thermal diffusivity of 3.1*10 cm /s for the rock/water system at depths
between 15 cm and 38 cm below the riffle bed. This value is applicable when
the water temperature in the rock section provides a stable density strati-
fication (channel water lighter than pore water in the rocks).
-------
Temperature, °C
6 8
Fig. 4.12. Vertical temperature profiles in
Riffle 17 during period of rock
temperature study.
30
-------
In the upper 15 cm of the riffle bed a mixed layer appears to exist.
Water temperature in that layer is very nearly the same as the channel water
temperature.
Under unstable conditions (pore water lighter than channel water) natural
convection takes place and thermal diffusivity is much increased.
4.4. OBSERVED OPERATIONAL WATER TEMPERATURE REGIME
Water Temperature Instrumentation
Temperature sensors (thermocouples) were installed at the upstream and
downstream ends and in the center pools and riffles of each of experimental
channels 1, 4, 5, and 8. The sensors provided continuous water temperature
information which was recorded on four strip charts for the four channels.
Three stations per channel were not sufficient, however, for the analysis of
stratification, surface heat transfer and mixing. Channel 1 was selected for
in-depth study of temperature regimes, and 21 continuous recording sensors were
installed in Channel 1 (schematic in Fig. 4.13). Ten of the sensors measured
water temperature 12.5 cm (5 in.) to 25 cm (10 in.) below the water surface.
There were two vertical arrays of four sensors each installed in Pools 6 and 12
which provided information on water temperature stratification and soil tempera-
ture. Three sensors were within 12.5 cm (5 in.) from the bottom of pools 4,
10, and 18, and one sensor 2 to 5 cm (1 to 2 in.) in the rocks of riffle 17.
All temperature sensors were connected to Leads-Northrup multi-channel color-
coded temperature chart recorders. Temperatures were recorded from November
1975 through December 1977.
Water Temperature Data Processing
Elements
The water temperature data processing described herein is for"water
temperatures measured at several stations in channels 1, 4, 5, and 8 of the
MERS. The data processing had the following elements.
•• Extraction from strip charts
• Review and estimation of missing or faulty data
• Calibration
• Storage
31
-------
Inflow from
Heot Exchonger
Channel Nos
87654321
Legend
i
V = vertical array of 3 probes placed at 5 and 10 in fronj channel
bottom and IS in. from water surface |
M = additional probe placed 1 and 2 in. in mud.
K = additional probe placed 1 to 2 in in riffle gravel. , ,
f
S = additional probes at 1,2, 3, 5, 10. and 15-20 in. in sojil(below
pool for non-continuous soil temperature profile
R = addtional probes at 1, 3, 6, 10, 15, and 20 in. in gravel below
riffle for non-continuous riffle gravel temperature profile. >
n = probe placed at 10 in from the gravel bottom rather than S
t
Placement of all other probes noted is 5 in. from bottom in both
riffles and pools
Outflow to
Mississippi
River
19
Fig. 4.13.' Thermal probe placement in MERS experimental channels (four channel pairs)
-------
These elements will be described in sequence.
Extraction of data
Water temperature data were extracted for four channels including two
hot channels, with up to 15 C above the Mississippi River temperature, one
warm channel, with up to 8 C above river temperature, and one ambient channel
which has no heat addition above river temperature.
For the statistical analysis, stations 2, 3, 9, 10, 17, and 18 of channels
4, 5, and 8, as well as stations 3, 4, 6, 12, 17, and 18 of Channel 1 were
selected. For each station, 3-hour water temperatures beginning at midnight
of each day were read off the strip charts and hand recorded on data sheets
(eight values per day).
Estimation of data
For periods during which no temperatures were recorded, estimation of
water temperature became necessary. This was done by one of the four follow-
ing procedures.
(a) Linear interpolation between water temperatures at upstream .
and downstream stations.
(b) When the period of missing data was only a few hours, and
when no water temperatures were recorded for all stations,
linear interpolation between the last water temperature
before the gap and the first water temperature after the
gap was used.
(c) When a temperature record was not available for the first
station (or the last station), linear extrapolation using
temperatures from the two immediate downstream (or upstream)
' stations'was used. • ' " "
(d) When no temperature record was available for all stations and
for a long duration, temperatures from another channel mul-
tiplied by a coefficient were substituted (linear intra-
channel transfer) .
Details on the periods of missing data, their causes, and the procedures used
in estimating the data were documented by Fu and Stefan (1979b) .
33
-------
Erroneous data were also found in the records. Some water temperatures
recorded^on the strip charts were not representative for several reasons:
(a) sensor might be in the air due to low flow, (b) sensor might be buried in
the sediment, and (c) sensor might not be recording the true water temperature
values. Owing to these reasons, estimation of corrected water temperature data
was also necessary. The procedure is explained in Fu and Stefan (1979b). Periods
during which estimates had to be made are also identified by Fu and Stefan (I979b)
Calibration
From January 1976 through June 1976, one complete calibration and numerous
spot checks of the recorded temperatures in Channel 1 were made. After July
1976, temperature calibrations were made approximately every month. Channels
4, 5, and 8 were included.
The corrections obtained from these calibrations were tabulated by Fu and
Stefan (1979b). Between the dates when calibrations were made, strip charts were
examined to identify discontinuities in the records. Calibrations were considered
valid only up to these points. Operational records and abbreviated remarks on
the strip charts were also consulted for this purpose.
Corrections of less than + 0.2 C were ignored. Corrections >|0.2| were
added to the hand recorded water temperatures. The entire water temperature
records for Channels 1, 4, 5, and 8 were stored on a magnetic tape. No cali-
brations were applied to the estimated water temperatures.
Water temperature data storage
All calibrated water temperatures were stored on a nine track magnetic
tape. The length of the temperature record of each station is listed in Table
4.3. The magnetic tape is divided into four files (multi-file tape). Each
file contains a temperature record of a single channel. Within each file,
temperature records are stored in the order listed in Table 4.3. The first
file contains temperature records of Channel 1, the second file those of
Channel 4, the third file those of Channel 5, and the fourth file those of
Channel 8. A computer program was used to read the raw data (without calibra-
tion) , to apply the calibration and to write the calibrated water temperatures
on the magnetic tape. More detailed information on the data tape can be
found in Fu and Stefan, (1979a).
34
-------
Ul
tn
1
1
TABLE 4.3. SUMMARY OF WATER TEMPERATURE INFORMATION
i
File , | Station
Sequence Set | Set No. in
No. Identifier Identifier Channel Order of
QN SI FI No. Storage
1 PFILED | CHANNEL 113
! 4
1 6
i 12
! 17
i 18
1
t
I
!
i
2 PFILED i CHANNEL 442
! 3
11
; ' 9
10
1 17
18
t 1 }
I •
3 PFILED CHANNEL 552
3
' 1 ' 11
9
10
17
18
i
i
' 1
i
STORED ON MAGNETIC TAPE (Cont'd)
i
Period
12/4/75
12/4/75
12/4/75
12/4/75
12/4/75
12/4/75
10/20/76
11/26/75
11/26/75
10/19/76
10/10/76
11/26/75
10/20/76
10/19/76
11/24/75
11/24/75
10/18/76
10/19/76
11/24/75
10/19/76
- 9/19/77
- 9/19/77
- 9/19/77
- 9/19/77
- 9/19/77
- 9/19/77
TOTAL
- 9/19/77
- 9/19/77
- 10/18/76
- 9/19/77
- 9/19/77
- 9/19/77
- 9/19/77
TOTAL
- 9/14/77
- 9/14/77
- 10/17/76
- 9/14/77
- 9/14/77
- 9/14/77
- 9/14/77
•
No.
of
Days
656
656
656
656
656
656
3,936
335
664
328
336
335
, 664
335
2,997
331
661
329
332
331
661
331
No. of
3 -hour
Temperatures
5,248
5,248
5,248
5,248
5,248
5,248
31,488
2,680
5,312
2,624
2,688
2,680
5,312
2,680
23,976
2,648
5,288
2,632
2,656
2,648
5,288
2,648
TOTAL
2,976
23,808
-------
TABLE 4.3 (Cont'd). SUMMARY OF WATER TEMPERATURE INFORMATION STORED ON MAGNETIC TAPE
File
Sequence
No.
QN
Set! Set
Identifier Identifier
SI ' FI
Station
No. in
Channel Order of
Mo. Storage
Period
No. No. of
of 3-hour
Days Temperatures
U)
-------
Water Temperature Data Analysis and Results
Elements
The water temperature data analysis as described herein included
• Development of a statistical analysis package
• Computation of a number of statistical parameters from 3-hour
data
• Determination of the seasonal cycle
• Determination of the diurnal cycle
• Analysis of longitudinal gradients
• Analysis of water temperature stratification in the pools
These elements will be described in sequence.
Development of a statistical analysis package
A computer program WTEMP1, described in detail by Fu and Stefan (1979a),
to extract from the 3-hour water temperature data statistical information
was developed including:
• Water temperatures at three probabilities of occurrence
specified by the user
• Maximum and minimum water temperatures in the period
analyzed
• Mean, standard deviation and dimensionless skewness
coefficient of the water temperature
WTEMP1 can handle a series containing .up to 19,200 individual water tempera-
tures, equivalent to 400 days of data from six stations in a single channel.
WTEMP1 can also nandle sliding statistic calculations for- sample periods of
up to 7 days. The-output is in tabular- form.
A second interactive program WTEMP2 for computing statistics on specified
data time intervals, composite station locations etc. from daily statistics was
also developed (Fu and Stefan, 1979a).
Computations
Computations have been carried out for individual stations and channel
composite using weekly sliding samples. The tabular output format from WTEMP1
is shown in Appendix H. The complete output comprises 151 pages for Channel 1
and 126 pages each for Channels 4, 5, and 8.
37
-------
Seasonal water temperature cycle
A graphical presentation of the statistical computer output is shown in
Figs. 4.14 and 4.15. The weekly mean, maximum and minimum water temperatures
have been plotted for two channels and several stations. Channel 1 was
artificially heated, and Channel 5 was not. The seasonal water temperature
cycle can be easily seen. The seasonal amplitude is on the order of + 15°C.
The annual mean differed by channel and station depending on the artificial
heating.
Water temperature fluctuations superimposed on the seasonal cycle are either
due to irregularities in the artificial heating of the inflowing water or due
to weather. Those due to heating are smoothed out along the channel.
Diurnal water temperature cycles
Diurnal water temperature cycles are shown in Figs. 4.16 and 4.17. Because
the inflowing water is taken from the Mississippi River, upstream stations in
the field channels reflect conditions found in the river, modified (or not) by
artificial heating. Because the field channels are shallower than the river,
diurnal water temperature amplitudes at the downstream end of the channels
are larger than upstream . Amplitudes of diurnal water temperature cycles
also vary strongly with season as seen in Fig. 4.16. Winter values have been
on the order of 1°C upstream and 2 to 3 C downstream. In the summer 2°C
diurnal amplitude upstream and 5 C amplitude downstream are not uncommon.
Longitudinal water temperature gradients
Water temperature differentials between the upstream and the downstream
ends of a channel are of interest for fish studies in the field channels. Such
longitudinal water temperature gradients are dependent on several other hydro-
thermal parameters as will be shown by theoretical analysiVin a later section.
It may; suffice^.here. to indicate that the- longitudinal_ AT can be expected to
rise with the addition of artificial heat. A greater flow rate will decrease
the longitudinal AT. Since the geometry of the MERS channels is fixed in
terms of length and width, the two most important controllable parameters are
flow rate and heat input. Actual values of longitudinal AT's will then depend
on the amount of heat rejected to the atmosphere, which is a weather depen-
dent process.
The daily minimum and maximum temperature difference between the inflow
and outflow of the artificially heated Channel 1 was recorded over a one-year
38
-------
uj a-
o
UJ
8-
YERR(S)
STflTION 4
s-
= 750
8-1
9
8-
Si
itn
n
in
LU-
a
8-
i TSO i.aa
YEHR(S)
STRTION S
i MO i ru
rEfirHS!
Fig. 4.14. Seasonal water temperature plots of Channel 1. (Cont'd)
39
-------
tfl
UJ
o
8-
CHflNNEL 1
STflTION 12
t sa
—i —r-
I DOC I.SO
8-1
in
in
UJ g.
C3
UJ
a
8-
1.73) 2 000
YEflR(S)
STflTION 17 •. Y'/W*
IM75 - a
1 * 1
°7SC 1. 000 1.250
8-,
*
' * ,"» '..• «"«••
1 1 [ ~ 1
1 SOD l.TSO 2 ODD 2 ZSO
YEflR(S)
1 i
2 SOD 2 TSO
!3
CO
UJ
a
s-
STflTION 18
YEflR(S)
Fig. 4.14. (Cont'd). Seasonal water temperature plots of Channel 1.
40
-------
CHRNNEL 5
3-
8
3
in
Bg-
UJ
Q
8-
YEflRtS)
STflTION 9
YERHt S)
Fig. 4.15. Seasonal water temperature plots of Channel 5.
(Cont'd)
41
-------
B-i
;
8-
R
3
IT)
e
a
8-
CHflNNEL 5
STflTION 10
—T~
i oao
—r-
1 23)
8-1
a
Si
^
tn
_j
LU g.
u 9
O
UJ
Q
8-
1 7S) J 000
YEflR(S)
STflTION 11
°TSQ I 000 1.291 I 500 l.TSO
8-
9
YEflR(S)
1
2 000
2.250
I
2.500
1
Z.7SO
I 750 I 000
YERR(S)
Fig. 4.15. Seasonal water temperature plots of Channel 5. (Cont'd)
42
-------
in
^>
to
C
UJ
Q
8-
CHflNNEL 5
STRTION 18
—r~
2 UO
1
2 73
—r~
1.2SO
I.TSD
YERR(S!
Fig. 4.15. (Cont'd). Seasonal water temperature plots of Channel 5.
43
-------
SIllllQll 3
STHIIUN 12
ClintltlEL I
1/13/76 fO l/2b/76
-Leu
• Ol • Ol
onus i
S1RIION 17
tOt too
OflTIS)
SlflllUII IU
utmsi
OHMSI
Fig. 4.16. Diurnal water temperature plots of Channel l.(Cont'd)
1/13/76 to 1/26/76
-------
OHYISI
Fig'. 4.16. Diurnal water temperature plots of Channel 1.
4/25/76 to 5/8/76
(Cont'd)
-------
.£>.
a\
CIIHtllllL I
II/U//6 10 11/22/76
unrisi
Fig. 4.16. Diurnal water temperature plots of Channel 1. (Cont'd)
11/9/76 to 11/22/76
-------
binilUII 3
ofmsi
SIMIOH S
SlflllUN
I •> 4 flu
SlflHON 17
CHHNNLL I
6/17/V/ TO b/JO/77
• oo ti to
oim si
DHYI5I
tta »»
onrisi
Fig. 4.16 (Cont'd). Diurnal water temperature plots of Channel 1.
6/17/77 to 6/30/77
-------
ClUUltttL S
1/19/77 10 2/1/77
tf a
DfiTIS)
smnmj is
oiirisi
<*ai *,o>
Fig. 4.17. Diurnal water temperature plots of Channel 5. (Cont'd)
1/19/77 to 2/1/77
-------
51 fill UN 2
40) ta
OHYISI
SiniiON 3
onris i
onrisi
sinnuii ID
STfHION 17
CIIIIMKEL S
1/19/77 TO S/2/7/
JOU <0u t OO tOO
onnsi
• OU !• U
DfiYI SI
SUII ION 10
onrisi
Fig. 4.17. Diurnal water temperature plots of Channel 5. (Cont'd)
4/19/77 to 5/2/77
-------
LH
O
•la,
8
bIHIIUII 10
siniioii
ClllVMIH 5
11/9/Vb 10 II/2.WO
DI1TISI
sL^/v
u OJ
onrisi
SIIIIIUII 18
UHllbl
I-1
Fig. 4.17. Diurnal water temperature plots of Channel 5. (Cont'd)
11/9/76 to 11/22/76 , '
-------
5IIII1UN 2
• 04 «0
bullion 3
bullion U
OHTISI
P. a-
H 8
y n
SUIT HIM 10
CIIIINIIlL 5
//I3/7V 10 7/2b/7/
oj i at
oiirisi
simian i/
omrisi
SIH1IUII 18
\_.
iiirrisi
Fig. 4.17 (Cont'd)
Diurnal water temperature plots of Channel 5,
7/13/77 to 7/26/77
-------
period. Complete results are given in Hahn et al, 1978a. Table 4.4 summarizes
the monthly mean and standard deviation of the longitudinal temperature differ-
ences.
Daily diurnal water temperature variations at the outflow of Channel 1
were also recorded over a one-year period (Hahn et al, 1978a). The monthly
mean and standard deviation of these values are also given in Table 4.4. It
should be emphasized that Table 4.4 gives the characteristics of a heated
channel, meaning that inflow temperatures were increased up to 15 C above the
temperature of the Mississippi River water. Longitudinal temperature gradients
and diurnal variations are not as large in a channel with ambient unheated
inflow.
Water temperature stratification in pools
It was found that the pool sections would frequently become thermally
stratified. Spot measurements of vertical water temperature stratification
in Channels 1, 5, and 8 were made using six thermistors attached vertically to
a point gauge and connected to a telethermometer. The point gauge was attached
to a cross bar which spanned the width of the channel. Vertical adjustment of
the point gauge enabled measurement of a complete temperature profile. An
example of the vertical stratification measurements is shown in Fig. 4.18.
Typically a thermocline would form at about 60 to 70 cm from the surface. The
lowermost 15 to 20 cm would have strong temperature gradients. Additional
measured water temperature profiles can be found in Appendix B of Hahn, 1978.
A frequency analysis of temperature stratification in pools 6 and 12 of
Channel 1 was made for the calendar year 1976. The continuous temperature
records of a vertical array of thermocouples were analyzed on a day by day
basis. Daily tabulations of the number of hours of stratification in the
ranges'0.5°C < AT < 1.0°C, 1.0° < AT < 2.0°C,.and AT > 2.0°C were made. AT
is the vertical temperature differential. The magnitude of AT was determined
from the temperatures measured by the top and bottom thermocouples in each
pool. If stratification occurred for one hour or more in a given pool, the
day was classified as a day with stratification. If the total number of hours
of stratification in a particular range was less than one, zero hours were
recorded. Stratification was observed somewhat more frequently in the down-
stream pools than the upstream ones. Water reaching the downstream pools has
had a longer exposure to solar radiation than water in the upstream ones. An
annual summary frequency analysis of the stratification data is given in
Table 4.5. A breakdown by months is shown in Figs. 4.19 and 4.20.
52
-------
TABLE 4.4. MONTHLY MEAN AND STANDARD DEVIATION (a) OF DAILY MINIMUM
AND DAILY MAXIMUM LONGITUDINAL (SPATIAL) TEMPERATURE DIFFERENCE,
AND OF DIURNAL (TEMPORAL) TEMPERATURE DIFFERENCE AT THE
OUTFLOW FOR A HEATED CHANNEL 1, DECEMBER 1975 - NOVEMBER 1976
Month
— -
December (1975)
January (1976)
February
March
April
May
June
July
August
September
October
November (1976)
Minimum
Longitudinal AT
' -Mean
5?f -
6.0
4.6
3.2
1.9
1.8
2.7
1.7
0.9
1.7
0.9
1.0
"- <°c)
" 2.0
1.5
1.3
2.7
1.8
1.3
1.5
0.9
1.1
0.7
1.0
1.0
Maximum
Longitudinal AT
Mean
" 8.0
9.5
7.7
6.9
5.6
5.7
4.3
4.0
4.5
3.5
3.4
5.7
oa
2.5
1.8
1.7
3.1
1.6
1.4
1.6
0.9
2.0
0.9
2.6
3.9
Diurnal
AT at Outflow
Mean
1.9
2.3
2.5
3.4
5.3
4.4
4.0
3.5
3.7
2.9
2.9
2.5
a
0.9
1.0
1.3
1.6
1.9
2.2
0.0
0.8
1.5
1.2
1.2
2.4
TABLE 4.5. FREQUENCY ANALYSIS OF TEMPERATURE STRATIFICATION IN CHANNEL 1
(January 1 - December 31, 1976)
Pool
NO.
- % days with stratification. ._ _6__.
12
"-" %- hours' with- stratification- - 6" -
12
\
o
> 0.5 C
48
34
-13 - -
11
/ertical AT
o
> 1.0 C
. 29- __
26
6 -
8
o
> 2.0 C
9
16
_ 2 ' ,., "...
3
53
-------
25
20
£
o
-------
100
80
60
40
20
o
o
0
5 ioo
80
60
40
Pool 6
Pool 12
D J
1975
M J J
1976
Fig. 4.19. Per cent days with stratification AT > 1.0 CO
and AT > 2.0°C A in Channel 1, pool 6 (top)
and pool 12 (bottom). Average channel flow
rate was 0.0175 m3/s until August 1976 and
0.0317 m3/s afterwards.
55
-------
o
50 n
40
30
20
10
0
£ 50
40
30
20
10
I I
Pool 6
Pool 12
D J
•"- 1975 -
M
M
J J
197-6
0-
Fig. 4.20. Per cent hours with stratification T > 1.0 C O
and T > 2.0°C A in Channel 1, pool 6 (top)
and pool 12 (bottom) . Average channel flow
rate was 0,175 m3/sec until August 1976 and
0.0317 m3/sec afterwards.
56
-------
SECTION 5
A STUDY OF LONGITUDINAL DISPERSION IN THE MERS FIELD CHANNELS
5.1. BACKGROUND ~. - - . . _ _
The flow field in the MERS" channels-is vecy non-unS-form and" mean flew^—--"•
velocities are very small. Water moves through the MERS channels in a way
far different from the movement of a piston in a cylinder (plug flow). Water
may move faster in the center of a riffle or a pool than at its edges. Trans-
verse and/or vertical mixing takes place between the faster and the slower
water masses. Vegetation is present in pools and riffles in the summer.
There is a weak jet effect at the transition from the riffles to the pools.
In the pools water masses may be entrapped near the banks. There can be
wind and there is intermittent stratification in the pool. As a result of
all these interacting processes, heat or material contained in the water will
spread out longitudinally along the channel. It was believed that the effect
of longitudinal spreading had to be incorporated in the water temperature
model.
Since a detailed study of the very complex low velocity channel nydro-
dynamics was not part of the research objective, the water temperature simula-
tion could be based on two conceptually different models.
(a) a tanks in series model with an adjustable residence time
parameter
(b) a.uniform .channel flow model._with-a longitudinal dispersion .
term - - - - -
The latter approach was chosen arbitrarily because it was believed and
later verified by the model results that the ooserved channel water tempera-
ture dynamics could be well predicted by a diffusion equation. The term
longitudinal dispersion coefficient is used herein to describe a bulk coeffi-
cient which lumps the effects of all the processes mentioned above. It
might also have been termed a bulk or hydrodynamic longitudinal diffusion
57
-------
coefficient. To use separate coefficients for riffles and pools was
con-sadered but the frequent transitions and short lengths of these
elements made it advisable not to use this approach.
The use of a diffusion equation to describe longitudinal mass trans-
port in a channel is only justified after the initial convective period
as shown by Fischer (1967). The MERS channels begin with a pool into
which water is discharged from a weir producing a very well mixed up-
stream condition.
The mean residence time of the water is anywhere from 10 to 30
minutes depending on flow rate. The convective length estimated from
the relationship
with i 1.8 m for a pool
h 0.8 m for a pool
n 0.1 (0.06 to 0.3}
g 9.8 m/sec
gives a value of approximately 5 m or far less than the length of the
first pool (30 m). It can therefore be considered that the convective
length does not reach beyond the length of the first pool.
A study of longitudinal dispersion in the MERS field channels was
made for two reasons. (1) A longitudinal dispersion coefficient was
needed in the unsteady channel water temperature model to be discussed
in Section 7, and (2) the size of the MERS channels is intermediate
between laboratory flumes and full-size canals; size,is an important
factor for longitudinal'dispersion and no"studies had ever been made on
a system of that intermediate size.
Longitudinal dispersion was first studied analytically in pipe flow
by Taylor (1954) and was later applied to open channel flow by Elder (1959)
For a straight, two-dimensional channel Elder obtained
Dr = 5.93 h u. (5-1)
Lj *
58
-------
where h is channel depth and u^ the shear velocity. Experimental verifi-
cation of Eq. 5-1 by Elder (1959), Fischer (1968), Thackston and Krenkel (1967),
and Sayre and Chang (1968) has indicated values of the constant from 5.3 to
15.7 for various flow and straight channel conditions.
Elder's equation, however, is not generally applicable to natural streams
and larger channels in which dispersion coefficients are much larger than
those predicted by Eq. 5-1. Fischer proposed that in natural streams velocity
profiles and transverse mix-ing are dominant in determining longitudinal
dispersion. Thus, the longitudinal dispersion coefficient depends upon
shoreline configuration, bends, and bed formation, as well as the stream
velocity (see e.g. Day, 1975).
A considerable amount of research on longitudinal dispersion in laboratory
flumes and in natural streams has been reported since the above cited early
studies. A recent summary of Fischer et al (1979) lists some of the investi-
gators who have contributed analytical or experimental information. Prior
reviews were given by Ditmars (1974) and Fischer (1973), El-Hadi and Davar
(1976).
The prediction of longitudinal dispersion for a stream of given geometry
without some form of field measurements of velocity profile or dye dispersion
is still a difficult problem. Relationships based on geometrical and hydrologic
stream characteristics were proposed, among others, oy Bansal (1971), McQuivey and
Keefer (1974), Jain (1976), Fischer et al (1979), and Liu and Chen (1980).
These relationships are of a semiempincal nature, e.g.,
Liu and Chen: D. = (0.27 to 1.25) u_ B2/h
LJ *
-- Fischer: • — O -- 0.011 U 2 ~B2/hu. - -~ —"
-- ]_, - *
where" D = longitudinal dispersion coefficient, "" ' ~ ~
L
U = mean velocity,
B = stream width,
h = mean stream depth, and
u^ = shear velocity.
59
-------
The channels investigated in this study do not have the accentuated lon-
gitudinal mixing characteristics of natural streams caused by meanders and
large widths. They are also very different from laboratory channels. Abrupt
changes in cross sectional areas from "riffle" to "pool" sections occur every
30 meters (100 ft). The parabolically shaped pool sections are 4 m wide and
approximately 1 m deep. The repeated expansions and contractions in flow
would be expected to affect longitudinal dispersion.
-- At elevated velocities longitudinal mixing is enhanced by ]et effects of
the-flow-from-the raffle sections into the. pool sections. Flow streamlines
which could be expected from such }et mixing are-sketched in Fig. 5.1. Lomax
and Orsborn (1971) studied the effects of jet mixing in a small, experimental
pool of circular shape and constant depth. A schematic of a measured outflow
concentration curve is shown in Fig. 5,2. Inflow started at t = 0 with a
concentration C in the pool. Concentration of the inflow was C = 0.
The pools at the Monticello Field Station have a length of nine times
the width and very low flow velocity. This confines ;jet mixing effects to
the upstream portion of each pool and reduced its significance. Jet mixing
was therefore evaluated as part of the longitudinal dispersion coefficient.
Dispersion in the MERS channels can also be affected by wind mixing,
plant growth, stratification and/or natural convection in the pools because
the mean flow velocities in the pools are only on the order of 1 cm/sec. The
measurements and the analysis integrated all those effects.
5.2. METHOD OF LONGITUDINAL DISPERSION COEFFICIENT DETERMINATION FROM
TRANSIENT WATER TEMPERATURES
The contribution of longitudinal dispersion to longitudinal water tempera-
..ture profilss.-increases with the transient character of water temperature.
The effect of longitudinal dispersion "is'most evident "when" a tempeFature pulse""
or a front of elevated or depressed temperature moves downstream. In this
section, temperature fronts will be related to longitudinal dispersion theory
and used to estimate longitudinal dispersion coefficients in the MERS experi-
mental channels.
The longitudinal dispersion coefficient for natural streams is often
determined from an instantaneous infection of a conservative tracer. The
one-dimensional mass transport equation for the tracer is
60
-------
Pool
Riffle
Fig. 5.1. Schematic streamlines at transition from riffle to pool in
Monticello experimental channels without (above) and with
(below) transverse winds.
Ideal longitudinal dispersion curve
Observed jet mixing in pool
Fig. 5.2. Scnematic outflow concencration curve for net mixing
in experimental pool (after Lomax and Orsoorne, 1971)
61
-------
with solution for an instantaneous injection
where C =• tracer concentration (M/L ),
M = mass of tracer injection (M),
A = channel cross-sectional area (L ),
U = cross sectional mean flow velocity (L/T),
x = distance from point of injection (L),
t = time from moment of injection (T), and
D. = longitudinal dispersion coefficient (L /T).
L
The use of Eq. 5-2 assumes that the shape of the tracer cloud will be
Gaussian. By the method of moments, or by a routing method, the longitudinal
dispersion coefficient between two points may be found. More comprehensive
details of D determination from tracer pulses in natural streams are given
by Godfrey and Frederick (1970), among others.
In the MERS experimental channels water temperature was used as a tracer
to estimate longitudinal dispersion coefficients. The most readily obtainable
trace, however, was not a heat impulse but a temperature front. Such fronts
could be obtained by turning on or shutting off the heat exchanger in the feed
line to Channel 1. (See Fig. 5.3 for examples.)
A longitudinal dispersion coefficient can be- accurately determined by
"-routing the temperature front downstream and adjusting D to minimize the
L
error in measured and predicted water temperatures. The analytical solution
for a temperature front assuming no sources or sinks is
9=| Erfc (x-Ut)/(4DTt)"'"| (5-4)
62
-------
10
O1
uf
_ 6
0°
OJ
ex
I
I
I
50 100 150 200 250
Distance Downstream (m)
300
350
400
Fig. 5.3. Temperature fronts in MERS channel 1 on November 17, 1976. Front
! inception is at t=0.
-------
where 0 = „}_[ _ *}„( and (5-5)
.2
v^T f
' ~ J
Erfc(g) = 1 - rr I d"Q dQ
Erfc is the complementary error function. Measured values of 9(x,t) were
compared to values computed by Eq, 5-4., The D value which minimizes the
error between computed and observed Q is the "measured" D for the test.
Heat transfer may also affect the longitudinal dispersion coefficient
calculated from a temperature front. The order of magnitude of this error was
estimated by computing a temperature front with and without heat transfer using
typical values of channel depth and width, flow rate, and longitudinal disper-
sion coefficient, with upper limiting values for the bulk surface heat transfer
coefficient, and the upstream temperature front (Gulliver, 1977). The computed
results indicate that surface heat transfer has no appreciable effect upon the
longitudinal dispersion coefficient determined by the routing method. Fig, 5.4
illustrates why this should be expected. Longitudinal dispersion can be di-
rectly visualized as the spread between the 0.16 and 0.84 values of dimension-
less concentration Q. Surface heat transfer reduces the values of 0 but not
appeciably the longitudinal spread.
5.3. FORMULATION OF DIMENSIONLESS LONGITUDINAL DISPERSION NUMBERS
Many investigators have expressed longitudinal dispersion by a dimension-
less number, D /u^ h , where u^ is bottom shear velocity. For the MERS
channels flow rate, rather than water surface slope, is more accurately mea-
sured. With
and
where f. = bottom friction factor for entire channel,
b
h = mean depth of channel, and
n = Manning's coefficient for the channel,
64
-------
084 —
Q
en
ui
with surface
heoi tronsfer
without surface
heot trorisfer
LM 04 m/sec
2
D(_- 28m /sec
x= 427m
-2 -10 -I
Ks=45 col cm doy C
T(CO)-T(o) = 20°C(at x=0)
8
CD
016 -
234
Time from Pulse Instigation (hr)
Fig. 5.4. Computed temperature fronts (normalized) which occur at
location 17 in Channel 1 with and without heat transfer.
-------
one derives for the dimensionless expression
\g n /
where Q = flow rate and B = mean surface width. For constant mean channel
depth h and constant Mannings 'n'
Dt, B DL
D^=_L_ « -_£-. _ (5_7)
For a stream or channel of variable cross section, Bansal (1971) proposed
a dimensionless local dispersion coefficient
V = I ITI - An L" (5-8,
B o h U AQ
where U = local cross sectionally averaged flow velocity,
U = average flow velocity, U, over stream reach,
h = average depth over stream reach,
B = average width over stream reach,
A = average cross sectional area over stream reach,
A = local cross sectional area, and
Q = flow rate.
D* and D * represent different expressions of the dimensionless
5
longitudinal dispersion coefficient for a channel of variable morphology.
5.4. RESULTS OF TEMPERATURE ROUTING TESTS
Five recorded temperature fronts were"routed through the MERS" channels
to determine the longitudinal dispersion coefficient DT. These five events
L
were selected by three criteria: 1) the temperature drop or increase had to
be strong enough to be accurately traced as it moved downstream, 2) external
conditions such as stratification, ice cover, or algal mats on the water sur-
face which would alter the computed dispersion coefficient had to be absent,
and 3) large diurnal variations or other transient water temperatures except
for the temperature front had to be absent. The data for 0, channel morphology,
66
-------
and travel times are given in Appendix A. Temperature records from 5, 6, or
7 stations were available. Values of 0 computed from Eq. 5.4 were compared
to measured 0 values (for 0 > 0.05 only). Values of D. were chosen
— Li
successively until the error between computed and measured 0's was minimized.
Below this range estimates of 0(x,t) were considered to be more charac-
teristic of entrapment rather than longitudinal dispersion.
A dispersion coefficient, representative of the entire channel, was
calculated for each experiment using the temperature front data from all sta-
tions and the corresponding equation 5-4.
Each set of data was reduced for four D formulations: (1) D. = constant,
— — " _ "
(2) Dr = D* Q/B, (3) D. = D* QA/(A B), and (4) D = D*Q/B . Results are given
L Ll D Ll
in Table 5-1.
The D formulation which is most applicable may be selected by computing
LI
and comparing ratios of the standard deviation of the fitted coefficients to
the mean of the five tests. These ratios represent the relative spread in
fitted coefficients. Numerical values are shown in Table 5-1, last line. The
dimensionless coefficients D* and D* represent about equally well the
a
longitudinal dispersion for variable geometry and both are better than the
constant DT. The choice between D* and D* remains purely arbitrary. The
L o
reason for both coefficients giving equally good results will be examined in
the next section.
Finally, the small standard errors and the consistency of the results
indicate that the temperature front approach to determining longitudinal
dispersion coefficient is valid and practical.
5.5. EFFECTS OF VARIABLE CROSS SECTIONAL AREA ON LONGITUDINAL DISPERSION
The subdivision of a channel into a series of-pools of long residence
time connected by riffles of short residence time makes the longitudinal
dispersion almost entirely dependent upon the mixing processes in each pool
and the number of pools. The processes are complicated and include:
« jet-like mixing at the transition from riffle to pool
» flow through the pool guided by vegetation and wind
« stratification effects.
It was not part of this project to analyze the motions in a pool in detail.
67
-------
TABLE 5-1. LONGITUDINAL DISPERSION IN MERS-CHANNELS FROM TEMPERATURE FRONT DATA
o>
co
Date
of
Experiment
11-17-76
11-18-76
11-22-76
12-06-76
12-14-76
Mean
S.D.
S.D.
Mean
*Standard error
Date
of
Experiment
11-17-76
11-18-76
11-22-76
12-06-76
12-14-76
Mean
r\
Flow L
3 2
(m /sec) (m /sec)
.0349 0.100 (.0142)*
.0382 0.095 (.0084)
.0373 0.100 (.0152)
.0394 0.124 (.0172)
.0428 0.128 (.0106)
0.109 (.0131)
0.015
0.140
•
TABLE 5-2. AVERAGE POOL AND RIFFLE
D. (l)=Const DT (2) =
L Li
DL (pools)
2 2
m /sec m /sec
0.100 0.099
0.095 0.093
0.100 0.100
0.124 0.124
0.128 0.124
.109 .100
D B
L
Q
9.3 (.0143)*
8.0 (.0082)
8.7 (.0153)
10.3 (.0168)
9.6 (.0108)
9.2 (.0131)
0.88
0.095
D FROM TEMPERATURE
L
D* Q/B
5 (riffles) D
L
2
m /sec
0.174
0.164
0.176
0.219
0.217
.190
D A B D
~ * L
B fiA
4.8 (.0136)* 7
4.4 (.0095) 6
4.9 (.0153) 7
5.6 (.0182) 8
5.3 (.0103) 7
5.0 (.0134) 7
0.46 0
0.093 0
FRONT ROUTING TESTS
JL
DL(3) = DB QA/(A B)
v (pools) Dr (riffles)
Ll Li
2 2
m /sec m /sec
0.102 0.026
0.102 0.026
0.114 0.029
0.135 0.034
0.139 0.035
.118 .030
B
L
Q
.55
.56
.07
.30
.88
.47
.68
.091
-------
The measured longitudinal dispersion coefficients D* (Table 5.1) are
of the same order as those reported by Fischer (1973) for laboratory channels.
This suggests that the combined effects of pool mixing and pool number were
comparable to the effects of advective motion and transverse mixing in a
straight channel.
That the pools dominate the effective longitudinal dispersion can be
seen from an expression for the flux T of material in an open channel.
3C
T = 0° - ^L i^ <5^9)
The advective flux QC in riffles and pools are comparable since Q and C are
the same. The dispersive flux is, however, a function of cross sectional
area A and D . The cross sectional areas of the pools are three to four
L
times those of the riffles (Table 4.2). The dispersive transport in the
pools is therefore dominant over the dispersive transport in the riffles.
It is for this reason that in Table 5.2 for all five experiments, the disper-
sion coefficients in the pools are comparable, whereas the riffle values are
not. (The values in Table 5.2 are computed from D* and D* values in
B
Table 5.1 using riffle or pool values of B and A.)
5.6. COMPARISON WITH OTHER LONGITUDINAL DISPERSION MEASUREMENTS
If the mean flow velocity U = Q/A in D* is replaced by the shear
velocity u^ , a comparison of these results with the dimensionless longitudinal
dispersion coefficients, D/(u.h), cited by Fischer (1973) is possible. Field
L *
measurements on the MERS experimental channels have indicated that u*/u=0.55
for these channels. The hjgh value of u^ /U is misleading, however, when
compared to the longitudinal dispersion coefficient. In the MERS channels
~ffio'st~of the «ater~surface-"drop-oecurs- in the r if fie-sections.--In fact, - - -•
water surface slope in the pools is so-smalt that it-cannot- be- measured.--- - - - .
As shown in the previous section, longitudinal dispersion occurs predomi-
nently in the pool sections. Thus the dimensionless coefficient, D /u^h ,
compares u representing surface slope in the riffle sections' and D
* L
representing longitudinal dispersion in pool sections. One would expect
this value to oe smaller than for most natural channels. For the mean value
of D* = 9.2,
69
-------
This value can be compared with others cited by Fischer (1973) . It is found
that the value given here is in the range of rectangular laboratory channels
with smooth sides, in which longitudinal dispersion is much less than that
found in natural channels. This is noteworthy since the geometry of the MERS
channels is very different from a laboratory flume.
5.7. EFFECT-- OF ROUGHNESS- (MANNING-1 S 'n-J-)- ON- LONGITUDINAL DISPERSION- IN _ __
THE MERS CHANNELS ..... ~
The mixing processes in the pools of the MERS channels are not very de-
pendent on Mannings 'n1 values. Vegetation may be a factor, but it is doubtful
that the relationship for shear flow in channels have much significance. Changes
in Mannings 'n1 caused by seasonal variations in vegetation are therefore not
expected to have much significance on the longitudinal dispersion coefficients.
The longitudinal dispersion coefficients reported in Table 5.1 are for
late fall/early winter conditions when vegetation was still present, although
not as extensive as in mid-summer. Mannings 'n1 values in November/December,
when the measurements were made, are probably close to a seasonal average.
70
-------
SECTION 6
HEAT TRANSFER ACROSS THE WATER SURFACE IN THE MERS FIELD CHANNELS
6.1. BACKGROUND
Net heat transfer across an air-water interface is the sum of several
components. The source term in Eq. 1-1 may'be written as
S = H - H. - H -H (6-1)
S I e c
where S = net heat flux per unit surface area,
— 2 —1
H = net shortwave solar radiation entering the water surface (EL T ),
H5 = net long wave radiation leaving the water surface,
H = energy leaving the water surface due to evaooration, and
e
H = energy convection from the water to the air.
c
Many different empirical, semi-empirical, and deterministic relationships to
determine the terms of Eq. 6-1 have been developed. Paily, Macagno, and
Kennedy (1974a)provided a detailed summary. Only the equations to be used in
this application will be given here. The parameters required to compute the
net heat transfer are temperatures (water, air, and dew point), barometric
pressure, solar radiation, cloud cover, and wind velocity. In this study
the effect of wind direction on heat transfer in the MERS channels was also
investigated.
A separate analysis of each term in Eq. 6-1 results in. a net heat flux
with nonlinear dependence on water temperature. In order to reduce complexity
in analytical derivations and numerical computations, methods of computing heat
flux as a linear function of water temperature have been proposed by Edinger,
Duttweiller, and Geyer (1968), Brady, Graves and Geyer (1969), Dingman and Assur
(1967), TVA (1968), Yotsukura, Jackman, and Faust (1973), Paily, Macagno, and
Kennedy (1974b), Jobson and Yotsukura (1972) and others. The main drawback
- 71 "
-------
to the.linear models is that they are only valid over a limited water tempera-
ture range. Each of these methods may be expressed in a formulation first
proposed by Edinger, Duttweiler, and Geyer (1968). Edinger, et al introduced
a hypotehtical equilibrium temperature, T , and proposed that the net heat
flux across an air-water interface may then be expressed as
S = -Ks (Ts - TE) (6-2)
where. T is water surface temperature and K is a bulk heat transfer
S S
coefficient at the watex surface. Equation 6-2 indicates that net heat
transfer will be small if (a) water surface temperature is close to equilibrium
temperature and/or (b) K is small (calm wind, high relative humidity, etc.).
S
T and K are indicative of water temperature response to weather conditions.
D S
The MERS experimental channels are typical of many small streams and
canals in that large diurnal variations in temperature occur. These diurnal
variations are of interest for three reasons.
(a) The diurnal temperature maxima and minima affect the channel eco-
system studied by other investigators.
(b) Currently available water temperature prediction models use
coefficients and time scales appropriate for deeper waters
(typically one day). Water temperature analysis at much
shorter time scales (3 hour) has not been much explored.
(c) Wind sheltering of the water surface in the MERS channels
by the banks is a strong possibility and needs to be
explored.
Thus, predicted temperature on less than daily time scales are required, and
the heat transfer equations and"tlme scale- of. daia input are chosen accordingly„
6.2. NET SHORTWAVE (SOLAR) RADIATION
Essentially all shortwave radiation is from the sun to the water body
(either direct or diffuse). In this study solar radiation is continually
measured as a weather parameter and net shortwave radiation is expressed as
the difference between (measured) incoming and (estimated) reflected radiation,
72
-------
or
H = H - H (6-3a)
s si sr v '
H = H (1-r) (6-3b)
O O J.
where H = incoming solar radiation (EL T ),
H = reflected solar radiation, and
r = total reflectivity o£ the water surface.
Anderson (1954) proposed an often used equation for total reflectivity of a
water surface,
r = A a8 (6-4)
where a = solar angle in degrees, and
A and B = constants depending upon cloud cover.
If time increments of one day or greater are of interest, the relations of
Koberg (1964) for H are more appropriate. Brady, Graves, and Geyer (1969)
calculated A and B as
Q *
A = 2.20 + -7— (C °'7 - 0.4)2/0.16 (6-5)
4 • U C
C °'7
B = -1.02 + !; , • + (C °'7 - 0.4)2/0.64 (6-6)
lo .0 r
where C = 1 - H /H = "Cloudiness ratio/" and
r si sm
H = solar radiation which would occur with clear skies.
sm
Ttie-se- ate-obv-iously -empirical-equations.— The -solar-angle- may-be .-determined
from, the.,following equations .as-J-is.ted,~in-PaiLy-_et..al.-il_^74^._._-._ ..—
sin a = sin $ sin fi + cos cos fi cos(h) (6-7)
where $ = geographic latitude (radians)
iS = declination of the sun (radians)
23.45 TT \2 IT ,,__ 1 ,. Bi
cos TTT (172-D) (6-8)
180 365
73
-------
D = day of the year (Jan. 1=1), and
,h = hour angle of the sun (radians)
= n (day hour - 12)/12 (6-9)
For this particular application, however, the average value of sin a over a
time increment is desired, rather than an instantaneous value. Thus, the
solar angle which will be used is
ct(rad) = sin (sin a) (6-10a)
where /. h2
sin a = . . / sin a dh
Vhl
... . . f sin h. sin h, }
= sin $ sin fi + cos $ cos 6 \ 2 - 1' ,,,„,_,
-- - - - - (6-lOb)
The value, of H may be estimated from the equations which are listed
Sfu
in Paily et al (I974a) .
H = H (1-C ) [a" + 0;5(l-a')l " (6-11)
sm so s '-
where C = a constant equal to the approximate value of some less significant
parameters (such as dust depletion) which are difficult to
determine,
= 0.01,
H = solar radiation incident on the atmosphere,
= I sin ct/R ,
° 2 2
I = solar constant = 2 cal/cm /min=2880 cal/cm /day, and
R = ratio of the actual to the mean distance from the sun to
the- earth. . - . .
1 + 0.017 cos|-^| (18.6-D)J
(6-13)
L •><•>-' J
a1 and a" are coefficients.
a1 = exp J m [- (0.465+0.130w)] [o.179+0.421 e~°*721jn] J (6-14)
a" = exp | m [- (0.465+0.134w)] [o.129+0.171 e~°'88m ]} (6-15)
74
-------
where m = optical air mass
P .P
a/ o
sin ct + 0.15 (a + 3.885)
a = solar angle in degrees ,
P /P = ratio of air pressure at location altitude to sea level air
3 O
pressure,
288 - 0.0065 Alt \5'256
288 j 15'17)
Alt = altitude in meters.
w = atmospheric moisture content (cm)
= 0.85 exp 0.110 + 0.0614 T , and (6-18)
o
T, = dewpoint temperature ( C) .
d
Although Eq. 6-4 through 6-18 are a rather long list to compute a simple
parameter, r, to input a typical value of reflectivity can create large errors
when time increments of less than oae day are used. For the MERS experimental
channels, r values from .04 to .2 have been computed over a day, and when
the importance of the H term is considered, these variations are too
large to be approximated accurately by "typical" or "average" values.
6.3. NET LONGWAVE RADIATION
Net longwave radiation is expressed as
H4 -°ew (Ts"- eaTa ' (6'19)
where T = water surface temperature ( K),
T = air temperature ( K),
cl
e = long-wave emissivity of the water surfaca (assumed equal
to water surface long-wave reflectivity),
= 0.970
e = emissivity of the atmosphere, and
3
a = Stefan-Soltzman constant.
75
-------
For atmospheric emissivity without cloud cover, £ , the Idso and Jackson
clC
(1969) formula will give accurate results for air temperatures above and
below freezing point.
= 1 - 0.261 exp -0.74xlO~4 T (°C)2 ] (6-20)
L Si J
ac
The Bolz formula is then used to find e <•
a
- . . ~*F~ = e " (KK. C-2) - (6-21)
«v~. ac— G-
wfieFe C" "= 'fraction cloud cover', and — --- .-"—---
c
K = a coefficient which depends upon cloud height.
The coefficient K varies between 0.04 and 0.25. A TVA (1968) study
recommends an average value of K = 0.17,
6.4. NET EVAPORATIVE AND CONVECTIVE HEAT TRANSFER
Evaporative heat transfer from a water surface may be expressed by the
relation,
HQ = pL(Wftn)z (esw - e^) (6-22)
where e = vapor pressure of the air at height z,
32
e = saturated vapor pressure at water surface temperature,
sw
Wftn = a wind function using wind velocity at height z,
L = latent heat of vaporization for water (E/M), and
p = density of water (M/L ).
...Saturation vapor pressure a_t. any air temperature. 3L(_JC)__may;Jpe:.comp_uted_pyer- _
water.. by_t_he .Magn.us-Tetons formula, (see. _Murray, 1967).
<*, - ..1.7, «XP ["•26r!£.;62"-16>]
Atmospheric vapor pressure is computed from relative humidity, RH
ea = IbT esa (6'24>
76
-------
The latent heat of vaporization is
L = 597.31 - 0.5631 T (6-25)
with L in cal/g and T in C.
The convective heat transfer from an air-water interface, when evaluated
according to Bowen (1926) may be expressed as
P (mb)
Hc = °'61 "WOT pL Wftnz (Ts - **> (6-26>
where T is air temperature at a height z above the water surface, p is
in mb and Wftn is the same as that for evaporative heat transfer.
Z
A number of empirical wind function formulas(Wftn) have been developed
z
for various conditions. A formula used by many investigators, e.g. Marciano
and Harbeck (1954), for natural water bodies is
(Wftn) = a + b W (6-27)
Z Z
where W = wind velocity at elevation z above the water surface, and
a and b = empirical constants.
Brady, Graves, and Geyer (1969) used the formula,
(Wftn) = a + b W 2 (6-28)
for heat loss from power plant cooling ponds. Shulyakovskyi (1969) incor-
porated a natural convection term in Eq. 6-27 to describe accentuated heat
loss from thermally loaded water bodies.
(Wftn)z = a + b Wz + C(A6 )1/3 , (6-29)
where 8 = virtual temperature or the temperature of dry air at the
same density as under tne given conditions, and
A8 = the difference in a virtual temperature between air at the
water surface and at 2 m height.
The difference in virtual temperature is, in effect, a theoretical tempera-
ture difference which represents differences in density. The relation for
A8 is
v
77 .
-------
A9 = T (1 + 0.378 e /p ) - T (1 + 0.378 e /p ) (6-30)
v s s 3 a a a
Ryan and Hacleraan (1973) excluded the constant term in Eq. 6-29 because of
the large number of studies at natural water temperatures which found a = 0.
They fitted laboratory natural convection and cooling pond data to the formula
(Wftn) = b W + c(A6 )1/3 (6-31)
Z Z V
Constants for the above equations have been determined for large water
bodies. For smaller water bodies like the MERS channels, the forced
convection constants (a and b) need to be evaluated separately. In the next
section the constants in Eq. 6-27, 6-28, 6-29, and 6-31 will be determined
from water temperature and weather data at the MERS and compared to results
for water bodies with large surface areas.
6.5. DETERMINATION OF WIND RELATED HEAT TRANSFER COEFFICIENT FROM MEASURED
STEADY-6TATE LONGITUDINAL TEMPERATURE PROFILES
General Procedure
With constant inflow rate and inflow temperature and quasi-steady
weather conditions, the longitudinal temperature profile in the channels
will approach steady state. Conversely, when water temperature records
indicate steady-state and weather conditions are nearly constant, the
longitudinal water temperature profile in the channels may be used to deter-
mine coefficients in the steady-state equation of the longitudinal water
temperature profile. If longitudinal temperature gradients are not large,
and if Eq. 6-2 is used for the source term, water temperature is predicted
by the solution to Eq. 1-1 for zero dispersion.
f - - f - . - -
8 T(x) ~ TE
where — = — —
9o o ' TE
T = equilibrium temperature,
£•
Io = T(x=0) ,
U = bulk channel flow velocity,
78
-------
m = Ks/Pcp h '
K = bulk surface heat transfer coefficient,
S
h = channel depth, and
c = specific heat of water.
From a measured longitudinal temperature profile and a computed equilibrium
temperature, the value of K may be determined by Eq. 6-32 as
Q
- - -- - - - - c~ Ti U ~ ~ ~~
- - . KS = -P-^T- lnf (
o
S Q e
or K = -p -^— In f- (6-34)
s IT y
Bx o
where Q = channel flow rate, and
B = mean channel surface width.
Errors in Computed Heat Transfer Coefficients
For a large measured longitudinal temperature gradient the contribution
of longitudinal dispersion, D , to water temperature gradients may be
u
significant. For this case the following solution to Eq. 1-1 with D > 0
Lj
must be applied.
(6-35)
The error Y(%) in determining K from Eq. 6-34 by assuming D =0 is
S L
In 9/9 (D =0) - In 9/9_ (D > 0)
*'
. For-the calculations-it is convenient to-rewrite Eq. 6—35 as,
In p2 (DT > 0) - - - I (1 + K D)*' - - 1 I (6-37)
V Ll """ ' *
O D
where D = 2 D/U
L
K = 2 m/U .
79
-------
For a given per cent error Y, equations 6-34, 6-36, and 6-37 may be
combined to give
pL\l JL . O _ o ,6 3g
UOO/ 100 4 ( JB
Eq. 6-38 may be solved for Y. Typical input parameters for K and D at
the MERS channels are
0 < D < 1.0 m/sec
Lt
9 x 10~5 < m < 2 x 10~4 sec"1
.05 < U < .2 m/sec
It is found that Y does not exceed 2 per cent. Longitudinal dispersion
therefore will not significantly affect the surface heat transfer coefficients
computed from steady-state longitudinal water temperature profiles, while
ignoring longitudinal dispersion.
Equation 6-34 also implies that the longitudinal temperature profile has
achieved steady-state. True steady-state can be approached but is rarely,
if ever, achieved precisely. Therefore an estimate of the unsteady character
of the longitudinal temperature profile and its effect on computations of K
s
is also needed. With D =0, Eq. 1-1 may be written as,
L
If ^BbS-
or
NS + ss
- m f- (6-40)
9o
where NS = 3T/3t is the unsteady contribution and
SS = (Q/Bh) 3T/3x is_ the steady-state contribution.
The per cent error Z in 9/9 created by assuming NS = 0 is
__ _ ,
100 ~ \ SS+NS
80
-------
If, for each assumed steady-state Q/Q , the location with the maximum 3T/3t
is used to compute NS, and SS is taken as the average over channel length,
Eq. 6-41 will give an estimate of the largest possible error in K computa-
tions.
Determination of Wind Function Coefficients
The bulk surface heat transfer coefficient K (cal cm" day" °C~ ) is
given by Edinger and Geyer (1965) as
KS = 9.256 + pL(Wftn){& + 0:61). . (6-42)
where Wftn = wind function (cm mb day" ), L = latent heat of
vaporization (cal/g) , and g (mb/ C) -• mean slope of the saturated
vapor pressure versus temperature curve between water surface and dew point
temperatures. Brady, Graves, and Geyer (1969) suggested the formula
3 = .4604 + .0197 T + .001585 T2 (6-43)
where T = (T + T )/2 (°C)
T = water surface temperature (°C), and
T, = dew point temperature (°C) .
Equations 6-37 and 6-29 may be combined to give the wind function for a steady-
state longitudinal temperature profile.
= pL(Wftn) = — - ~ - cal cm~2 day"1 (6~44)
Equation 6-44 is not dependent upon channel depth, indicating pool stratifi-
cation does not affect computations of F . This is so because stratifica-
tion decreases' not only the effective channel depth, but also residence time
as long as the flow rate is constant.
Using Eqs. 6-34 and 6-44, K and F were computed for 47 cases of
near steady-state longitudinal temperature profiles. The parameter NS/SS
and the maximum per cent error in K_ due to unsteady water temperatures
were also computed for each ease. The data and computational- results of all
cases is given in Appendix B, Table B-l. A sample calculation is also given
in Appendix B.
81
-------
A least-squares fit to determine the coefficients in Eqs. 6-27, 6-28,
6-29, and 6-31 was undertaken with all 47 cases. In addition, the data
were also fitted to the equation
F = nl,(Wftn) = b W + 5.52(Ae ) (6-45)
w z z v
z
where c = 5.52 was found by Ryan and Harleman (1973) . Equation 6-45 was
included because natural convection should not depend strongly on water
surface area, and Ryan and Harleman's (1973) results may be applicable to
the MERS channels.
The following five equations result from a least squares fit on computed
F values from the MERS field channels:
1. F (cal cm"2 day'1 mb'1) = a + b W (m/sec) (6-46)
w y
with a = 17.5
b =• 2.40
and a standard error = 3.79
2. F = a + b Wn2 (6-47)
w 9
with a = 21.18
b = .316
and a standard error = 3.74
3. F = a + b W. + c(A8 )1/3 (6-48)
w y v
with a = 9.17
...... 7 ." b ~= 2768 "~" ".""'_";• ---- •• ----- -------
c = 2.78 ' ' ------ ' - ---- • -
and a standard error = 3.72 " "" ......... "
4. Fw = b W9 + c(A8v)1/3 (6-49)
with b = 3.10
c = 5.68
and a standard error = 3.81
82
-------
5. FW = b Wg + 5.52 (A9V)1/3 (6-50)
with b = 3.21
and a standard error =3.77
Equations 6-46 and 6-47 are compared to experimental (field) values of
F in Figure 6-1. Three ranges of the parameter NS/SS are indicated by
different symbols. The scatter in Fig. 6-1 is not uncommon, as indicated
e.g. by similar plots by Brady, Graves, and Geyer (1969). Equations 6-48,
- 1" /O
6^49, and 6-50 are compared_to computed values of F - 2.78(A9 ) _ ,
FW - 5.68(A9V)1/3 , and FW - 5.52(A9v)1/3 in Figs. 6-2, 6-3, and" 6-4,
respectively. Comparison of the standard errors of Eqs. 6-46 through 6-50
and comparison of the individual curves with F data indicates that all
five equations give a remarkably similar fit to the data. This indicates
that for a given type of water body, in this case narrow open channels with
constant thermal loading, any of the five equations may be applied as a wind
function.
Paily, Macagno, and Kennedy (1974a)compiled a list of wind function
formulas given by various authors for different water bodies. Some of the
wind function formulas which compare with those used here are given in Table
6-1. Table 6-1 gives (Wftn), = F /(pL) as a function of wind velocity at
/ W2
2 m above the water surface. Equations 6-46 through 6-50 were converted
to the form of Table 4 and added on the bottom.
Wind velocities at 9 m above the water surface were related to wind
velocities at 2 m, W , by the formula
W2 = Wg(2/9)'3 (6-51)
as suggested by Paily,--Ma'eagno,~ and"Kennedy (1974a)for wind over land.
Because of the- small size of. jbhe Monticello channels',, the ve.rtical wind
velocity profile is essentially that over land. A comparison of the relation-
ships found in the MERS-channels and from previous investigations as shown in
Table 6-1 can now be made.
A number of the investigators (Kohler, Marciano and Harbeck, Harbeck,
Hughes, Turner, and Ficke).listed in Table 6-1 used a wind function formula
83
-------
Q.
O 0
-------
oo
U1
30
o
T3
E
o
8 20
CD
c\j
i
5
10
p NS/SS=0
O 0
-------
2P
24
-t
CO
. ' 0
""TT
n NS/SS=O
O 0
-------
CD
-J
o
XI
o
u
28
24
20
'6
CD
<
OJ
m
IT)
0
D NS/SS=0
O 0
-------
TABLE 6-1. WIND FUNCTION FORMULAS DETERMINED BY VARIOUS
INVESTIGATORS*
Investigator
Kohler (1954)
Marciano and
Harbeck (1954)
Harbeck (1958)
Rymsha and
Doschenko (1958)
Hughes (1964)
Turner (1966)
Wind Function Formula (Wftn)
(cm/day-mb)
0.005175+4.995x10 W_
0.01134
0.01309 W.
0,0209+9.107x10 (T-Ta2)+0.01018 W
0.008864
0.02045
Water Body
Lake Hefner,
barge and
midlake
stations
Lake Hefner
Lake Mead
Winter Time
Rivers
Salton Sea,
California
Lake Michie,
No. Carolina
Brady, Graves,
and Geyer (1969)
Shulyyakovskyi (1969)
Ficke (1972)
Ryan and Harleman
(1973)
0.0239+0.00147
1/3
0.015+0.0094 (A0v)+0.0112
0.0125 W,
0.00934 (A6v)1/3+ 0.01107
Power Plant
Cooling Ponds
Pretty Lakes,
Indiana
Cooling Ponds
and Laboratory
Data
Stefan, Gulliver
Hahn & Fu
0.0296+0.00637
0.0358+0.00132
1/3
0. 0155+0. 0047(A9v)+0. 00711
0.0096 (A9v)1/3+0. 00832 W2
0.00934
. 00852
MERS channels
6-46
6-47
6-48
6-49
6-50
*Listed by Paily, Macagno, and Kennedy (1974a).
88
-------
like Eq. 6-46 for lakes with natural water surface temperatures. These inves-
tigators found a very small or zero constant term. The large constant term
in Eq. 6-46 is believed to reflect natural convection due to the artificially
elevated water surface temperatures in the Monticello Field channels studied
herein.
Brady, Graves, and Geyer (1969) used a wind function similar to Eq. 6-47
forT^at loss- from power plant cooling" ponds, given in Table 6-1. Their
~v~ _ ^ — — - .. _ .
constant term is smaller than, that of Eq.-6-47.- The explanation for this is
similar to that for^tfie higher constant term in Eq. 6-46. The-virtual
temperature difference (A6 ) in the cooling ponds studied by Brady et al
varied from 5 to 12°C. The range of A6 in the open channel studied here
was 7 to 35 C, resulting in stronger natural convection effects.
Equation 6-48b has the same form as that used by Shulyakovskyi (1969) .
The constant terms in the two equations are strikingly similar. However, the
natural convection term in Eq. 6-48 is approximately half that of Shulyakovskyi's.
With three fitted constants, this equation should give the best fit, but the
standard error is only slightly smaller than with two fitted constants.
The natural convection constant in Eq. 6-49 is very similar to that
proposed by Ryan and Harleman (1973), which was found to be satisfactory in
a numoer of evaporation experiments at the laboratory scale and the field
scale. The Monticello data essentially confirm Ryan and Harleman's natural
convection term. In fact, the results are so close that a least-squares
fit using Ryan and Harleman's (1973) natural convection constant, Eq. 3-50,
gave a smaller standard error than Eq. 6-49.
As expected, Eqs. 6-46 through 6-50 have smaller forced convection terms
than- the corresponding equations in Table 6-1. One possible cause is that
the Monticello channels experience a land surface rather than a water surface
wind boundary layer. More important yet is the observation that, except for
a narrow range of wind directions, the relatively high channel sides at
Monticello create a separated wind ooundary layer over the channel water
surface. The reduced wind velocities near the water surface in a separated
boundary layer would cause less forced convection than in an unseparated
boundary layer.
89
-------
Since any of the equations 6-46 through 6-50 will give similar accuracy
in determining a wind function for the channel studied, criteria for selec-
tion of a wind function equation will be based upon how well the equation
describes physically real processes and whether the equation may be used
for channels which do not have heat additions. Equation 6-49 fits both of
these criteria. Equations 6-46, 6-47, and 6-48 are not as good as Eqs.
6-49 or 6-50 because they contain an unexplained constant rather than a
natural convection term. In addition, for natural water"surface temperatures
where A6 is small, Eq. 6-49 will reduce to a form similar to wind functions
determined from lake evaporation in Table 6-1. Finally, the natural convec-
tion term determined by least-squares fit from the MERS field channel
data is very close to the one derived from laboratory data by Ryan and
Harleman (1973). This gives further credibility to the wind function in
Eq. 6-49 , which is selected as the most appropriate relationship.
90
-------
SECTION 7
DYNAMIC WATER TEMPERATURE PREDICTION IN OPEN CHANNELS
7.1. GENERAL OBJECTIVE
To describe water temperature dynamics in an open channel under all
types of weather conditions, a numerical solution of Eq. 1-1 is necessary.
The most important terms to be described numerically are the source terms
which incorporate heat transfer across the air-water interface and longitudinal
dispersion. The source terms are highly non-linear and require an iterative
scheme to accurately predict air-water heat transfer. Linearized source terms
may be used to solve Eq.l-r analytically for temperature distribution but the
solutions with linearized equilibrium temperature relationships are limited
to specific upstream water temperature and steady-state weather conditions. A
one-dimensional unsteady finite difference model for stream temperature
prediction titled MNSTREM for Minnesota Stream Water Temperature Prediction
Model has been developed. MNSTREM will be described herein and applied to
water temperature prediction in the MERS field channels. MNSTREM has unusually
high time resolution, meaning that it can predict rapid water temperature
changes in very shallow water.
7.2. FORMULATION OF ONE-DIMENSIONAL FINITE DIFFERENCE EQUATIONS FOR WATER
TEMPERATURE PREDICTION
The finite difference equation for water temperature with convection,
longitudinal dispersion, and air-water heat exchange will be formulated_in a
Crank-Nicholson implicit scheme. The basic transport equation is
where T = temperature ( C) ,
t = time (T) ,
x = distance (L) ,
2
A = cross-sectional area (L ) ,
Q = flow rate (L3/T) ,
D = longitudinal dispersion coefficient (L /T) ,
91
-------
B = channel surface width (L) ,
p = density of water (M/L ),
c = specific heat (E M~ °C~1) , and
P 2
S = net rate of heat transfer through air-water interface (E/L )
The equations will be formulated in the control volume , shown in Fig. 7.1
about the i»]+J5 location. Then
A - A
A ~
At .
and
BS BiSi,-j+% VSi.
pCp PCP 2 pCp
since Crank-Nicholson implicit differences are being used. Formulation of the
x-differences depends upon- the grid Peclet number.
Pe - (7.2)
L
If Pe is large/ an upstream difference scheme is appropriate (as in QUAL I) ,
and if Pe is small a central difference scheme is appropriate. For many open
channel conditions where D is significant, Pe is in an intermediate region.
ij
For the MERS experimental channels,
Q = .05 m /sec
A = 1.5 m3
D =: .1m /sec
Lt
AX = 3 •+ 8 m, and
Pe a 1 ->• 3
which is in-the range where neither the upstream nor central difference formu-
lations' are accurate (the-greatest error occurs- at-Pe=2).- The exact solution
to the steady-state transport equation without a source term gives (Chien, 1977)
. 31
3x VL ax
- T
i+l
exp(Pei+1) -
1 ^ *P T
jxp(Pei_1 ) - 1J
. (7-3)
exp(Pe.
92
-------
•Ax
j-H
i.j-M/2
] o
i-l
At
O
i-H
i= distance , node
j.— time node
x= distance
t= time
Ax= distance increment
At=time increment
Fig. 7.1. Control volume in time-distance coordinates for fornjulation
of finite differences.
-------
At large values of Pe, however, the exponential terms are expensive to compute.
A power 'law approximation to Eq. 7-3 which has been developed by Patankar and
Baliga (1978) is
«s
3x \ L 3x
Ax
Ti *
5 1
5 -,
(7-4)
In Eq. 7-4 the terms in brackets -signify the largest of the individual terms.
Because one temperature is assumed to be valid over a whole control volume,
the Peclet numbers between matrix locations use the hormonic mean of D .
L
Pe
(7-Sa)
(7-Sb)
Because the finite difference equation is to be Crank-Nickolson implicit
(rather than fully implicit) the convective-dispersive term in Eq. 7-4 should
be:
r 1
- S 22
~ 2 3x
JL
ax
3x
3x
(7-6)
Combining Eq. 7-1, 7-4, and 7-6 gives the finite difference equation:
a T . ,-bT ,+cT
1-1,3+1
+Z =0
1,3
(7-7)
where a =
Q At
2 A Ax
1 +
t1-0-1!^
94
-------
.At
b = 1 + a + c , and
At B. S. .
Z. . = (l-a-c)T. . + a T. , + c T. , . + - - — *fj *
l -
. . . . . , . , .
1,3 i,: i-l,D i+l, J pcpAi
To solve for the longitudinal temperature profile at t=j+l, Eq. 7-7 is
first rewritten as
Ti,3+l + Wi Vl, = -
were at all locations except the boundaries.
W
i b + a W.
and a G + Z
At 1=1, T is known and
1 ,
Z + a T
= -c/b and G = -3 - (7-11)
At the lower end of the channel reach to be studied (i=I), a constant gradient
of temperature will be assumed, or
T = 2 T - T
1+1, 3+1 I,j+l 1-1, 3+1
Then,
and finally
W- =-0 - - -
(a-c) S-l + ZI
I ~ b-2c+(a-c)W
= GI '
-------
The longitudinal temperature profile is then found through back substitu-
tion from i=I-l to i=2 through a rearranged form of Eg. 7-8.
T = G - W T
i, 3+1 i i i+1, j+1
Eq. 7-8 through 7-12 comprise a tri-diagonal matrix algorithm (TDMA) ,
The implicit method of solution is always stable and convergent (Mickley
et al, 1957). However, because a Crank-Nicholson implicit scheme, rather
than" fully implicit, is used oscillation matrices in the solution for tempera-
Lures can develop.- "A Crank-Nicholson-scheme- is used because it is more- •-
accurate when the transition to the final steady values of the dependent
variable are of interest. For finite difference procedures which are both
convergent and stable, a reduction of increment step size gives a good indica-
tion of result reliability (Mickley et al, 1957). To assure that no oscilla-
tion occurs in the matrix, and to test the significance of numerical damping,
the numerical solution is computed with half the original time step. If no
significant difference exists between the three runs, then no oscillation or
numerical damping problems are assumed to exist.
The source term in Eq. 7-1 is nonlinear ly dependent upon water temperature.
Therefore, the term S. cannot be included in the TDMA' scheme and is
included through iteration. First, it is assumed that S, . = S. . Then
i ,3+1 i r 3
S is computed from T. determined from the TDMA. Finally, S.
1,3+1 1,3+1 1,3+1
is used to solve new T. values. These two steps are repeated until
i , 3+-I-
there is no significant change in the T. . values. For open channels,
i , 3+-1-
such as the MERS field channels, where heat exchange is important this type
of iteration can significantly reduce errors in the prediction of water tem-
perature.
7_.3Y DESCRIPTION OF NUMERICA^FINITE" DIFFERENCE" COMPUTER" PROGRAM~FOR WATER-'" "
-"---"TEMPERATURE PREDICTION-- ../•.. ..- - . .
The water temperature prediction program MNSTREM begins with a channel
reach of specified length, distance increment, time increment, and initial
water temperature in each segment. The channel reach is divided into a
number of segments. Each segment is assumed to have a constant surface
width and cross sectional area and may consist of one or more distance increments.
96
-------
MNSTREM reads cross sectional area and surface width for each segment
as well as dimensionless longitudinal dispersion coefficient for the channel
and latitude and longitude of the study site. Upstream water temperatures,
flow rates, and the required weather parameters are read by MNSTREM to con-
tinuously compute water temperatures. Measured water temperatures may be
input to MNSTREM for comparison with computed water temperatures. Individual
residuals, the residual sum of squares, and the RMS residual are then computed.
A more detailed description of all required input parameters is given in
Appendix C (User's guide). A general flow chart for MNSTREM is given in
Fig. 7.2. The program listing, given in Appendix D, provides a more detailed
description of the required computations. A sample input is given in Appendix
E. Appendix F gives a partial sample output which includes a comparison with
measured water temperatures and computed residuals.
The most important input parameters for accurate water temperature pre-
diction with MNSTREM are:
a. Flow rate and cross sectional areas since water temperature
prediction is very sensitive to residence time,
b. Solar radiation, relative humidity, air temperature, and
wind velocity measurements for accurate air-water heat
transfer prediction, and
d. If temperature fronts or pulses are being routed D ,
At, and Ax must be chosen carefully to match real
longitudinal dispersion with the combination of D
I*
and numerical dispersion in the model.
1 \
7.4. EXAMPLES OF WATER TEMPERATURE PREDICTIONS IN THE MERS FIELD CHANNELS
\ •
The wind function (Eq. 6A49 in Table 6-1) and the dimensionless longi-
A
tudinal dispersion coefficient.-1 D* = D B/Q = 7. 47 (Section 5.4) are used
I
in MNSTREM to predict water temperatures in the MERS field channels. MNSTREM
is applied to water temperature prediction over four different time periods
and compared to recorded water temperatures.
a. January 31 through February 4, 1976, with weather data averaged over
1 hour time increments (Run A) .
97
-------
Program MNSTREM
'Read input data describing channel
reach, initial water temperatures,
flow rates, and measured temperatures
1
SUBROUTINE TEMP
CALL PRNOUT
SUBROUTINE PRNOUT
Prints intermediate statistics,
water temperature data, final
water temperature predictions,
measured water temperatures, and
residuals
Test increment size
if desired
Fig. 7.2. General flow chart for MNSTREM, the Minnesota Stream
Water Temperature Prediction Model. (Cont'd)
98
-------
SUBROUTINE TEMP
SUBROUTINE ITEMP
Sets up the initial water
temperature profile at each
i node from input data for
each segment.. _ -
Main Loop
One loop for each input
of weather data
CALL HFLUX
assuming T ,=T.
13
1
SUBROUTINE HFLUX
Computes heat flux source, S, and
longitudinal dispersion coefficient
for each channel segment
Second Loop
One loop for each hour
2.
SUBROUTINE QRATE
Computes flow rate for each hour
Tri-diagonal matrix algorithm
solution to determine"water
.~temg_epa.-tur-e£L at—time, step
Fig. 1.2. General flow chart for MNSTREM, the Minnesota Stream
Water Temperature Prediction Model. (Cont'd)
99
-------
I
Are further iterations \
on source term needed? f~
No
Yes
CALL HPLUX
Using new values
of T . ,
if D+l
SUBROUTINE HFLUX
Compute Residuals
Fig.
7.2 (Cont'd) .
General flow chart for MNSTREM, the Minnesota
Stream Water Temperature Prediction Model.
100
-------
b. April 3 through April 12, 1976, with weather data averaged over
3 hour time increments (Run B).
c. June 20 through July 19, 1977, with weather data averaged over
3 hour time increments (Run C).
d. November 16 through December 6, 1976, with weather data averaged
over 3 hour time increments (Run D).
Finally, the accuracy of using instantaneous rather than time-averaged readings
of weather data is evaluated by comparison of water temperature predictions
with the two types of data input to the model. In addition, the predictive
capabilities of time-averaged weather data over various averaging time periods
is analyzed.
In order to predict water temperatures the input data listed in Section
7.3 and Appendix C must be obtained. The water data acquisition system has
been described in Section 4.4. The MERS weather station is described in
Appendix G. Hourly inflow water temperatures and channel water temperatures
were read on water temperature strip charts and adjusted according to calibra-
tion measurements as described in Section 4.4. A weather station to measure
and record weather parameters was installed near the center of the MERS. Solar
radiation was measured with a 50 ^unction Epply pyranometer at 2 m height. Wind
velocity was measured at 9 m height with a Science Associate's anemometer.
YSI 80 recorders were calibrated and used for recording both wind velocity and
solar radiation. Air temperature and relative humidity were recorded on a
Belfort spring wound hygrothermograph, with a bimetal temperature sensor and
hair hygrothermograph. The hygrothermograph was calibrated weekly. When in-
strumentation breakdown occurred, the required weather data were obtained from
Morthern States Power Company records at the nearby Monticello Power Plant.
Fraction cloud cover was estimated from NOAA Minneapolis-St. Paul and St. Cloud
local weather data and from the Epply solar radiation records. Finally, cross
sectional area and surface width at each channel segment were determined from
field measurements as in Section 4.1.
When a convective-diffusive equation such as Eq. 6-1 is modelled through
finite differences, care must be taken to insure that the numerical dispersion
caused by the finite size of the well-mixed control volumes is small when
compared to the physical dispersion. Longitudinal dispersion is included
-------
in Eq. 6-1 as a diffusive terra. Banks (1974) studied a one-dimensional un-
steady mixed cell model with a concentration front similar to the temperature
fronts in Section 5. He found that at an infinite number of cells the model
is equivalent to Eq. 6-1 without the source term where the term Q
replaced the longitudinal dispersion coefficient (D ) . If we define an
L
approximate numerical longitudinal dispersion coefficient,
then the grid Peclet number may be expressed as
Pe = 2 "nun/0!, (7'14)
Although D is only truly equivalent to a diffusion coefficient at an
infinite number of cells, it gives a valid estimate of existing numerical
diffusion with a finite number of cells. It is interesting to note that
the greatest difference between central difference or upstream difference
formulations as described in Section 7.2 occurred at D = D . For the
num L
hybrid formulation used here (Eq. 7-4) a small grid Peclet number must be
used if longitudinal dispersion is important. The use of various distance
increments indicated that in the four water temperature prediction, Runs A
through D, longitudinal dispersion was crucial only to Run D, where two
temperature fronts were routed through the channel.
Run A: 1800 h January 31 through 2200 h February 4, 1976, with weather data
averaged over 1 hour time increments.
... This" four-dayper rod- -rnciudes~a- temporary- decrease—in-up-stream- water
temperature-and a severe cold-front. Upstream-water-temperature was mostly •
"maintained'at 15 C through the heat exchangers. Because of the artificially '
high water temperature, air-water heat transfer was large. The mean channel
flow velocity was constant at 0.016 m/sec. Air and dew point temperatures
were obtained from records at the NSP Monticello weather station. Computed
water temperatures were compared to observed hourly water temperatures at
channel locations 3, 7, 9, 13, 15, and 17. The standard error was computed
to be 0.22 C. When one considers that a decrease in water temperature of
102
-------
7°C along the channel was common, this standard error can be considered very
low. Predicted and observed water temperatures from channel locations 3, 9,
and 17 are shown in Fig. 7.3.
Run B: 0000 h April 3 through 2400 h April 12, 1976, with weather data
averaged over 3-hour time increments
For these ten days in early spring the relatively high solar radiation
and low relative humidity are in-contrast to Run A for the winter period. Air
temperature was mild, between 20 and -3 C. The mean channel flow velocity
during this period was approximately 0.013 m/sec. With the large solar radia-"
tion values, a small error in measured solar radiation or predicted water
surface reflectivity can cause significant' errors in computed water tempera-
tures at 1200, 1500, and 1800 hr. This is shown to occur periodically in
Fig. 7.4, where computed and observed water temperatures for channel locations
3, 12, and 17 are compared. While longwave, evaporative, and conductive
heat transfer were all of similar values for the winter Run A, the larger
air temperatures and lower relative humidity of Run B decreased the importance
of conductive heat transfer and increased the importance of evaporative heat
transfer. Residuals of predicted and observed water temperatures were taken
every 3 hours at channel locations 3, 6, 12, and 17, and the standard error
was determined to be 0.22 C. This standard error compares very well to the
daily variation of 7-8 C at location 17.
Run C: 1500 h June 20 through 2100 h July 19, 1977. with weather data
averaged over 3 hour increments
This run encompasses a full month in midsummer when solar radiation is
of the greatest importance. Large variations in cloud cover, wind velocity,
and solar radiation occurred over the month period. In declining importance,
most of the air-water heat transfer was due to solar radiation, evaporation,
and longwave radiation. Convective heat transfer-provided a very small.
contribution. The mean channel flow velocity throughout the period for
Run C (Period C) was approximately 0.020 m/sec. Computed and observed water
temperatures for channel locations 3 and 17 are compared in Figs. 7.5a, 7.5b,
and 7.5c. As with Run A and Run B, upstream water temperature was maintained
10-15 C above ambient so that the heat transfer is about what is expected
from cooling ponds. The temperature regime which would occur under ambient
103
-------
O
o
tu
Q.
E
O Location 3
A Location 9
Q Location 17
Jon 31
Feb I
Feb 2
Feb 3
Feb 4
Fig. 7.3. Run A. Hourly water temperatures at three stations in channel 1 from January 31
through February 4, 1976. Station locations are identified in Fig. 4.13. Symbols
are recorded values; solid lines are values simulated with MNSTREM.
-------
22
20
18
o
<° I/I
t- 0.1^
12
10
8
-] 1 p r
A Location 3
D Location 12
O Location 17
12
Fig. 7.4. Run,B. Three-hour water temperatures at three stations in Channel 1. Station
locations are identified in Fig. 4.13 (3 is upstream and 17 is downstream).
Symbols are recorded values; solid lines are values simulated with MNSTREM.
-------
o
en
A Location 3
Location 17'
Fig. 7.5a. Run C. Three-hour water temperatures at Station 3 (upstream) and Station 17 (downstream)
in Channel 8. Symbols are recorded values; solid lines are values simulated with MNSTREM.
-------
A Location 3
O Location 17
26 -
8
Fig. 7.5b. Run C. Three-hour water temperatures at Station 3 (upstream) and Station 17
(downstream) in Channel 8. Symbols are recorded values; solid lines are
values simulated with MNSTREM.
-------
o
CD
42
40
o
0)
36
34
32
30
28
A Location 3
O Locatiqn 17
II
12
13
14
15
July 1977
16
17
18
19
Pig. 7.5c. Run C. Three-hour water temperatures at Station 3 (upstream) and Station 17
(downstream) in Channel 8. Symbols are recorded values; solid lines are
values simulated with MNSTREM.
-------
conditions is shown during the period June 24 through June 26 when no addi-
tional heat was supplied to inflowing water upstream. Residuals of predicted
and observed water temperatures were taken every 3 hours at channel locations
3, 9, and 17, and the standard error was determined to be 0.29°C.
Run D: 0000 h November 16 through 2400 h December 5, 1976, with weather
data averaged over 3 hour increments
This late fall period includes four temperature fronts (2 increases and 2
decreases), a steady decrease in air temperature from approximately 5°C to -20°C,
and-consistently^large wind velocities. Upstream- water temperatures wer-e
o
maintained around 10 C above ambient. During the temperature front periods,
the computational time increment was reduced from 3600 sec to 400 sec to
prevent oscillations in the solution matrix. This is usually required when
3T/3x is very large, as in a temperature front or pulse. Grid Peclet number
was approximately 1 in the pools, so D YD = 1/2 during this run and D*
nuin Lt
was adjusted accordingly. The mean channel flow velocity during this period
was approximately 0.015 m/sec. As with Run A, air temperature and dew point
temperature were obtained from weather station records for the NSP Monticello
Power Plant. The predicted and observed water temperature, from channel
locations 3 and 17 are compared in Fig. 7.6a and 7.6b. During much of the
period predicted downstream temperatures are larger than observed. This
could be because the wind direction was predominantly parallel to the down-
stream channel sections, and wind direction is not included in the wind
function determined in section 6.5. In tne later half of Period D, low air
temperatures and a variable wind velocity created a highly dynamic temperature
regime with large longitudinal gradients. As with Run A, evaporative, con-
vective, longwave radiation and solar radiation heat transfer were all of
-equal—importance. -Residuals-of- predicted.and, observed..water, .temperatures . _
were taken-every--3 hours, at. ctiamiel,J-ocationaJJ^..6J.,12,_._and_.17_, and the "standard
error was determined to be-0-.32 CT --- -' - '"._.. .. .._.'_
The water temperature predictions were chosen to cover all four seasons
of the year. Runs A, B, C, and D indicate that MNSTREM accurately predicted
water temperatures which were both ambient and artificially heated for each
season of the year. Upstream water temperature varied from 0 to 40°C and
air temperature varied from -25 C to +35 C. With this wide range of conditions,
109
-------
I ' 1 '
A Location 3
O Location 17
20 21
November 1976
24
25
Fig. 7.6a. Run D. Three-hour water temperatures at Station 3 (upstream) and Station 17 (down-
stream) in Channel 1. Symbols are recorded values; solid lines are values simulated
with MNSTREM.
-------
A Location 3
O Location 17
27 28
November 1976
3 4
December 1976
Fig. 7.6b. Run D. Three-hour water temperatures at Station 3 (upstream) and Station 17 (down-
stream) in Channel 1. Symbols "are recorded values; solid lines are values simulated
with MNSTREM.
-------
the maximum standard error for any run was 0.316 C, which compares to a calibra-
tion error of 0.10 to 0.2 C for the thermistors which measured water temperature.
The low standard errors are a good indication of MNSTREM's capability for pre-
dicting water temperature, especially when one considers that MNSTREM has no
coefficients determined from the data of any of the predictive Runs A, B, C,
or D.
When predicting water temperatures over a stream reach, an investigator
must decide- what~time peraod~ts* adequate for averaging measured weather para-
meters. Since hourly averages of weather data were used for Run.A, this period
was used to compare the accuracy of water temperature predictions with weather
data averaged over 1, 3, 6, and 12 hour periods. The results, given in Table
7.1, indicate that for the period in Run A, averaging weather parameters over
1, 3, and 6 hour periods gave equally good predictions.
TABLE 7.1. STANDARD ERROR OF WATER TEMPERATURE
PREDICTIONS WHEN WEATHER PARAMETERS
ARE AVERAGED OVER 1, 3, 6, AND 12 HR PERIODS
Time
Averaging
Period for
Weather Parameters
(hr)
1
3
6
12
Standard Error of
Prediction
C°0
0.222
0.243
0.239
0.407
It should be noted that Run A was for a winter period when daily variations
of solar radiation are small. " Other weather parameters such" as wind"velocity
and air temperature, however, _were highly variable in Run A.
Instantaneous values of weather data are more easily obtained than time
averaged values. Much of the NOAA Local Climatological Data, for example,
is given as instantaneous values every three hours. To test the accuracy of
water temperature prediction with instantaneous rather than time-averaged
readings of weather data, water temperatures were predicted with 3, 6, and 12
hour instantaneous weather data for Period A. In addition, the effect of
112
-------
lead time in the input of the instantaneous data to the computer model was
studied. The results are given in Table 7.2. For six-hour instantaneous
weather data, a lead time of 1.5 hours means the weather data is applied 4.5
hours before the reading and 1.5 hours after the reading time. The results
of Table 7.2 indicate that, for example
1) a three-hour instantaneous reading of weather data taken at
1200 would give the best prediction if applied from 1030 to
1330, based on the lowest standard error of prediction
2) a six-hour instantaneous reading at 1200 should be applied
from 0830 to 1430, and
3) a 12-hour instantaneous reading at 0100 should be applied
from 0230 to 1430.
The results given in Table 7.2 are only strictly applicable to the MERS
channels for Period A. They give an indication of what can be expected in
a shallow channel.
113
-------
TABLE 7.2. STANDARD ERROR OF WATER TEMPERATURE PREDICTIONS WITH
INSTANTANEOUS READINGS OF WEATHER DATA TAKEN AT 3, 6,
AND 12 HOUR INCREMENTS*
Lead Time
in Input of Re
Weather Data e
(hr)
0.5
1.5
2.5
3.5
4.5
5.5
6.5
7.5
8.5
9.5
10.5
11.5
,. Standard Error
very: 3 (hrs) 6 (hrs)
<°0 (°C)
0.265 0.379
0.248 0.271
0.278 0.261
0.308
0.420
0.364
*
12 (hrs)
<°C)
0.678
0.474
0.297
0.344
0.534
0.402
0.539
0.649
0.661
0.705
0.793
0.845
*A lead time of 2.5 hours means the instantaneous data is applied for 2.5 hours
after the reading time, and N-2.5 nours before the reading time, where
N = 3, 6, or 12.
114
-------
REFERENCES
Anderson, E. R.' (1954) "Energy-Budget Studies," Water-Loss Investigations,
Lake Hefner Studies, U. S. Geological Survey Professional Paper No. 269.
Banks, R. B. (1974) "A Mixing Cell Model for Longitudinal Dispersion in Open
Channels," Water Resources Research, Vol. 10,No. 2.
Bansal, M. K. (1971) "Dispersion in Natural Streams," Proc. ASCE, Vol. 97,
No. HY11.
Barnes, Jr., H. H. (1967) "Roughness Characteristics of Natural Channels,"
U. S. Geological Survey Water Supply Paper No. 1849.
Bowen, I. S. (1926) "The Ratio of Heat Losses by Conduction and by Evapora-
tion from any Water Surface," Physical Review, Vol. 27.
Brady, D. K., W. L. Graves, and J. C. Geyer (1969) "Surface Heat Exchange at
Power Plant Cooling Lakes," The Johns Hopkins University, Dept. of
Geography and Environmental Engineering, Report No. 49-5.
Chien, J. C. (1977) "A General Finite-Difference Formulation with Application
to Navier-Stokes Equations," Computers and Fluids, Vol. 5(1):15.
Day, T. J. (1975) "Longitudinal Dispersion in Natural Channels," Water
Resources Research, Vol. 11, No. 6, p. 909.
Dingman, S. L. and A. Assur (1967) "The Effect of Thermal Pollution on River-
Ice Conditions—Part II, U. S. Army Cold Regions Research and Engineering
Laboratory, Hanover, New Hampshire.
"D~£tm~arsT~JY~ir.~~{r9T47—"Mixing- and-Transport— (lifcer-a-ture -review-)-r-Journal. . _
of the Water Pollution Control Federation, Vol. 46, No. 6, pp. 1591-1604.
Sdinger, J. E. and""J. C. Geyer (1965) Heat "Exchange in the_Environment,•- —.
The Johns Hopkins University, R.P. 49-2.
Edinger, J. E., D. W. Duttweiler, and J. C. Geyer (1968) "The Response of
Water Temperatures to Meteorological Conditions," Water Resources Research,
Vol. 4, No. 5.
Elder, J. W. (1959) "The Dispersion of Marked Fluid in Turbulent Shear Flow,"
Journal of Fluid Mechanics, Vol. 5, pp. 554-560.
115
-------
El-Hadi, N.D.A. and K. S. Davar (1976) "Longitudinal Dispersion for Flow
Over Rough Beds," Journal of Hydraulic Div., ASCE, Vol. 102, HY4, pp. 483-
493.
Ficke, J. F. (1972) "Comparison of Evaporation Computation Methods, Pretty
Lakes, Lagrange City,Northwestern Indiana," U. S. Geological Survey
Professional Paper No. 65.
Fischer, H. B. (1967) "The Mechanics of Dispersion in Natural Streams,"
Proc.. ASCE. Vol. 93, HY6.
Fischer, H. B. (1968) "Dispersion Prediction in Natural Streams," Proc.
ASCE, Vol. 94, SA5.
Fischer, H. B. (1973) "Longitudinal Dispersion and Turbulent Mixing in Open-
Channel Flow," Annual Reviews of Fluid Mechanics, Vol. 5, pp. 59-78.
Fischer, H. B. et al (1979) Chapter 5 in Mixing in Inland and Coastal Waters,
Academic Press.
Fu, A. Y. and H. Stefan (1979a) "User's Manual for Operational Water Tempera-
ture Statistics Computer Programs WTEMP1 AND WTEMP2," St. Anthony Falls
Hydrulic Laboratory, University of Minnesota, External Memorandum No. 162.
Fu, A. Y. and H. Stefan (1979b) "Water Temperature Data Processing for the
Experimental Field Channels at the USEPA Ecological Research Station in
Monticello, Minnesota," St. Anthony Falls Hydraulic Laboratory, University
of Minnesota, External Memorandum No. 167, 32 pp.
Godfrey, R. G. and B. J. Frederick (1970) "Stream Dispersion at Selected
Sites," U. S. Geological Survey Professional Paper 433-K.
Gulliver, J. S. (1977) "Analysis of Surface Heat Exchange and Longitudinal
Diapersion in a Narrow Open Field Channel with Application to Water
Temperature Prediction," M.S. Thesis, University of Minnesota, Department
of Civil and Mineral Engineering, 187 pp.
Gulliver, J. S. and H. Stefan (1980) "Soil Thermal Conductivity and Tempera-
ture Prediction in the Bed of the Experimental Field Channels at the
USEPA Ecological Research Station in Monticello, Minnesota," St. Anthony
Falls Hydraulic Laboratory, University of Minnesota^ External Memorandum
No. 165. ... . "
Hahn, M. G., J. S. Gulliver", and H. Stefan (1978a) "Operational Water Tempera-
ture Characteristics in Channel No. 1 of the USEPA Monticello Ecological
Research Station," St. Anthony Falls Hydraulic Laboratory, University
of Minnesota, External Memorandum No. 151, 46 pp.
Hahn, M. G., J. S. Gulliver, and H. Stefan (1978b) "Physical Characteristics
of the Experimental Field Channels at the USEPA Ecological Research Station
in Monticello, Minnesota," St. Anthony Falls Hydraulic Laboratory,
University of Minnesota, External Memorandum No. 156, 36 pp.
116
-------
Hahn, M. G. (1978) "Experimental Studies of Vertical Mixing in Temperature
Stratified Waters," M.S. Thesis, University of Minnesota, Minneapolis,
Minnesota, 179 pp.
Harbeck, G. E. (1958) "Water Loss Investigations, Lake Mead Studies," U.S.
Geological Survey Professional Paper No. 298.
Hughes, G. H. (1964) "Analysis of Techniques Used to Measure Evaporation
from Salton Sea, California," U. S. Geological Survey, Professional Paper
No. 272-H.
Idso, S. B. and R. D. Jackson (1969) "Thermal Radiation from the Atmosphere,"
Journal of Geophysical Research, Vol. 74, No. 23.
Jain, S. C. (1976) "Longitudinal Dispersion Coefficients for Streams," Journal
of the Environmental Engineering Division, ASCE, Vol. 102, No. EE2, April
1976, pp. 465-474.
Jobson, H. E. and N. Yotsukura (1972) "Mechanics of Heat Transfer in Non-
Stratified Open Channel Flows," Institute of River Mechanics Paper,
Colorado State University.
Koberg, G. E. (1964) "Methods to Compute Long-wave Radiation from the
Atmosphere and Reflected Solar Radiation from a Water Surface," U. S.
Geological Survey, Professional Paper 272-F.
Kohler, M. A. (1954) "Lake and Pan Evaporation," Lake Hefner Studies, U. S.
Geological Survey Professional Paper No. 269.
Limerinos, J. T. (1970) "Determination of the Manning Coefficient from
Measured Bed Roughness in Natural Channels," U. S. Geological Survey
Water Supply Paper 1898-B.
Liu, H. and A.H.D. Cheng (1980) "Modified Fickian Model for Predicting
Dispersion," Journal of the Hydraulics Div.,ASCE, Vol. 106, No. HY6,
June 1980, pp. 1021-1040.
Lomax, C. C. and J. F. Orsborn (1971) "Flushing of Small Streams," USEPA
Water Pollution Control Research Series, 16010, DMG.
Marciano, J. J. and G. E. Harbeck (1954) "Mass Transfer Studies," Lake
Hefner Studies, U. S. Geological Survey Professional Paper No. 267.
McQuivey, P. S. and T. N. Keefer (1974) "Simple Method for" Predicting-
Dispersion in Streams," Journal 3nv. Sngrg. Div., ASC2, Vol. 100, EE4,
pp. 997-1011.
Mickley, H. S., T. K. Sherwood, and C, E. Reed (1957) Applied Mathematics
in Chemical Engineering, McGraw-Hill, New York.
117
-------
Murray, F. W. (1967) "On the Computation of Saturation Vapor Pressure,"
Journal of Applied Meteorology, Vol. 6, No. 1.
Olsen, R. M. (1973) Essentials of Engineering Fluid Mechanics!, Intext
Educational Publishers, New York.
Paily, Er. P., E. 0. Macagno, and J. F. Kennedy (1974a) "Winter-Regime Surface
Heat Loss from Heated Streams," Institute of Hydraulic Research, Report
No. 155, University of Iowa.
Paily, P.., P..,. E. 0. Macagno, and J. F.Kennedy (1974b) "Winter Regime Thermal
Response of Heated Streams," Proc. ASCE, Vol."100, HY4.
Patankar, S. V. and B. R. Baliga (1978) "A New Finite Difference Scheme
for Parabolic Differential Equations," Numerical Heat Transfer, Vol. 1:
27-37.
Ryan, P. J. and D.R.F. Harleman (1973) "An Analytical and Experimental Study
of Transient Cooling Pond Behavior," Parsons Laboratory for Water Resources
and Hydrodynamics Report No. 161, Massachusetts Institute of Technology
Rymsha, V. A. and R. V. Doschenko (1958) "The Investigations of Heat Loss
from Water Surfaces in Winter Time," (in Russian), Turdy GG1, No. 65.
Sayre, W. W. and F.. M. Chaag (1968) "A Laboratory Investigation of the
Open Channel Dispersion Process for Dissolved, Suspended and Floating
Dispersants," U. S. Geological Survey Professional Paper 433-E.
Shulyakovskyi, L. G. (1969) "Formula for Computing Evaporation with Allowance
for the Temperature of Free Water Surface," Soviet Hydrology - Selected
Papers, No. 6.
Stefan, H. and J. S. Gulliver (1980) "Pore Water Temperatures and Heat
Transfer in a Riffle (Rock) Section of the Experimental Field Channels
at the USEPA Ecological Research Station at Monticello, Minnesota,"
St. Anthony Falls Hydraulic Laboratory, University of Minnesota, External
Memorandum No. 166.
Taylor G. I. (1954) "The Dispersion of Matter in Turbulent Flow Through a
Pipe," Proc. Royal Society of London, No. 223-A, pp. 446-468.
., Tennessee Valley Authority (1968) "The Water Temperature Regime of Fully
Mixed Streams,""Water Resources"Research Laboratory, Report No. 15.. ~ " •
Thackston, E. L. and P. A.Krenkel (1967) "Longitudinal Mixing in Natural
Streams," Proc. ASCE, Vol. 93, SA5.
•
Turner, J. F. (1966) "Evaporation Study in a Humid Region—Lake Michie,
North Carolina," U. S. Geological Survey Professional Paper No. 272-G.
Yotsukura, N., A. P. Jackman, and C. R. Faust (1973) "Approximations of Heat
Exchange at the Air-Water Interface," Water Resources Research, Vol. 9,
No. 1.
118
-------
APPENDIX A
TEiMPERATURE FRONT DATA FOR DETERMINATION
OF LONGITUDINAL DISPERSION COEFFICIENT
119
-------
TABLE A-l. CHANNEL 1 CROSS-SECTIONAL AREAS AND SURFACE WIDTHS
FOR TEMPERATURE FRONT EXPERIMENTS
Section Nr
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
Cross-sectional
Area (ra2)
1.17
0.46
1.44
0.48
1.62
0.37
1.64
0.37
1.77
0.37
1.76
0.39
1.83
0.38
1.37
0.38
Surface Width (m)
3.05
2.01
3.35
1.89
3.35
2.04
3.35
1.95
3.35
1.80
3.35
1.80
3.26
1.80
3.20
1.65
TABLE A-2. RESIDENCE TIMES TO EACH STATION IN TEMPERATURE
FRONT EXPERIMENTS (HRS)
Date of
Experiment
11/17/76
11/18/76
11/22/76
12/06/76
12/14/76'
3
0.31
0.37
0.31
0.15
0.29
5
0.76
0.82
0.74
0.56
0.67
7
1.25
1.28
1.20
1.06
1.08
Stations
9
1.72
1.68
1.65
1.50
1.47
1
2.22
2.14
2.06
1.98
1.88
13
2.83
2.58
2.64
2.50
"2.30
15
3.64
3.04
3.13
2.99
2-. 86
17
4.15
3.41
3.58
3.46
- —
120
-------
TABLE A-3. MEASURED 0 VALUES FOR 11/17/76 TEMPERATURE FRONT (Cont'd)
Time Since
Front Inception
(ht)
.25
.31
.38
.44
.50
.56
.63
.69
.75
.81
.88
.94
1.00
1.13
1.25
1.38
1.50
1.63
1.75
1.88
2.00
2.13
2.25
2.38
2.50
2.63
2.75
2.88
3.00
3.13
3.25
3.38
3.50
3.63__
3.75
4.00 - -
4-.25~- _..
4.50
4.75
5.00
Stations
3
.83
.49
.33
.23
.15
.10
.06
.03
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—_
-
--
.
—
—
5
1.00
1.00
1.00
1.00
.96 -
.94
.80
.65
.52
.39
.31
.23
.15
.10
.05
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
7
1.00
II
It
n
1.00
1.00
1.00
1.00
1.00
—
—
.99
.94
.76
.49
.32
.21
.14
.09
.05
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
9
1.00
II
n
n
n
ii
it
it
n
1.00
1.00
1.00
1.00
—
—
.97
.83
.63
.46
.32
.21
.15
.11
.06
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
11
1.00
n
it
n
it
n
n
n
it
n
n
n
n
1.00
1.00
1.00
1.00
1.00
—
.90
.78
.60
.47
.37
.29
.22
.18
.14
.11
.08
—
—
—
—
—
—
—
—
—
13
1.00
n
n
it
n
n
II
11
U
n
II
n
n
n
ii
ii
it
n
1.00
1.00
1.00
1.00
.97
.89
.76
.67
.55
.47
.40
.35
.31
.27
.24
—
—
—
—
—
—
15
1.00
It
It
II
II
H
II
It
It
II
II
n
n
u
n
u
n
ti
11
ii
11
u
1.00
1.00
1.00
—
.96
.89
.80
.73
.66
.60
.55
.51
.47
.40
.36
.33
.30
.23
Note: 0 is defined by Eg. 5.5.
0 = [T(t) - T(0)]/[T(oo) - T(0)]
121
-------
TABLE A-3. MEASURED 0 VALUES FOR 11/18/76 TEMPERATURE FRONT (Cont'd)
Time Since
Front Inception
(hr)
.25-_
.31
.38
.44
.50
.56
.63
.69
.75
.81
.88
.94
1.00
1.13
1.25
1.38
1.50
1.63
1.75
1.88
2.00
2.13
2.25
2.38
2.50
2.63
2.75
2.38
3.00
3.13
""• "3V2S •"
3-jJ5-
" .3.50 .. . .
"" 3763 " """" " ~
3.75
3
.89
.66
.51
.38
.28
.20
.15
.10
.08
.07
.05
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
— "
—
5
1.00
"
1.00
—
—
.95
.86
.73
.61
.50
.40
.31
.23
.13
.09
.07
. —
—
—
—
—
—
—
—
—
—
—
—
—
—
--
" —
—
7
1.00
"
n
It
it
it
n
it
1.00
1.00
—
—
.91
.72
.53
.45
.22
.13
.08
.07
—
—
—
—
—
—
—
—
—
—
—
~-
- ---
—
Stations
9
1.00
n
ti
- ,«*-
fi
n
n
n
n
n
1.00
1.00
—
—
—
.90
.75
.57
.42
.31
.21
.13
.08
—
—
—
—
—
—
—
- — — — -
—
—
-" -' --
—
11
1.00
ri*
it
~ .» -
n
n
n
n
it
it
n
it
1.00
1.00
1.00
—
—
.95
.88
.76
.64
.52
.37
.27
.17
.10
.05
—
—
—
—
—
—
--
—
13
1.00
, . n
n
*"
n
n
n
n
it
n
n
n
H
ft
It
1.00
1.00
—
—
—
.93
.87
.77
.68
.57
.45
.32
.21
.15
.10
.06
—
—
—
—
15
1.00^
II
II
II
t|
11
"
11
II
II
II
n
n
n
it
n
ii
1.00
1.00
1.00
—
—
—
.92
.87
.79
.68
.62
.52
.42
- -- .34-
.-22"
.16
.10 "
.07
122
-------
TABLE A-3. MEASURED 8 VALUES FOR 11/22/76 TEMPERATURE FRONT (Cont'd)
Time Since
Front Inception
(hr)
.19
.25
.31
.38
.44
.50
.56
.63
.69
.75
.81
.88
.94
1.00
1.13
1.25
1.38
1.50
1.63
1.75
1.88
2.00
2.13
2.25
2.38
2.50
2.63
2.75
2.88
3.00
3.13
.3.25
3.38
3.50
3.63
3.75
3.88
" 4.00
Stations
3
.99
.76
.46
.30
.18
.12
.07
.05
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
5
1.00
11
11
it
1.00
.96
.86
.75
.57
.47
.37
.28
.20
.14
.09
.05
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
7
1.00
tl
II
II
If
II
tl
It
1.00
1.00
.99
.96
.92
.84
.59
.43
.28
.20
.13
.08
.05
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
9
1.00
rt
n
it
n
n
it
n
n
n
tt
n
»
1.00
1.00
.96
.85
.69
.51
.39
.29
.22
.18
.14
.10
.08
.07
.05
—
—
—
— - -
—
—
—
—
—
--
11
1.00
n
n
ii
it
n
tl
it
n
it
n
n
tt
n
n
n
1.00
.98
.96
—
.75
.62
.46
.40
.33
.26
.20
.17
.14
.12
—
.09
—
—
.06
—
—
—
13
1.00
n
ti
n
n
ii
n
n
it
n
ft
tt
11
ii
n
n
it
it
it
ii
1.00
.98
.91
.84
.72
.60
.48
.42
.36
.30
.25
.21
.17
.15
—
.12
—
.09
15
1.00
n
n
n
n
n
ti
n
n
n
it
it
n
n
tt
ii
ti
ti
ti
. H
n
n
1.00
1.00
.97
.92
.85
.77
.68
.59
.50
.42"
.36
.31
.28
.23
.20
.17
123
-------
TABLE A-3. MEASURED 0 VALUES FOR 12/06/76 TEMPERATURE FRONT (Cont'd)
Time Since
Front Inception
(hr)
.13
.19
.25
.31
.38
.44
.50
.56
.63
.69
.75
.81
.88
.94
1.00
1.13
1.25
1.38
1.50
1.63
1.75
1.88
2.00
2.13
2.25
2.38
2.50
2.63
2.75
2.38
3.00
3.13
3.25
3.38
3.50
3.75
3.88
4.00
Stations
3
.64
.37
.24
.15
.10
.07
—
—
—
—
—
—
—
—
—
—
—
— •
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
5
1.00
1.00
1.00
—
.95
.84
.69
.50
.36
.27
.21
.15
.12
.08
.07
—
—
— '
—
—
—
--
—
—
—
—
—
—
—
—
—
. —
. —
—
—
—
—
—
7_
1.00
if
n
n
ti
"
ti
n
1.00
-~
.97
.93
.87
.77
.63
.38
.26
.18
.12
.08
.05
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
9
1.00
11
11
n
ti
n
n
n
n
ti
1.00
ii
11
1.00
.95
.86
.70 •
.50
.38
.27
.20
.16
.12
.09
.07
--
—
—
—
—
—
—
—
—
—
—
—
11
1.00
II
n
n
ti
n
ti
n
n
n
tr
11
it
n
n
it
it
1.00
—
.91
.77
.62
.47
.35
.27
.21
.14
.11
—
.08
—
—
—
—
—
—
—
--
13
1.00
ft
n
n
n
n
n
"
ii
n
n
it
n
11
ti
ti
ti
11
it
11
1.00
—
.95
.85
.73
.61
.50
.41
.35
.29
.24
.21
.18
.15
—
.11
—
.09
15
1.00
n
n
n
n
n
n
n
n
H
ii
n
n
n
n
n
n
n
n
n
n
ii
n
1.00
1.00
.97
.90
.79
.69
.60
.48
.42
.35
.29
.21
.18
.15
.13
124
-------
TABLE A-3(Cont'd). MEASURED 9 VALUES FOR 12/14/76 TEMPERATURE FRONT
Time Since
Front Inception 3
(hr)
.25 .72
.31 .42
- -.--..38 ~- -. .26
-- -V44 ~ — ~ - .15'
.50 -.09
- ~ .56 ' * 'rO^S
.63
.69
.75
.81
.88
.94
1.00
1.13
1.25
1.38
1.50
1.63
1.75
1.88
2.00
2.13
2.25
2.38
2.50
2.63
2.75
2.88
3.00
3.13
3.25
3.38
" 3-.50 " -- "
3.63. - .. ". ._—
3.75-
4.00 " —
5
1.00
n
11
Q-G—
- .ati
.61
.45
.34
.25
.21
.16
.13
.08
.05
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
— —
— __ _
—
~ — "
7
1.00
11
"
n
«
n
—
—
.95
.88
.79
.65
.43
.29
.19
.10
.07
—
—
—
—
—
—
—
—
—
—
—
—
—
—
— —
— .
. —
—
Stations
9
1.00
n
n
- ii-
n
n
n
it
it
n
n
1.00
.95
.82
.64
.47
.35
' .24
.16
.11
.06
—
—
—
—
—
—
—
—
—
—
—
...
—
—
11
1.00
II
It
II-
"ll
II
II
n
it
it
it
ii
1.00
—
.95
.90
.79
.66
.51
.39
.29
.21
.15
.11
.10
—
—
—
—
—
—
— —
—
—
13
1.00
n
n
- it
11
n
it
it
n
n
n
n
n
ft
it
11
1.00
—
.95
.88
.81
.69
.55
.45
.36
.29
.24
.19
.16
.12
.10
.09
— —
...07._.
—
~~—
15
1.00
If
tl
II
It
It
It
It
It
II
tt
n
it
ti
ti
it
it
n
1.00
—
.97
.90
.84
.76
.66
.56
.48
.44
.37
.32
.29
. 26
— _
.19
.13
125
-------
APPENDIX B
DATA AND COMPUTATIONS OF BULK SURFACE HEAT
TRANSFER COEFFICIENT AND WIND FUNCTION PARAMETERS
FROM STEADY-STATE LONGITUDINAL WATER TEMPERATURE PROFILES
126
-------
Sample Calculation of Bulk Surface Heat Transfer Coefficient
and Wind Function Parameters from Steady State
Longitudinal Temperature Profile
A near steady-state longitudinal temperature profile was observed on
strip chart records of water temperature on November 29, 1976. Weather
o
records indicated that T = -14 C, W = 4.5 + 1 m/sec, and relative humidity
a y —
(RH)was 91 per cent. Using psychometric tables,- dew point was determined to
be -14 C for the given T and RH.--This and all-other-steady-state--per igdsT
Cl "*
occurred during non-daylight hours, so T = T was a good approximation.
£ U
Water temperatures were adjusted by the corresponding temperature cali-
bration and values of 6/6 = [T(x)-T ]/(T -T ) were computed and plotted on
semi-log paper as shown in Fig. B-l. A best fit line through all points
gave 6/9 (x=487.7 m)*= .796. All points, however, are enclosed by maximum
and minimum 9/9 (x=487.7 m) values of .865 and .760, respectively.
Sample Computation:
For 8/eQ = .796, x = 1600 ft = 487.7 m, Q = 473 gpra = 29.8 2,/S
and b = 9.5 ft = 2.9 m, K from Eq. 6-34 = 41.66 cal cm"2 day ~l °C~1
5
T =8.0 (mean channel water temperature) C
w
Then
3 = from Eq. 6-43 = .4156 mb/ °C .
and
F from Eq. 6-42 = 31.59 cal cm" day mb
w
p = air pressure = 1013 mb '(assumed)
3 ^
_.. e = saturation vapor pressure at T = 8.05-mb -
e = saturation vapor pressure at T = 1.56 mb .
clZ Q
Then
A6v = from Eq. 6-30 = 22.95°C
and
FW - 5.88(A6y) / from Eq. 3-45 = 15.91 cal cm"1 day"1 mb"1
* 1600 ft.
129
-------
U)
O
fi/
fa
1.0
0.9
0.8
0.7
0.6
0.5
0.4,
O Observed water temperatures
—— "Best fit" steady state profile
I
I
I
100
200 300 400
Distance Downstream (m)
500
Fig. B-l. Steady state water temperature profile and "best fit" curve for
November 29, 1976.
-------
The maximum AT/At of all channel locations was 2.0 C/day, and the
temperature difference between inflow and outflow temperatures was AT = 5.0°C.
Thus the steady state parameter was
N§ _ AT/day
=
SS 30.39 Q(gpm)(AT)
where Q = 473 gpm. The maximum per cent error in T(x) from non-steady-
3 = T _ T
o o E
state temperatures is 100(1 - 1/1.05) = 5%. Since 9 = T - T_ = 15° + 14°
Q O £»
29°C, from Fig. B-l
T{x = 1600 ft) = .796 9 + T_ = 9.08°C
O Ct
and the maximum per cent error Z in 9/0 (x = 1600 ft)
0/9 -[.95T)x = 1600) - T ]/8
Z = 100 —2 - _ - E_o . 2%
Then the maximum error in K calculations is
s
In 9/9 - In (1.0 2 9/9 )
%(NS/SS) error in K = + 100 - ." , ta - = + 8.7%
s — xn y/o —
At a flow of 473 gpm, the accuracy of flow rate measurement in the MFS
experimental channels is *• 4 per cent. Equation 6-34 indicates that this
will transfer directly to a +4 per cent error in K computations. Thus,
™ 5
there are three significant contributions to errors in K .
Equations 6-42 and 6-47 indicate that the magnitude of the maximum error
in K , rather than the per cent error, is carried to both F and
s w
_ I-/*! - »- - _ _-,_.._ - __ _ ___ - -_ — ___
F - 5.52(A9 ) -- . The maximum and minimum values of K , F , and- -F-- - 5.52
wv _ sw.w
(A6T ) which would occur due to each .component of error and with a combi-
nation of all three errors are given in Table B-2. With such a large
potential error, the scatter in Figs. 6-1 through 6-4 is understandable.
131
-------
TABLE B-2^ MAXIMUM AND MINIMUM VALUES OF K , F , and F -5.52(A9 ) 1f IN SAMPLE
S W W V
i' COMPUTATION DUE TO MAXIMUM POSSIBLE ERRORS IN DATA
A i
Best estimate
Ranges of estimates:
K (cal cm"2 day"1 °C~1) F (cal cm"2mb~T)
41.66
31.59
Fw-5.52(A8v)
1/3
—2 —1 —1
(cal cm day mb )
15.91
Including estimate of 6/0
o
Maximum
Minimum
Including errors
steady character
temperature ' '
Maximum 3
Minimum « .
Including errbr^
measurement . i
' l
Maximum »'
Minimum --
' 1
i
)
di|e to un-
of water
i
•>
t
in flow
i
1
i
50
26
45
38
43
40
.11 39.83 !
.48 16.79
.28 35.12
.04 28.07
.33 32.22
.0 29.98
24.15
1 1.11
i
19.44
12.38
1
17.54
14.29
Including corabinatipn of
all sources of error
Maximum '
t
Minimum /
1
56
23
.65 46.21
.21 13.61
30.52
- 2.08
-------
APPENDIX C
USER'S GUIDE TO THE MINNESOTA
STREAM WATER TEMPERATURE PREDICTION MODEL (MNSTREM)
As noted in Section 7.3, the most important parameters for accurate
water temperature prediction with MNSTREM are flow rate, cross-sectional area,
solar radiation, air temperature, relative humidity, and sometimes longitudinal
dispersion coefficient. There are also a number of other parameters necessary
in running the model. This appendix gives a line by line description of the
input file with format specifications given in parenthesis. Sample data from
Appendix E are given in brackets.
133
-------
Line 1 (8P 10.2)
DSTAR, TSTART, DELX, DTIME, DAY1, TLNGTH , LAT, ALT
DSTAR = Dimensionless longitudinal dispersion coefficient = D. B/Q [7.47]
L
TSTART = Starting hour of day (hours.minutes) [00.00]
DfiLX = Distance increment (ft) for computations [12.50]
DTIME = Time increment (sec) for computations [3600]
DAY1 = Day of month for starting day [16]
TLNGTH = Total length of channel reach (ft) [1700]
LAT = Latitude of study site (degrees) [45]
ALT = Altitude of study site (ft) [1000]
Line 2 (8F 10.2)
DATE, CHDATE
DATE = Month.Day of starting day [11.16]
CHDATE = Number of days in month/100 [.30]
Line 3 (915)
NHOUR, MHOUR, NRES, NP, NSECT, IFLAG, JFLAG, KFLAG, LFLAG
NHOUR = Total hours of time period for computer run [480]
MHOUR = Number of hours between weather data input [3]
NRES = Number of channel locations with observed water temperatures [3]
NF = Number of flow rates read in [11]
NSECT = Number of constant cross-sectional segments in channel reach [17]
IFLAG = 1 if observed water temperatures are read in
JFLAG = 0 if the longitudinal dispersion model is used
JFLAG = 1 if computations are by upstream differences (numerical
dispersion only)
KFLAG =1 if D , rather than DSTAR, is specified at each segment
Jj
LFLAG = 1 if test of increment size is called
Set 4 (9F5.1) Two lines
T0(l) (TI(I), 1=1, 17)
T0(l) = inflow temperature for first time step (°C) [11.3]
TI(I) = initial temperature at channel segment I [TI(1)=10.7,TI(16)=8.1]
'Line 5 (10A8)
HEAD (I), 1=1, 10 . . . .
HEAD = Heading for output of weather data which includes date and
time of run
134
-------
Set 6 (7F 10.4). NHOUR/MHOUR lines
TA(J), RH(J), RAD(J), W(J), WDIR(J), P(J), CC(J), TIME
Set 6 is repeated for J=l to J=NHOUR/MHOUR
TA(J) = Air temperature (°C) [TA(1)=-10]
RH(J) = Relative humidity (%) [RH(1)=100]
RAD(J) = Total solar radiation (cal cm-2 min~l) [RAD(3)=.09]
W(J) = Wind velocity (m/sec) [w(4)=3]
WDIR(J) =s Wind direction (degrees) [Blank column in Appendix E]
P(J) = Air Pressure (Bars) [PA(1)=979]
CC(J) = Fraction cloud cover [CC"fl3)=0.5]
TIME = Time of day (hours), [03 in first row-of Set 6] -
The last column gives data>-which is "not r.ead-into-
If IFLAG = 0, the next set of lines is skipped.
Set 7(18F 6.2). NHOUR/MHOUR lines
TO(J), (TT(IT+1, IT+NSECT)
Set 7 is repeated for
J=l to J=NHOUR/MHOUR, where IT=>(J-1) *NSECT
TT(I) = Observed water temperatures. Use 0.0 if no observed water
temperature is input at segment. [Not in Appendix E]
Temperature measurements taken at a given time would occupy one row.
Set 8(12F 5.1). NHOUR/12+1 lines
TO(J), J-l, NHOUR+1 [TO(1)=11.1, TO(13)=10.8]
TO(J) = Upstream water temperature which is read for each hour
Date and hour given in the last two columns is not read into MNSTREM
Set 9 (10F 7.0). NF/10+1 lines
Q(I), 1=1, NF
Q(I) = Flow rate (gpm) [Q{l)-473, Q(11)=456]
NF = Number of input flow rates
Set 10 (1015). NF/10+1 lines
---" IF(I)y-1=1, NF - -- - - ----_-_ -- -- -.-.-.
"IF(I) =~TJumber of hours Q"(I) is "applied [lF(l)=8i 'IF(II)=192]
Set 11(10F 7.2). MSECT/10+1 lines
A(I), 1=1, NSECT
A(I) = Cross-sectional area of segment I(ft ) [A(l)=14.2, A(17)=20]
Set 12 (10F 7.2). NSECT/10+1 lines
B(I), 1=1, NSECT
2
B(I) = surface width of channel segment I(ft ) [B(l)=10.8, B(17)=10.5]
135
-------
Set 13 (10F 7.2). NSECT/10+1 lines
DL(I), 1=1, NSECT
DL(I) = Longitudinal dispersion coefficient for channel segment I
(m2/sec) [Not in Appendix E]
136
-------
APPENDIX D
PROGRAM LISTING FOR THE MINNESOTA
STREAM WATER TEMPERATURE MODEL
MNSTREM
137
-------
00001
00002C
00003C
00004C
00005C
00006C
00007C
00008C
00009C
00010C
00011C
000 12C
00014C
00015C
0001 AC
00017C
0001 8C
0001 9C
00020C
00021C
00022C
00023C
00024C
00026C
00027C
0002BC
0002°C
00030C
0003 1C
00032C
00033C
00034C
00035C
00036C
00037C
00038C
00039C
00040C
00041C
00042C
00043C
00044C
0004SC
00046C
00047C
00048C
00049C
" 00050C
00051C
00052C -
00053C
00054C
00055C
OOOSoC
00057C
ooosac
00059C
00060C
00061C
00062C
00063C
00064C
00065C
PROGRAM MNSTREM
ALT = ALTITUDE(M)
B(K) = SURFACE WIDTH (M)
CC(I) = FFACTION CLOUD COVER FOR TIME INCREMENT I
CHDATE = THE NUMBER OF DAYS IN THE STARTING MONTH /100
DATE = MONTH DAY
DAY1 = STARTING DAY OF YEAR
DELX = DISTANCE INCREMENT FOR COMPUTATIONS
DL = LONGITUDINAL DISPERSION COEFFICIENT
DSTAR = DIMENSIQNLESS LONGITUDINAL DISPERSION COEFFICIENT
DTIME = TIME INCREhENT FOR COMPUTATIONS (SEC)
.EKRttAX..= J1AXIH.UJi_£RRa&.IiUC-TJl CM*,NG,e IN INCFEr<5NT SIZE
ERRMEAN- = MEAN ERROR DDUE VO LHANGc IN- INCREMENT SIZE""
FI s-COEFFIUENT UHZCH. Sr£C:?IES..TJCtl£. STEP METHOD ,.
FI =-05 FOI\ CRANtv-NieOLSGh' Irt.'-LICIT hETHOD
FI- = IX) FOR FULLY IMPLICIT METHOD - - . " :;_."._.
HEAD = HEADING FOR WEATHER DATA OUTI-UT
ITEST = 0 IF NO INCREMENT SIZE TEST HAS OCCURED
LAT = LATITUHE(RADIANS)
M= TOTAL NUMBER OF RESIDUALS
MHOUR= NUMBER OF HOURS BETWEEN WATER TEMPERATURE AND
UEATHER DATA INPUT
NITER = NUMBER OF ITERATIONS ON SOURCE TERM FOR EACH TIME STEP
Nf= NUMBER OF FLOU RATES READ IN
NHOUR= TOTAL HOURS OF TIME PERIOD
NSECT = NUMBER OF SECTIONS IN lOTAL CHANNEl LENGTH
NT = TOTAL NUMBER OF TIME INCREMENT STEPS
F= AIR ^RESSUFE (MB)
FRNOUT = A SUBROUTINE FOR OUTPUT OF WATtR TEMPERATURES
NSFTN
138
-------
\
00066C Q = CHANNEL FLOU RATE
00067C RAD= INCOMING SHORTWAVE RADIATION , RESID(4000)
00081 COMMON/C/ TSTART.TLNGTH, IFLAG, JFLAG.M-LAG.JTIME,NSECT,NDELEAT NSFTN
00082 COMMON/D/ T0(710),TI(18)
00083 COMMON/E/ Q<30),IF(30),NF
00084 COMMON/F/ TA(250) , RH(250>. RAD(250).U<250),UDIR(250),P(250>,CC<250)
00085 COMMON/G/ S(17),SIGMA, DAY1, LAT,ALT,TOLIK 150)
00086 COMMON/H/ A(17).B(17),AK(17),BK(17),CK(17),DL(17), CELX,
00087 tDTIME ,FE(17),FC(17),FI
00088 DIMENSION HEAD<10)
00089 REAL LAT
00090 EXTERNAL TEMP.FRNOUT
00091C READ PARAMETERS AND FLAGS
00092 READ(7,800) DSTAR, TSTART, DELX, DTIME,DAY1,TLNGTH,LAT,ALT
00093C DELX SHOULD BE CHOSEN SUCH THAT TLENGTH/NSECT/DELX = EVEN NO
00094 READC7.800) DATE,CHDATE,FI
00075 LAT-LATH3 14/180
00096C CONVERT TO METRIC
00097 OELX=DELX* 3043
00098 TLNGTH=TLNGTH* 3048
00099 READ(7,810) NHOUR, MHOUR. NF, NSECT.IFLAG.JFLAG,KFLAG.LFLAG.NITER
00100C IF IFLAG=0, NO OBSERVED WATER TEMPERATURES ARE READ
00101C IF JFLAG=0, THE LONGITUDINAL DISPERSION MODEL IS USED ' •
00102C IF JFLAG = I/ COMPUTATIONS ARE BY UPSTREAM DIFFERENCES
00103C (NUMERICAL DISPERSION ONLY)
00104C IF KFLAG=0, CSTAR, RATHER THAN DL, IS SPECIFIED
00105C IF LFLAG=1, TEST OF INCREMENT SIZE IS CALLED
00106 SIGMA=1 17E-7
00107 NT=NHOUR/MHOUR
00108 ITEST=0
00109 ITER=0
00110 NDELEAT=0 NSFTN
00111C READ INITIAL UATER TEMPERATURES
00112 READ(7,320> T0(1),(TI(I),1=1,17)
00113C READ AND CONVERT UNITS OF WEATHER DATA
00114 READ(7,030) (HEAD(I),1=1,10)
00115 PRINT 835, (HEAD(I),1=1,10) NSFTN
00116- PRINT 840_ - ._.-"- ..-..•-' NSF-TN
00117 . TIME=TSTART+FLOAT(MHOUR>
00118 DO 4 J=1,NT "
00119 NH=J*MHOUR
00120 ITER=ITER+1
00121 IFdTER EQ 87) PRINT 840 NSFTN
00122 IF(TIME GE 24 ) TIME=TIME-24
00123 READ(7,850) TA(J), RH(J), RAtK J),UiJ),UDIR(J>.P(J),CC.TIME NSFTN
00126 TIME=TIME+FLOAT(MHOUR>
00127 4 CONTINUE
00128 IF(IFLAG EQ 0> GO TO 7
00129C READ UATER TEMPERATURE TEMPERATURE DATA
00130 IT=0
.00131 DO 5 J=1,NT..
139
-------
APPENDIX G
MERS WEATHER STATION
The MERS channel water temperatures are controlled by inflow and
weather conditions. A station to measure and record weather parameters was
therefore installed near the center of the field station. Solar radiation
and wind velocity were recorded since December 3, 1975. A 50 - junction Epply
pyranometer and a Science Associate's anemometer were mounted at 2 m and 9 m,
respectively, above the ground. A Belfort hydrothermograph, with a bimetal
temperature sensor and hair hydrometer, were mounted in a weather shelter
and recorded air temperature and relative humidity since January 29, 1976.
A Texas Electronics wind direction vane and transmitter, mounted at 4 m, with
power supply unit and Rustrak recorder was put in operation on April 2, 1976.
Sample strip chart records are shown in Figs. G-l and G-2. The weather station
is shown in Fig. G-3. A one-year cycle of measured air temperatures is shown
in Fig. G-4.
157
n ' co
CD ip Cj3 C 5 CD CJ> CJ
J^fcHt-F-rrT I ' i . i—j
! * *' i •' I I ' i ' I i I
^trlJ^jlTm
^itr! rrtTTTrrr
ir^rmTnlTTi
b'lJJj.vJ-L.L
'3!-JX-{J4-
I ! !'
V,1 I I i i'
IXTI^!.! _[:.
j). ->_,_1'_1.^
• T_"'
1 *»':•
. •<.<*•
——>x-»i-
> • .V-
,• ;Jr"'
' II
_LH
o. --r--:-;-1--
* 11 '.^,'. '
O1 ' ' I , ! __.T-'}''^ x' ! ' I I ' II
, „ ;.,,vvi—K. 1 —I'r-
o'-iiilL-^"1.'-!:1..!:..'::!'].^
**l i I _ «V.i-,- i ii ii '
Q^II'!^ i''1.^'.1!.:..'.^.!1:.1']!1.
^iT?l'r:'!ii'^"!!^^
ol-^ll1!!^
-------
' ! MONDAY v TUESDAY v WEDNESDAY v THURSDAY v
1 ' 4- S 8 lO^IJa 4-68 lO^Jta 4-68 lOXIf? 468 lofjl 2 468 /qXJfS 468 10WS f 6 Q IQXffS 4-68 10^2 » 6 t
^ —~—————~———————————————~~ ———• '~
si
Fig. Gr2 - Sample hygrothermograph strip chart for degrees £e1c'ius air temperature
(top) and percent relative humidity (bottom). *
-------
r
Fig. G-3. Weather station with hygrothermograph and shelter,
pyranometer, wind vane, and anemometer. Top - complete
weather station. Bottom - hygrothermograph, small
weather shelter, and pyranometer.
160
-------
APPENDIX H
SAMPLE OUTPUT FROM PROGRAM WTEMP1
WATER TEMPERATURE- STATISTICS
********
WATER TEKPERATUPE STUDY
US EPA NONTICELLO ECOLOGICAL RESEARCH STATION
"CONDUCTED 8Y~"
UNIVERSITY Of (IIMNSSOTA
ST. ANTHONY FALLS HYDRAULIC LABORATORY
STATISTICAL ANA'LYS'IS OF "OPERATIONAL" MATE'R TENPERATURES FOR "
CHANNEL *» STATIQN~~2
NOTES»
»ERICD CT ANALYSIS IS FROM 10/20/76 TO 9/19/77. ALL STATISTICS
DECOROED ARE SLIDING 7-DAILY PARAMETERS. COMPUTED ft OH THREE-HOURLY
DiTt. ALL STATISTICS, EXCEPT SKEVNESS* ARE RECORDED IN DEGREES CEL-
CIU5. SKEVNESS IS OInENSIONLESS. ALL CALENDAR WEEKS START ON SUNDAY.
162