URBAN STORM WATER RUNOFF
DETERMINATION OF VOLUMES AND FLOWRATES
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
NERC - CINN
EDISON, N.J.
-------
URBAN STORM WATER RUNOFF
DETERMINATION OF VOLUMES AND FLOWKATES
BY
Van Te Chow and Ben Chie Yen
Department of Civil Engineering
University of Illinois at Urbana-Champaign
Urbana, Illinois 61801
Contract No. 68-03-0302
Program Element No. 1BB034
Project Officer
Chi Yuan Fan
Storm and Combined Sewer Section (Edison, N.J.)
Advanced Waste Treatment Research Laboratory,.
National Environment Research Center
Cincinnati, Ohio 45268
NATIONAL ENVIRONMENTAL RESEARCH CENTER
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
CINCINNATI, OHIO 45268
-------
FOREWORD
A prerequisite for effective control of storm runoff pollution
is a reliable method to predict the ..quantity of the storm runoff. The
time distribution of storm runoff from an urban drainage system depends
on the areal and temporal distributions of the intensity of the rainfall,
the frequency of the rainstorm, and the physical characteristics of the
drainage system. Numerous methods have been proposed to evaluate urban
runoff from rainfall. Many have been accepted for engineering
applications while others need yet to be tested and verified. Therefore
an investigation to evaluate the methods on a common basis would be a
significant contribution to the recent efforts on pollution control.
11
-------
ABSTRACT
An investigation is made to (a) develop a method of depth-
duration-frequency analysis for precipitation events having short return
period (high frequency) for urban storm water runoff management and control
purposes; (b) develop a new high accuracy urban storm water runoff de-
termination method which when verified, can be used for projects requiring
high accuracy detailed runoff results and can also be used as the
calibration scale for the less accurate urban runoff prediction methods;
and (c) compare and evaluate selected urban storm water runoff prediction
methods. The eight methods evaluated are the rational method, unit
hydrograph method, Chicago hydrograph method, British Road Research
Laboratory method, University of Cincinnati Urban Runoff method, EPA
Storm Water Management Model, Dorsch Hydrograph Volume method, and
Illinois Urban Storm Runoff method. The comparison and evaluation is
done by using four recorded hyetographs of the Oakdale Avenue Drainage
Basin in Chicago to produce the predicted hydrographs by the methods and
the results are compared with recorded hydrographs. The relative merits
of the methods are discussed and recommendations are made.
iii
-------
TABLE OF CONTENTS
Page
Foreword ii
Abstract iii
List of Figures vi
List of Tables viii
Acknowledgments ix
I. Summary and Conclusions 1
II. Recommendations 3
III. Introduction 7
1. Problems of urban storm water runoff 7
2. Objectives and scope of study 9
IV. Precipitation Analysis 11
1. Rainfall depth-duration-frequency analysis for urban runoff 11
2. Infiltration and other abstractions 47
3. Snow melt 51
V. Characteristics of Oakdale Avenue Basin 53
1. Surface drainage pattern 53
2. Sewer system 59
VI. Surface Runoff Model 66
1. Runoff in sub catchments 66'
2. Gutter flow routing 75
3. Inlets . 77
4. Program description and data preparation . 81
VII. Sewer System Routing Model 93
1. Sewer network representation 93
2. Method of solution 95
3. Computer program description 100
IV
-------
Page
VIII. Water Quality Model 105
1. Water quality model formulation 106
2. Program description and data preparation 109
IX. General Description of Other Methods Evaluated 115
1. The rational method 118
2. Unit hydrograph method 118
3. Chicago hydrograph method 120
4. Road Research Laboratory method 121
5. Cincinnati urban runoff model 123
6. EPA storm water management model 126
7. Dorsch hydrograph-volume-method 128
X. Evaluation of Methods 130
1. Hydrographs for the Eight Methods 131
2. Comparison of the methods 148
XI. References 157-161
XII. Notation 162-163
Appendices
A. Listing of computer program for frequency 164
analysis of hourly precipitation data
B. Listing of Computer Program for the Illinois 202
Surface Runoff Model
C. Listing of Computer Program for the Illinois - 219
Sewer System Water Quality Model
-------
FIGURES
No. Page
1 Example Hyetograph 17
2 Trapezoidal Representation of Hyetograph 17
3 Probability Distribution of Rainstorm Depth 21
4 Probability Distribution of Rainstorm Duration 22
5 Probability Distribution of Average Rainstorm 23
Intensity
6 Probability Distribution of Elapse Time Between 24
Rainstorms
7 Conditional Distributions of Average Rainstorm 27-31
Intensity
8 Gamma Density Function Parameters as Functions 38
of Rainstorm Duration
9 Exponential Density Function Parameter as 39
Function of Rainstorm Duration
10 Non-Exceedance Probabilities for Rainstorm Depth 42-43
and Elapse Time Between Rainstorms
11 Conditional Distributions of Rainstorm Duration 44-45
and Hyetograph First Moment Arm
12 Oakdale Basin Location Map 54
13 Drainage Pattern of Oakdale Basin 55
14 Schematic Drawing of Gutter Cross Section , 58
15 Details of Grate Inlets 60-61
16 Circular Sewer Flow Cross Section 63
17 Evaluation of Weisbach Resistance Coefficient • 70
18 Computational Grid for Semi-Implicit Four-Point 73
Backward Difference Scheme
19 Flow Chart for Illinois Surface Runoff Model 82
Computer Program
20 Composition of Illinois Surface Runoff Model 83
Computer Program
VI
-------
No.
21 Identification of Types of Gutter and Inlet 88
22 Example of Illinois Surface Runoff Model 91-92
Data Preparation
23 Solution by Method of Overlapping Y-Segments 99
24 Flow Chart for ISS Model Computer Program Flow 101
Prediction Option
25 Flow Chart for Illinois Sewer System Water 112
Quality Model Computer Program
26 Hyetograph and Hydrographs for May 19, 1959 132
Rainstorm
27 Hyetograph and Hydrographs for July 2, 1960 133
Rainstorm
28 Hyetograph and Hydrographs for April 29, 1963 134
Rainstorm
29 Hyetograph and Hydrographs for July 7, 1964 135
Rainstorm
30 Ten-Min Unit Hydrograph for Oakdale Basin 139
31 S-Curve and One-Min Unit Hydrograph for Oakdale 140
Basin
32 Sensitivity of Computed Hydrograph to Horton's 151-152
Infiltration Parameters
vii
-------
TABLES
No. Page
1 Statistics of Parameters for Summer Rainstorms 19
at Urbana, Illinois
2 Conditional Statistics of Rainstorm Parameters 32-36
for Different Durations
3 Dimensions of Gutters of Oakdale Avenue Drainage 64
Basin
4 Dimensions of Alleys of Oakdale Avenue Drainage 65
Basin
5 Dimensions of Sewers of Oakdale Avenue Drainage 65
Basin
6 Urban Runoff Prediction Methods Evaluated 116-117
7 Rational Method Computation 136
8 Comparison of Urban Storm Runoff Methods 137
9 Headings for Computation of Unit Hydrograph 142
Method
10 Values of Infiltration Parameters Used for Illinois 146
Urban Storm Runoff Method shown in Figs. 26 to 29
viii '
-------
ACKNOWLEDGMENTS
)
During the period of this project, many people have contributed
to the progress of the work. Many students and assistants at the University
of Illinois at Urbana-Champaign have helped in various phases of data
collection, card punching, data analysis, and numerous other computations.
Many engineers in various federal and local government agencies, consulting
firms, research institutes and universities provided indispensable help in
supplying data, information exchanges, and discussions, although only
selected data are used in this report and not all the ideas provided by
these people are presented here. The names of these people are too numerous
to be listed here and their contributions are greatly appreciated. Particular
thanks are due Professor Wayne Huber of the University of Florida, Dr. Albin
Brandstetter of Battelle Pacific Northwest Laboratories, and Dr. P. Wisner of
James F. MacLaren, Ltd. for their help.
Mr. A. Osman Akan, Graduate Research Assistant in the Department
of Civil Engineering, University of Illinois at Urbana-Champaign, contributed
most to the successfulness of this research project. He programmed the
Illinois Surface Runoff Model and the Water Quality Model, provided the
interface between these models and the Illinois Storm Sewer System Simulation
Model, performed the tedious and painstaking testing and computer computations,
participated in the field survey for data collection, and wrote the section on
Program Description and Data Preparation for the Illinois Surface Runoff
Model. Without his help, this report would not be materialized. Mr. T. A.
Ula, Graduate Research Assistant, performed the frequency analysis in addition
to his contributions in data collection and various aspects of computations.
Throughout this research project, Mr. Chi-Yuan Fan, Project Officer,
and Mr. Richard Field, Chief, Storm and Combined Sewer Section, EPA Advanced
Waste Treatment Research Laboratory, provided smooth and fruitful cooperation.
ix
-------
Their help in supplying information on data sources and coordination with other
researches is indispensable and deeply appreciated.
Thanks are also due Mrs. Norma Barton for her typing of the report
and other materials for the project.
x
-------
I. SUMMARY AND CONCLUSIONS
The purpose of this report is to provide information which is
useful for operation and management of urban storm and combined sewer
runoff and also useful for design and operation of overflow treatment
and other control facilities. This report is divided into three major
parts. The first part contained in Chapter IV is precipitation analysis
on depth-duration-frequency analysis of short return period (high
frequency) rainstorms. Conditional probability is utilized and method
of application on urban runoff problems is demonstrated using the hourly
precipitation record available from the U.S. Environmental Data Service,
National Climatic Center.
The second part of the report consists of Chapters VI, VII,
t
and VIII describing the Illinois Urban Storm Runoff method. It
includes the development of the Illinois Urban Surface Runoff model to
couple with the existing Illinois Storm Sewer System Simulation Model
and the formulation of a non-reactive water quality model to compute
the concentration of the pollutants of urban runoff. In the surface
runoff model, Morton's formula is used to evaluate infiltration. The
overland flow is computed by using the kinematic wave method together
with Darcy-Weisbach's formula to estimate the friction slope. The
gutter flow is computed by using the kinematic wave method together
with Manning's formula to estimate the friction slope. Inlets to
catch basins are classified in types according to their geometry, and
weir and orifice formulas are used to estimate the discharge. The
sewer flow routing utilizes the complete St. Venant equations
accounting for the junction backwater effects.
-------
The third part of the report consisting of Chapters V, IX,
and X is a comparative evaluation of eight selected urban runoff
prediction methods using the data from the Oakdale Avenue Drainage
Basin in Chicago. The methods evaluated are the rational method, unit
hydrograph method, Chicago hydrograph method, British Road Research
Laboratory method, University of Cincinnati Urban Runoff method, EPA
Storm Water Management Model, Dorsch Hydrograph Volume Method, and
the Illinois Urban Storm Runoff method. This part of study is believed
to be of particular interest to practicing engineers. The
evaluation is made by using the recorded hyetographs of four rainstorms
applying the methods to compute the predicted runoff hydrographs and
then compare the results with the recorded hydrographs. This
comparative study suggests that the most suitable method to be used
for an urban runoff problem depends on the accuracy required for the
project. If only a quick simple approximate result of peak runoff
rate is needed, the rational method is quite satisfactory, whereas for
a project involving a large amount of money and high accuracy and
details of the temporal and spatial distribution of the runoff are
required, the Illinois Urban Storm Runoff method will be a suitable
choice, and the Dorsch method and EPA SWMM may be the alternatives
if backwater effects are not important. For in-between accuracy,
the unit hydrograph method is recommended, if possible. When the
unit hydrograph for the drainage area is not available, in most
cases the Road Research Laboratory method appears to be superior to
the Chicago and University of Cincinnati methods.
-------
II. RECOMMENDATIONS
Based on the results of this investigation, the following
recommendations related to the determination of flow rate of urban
storm water runoff are made:
1. The most suitable method to be used for determination of
urban storm runoff depends on the objective and size of
the project, the required accuracy and detail of the
flow rate determination, and the available data.
2. For the purpose of storm water runoff control and
management, when required results are extensive (such as
discharge and depth or velocity at different times and
at many locations of the drainage basin) and the required
accuracy is high such as in the case of a high-valued
urban area or a high-cost project, and the detailed basin
data is available, the Illinois Urban Runoff method appears
to be a suitable choice. The Dorsch hydrograph-volume
method may also be used. If the backwater effect of
sewer junctions is not important, the EPA SWMM is a
good alternative.
3. For a cheap, simple and quick estimation of the peak
runoff rate without requiring the entire runoff hydro-
graph, the rational method can often be used
satisfactorily. However, one should always bear in
mind the limited accuracy the rational method can provide -
a sacrifice of accuracy for the sake of simplicity.
4. For projects or problems requiring moderate accuracy
and determination of the entire or part of the runoff
hydrograph, the unit hydrograph method is the simplest
-------
and cheapest to use if the unit hydrograph is available or
can be reliably derived. In using the unit hydrograph,
one should check that there should be no significant
change in basin characteristics. If the unit hydrograph
is not available, the order of choice would be the British
Road Research Laboratory method, Chicago hydrograph method,
and University of Cincinnati Urban Runoff Model.
5. The comparative study is made on only eight methods and is
neither exhaustive nor exclusive. Other methods not .
mentioned here may also be .useful.
6. Existing and available data on urban rainfall, runoff, and
basin characteristics are generally inadequate for a re-
liable accurate evaluation of the methods. Neither are
the data adequate for engineering operation, management
and design purposes. This inadequacy is in both details
and accuracy. For example, the best readily available
rainfall data, the hourly data from the Environmental
Data Service, National Climatic Center, is inadequate in
view of many urban drainage basins having a time of
travel much shorter than an hour.
7. More study on urban infiltration would be desirable. This
includes the determination and listing of infiltration
parameters such as f and k in Horton's formula for"
typical urban surface conditions, and data on informa-
tion such as anticedent moisture condition related to
infiltration.
8. A detailed accurate evaluation of the surface runoff
part of the methods is desirable. Because of the lack
-------
of sufficient detailed geometric (percentage and
distribution of roofs, lawns, sidewalks, driveways, etc.)
and infiltration data, it is believed that the effects
of the surface runoff on the methods have not been
adequately evaluated. In addition, surface runoff hydro-
graphs at certain urban locations will also be useful
for engineering purposes.
9. Further study on hydraulic characteristics of inlets is
necessary. There is no use to have a highly sophisticated
surface runoff model if its downstream control, the
inlet, cannot be accurately modeled. Existing informa-
tion is inadequate to represent the large number of types
of inlets now being used.
10. Improvement on the manner to handle surcharge and
supercritical flow in sewers would be useful and desirable.
11. Further study on differences between design and flow
prediction purposes on modeling and data requirement may
be useful for urban storm runoff management.
12. In view of the probabilistic nature of the physical
conditions of the drainage system, such as clogging of
inlets and gutters, interference of parked vehicles on
street and gutter flows, change of roughness of sewers,
etc., development of a probabilistic method to account
for such uncertainties will be useful for both design and
operation purposes.
13. It is of course possible to further improve the hydraulic
aspects of the existing sophisticated methods. For
instance, the kinematic wave routing of the overland and
-------
gutter flows can be replaced by solving the more accurate
St. Venant equations. However, at present such im-
provement appears to be unfruitful and immature because
of the uncertainties involved in the basin physical
properties and the detail and reliability of the data.
Furthermore, such improvement would require considerably
more computer time, thus making the method impractical.
14. An advanced stochastic approach to analysis and predicted
high frequency rainfall as an alternative to the method
proposed in Sec. IV-1 may be useful.
-------
III. INTRODUCTION
III-l. Problems of Urban Storm Water Runoff
Metropolitan areas, in the United States as well as elsewhere in
the world, are growing at an unprecedented rate. One result often associated
with urbanization is the deterioration of the living environment due to
either lack of comprehensive planning or incapability of the cities to keep
pace with the growth. Among the vital facilities in preserving the living
environment, urban sewer drainage systems affect directly the quantity
and quality in disposing urban waste water. Large amounts of money and
resources are involved in the design, construction, modification, operation
and maintenance of urban sewer systems. There are two major types of
problems related to urban storm water runoff: (1) flooding due to in-
adequate sewer capacity causing damage of properties and disruption of
traffic and other human activities; and (2) pollution due to storm runoff.
The flooding problem is a design problem involving design rainstorms having
return period of once in several years. Many methods have been
developed for sewer design purpose in the past 130 years (Chow, 1962, 1964).
The common objective of these methods is to provide a design flow for
sizing the sewers of a new storm drainage system or an existing system
with adequate capacities to dispose of this once-in-many-year design flow
without flooding.
Conversely, the storm runoff pollution problem is an operation
and management problem. It involves consideration of rainstorms with
frequencies of several times in a year. For an existing drainage system
the problem involves the management of the time distribution of the
quantity and quality of storm water so that the runoff would not overload
-------
the treatment facilities or unacceptably pollute the receiving water bodies.
This may require some modification of the drainage facilities. For new drainage
systems a balanced design considering both the pollution control and adequate
capacity to avoid flooding would require an optimization analysis.
During the past decade the public as well as the engineering
profession has been greatly alarmed by the pollution aspects of urban storm
runoff. Many investigations have been sponsored by EPA and by other
agencies on the extent of urban runoff pollution (e.g., see American Public
Works Association, 1967, 1969; Engineering-Science, Inc., 1967; Envirogenics
U"'" L—' I
Company, 1971; Hawkins and Judd, 1972; Sartor and Boyd, 1972; Weibel et al.,
*•"•
1964). Recent studies on urban storm runoff pollution and efforts on its
control and management have been summarized in two excellent literature
i— u^-
reviews (Field and Weigel, 1973; Field and Szeeley, 1974).
One fundamental prerequisite for an efficient and economic design
and operation of urban sewer systems is a reliable method to predict the
quantity and quality of water handled by the system, particularly the time
distribution of runoff due to rainfall at important locations such as
junctions and overflow facilities. Numerous methods have been proposed
to estimate rainstorm runoff. Many of these methods are of regional
nature whereas many others are applicable only to rural or similar natural
drainage basins. For those methods applicable to urban areas, some treat
i
the drainage area as a "black box" without considering the time or space
distribution of the runoff in the drainage system; others treat the
drainage system as sequential overland and channel flows without considering
the detention storage in the drainage system due to backwater effects of
the junctions. Furthermore, many of the improved methods which have been
proposed recently are not yet widely accepted mainly because their relative
-------
merits have not been assessed, nor have they been compared in a comprehensive
manner on a common basis to the previously developed methods.
However, in order to comply with recent pollution control laws,
attempts have been made to utilize the storage capacity of a sewer system
to detent storm runoff and to control overflow in order to reduce the cost
in handling urban waste water problems (Poertner, 1972; Anderson et al., 1972;
Field and Struzeski, 1972). Examples of such attempts are the automatic
control systems built or proposed in Minneapolis-St. Paul (Anderson, 1970;
Tucker, 1971), Seattle, San Francisco, and Detroit (Field and Struzeski,
1972). Obviously, a reliable storm runoff prediction method would be
particularly useful and beneficial for such efforts. In fact, without a
reliable runoff quantity prediction it is most unlikely that a satisfactory
runoff quality prediction can be achieved.
III-2. Objectives and Scope of Study
As discussed in the preceding section, a "demonstration" type
investigation to evaluate quantitatively on a common basis the relative
merits of the conventional as well as recently developed runoff prediction
methods is needed for improvement in urban storm water management and
pollution control. Since many of the urban storm runoff models are being
,-
compared quantitatively in an EPA project at Battelle Pacific Northwest
Laboratories using hypothetical data, it is not the intention of this
investigation to evaluate all of the existing urban storm runoff models.
Accordingly, the major objectives of the present study are: (a) to
develop a surface runoff model, coupled with an existing sewer routing model
previously developed at the University of Illinois, the Illinois Storm Sewer
System Simulation Model, to form an urban storm runoff method; and
(b) to evaluate this new method using actual field data and also compare
-------
the results quantitatively with other relatively popular methods using the
same field data from the Oakdale Avenue Drainage Basin in Chicago. The
methods compared in this study include the rational method, unit hydrograph
method, Chicago hydrograph method, British Road Research Laboratory method, EPA
Storm Water Management Model, University of Cincinnati method, Dorsch Hydrograph-
Volume method, and the Illinois Urban Storm Runoff method. The description of
the last method is given in Chapters VI and VII. Brief descriptions of the
other seven methods are given in Chapter IX.
2
In addition to the Oakdale Drainage Basin (0.052 km or 12.9 ac)
which was chosen to test the eight chosen methods because of its available
data and previous studies, a much bigger basin, the Boneyard Creek
2
Drainage Basin in Champaign-Urbana, Illinois (9.3 km or 3.0 sq mi) was
also used to test the applicability of the Illinois Urban Storm Runoff
method on large basins. However, because of the enormous cost, time, and
manpower that would be involved if the other methods were also tested on
the Boneyard Basin, and the result would probably produce no additional
information than that from the Oakdale Basin, the other seven methods were
not tested on the Boneyard Basin. Furthermore, because of the large amount
of data for the Boneyard Basin due to its size, its inclusion here would
»-
make this report voluminous. Therefore, its physical properties together
with the results to demonstrate the applicability of the Illinois Urban
Storm Runoff method on large basins may later be presented as a separate
supplementary report.
Also, an auxiliary objective to objectives (a) and (b) mentioned
above is (c) to develop a method of rainfall depth-duration-frequency
analysis suitable for urban storm runoff management purposes.
10
-------
IV. PRECIPITATION ANALYSIS
Not all the rainwater falling on an urban area becomes runoff.
There are infiltration and other losses called abstractions. The total
precipitation subtracted by the abstractions is called precipitation excess.
Since both precipitation and abstractions are statistical quantities,
precipitation excess is also of statistical nature. In this chapter
precipitation and abstractions are discussed in view of urban storm runoff.
IV-1. Rainfall Depth-Duration-Frequency Analysis for Urban Runoff
Natural rainfalls have finite duration and areal and temporal
variation of their intensities. It has yet to be found that two rainfalls
are identical. Thus rainfall data is analyzed statistically to be useful
for engineering purposes. The intensity of rainfall is a function of its
duration, frequency, and area, which has been discussed extensively
elsewhere (e.g., U.S. National Weather Service, 1961; Chow, 1964). There
are two types of rainfall information needed for urban storm runoff manage-
ment purposes. The first is the duration and maximum intensity for rainfalls
having long return periods of a number of years to be used for design and
safety considerations. The other is the information on high frequency
rainstorms with return periods less than a year to be used for operation
and pollution control purposes.
Because of the large number of rainfalls involved for a given
location, conventionally in engineering hydrology as well as in meteorology
only the maximum values in the form of partial duration series or annual
maximum series are analyzed to establish the rainfall intensity-duration-
frequency relationship for long return period events. This is done so
because usually for the purpose of design of drainage facilities based on
catastrophic failure concept (Yen and Ang, 1971) rainfall of small return
11
-------
period is usually of no significance. However, this is not the case from
an operational viewpoint for the control of pollution due to urban storm
runoff for which every rainstorm contributes to the problem. Because of
the lack of dillution effect due to their small volume, the low return
period rainstorms probably contribute considerably more to pollution
per unit volume of water than the long return period ones. Treatment
plants, detention storages, overflows and other sewer facilities designed
for small capacities corresponding to most frequent rainfalls would not
be able to control the runoff from less frequent rainstorms. Contrarily,
such facilities designed to operate at full capacity only once in every
several years would be costly and unjustified if the safety consideration
does not require it.
Obviously, the volume and time distribution of storm runoff
quantity and quality from an urban drainage system depend on the areal
and temporal distributions of the rain falling on the urban basin.
Consequently, the depth, duration, and frequency of the rainfall and
other parameters defining the internal pattern of the rainstorm should
all be considered. In fact, the time elapsed between successive rain-
falls is also an important factor in determining the quantity and
especially the quality of storm runoff. Conceivably, most of the water
from a rainstorm followed soon after an earlier heavy rainfall would
become runoff and the quality of the water would be relatively better.
Information on long return period rainfall useful for urban
runoff studies have been well established and can readily be found
(e.g., Chow, 1953; U.S. Weather Service, 1961). Unfortunately,
analytical information on high frequency rainstorms with return periods
less than a year which is useful for urban engineering purposes is
practically nonexistant. This lack of information is due mainly to the
12
-------
large amount of data needed to be analyzed and to the difficulty in
defining precisely the duration and intensity of a rainstorm. The purpose
of this part of the study is to provide a practical and realistic model
on a probablistic basis for rainfall input for urban storm runoff studies
and operations.
The frequency analysis methods for rainfalls of long return
period can be extended to short return period rainfalls. However, unlike
the long return period case for which data can be selected to form a
partial series for the analysis, for short return period analysis the
entire set of data, i.e. , all the rainstorms recorded, are utilized.
Grayman and Eagleson (1969) applied this concept on hourly rainfall
data for 546 rainstorms in a 5-year period at Boston.
Ideally, the precipitation data used for short return period
frequency analysis should be a continuous record of hyetograph (rainfall-
time curve). For most precipitation recording-gaging stations in the
United States, the most detailed precipitation data readily available
from the Environmental Data Service, National Climatic Center* are
hourly records. Since many urban drainage basins have the time of travel
of surface runoff less than an hour, the hourly rainfall data is
obviously inadequate and unsatisfactory. Data of 5- oj: 10-min
intervals would be much more satisfactory. Unfortunately, rainfall
records having time intervals shorter than one hour may. only be
obtained by the user from the original record, if available, which is
a very time consuming process and most unlikely to be undertaken by
*Hourly rainfall data for recording raingages in the U.S. Weather Service
system are available at cost on punched cards from U.S. Department of
Commerce, National Oceanic and Atmospheric Administration, Environmental
Data Service, National Climatic Center, Federal Building, Asheville,
N.C. 28801. , .
13
-------
practicing engineers. Therefore, from a practical viewpoint, this part
of research on rainfall frequency analysis method is developed by using
hourly rainfall records. Nevertheless, the methodology is equally
applicable to data of other time intervals.
The method of analysis which has been written in a computer
program (Appendix A) will be described in the following in this section.
The data used to demonstrate the application of the method are the
point hourly precipitation data on punched cards as provided by the
National Climatic Center for the Morrow Plot raingage at Urbana, Illinois,
(11-8740 Urbana), covering 14 years from 1959 to 1972. Because of the
seasonal characteristics of rainstorms, the analysis is carried out separately
for different seasons. The analysis for the record of June, July, and
August consisting of 455 rainstorms is presented here as an example.
(A) Identification of rainstorms. - A rainstorm is defined here as a
period of continuous non-zero rainfall. However, since the input data
provided on punched cards do not identify traces (rainfall less than
0.01 in./hr or 0.25 mm/hr), there is no way to differentiate traces
from hours of no rainfall in the data. The duration of a rainstorm,
t, in hours, is defined as the length of time between the beginning and
the end of a continuous non-zero rainfall. The total volume of a
rainstorm as customarily expressed in depth, D in in. or mm covering the
entire area considered, is equal to the sum of the depth for each time
interval xrithin the duration of the rainstorm, i.e.,
D - ? d (1)
j=l J
in which d. is the depth for the j-th time interval and n is the number
of time intervals of the rainstorm. In engineering practice, as a matter
of convenience, equal time interval At is usually used and the standard
14
-------
time interval used here is one hour according to the standard U.S.
Weather Service data, although other time intervals can also be specified.
The average intensity of the rainstorm, i, is defined as equal to D/t,,
and usually expressed in in./hr or mm/hr. The elapse time between
successive rainstorms, t, , is defined as the time between the end of a
b
rainstorm and the beginning of the next rainstorm. The computer program
traces the entire record, identifies the rainstorms and determines the
values of D, i, t,, and t, .
a b
(B) Calculation of rainstorm parameters. - A schematic drawing of an
example hyetograph is shown in Fig. 1. The other parameters of the rain-
storms essential for a statistical analysis are computed as follows. The
standard deviation of the rainstorm depth, a, in in. or mm, is computed
as n „
I (d,-d)2
°d ' 1^=^—1 1/2 (2)
where the average depth per time interval, d in in. or mm, is
n
I d.
n
(3)
The first moment arm of the hyetograph with respect to the beginning
time of the rainstorm, t in hr, is
_ n n
t = At [ I (j-0.5) d.]/ Id (4)
2
and the corresponding second moment arm, G in hr , is
G=(At)2[£ (j-0.5)2d. +~ I d ]/ ? d (5)
L 3 3
15
-------
(C) Nondimensional hyetograph. - In order to describe in more general
terms the time distribution of the rainfall of a rainstorm, the
hyetograph is nondimensionalized by using the rainstorm depth D and
duration t_. as the nondimensionalizing parameters. Therefore, for a time
interval, the nondimensional depth d. = d./D where the superscript (o)
J J n
represents the nondimensional quantity. Accordingly, D = £ d. = 1
j=l 3
and t = 1. The nondimensional average intensity i = D /t = 1.
Similarly, the nondimensional average depth is
n
The nondimensional standard deviation of the depth .is
n 0
«°i - 5°> 1/2 a
The first moment arm of the nondimensional hyetograph is
t° = f- (8)
d
The second moment arm of the nondimensional hyetograph is
G° = -| (9)
fcd
(D) Shape of hyetographs. - The shape of the hyetographs may be approxi-
mated by some simple geometric figures. Assuming that the hyetographs
can be represented by trapezoids shown in Fig. 2,
t, = a + b + c (10)
d
16
-------
E
fc
o.
o>
o
o
o:
At
Time, t, in hr
Fig. 1. Example hyetograph
o.
o>
O
o
o:
Time, t
Fig. 2. Trapezoidal representation of
hyetograph
17
-------
t,(t, + a + c) + c(2a + c)
- 3(td+.)
and
j + (a+c) tj + (a+c)2 (t,+c) + ca(2a+c)
d d d
6(td + c)
Equations 10 through 13 can be solved simultaneously for a, b, c, and h
so that the trapezoid representing the hyetograph can be determined.
To demonstrate the methodology in application, it is assumed that a
special case of the trapezoidal shape with c = 0, i.e. triangles, can be
used to approximate the Urbana rainfall data. For this triangular case,
solving Eqs. 10, 11, and 12 yields
a = 3t - t, (14)
d
b = 2t. - 3t (15)
d
and
h = 2d (16)
For nondimensional hyetograph,
a° = a/t, = 3t° - 1 (17)
d
b° = b/t. = 2 - 3t° (18)
d
and
h° = h/D = 2d/D = 2/n " . (19)
The statistics of the parameters for the 455 summer rainstorms
at Urbana, Illinois are given in Table 1.
(E) Frequency analysis of rainstorm parameters. - With the rainstorm
parameters computed for every rainstorm in the record, a one-way
frequency analysis can subsequently be made for each parameter. The computer
18
-------
Table 1. STATISTICS OF PARAMETERS FOR SUMMER
RAINSTORMS AT URBANA, ILLINOIS
Parameter
tb , hr
td , hr
in.
D ,
Tnni
in./hr
i ,
mm/hr
in.
Q
Htfrt
t , hr
G , hr2
a , hr
b , hr
in.
h ,
o
t°
G°
0
a
b°
h°
Mean
59.4
2.45
0.31
7.87
0.11
2.79
0.06
1.52
1.14
2.90
0.95
1.50
0.21
5.33
0.13
0.47
0.30
0.42
0.58
1.23
Standard
Deviation
81.9
2.08
0.51
13.0
0.15
3.81
0.11
2.79
1.01
7.10
1.39
1.80
0.30
7.62
0.15
0.11
0.11
0.32
0.32
0.65
Min
1
1
0.01
0.25
0.01
0.25
0
0
0.50
0.33
-5.00
-2.10
0.02
0.51
0
0.13
0.05
-0.62
-0.40
0.14
Max
744
14
3.52
89.4
1.05
26.7
0.87
22.1
7.68
71.2
10.0
13.8
2.09
53.1
0.48
0.80
0.68
1.40
1.62
2.00
Number of Rainstorms = 455
19
-------
program (Appendix A) has a one-way frequency analysis subroutine which
tabulates for any given parameter the frequencies (number of observations
over given intervals), relative frequencies (frequency divided by the
total number of observations), probability densities (relative frequency
divided by the interval size) and non-exceedance probabilities
(cumulative relative frequencies). The mean and standard deviation of
the parameter are also calculated and the maximum and minimum values are
found (Table 1). Histograms of the probability densities for the rain-
storm parameters can then be plotted as shown in Figs. 3, 4, 5, and 6
for the rainstorm depth, duration, intensity, and elapse time between
rainstorms, respectively.
(F) Fitting of probability density functions. - The next step is to fit
some probability density functions to the histograms of the rainstorm
parameters. Exponential distribution and gamma distribution have
a non-negative range and have been applied in the analysis of
non-negative valued rainstorm parameters. The probability density
function for the exponential distribution is
f(x) = J e~X/B x ^ 0 and B > 0 (20)
D
where both the expected value and the standard deviation of the distri-
s
bution are given by B. This single-parameter distribution can easily
be fitted to data (Grayman and Eagleson, 1969) but the fact that it has
the same expected value and standard deviation makes its use difficult
to justify in most applications. The exponential density functions
fitted in Figs. 3 through 6 all assume a value of B equal to the
observed mean. Use of the mean rather than standard deviation for B im-
plicitly implies that the mean is more "meaningful" than the standard de-
viation. The probability density function for the gamma distribution is
20
-------
Rainstorm Depth, D, in mm
50
I I I I
Gamma Distribution with C = 0.37, 8=0.84
Exponential Distribution with 8 = 0.31
- „ I -. L
I 2 3
Rainstorm Depth, D, in in.
100
0.35
0.30
0.25
>>
u
0.20 O-
£
o>
_>
0.15 2
0>
or
o.io
o.os
Fig. 3. Probability distribution of rainstorm depth
21
-------
0.4
0.4
6 8 10 12
Rainstorm Duration, td, in hr
14
— 0.3
i i i i i i i i i i i i i i
Gamma Distribution with C = I139, 8=1,76
Exponential Distribution with B=2,45
u
c
a>
3
a>
>
a>
DC
16
Fig. 4. Probability distribution of rainstorm duration
22
-------
Ill
10
20
30
o>
o
JO
O
.O
O
Gamma Distribution with C = 0,54, B = 0,20 —
Exponential Distribution with B = O.II
n I In rrl m I
0.22
0.20
0.18
0.16
0.14
3
CT
a>
0.10 >
a>
Od
O.O8
0.06
0.04
0.02
0.2 0^ 0.6 0.8 1.0
Average Rainstorm Intensity, i, in in./hr
1.2
Fig. 5. Probability distribution of average rainstorm intensity
23
-------
•
0.03 -
^
«-»
*~
>»
1 1
s
\
1
L
y
i
I
s
h
m
£
f
JK
lflk»— nL 1 _ 1
200 400
J_
tb, in hr
600
1 „
0.3
^
u
c
a>
er
a>
U.
0.2
a>
•w
O
0
tr.
O.I
r>
800
Elapse Time Between Rainstorms, tb, in hr
Fig. 6. Probability distribution of elapse time between rainstorms
24
-------
1 C-l ~X/B
f(x) = - - x e x £ 0 and B, C > 0 (21)
F(C)BU
where the gamma function
r(C) = z" dz for all C > 0 (22)
J0
The expected value
E(x). = CB (23)
and the variance
V(x) = CB2 (24)
The exponential distribution is actually a special case of gamma dis-
tribution with C = 1. The gamma distribution being a two-parameter
distribution provides an extra degree of freedom compared to the
exponential distribution in fitting data, and both mean and standard
deviation of the observed data can be preserved by adjusting its para-
meters C and B. In Figs. 3 through 6, the values of C and B for the
fitted gamma probabity density function are computed by using the
observed values of mean and standard deviation:
2
" V(x)
. B - =£0. (26)
s
It can be seen from Figs. 3 through 6 that in general the gamma dis-
tribution provides a better fit to the data than the exponential
distribution.
(G) Conditional frequency analysis of rainstorm parameters. - As
the rainstorm parameters are not truly independent, conditional
probabilistic analysis is included to provide information useful for
solving urban storm runoff problems. The computer program has a
two-way frequency analysis subroutine which tabulates, for any two given
25
-------
parameters, a two-dimensional table of frequencies, relative frequencies,
and probability densities. The results can then be used to plot
three-dimensional histograms of joint probability densities for pairs
of rainstorm parameters. However, it is a formidable task to fit
bivariate joint density functions to three-dimensional histograms. Mainly
for this reason and also from a practical viewpoint, in this study joint
distributions are dealt with through the use of conditional distributions.
The computer program has a sorting subroutine which can sort the values
of a rainstorm parameter in an ascending order and can rearrange
simultaneously the corresponding values of another parameter. By using
this subroutine together with the one-way frequency analysis
subroutine, conditional frequency analysis can be carried out for pairs
of rainstorm parameters.
The example described here is to find the conditional distribu-
tions of average rainstorm intensity, i, for different rainstorm durations.
The same procedure can be applied to any other pair of dependent rainstorm
parameters. By using the sorting subroutine, values of t, are sorted in
an ascending order (here they vary from 1 to 14 hr), and corresponding
values of i are rearranged simultaneously. For a given value or a given
range of values of t,, the corresponding i values are picked up for a
one-way frequency analysis as described in (E). By repeating this one-way
frequency analysis of i for different values of t,, a set of histograms
of conditional probability densities of i are obtained as shown in
Figs. 7a to 7e. Exponential and gamma density functions were then
.fitted to these histograms with parameters evaluated based on the
observed conditional means and standard deviations. Tables 2 (a) to (e)
give the conditional means, standard deviations, minimum and maximum
values of each rainstorm parameter for different values of t,.
26
-------
20
Average Rainstorm Intensity, i, in mm/hr
10 20
T
T
Gamma Distribution with C = 0.22, 8=0.32
Exponential Distribution with 8 = 0,07
15
co
c
-------
CVJ
CO
C
0>
o
"o
c
o
o
o
Average Rainstorm Intensity, i, in mm/hr
10 20
Gamma Distribution with C =0,5, B = 0,24
Exponential Distribution with B = 0,I2
30
0,32
0,28
0.24
0.20
0,16
0.12
c
0>
3
er
CD
CD
>
CD
(E
0.08
0.04
ni In I n i
0,2 0.4 0.6 0,8 1.0
Average Rainstorm Intensity, i, in in./hr
—'o
1.2
Fig. 7. (b) t, = 2 hr
28
-------
ro
u
•o
in
c
0>
o
o
c
o
TJ
C
o
o
10
Average Rainstorm Intensity, i, in mm/hr
10 20 3O
Gamma Distribution with C = l!l9, B = O.IO
8|— Exponential Distribution with B = 0.12
0.2 0.4 0.6 0.8 1.0
Average Rainstorm Intensity, i, in in./hr
0.20
0.16
0.12
0.08
u
0>
cr
o>
0.04
1.2
Fig. 7. (c) t, = 3 hr
29
-------
Average Rainstorm Intensity, i, in mm/hr
10 20
Average Rainstorm Intensity, i, in in,/hr
Fig. 7. (d) t = 4, 5 hr
30.
•
in
^ a •
i y
« i
0 II
ll
1 '
_o
•5
o
0 2
0
1
Gamma Distribu
Exponential Dis
i'
^y
-
'
.
K
•
y
s
0.2
»
If
»
1 ' 1 ' |
tion with C = I.I5, B=0,I3 — _
— •
—
*-_J., , I ,"
0.4 0.6 0.8 1.0 1
0.16
U
a>
0.12 3
cr
a>
U_
a>
0.08 ^
IT
0.04
0
2
30
-------
Average Rainstorm Intensity, i, in mm/hr
14°
12
~
T"
(O
U—
— ^ 10
H-
*J
c 8
OJ
o
"o
—
—
— •
c I
-S 1
.t: 6&-
•o
o
0
4
2
„
\
•\
^f
Nf
-!A
r
_
f
1
\
\\
-
1
Gamma
10
1
Distribution
20 30. „
1 1
1 1
—
with C = 2.56, B = 0,063
Exponential Distribution with B = 0
t
16 — — —
—
—
—
_
—
—
\
»
s
\
S
t
\
^
\
*
1
i
^^
m
9
str-l
—
—
r i i
0.24
0.20
O
o
0.16 §.
0>
0>
—
0.12 0
-------
Table 2. CONDITIONAL STATISTICS OF RAINSTORM
PARAMETERS FOR DIFFERENT DURATIONS
(a) td = 1 hr
Parameter
t , hr
td » hr
in.
D
mm
in./hr
i ,
mm/hr
in.
°d »
mm
t , hr
G , hr2
a , hr
b , hr
in.
h ,
mm
o
ad
t°
G°
o
a
b°
h°
Mean
55.6
1
0.07
1.78
0.07
1.78
0
0
0.50
0.33
0.50
0.50
0.14
3.56
0
0.50
0.33
0.50
0.50
2
Standard
Deviation
94.5
0
0.15
3.81
0.15
3.81
0
0
0
, 0
0
0
0.29
7.37
0
0
0
0
0
0
Min
1
1
0.01
0.25
0.01
0.25
0
0
0.50
0.33
0.50
0.50
0.02
0.51
0
0.50 '
0.33
0.50
0.50
2
Max
744
1
1.03
26.2
1.03
26.2
0
0
0.50
0.33
0.50
0.50
2.06
52.3
0
0.50
0.33
0.50
0.50
2
Number of Rainstorms = 176
32
-------
Table 2. (b) t = 2 hr
d
Parameter
*b ' hr
td ' hr
in.
D ,
mm
in./hr
i ,
mm/hr
in.
Q
nun
t , hr
G , hr2
a , hr
b , hr
in.
h
o
°d
t°
G°
o
a
b°
h°
Mean
67.8
2
0.24
6.10
0.12
3.05
0.07
1.78
0.92
1.18
0.76
1.24
0.24
6.10
0.23
0.46
0.29
0.38
0.62
1
Standard
Deviation
76.1
0
0.34
8.64
0.17
4.32
0.12
3.05
0.27
0.54
0.81
0.81
0.34
8.64
0.15
0.13
0.13
0.40
0.40
0
Min
1
2
0.02
0.51
0.01
0.25
0
0
0.53
0.40
-0.41
-0.43
0.02
0.51
0
0.27
0.10 -
-0.20
-0.22
1
Max
327
2
2.09
53.1
1.05
26.7
0.87
22.1
1.48
2.29
2.43
2.41
2.09
53.1
0.48
0.74
0.57
1.22
1.20
1
Number of Rainstorms = 134
33
-------
Table 2. (c) t . = 3 hr
d
Parameter
D
Q
D
i
°d
t
G
a
b
h
o
CTd
t°
G°
o
a
b°
h°
, hr
, hr
in.
' ram
in./hr
mm/hr
in.
9
mm
, hr
,hr2
, hr
, hr
in.
'• mm
Mean
54.3
3
0.36
9.14
0.12
3.05
0.09
2.29
1.38
2.58
1.14
1.86
0.24
6.10
0.23
0.46
0.29
0.38
0.62
0.67
Standard
Deviation
70.2
0
0.33
8.38
0.11
2.79
0.10
2.54
0.43
1.28
1.29
1.29
0.22
5.59
0.11
0.14
0.14
0.43
0.43
0
Min
1
3
0.03
0.76
0.01
0.25
0
0
0.56
0.48
-1.34
-1.18
0.02
0.51
0
0.19
0.05.
-0.45
-0.40
0.67
Max
304
3
1.44
36.6
0.48
12.2
0.53
13.5
2.40
5.97
4.18
4.34
0.96
24.4
0.45
0.80
0.66
1.40
1.45
0.67-
Number of Rainstorms = 64
34
-------
Table 2. (d) t = 4 and 5 hr
d
Parameter
tb , hr
td , hr
in.
D ,
mm
in. /hr
i
mm/hr
in.
°d '
mm
t , hr
G , hr2
a , hr
b , hr
in.
h
mm
o
t°
G°
o
a
b°
h°
Mean
55.7
4.32
0.65
16.5
0.15
3.81
0.15
3.81
1.97
5.17
1.58
2.74
0.30
7.62
0.20
0.46
0.28
0.37
0.63
0.47
Standard
Deviation
64.9
0.47
0.61
15.5
0.14
3.56
0.15
3.81
0.60
2.57
1.73
1.77
0.28
7.11
0.08
0.13
0.13
0.40
0.40
0.05
Min
1
4
0.08
2.03
0.02
0.51
0.01
0.25
0.66
0.75
-2.02
-1.30
0.04
1.02
0.06
0.17
0.05
-0.51
-0.33
0.4
Max
292
5
2.86
72.6
0.72
18.3
0.53
13.5
3.42
12.4
5.30
6.29
1.43
36.3
0.39
0.78
0.65
1.33
1.51
0.5
Number of Rainstorms = 50
35
-------
Table 2. (e) t, = 6-14 hr
d
Parameter
tb , hr
td ' hr
in.
' mm
in. /hr
i ,
mm/hr
in.
mm
t , hr
G , hr2
a , hr
b , hr
in.
' mm
o
t°
G°
o
a
b°
h°
Mean
61.5
8.55
1.35
34.3
0.16
4.06
0.15
3.81
3.83
22.0
2.93
5.61
0.32
8.13
0.12
0.44
0.27
0.32
0.68
0.25
Standard
Deviation
77.8
2.62
0.85
21.6
0.10
2.54
0.13
3.30
1.70
17.6
3.42
3.31
0.21
5.33
0.06
0.13
0.13
0.39
0.39
0.07
Min
1
6
0.30
7.62
0.05
1.27
0.03
0.76
1.00
2.54
-4.99
-2.10
0.10
2.54
0.04
0.13
0.05
-0.62
-0.30
0.14
Max
324
14
3.52
89.4
0.53
13.5
0.74
18.8
7.68
71.2
10.0
13.8
1.07
27.2
0.27
0.77
0.68
1.30
1.62
0.33
Number of Rainstorms
31
36
-------
The above procedure gives a set of fitted conditional density
functions for i corresponding to different values of t,. Although it is
not necessary from a practical viewpoint, it is often convenient for easy
application to obtain a single expression as the conditional density
function of i given t • i.e. f(i|t,). For the present example, this is
done by using the exponential and gamma density functions and by expressing
the parameters of these density functions as functions of the storm
duration, using the values fitted in Figs. 7a to 7e. For the gamma
distribution, as shown in Fig. 8, assuming an exponential relationship
between B or C and t , the fitted expressions are
B = 0.33td~°-79 (27)
and
C = 0.31t, (28)
d
Thus, from Eqs. 21, 27 and 28, the conditional gamma density function is
0.79, 0<31td .0.31t,-l , ... 0.79,
> x d
d
Likewise, for the exponential distribution, a plot of B against t, (Fig. 9)
yields
B = O.OSt,0'37 (30)
d
and from Eqs. .20 and 30, the conditional exponential density function is
.) = 12.5t, °'37 exp(-12.5it, °'37) (31)
d d d
37
-------
0.5
I I I I I I I I
0.4
•B = 0,33 t.
-0,79
0,3
CD
0.2
O.I
468
Rainstorm Duration, t., in hr
10
12
i i
i i i i
C 2
O
I
C = 0,3I t
466
Rainstorm Duration, td, inhr
10
12
Fig. 8. Gamma density function parameters as functions
of rainstorm duration
38
-------
0.20
I I I I I
CO
Rainstorm Duration, t., in hr
a
Fig. 9. Exponential density function parameter as function
of rainstorm duration
39
-------
The conditional probability analysis can be extended to obtain
joint density functions for practical uses. For instance,
f(i, td) = f(i|td) f(td) (32)
By using the exponential density functions, Eqs. 31 and 32 together with
Eq. 20 (with value of B = 2.45 for f(t.) given in Fig. 4) yield
f(i, t,) = 5.1tJ~°'37 exp(-12.5it ~°'37 - O.Alt,) (33)
d d d d
Similar analyses can be performed on other parameters, and the
procedure can be extended to trivariate cases, e.g., f(t|D, t,). However,
d
trivariate frequency analysis is rather tedius, requiring large amount of
data, and at present is unlikely to be undertaken by most engineers.
(H) Application procedure. - in urban storm runoff problems, two types
of application often arise in connection with the statistical analysis
just described. The first is for design of certain facilities such as
treatment plants and overflow devices. The second is for operational
purposes involving the prediction of the time of occurrence, depth and
duration of the next rainstorm after a rainstorm of a given depth and
duration has just occurred.
To illustrate the application to design, assume that the storm
f'
runoff quantity is the controlling factor in determining the capacity of
a waste water treatment plant and that the plant is to be designed with
a capacity which will be exceeded on the average at most twice a month.
For the 3-month summer rainstorm data over the 14-year period at Urbana,
the average number of rainstorms for a 3-month summer period is 455/14 =
32.5. Assuming that the Urbana data is applicable for the design under
consideration, the given exceedance frequency of twice a month (6 times
in 3 months) corresponds to an exceedance probability of 6/32.5 = 0.185
40
-------
and a non-exceedance probability of 0.815. It is obvious that among
rainstorms of the given frequency those with large depth of rainfall
and short duration are most critical to the design. Assuming that
depth is the most significant one among all the rainstorm parameters
considered, from the non-exceedance probability curve shown in Fig. lOa
that for the given design frequency the design rainstorm depth is
12.9 mm (0.51 in.). Subsequently, from the conditional probability
density function of t,, f(t, |D), the most frequently occurred
duration for this depth, i.e., the mode of f(t |D), can be found and
used as the design duration. For the Urbana summer data for D = 12.9 mm
(0.51 in.), based on the 49 rainstorms with D between 8.9 and 16.5 mm
(0.35 and 0.65 in.), the mode of t is 2 hr (Fig. lla). Likewise, the
mode of t can be found from the conditional probability density
function f(t|D), or alternatively, f(t|t ) or f(t|D, t ). Because t
appears to be more sensitive to t than to D, and because the trivariate
conditional probability density function is difficult to obtain, the
mode of f(t|t ) is adopted as the design value of t. For the example
Urbana data (Fig. lib) the value of t used for the design is 1.025 hr.
Substituting the values of D, t , and t into Eqs. 14, 15, and 16 yields
the design values of the triangular hyetograph shape factors a, b, and
h, being 1.075 hr, 0.925 hr, and 12.9 mm (0.51 in.), respectively. The
design hyetograph thus determined can subsequently be routed through
the surfaces and sewers of the drainage basin using the routing
methods that will be described later to give the design discharge or
hydrograph for the treatment plant.
For the case of application to storm water runoff control,
the problem is of the nature of flow prediction for management purposes.
For instance, for an existing drainage system, when a rainstorm comes,
41
-------
.0
o
2
a.
o>
o
0>
0)
o
X
UJ
o
z
Rainstorm Depth , D, in mm
50
I 2 3
Rainstorm Depth, D, in in. -
100
Fig. 10. Non-exceedance probabilities for rainstorm depth and
elapse time between rainstorms
(a) Rainstorm depth
42
-------
£ 0.4
u
x
UJ
I
g
Z 0.2
i r i i
I I I
°0 24 6 B 10
Elopse Time Between Roinstorms, f^, in hr
200 40O 6OO
Elapse Time Between Rainstorms, t^, in hr
800
Fig. 10. (b) Elapse time between rainstorms
43
-------
1 I \ I I I I
i i I I i r
Gamma Distribution With C = 3.00 , B = 1.08 ——
Exponential Distribution With B = 3.24
to
(£)
6
V
a
VI
"in
to
O
en
a>
O
a
c
o
O
o
0.4
0.3
0.2
o
c
a>
3
cr
O)
>
a>
o:
O.I
12
14
16
Rainstorm Duration, td, in hr
Fig. 11. Conditional distributions of rainstorm duration
and hyetograph first moment arm
(a) Rainstorm duration for 0.35" ^ D < 0.65"
44
-------
I I I I I I I I I I I I I
0.15
CM
ii 2
C
o>
O
o
C
o
~ I
TJ
C
O
o
0.10
o
c
0)
0)
w.
u_
o>
or
0.05
0.2 0.4 0.6 0.6 1.0 1.2 1.4
First Moment Arm Of Hyetogroph , T , in hr
1.6
Fig. 11. (b) First moment arm of hyetograph for t, = 2 hr and
for all D
A5
-------
it is desirable for operational purposes to know the time of occurrence,
depth and duration of the next rainstorm so that decision can be made
on the utilization of in-line storage and other control facilities.
Consider the case when a rainstorm of depth D and duration t_. has just
occurred. Data from Urbana and other locations indicate that the
elapse time between rainstorms is nearly independent of the depth and
duration of the preceding rainstorm. Therefore, assuming that t, is
b
independent of D and t , , the most probable time of occurrence of the
next rainstorm can be estimated as the mode of the probability density
function such as the one shown in Fig. 6. The depth of the next
rainstorm can be evaluated as the mode of the conditional distribution
of the depth of a rainstorm given the depth of the previous rainstorm,
f(D«|D ). Subsequently the duration of the next rainstorm can be
determined from the mode of the conditional probability density function
f(t |D ) and the shape of the hyetograph for the next rainstorm from
f(t|t,_) as described in the design application.
Of course, many refinements and improvements can be made on
the procedures for runoff control and for design just described. For
example, in the flow prediction for storm water runoff control, the
elapse time t, can be estimated by using a. non-exceedance probability
b
function for t. (Fig. lOb) with an assumed or selected non-exceedance
b
probability. The depth of the next rainstorm can be estimated by using
f(D2]D1, tdl) instead of fCDD.^, and t by f(t|td2> D2) instead of
f(t|t,2). The corresponding risks of the prediction can be evaluated
accordingly using the joint density functions of the parameters involved.
However, such refined methods are rather tedious and complicated, and
available data are often inadequate to establish the needed conditional
-------
probabilities. Therefore, they are suggested to be considered only in
the future after significant information based on simpler procedures
are obtained and when adequate data are available.
IV-2. Infiltration and Other Abstractions
Not all of the rainfall produces runoffs. In hydrology the
losses that do not produce surface runoff are called abstractions.
These losses consist of interception, evaporation, transpiration, in-
filtration, and depression storage. Interception is the amount of
rainwater being intercepted by trees, vegetation, posts and buildings
that never reach the ground surface. For urban areas the rainfall on
roofs which is drained to the surface or directly to the sewers is
not considered an interception. The relative importance of inter-
ception on runoff depends on the intensity and duration of rainfall.
In urban areas usually there is no dense woods or vegetation, the
amount of interception is no more than a fraction of an inch (a few mm)
and mostly occurs during the beginning of the rainfall. Therefore,
for relative heavy rainfall of short duration, which would be of
importance for design or pollution control because of overflow, the
amount of interception is less than a few percent of the runoff
volume and can be neglected without causing serious accuracy problems.
Evaporation and transpiration are often considered
simultaneously for obvious practical reasons. Evapotranspiration may
be important when the water balance over a long period is considered.
But as shown by Shen et al. (1974), it is negligible when heavy
rainfall over short duration is considered, particularly in view of
the vegetation and tree situation in most urban areas. Actually, there
is no difficulty to include interception and evapotranspiration losses
47
-------
when reliable formulas become available. At present the amount of these
losses is assumed negligible.
Infiltration loss is a major factor affecting surface runoff.
Infiltration is defined as the process of water flowing through the
ground surface, i.e., the interface between the fluid environment and
the soil environment below. Various theoretical approaches and empirical
formulas have been proposed to estimate infiltration. At the beginning
of this research project a study was made to use Philip's (1969) theory
to derive a four-parameter method to account for infiltration. Un-
fortunately, available field data are 'inadequate to substantiate this
approach and hence the simpler and popular Horton's formula is used. In
fact, even for Horton's formula which is a three-parameter function,
there are difficulties to establish the values of the parameters based
on existing data. Philip (1969) also proposed a two-parameter
approximate infiltration equation. However, this equation has not been
adequately tested nor has it been widely accepted by Civil Engineers.
Horton's formula is
f = fc + (f0 - fc) e"kt (34)
in which f is the instantenous infiltration capacity; f and f are the
initial and final infiltration capacities, respectively; t is time; and
k is an exponent accounting for the decay rate of infiltration. For
Horton's equation to apply, the water supply rate (rainfall and water
stored on land surface) must be equal to or greater than the infiltra-
tion capacity at that instant. Otherwise, the entire amount of water
is assumed infiltrated.
The difficulty in applying Horton's formula to actual drainage
basins arises partly from the fact that in experimental and field
48
-------
rainfall-runoff studies, the soil properties and surface moisture condi-
tions are often inadequately recorded, and partly because of a natural
drainage basin, the soil condition is inevitably nonhomogeneous and
the values of f , k and f are different for different areas within the
c o
drainage basin. The measured basin runoff hydrograph merely reflects the
integrated effects of infiltration and other factors, and there is
actually no single set of representative values of f , f , and k for
the entire basin.
Theoretically, strictly speaking, for a given soil none of
the values of k, f , and f is constant. They depend on the soil type,
fluid properties, moisture condition of the soil, and water pressure
(usually depth) on the ground surface. The initial infiltration
capacity, f , obviously depends heavily on the initial soil moisture
condition. The final infiltration rate, f , and the exponent
expressing the decay rate, k, are the soil properties and should be
constant if secondary effects such as those due to changes in soil
flow potential near the ground surface, in water depth, in fluid
properties and seasonal effects are neglected. Philip (1969, Figs. 16
and 17) has shown that for a given soil and liquid, f is essentially
constant but the decay rate decreases with decreasing initial moisture
and with increasing overlaying water depth, i.e., k decreases with
decreasing initial moisture content or increasing water depth. However,
for a given surface in the practical range of conditions the variation
of k is relatively small. In other words, from a purely theoretical
viewpoint considering a liquid entering a porous medium, the values
of f and k are not constant, but from a practical viewpoint in using
Horton's formula, the values of f and k can be treated as essentially
49
-------
constants. The approximate constancy of f for a given soil surface
has generally been accepted. The relative constancy of k is still a
matter of debate. Many researchers using experimental data tried to
show that k varies considerably with initial soil moisture condition
and other factors. However, a small error in measurement in f would
c
easily give a varying k for the same soil. Of course, the dif-
ferences may also reflect the seasonal effects. Unfortunately, for the
purpose of the present study on storm runoffs this uncertainty on
infiltration imposes a serious problem on the accuracy of the results,
making a reliable comparison of the runoff prediction methods difficult.
Should Horton's formula be used with f and k treated as constants, it
is suggested that different sets of values be used for different
seasons.
A simple one-parameter approach, the -index method has also
been used for rainfall-runoff studies. The method assumes a constant
infiltration rate over the period of rainfall. This method is com-
patible with the requirement for the more sophisticated urban runoff
methods evaluated in this study.
It should be noted here that not all the infiltrated water
is necessarily lost because some water may find its way through sub-
surface flow to contribute to the basin runoff. In urban basins this
may occur as infiltrated water entering sewers through joints and
other leakages. However, such a case is not considered in this study.
The amount of loss due to depression storage depends on how
the term is defined. Loosely it is usually defined as the water to
fill the ground depressions before surface runoff starts. The amount,
obviously, is a function of the surface texture. Actually, as it is
50
-------
commonly defined, the depression storage includes a thin layer of water
held by surface tension before surface runoff starts. Expectedly
the depression storage is of statistical nature. The most commonly
used formula for depression storage supply rate (in./hr or mm/hr) is
(Linsley et al. , 1949)
8 = (i - f) exp FlS-] (35)
c
in which SG is the depression storage capacity expressed in depth (in.
or mm); P is the cumulative rainfall in depth; and F is the cumulative
infiltration in depth.
IV- 3. Snow Melt
Runoff from snow melt in urban areas differs from that of rural
areas in two major aspects. First there is more heat available for snow
melting in urban areas than in rural. Second and most important, the
intense human activities and interferences in urban areas would hasten
the melting process. Although runoff from snow melt is never a problem
in urban storm sewer design because its magnitude is smaller than the
flash flood of urban runoff due to heavy rainstorms, it is of considerable
importance in pollution control because of the quality of the melted
water, particularly when de-icing additives are used to speed up the
melting process.
The energy needed for snow melting comes mainly from three sources:
the radiant heat from the sun, the conduction heat from the environment, and
the latent heat of vaporization released by the condensation of water vapor.
The first two are the major ones to be considered for urban snow melting.
Many empirical and semi-empirical studies have been made on snow melt in
51
-------
rural areas (Chow, 1964, Section 10). However, these studies consider a
time of melting usually much longer, than that which would be interesting
and useful for urban settings, and the interferences of human activities
are not included. Snow melt in urban areas is a topic practically
untouched. A possible approach is to consider the energy budget of snow
melting and the thermodynamic processes involved, such as the idea outlined
by Eagleson (1970, Chap. 13). But such an idea has not been extended or
developed into any form nearly adoptable in practice, and to undertake
such a research is beyond the scope of this study.
From an engineering viewpoint the worst condition in terms of both
runoff quantity and quality for snow melt is melting of snow under warm rain.
For this situation the following daily snow melt formula recommended by
the U.S. Army Corps of Engineers (1960) is tentatively adopted in this study:
M = 0.007 P (T - 32) (36)
r a
in which M is the daily snow melt in in., ?r is the daily rainfall in in.,
and T is the mean daily temperature of saturated air at 10-ft above
a
ground in °F. If M and P are in mm and T in °C, the formula becomes
M = 0.013 P T (37)
1C Si
52
-------
V. DRAINAGE SYSTEM CHARACTERISTICS OF OAKDALE AVENUE BASIN
The Oakdale Avenue Drainage Basin in Chicago was selected to
verify and evaluate the urban storm runoff models. It is one of the
very few drainage basins for which relatively compatible and reliable
data are available. The basin is located in a residential section in
the city of Chicago (Fig. 12). It is approximately 2 1/2 block long
2
by 1 block wide and has a drainage area of 0.052 km or 12.9 acres. The
basin consists entirely of residential dwellings, and the drainage char-
acteristics are relatively uniform. The street pattern may be considered
as being typical of many cities in Illinois and the United States.
V-l. Surface Drainage Pattern
The drainage pattern of the land surface of the Oakdale Basin
is shown in Fig. 13. There are 30 inlet catch basins, each delivering
its water to a sewer junction and each receiving water from a gutter
except inlets 7, 14, 15, 21, 29, and 30 (Fig. 13) which receive water
from two gutters for each of these 6 inlets. Each of the 36 gutters is
contributed by water from one or more sub catchments as shown by the
dashed lines in Fig. 13. The subcatchment area of the Oakdale Basin
consists of four types of surfaces: roofs, lawns, paved.sidewalks, and
street pavements. The relative percentages of size of these four types
of surfaces vary from subcatchment to subcatchment, and as one would
expect, also change with time as a result of change of land uses and
structures. Detailed data on the distribution of these four types of
surfaces for each of the subcatchments in the Oakdale Basin is not
available.
53
-------
Lake
Michigan
N
I
o
E
o
w
o
_J
z*
n
»
f
, . i
"n i
o>
•^--1
ZL
W.Wellington Ave.
L— ULF"1 r--
W. Oakdale Ave. |
o
0>
1 _g.
i W. George St. z
I
n ^
*
• — ' °
z
1
0
I-
250
5OO ft
50 100 ISOm
Fig. 12. Oakdale Basin location map
54
-------
m
r — F
— M<
1
ft
1
t
t '
1 Ifl
1
ri
t i n~«—
115
in
1
1
I9~ I3/
' .n
1 i
14 lOj,
I
n , |
1 '
9J. el
1 »
J
• 113
1110
l-n
ft! ! — '
L- .^
30" 28N
t 26([
! ,
1
L
(29 251 24T 23T . 22X
I : t i i i I •*
_J J_ _L _
U t i *
i J 1 * j 1 *
i
j
t
i
21 , 1ST , I7f
"•i ! 1 !
I i. _
t f
— — 1
L * J
o
a
—i—i
ij I Jl
^ U_ 'X j_
§108 .^±•104 §|Q3 IJI02 [•IOIB3IOO
^ » iH
UJ
i
UJ
3
0 80 Gutter Inlets
• Sewer Junctions
~*~ Surface Flow Direction
Scale
250 ft
75 m
Fig. 13. Drainage pattern of Oakdale Basin
-------
In most metropolitan areas, data on such detailed land surface
uses are generally nonexistent. For a small drainage basin like the
Oakdale Basin, it is possible to conduct a detailed survey to actually
measure the physical properties of each of the subcatchments. However,
even for such a small basin, at least several man-months of work is
needed to obtain the data. Atop of this difficulty there is always the
problem of getting permission to survey in private properties. One
simpler, faster, and less costly method is to use aero-photos or satellite
pictures. Because of the time limitation of the research project, this
photogrammetric method was not undertaken. Nevertheless, even if such
detailed surface use data is available, it would be extremely tedious
and costly to actually use it to route the rainwater through all
the surfaces to the gutters. In view of the practical consideration
of the costs involved in obtaining and using the detailed subcatchment
surface data and the seasonal changes of the surface characteristics, it
appears most unlikely that any practical urban storm runoff simulation
model would require to use such detailed subcatchment information. For
the case of the Oakdale Basin used in this research, although the instru-
mented survey was not extended to cover the details of the subcatchments,
a thorough visual survey was conducted in order to provide reliable
information for the evaluation of the urban storm runoff methods.
The streets in the Oakdale Basin are 8.5 m (28 ft) wide, paved
with asphalt having a cross-slope of 2.4 to 3.0% on both sides of the street
crown. The longitudinal slope of the streets varies as listed in Table 3.
The street slopes are typically flat as for most of the Midwest cities.
Between Leclaire Street and the outlet of the basin Oakdale Avenue is
actually sloping down towards west although the sewer line underneath
has its slope in the opposite direction.
56
-------
The street gutters are all triangular in cross section formed
by cast-in-place concrete. The curbs are essentially vertical and the
curb height varies along the gutters and of-ten interrupted by driveways.
The most frequently observed curb height is 0.25 m (10 in.). The lateral
angle, 9, between the gutter bottom and the vertical is 1.54 radians. There
is a break of lateral slope where the gutter bottom joins the street
pavement. Consequently, assuming that the water surface of the gutter
flow is horizontal along the lateral direction, the relationship between
the flow cross sectional area, A, and depth measured from the apex of the
gutter, h, (Fig. 14) is
v,2
A = — for h <; W cote (38a)
2 cote
2
A = Wh - W—£0 for h > W C0t9 (38b)
in which W is the gutter width. The hydraulic radius R is
sin6 (39a)
uya;
2TT+ cose)
in, w2 cote
Wh - 7.
R = L - (39b)
h H
sin6
and the water surface width b is
b = h/cote for h <, W cot6 (40a)
b = W for h > W cote (40b)
57
-------
Fig. 14. Schematic drawing of gutter cross section
58
-------
The longitudinal slopes of the gutters are the same as those for the streets.
The gutter width measured horizontally is 0.3 m (1 ft). The geometric
dimensions of the 36 gutters are listed in Table 3. In most of the time the
gutters are kept reasonably clean although some debris have been observed.
Accordingly Manning's roughness factor n is estimated to be 0.013 for the
gutters.
A certain amount of rainstorm water is discharged directly from
subcatchments into the alleys between the streets as shown in Fig. 13.
Hydraulically these alleys act like wide shallow channels to transport the
water into inlets or gutters. Most of the alleys have concrete surface with
uneven joints and cracks and their estimated Manning's roughness factor n
is 0.016. The length, width, and slope of the alleys are listed in Table 4.
The 30 inlets in the basin are grate inlets either circular or
rectangular in shape as listed in Table 3. The details for the circular
grate inlets are shown in Fig. 15a and those for the rectangular in Fig. 15b.
The approximate locations of the inlets are identified by the inlet numbers
in Fig. 13. The distances between the inlets are given in Table 3. Some
of the inlets do not start from the curb line but offset slightly and
extend beyond the gutter proper into the street pavement. Such irregularity
occurs mostly for replacing inlets with clogged inlet catch basins.
Apparently some of the inlet catch basins have the clogging problem. There
is no record to identify whether the inlet catch basins surveyed now in
1973-74 are the same as those a decade ago.
V-2. Sewer System
The combined sewer system of the Oakdale Avenue Drainage Basin
consists of 18 circular sewer pipes and 18 junctions or manholes plus the
sewer system outlet. The diameter of the concrete pipes ranges from 0.25 m
59
-------
SECTION B-B
SECTION A-A
SECTION D-D
y_Jfe222
SECTION E-E
--
HT
N.
— L -
SECTION F-F
SECTION H-H
SECTION K-K
SECTION L-L
Fig. 15. Details of grate inlets
(a) Circular grate inlets
60
-------
E 1
D
1
0
D
—
1
I
:
—
i
3
!
]
D
I
J
PLAN
SECTION E-E
7Ti£j
i/r I/T I/T i/r
ir'-l
SECTION D-D
INLET GRATE
iok Casflron
'*:'8'Wo-fer dound Macadam
i/T^-' or &'Concrete <5
ASSEMBLY
INLET CASTINGS
Fig. 15. (b) Rectangular grate inlets
61
-------
(10 in.) to 0.76 m (30 in.) and the pipe roughness is estimated to be 0.01
ft or 3 mm. The junctions or manholes are marked as 3-digit numbers (e.g.,
101 to 118) in Fig. 13 and the dimensions of the sewers are listed in
Table 5. Most of the manholes are of flow- through type with a half-cut
pipe embedded at the manhole bottom connecting the upstream and downstream
sewers to induce smooth flow at low discharge. None of the junctions or
2
manholes has a horizontal cross sectional area bigger than 2 m (20 sq ft)
and hence their storage capacity is relatively small.
For circular sewers such as those used in the Oakdale Basin, when
the pipe is flowing partially filled with a depth h and central angle
6 = cos [1 - (2h/D)] as shown in Fig. 16, the flow area A and the
corresponding hydraulic radius R and water surface width b are
A = 5- (6 - sinG) (41a)
o
R = (1 _ sine.) (41b)
(41c)
siny
for 0 ^ 6 ^ 2ir and D is the pipe diameter.
At the outlet of the Oakdale Avenue Drainage Basin a Simplex 0.76 m
(30 in.) Type "S" parabolic flume is placed in a vault at the corner of
Oakdale and Lamon Avenues to measure and record the basin runoff. This runoff
measurement and recording system together with a tipping bucket recording
rain gauge located at one block north of the basin has been in operation
since 1959 measuring rainfalls and runoffs. Details of these measuring
devices and the data collected can be found elsewhere (Tucker, 1968) and are
not presented here.
62
-------
Fig. 16. Circular sewer flow cross section
63
-------
Table 3. DIMENSIONS OF GUTTERS OF OAKDALE AVENUE DRAINAGE
BASIN
Gutter
From To
; Inlet Inlet
.
: - 1
! 1 2
; 2 3
3 7
4
1
i - 5
! 4 7
i 5 6
! 6 8
! 8 9
i 9 10
i 10 14
1 - 11
: - 12
; 11 13
i
i 12 14
i 13 15
! - 15
I - 16
i 16 17
i 17 18
! 18 21
! - 19
1
i - 20
i 19 21
; 20 22
' 22 23
23 24
• 24 25
: 25 29
j - 26
': ~ 27
' 26 28
27 29
; 28 30
30
Gutter |
Length
ft m
58 18
201 61
192 59
120 37
104 32
104 32
•
46 14
42 13
117 36
194 59
200 61
114 35
100 30
100 30
192 59
192 59
96 29
96 29
58 18
201 61
192 59
120 37
100 30
100 30
42 13
42 13
117 36
194 59
200 61
114 35
151 46
151 46 .
128 39
128 39
96 29
96 29
Size of
Contributing
Sub catchments
2 '
ac m
0.17 690
0.53 2140
0.51 2060
0.32 1290
0.04 ., 160
(0.54)* (2190)*
0.04 160
(0.64)* (2590)*
0.03 120
0.02 80
0.32 1290
0.51 2060
0.53 2140
0.31 1250
0.02 80
0.02 80
0.07 M 280
(0.50)* (2020)*
0.07 280
(0.68)* (2750)*
0.26 1050
0.26 1050
0.17 690
0.53 2140
0.51 2060
0.32 1290
0.01 ., 400
(0.-54r (2190)*
0.01 400
(0.78)* (3160)*
0.02 80
0.02 80
0.32 1290
0.51 2060
0.53 2140
0.31 1250
0.13 530
(0.64)* (2590)*
0.13 „ 530
(0.41)* (1660)*
0.05 200
0.05 200
0.26 1050
0.26 1050
; Type of Grate
] Inlet at Down-
Longitudinal stream End of
Slope Gutter*
0.0012 C
0.0012 1 R
0.0012 i C
0.0012 1 R
0.0010 C
i
0.0010 i C
0.0010 R
0.0010 i R
0.0027 ; R
0.0027 • R
0.0027 i C
0.0027 R
0.0010 C
o.ooio r c
0.0010 i R
i
0.0010 R
i
0.0010 C
0.0010 i C
0.0012 ' C
0.0012 R
0.0012 i C
0.0012 i C
0.0010 i C
(
0.0010 I C
;
0.0010 ; R
o.ooio i „ c
0.0027 : ' C
0.0027 ! R
0.0027 • C
0.0027 -R
0.0010 C
0.0010 C
0.0010 i R
o.ooio ; R
0.0010 C
0.0010 i C
Contribution
to Sewer
Junction
102
104
106
109
107
107
109
109
110
112
114
117
115
115
117
117
118
118
102
104
106
109
108
108
109
109
110
112
114
117
116
116
117
117
118
118
Type of grate inlets: C = circular, R = rectangular
^'Contribution by alleys
Gutter width W = 1.0 ft (0.30m)
Gutter bottom inclined 88° 15' from vertical curb
Manning's n for gutters = 0.013
64
-------
Table 4. DIMENSIONS OF ALLEYS OF OAKDALE AVENUE
DRAINAGE BASIN
Lo
Alleys between
Wellington and
Oakdale
Alleys between
Oakdale and
George
nation
West of Leclaire
East of Leclaire
West of Leclaire
East of Lavergne
West of Leclaire
East of Leclaire
West of Lavergne
East of Lavergne
Length
ft m
395 120.4
295 89.9
335 102.1
295 89.9
396 120.7
236 71.9
380 115.8
290 88.4
Width
ft m
15.5 4.7
15.5 4.7
15.5 4.7
15.5 4.7
15.5 4.7
15.5 4.7
15.5 4.7
15.5 4.7
Slope
0.0047
0.0049
0.0042
0.0053
0.0043
0.0053
0.0040
0.0054
Contributing
to inlet
13
14
5
4
26
27
20
19
Table 5. DIMENSIONS OF SEWERS OF OAKDALE AVENUE
DRAINAGE BASIN
' Sewer
From
•Node
118
: 115
116
. 117
i 114
• 113
112
111
110
107
108
109
106
105
104
103
102
101
To
Node
117
117
117
114
113
112
111
110
109
109
109
106
105
104
103
102
101
Length
ft
108
170
105
134
34
168
158
38
131
50
45
153
39
156
156
61
73
100 32
m
32.9
51.8
32.0
40.8
10.4
51.2
48.2
11.6
39.9
15.2
13.7
46.6
11.9
47.6
47.6
18.6
22.3
9.8
Slope
0.72
0.71
1.08
0.45
0.45
0.45
0.40
0.40
0.40
3.78
4.20
0.35
0.35
0.35
0.30
0.30
0.30
0.30
Diameter
ft
1.00
0.83
; 0.83
1.25
1.25
-1.25
1.50
1.50
1.50
0.83
0.83
; 1.75
1.75
• 1.75
2.00
2.00
2.00
2.50
m
0.30
0.25
0.25
0.38
0.38
0.38
0.46
0.46
0.46
0.25
0.25
0.53
0.53
0.53
0.61
0.61
0.61
0.76
65
-------
Vi. ILLINOIS SURFACE RUNOFF MODEL
The Illinois '."rban St.orm Runoff method actually consists of two
parts: the surface runoff model and the sewer system routing model. The
input into the surface runoff model is the hyetograph and the output is the
inlet hydrographs which constitute the input into the sewer system routing
model. The sewer routing model is the Illinois Storm Sewer System Simulation
Model (Sevuk et. al. , 1973) and will be described briefly in the following
chapter. The Illinois surface runoff model is a recent development and will
be discussed in this chapter. It should be noted here that improvement and re-
finements are continuously being made on both surface and sewer models and those
reported here are the most up-to-date versions at the time of writing this report.
VI-1. Runoff in Subcatchments
The surface runoff is subdivided into two subsequent parts: the
subcatchments which consists of only strips of overland flows receiving
rainfall as the input; and the gutters which receive water from the sub-
catchments as well as from direct rainfall and deliver the water into inlet
catch basins to produce inlet hydrographs. The overland surface of a
drainage basin can be approximated by a number of equal-width rectangular
strips of different lengths. A large number of such strips of narrow
width will closely approximate the actual overland surface. But this
will require a large amount of computations without significant improve-
ment in accuracy. Contrarily, to'o few strips would approximate the
actual geometry poorly.
Time varying free-surface flow including overland flows can be
described mathematically by a pair of partial differential equations called
the St. Venant equations (Chow, 1959; Yen, 1973a, 1973b; Sevuk et al., 1973)
66
-------
<«>
^ TT *\ TT 'M, 1
o V d V on 1 i
- + V- + g cosO — = g(So-Sf) + ^ (UrV)q da (43)
in which x is the direction of the flow measured along the bed; t is time;
A is the flow cross sectional area; b is the width of the free surface;
D = A/b is the hydraulic depth; V is the cross sectional average flow velocity;
h is the depth of the flow above the invert; 6 is the angle between the channel
bed and the horizontal; S = sinQ is the bed slope; Sf is the friction slope;
a is the perimeter bounding A; q is the lateral discharge per unit length of a
having a velocity component U along the x-direction when joining or leaving
the flow; and g is the gravitational acceleration. The first equation is the
equation of continuity and the second the momentum equation.
With the appropriate initial and boundary conditions, these
two equations can be used to solve numerically using digital computer
for the overland flow on sub catchments. However, the solution requires
considerable amount of computer time and in view of the usually large
number of sub catchments (overland strips) such an approach is feasible
but impractical. In field conditions, the accuracy of data on overland
geometry and rainfall input usually does not render the accuracy that
the St. Venant equations can provide. Several approximations of Eqs. 42
and 43 are possible. Often the overland flow is approximated by using
the Manning's formula or the Izzard's method which essentailly assumes
the flow or the rainfall to be steady. Other approximations of the
St. Venant equations include the kinematic wave model and diffusion
wave model (Yen, 1973a). From past experience the non-linear kinematic
wave model was found to be most suitable for solving overland flows
67
-------
because it does not require a downstream boundary condition and hence
considerably reduces the computer time and difficulties and yet its accuracy
is substantially better than that given by Manning's formula. Those who are
interested in the relative accuracy of the different approximate models can
refer elsewhere (Sevuk, 1973, Yen, 1973a) .
For the kinematic wave approximation, the inertia and pressure terms
of the momentum equation (Eq. 43) are neglected; thus,
SQ = Sf (44)
The friction slope S can be estimated by using the Darcy-Weisbach formula
2
in which f is the Weisbach resistance coefficient given by the Moody diagram
and R is the hydraulic radius, by the Manning formula
Sf = 2^2 V2 R-A/3 (46)
where n is the Manning's roughness factor, or by the Chezy formula
V2
sf = x
in which c is the Chezy factor. The continuity equation (Eq. 42) and Eq. 44,
together with the initial condition and one upstream boundary condition,
can be solved numerically for the unsteady flow. Equation 46 is for V in
fps and R in ft; if V is in m/sec and R in m, the coefficient is unity
instead of 2.22.
In selecting the resistance formula to approximate the friction
slope, the Weisbach coefficient has the advantage of being dimensionless and
having better theoretical justification, whereas Manning's n has the advantage
68
-------
of being nearly constant independent of flow depth for flows over rough
boundaries with sufficiently high Reynolds number. However, for overland
flows the depth is usually so shallow and the Reynolds number of the flow not
sufficiently high that it would be erroneous to consider n to be constant
(Chen and Chow, 1968; Yen, 1975). Therefore, the Weisbach formula
(Eq. 45) is adopted in this model to evaluate the overland flow.
In Eq. 45, the value of f is given by the Moody diagram which
can be found in standard hydraulics reference books (e.g., Rouse, 1950;
Chow, 1959). For the case of overland flow under rainfall, limited in-
formation was given by Yen et al. (1972) and Shen and Li (1973). Based on
the available information, the Weisbach f is computed as
f=| (48)
for laminar flow, in which H = VR/v is the Reynolds number of the flow
where v is the kinematic viscosity; and the coefficient C is
C = 24 + 101 i°'4 (49a)
for i in mm/hr, or
C = 24 + 27 i°'A (49b)
for i in in./hr. Since the surface of natural overland is inevitably
rough, for turbulent flow, f is constant as
1 9R
— = 2 log f^+ 1.74 (50)
/f k
where k is a length measure of surface roughness. The transition between
Eqs. 48 and 50 is shown schematically in Fig. 17. The critical Reynolds
number ]R determining which equation should be used is
69
-------
VJ
O
o»
O
2 Log •*£- + 1.74
Log R
Fig. 17. Evaluation of Weisbach resistance coefficient
-------
Kc = C (2 log ^+ 1.74)2 (51)
When the Reynolds number of the flow E. < TR , Eq. 48 applies. Otherwise,
Eq. 50 is used. It should be noted that actually there is a transition
between laminar (Eq. 48) and fully developed rough turbulent flow (Eq. 50).
This transition is neglected here and the steady uniform flow values of f
are used for unsteady cases. This approximation can of course be improved
when more information on f becomes available.
The water input onto the subcatchment surface to produce the
flow is the lateral flow q in Eqs. 42 and 43. The value of q is equal to
the rainfall minus infiltration. The infiltration is estimated by using
Morton's formula (Eq. 34). When the rainfall rate is smaller than the
infiltration capacity, the deficiency is supplemented by the water on the
surface, if any.
In solving Eqs. 42, 44 and 45 numerically, the initial condition
to start the solution cannot be zero depth and zero velocity because this
condition will impose a mathematical singularity. In reality, when rain
falls on a dry overland surface, there is indeed an initial wetting
process before runoff starts. The surface tension will hold a
small amount of water without producing runoff. Therefore, the initial
condition for the overland runoff from the subcatchments can be assumed
as a small finite depth with zero velocity. In other words, immediately
following the commencement of rainfall, after infiltration and other
.losses are subtracted, the water left on the overland surface simply"
accumulated without producing runoff until the initial depth is reached.
This initial depth depends on the slope and nature of the overland surface.
Future studies will provide more information on this initial depth. It
suffices at present to assume the initial depth to be 0.0012 in. or
0.03 mm. It has been found that the final solution is practically unaffected
71
-------
by the value of the initial depth so long as it is assumed within a
reasonable realistic range.
Several numerical schemes can be used to obtain the solution
(Sevuk and Yen, 1973). A 4-point, noncentral, semi-implicit scheme is used
to solve the equations because of its independent selection of the time and
space increments (At and Ax) in the computations without stability problems
and consequently saves computer time.
A more accurate and convenient form of Eq. 42 for the purpose of
numerical solution is
where Q is the flow rate at any flow section, A is the flow cross-sectional
area and q = qda is the lateral flow per unit length of flow in
'0
x-direction, being positive for inflow. Applying the chain rule of
differentiation to Eq. 52, and letting G(h) = 9Q/9h and b(h) = 9A/9h, one
obtains
9h
For the semi-implicit four-point backward difference scheme adopted, re-
ferring to the fixed rectangular grid in Fig. 18, the coefficients and
partial differential terms of Eq. 53 may be approximated by the "following
expressions
G = (GD + Gc) (54a)
b = (bD + bc) (54b)
f = fe (hc - V (55a)
f =(h+h-h" (55b)
72
-------
A
0)
E
U)
2 i-l
Space,x
B
i + l
Fig. 18. Computational grid for semi-implicit four-point
backward difference scheme
-------
Substitution of Eqs. 54 and 55 into Eq. 53 yields
(GC + V (hC - V + 4A7 (bD + V % + hC - hA - V = % (56)
The flow parameters at grid points A and B are known either from the initial
conditions or from previous time step computations, and the flow parameters
at point D are known either from the upstream boundary condition or from
previous computations. Therefore, with Gr and b_ being specified functions
C* L>
of h , the only unknown in Eq. 56 is h .
C 0
Surface runoff on subcatchments usually occurs in the form of
open-channel flow in wide channels. Consequently the runoff problem can
be simplified by solving for the discharge per unit overland width, Q .
Hence, in Eq. 56, b. = b,, = b_, = b = 1. For laminar flow in a wide
A JJ L. 1)
rectangular channel, combining Eqs. 44, 45 and noting that Q = VR and
R = h, one obtains
8gSo 3
Hence
8Q 24gS
r - u - Q u
G ~ 3h~ ~ ~CT~ h
For the case of turbulent flow, Eqs. 45 and 50 yield
Q = v/8is~ (2 log ^ + 1.74) h372 " (59)
U O K
and
G = /8gS^ (3 log - + 3.47) h (60)
O K.
74
-------
The depth of water h at any subcatchment flow section is obtained by
using Eqs. 56 and 58 for laminar flow, and Eqs. 56 and 60' for turbulent
flow. Newton's iteration technique is used for the numerical processes.
Knowing the flow depth, the discharge per unit width is evaluated by using
Eq. 57 and Eq. 59 respectively, for laminar and turbulent flows. After
the flow parameters at grid point C is computed, the computations for the
next downstream station at the same time level (grid point E in Fig. 18)
can be performed. After the flow parameters at all the stations at a given
time level are evaluated, the computations is advanced to the next time
level starting from the upstream end.
VI-2. Gutter Flow Routing
Street gutters and surface runoff in defined channels receive water
from overland runoff of subcatchments, from upstream water sources, if any,
and directly from rainfall. The input water is transported through the gutter
or channel into the inlet catch basin to produce the inlet hydrographs for
the sewer runoffs. Horton's formula (Eq. 34) is assumed applicable to account
for infiltration. Theoretically, the gutter flow is also described mathe-
matically by the St. Venant equations (Eqs. 42 and 43). Again, using these
equations to solve for gutter flows is feasible but impractical in view of the
large number of gutters for a drainage basin and the accuracy, of the input
data. In field conditions gutters are rarely prismatic channels because of
poor control in construction and interruption of local facilities such
as driveways and other intersections. Furthermore, in actual operation,
gutters are often obstructed by debris, parked cars and the like. Such
obstructions are time varying and random in nature. It is possible to describe
the b.oundary condition .of the gutters precisely. for any particular ..runoff. .. .
75
-------
considered. Therefore, solving the gutter runoff by using the St. Venant
equations, which require large amount of computer time and detailed geometry
data, cannot be justified. Consequently the kinematic wave approximation is
adopted for gutter flow routing as for the case of overland flow in
sub catchments .
The differential equations of the kinematic-wave model for gutter
flows are the same as those for overland flows, i.e., Eqs. 42, and 44.
Therefore, the same solution technique can be used for both cases, and
Eq. 56 is the finite difference equation representing also the gutter flow.
However, the flow conditions in gutters are often within the range where
Manning's formula (Eq. 46) can be used to approximate the friction slope,
S . Use of Manning's formula instead of Darcy-Weisbach' s (Eq. 45)
simplifies the computation as it is no longer needed to check the in-
stantaneous flow Reynolds number in order to estimate Sf. Thus, from
Eqs. 44 and 46,
0 = i s 1/2 AR2/3 (61)
n o
and
(f AT^ f * ^ f >
where C = 1 in SI system and 1.49 in English system. The terms A, R,
3A/3h and 3R/8h in Eq. 62 should be evaluated from the cross-sectional
shape of the gutter under consideration. For instance, for the
triangular gutters in the Oakdale Avenue Drainage Basin these terms can
be evaluated by using Eqs. 38 and 39. Specifying b in Eq. 56 as a
function of h from the gutter geometry, it is possible to solve Eqs. 56
and 62 simultaneously for the flow depth h at grid point C by use of
L*
Newton's iteration technique. The discharge then is evaluated from
76
-------
Eq. 61. The computation progresses in the downstream direction for each
time level as described for subcatchment runoffs.
The initial condition for the gutter flow routing is essentially
the same as that for runoff in subcatchments. As to the boundary
conditions, kinematic wave model requires only the upstream conditions
be provided. In the present gutter routing model the upstream boundary
condition is provided as specified flow depths at the upstream end of the
gutter at each time level. These flow depths are evaluated within the
model using Manning's formula (Eq. 46) and they correspond to the carry-over
of water from the upstream inlets. When there is no such carry-over at the
upstream end of a gutter, the depth of water at the upstream flow section
is assumed to be always equal to the initial depth.
Vl-3. Inlets
Inlets are one of the most important components of urban drainage
systems to determine the time distribution of urban storm runoffs. They
control the amount of water to flow from gutters into sewers. The hydraulic
characteristics of an inlet depend on the geometric properties of the inlet.
Unfortunately, despite the large number of inlets used in streets and highways,
the geometries of inlets have never been standardized. Furthermore,
in their operation, inlets are seldom kept clean to be free from foreign
materials partially clogging the inlet.
Inlets can be classified as curb type and grate type. A combina-
tion of the two is also used. Hydraulically they can be described by the
weir formula
Q = CdbH3/2 (63)
or the orifice formula
Q = CdAH (64)
77
-------
in which C, is a discharge coefficient; A is the cross-sectional area of the
orifice opening; b is the length of the weir; and H is the available head.
The value of H depends on the gutter or surface flow depth near the inlet.
The difficulty in using Eqs. 63 and 64 for inlet flow computation is the
wide range of variations of the values of C, and A, particularly for unclean
inlets. Also, the determination of the range of application of the weir and
orifice formulas is a matter of debate.
The inlet imposes a backwater effect on the gutter flow. For a
supercritical flow in the gutter, disturbance waves cannot propagate upstream
and hence the numerical solution of the St. Venant equations or its nonlinear
kinematic-wave approximation can proceed forward from upstream without de-
pending on the downstream boundary conditions.
For a subcritical flow the gutter flow is directly affected
by the hydraulic conditions at its downstream end, i.e., the inlet. Conse-
quently the inlet flow condition, which is by itself unknown and yet to be
solved, becomes the necessary downstream boundary condition for the numerical
solution of the St. Venant equations. Contrarily, for the kinematic-wave
approximation, no downstream boundary condition is required. Consequently,
the gutter flow can be solved without requiring simultaneous solution of the
yet unknown inlet flow conditions, and hence the solution technique can be
simplified and the required computer time greatly reduced. The inlet flow
can subsequently be computed as will be described later. Such an approxi-
mation neglecting the backwater effect due to the inlet of course differs
from the reality. However, in view of the uncertainties on the physical
conditions of the gutters and inlets, it appears to be justified from a
practical viewpoint that the gutter flow is computed by using the nonlinear
kinematic-wave approximation and the inlet"f low is "computed independently.
78
-------
A more accurate and potentially practical approach to gutter-inlet
flow solution is to use generalized nondimensional curves describing inlet
runoff hydrographs for different input and geometry conditions (Akan, 1973).
However, this approach requires at least a certain degree of standardization
tion of the inlets in order to avoid a large number of nondimensional graphs
and hence it is not adopted here, although it may be used in the future for
the refinement of the surface runoff model.
In the Illinois surface runoff model, the average depth plus the
velocity head at the end of the gutter is used as the value of H in Eqs. 63
and 64 for the calculation of the inlet discharge. The discharge coefficient
C, in these equations is assigned different values according to the type of
inlet under consideration. For instance, in Eq. 63, C, is assumed to be
equal to 3.0 for grate inlets with longitudinal bars and for combined in-
lets, 2.4 for grate inlets with diagonal bars, 2.7 for grate inlets with
cross bars, and 1.2 for curb openings. The corresponding C, values in
Eq. 64 are 0.60, 0.48, 0.54, and 0.30, respectively. The inlet discharge
is first computed by using Eq. 63 until this equation gives a discharge
greater than the discharge of the approaching gutter flow. From then on,
it is assumed that the flow around the inlet has the characteristics of
orifice flow, and the inlet discharge is computed by using Eq. 64. During
the recession of the gutter runoff, when the computed inlet discharge
using the weir formula is smaller than the approaching gutter flow, the
inlet discharge is assumed to be computed again by using the weir formula,
Eq. 63.
When the approaching gutter flow is greater than the inlet
discharge, the excessive water is assumed carried over the inlet to continue
on as the input flow into the next gutter immediately following. Should
79
-------
there exist more than one downstream gutter such as at an intersection,
the model allows a distribution of the carry-over flow among these
downstream gutters. However, the distribution factors should be provided
on the program data cards. This carry-over of excessive flow from inlet
is assumed to continue until the flow reaches a low point such as
Junctions 109 and 117 in Fig. 13 where no further carry-over can reason-
ably be assumed and a reservoir storage routing is performed for the
discharge through the last inlet and the storage around it. :
The assumptions on the distribution of carry-over flow, on the
transition between the weir and orifice flows, and on the values of C.,
d
are not precise as the reality. Improvement and refinements on these
aspects can be made in the future when more reliable and useful laboratory
and field data become available.
80
-------
VI-4. Program Descjciy_ition_and_ Data Preparation
The Illinois Surface Runoff Model is programmed in Fortran IV
language for computer solutions. The input into the computer program
consists of the geometric characteristics of subcatchments, gutters,
inlets and the identification of the sewers joining the inlet catch
basins, and also the rainfall hyetographs. The output is the inlet
hydrographs which serve as the input into the sewer system. The program
also performs water quality computations of the runoff to produce inlet
pollutographs. The formulation and details of the water quality model
will be given in Chapter VIII.
(A) Program Description. - The computer program of the Illinois
Surface Runoff Model as listed in Appendix B allows the consideration of
a maximum number of 100 gutters at a time. Along each gutter, the
subcatchments can be approximated by as many as 10 rectangular strips.
These strips of overland areas may have different lengths, slopes, and
surface and infiltration properties, but should be equal in width. Two
different pollutants are considered at a time for each gutter and
subcatchment strip. The program allows for the entire basin a maximum
of five zones of rainfall with different hyetographs. The computational
logic is shown schematically in Fig. 19. The computer storage requirement
for the program in its present form is about 400K. If more st-orage is
available, the program can easily be modified to consider larger basins.
This modification can be achieved by simply changing the arrays in the
dimension statements.
The computer program consists of one main program and six sub-
routines. The relationship between the main program and the subroutines
is shown schematically in Fig. 20. A brief description is as follows:
MAIN PROGRAM: It reads and stores data for the entire basin. It performs
81
-------
CONSIDER THE MOST
UPSTREAM GUTTER
SET FLOW CONDITION
EQUAL TO INITIAL CONDITION
CONSIDER NEXT GUTTER
TIME = TIME + AT
IS THE TIME
GREATER THAN ESTIMATED
INLET HYDROGRAPH
DURATION?
HAVE ALL THE
GUTTERS IN THE SYSTEM
EEN CONSIDERED?
COMPUTE UPSTREAM
BOUNDARY CONDITION
(SUBROUTINE UPSBO)
PERFORM STORAGE
ROUTING COMPUTATIONS
(SUBROUTINE STROUT)
IS THE
GUTTER IMAGINARY?
COMPUTE DIRECT INFLOW
HYDROGRAPHS AND
POLLUTOGRAPHS FOR
EACH OF THE SEWER NODES
(SUBROUTINE SWRJNT)
i
r
PRINT OUT AND PUNCH
RESULTS ON CARDS
COMPUTE QUANTITY AND
QUALITY OF LATERAL
INFLOW FROM SUBCATCHMENTS
(SUBROUTINE OVLFLO)
COMPUTE DIRECT LATERAL
INFLOW FROM RAINFALL
PLUS SNOWMELT MINUS
INFILTRATION
(SUBROUTINE RASNIN)
COMPUTE DEPTH AND RATE
OF FLOW AT EACH STATION
ALONG THE GUTTER
STOP
COMPUTE WATER QUALITY
OF GUTTER FLOW
COMPUTE FLOW INTO INLET
AND CARRY-OVER FROM
INLET AT THE END
OF GUTTER
(SUBROUTINE DOWBOU)
Fig. 19. Flow chart for Illinois surface runoff model computer program
82
-------
eC
O£
O
O
C£
Q-
2!
1— i
=a:
•\
/
UJ
C g
t— <
fe
O
Ol
X f/> V
^ _, J 7
V -^
X X.
RASNIN
OVLFLO
UPSBO
DOWBOU
STROUT
SWRJNT
Fig. 20. Composition of Illinois surface runoff
model computer program
83
-------
the gutter routing computations based on nonlinear kinematic wave and
Manning's equations. Newton's iteration technique is used to solve
Eqs. 56 and 62, and discharge is computed by using Eq. 61. The
computations are made starting from the most upstream gutters and
proceeding towards downstream. A mechanism is built in the program to
decide which gutters should be considered first. This allows the
gutters in the drainage basin to be numbered arbitrarily from 1 to 100
while preparing data. However, the user should specify the flow direction
in each gutter. The water quality computations for gutter flows are also
performed in the main program.
SUBROUTINE OVLFLO: Flow in subcatchment strips is computed in this
subroutine. Newton's iteration technique is used to solve the nonlinear
kinematic wave equations with Darcy-Weisbach's formula; i.e. Eqs. 56 to
61. An approximate form of the Moody diagram (Fig. 17) is built in the
subroutine to estimate the resistance coefficient. This subroutine is
called from the main program while the computations are being done for
each grid point along the gutter, unless the strip characteristics are
identical for which the preceding values can be used. The water quality
computations for each subcatchment strip are also performed in this
subroutine as will be described in Chapter VIII. The output from sub-
routine OVLFLO provides part of the lateral inflow for the gutter routing
in the main program.
SUBROUTINE RASNIN: -This subroutine computes from rainfall the rate of
lateral inflow for the subcatchments and part of the lateral inflow for
the gutters. The inflow is evaluated as rainfall plus snowmelt minus
infiltration. This subroutine is called from the main program and
subroutine OVLFLO.
84
-------
SUBROUTINE UPSBO: This subroutine computes the upstream boundary
condition for each gutter. It is called by the main program when com-
putations are made for the upstream end of each gutter at each time
level. The carry-over water from all the immediately upstream inlets
are summed up and a corresponding flow depth is computed using Manning's
formula. This computed depth is accepted as the upstream flow depth
for the gutter being examined at the time level being considered.
SUBROUTINE DOWBOU: Knowing the gutter outflow as computed by the main
program at each time level, this subroutine is called upon to evaluate
the inflow into the inlets and the carry-overs. The flow into the
inlets is computed by using the weir or orifice formulas (Eqs. 63 and
64) as explained in Sec. VI-3.
SUBROUTINE STROUT: This subroutine provides a storage routing procedure
around the inlets where there exist no immediate downstream gutters and
hence there is no carry-over. The inflow to the storage area consists
of the outflow from the upstream gutters. The outflow from the storage
area is the flow into the inlet computed by using the orifice formula
(Eq. 64).
SUBROUTINE SWRJNT: This subroutine prepares the output (inlet hydrographs)
in data cards and provides the link between the Illinois Surface Runoff
Model and the Illinois Storm Sewer System Simulation Model. After all
the flow hydrographs into the inlets are computed, this subroutine is
called by the main program. For the sake of convenience in linking
the surface runoff and sewer routing models, the output hydrographs of
the surface runoff model are identified by the sewer nodes, i.e., sewer
manholes or junctions, instead of the corresponding inlets if they
carry different identification numbers. When there are more than one
gutter inlets discharging into the same catch basin or sewer node, the
ordinates of those inlet hydrographs are summed up. The computed
85
-------
hydrographs for each of the sewer nodes are provided on computer cards
as a part of the surface runoff program output. These cards are in a
format compatible to the input data card requirements for the Illinois
Storm Sewer System Simulation Model and can be used as part of the data
deck for sewer routing. The ordinates of the sewer node inflow
hydrographs and pollutographs are also printed out from this subroutine.
(B) Data Preparation. - Detailed information on basin characteristics
is needed for the Illinois Surface Runoff Model and hence a number of
data cards are required for the model. The data deck consists of the
following sets in the order of presentation:
(1) General description of drainage basin: This set consists of two
cards. The first card of the first set in the data deck specifies
whether the data is provided in English or metric system of units. If
the English system is used, the integer number 1 should be punched in
the first column of the card. If the metric system is used, the integer
number 2 should be punched in the first column. The second card of the
set specifies the following information: the total number of gutters in
the system; the total number of sewer nodes in the system; the total
number of rain-zones considered; an integer number that indicates the
frequency of printed output (e.g. when the number is equal to 1, the
output is printed out at every time level); the time interval of com-
putation in min; time in min when the execution should stop for each
gutter corresponding to the estimated duration of the gutter outflow
2 2
hydrographs; the gravitational acceleration = 32.2 ft/sec or 9.81 m/sec ;
the average daily temperature on the day of rainstorm (in F or C); the
2 2
kinematic viscosity of water in gutter (in ft /sec or mm /sec); and the
constant C in Eq. 48. The first four quantities must be punched in 15
formats and the remaining six quantities must be punched in 15
format. The time interval and the stop execution time values should be
86
-------
selected such that the latter must be an exact multiple of the former,
and there should be no more than 100 time steps of computation. When no
snowmelt is involved the space for the daily temperature is left blank.
(2) Hyetographs: This set of data consists of several subsets of cards.
Each subset corresponds to a rainfall zone. The first data card in
each subset gives time in min at which the rainstorm starts (F10.0);
time in min the rainstorm stops (F10.0); the total daily rainfall on
the day of rainstorm in in. or mm (F10.0); and an integer number to
specify the number pairs of time and rainfall intensity values used to
describe the hyetograph (15). The other cards following in the same
subset give the hyetograph ordinates. There should be no more than 8
pairs of values on each card corresponding respectively to time in min
(F5.0) and rainfall in in./hr or mm/hr (F5.0). There should be no more
than 100 pairs of time-intensity values to describe a hyetograph.
(3) Gutter, Inlet and subcatchment descriptions: This set contains a
number of subsets. Each subset corresponds to a gutter considered in
the system. Each gutter is given the same number as the inlet at its
downstream end. When there exists no inlet, an imaginary inlet should be
assigned. The types of gutters and inlets considered in the program are
represented by a number as shown in Fig. 21. When an imaginary inlet
(type =0.0) is introduced, the flow into the inlet is always zero, and
the gutter outflow is equal to the carry-over from the imaginary inlet.
When storage routing is required for an inlet, it is necessary to -use
the concept of imaginary gutter (type = 0.0). When there is an imaginary
gutter, the program does not perform any gutter routing but it calls the
subroutine STROUT to perform storage routing. The data cards required
for each subset contains the following four or more cards: (a) The
first card of the subset specifies the gutter number (15); number of grid
points to be considered for gutter routing (15) ; rain-zone number to
87
-------
Gutter Type = 1.0
Gutter Type = I.O
v/yyyyyyyiwvw'''
Gutter Depth
Gutter Type = 2X)
Gutter Type = 2.0
(a) Gutter Types
I /A i
/
/.
f f f
/ /. _ ./
/ V V V N x/ /- - /• -t
A x x x x y / / /
A\\\\7 /— / /-
' f
i
Inlet Type =1.0 Inlet Type=2.0 Inlet Type=3.0 Inlet Type =4.0 Inlet Type =5.0
Imaginary Inlet , Inlet Type = 0.0
(b) Inlet Types
Type = 2.0
Type = 2.0
Type = 2.0
Type = 2.0
(c) Sewer Node Types
Fig. 21. Identification of types of gutter and inlet
88
-------
which the gutter belongs (15); type of gutter (F5.0): length of
gutter in ft or m (F5.0); width of gutter in ft or m (F5.0);
longitudinal slope of gutter (F5.0); depth of gutter in ft or m for
rectangular gutters (F5.0); the angle between the gutter plane and the
vertical in radians (F5.0); inlet type at the end of the gutter (F5.0);
width of the inlet in ft or m (F5.0); length of the inlet in ft or m
(F5.0); ratio of total area of the openings to the total area of the inlet
(F5.0); width of street pavement measured from the crown to the gutter in
ft or m (F5.0); the uniform initial depth of flow along the gutter in ft
or m (F5.0); and Manning's roughness factor for the gutter (F5.0). (b) On
the second card of the subset are the initial infiltration capacity of
gutter surface in in./hr or mm/hr (F5.0); the final infiltration capacity
of gutter surface in in./hr or mm/hr (F5.0); Horton's constant of decay
rate of infiltration, k, for gutter surface in hr~ (F5.0); initial
infiltration capacity for street pavement in in./hr or mm/hr (F5.0);
final infiltration capacity for street pavement in in./hr or mm/hr (F5.0);
Horton's constant of decay rate of infiltration, k, for street pavement
in hr (F5.0); initial concentration of the first pollutant associated
with gutter flow in ppm or mg/1 (F5.0); and initial concentration of the
second pollutant associated with gutter flow in ppm or mg/1 (F5.0).
(c) On the third card of the subset are the numbers of six immediately
upstream inlets (615) ; and the proportions of carry-over from these
inlets that go into the gutter (6F5.0). Note that when there are no
upstream inlets, this card should still be there but with nothing
punched on. (d) On the fourth card of the subset are the length of a
subcatchment strip in ft or m (F5.0); the slope of the strip (F5.0);
surface roughness in ft or m (F5.0); initial capacity of infiltration in
in./hr or mm/hr (F5.0); final infiltration capacity in in./hr or mm/hr
(F5.0); Horton's k for infiltration in hr (F5.0); uniform initial
89
-------
depth of flow along the strip in ft or m (F5.0); initial concentration
of the first pollutant associated with the subcatchment flow in the
strip in ppm or mg/1 (F5.0); initial concentration of the second
pollutant in ppm or mg/1 (F5.0); and number of computation grid points along
the strip (15). This card (d) should be repeated for. each of the
subcatchment strips starting from the one at the upstream end of the
gutter. The number of subcatchment strips is equal to the number of
computation grid points along the gutter minus one.
When considering an imaginary gutter, the area of the storage
2 2
surface is punched as the fifth quantity on card (a) in ft or m
instead of the gutter length. The other gutter properties can be assigned
any values since they will not be used. The number of grid points
should be assigned the value 2. Then the second (c) and the fourth (d)
cards each can be replaced by a blank data card.
(4) Sewer node description: In this set of data one card is needed to
describe each sewer node. There are two different types of sewer nodes
to be considered. The type 1.0 represents the junctions of sewers in the
layout. The upstream nodes without any incoming sewer pipes are classified
as type 2.0. On each card of this data set the following information is
*)
required: the sewer node number (15); the type of sewer^node (F5.0);
3
the base flow for the sewer node in cfs or m /sec (F5.0); the concentra-
tion of the first pollutant associated with the base flow in ppm or fflg/1
(F5.0); the concentration of the second pollutant in ppm or mg/1 (F5.0)
and the inlet identification numbers of up to ten gutter inlets discharging
into the sewer node under consideration (1015). When there are less than
ten gutter inlets discharging into the sewer node, the excess space
should be left blank.
A conceptual simple drainage system and the corresponding data
representation is shown in Fig. 22 as an example.
90
-------
I O
m
8
1-
I o
i in
f
8
B
| M_
75'
__j p
111!
.s
BB
15
i S = 0.01
>-+ m
LB
|*
•
L..
i 1
|
it
i
50' 50'
i
"in
CVJ
50'
n
7 5 j 8
K
ro
30'
^n^
ii1
in
Sewer Junctions
Identified By 15 And 16
Inlets 6 And 8 Are
Imaginary
Gutter 5 Is Imaginary
Subcatchments I
f o « LO in./hr; fe • 0.3 in./hr,
k = 4.0 , S • 0.06 ,
Surface Roughness B QOI ft
Gutters And Street Pavements !
fo=O.O5inyhr; fc«O.OI in./hr,
k = 15.0, S = O.OI, n *OjOI3
Grate Inlets '.
Width = 1.0 ft,
Length = 2.0ft
Opening Ratio =0.60
Sewer Line
2' 14'
14' 2'
1.5 rod. 15 rod?1
Section AA
14' 2'
Section BB
I, in7hr
i
Rainfall Hyetograph
(A Simple Zone)
Fig. 22. Example of Illinois surface runoff model data preparation
(a) Layout of conceptual basin
91
-------
JUJ. J. j-
• 0
SET 2 ! o.S*
60.0
2
0.05
1
75.
0.
1
0.05
75.
75.
3
0.05
i . 90.
| 90.
: 90.
1 90.
4
0.05
1
SET 3 200'
50.
50.
7
0.05
75.
75.
100.
: 6
; 0.05
0.
0.
8
0.05
37.5
; 75.0
75.0
5
4
15
SET 4 ,6
2
0.0
0.0
3
0.01
0.06
0.
3
0.01
0.06
0.06
5
0.01
0.06
0.06
0.06
0.06
4
0.01
O.C6
0.06
0.06
4
0.01
0.06
0.06
0.06
3
0.01
0.
0.
4
0.01
0.06
0.06
0.06
2
7
2."
1.
1
60.0
5.0
1
15.0
0.01
0.
1
15.0
0.01
0.01
1
15.0
0.01
0.01
0.01
0.01
1
15.0
0.01
0.01
0.01
1
15.0
0.01
0.01
0.01
1
15.0
0.
0.
1
15.0
0.01
0.01
0.01
1
2
0.2
0.2
2
3.0
1.0
0.05
1.0
0.
1.0
0.05
1.0
1.0
1.0
0.05
1.0
1.0
1 .0
1.0
2.0
0.05
1.0
1.0
1.0
2.0
0.05
1.0
1 .0
1.0
2.0
0.05
0.
0.
2.0
0.05
1.0
1.0
1.0
0.0
3
0.5
0.5
1
3.5
LO.O
100.
o.oi
0.3
0.
100.
0.01
0.3
0.3
200.
0.01
0.3
0.3
0.3
0.3
150.
0.01
0.3
0.3
0.3
150'.
0.01
0.3
0.3
0.3
VO.
0.01
0.
0.
VO.
o.oi
0.3
0.3
0.3
324,
6
0.1
0.1
.0
3.0
2.0
15.0
4.0
0.
2.0
15.0
4.0
4.0
2.0
15.0
4.0
4.0
4.0
4.0
2.0
15.0
4.0
4.0
4.0
2.0
15.0
4.0
4.0
4.0
2.0
15.0
0.
0.
2.0
15.0
4.0
4.0
4.0
8
1
2
90.
9
2Q.O
0.01
0.5
1.0
.0001
0.
0.01
0.5
.0001
• oooi
o.oi
0.5
.0001
.oooi
.0001
• oooi
n.ol
0.5
.0001
.0001
.0001
o.oi
0.5
.0001
• oooi
.oooi
o.oi
0.5
0.
0.
0.01
0.5
.0001
.0001
.oooi
1.0
3
0
3.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
1.
0
0
1
5
0
1
5
5
0
1
5
5
5
5
1
1
5
5
5
1
1
5
5
5
1
1
1
1
5
5
5
0
4
30
1.
0.
0.
1.
0.
0.
1.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
1.
32
.0
5
1
5
1
1
5
1
1
1
1
1
1
1
1
1
1.
1
1
1
0
6
.2 34.0 0.000011 24.0
3.0 40.0 3.0 50.0 3.0 55.0 3.0
1.0 1.0 2.0 0.6 14.0.00010.013
3
1.0 1.0 2.0 0.6 14.0.00010.013
3
3
1.0 1.0 2.0 0.6 14.0.00010.013
3
3
3
3
1.0 1.0 2.0 0.6 14.0.00010.013
6
3
3
1.0 1.0 2.0 0.6 14.0.00010.013
3
3
4
0.0 0. 0.0 0. 14.0.00010.013
0.0 0. 0. 0. 1«. 0.00010.013
3
3'
3
1.0 1.0 2.0 0.6
1.0 1.0 1.0
7 8
Fig. 22. (b) Data representation
92
-------
VII. SEWER SYSTEM ROUTING MODEL
As mentioned previously, the sewer routing model of the Illinois
Urban Storm Runoff Method is the Illinois Storm Sewer System Simulation
Model (ISS Model). The ISS Model actually consists of two options: the
design option for which the size of the sewers are to be determined, and
the flow prediction option for which the size of the sewers and junctions
are known and the objective is to compute the runoff hydrographs for
given inputs. Since the ISS Model haa been reported in detail
elsewhere (Sevuk et al., 1973) , and the objective of this research is to
investigate methods of prediction for urban storm runoffs for the
purposes of pollution control and management, only the flow prediction
option of the ISS Model is briefly summarized in this chapter. Those
interested in design of storm sewer networks are recommended to refer
to a comparative study on using the ISS and other methods (Yen and
Sevuk, 1975).
VII-1. Sewer Network Representation
One of the most important aspects in solving sewer flow
problems in a network of sewers is to properly and systematically
represent the geometric sequences of the sewers. This is particularly
important if computer solution is used for which a logical means of
selecting the proper order of sewers is a prerequisite of solution. For
a small network consisting of a few sewers, it is not difficult to
assign specifically the sequence that the computation should follow.
For a large system consisting of many sewers, and particularly with the
possibility of alternation of the sewer connecting pattern, it is more
desirable and practical to set up some rules that the computer can
follow to select the sequences. A computer program written for
93
-------
specific sequence can be used only for those networks with patterns
following that sequence, and it is necessary to number the sewers
precisely as the sequence requires. Contrarily, a computer program
written for a certain sequence selection rule allows arbitrary patterns.
In view of the variability of urban sewer systems and the large number
of sewers involved in each of the system, the latter approach of setting
a rule for sequencing to allow arbitrary numbering of the sewers is
adopted in the ISS Model.
In this approach the members of a sewer system are represented
by a node-link system commonly used in network analyses. The nodes are
the junctions or manholes that join the sewers and each is assigned
arbitrarily a number, as shown in Fig. 13 and Table 5 for the Oakdale
Avenue Drainage Basin. The sewers are the links in the network and they
are represented by two numbers, the first being the number of the
upstream node and the second the downstream node. The outlet of sewer
system is the root node and is joined by only one sewer. Thus, two
sewers that join at the same node each will have one of its two
identification numbers identical to the others. The sequence becomes
a systematic search of marked numbers, and the computer can easily
determine the connectivity, i.e. , the pattern, of the network.
In solving for the flow in a sewer network, the solution
obviously proceeds from upstream sewers toward downstream, no matter if
the backwater effect is accounted for or not. To determine which .
upstream sewers should be solved first and the sequence of the sewers
to be solved, a systematic searching method is adopted. The search
starts from the root node, i.e., the outlet of the sewer system, to
detect if the sewer connected to this node has already been solved.
If this connected sewer is not yet solved, then the search moves to its
94
-------
upstream node and the process is repeated again. For nodes joining two
or more branches, the process can be proceeded one by one following
certain order, e.g., following the relative order of the branches stored
in the computer, or following the order, say, from left to right of
the branches connected to the node. For example, if the latter rule of
left to right is used, the search will start from the root node. Since
the flow for the last sewer connecting to the root node has not yet
been solved, the search will move to the upstream node of the last
sewer. Assuming that there are two other sewers joining to this node,
the search will first look if the left one was solved. If it has
not been solved, the search will automatically move up to the
upstream node of this left sewer and repeat the entire process again.
If the left sewer has already been solved, the search will then move
to the sewer at the right. If this right sewer has not yet been solved,
the search will move to its upstream node. If the right sewer has
already been solved, this implies that the sewer with its upstream end
joining the node is the only one to be solved. After the solution for
this sewer is obtained, the search returns to its downstream end node.
In this manner, the solution is obtained systematically, branch by
branch, from upstream towards downstream.
VII-2. Method of Solution
In the ISS Model, the flow in each of the sewers is determined
by solving the St. Venant equations (Eqs. 42 and 43). The friction
slope, Sf, is evaluated by using the Darcy-Weisbach formula (Eq. 45).
The Weisbach resistance coefficient f is estimated by using a simplified
form of the Moody diagram. Since laminar flow rarely occurs in sewers,
only turbulent flow is considered. For fully developed turbulent flow
95
-------
in hydraulically rough conduits, Eq. 50 applies. For hydraulically
smooth conduits, the Blasius formula is
for H < 4 x 10 . In sewers hydraulically smooth boundary flow seldom
occurs with H > 10 . Hence if the Blasius formula is assumed to apply
up to slightly higher Reynolds number and the transition between smooth-
and rough-surface flows is neglected, the threshold Reynolds number, H*, is
pTj Q
P* = 0.633 (log ^ + 0.87) (66)
The value of f is computed from Eq. 65 or Eq. 50 depending on whether
]R is less or greater thanlR*.
One initial condition and two boundary conditions are needed
to solve the St. Venant equations. For supercritical flow, the two
boundary conditions are furnished by the flow conditions at the
upstream node of the sewer. This imposes no computational problem
since the solution is proceeding towards downstream. However for
subcritical flow, one boundary condition is furnished from the upstream
node and the other should be from the downstream node. This downstream
boundary condition physically represents the backwater effect from
the junction to the sewer and usually it is an unknown to be solved.
The junction flow condition, in turn, is determined by not only its
physical properties but also the flow conditions of all the sewers
joining to the junction. Therefore, to solve for the flow in a sewer
network, it is necessary either to solve simultaneously all the
equations describing mathematically the flows in all the sewers and
junctions of the network, or to subdivide the network into components
for solution by successive approximations. The first approach is
96
-------
possible and practical for small systems consisting of a few sewers and
junctions. For large systems this simultaneous solution method would
easily become immanageable. Therefore the second approach is adopted
for the ISS Model using a technique called overlapping Y-segments.
In the ISS Model the sewer systems are considered as a tree
type network, each consisting of branches formed by a number of connected
Y-segments. Each Y-segment contains three sewers joined by a common
junction. The hydraulic condition of a junction is accounted'for by a
dynamic equation in addition to the continuity equation commonly used.
The continuity equation is
<>i + Q2 + *, - Q3 - iT (67)
in which s is the storage in the junction; Q. is the direct inflow into
the junction; and the subscripts 1 and 2 represent the inflow sewers
and 3 the outflow sewer from the junction. If the storage of the
junction is negligible, the right-hand side of Eq. 67 is equal to zero.
The dynamic equation for a junction with large storage, i.e.,
reservoir type junction, is
V 2
Zl + hl = Z2 + h2 = Z3 + h3 + 2i~ (68>
in which h is the depth of sewer flow at the junction and z is the
elevation of the sewer invert above a reference horizontal datum. For
a point-type junction with negligible storage, the velocity head term
2
(V,. /2g) in Eq. 68 is assumed equal to zero. Furthermore, if the
inflowing sewer has a drop producing a free-fall of the flow, the depth
for that sewer is equal to the critical flow depth corresponding to the
instantaneous discharge of the sewer.
97
-------
The St. Venant equations (Eqs. 42 and 43) and the junction
equations (Eqs. 67 and 68) can be solved simultaneously for a Y-segment
with known upstream boundary conditions of the inflowing sewers and
assumed (if unknown) downstream boundary condition for the outflowing
sewer. Since the downstream boundary condition is assumed, the solution
is only approximate and a successive overlapping Y-segment technique
is adopted to improve the accuracy of the solution. The technique is
shown schematically in Fig. 23. Numerical solutions are first obtained,
branch by branch, for those Y-segments whose inflowing sewers are
connected to the inlet catch basins. The prescribed inlet hydrographs
(or the outflow hydrographs from the Illinois Surface Runoff Model) and
the compatibility conditions at the junction are used as the boundary
conditions. In addition, if the downstream condition of the Y-segment
is unknown, the forward differences are used as a substitution for the
downstream boundary conditon. The solution is obtained by applying
the St. Venant equations to each of the three sewers and Eqs. 67 and
68 to the junction and solving these equations simultaneously for the
Y-segment using a first-order characteristic method (Sevuk et al., 1973).
After the computation is completed, the first trial solution for the
outflowing sewer is discarded but the "true" solution for the inflowing
sewers is retained. Thus the inflow hydrograph into the junction of
the current Y-segment is obtained. This junction will serve as one of
the two inlets of the next Y-segment and the original outflowing sewer
will become an inflowing sewer for the advanced new Y-segment. This
procedure is repeated until the entire network is solved. For the
last segment of the system, the prescribed boundary condition at its
downstream end, the outlet of the system, is used and thus the numerical
solution over the entire network is completed.
98
-------
b.
a. Complete solution domain
b b
b. Solution domains for first-
order sewers (1 throuqh 5)
c. Reduced solution domain
d. Solution domain for second-
order sewers ( .6 and 7)
e. Final solution domain for
the sewers 8, 9., and 10
Fig. 23. Solution by method of overlapping y-segments
99
-------
The ISS Model, in its present form, can simulate the flow with
regulating and operational devices if they are located at the system
outlet. From the hydraulics viewpoint, such control facilities can be
expressed mathematically as stage-time relationship h = f(t); velocity-
depth relationship V = f(h); discharge-depth relationship Q = f(h); and
discharge-time relationship Q = f(t).
Similar to the case of solving numerically the overland and
gutter flows, the initial condition for the sewer routing cannot be
dry-bed; i.e., zero .depth and zero velocity will impose a computational
singularity. For combined sewers, the initial condition can be
evaluated from dry-weather flow. For storm sewers without initial base
flow, a small and negligible base flow is assumed to start the
computation.
In using the overlapping Y-segment technique, only three
sewers are considered at a junction. For a junction with more than three
joining sewers, only three can be considered for direct backwater
effects. Others, preferably those with small backwater effects from
the junction, can be treated as direct inflows, i.e., as Q. in Eq. 67.
For a junction joined by only two sewers, the third sewer of the
Y-segment can be considered as imaginary with zero length.
VII-3. Computer Program Description
The ISS Model just described has been programmed for computer
solution. A macro flow chart showing the logic of solution for the
flow prediction option of the computer program is given in Fig. 24. The
program begins its execution by reading a set of control and data cards
which describes the run control specifications (user commands), sewer
system layout and physical characteristics of system components and
100
-------
READ RUN CONTROL
SPECIFICATIONS
READ SYSTEM LAYOUT AND
PHYSICAL CHARACTERISTICS
OF SYSTEM COMPONENTS
TES
READ INFLOW HYDROGRAPHS
FIND Y-SEGMENT WITH
UPSTREAM SEWER(S)
EMANATING FROM INLET
CATCH BASIN(S)
TREAT JUNCTION OF
COMPUTED SEGMENT
AS INLET OF PRUNED
NETWORK
*
TRIM UPSTREAM
SEWERS OF COM-
PUTED Y-SEGMENT
IS
THERE
ANOTHER SEWER
SYSTEM
IS
NO/THIS LAST
SEGMENT OF
YSTEM
COMPUTE INITIAL
CONDITIONS
COMPUTE TIME INTERVAL
AND ADVANCE TO NEXT
TIME STEP
COMPUTE FLOW CONDITIONS
AT UPSTREAM BOUNDARY
STATIONS
COMPUTE FLOW CONDITIONS
AT INTERIOR STATIONS
COMPUTE FLOW CONDITIONS
AT DOWNSTREAM BOUNDARY
STATION
COMPUTE FLOW CONDITIONS
AT JUNCTION STATIONS
HAS
FLOW AT
ENTRANCE OF
DOWNSTREAM SEWER
APPROACHED BAS
FLOW
Fig. 24. Flow chart for ISS model computer program flow
^prediction option
101
-------
inflow hydrographs. After verifying the accuracy and completeness of
this information, the program starts to execute the numerical
simulation phase. Output from the numerical simulation consists of a
tabular print out of the computed or specified sewer diameters, time
variations of flow rate, velocity, and depth at sewer entrances, and
space variations of flow rate, velocity, and depth along the sewers
at specific time intervals. In addition to this default print out
(the amount of which is under the control of user through various
control card options) , optional data capture facilities are provided
which allow the user to produce his own plots or graphic displays
of flow behavior at the inlets, junctions, and outlets of the network.
After completion of the flow simulation for one sewer network, the
program proceeds on to the next network, if any, until all such
networks are processed.
At its present form, the computer program of the ISS Model
cannot account for moving hydraulic jumps and surges within the sewers
and only circular sewers can be considered. The latter restriction
may be removed, since for sewers with noncircular cross sections,
simulation can still be made by using equivalent diameters of the
conduits based on best fit of hydraulic radius or cross sectional area.
The computer program of the ISS Model consists of around 3,000
statements written in PL/1-source code, as well as an assembler
language subroutine. The PL/1 portion includes several short external
subroutines, a main controlling section with network control card input
routines, and a large set of numerical computation routines. The
program is written to be composed of modules, whenever possible, so that
most of the routines are separately compiled and then linkage edited
together with the assembler language routine to form the executable load
module.
102
-------
The program can be adapted for execution on most large IBM
360 or 370 models running under OS/360 (preferably OS/MVT) operating
system or its VS equivalent. A memory storage of 300K bytes is
desirable although under certain conditions 100K bytes may be
acceptable. The machine must have the floating point instruction
set. In implementation, the program resides on a direct-access storage
device, such as a disk, a drum, or a data cell drive. This enables
the program to be loaded quickly and efficiently into the main storage
without entailing the added overhead of recompilation each time the
program is executed.
Despite the sophistications of the theory and programming
techniques of the ISS Model, the computer program was written for easy
adoption by users having only elementary computer knowledge. Pro-
gramming experiences in PL/1 language will be helpful but not essential
for the implementation of the program. The program can be implemented
through the use of either a distribution tape or the program listing.
The former approach is recommended, particularly for those users who
are not technically oriented on computer operations. Prospective
users may obtain the distribution tapes at cost of the magnetic tape,
handling, postage, and computer time for duplication.
Standard distribution tapes are provided on 1600 bytes per inch,
2400 ft magnetic tapes. When requesting for a distribution tape the
following information should be provided:
.(a) Model of computer system, e.g., IBM System 370, model
168 MP.
(b) Operating system in .use, e.g., OS/MVT with HASP.
(c) Region size available for the program, e.g., 220K.
103
-------
(d) Type of disk or other direct-access device to which the
program will be transferred from the distribution tape,
e.g. , 2314, 2319, 3330.
(e) Compilers available at the installation, e.g., either,
neither or both of the OS PL/1 (F) and PL/1 Optimizing
compilers.
However, users who prefer to use the program listing for
implementation or modification should have adequate understanding of the
following IBM systems 360 and 370 concepts:
(a) creation, use, and modification of partitioned data sets
for storage of source and object decks, as well as load
modules;
(b) use of the linkage editor to produce executable load
modules, and method of executing the load modules thus
created;
(c) use of other IBM system utilities, such as IEHMOVE,
IEBUPDTE, etc;
(d) use of tapes, disks, and other direct access devices; and
(e) use of IBM job control language (commonly referred to as
JCL) .
Those who need to use the program listing and yet unfamiliar with computer
applications are advised to seek assistance from experienced programmers.
The listing of the computer program, input and output data format,
program structures, operational procedure and the finite differences equa-
tions used for numerical solution have been reported in details by Sevuk
et al. (1973) and are not repeated here for brevity. Those interested in
details should refer to that report.
104
-------
VIII. WATER QUALITY MODEL
As mentioned in Chapter III, Introduction, storm runoff has been
known as a significant source of pollution because of its capability to
wash and carry pollutants on its course of flow. With recent stringent
requirements on water quality for the control of pollution, the necessity,
desirability, and economical feasibility of treatment of storm runoff is
often an unavoidable concern. For combined sewer systems, the problem is
further complicated by the variable quality as well as the quantity of
the dry weather flow.
For pollution control purposes, ideally the time and spatial
distributions of the quality of runoff should be known, at least at
certain key locations like sewer outlets and overflows. This knowledge
is particularly useful, for example, for a selective withdraw and treat-
ment of storm runoff. However, such detailed information requires, from
an experimental viewpoint, the time-consuming, tedious and expensive
measurements and laborous analyses of the data, and from a theoretical
viewpoint, the precise theories concerning diffusion and dispersion of
pollutants and the chemical and thermal processes involved. Furthermore,
a prerequisite for an accurate water quality analysis is a reliable
quantitative prediction of water.
The water quality model introduced in this chapter is an attempt
to provide a means to evaluate water quality of storm runoff as a
supplement to the quantitative evaluations of the Illinois Surface Run-
off Model and the ISS Model. The model is a relatively simple one by
using a one-dimensional approach considering the time and spatial varia-
tions of the pollutant expressed in terms of cross-sectional averaged
concentration. Only the equation of conservation of mass is used.
105
-------
Implicitly this involves the assumption of instant mixing within the
flow cross section. The dynamic equation expressing the gravity effect
(buoyancy or settling) is not used in the model. Neither are the effects
of biological and chemical reactions accounted for. Improvements on
these aspects, of course, can and should be made in the future. However,
this is to be done after the simple model has been adequately tested
with data. At present available field measurements are mostly for limited
number of locations in a drainage basin and do not provide enough data on
spatial and temporal variations of storm runoff water quality to verify
the model sufficiently.
VIII-1. Water Quality Model Formulation
The water quality model described here utilizes the computed
discharge hydrographs and flow area-(or depth-) time relationship at
desired locations of the drainage system from the surface runoff and ISS
Models together with the initial pollutant distribution to evaluate the
transport of pollutant by storm runoff. The output is a set of polluto-
graphs (pollutant concentration-time graphs) at the desired locations.
The mass conservation equation of a pollutant expressed in
concentration c for the flow with a discharge Q and flow cross sectional
area A
<69>
in which c., is the concentration for the lateral flow q.. The term
c0qfl includes concentrated pollutant dosage in which case the value of
Xf JC
c^q,, is simply replaced by the value of the dosage. Substitution of the
continuity equation (Eq. 52) into Eq. 69 yields
rc> (70)
106
-------
By using the computational grid shown in Fig. 18, Eq. 70 can be written
in finite difference form as
Ac + 'c = «* < Vc>
However, because of the details of the output of the Illinois
Surface Runoff Model and the ISS Model, the computational procedure for
water quality for these two parts of runoff are different. For the
surface runoff part, the water quality computation is programmed in the
Illinois Surface Runoff Model and the quantity and quality computations
are done concurrently. In the sewer system runoff part, the water
quality model is programmed separately using the output from the ISS
Model as the input.
(A) Water quality computation for surface runoff. - Water quality for
surface runoff is computed within the Illinois Surface Runoff Model. In
Eq. 71, if the computational grid points D and C represent respectively
the upstream and downstream ends of a gutter or subcatchment strip, the
flow parameters for D at each time level are. known from the upstream
boundary conditions. The values of A and Q are supplied from the
u u
water quantity computation of the model. For gutter flows, the values
of c and q are known from the subcatchment runoff. For subcatchment
Xj JG *""
flows, the values of c and q are known from rainfall (usually with
Jo Jo
c = 0) and other known lateral flows, if any. Hence, from Eq. 71 the
concentration c can be solved explicitly for each time level as
o
, CBAC CDQC
q£°£ At Ax ,-,.,.
CC = Ac Qc (72)
AT + Ax" + q£
(B) Water quality model for sewer system runoff. - In view of water
quality routing, a sewer system can be considered to consist of two
107
-------
different types of elements; namely, the sewer conduits and the junctions
with substantial storage capacities. For a sewer with no lateral flow,
Eq. 71 yields for the downstream end C
(AF + Ax"} °C ~ Ax' CD = ~AT~ (73)
and for the upstream end
L 4. f
Ax °C lAt Ax; °D At (74)
With the discharge and flow cross-sectional area given from the output of
the ISS Model, and c, and c, known from the initial condition or previous
o A
time computations, Eqs. 73 and 74 can be solved for c^, and c at each
**s . i)
time level using Cramer's rule. This procedure provides the pollutographs
at both ends of sewers which can also be used for water quality com-
putation at sewer junctions.
For a sewer junction, conservation of mass of the pollutant
gives
(75)
in which c is the pollutant concentration for the volume of water s in
the junction; c. is the concentration of discharge Q. from the i-th
sewer into the junction, Q being positive for inflow and negative for
outflow; and c. is the concentration for the j-th direct inflow Q..
Writing Eq. 75 in finite difference form and solving for the con-
centration at the present time level one obtains
1 , , At ,r v ~ \ i /-,^N
r— [c h + — (I c.Q. + i c.Q.)J (76)
P°° i X X j 3 3
108
-------
in which A is the constant horizontal cross sectional area of the
junction being considered; h is the depth of water in the junction; and
the subscripts (p) and (o) denote respectively the present and previous
time levels.
VIII-2. Program Description and Data Preparation
As mentioned in the preceding section, the water quality
computations for the surface runoff and sewer system are handled separately.
The computation for the surface runoff quality is done as a part of the
surface runoff model whereas that for the sewer is done through a
program supplement to the ISS Model. Since the Illinois Surface Runoff
Model Program has been described in Sec. VI-4, only the quality part of
the program is described here.
(A) Program description and data preparation for surface runoff quality
computation. - Because the discharges and flow cross-sectional areas of
the gutter and subcatchment flows required by Eq. 72 for quality com-
putation are computed in the surface runoff model but not provided as a
part of the surface runoff output, it is more advantageous to integrate the
quality computation into the surface runoff program than to separate it.
The computer program for the surface runoff allows the consideration of
two pollutants at a time. The- quality of subcatchment flow is computed
in subroutine OVLFLO using Eq. 71. At the upstream end of a -sub-
catchment strip, the concentration of each pollutant is assumed to remain
the same as the initial concentration. Knowing the flow from the
subcatchments, the water quality computations for gutter flow are made
in the MAIN program. In Eq. 71, the total lateral flow into a gutter from
all the components (subcatchments, direct rainfall etc.) is computed at
every time level and averaged over the length of the gutter. The
concentration in lateral inflow is evaluated as the average of
109
-------
component flow concentrations weighed with respect to flow rate. The
upstream boundary condition for a gutter is computed in subroutine
UPSBO. The concentration of each pollutant at the upstream end of the
gutter is evaluated at each time level as the average of concentrations
weighed over the flow rate of the carry-overs from the immediately
upstream inlets. When there is no upstream inlets, the concentration
at the upstream end is equal to the initial value at all time levels of
computation. The quality of flow into an inlet is assumed to be the
same as the corresponding gutter outflow. When there are more than one
gutter inlets discharging into the same catch basin or sewer node, an
average pollutograph for each pollutant weighed with respect to the inlet
flow rates is computed for the sewer node in subroutine SWRJNT.
The data required for the quality part of the surface runoff
model consists of the initial concentration of each pollutant in
subcatchment strips and gutters. The input format and the preparation
of data cards has been described in Sec. VI-4. The output from the computer
program includes the ordinates of direct inflow hydrographs and the
corresponding pollutographs for all the sewer nodes printed out at equal
time intervals.
(B) Description of sewer system water quality model. - The sewer system
water quality model is programmed in Fortran IV language as a supplement
to the ISS Model. The input to the computer program includes the data on
the depth and discharge at the entrance and exit of each sewer and the
volume of water at each storage junction at given times as provided by
the output of the ISS Model. In addition, the pollutographs as obtained
from the surface runoff quality computation, and the direct inflow
hydrographs, if any, and the sewer system layout are also input into the
110
-------
program for runoff quality routing in the sewer system, using Eqs. 73,
74, and 76. The output from the computer program consists of the
pollutographs representing the time variation of pollutant concentration
at the sewer system outlet, at the storage junctions, and at the entrance
and exit of all sewers.
The computer program of the Illinois Sewer System Water Quality
Model allows the consideration of two different pollutants at a time.
The sewer system may consist of as many as 100 sewers. Arbitrary
identification numbers can be assigned to the sewer nodes (junctions
and manholes). Because the backwater effect is already accounted for
in the quantity computation, the sewer system is not restricted to
tree-type networks. The computer program consists of approximately
200 statements and the storage requirement is 300 K. When a sewer
system consisting of more than 100 sewers is to be considered, the
program should be modified by simply changing the DIMENSION statements.
This may cause an increase in storage requirement. The computer program
is listed in Appendix C and the computational logic is shown
schematically in Fig. 25.
(C) Data preparation for sewer system water quality model computer
program. - The data deck for the sewer system water quality model
s
consists of three sets of cards. The first set is a simple card
describing the general information on the sewer system and input data.
The order and format of this information are as .follows: The total
number of sewers in the system (15) ; total number of reservoir-type
junctions in the system (15); number of points used to describe each
of the discharge and stage hydrographs at the entrance and exit of
each sewer (15) ; number of points used to describe direct inflow hydro-
graphs at each of the sewer nodes (15) ; identification number of the
111
-------
SEWER SYSTEM
CONSIDER ARBITRARILY
SELECTED FIRST SEWER
_L
READ SPECIFIC DATA FOR
THE SEWER CONSIDERED
TIME «= ZERO
SET CONCENTRATION OF EACH
POLLUTANT EQUAL TO INITIAL
VALUES AT ENTRANCE AND EXIT
-»JTIME " TITC + AT
GREATER THAN INPU
HYDROGRAPH
URATIONS
HAVE ALL SEWE
IN THE SYSTEM BEEN
CONSIDERED?
NO
COMPUTE CONCENTRATION OF
EACH POLLUTANT AT EXIT
AND ENTRANCE AT NEW TIME
LEVEL AND STORE VALUES
PRINT OUT TIME, DISCHARGE
AND POLLUTANT CONCENTRATIONS
AT ENTRANCE OF EXIT
REPLACE OLD CONCENTRATION
VALUES BY NEW VALUES
CONSIDER ARBITRARILY SELECT-
ED FIRST STORAGE JUNCTION
YES
READ DATA FOR THE
JUNCTION CONSIDERED
I TIME "ZERO
SET CONCENTRATIONS EQUAL
TO INITIAL VALUES
CONSIDER
NEXT JUNCTION
I
NO
IS TIME
GREATER THAN INPUT
HYDROGRAPH
DURATIONS?
HAVE ALL
JUNCTIONS BEEN
CONSIDERED?
COMPUTE CONCENTRATION OF EACH
POLLUTANT AT NEW TIME LEVEL
T
PRINT OUT VALUES OF TIME
STORAGE AND POLLUTANT
CONCENTRATIONS
*
REPLACE OLD VALUES OF EACH
QUANTITY BY NEW VALUES
STOP
Fig. 25. Flow chart for Illinois sewer system water
quality model computer program
112
-------
sewer node representing the system outlet (15); and the time interval
at which ordinates of the direct inflow hydrographs at the sewer nodes
are provided as input (F5.0).
The second set of cards are for sewer data. This set consists
of as many subsets as the total number of sewers, each subset corres-
ponding to a sewer. The order of these subsets is arbitrary within the
set. Each subset contains the following information: (a) The first
card in each subset provides the description of the corresponding sewer:
The identification number of the junction node at the upstream end of
the sewer (15); the identification number of the junction node at the
downstream end of the sewer (15); length of the sewer (F5.0): diameter
of the sewer (F5.0); initial concentration of the first pollutant at
the upstream end (F5.0) and at the downstream end (F5.0); initial con-
«
centration of the second pollutant at the upstream end (F5.0) and at
the downstream end (F5.0). (b) The data cards following the first one
in each subset contain the ordinates of the storage and flow hydrographs
at the entrance and exit of the sewer. These values need not be
provided at equal time intervals because the output from the ISS Model
on the depth and discharge at the entrance and exit of each of the
.sewers are at irregular time intervals. Each card contains three
groups of data punched in sequence and each group consists of the time
(F6.0), discharge at the entrance (F5.0), flow depth at the entrance
(F5.0), discharge at the exit (F5.0), and flow depth at the exit (F5.0).
The third set of cards are for the data for reservoir-type
junctions, containing as many subsets as the total number of
However, the direct inflow into the sewer junctions and the water stage
in the junctions can be obtained respectively from the surface runoff
and sewer routing models at constant time intervals. Therefore these
values are provided as input at regular time intervals.
113
-------
reservoir-type junctions in the sewer system. In each subset representing
a junction, (a) the first card contains in sequence the information on
the identification number of the junction node (15), cross sectional
area of the junction (F5.0), and the initial concentration of the first
pollutant (F5.0) and second pollutant (F5.0) in the junction; (b) the
data cards following the first card of the subset contain the ordinates
of direct inflow hydrographs into the junction, the ordinates of the
corresponding pollutographs and the water stage of the junction. All
*
values should be given at equal time intervals. Each card contains
three.groups of data punched in sequence and each group consists of:
time (F6.0), rate of direct inflow (F4.0), concentration of first
pollutant of the direct inflow (F4.0), that of the second pollutant (F4.0);
and the water stage in the junction (F4.0).
The output of the sewer quality model includes the time,
discharge, and concentration of the two pollutants at the sewer outlet
and at the entrance and exit of all the sewers, printed out at equal time
intervals under the appropriate headings, as well as the storage and
concentration of the two pollutants in each of the reservoir-type
junctions of the system at the same constant time intervals.
"See footnote in the preceding page.
114
-------
IX. GENERAL DESCRIPTION OF OTHER METHODS EVALUATED
Numerous rainfall-runoff "models" have been proposed by previous
investigators. Most of these models were developed for rural areas. These
include the Burkli-Ziegler and other similar formulas (Chow, 1962), the
monograph methods such as the methods of ARS-SCS (U.S. Soil Conservation-
Service, 1971), BPR (Potter, 1961), California (1953), Chow (1962), and
Cook (Hamilton and Jepson, 1940), and many of the hydrograph methods.
Among those models applicable to urban areas, some are strictly
for overland runoff prediction. These include Izzard's (1946) and Horton's
(1938) methods. The models which consider both overland and sewer flows
can be divided into two groups. The lumped system group, including the
rational method and the unit hydrograph method, treats the drainage basin
as a black box producing output (basin runoff) from given input without
considering what is happening to the flowing water within the basin. The
distributed system approach which includes most other urban runoff models,
routes the rainfall excess through the overland surface and sewers to produce
the runoff hydrographs.
It would be tedious to list and costly to compare in this study
all the methods applicable to urban drainage systems. Therefore, only eight
methods which are either most well known and widely adopted or having great
application potentials are evaluated in this report. The Illinois Urban
Storm Runoff method is described in Chapters VI and VII. The other seven
methods are briefly described in this chapter, following roughly their
relative order-of-complexity. Major features of these eight methods are
listed in Table 6. Presumably the most authoritive description of the
procedures in using a method is that by the original method developer.
Therefore, the steps for application of these methods are not repeated here.
115
-------
Table 6. URBAN RUNOFF PREDICTION METHODS EVALUATED
Model
Rational
Unit
hydrograph
Chicago
hydrograph
RRL
UCUR
Input data
Rainfall
Average In-
tensity over
the duration
Hyetograph
Hyetograph
Hyetograph
Hyetograph
Abstractions
Accounted by
runoff
coefficient
Infiltration by
ij> index or
Morton's
formula; other
abstractions,
if accountable
Infiltration by
Morton's for-
mula and depres-
sion storage by
an exponential
function
Pervious areas
produce no
runoff and all
rainfall on
impervious
areas becomes
runoff
Infiltration
from rainfall
only by
Horton's
formula and de-
pression storage
by an exponen-
tial function
Basin properties
iasin size
Jase flow
Overland surfaces;
lengths, slope, cross-
sectional dimensions
and roughness of
gutters and sewers
Areas of directly
contributing im-
pervious surfaces;
time of travel of
impervious areas
Length, slope, and n
for overland surfaces;
length of gutters;
diameter, slope and n
of sewers
Runoff routing
Overland
Izzard's
method
Gutters
Linear kinematic
wave storage
routing with
Manning's formula
Flow time - area method
Manning's
formula and
empirical de-
tention
storage
function
Continuity equa-
tion of steady
spatially varied
flow
ewers
Linear kinematic
wave storage
routing with
Manning's formula
or time offset
method
Reservoir routing
lagged by time of
travel in sewers
No routing,
lagged by time
of travel in
sewers
Output results
Peak discharge
Basin runoff
lydrograph
Basin runoff
hydrograph
Basin runoff
hydrograph
Runoff
hydrograph
Selected references
Chow, 1962;
ASCE and WPCF, 1969
Chow, 1964
Tholin and Keifer,
1960
Watkins, 1962;
Terstriep and Stall,
1969
Papadakis and Preul,
1972; Univ. of
Cincinnati, 1970
-------
Table 6. (continued)
Model
SWUM
Dorsch
Illinois
Input data
Rainfall
Hyetographs,
allows areal
variation
Hyetqgraphs
Hyetographs,
allows areal
variation
Abstractions
Infiltration
by Morton's
formula and
depression
storage .values
Infiltration
by Horton's
formula and de-
pression
storage
Infiltration
by Horton's
formula and
initial deten-
tion storage
Basin properties
Overland surface
length, width, rough-
ness n, slope and
percent imperviousness ;
length, slope, cross-
sectional dimensions
and roughness of
gutters and sewers
Overland surface and
gutter length, slope,
and roughness ', sewer
size, length, slope
and roughness
Length, width, slope,
and roughness of
overland surface
elements; length,
roughness, cross-
sectional dimensions
and slope of gutters;
type and dimensions
of inlets: length,
slope, roughness and
diameter of sewers,
size of manholes and
junctions
Runoff routing
Overland
Manning's formula
with uniform depth
Kinematic wave
model
Nonlinear kine-
matic wave with
Darcy-Weisbach's
eq.
Gutters
Linear kinematic
wave model,
storage routing
with Manning's
formula and
continuity
equation
Kinematic wave
model
Nonlinear
kinematic
wave with
Manning's
formula
Sewers
Improved non-
linear kine-
matic wave
model
St. Venant
eqs. with
partial back-
water effects
St. Venant
eqs. with back
water effects
Output results
Hydrographs of
runoff quantity
and quality, also
depth of flow
Runoff hydro-
graphs and depth
Runoff hydro-
graphs, also
depth and
velocity
Selected references
Metcalf & Eddy, Inc
et al., 1971
Heaney et al. , 1973
Klym et al. , 1972
Sevuk et al., 1973
-------
Those who want to use the particular methods should refer to the original
reports for detailed procedures.
IX-1. The Rational Method
The rational method is the oldest, simplest, and most widely adopted
method for storm runoff estimation (Chow, 1962). In the rational method, the
peak rate of storm runoff, Q , is estimated as
Q = CiA (91)
in which C is a dimensionless runoff coefficient; i is the rainfall intensity
and A is the size of the drainage area. The infiltration and other abstractions
from the rainfall is implicitly accounted for by the runoff coefficient. The
value of i is equal to the average rainfall intensity over a duration equal to
the so-called time of concentration. Details of application of the rational
method to urban storm sewer design can be found elsewhere (e.g. ASCE and WPCF,
1969; Yen et al., 1974).
The drawbacks of the rational method have been discussed by many
investigators (e.g., see Chow, 1964; McPherson, 1969). For urban storm runoff
quantity and quality control, the most serious drawback of the rational formula
is that it gives only the peak discharge, Q , and provides no information on
the time distribution of the storm runoff.
IX-2. Unit Hydrograph Method
Since Sherman (1932) proposed the concept of unit hydrograph, it
has been used to study the rainfall-runoff relationship for rural as well
as for urban areas. A unit hydrograph for a drainage basin is defined as
the discharge-time graph (hydrograph) of a unit volume of direct runoff
(usually expressed as unit depth) from the basin produced by an areally
and temporally uniformly distributed effective rainfall of a specified unit
118
-------
duration. The unit hydrograph for a drainage basin is obtained by reduc-
tion from previous rainfall and runoff data, by synthetic means, or by
transpose of the unit hydrograph from a neighboring basin of similar
physical characteristics. The unit hydrographs for different durations
can be obtained by direct derivation or by the S-hydrograph method. The
procedures to derive the unit hydrograph for a drainage area and to apply
it can be found in standard hydrology reference books (e.g., Chow, 1964).
With the unit hydrographs of different durations for a given
drainage basin known, the procedure to produce runoff hydrographs for given
rainstorms is as follows:
(a) The rainfall excess is first computed by subtracting abstrac-
tions from the total rainfall. For urban storm runoff
studies among the different abstractions usually only infiltra-
tion is considered although other abstractions can also be
included without causing much difficulty. Often the 4>-index
method of constant infiltration rate is used because of its
simplicity to estimate the infiltration, although Horton's
formula (Eq. 34) and other methods can also be used.
(b) Subdivide the rainfall excess obtained in (a) into a number
of small rainstorms of various durations each of which has
nearly uniform distribution of intensity of rainfall excess
over its duration. Determine the amount (depth) of rainfall
excess, duration, and the beginning time for each of- these
subdivided rainstorms.
(c) For each of these subdivided rainstorms, apply the unit hydro-
graph with a duration closest to the rainstorm's duration.
The ordinates of the unit hydrograph are multiplied by the depth
of the rainfall excess to give the runoff hydrograph for that
119
-------
component rainstorm. This procedure is repeated in sequence
for all the rainfall excesses of different depths and dura-
tions occurring at different times. The resulted runoff
hydrographs for all the component rainstorms are then added
linearly with appropriate time shifting (to account for the
different times for different rainfall excesses) to give the
combined runoff hydrograph due to the rainstorm.
(d) A base flow, such as dry weather flow, is added to the
combined runoff hydrograph obtained in (c) to give the total
runoff hydrograph.
For urban drainage basins which are often of the size smaller than
2
several km (several sq mi), reliable unit hydrograph usually cannot be
established because of lack of data. Also, for the unit hydrographs obtained
based on past records to be applicable, the drainage basin characteristics
must be time invariant, i.e., no significant changes of the physical proper-
ties of the basin in time. Furthermore, from the fluid mechanics viewpoint
the unit hydrograph theory may not be reliable for small drainage areas of
a few acres (hecters) or smaller.
IX-3. Chicago Hydrograph Method
The Chicago method is a steady-flow hydrograph routing method to
determine the time distribution of the quantity of storm runoff (Tholin and
Keifer, 1960). The method takes into consideration storages in gutters and
sewers. For a given drainage area and rainstorm, the infiltration is
computed by using Horton's formula (Eq. 34) and the surface depression
storage is computed by using an empirical function (Eq. 35). The computed
rainfall excess is then routed through the overland surface using a modified
Izzard's method. The gutter flow is routed using a storage routing method
120
-------
with Manning's formula. The same routing method is applied to the sewer
laterals and mains. To simplify the computation, a time-offset method
to subdivide the hydrographs for sewer routing was also proposed. However,
this time-offset method has no theoretical basis and is not used in the
present study.
Detail procedure of the Chicago method is described in a paper
by Tholin and Keifer (1960). Hand computation of the method is tedious and
time consuming. The Department of Public Works of the City of Chicago has
a computer program of the Chicago method.
IX-4. Road Research Laboratory Method
The British Road Research Laboratory (RRL) method is another
hydrograph routing method developed specifically for urban areas (Watkins,
1962; Terstriep and Stall, 1969). The method was developed for the purpose
of determination of design runoff hydrographs although it can also be used for
flow prediction purpose. The method assumes that pervious areas and
impervious areas not directly connected to the drainage system produce no
runoff, and all the rainfall on impervious areas directly connected to the
drainage system becomes runoff. A linear flow time-area concept similar to
that adopted in the development of Chow's method (1962) is used to establish
»-
the hydrographs for the impervious areas. A time of entry, similar to the
time of concentration for the rational method, is estimated by experience
as the time required for the directly connected impervious area to contribute
to the flow in.to the inlet catch basin. Terstriep and Stall (1969) proposed
to use a formula suggested by Hicks to compute the time of entry for overland
flow and Izzard's modification of Manning's formula for gutter flow. Thus,
in the original version the required overland and gutter data consist only
of the area of the directly connected impervious subcatchments and their time
121
-------
of entry, whereas for the Terstriep and Stall version the slope and length
of the overland surfaces and the slope, shape, and Manning's n for the
gutters are required. The inflows from different contributing areas of an
inlet are combined linearly with appropriate tine lag to give the inlet
inflow hydrograph.
The inlet hydrographs are then routed through the sewers using a
reservoir routing technique. The time of travel in a sewer is computed as
t = L/V, where L is the length of the sewer and V is the full-pipe flow
velocity computed by using the Chezy formula (Eq. 47) in Colebrook-White
form with c in m /sec given by
c = 17.72 log
148000
+ 8.04 x 10
/RS
R
(92)
where the surface roughness k is in mm and hydraulic radius R in m.
The inflow hydrograph of a sewer is the combination of the outflow
hydrographs, with appropriate time lag as computed by the time of travel, of
the inlets and sewers joining at the upstream end of the sewer being considered.
This inflow hydrograph is then routed through the sewer by using the con-
tinuity equation
li + Z2 Qi + Q2 S2 - si
2 2 At (93)
in which I is the inflow rate as given by the inflow hydrograph., Q is the
outflow rate of the outflow hydrograph, s is the storage in the sewer, and the
subscripts 1 and 2 refer to the beginning and end of the time interval At, respec-
tively. A storage-discharge relationship is then supplemented to Eq. 93 to give
the outflow rate. Originally, Watkins suggested to use the recession part
of recorded runoff hydrograph to establish the storage-discharge relationship.
In a later version, it was suggested to approximate this relationship using
Chezy's formula (Eq. 47) with c given by Eq. 92 assuming instantaneously the
sewer flow is steady and uniform with a slope equal to the sewer slope.
122
-------
A linear interpolation between the values of hydraulic radius and flow area was
suggested to avoid time consuming iterative solution.
Obviously, the RRL method would be most successful for drainage
basins of nearly equal in sizes of impervious and pervious areas and for
rainstorms of moderate intensities and durations. From the theoretical view-
point, the criticisms on the RRL method are numerous (Keeps and Mein, 1973).
For instance, in considering the sewer flow, the velocity is computed using
Eqs. 47 and 92 assuming a full-pipe flow, whereas in computing the sewer
storage, the sewer is assumed to have unlimited storage capacity as required
by Eq. 93. In general, the method does not take into account the actual
physics of storm water flow on the land surfaces and in sewers. The method
should be considered as empirical rather than theoretical.
IX-5. Cincinnati Urban Runoff Model
The University of Cincinnati Urban Runoff (UCUR) Model is basically a
hydrograph routing model developed under a grant from the U.S. Environmental
Protection Agency (University of Cincinnati, 1970; Papadakis and Preul, 1972).
The model has many similarities to the Chicago method and hydraulically it
is essentially a linear kinematic-wave model. As in the Chicago method,
infiltration is subtracted from rainfall using Morton's fprmula (Eq. 34)
(Eq. 34) with the aid of Jens' (1948) curves, and the depression storage
by using the same exponential function (Eq. 35) as in the Chicago method.
No infiltration is allowed from water stored on the land surface. Intercep-
tion and evapotranspiration are neglected. The resulted hyetograph of
rainfall excess can then be used for routing. The basin is subdivided
into a number of subcatchments each having homogeneous infiltration
characteristics. At any instant the rainfall is assumed uniformly distri-
buted over the entire basin but the abstractions are different for different
subcatchments.
123
-------
In routing the storm water, the overland flow is computed by using
Manning's formula (Eq. 46) coupled with an empirical detention storage function.
Assuming uniform flow with depth equal to the hydraulic radius, Manning's
formula combined with the continuity equation gives a relationship for detention
storage per unit width of the overland strip under equilibrium condition,
de'
d = C
nL)°'6L
e e S0.3
in which i is the overland flow supply rate which is equal to the intensity
of rainfall excess; n is the Manning roughness factor; L and S are the length
and slope of the overland strip, respectively; and C is a dimensional co-
efficient, equal to 0.625 sec/m " in SI system with length in m and time in
sec, and equal to 0.437 sec/ft ' in English system or 0.00073 if i is in
in./hr. The empirical detention storage function to be used together with
Eq. 46 to solve for the overland flow hydrograph is
d 0.3
h = d + 0 .6d <-) (95)
in which d is the detention storage and h is the flow depth at the exit of
the overland strip. During recession the ratio d/d is assumed to be unity.
It should be noted here that in UCUR method the overland flow depth does
not include the depth of water for depression storage.
The gutter outflow is assumed equal to the sum of the upstream
inflow, Q , and lateral inflow from overland at the same time increment; i.e.,
Q = Q +q?L (96)
U X/
in which q is the overland flow supply per unit length of gutter and L is the
A/
124
-------
gutter length. No time lag of gutter flow is considered. Thus, despite the
required detailed data for overland surfaces (length, slope, and Manning's n
of the overland surfaces), the required gutter data is rather simple: only
the length of the gutter. Neither the slope nor the cross section of the
gutter is needed. In Eq. 96 it is implicitly assumed that the overland
surface connected to the gutter is approximately rectangular in shape.
Because the effects of the slope, shape, and surface roughness of the gutter
are not accounted for, it can be expected that the UCUR method is not re-
liable when the gutter storage is important, such as the cases with small
gutter slope and street cross slope and with heavy rainstorms.
For sewer flows, no pipe storage is considered. The hydrographs
are simply lagged, without changing shape, by an average flow velocity, V ,
ciV
which is the weighed velocity of the flow velocity of the sewers computed
by using Manning's formula; i.e.,
The procedure has been programmed for digital computer solutions.
The method is rather tedious and complicated in view of the assumptions made
to simplify the hydraulic routing of storm runoffs and the accuracy that
the method can provide. Keeps and Mein (1973) modified the program to reduce
the computer time and to allow defined sewer networks. They used Newton-
Raphson's iteration instead of trial-and-error to solve the overland and
sewer flow equations and omitted the weighed average velocity computation
(Eq. 97). Their modifications substantially improve the applicability of
the method.
125
-------
IX-6. EPA Storm Water Management Model
The Storm Water Management Model (SWMM) was first developed jointly
by Metcalf & Eddy, Inc., University of Florida, and Water Resources Engineers,
Inc. (1971) under the sponsorship of the U.S. Environmental Protection Agency.
Subsequent separate modifications and improvement of SWMM by the University
of Florida researchers and Water Resources Engineers resulted in the EPA
SWMM version and the WRE SWMM version, respectively. The version compared
in this report is the nonproprietary EPA SWMM (Heaney et al., 1973).
SWMM is a comprehensive urban storm water quantity and quality
runoff prediction and management simulation model. The surface runoff is
routed by using a linear kinematic approximation and the sewer routing is
an improved nonlinear kinematic wave model incorporating some features of
the nonlinear quasi-steady dynamic-wave model. Because SWMM is probably
the best documented one among the recently developed models, details of the
EPA SWMM can easily be found in reports by Metcalf & Eddy, Inc. et al.
(1971) and Heaney et al. (1973).
In the five methods previously described in this chapter, at a
given time only one hyetograph is permitted for a drainage basin. In other
words, the rainfall is assumed uniformly distributed over the entire basin
and no areal variation is permitted. In SWMM several hyetographs applied
to different subcatchments at any given instant can be specified. Like
other methods, interception and evapotranspiration are neglected.
Infiltration is accounted for by using Morton's formula (Eq. 34).
Different degrees of permeability and different infiltration parameters for
Horton's formula may be applied to different subcatchments. Infiltration
is allowed for water stored on the surface in addition to rainfall. Overland
flow is assumed not to occur until depression storage is filled. Also, if
126
-------
not specified, by default one-quarter of the impervious surface area is assumed
to be of zero depression storage.
Overland flow is assumed to be one- dimensional with its depth being
constant along its length at any given instant. A quasi-steady flow routing
is applied to the overland flow using the continuity equation (Eq. 93) and
Manning's formula (Eq. 46) with depth equal to hydraulic radius. Within each
time interval the flow is assumed steady and uniform. The subcatchment data
required include the length, width, slope, surface roughness, and infiltra-
tion parameters and percent imperviousness.
The gutter flow is routed using the quasi-steady storage routing
approach with the storage continuity equation and Manning's formula. Each
gutter is assumed to be fed uniformly along its length by the lateral flow
from overland surface.
The sewer routing is by a modified nonlinear kinematic-wave scheme.
The continuity equation (Eq. 42) and Manning's formula (Eq. 46) are used
with the slope assumed equal to friction slope, and the flow is assumed to
be steady within each time interval. The sewer can be of any cross sectional
shape. The continuity equation is put into a finite difference form as
follows
in which Q is the discharge; A is the flow cross sectional area; L is the
sewer length; the subscripts u and d denote respectively the upstream and
downstream conditions; and the subscripts 1 and 2 represent respectively
the conditions at the beginning and end of the time interval At. The time
derivative is weighted W at the downstream station and the spatial deri-
vative is weighted W at the end of At. The backwater effect is considered
127
-------
by including the convective acceleration term in the equation of motion and
the friction slope is computed by
h - - h,, V2, - V2
s =
S
f o L 2gL
To simplify the computations for various conditions, these equations are
normalized using values of flow rate and area for conditions of the conduit
flowing full. .
If the flow is supercritical, no routing is attempted and the
sewer outflow is assumed equal to its inflow. If the backwater effect is
expected to be small and the sewer is circular in cross section, the gutter
flow routing method can be used as approximation to the nonlinear kinematic
wave routing.
IX- 7. Dorsch Hydrograph-Volume-Method
The Dorsch Hydrograph-Volume-Method (Klym et al. , 1972
relatively unknown in the U.S. until recently. It is a flow prediction
model and in most aspects it is more sophisticated than any one of the
preceding six models briefly described in this chapter. It takes hyetograph
as its rainfall input and a statistical method can be adopted for determina-
tion of frequent urban rainstorms. Interception and evap'otranspiration are
neglected and infiltration is computed by using Horton's formula (Eq. 34),
allowing different parameters for different sub catchments.
The surface flow is routed by a kinematic-wave scheme (Eqs. 42
and 44) similar to the Illinois surface runoff model but Manning's formula
(Eq. 46) instead of Darcy-Weisbach's formula (Eq. 45) is used to give the
friction slope. No consideration is given on changes of n for shallow
depth and due to rainfall. Sewer flows are routed by using the St. Venant
equations (Eqs. 42 and 43) with Manning's formula (Eq. 46) to give the
128
-------
friction slope. Backwater effect of junctions are partially accounted for.
The large number of the sewer and junction flow equations are solved
simultaneously using an implicit finite difference scheme. Because the
Dorsch HVM is a proprietary model, users should refer to the model developer
for procedural details. Because of its relative sophistication and pre-
sumably better accuracy, the data requirements for the Dorsch method are
considerably more extensive than those for the preceding six models.
129
-------
X. EVALUATION OF MODELS
The seven urban storm runoff simulation models described in Chapter IX
and the Illinois Urban Storm Runoff Model are compared using four rainstorms on
the Chicago Oakdale Avenue Drainage Basin. Based on the measured hyetographs
the runoff hydrographs are predicted by using the eight simulation models (the
rational method gives only the peak discharge) and the results are compared
with the measured hydrographs. Such an evaluation using only a few measured
rainstorms of course is limited in scope and one should be cautioned on the
generalization of the conclusions.
Theoretically, it would be desirable also to compare the different
models on drainage basins having their sizes an order-of-magnitude bigger
than that of the Oakdale Basin. However, in view of the requirements
of data details and accuracy for the more sophisticated models and
the computer time required for large basins, comparison using many rain-
storms on large basins would be extremely costly and time consuming. For
applications some errors on the data and results may be acceptable, where-
as for the purpose of model evaluation and comparison excessive errors
are intolerable. Furthermore, from the hydrodynamic viewpoint it is more
desirable to evaluate the models with smaller basin size with sufficient
*•
and accurate details than to larger basins, because for the larger basins
the local effects tend to be averaged out and not as clearly reflected in
the basin outflow hydrographs as for the smaller basins. For the Oakdale
Avenue Basin, because of its small size and proximity, the writers were
able to check in detail the field situation needed for the more complex
models within the time and budget allowance. With the presently
available data, comparison using a larger basin would not provide additional
new conclusions. In fact, as will be discussed later, the testing on the
Oakdale Basin suggests that a testing of the surface runoff part of the
models on a smaller area with accurate data is highly desirable.
130
-------
X-l. Hydrographs for the Eight Methods
The four rainstorms together with their respective measured runoff
hydrographs for the Oakdale Avenue Drainage Basin in Chicago selected for the
comparison study of the eight methods are the rainstorms of May 19, 1959,
July 2, 1960, April 29, 1963 and July 7, 1964. These four rainstorms are
chosen based on the following reasons: (a) The measurements of rainfall and
runoff are supposedly reliable as recommended by Tucker (1968). (b) The
rainfall is relatively heavy and hence presumably the measurements are
relatively accurate. (c) The rainstorms have been used by previous investi-
gators for establishment of their methods or for testing. The measured hyeto-
graphs and hydrographs are taken from a report by Tucker (1968) and reproduced
in Figs. 26 to 29. The physiographic characteristics of the Oakdale Drainage
Basin have been described in Chapter V and supplementary information can be
found elsewhere (e.g. Tucker 1968).
(A) Rational method. - The computation for the rational method is shown in
2
Table 7. The time of concentration, t , for the 0.052 km (12.9 ac)
drainage basin is 23 min, determined using the monograph by Kerby (1959) for
the flow on the subcatchment draining into the sewer Junction 117 (Fig. 13)
plus the sewer flow time from Junction 117 to the basin outlet which is com-
puted by using Manning's formula assuming barely filled gravity pipe flow.
The actual time of concentration is probably shorter as indicated.by the
recorded hyetographs and hydrographs. Therefore, a time of concentration of
20 min is adopted here. The rainfall intensity i used in the rational formula
is the average intensity of the recorded rainfall over a duration equal to t .
The runoff coefficient C used is 0.60 which is the value one would adopt from
standard tables (Chow, 1962, 1964; ASCE and WPCF, 1969) corresponding to the
surface condition of the Oakdale Basin. As customarily done for the rational
method, no abstraction is made from the rainfall and no adjustment of the C
131
-------
a>
a
CC
May 19, 1959
T.
U
12
10
in
>»—
o
o
.c
U
-UCUR
Chicago
Dorsch
ILL
Unit Hydrograph
100
80
E
_c
60 ^
'in
c
Q>
40 —
"o
a
20 CC
0.2
in
°E
_c
0>
w
o
O.I
10
20 30 40
Time in min
50
Fig. 26. Hyetograph and hydrographs for May 19, 1959 rainstorm
132
-------
.= 3 —
c
0>
2 —
o
cr
—
—
Ill II 1
I
1
I
1
1
1
1
1
I
July 2, I960
—
—
111 IT i
— 80
— 60
— 40
— 20
E
v>
c
O)
o
or
m
«•-
o
ID
O»
w
Q
— 0.5
o
o>
O)
o>
in
O
20
40
Time in min
60
Fig. 27. Hyetograph and hydrographs for July 2, 1960 rainstorm
133
-------
tn
c
0>
^
_ I
o
«^-
c
a
IT
L
a>
o>
o
CO
April 29, 1963
.J
UJUI
60
E
20 o
c
'o
(E
o
a>
o»
w
a
.c
o
10
20 30 40
Time in min
Fig. 28. Hyetograph and hydrographs for April 29, 1963 rainstorm
134
-------
— 120
Unit Hydrograph
-Recorded
)cA\\
Fig. 29. Hyetograph and hydrographs for July 7, 1964 rainstorm
135
-------
Table 7. RATIONAL METHOD COMPUTATION
Area = 12.9 ac or 0.052 km2, Runoff Coefficient C = 0.60
Rains torm
Duration, min
in./hr
mm/hr
cfs
m /sec
1 May 19, 1959
i
±, from j. 0
to i 20
20
1.17'
29.7
9.1
0.26
July 2, 1960
15
35
20
1.77
44.9
13.7
0.39
April 29, 1963
13
33
20
0.75
19.0
5.8
0.16
July 7, 1964
18
38
20
1.29
32.8
10.0
0.28
136
-------
Table 8. COMPARISON OF URBAN STORM RUNOFF METHODS
^\. Rains torm
Method ^^\^
-^
Rational
Unit hydrograph
Chicago
RRL
UCUR
: EPA SWMM
DORSCH HVM
Illinois
Recorded
May 19, 1959
0 t
VP P
3
,- m
cfs — mm
sec
9.1 0.26
6.9 0.20 17.0
9.3 0.26 20.0
7.3 0.21 18.7
10.2 0.29 16.5
7.6 0.22 17.0
8.0 0.23 13.5
7.4 0.21 16.0
7.2 0.20 16.9
July 2, 1960
Q t
P P
3
cfs — min
sec
13.7 0.39
15.3 0.43 36.2
15.5 0.44 38.8
14.2 0.40 36.2
18.2 0.52 38.0
15.6 0.44 37.2
18.5 0.53 32.0
15.4 0.44 33.9
17.5 0.49. 32.2
April 29, 1963
Q t
P P
3
cfs — min
sec
5.8 0.16
6.0 0.17 31.7
8.2 0.23 34.0
4.4 0.12 33.8
7.8 0.22 31.5
5.7 0.16 31.5
6.8 0.19 29.0
6.7 0.19 30.3
6.7 0.19 30.0
— — — — — ^— ^— — — — ^— — —
July 7, 1964
Q t
P P
3
cfs — min
sec
10 ..0 0.28
10.7 0.30 38.0
12.2 0.35 35.2
11.5 0.33 38.8
8.8 0.25 36.0
:
9.9 0.28 36.6'
9.6 0.27 36.5
-------
value is made to account for the preceding rainfall or antecedent surface
wetting conditions. As shown in Table 8, the computed peak discharges by the
rational method can be considerably different from the measured values.
(B) Unit hydrograph method. - As discussed in the preceding chapter, only
in rare cases that sufficient data are available for urban drainage basins to
establish their respective unit hydrographs. Fortunately, the Oakdale Basin
is one of those few that its unit hydrographs can be obtained. From the
recorded data recommended by Tucker (1968), rainfalls of 10-min duration
were selected together with the corresponding hydrographs to establish the
10-min unit hydrographs. The eight 10-min unit hydrographs so obtained are
shown in Fig. 30. Based on these eight unit hydrographs the average 10-min
unit hydrograph is plotted as shown in Fig. 30 and is used in this study
for runoff prediction. The 1-min unit hydrograph is subsequently obtained
by using the S-hydrograph method (Chow, 1964) as shown in Fig. 31 and checked
against recorded data. The 2-min unit hydrograph can then be derived by
adding one 1-min unit hydrograph to another 1-min unit hydrograph with a
1-min time lag. It should be mentioned here that in Fig. 30 the deviation
among the eight hydrographs for the eight recorded 10-min rainfalls may
reflect, in addition to data accuracy and linear approximation of a nonlinear
physical phenomenon, the effect of changes in basin characteristics with
season and time.
In applying the unit hydrographs of different durations to the
four rainstorms to obtain the runoff hydrographs, the abstraction made from
the total rainfall to give the rainfall excess is a rather subjective matter.
As discussed in Sec. IX-2, in this investigation for the sake of consistency,
Morton's formula (Eq. 34) is used to estimate the infiltration with the values
of the parameters roughly the same as those used for the Illinois surface
138
-------
Discharge in cfs
8
s
s
I I I I
.o rr ro —
oo — *° to
Discharge in m3/sec
Fig. 30. Ten-min unit hydrograph for Oakdale Basin
139
-------
Discharge in cfs
Discharge in m3/sec
Fig. 31. S-Curve and one-min unit hydrograph for Oakdale Basin
140
-------
runoff model. Also, the same base flow is added to both methods.
With the rainfall excess and base flow known, the standard
procedure for using unit hydrographs is applied. The headings of the
table used for the computations is shown in Table 9 and the results are
shown in Figs. 26 to 29 respectively for the four rainstorms tested.
(C) Chicago hydrograph method. - As mentioned in Sec. IX-3, the Chicago
method actually has two versions. The time-offset version is not adopted
in this study because the sewer layout of the Oakdale Basin is not favorable
for using this version. Therefore, the storage routing version is used.
The hydrograph for the July 7, 1964 rainstorm by the Chicago method shown
in Fig. 29 was taken from the University of Cincinnati (1970) report.
Because the computer program for Chicago method could not be released by
the city of Chicago and only three other rainstorms were to be evaluated,
it is not worthwhile in this study to write a computer program and hence
the hydrographs for the other three rainstorms were computed by hand
calculations. The computed results are plotted in Figs. 26 to 28 for
comparison.
(D) Road Research Laboratory method. - The hydrographs by the British
Road Research Laboratory method for the four rainstorms tested have been
computed by other investigators. The July 7, 1964 rainstorm was taken
from Terstriep and Stall (1969) and plotted in Fig. 29. The July 2, 1960
rainstorm was taken from the University of Cincinnati (1970) report and
shown in Fig. 27. For the other two, the May 19, 1959 and April 29, 1963
rainstorms, the hydrographs were taken from James F. MacLaren (1974) and
replotted in Figs. 26 and 28. Of course the results would be different
if the directly contributing impervious areas are taken differently from
the values assumed by those investigators.
141
-------
Table 9. HEADINGS FOR COMPUTATION OF UNIT HYDROGRAPH METHOD
Time
Recorded
Rainfall
Infiltration
Rate
Average Rate
During Period
Rainfall
Excess
Component Rainstorm
Duration
Excess
Runoff = Rainfall Excess x Unit Hydrograph Ordinates
Component
Rainstorm 1
Component
Rainstorm 2
Component
Rainstorm 3
Base
Flow
Total
Runoff
N3
-------
(E) University of Cincinnati Urban Runoff Model. - The hydrographs of
rainstorms of July 2, 1960 and July 7, 1964 for the UCTJR model were
originally computed by the model developers (University of Cincinnati,
1970) and replotted here in Figs. 27 and 29. The hydrographs for the
other two rainstorms, May 19, 1959 and April 23, 1963, were computed by
James F. MacLaren (1974) for another comparative study and borrowed for
this study as shown in Figs. 26 and 28.
(F) EPA Storm Water Management Model. - The hydrograph of the July 2,
1960 rainstorm for SWMM shown in Fig. 27 was taken from the report by the
model developers (Metcalf & Eddy, Inc. et al., 1971). It was calculated
using the original version of SWMM. However, the quantity part of the
original SWMM is essentially the same as the later modified version called
EPA SWMM. Hence this hydrograph can be used without recalculation. The
hydrographs for the rainstorms of May 19, 1959 and April 29, 1963 were
taken from the study by James F. MacLaren (1974) and replotted in Figs.
26 and 28.
(G) Dorsch Hydrograph Volume Method. - The Dorsch HVM is the only
proprietary model compared in this study because the hydrographs for the
May 19, 1959, July 2, 1960, and April 29, 1963 rainstorms are readily
available (James F. MacLaren, 1974) and replotted in Figs. 26 to 28. For
the same reason of proprietary no attempt was made to use the Dorsch HVM
to calculate the hydrograph for the July 7, 1964 rainstorm.
(H) Illinois Urban Storm Runoff Method. - The Illinois Urban Storm Runoff
Method was used to compute the hydrographs for all the four rainstorms
tested and the results are plotted in Figs. 26 to 29. Since details of
this method are not available elsewhere they are described as follows.
143
-------
The measured rainfall recorded as the hyetograph is subtracted
by infiltration to give the rainfall excess. The rainfall excess is then
routed through the subcatchments using the nonlinear kinematic-wave
approximation. When the rainfall rate recorded on the hyetograph is smaller
than the infiltration capacity, the depth of water detented on the surface
supplements the supply to infiltration.
As mentioned in Chapter V, the subcatchments of the Oakdale Basin
consists of four types; roofs, lawns, sidewalks (including drive-ways), and
street pavements. The asphalt street pavements are 4.27-m (14-ft) wide (one-
half of street width) having an average cross slope of 2.7%. The infiltration
through the street pavement is rather small and the values of f and k for
Horton's formula (Eq. 34) used are 0.01 in./hr or 0.25 mm/hr and 10 hr ,
respectively. Without adequate information on the initial wetting conditions
the value of the initial infiltration f for the four rainstorms is assumed
o
to be 0.02 in./hr or 0.5 mm/hr. It was found that because of the small
infiltration through the pavement, the results are not sensitive to the
values of f and k used.
o
Theoretically, the rainwater falling on the street pavements can
be routed through the length of the pavement strips as described in Sec. VI-1.
However, because of the short pavement length perpendicular to the gutter
(4.3 m or 14-ft), the pavement flow time is less than half a minute and much
less than the time intervals used for gutter flow computations. Consequently,
such routing of flow through the street pavement would substantially increase
the computer time and cost with essentially no improvement in the accuracy
of the results. Furthermore, the kinematic wave routing method described
in Sec. VI-1 does not account for the backwater effect from the gutter and
hence the solution is only approximate. Therefore, it was decided not to
144
-------
perform the street pavement flow routing for the Oakdale Basin and the com-
puted rainfall excess falling on the pavement is assumed to be transformed
into lateral flow to the gutter directly.
As discussed in Sec. V-I, it is possible but impractical to con-
sider in urban runoff prediction the details of distributions of the areas
for roofs, lawns, and sidewalks separately and accurately. Therefore, in
the present sutdy for the Oakdale Basin for the purpose of flow routing these
three types of surfaces are treated together as homogeneous overland surfaces
and an average slope of 6% and hydraulic surface roughness of 0.01 ft or
3 mm is assumed.
As discussed in Sec. IV-2, theoretically, the final infiltration
capacity f and the exponent k in Morton's formula (Eq. 34) should be nearly
constant for the pavements as well as for the roofs, lawns and sidewalks.
The only reasons that these values would change are the seasonal changes
and changes of the drainage basin physical characteristics in the span of
5 years. However, sufficient data were not available to establish the
values of f arid k for the Oakdale Basin. Consequently a costly and time
consuming, trial-and-error process based on the recorded hydrographs was
adopted. This approach is neither accurate nor foolproof and it is highly
s
undesirable. This difficulty necessarily points out the importance of
the need for detailed accurate data on land surface conditions in storm
runoff prediction. For the four rainstorms tested, no suitable common
values of f and k were found. In fact, the computed runoff hydrograph
is sensitive to the values of f , k, and f assumed as shown in Fig. 32
for rainstorms of July 2, 1960 and July 7, 1964. Values of Horton's
infiltration parameters used for the computed hydrographs of the Illinois
Urban Storm Runoff Method shown in Figs. 26 to 29 are listed in Table 10.
145
-------
Table 10. VALUES OF INFILTRATION PARAMETERS USED FOR ILLINOIS
URBAN STORM RUNOFF METHOD SHOWN IN FIGS. 26 TO 29
'"•--<
"••--,. Rainstorm
Parameter \^
f in./hr
0
f mm/hr
0
fc in./hr
f mm/hr
c
k hr'1
May 19, 1959
1.00
i .,
25
0.15
3.8
10
July 2, 1960
0.55
•
14
0.05
1.3
15
April 29, 1963
0.80
I
20
2.5
15
July 7, 1964
1.00
25
0.5
12.7
6
146
-------
The variations of f and k estimated from recorded hydrographs for the
four rainstorms tested due in addition to the error of using average values
to represent the integrated effect over the entire basin of infiltration of
individual component areas of different types, may also actually reflect the
effects of movement of the rainstorms and nonuniform areal distribution of the
rainfall intensity, as well as differences in roof-top and other surface retentions.
In routing the runoff, the direction of the subcatchment flow is as-
sumed perpendicular to the gutter. The overland surface (particularly the street
pavements) of the Oakdale Basin, as in many American cities, is reasonably
homogeneous and that only limited number of subcatchment strips need to be
computed and other strips can simply use the result obtained. This simplifi-
cation considerably reduces the computer time and costs without sacrificing
the accuracy of the results. The initial condition for the overland flow
as well as for gutter flow is a depth of 0.0012 in. or 0.03 mm of water
with zero velocity.
The subcatchment. runoffs together with the direct rainfall consti-
tute the input into the gutter flow. The gutter infiltration is evaluated
using Horton's formula (Eq. 34) with the same values of f , k and f as for
c o
the street pavement since no better information is available, although the
model allows different sets of values for the pavement and gutter. The
gutter flows are routed from the upstream, i.e., from the gutter with the
highest elevation of the basin towards downstream, as shown by the arrows
in Fig. 13. The westward flow direction of Oakdale Avenue east of Junction
118 is opposite to the direction of the slope of the sewers underneath. At
the downstream end of a gutter, water flows into the inlet with a capacity
given by Eqs. 63 and 64, in which the discharge coefficients C, used are 3.0
and 0.6, respectively. Not all of the inlets are placed immediately from
the curb line and some may extend beyond the gutter proper at the pavement
147
-------
side. Such irregularities are not accounted for in the model. Also, the
circular grate inlets are approximated by rectangular inlets of equivalent
area in the computation. Moreover,' as mentioned in Sec. VI-3, should the
inlet discharge be smaller than the gutter discharge, the excessive flow
is assumed carried through the inlet continuing into the next gutter
immediately following.
The outflow hydrographs obtained from the surface runoff model
are then used as the input into the sewer system for routing using the ISS
Model to give the drainage basin runoff hydrographs. Since the sewer
junctions in the Oakdale Basin have small cross-sectional areas, they are
considered as point-type junctions in the ISS Model. Also, the flow
measurement flume located at the outlet of the Oakdale Basin is approxi-
mated by a 0.76-m (30-in.) diameter pipe in calculations. The computed
hydrographs for the four rainstorms tested using the Illinois Urban Storm
Runoff Method with the values of infiltration parameters for subcatchments
listed in Table 10 are plotted in Figs. 26 to 29 for comparison with the
recorded hydrographs and computed hydrographs by other methods.
X-2. Comparison of the Methods
The evaluation of the eight methods are made against the recorded
hydrographs for the four rainstorms tested. Of course there is no guarantee
on the measurement accuracy of the recorded results. In fact, the accuracy
of rainfall and runoff measurements for the Oakdale Basin data has not been
adequately established. For example, the sharp drop at the peak discharge
of the July 7, 1964 data and the strange shape of the recorded hydrograph
in the next five min appear to be somewhat suspicious. For the July 2, 1960
rainstorm the peak discharge occurred at the same instant as the end of the
maximum rainfall intensity, leaving no time lag between the two, is also
questionable. Nevertheless, in view of the general situation of poor quality
of field data for urban storm runoffs, the Oakdale data should be considered
as one of the best sets that one can utilize.
148
-------
The comparison of the methods can be made from two aspects. First
is the ability of the methods to reproduce the recorded hydrographs based on
the recorded rainfalls. Second is the computational time and cost required
for the methods. For the first, in addition to the hydrographs shown in
Figs. 26 to 29, the peak discharge Q and its time of occurrence, t , for the
four rainstorms tested are listed in Table 8.
As mentioned in Chapter IX the rational method gives only the peak
discharge and hence its comparison with other methods is limited. The rational
method is extremely simple and gives a reasonable accurate estimation of Q if
the runoff coefficient C and the time of concentration, t , used for rainfall
c
intensity determination are properly chosen. However, the choice of C and
t is more an art of judgement than a scientific precision determination, and
as customarily used neither C nor t takes into account the preceding surface
moisture condition or the intensity and areal distribution of rainfall.
The unit hydrograph method provides reasonable though not very
accurate runoff hydrographs, and it is simple, fast, and straight forward if
the unit hydrograph for the drainage basin is available. Because the unit
hydrograph theory itself involves the assumption of linearity (Yen et al. , 1969,
Yen et al., 1973), one cannot expect high accuracy prediction from unit hydro-
graphs, particularly for small urban drainage basins of a-few acres (hectares).
From the practical viewpoint, the biggest problem of using unit hydrograph is
its availability. For most urban drainage basins there are no data to es-
tablish the unit hydrographs. Occasionally, transposing the known unit hydro-
graph of a neighboring drainage basin of similar nature with appropriate
size and other adjustments may be used. This is indeed a practical approach
but one should not expect high accuracy from it.
149
-------
None of the methods evaluated consistently reproduces the re-
corded runoff hydrographs faithfully. In general the three most sophisticated
methods, namely the Illinois Urban Storm Runoff Method, the Dorsch HVM, and
the EPA SWMM usually give better results than the other methods. The Chicago
method, RRL method,aid UCUR method may give results considerably different
from the recorded hydrographs and from those predicted by the more sophisticated
methods. However, a more precise, detailed, and meaningful comparison of
accuracy is not possible at present because of the uncertainties involved in
the amount and area and time distributions of infiltration as discussed earlier.
As shown in Fig. 32, the predicted runoff hydrographs are quite sensitive to
the values of infiltration used. Although attempts had been made in this
study to use the same infiltration function as much as possible, differences
still exist for different rainstorms and for different methods. As pointed
out by Torno (1974), infiltration of an amount different from that for other
methods was used for the UCUR method for the July 2, 1960 and July 7, 1964
rainstorms. Presumably the agreement of the UCUR predictions with the re-
corded hydrographs would be poorer if the same infiltration is used.
Furthermore, at least for the Illinois Urban Storm Runoff Method
and presumably also for other sophisticated methods, the accuracy on the
details of the basin geometry have profound influence on Ithe shape of the
hydrograph as shown in Fig. 32a. During the investigation of the effect
of infiltration on hydrographs, it was thought that since the Oakdale Basin
is reasonably symmetric with respect to Oakdale Avenue, computer time and
cost can be saved by computing the surface runoff for one half of the basin
and then double the result for sewer routing as a first approximation
before further refinement is made. The rainstorm of July 2, 1960 was tested
with the actual geometry and symmetric approximation for infiltration k = 20,
150
-------
O)
"c
O I
•*- '
o
a:
o
20
July 2 , I960
I \
100
80 £
E
60
40
20
Q)
o
tr
15
to
•*-
u
c
O)
o
CO
10
fo = 0.55
f c = 0.02
k = 20.0
fo = 0.55
fc = 0.02
k = 20.0
Symmetric
f o = 0.75
f c = 0.50
k = 6.0
fo = 0.55
fc = 0.05
k = 15.0
0.5
0.4
u
o>
10
0.3
E
c
o
(A
0.2
O.I
20
40
Time in min
60
eo
Fig. 32. Sensitivity of computed hydrograph to Horton's infiltration parameters
(a) Rainstorm of July 2, 1960
-------
r 3
'in
c
0)
a
cr i
o
14
12
10
6
o
CO
July 7, 1964
120
IOO
8O •=
6O
v>
c.
o>
40 "5
2O
o
cr
n n n n n n
fo = 0.80
fc = 0.30
k = 2.0
fo = 1.00
f c = 0.5O
k = 6.0
0.3
o
a>
IQ
f0 = 0.55
fc = P.02
k = 20.0
0.2
£
c
a>
o>
O
-------
f = 0.5 mm/hr (0.02 in./hr) and f = 14 mm/hr (0.55 in./hr). The
difference between the two resulted hydrographs of assumed symmetric and
nonsymmetric cases are surprisingly large, 17% difference in the peak
discharge. Presumably when the amount of infiltration increases the
difference would decrease. Nevertheless this result indicates that for
the sophisticated models for overland flow routing the accuracy of the
basin geometry is one important factor to produce reliable hydrographs.
The.difficulty in obtaining reliable comparison for the different
methods because of the uncertainties in infiltration suggests that further
study should be done on this line and that an even smaller drainage basin
than the Oakdale Basin with accurate detailed data should be used for
testing, and perhaps a separate testing on the surface runoff and sewer
runoff be conducted. This separate testing of surface and sewer runoffs
would more positively identify the relative merit of the different methods.
At present, it is intended in this report to provide from a practical
viewpoint some useful information for an engineer to decide which method
he should choose to use for his particular project.
As mentioned earlier, the computation time and effort for the
rational method is minimal, and that for the unit hydrograph method is
also small provided the unit hydrograph is available. "For the other six
methods use of computer is highly recommended, and in fact the three most
sophisticated methods cannot be done without using a computer. The rela-
tive amount of computer time for the six methods varies depending on the
rainstorm and drainage basin. In addition to the experience gained in
this study, some information on computer time for certain methods have
been reported by Keeps and Mein (1973) and James F. MacLaren, Ltd. (1974).
Roughly for a drainage basin like the Oakdale with the rainstorms similar
153
-------
to the four tested, the process time for RRL and Chicago methods are of
the order of a tenth of a min on the University of Illinois IBM 360/75
system, with the latter slightly longer. The UCUR method requires a
computer process time more than one order of magnitude longer than the
RRL method with no obvious improvement in accuracy. The computer time
for EPA SWMM is about the same as that for UCUR method and the Dorsch HVM
is slightly more. The Illinois Urban Storm Runoff Method requires the
longest computer time among the eight methods evaluated, being about twice
as much as that for the EPA SWMM, or somewhat more than 20 min for each
of the four Oakdale rainstorms tested.
The evaluation made in this as well as other similar studies on
urban runoff methods is definitely of limited scope and not exhaustive,
and one should interpret the results with caution. Nevertheless, it can be
concluded that if only the peak discharge is required and a quick result
is expected without high accuracy, the rational method is the most suitable
one to use. However, if the runoff hydrograph is needed such as for flow
regulation and storm runoff pollution control purposes the rational method
is unacceptable and other methods should be sought. Among the seven methods
evaluated that produce runoff hydrographs, the Illinois Urban Storm Runoff
Method most likely will give the most accurate result and is recommended
if high accuracy is required and no restrictions on computational costs.
The drawbacks of the Illinois method are that (a) it requires a large
amount of computer time and hence costs; (b) it requires detailed data
on basin physical characteristics that are not required by the other methods;
and (c) at present the sewer routing model (ISS Model) allows only circular
pipes and for sewers having other cross sectional shapes the equivalent pipe
diameter giving similar depth-area or depth-hydraulic radius relationship
should be first determined.
154
-------
The Dorsch HVM, though perhaps less accurate than the Illinois
method, also requires less computer time, so it can be used if the program
is readily available. Nevertheless, its data requirement is also as
detailed and demanding as the Illinois method.
The EPA SWMM can be operated with much less required data on
basin physical characteristics. Although it may not produce consistently
relatively reliable hydrographs as the Illinois or Dorsch methods, it is
much cheaper to use and is well documented and relatively most well known.
It also has many other features that the Illinois and Dorsch methods do not
have. It is recommended as another useful practical method.
When less accurate result is acceptable but the entire hydrograph
is required, the unit hydrograph method should be used whenever it is
possible. It is relatively simple, cheap, and fast, and it offers an
accuracy at least comparable, if not better, than the Chicago, UCUR, or
RRL methods. It does not require the use of computer to give the hydrograph,
although using computer would save time and programming is rather simple.
However, if synthetic means or basin transposition has to obtain
the unit hydrograph, one should be careful on the reliability and accuracy.
Should the unit hydrograph not be available for the drainage basin
of interest and the required accuracy of the hydrograph is not very high,
the RRL method appears to be the next choice because its data requirement
is not very high and the computer time is rather short. The Chicago method
does not seem to offer anything better than the RRL method, in fact, it
often over-predicts the peak discharge, and it requires more basin data
and computer time. The UCUR method differs little from the Chicago
method and yet requires considerably more data and longer computer time.
So it should be considered only as the last resort if other methods can-
not be used.
155
-------
Finally, it should be mentioned here that both the ISS Model and
EPA SWMM have built-in mechanism to consider the effectiveness of using
in-line storage for runoff quantity and quality control purpose should it
be desired to do so. Presumably the Dorsch HVM can also consider this
storage effect. The total runoff volume of course simply corresponds to
the total area under the hydrograph. The shape of the hydrograph, i.e.,
the time distribution of the runoff, will be altered with different
in-line storage. An example of alternating a storm sewer system design
by including in-line storage utilizing the ISS model has been reported
by Yen and Sevuk (1975).
156
-------
XI. REFERENCES
Akan, A. 0. Unsteady Gutter Flow into Grate Inlets. M.S. Thesis, Dept. of
Civil Eng., Univ. of Illinois at Urbana-Champaign, August 1973.
American Public Works Association. Problems of Combined Sewer Facilities
and Overflows. Water Pollut. Control Res. Ser., No. WP-20-11, Federal Water
Pollution Control Administration, December 1967.
American Public Works Association. Water Pollution Aspects of Urban Runoff.
Water Pollut. Control Res. Ser. No. WP-20-15, Federal Water Pollution Control
Administration, January 1969.
American Society of Civil Engineers and Water Pollution Control Federation.
Design and Construction of Sanitary and Storm Sewers. ASCE Manual No. 37,
1969.
Anderson, J. J. Real-Time Computer Control of Urban Runoff. Jour. Hyd. Div.,
ASCE 96(HY1):153-164, January 1970.
Anderson, J. J., R. L. Gallery, and D. J. Anderson. An Investigation of the
Evaluation of Automation and Control Schemes for Combined Sewer Systems.
Tech. Report No. 4, OWRR Proj. C-2207, Dept. of Civil Eng., Colorado State
Univ., Fort Collins, Colorado, January 1972.
California Department of Public Works, Division of Highways. California
Culvert Practice. 2nd ed., 1953.
Chen, C. L. and V. T. Chow. Hydrodynamics of Mathematically Simulated
Surface Runoff. Civil Eng. Studies Hyd. Eng. Ser. No. 18, Univ. of
Illinois at Urbana-Champaign, 111., 1968.
Chow, V. T. Frequency Analysis of Hydrologic Data with Special Application
to Rainfall Intensities. Eng. Expt. Sta. Bull. No. 414, Univ. of Illinois
at Urbana-Champaign, 111., 1953.
Chow, V. T. Open-Channel Hydraulics. McGraw-Hill Book Co., New York, 1959.
^
Chow, V. T. Hydrologic Determination of Waterway Areas for the Design of
Drainage Structures in Small Drainage Basins. Eng. Expt. Sta. Bull. No. 462,
Univ. of Illinois at Urbana-Champaign, 111., 1962.
Chow, V. T. ed. Handbook of Applied Hydrology. McGraw-Hill Book Co.,
New York, 1964.
Eagleson, P. S. Dynamic Hydrology. McGraw-Hill Book Co., New York, 1970.
Engineering-Science, Inc. Characterization and Treatment of Combined Sewer
Overflows. Report submitted by the City and County of San Francisco Dept.
of Public Works to the Federal Water Pollution Control Administration,
November 1967.
157
-------
Envirogenics Company (Division of Aerojet-General Corp.)- Urban Storm
Runoff and Combined Sewer Overflow Pollution. Water Pollut. Control Res.
Ser. No. 11024 FKM 12/71, U.S. Environmental Protection Agency, December
1971.
Field, R., and E. J. Struzeski, Jr. Management and Control of Combined
Sewer Overflows. Jour. Water Poll. Cont. Fed., 44 (_7): 1393-1415, July
1972.
Field, R., and P. Szeeley. Urban Runoff and Combined Sewer Overflow
(literature review). Jour. Water Pollut. Control Fed., 46(6):1209-1226,
1974.
Field, R., and.P. Weigel. Urban Runoff and Combined Sewer Overflow
(literature review). Jour. Water Pollut. Control Fed., 45(6}:1108-1115,
1973.
Grayman, W. M., and P. S. Eagleson. Streamflow Record Length for Modelling
Catchment Dynamics. Report No. 114 Hydrodynamics Lab., MIT, Cambridge,
Mass., 1969.
Hawkins, R. H., and J. H. Judd. Water Pollution as Affected by Street
Salting. Water Resources Bull., 8(6):1246-1252, December 1972.
Heaney, J. P., W. C. Huber, H. Sheikh, J. R. Doyle, and J. E. Darling.
Storm Water Management Model: Refinements, Testing and Decision-Making.
Report, Dept. of Environmental Eng. Sci., University of Florida, Gainsville,
Florida, June 1973.
Keeps, D. P., and R. G. Mein. An Independent Evaluation of Three Urban
Stormwater Models. Civil Eng. Research Report No. 4/1973, Monash Univ.,
Clayton Vic., Australia, 1973.
Horton, R. E. The Interpretation and Application of Runoff Plot Experiments,
with Reference to Soil Erosion Problems. Proc. Soil Sci. Soc. Am. 3:340-349,
1938.
Izzard, C. F. Hydraulics of Runoff from Developed Surfaces. Proc. Highway
Res. Board, 26:129-146, 1946.
James F. MacLaren, Ltd. Review of Canadian Storm Sewer Design Practice
and Comparison of Urban Hydrologic Models. Unpublished Draft Report for
Canadian Centre for Inland Waters, 1974.
Jens, S. W. Drainage of Airport Surfaces — Some Basic Design Considerations,
Trans. ASCE, 113:785-809, 1948.
Kerby, W. S. Time of Concentration for Overland Flow. Civil Engineering,
29O): 174 (issue p. 60), March 1959.
Klym, H. , W. Koniger, F. Mevius, and G. Vogel. Urban Hydrological Processes.
(Presented in the Seminar Computer Methods in Hydraulics at the Swiss
Federal Institute of Technology). Zurich, Dorsch Consultants, Munich,
Germany, February 17, 1972.
158
-------
Linsley, R. K. , M. A. Kohler, and J.L.H. Paulhus. Applied Hydrology,
McGraw-Hill Book Co., New York, 1949.
McPherson, M. B. Some Notes on the Rational Method of Storm Drain Design.
Tech Memo. No. 6, ASCE Urban Water Resources Program, January, 1969.
Metcalf & Eddy, Inc., University of Florida, and Water Resources Engineers,
Inc. Storm Water Management Model. Water Poll. Cont. Res. Ser. No. 11024
Doc, Vol. 1-4, U.S. Environmental Protection Agency, 1971.
Papadakis, C. N., and H. C. Preul. University of Cincinnati Urban Runoff
Model. Jour. Hyd. Div., ASCE,. 98(HY1Q):1789-1804. October, 1972.
Philip, J. R. Theory of Infiltration. In: Advances in Hydrosciences,
ed. by V. T. Chow, 5:215-296, 1969.
Poertner, H. G. Existing Automation, Control and Intelligence Systems for
Metropolitan Water Facilities. Tech. Report No. 1, OWRR Proj. C-2207, Dept.
of Civil Eng., Colorado State Univ., Fort Collins, Colorado, January 1972.
Potter, W. P. Peak Rates of Runoff from Small Watersheds. BPR Bull. Hyd.
Design Ser., No. 2, April 1961.
Sartor, J. D., and G. B. Boyd. Water Pollution Aspects of Street Surface
Contaminants. Environmental Protection Technology Ser. No. EPA-R2-72-081,
U.S. Environmental Protection Agency, November 1972.
Sevuk, A. S. Unsteady Flow in Sewer Networks. Ph.D. Thesis, Dept. of
Civil Eng., Univ. of Illinois at Urbana-Champaign, Illinois, 1973.
Sevuk, A. S., and B. C. Yen. A Comparative Study on Flood Routing Computa-
tion. Proc. Internat. Synrp. on River Mechanics, 3:275-290, Bangkok,
January 1973.
Sevuk, A. S., B. C. Yen and G. E. Peterson. Illinois Storm Sewer System
Simulation Model: User's Manual. Research Report No. 73, Water Resources
Center, Univ. of Illinois at Urbana-Champaign, 111., October 1973.
Shen, H. W., and R. M. Li. Rainfall Effect on Sheet FloW over Smooth
Surface. Jour. Hyd. Div.3 ASCE, 99(HY5.):771-792, May 1973.
Shen, Y. Y., B. C. Yen, and V. T. Chow. Experimental Investigation of
Watershed Surface Runoff. Civil Eng. Studies Hyd. Eng. Ser. No. 29, Univ.
of Illinois at Urbana-Champaign, Illinois, September 1974.
Sherman, L. K. Streamflow from Rainfall by the Unit-Graph Method.
Eng. News Rec., 108:501-505, April 7, 1932.
Terstriep, M. L., and J. B. Stall. Urban Runoff by Road Research
Laboratory Method. Jour. Hyd. Div., ASCE, 95(HY6):1809-1834, November
1969.
Tholin, A. L., and C. J. Kiefer. The Hydrology of Urban Runoff. Trans.
ASCE, 125:1308-1379, 1960.
159
-------
Torno, H. Discussion of: Methods for Determination of Urban Runoff by
C. N. Papakakis and H. C. Preul. Jour. Hyd. Div., ASCE, 100(HY8):1179-
1180, August 1974.
Tucker, L. S. Oakdale Gaging Installation, Chicago - Instrumentation and
Data. Tech. Memo. No. 2, ASCE Urban Water Resources Research Program,
1968; available from NTIS as PB 182787.
Tucker, L. S. Control of Combined Sewer Overflows in Minneapolis-St. Paul.
Tech. Report No. 3, OWRR Proj. C-2207, Dept. of Civil Eng., Colorado State
Univ., Fort Collins, Colo., October 1971.
U.S. Army Corps of Engineers. Runoff from Snowmelt. Engineering and Design
Manuals EM 1110-2-1406, January 5, 1960.
U.S. National Weather Service. Rainfall Frequency Atlas of the United States.
Weather Bureau Tech Paper No. 40, 1961.
U.S. Soil Conservation Service. SCS National Engineering Handbook, Section 4,
Hydrology. p. 10.1-10.24, January 1971.
University of Cincinnati Department of Civil Engineering. Urban Runoff
Characteristics. Water Pollut. Control Ees. Ser. No. 11024 DQU, U.S.
Environmental Protection Agency, October 1970.
Watkins, L. H. The Design of Urban Sewer Systems. Road Research Technical
Paper No. 55, Dept. of Sci. and Ind. Research, London, 1962.
Weibel, S. R., R. J. Anderson, and R. L. Woodward. Urban Land Runoff as a
Factor in Stream Pollution. Jour. Water Pollut. Control Fed. , 36:914-929,
1964.
Yen, B. C. Methodologies for Flow Prediction in Urban Storm Drainage Systems.
Research Report No. 72, Water Resources Center, Univ. of Illinois at Urbana-
Champaign, 111., August 1973a.
Yen, B. C. Open-Channel Flow Equations Revisited. Jour, Eng. Mech. Div.3
ASCE, 99(EM5):979-1009, October 1973b.
Yen, B. C. Discussion of: Numerical Model of St. Lawrence River Estuary
by D. Prandle and N. L. Crookshank. Jour. Hyd. Div.3 ASCE, 101
-------
Yen, B. C., and A. S. Sevuk. Design of Storm Sewer Networks. (To appear
in Proa.j ASCE, 1975.
Yen, B. C., W. H. Tang, and L. W. Mays. Design of Storm Sewers Using the
Rational Method. (To appear in Water & Sewage Works, 1974.)
Yen, B. C., H. G. Wenzel, Jr., and Y. N. Yoon. Resistance Coefficients for
Steady Spatially Varied Flow. Jour. Hyd. Div. 3 ASCE, 98(HY8):1395-1410,
August 1972.
161
-------
XII. NOTATION
A = area;
B = parameter for exponential density function (Eq. 20) or gamma density
function (Eq. 21);
b = water surface width;
C = parameter for gamma density function (Eq. 21); also, constant;
also, runoff coefficient;
C, = discharge coefficient;
c = Chezy's roughness factor; also, concentration;
D = rainfall depth; also, diameter; also, hydraulic depth = A/b;
d = detention storage;
d = average rainfall depth per time interval (Eq. 3);
d. = depth of rainfall of the j-th time interval;
J
E(x) = expected value of x
F = cumulative infiltration expressed in depth;
f = infiltration capacity; also, Weisbach resistance coefficient;
f = final infiltration capacity;
f = initial infiltration capacity;
G = second moment arm of hyetograph (Eq. 5); also, a function of flow depth, = 8Q/3h;
g = gravitational acceleration;
,*
H = available head of flow;
h = flow depth;
I = inflow rate;
i = rainfall intensity;
j = index number;
k = decay rate of infiltration in Horton's formula; also, surface roughness;
L = length
M = daily snow melt;
162
-------
n = Manning's roughness factor; also, a number;
P = cumulative rainfall in depth;
P = daily rainfall;
Q = discharge;
Q = peak discharge;
q = lateral inflow per unit length of a;
q = qda, lateral discharge per unit length of flow;
'a
R = hydraulic radius;
IR = VR/v, Reynolds number;
S = slope;
S = friction slope;
S = sinO, bed slope;
s = depression storage supply rate; also, storage;
s = depression storage capacity expressed in depth;
T = mean daily temperature of saturated air at 10-ft level;
3.
t = time;
t = first moment arm of hyetograph (Eq. 4):
t =elapsed time between the end of a rainstorm and the beginning of the
following rainstorm;
t , = rainfall duration;
U = x-component of velocity of lateral flow;
V = flow velocity;
V(x) = variance of x;
W = gutter width;
x = longitudinal direction;
z = elevation above horizontal reference datum;
6 = angle:
v = kinematic viscosity of water;
a = perimeter bounding flow area A: and
a, = standard deviation
d
163
-------
APPENDIX A
LISTING OF COMPUTER PROGRAK FOR FREQUENCY
ANALYSIS OF HOURLY PRECIPITATION DATA
The frequency analysis program is written in Fortran IV language.
The program in its present form requires about 160K bytes of storage. The
computer time required for the analyses described in Section IV-1 of this
report is about 90 seconds. Storage and time requirements may vary con-
siderably depending on the particular application.
Program descriptions, input and dimension requirements are
given at the beginning of the program as comment statements.
The listing of the computer program is given below.
164
-------
C*«*t±*i>«»:*«**e*«***a**X***tt************»* ******************************
C***:«**4*«f»***00*****«***«**-**»«
-------
C PROGRAM 2 ... SDRUNG VALUES OF RAINSTORM PARAMETERS.
C
C * THIS PROGRAM PERFORMS THE FOLLOWING OPERATIONS
C
C OPERATION 2-1 ...
C
C PRINTS THE VALUES OF THE RAINSTORM PARAMETERS
C FOR ALL RAINSTORMS IN THE RECORD.
C
C OPERATION 2-2 ...
C
C SCRTS THE VALUES OF THE RAINSTORM PARAMETERS IN ASCENDING
C ORDER AND PRINTS THE RESULTS.
C
C OPERATION 2-3 ...
C
C FOR ANY GIVEN RAINSTORM PARAMEfER,SORTS THE VALUES OF THAT
Z PARAMETER ANO PRINTS THE RESULTS TOGETHER WITH THE
C CORRESPONDING VALUES OF THE OTHER PARAMETERS.
C
C OPERATION' 2-4 .. .
C
C FOR ANY GIVEN RAINSTORM PARAMETER,SORTS A GIVEN SUBSET
C OF THE VALUES OF THAT PARAMETER AND PRINTS THE RESULTS
C TOGETHER WITH THE CORRESPONDING VALUES OF THE OTHEP
C PARAMETERS.
C
C * ALL OPERATIONS ARE OPTIONAL.
C
f1 •»*._«• _«-•_ •*_«._«__»_ __««._____ _ _ -j nr mr mm - -r -r iim w^^r - - _u_ -.- r — -a. i— -im i_ _i_- -f -m am j -wr mj — .. j_ — —. —
C
C PR03RAM 3 ... ONE-WAY FREQUENCY ANALYSIS OF THE VALUES OF
C . RAINSTORM PARAMETERS.
C
C * THIS PROGRAM PERFORMS ONE-WAY FREQUENCY ANALYSIS ON ALL VALUES
C OR DN A SUBSET OF VALUES (CORRESPONDING TO A GIVEN SUBSET OF
C VALUES OF ANOTHER PARAMETER) OF ANY GIVEN RAINSTORM PARAMETER.
C
C OPERATION 3-1 ...
C ONE-WAY FREQJENCY ANALYSIS OF ALL VALUES OF A GIVEN
C RAINSTORM PARAMETER.
C
C OPERATION 3-2 ...
C ONE-WAY FREQUENCY ANALYSIS OF A SUBSET OF VALUES
C (CORRESPONDING TO A GIVEN SUBSET OF VALUES OF ANOTHER
C PARAMETER) OF A GIVEN RAI'STORM PARAMETER.
C
C * THE OUTPUTS FROM THIS PROGRAM ARE THE TA3LES OF FREOUENCIES
-------
C OF OBSERVATIONS 3VE1 GIVEN CLASS INTERVALS), RELATIVE FREQUENCIES
C (FREQUENCY DIVIDED BY THE TOTAL NJMBER OF OBSERVATIONS) , PR3BJ BI-
C LITY DENSITIES(RELATIVE FREQUENCY DIVIDED BY THE INTERVAL SIZEI
C AND NON-EXCEEDANCE PROBABIL IT I ESI 2U^UL ATI VE RELATIVE FREQUENCY).
C IN ADDITICN, MINI MJM, MAXIMUM, MEAN AND STANDARD DEVIATICN OF THE
C VALJES CONSIDERED A* E PRINTED.
C
C * ALL OPERATIONS ARE OPTIONAL.
C
£ — — — — _____ __—-._,—.__________________..___-_.. _________-___-.____-__..____— — —,—____-._
C
C PROGRAM 4 ... TWO-WAY FREQUENCY ANALYSIS OF THE VALUES CF PAIRS CF
C RAINSTORM PARAMETERS.
C
C * THIS PROGRAM PERFORMS TWO-WAY FREQUENCY ANALSIS ON THE VALUES OF
C ANY GIVEN TWO RAINSTORM PARAMETERS.
C
C * THE OUTPUTS FRCM THIS PROGRAM APE THE TWO-WAY TABLES CF
C FRE3UENCIES(NUv.bE^ DF OBSERVATIONS OVER GIVEN CLASS INTERVALS),
C REHTIVE FREQUENCIES tFREOJENCY DIVIDED BY THE TOTAL NUMBER OF
C OBSERVATIONS) , AND PROBABILITY DENSITI ES IRELATI VE FREQUENCY DIVIDED
c BY THE INTERVAL AREA).
c
C * THE EXECUTION CF PROGRAM <* IS OPTIONAL.
C
c ---------------------------------------------------------------------
C
C SUBROUTINES ...
C
C * THE FOLLOWING SUBROUTINES ARE USED IN THE PROGRAMS
C
C SORT
C TAB1
C TAB2
C
C * SEE SUBROUTINE LISTINGS FOR PROGRAM DESCRIPTIONS.
C
C DIMENSION REQUIREMENTS.
C
: * X( ) IS USED TO STORE THE HOURLY PR ECIPITAT ICN VALUES FCR ALL
C HOURS DF A SEASON. ITS DIMENSION SHOULD BE EQUAL TO OP. GREATER
C THAN THE NUMBER OF HOURS WITHIN THE SEASON.
C
C * 0< ) IS USED TC STO*E THE HOURLY PRECIPITATICN VALUES FCR A RAIN-
C STORM.ITS DIMENSION SHOULD BE EOUiL TO 0?. GREATER THAN THE NUMBER
C OF HOURS FCR THE LONGEST DURATION RAINSTORM IN THE RECORD.
167
-------
c
C * P( ) IS USED TO STCRE THE VALUES DF THE RAINSTORM PARAMETERS
C FOR ALL RAINSTORMS IN A SEASON.THE POW DIMENSION SHOULD BE E3UAL
C TO 3R GREATER THAN THE MAXIMUM NUMBER CF RAINSTORMS EXPECTED IN
C ANY ONE SEASON.
C
C * PA< 1 IS USED TO STORE THE VALUES OF.THE RAINSTORM PARAMETERS
C FOR ALL RAINSTORMS IN THE RECORD.THE ROW DIMENSION SHOULD BE EQUAL
C TO OR GREATER THAN THE MAXIMUM NUMBER OF RAINSTORMS EXPECTED FCR
C THE WHOLE RECORD.
C
C * PAH ) SHOULD HAVE THE SAME DIMENSICNS AS PA( ).
C * A( ) AND B! ) SHOULD BOTH HAVE THE SAME DIMENSION AS THE ROW
C DIMENSION OF PA( I .
C
C * FRE3( ),PCT( J,DEN< ),CDEN< » AND XX( ) SHOULD ALL HAVE THE SAME
C DIMENSION AND IT SHOULD BE AT LEAST ONE MCP.E THAN THE MAXIMUM
C ' NUMBER OF CLASS INTERVALS THAT WILL BE USED FCR ANY ONE OF THE
C CNE-WAY FREQUENCY ANALYSES.
C
C
C INPUT REQUIREMENTS.
f
C * INPUT STATEMENTS TOGETHER WITH THE DESCRIPTIONS OF THE INPUT
C PARAMETERS ARE LISTED BELOW IN THE ORDER THEY APPEAR IN THE
C PROGRAMS.INPUT DATA SHOULD BE READ IN THE SAME ORDER.
C
f —______ _,^______—._.. _____„.__ _____.. — •_.___ ____———— — — __—_——— — — ____ — ______....
C
C PROGRAM 1 ...
C
C *** INPUT 1-1
C READ<5,57) IDIMX.IOIMPA
C 57 FORMAT(21*)
C IDIMX=DIMENSION OF VECTOR X( ).
C *** IDIMPA=ROW DIMENSION OF MATRIX PA( ).
C
C *** INPUT 1-2
C READ(5,1) IUNIT
C READ(5,1> IPRSI
C 1 FOR1AT< ID
C IUNIT=1 IF ALL INPUTS ARE IN ENGLISH UNITS.
C IUNIT=2 IF ALL INPUTS ARE IN SI UNITS.
C IP ;SI = 0 IF ALL INPUTS ARE IN SI UNITS.
C IPRSI=0 IF ALL INPUTS ARE IN ENGLISH UNITS AND IF THE RESULTS ARE
C TO 3E PRINTED ONLY IN ENGLISH UNITS.
168
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
• C.
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
-
c
c
c
c
c
IPRSI=1 IF ALL INPUTS ARE IN ENGLISH UNITS AND IF THE RESULTS ARE
***
***
2
3
***
***
4
***
***
5
***
***
13
12
15
***
***
TO BE PRINTED BDTH IN ENGLISH ANO SI UNITS.
THE ATOVE APPLY TO ALL PROGRAMS.
INPJT 1-3
READ(5,2) NUMMCN
FOF-MATU2)
NUMMON=NJMBER OF MONTHS IN THE SEASON ( 01-12).
READ(5,3) (MONTHU), 1=1, NUMMCN)
FORMAT(I2)
MCNTH( )=MCNTH NUMBERS IN ORDER 10 1- 12= JAN-DEC ).
ONE CARD FOR EACH MONTH.
INPJT 1-4
READ(5,4) NUMYR.LASTYR
FORM AT (2 12)
NUMYR=NUM8ER OF YEARS IN THE RECORD
LASTYR=LAST TWO DIGITS OF THE LAST YEAR IN THE RECORD.
INPUT 1-5
READ (5, 5) IPR1
FORM AT (ID
IPR1=1 IF.FCR EACH S EASON, HOURLY PRECIPITATION VALUES APE
PRINTED F3R EACH DAY AND HOUR OF THE SEASON TOGETHER WITH
VALUES OF THE RAINSTORM PARAMETERS FOR ALL RAINSTORMS IN
SEASON. OTHERWISE, I PR 1 = 0
INPJT 1-6
READ(5,13J IYR,IMP,IOY,
-------
C 80 FORMAT (II)
C IPP?=1 IF EXECUTION OF PROGRAM 2 IS DESIRED.OTHERWISEtIPR2=0 :
C *** IF IPR2 = 0 NO OTHER INPUT CARDS SHOULD BE USED FOF. PfrOGR/M 2.
C
C OPERATION 2-1 ...
C
C *** INPUT 2-2
C REAO(5,81) IPR21
C 81 FORIAT(Il)
C *** IPR21=1 IF OPERATION 2-1 IS TO BE PERFORMED.OTHERWISE,IPR21=0
C
C OPERATION 2-2 ...
C . . . .
C *** INPJT 2-3
C 90 REA3(5,82) IPR22
C 82 FOP1ATU1J
C *** IP*22=l IF OPERATION 2-2 IS TO BE PERFORMED.OTHERWISEfIPR22=0
C
C OPERATION 2-3 ...
C
C *** INPUT 2--4
C 98 REAO<5,85) IP323
C 85 FORMAT(Il)
C *** I?R23=1 IF OPERATION 2-3 IS TO BE PERFORMED.OTHERWISEtIPR23=0
C
C *** INPJT 2-5
C READ(5,99I NUMT
C 99 FORMAT(12)
C NUMT=NUM3ER OF TIMES OPERATION 2-3 IS TO BE REPEATED.
C *** IF IPR23=0 DISREGARD THIS INPUT.
C
C *** INPUT 2-6
C REAO(5,9^) NOPAR
C 94 FORMAT(I2)
: NOP*f\ = THE NUM3ER OF THE RAINSTORM PARAMETER FOR WHICH OPERATION
C 2-3 WILL 3E CARRIED OUT.
C (NUMBERING DF PAINSTORM PARAMETERS IS DESCRIBED AT THE END OF
C THIS SECTION.)
C ONE CARO FCP EACH RE°ETITION.
C THE TOTAL NUMBER OF CARDS IS EQUAL TO NUMT.
C *** IF IPR23=0 DISREGARD THIS INPUT.
C
C OPERATION 2-4 ...
C
C *** INPJT 2-7
C 104 REAO(5t126) IPR24
C 126 FQR^AT(Il)
C *** IPR24=1 IF OPERATION 2-4 IS TO BE PERFORMED.OTHERWISE,IPR24=0
170
-------
c
C *** INPJT 2-8
C READ(5,127) NUMT
C 127 FORMAT (12)
C NUMT=NUMBEP. Or TIMES OPERATION 2-4 IS TO BE REPEATED.
C *** IF IPR24=0 DISREGARD THIS INPUT.
C
C *** INPUT 2-9
C REAO(5f107) N3°AR,XLLtXUL
C 107 FOR1AT(I2,3X,2F10.0>
C NCP4R=THE NUMBER OF THE RAINSTORM PARAMETER FOR WHICH OPERATION
C 2-4 WILL BE CARRIED TUT.
C XLL=LOWER LIMIT FOR THE VALUES OF THE PARAMETER CONSIDERED.
C XUL=U°PE* LIMIT FOR THE VALUES OF THE PARAMETER CONSIDERED.
C CNLV THOSE VALUES .GE . XLL AND .LT. XUL ARE COSIDERFD.
C ONE CARD FOR EAZH REPETITION.
C THE TOTAL NUM3ER OF CARDS IS EQUAL TO NUMT.
C *** IF IPR24=0 DISREGARD THIS INPUT.
C
C""""*"""~~" ~* — ——— ——————————— -» _ _ -—.——.—_—. ——— ___—— _________________•—_-,
C
C PROGRAM 3 ...
C
C *** INPJT 3-1
C READ(5,150) IPR3
C 150 FCR1AT( ID
C IPR3=1 IF EXECUTION OF PROGRAM 3 IS DESIRED.OTHER WISE,IPR3=0
C *** IF IPR3=0 NO OTHER INPUT CARDS SHOULD BE USED FOR PROGRAM 3
C
C OPERATION 3-1 ...
C
C *** INPUT 3-2
C READ(5,152) IPR31
C 152 FOPMAT(Il)
C *** IPR31=1 IF OPERATION) 3-1 IS TO BE PERFORMED. OTHERWISE,I PR31 = 0
C
C *** INPUT 3-3
C READ<5,15<») NUMT
C 154 FQRMATU2)
C NUMT=NUMBER 3F TIMES OPERATION 3-1 IS TO BE REPEATED.
C *** IF IPR31=0 DISREGARD THIS INPUT.
C
C *** INPUT 3-4
C READ(5,155) NO»AR,U93I 1),UBO<3),U30(2)
C 155 FORMAT(I2f8Xf3F10.0)
C NOP4R=THE NUM3ER OF THE RAINSTORM PARAMETER FOR WHICH OPERATION
C 3-1 WILL BE CARRIED OUT.
C UBO(1I=LOWER LIMIT Or THE FIRST CLASS INTERVAL.
171
-------
uem 3)*'j?p,- ^ M-iT, - - -<;T CLASS INTERVAL.
UBGi 2»=NlMBER CF CLA^b iNitKVALS.
C IF J60C1)=UP013),THE PROGRAM USES MINIMUM AND MAXIMUM VALUES
w Or THE RAINSTORM PARAMETER AS UBOl1) AND UBC(3)tRESPECTIV ELY.
C UBQI2) MUST INCLUDE TWO CLASS INTERVALS FOR THE VALUES UNDER
C AND ABOVE LIMITS.
C INTERVAL SIZE IS COMPUTED AS FOLLOWS
C
C XINTS£=(UBO(3)-UBO<1))/(UB0(2)-2.)
C
C A COUNT IS CLASSIFIED INTD A PARTICULAR INTERVAL IF THE VALUE IS
C .GE. THE LOWE* LIMIT OF THAT INTERVAL BUT .LT. THE UPPER LIMIT
C OF THE SAM£ INTERVAL.
C ONE CARD FCR EACH REPETITION.
: THE TOTAL NUMdE* OF CARDS IS EQUAL TO NUMT.
C =?*•* IF IPR31=0 DISREGARD THIS INPUT.
C
C OPERATION 3-2 ...
C
C *** INPJT 3-:>
: 153 REAO<5,ra^» IPP.32
C 18* FORMAT< II)
C *** IPR32=1 IF OPERATION 3-2 IS TO BE PERFORMED.OTHERWISEtIPR32=0
C
C *** INPUT 3-6
C READC5.185) N'JMT
C 185 FORMAT(I2)
C NUMT = ;NUMBER OF TIMES OPERATION 3-2 I S TO BE REPEATED.
C *** IF IPR32=0 DISREGARD THIS INPUT.
C
C *** INPUT 3-7
C READ(5f187) NOPARB, N?PAR,XLL,XULf 'JBO 11J tUBO( 3 ) t UBO t 2 )
C 187 FCR1AT12I2,6X,5F10.0)
C ONE-WAY FREQUENCY ANALYSIS WILL BE PERFORMED ON A SUBSET OF VALUES
C OF RAINSTORM PARAMETER NOPAR CORRESPONDING TO THE SUBSET OF VALUES
C OF RAINSTORM PARAMETER NQPAR3 WHICH ARE .GE. XLL AND .LT. XJL
C UBO(1),USO(2),UBO(3) ARE AS DEFINED FOR INPUT 3-4
; ONE CARD FCR EACH REPETITION.
C THE T3TAL NUMBER OF CARDS IS EQUAL TO NUMT.
C *** IF IPR32=0 DISREGARD THIS INPUT.
C
C~~~~~~*~""*" ———————————————————————————————— • —————-.^———————
C
C PROGRAM <» ...
C
C *** INPJT <»-!
C READ (5,200) I
C 200 FORMAT(Il)
172
-------
C IPR4=1 IF EXECUTION OF PROGRAM A IS DESIRED.OTHERWISE.IPR4 = 0
C *** IF IPR4=0 NO OTHER INPUT CARDS SHOULD BE USED FOR PROGRAM A
C
C «•** INPUT 4-2
C READ(5,203) NUMT
C 203 FORMAT!12)
C *** NUMT=NUMdER OF TIMES PROGRAM 4 IS TO BE EXECUTED.
C
C *** INPJT 4-3
C REA3<5,205) NOV(1),J302(1,1),UB02(3,1)
C REAO<5,205) NOVl2),UB02(1,2),UB02<3,2)
C 205 FOPMAT(I2,8X,2F10.0)
C NOV(1)=THE NUMBER OF THE FIRST RAINSTORM PARAMETER TO BE
C CROSS-TABULATED.
C NOV(2)=THE NUMBER OF THE SECOND RAINSTORM PARAMETER TO BE
C CROSS-TABULATED.
C UB02litJ)=LCWER LIMIT OF THE FIRST CLASS INTERVAL FOR THE
C J TH VARIABLE, J=l,2
C UB02(3,J)=UPPER LIMIT OF THE LAST CLASS INTERVAL FOR THE
C J TH VARIABLE, J = l,2
; IF UOB2CliJ)=UBC2<3,J), THE PROGRAM USES MINIMUM AND MAXIMUM
C VALJES OF THE VARIABLE J AS UB02(1,J) AND UB02(3,J),RESPECT IVELY.
C
C UB02<2,J)=NUMBER OF CLASS INTERVALS FOR THE J TH VARIABLE, J=l,2
C UB02(2,J) MUST INCLUDE F3* EACH VARIA3LE TWO CLASS INTERVALS FOR
C THE VALUES UNDER AND ABOVE LIMITS.
C IN THIS PROGRAK.NUM3ER OF CLASS INTERVALS FOR BOTH VARIABLES IS
C EQUAL TO 20.THAT IS, JB021 2, 1) ='JB02( 2, 2) = 20. HOWEVER .DESIRED
C INTERVAL SIZES FOR VARIABLES CAN BE OBTAINED BY PROPER
C CHOICE OF CORRESPONDING UB02(1,J) AND UB02C3.J) VALUES.
C INTERVAL SIZE FOR EACH VARIABLE IS COMPUTED AS FOLLOWS
C
C (UB02(3,J)-UB02(l,J)»/(UB02<2,J)-2.) J=I,2
C
C FOR EACH VARIABLES COUNT IS CLASSIFIED INTO A PARTICULAR INTERVAL
C IF THE V4LUE IS .GE. THE LOWER LIMIT OF THAT INTERVAL BUT .LT. THE
C UPPER LIMIT OF THE SAME INTERVAL.
C TWO CARDS FOR EACH REPETITION.
C *** THE TOTAL NUMBER OF CARDS IS EQUAL TO 2*NUMT.
C
^ __—.____-_* ^ — —H •«.— — —-—-— 1 •- •- _ -im -• — __»_«_••»—_— -• .— _i_ _•• • — r j- mji —-— _ -I-- "M. • --— —I _~ I ~M_ — _
C
C PREPERATICN OF THE CARD DECK OF HOURLY PRECIPITATION DATA.
C
C * THE MAIN INPUT TO THE FREQUENCY ANALYSIS PROGRAM IS THE CARD DECK
C OF HOURLY PRECIPITATION DATA.
C
C * THE DECK SHOULD BE PREPARED AS FOLLOWS
173
-------
c
C (1) FOR EACH DAY WITH PRECIP ITATION,2 CARDS ARE PUNCHED
C AS FOLLOWS
C
C CARD NO. 1
C COLUMNS
C
C 1-i LEAVE BLANK OR USE FOR STATION IDENTIFICATION PURPOSES.
C IN U.S.A. NATIONAL CLIMATIC CENTER FORMAT,CCLS. 1-2 ARE USED
C TO PUNCH STATE CODE AND COLS.3-6 ARE USED TO PUNCH STATION
C NUMBER.
C 7-S LAST TWO DIGITS OF THE YEAR (55=1955).
C 9-10 MONTH N3. (01-12=JANUARY-DECEMBER•.
C 11-12 DAY OF THE MONTH (01-31J.
C 13 CARD NUMBER (=1).
C 14-16 PRECIPITATION AMOUNT FOR HOUR ENDING 01 LOCAL STANDARD TIME.
C VALUES ARE PUNCHED AS INTEGERS TO ONE-HUNOREOTHS OF AN INCH
C OR TO ONE-TENTHS OF A MILLIMETER (021=0.21 IN. OR 2.1 MMI.
C FOR NO PRECIPITATION (THAT IS LESS TH«N 0.01 IN. OR 0.1 MM),
C LEAVE BLANK OR PUNCH '-BE',WHERE B INDICATES BLANK.THE
C LATTER IS USA-NCC PRACTICE.
C
r
Of •
c
C 47-49 PRECIPITATION AMOUNT FOR HOJR ENDING 12 LST.
C 50-80 BLANK.
C
C CARD NO. 2
C COLUMNS
C
C 1-6 SAME AS IN CARD 1.
C 7-3 SAME AS IN CARD 1.
C 9-10 SAME AS IN CARD 1.
C 11-12 SAME AS IN CA*D 1.
C 13 CARD NUMBER (=2).
C 14-16 PRECIPITATION AMOUNT FOR HOUR ENDING 13 LST.
C
C
C
C 47-49 PRECIPITATION AMOUNT FOR HOUR ENDING 24 LST.
C 50-53 LEAVE BLANK OR PUNCH DAILY TOTAL.LATTER IS USf-NCC PRACTICE.
C 54-57 LEAVE BLANK OR PUNCH MONTHLY TOTAL ON THE LAST DAY OF THE
C MONTH WITH PRECIPIT AT ION.LATTER IS USA-NCC PRACTICE.
C 58-78 BLANK.
C 79-80 NEXT DAY WITH PRECIMITATION.FOR THE LAST DAY OF THE MONTH
C WITH PRECIPITATION,THIS IS '01'.
C
C (2) FOR FIRST DAY 3F THE MONTH, THE TWO CARDS SHOULD ALWAYS .BE
174
-------
c
c
p
c
c
c
c
f
V
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
(3)
(V»
PUNCHED. IF THERE IS NO PRE: I PI TATION FOR THAT DAY, DAILY
TOTAL CAN BE PUNCHED AS '-BBB'.IF THERE IS ND PRECIPITATION
FOR THAT MONTH, MONTHLY TCT6L CAN BE PUNCHED AS «-BBB' .
THESE TWD ARE USA-NCC PRACTICES.
HOURLY PRECIPITATION DATA FOR RECORDING RAINGAGES IN THE
U.S. WEATHER SERVICE SYSTEM ARE AVAILABLE AT COST ON PUNCHED
CARDS FROM U.S. DEPARTMENT 3F COMMERCE, NATIONAL OCEANIC AND
ATMOSPHERIC ADMI NISTRAT ION, ENV IRCNMENT AL DATA SERVICE,
NATIONAL CLIM4TIC CENTER t FEDERAL 3UI LDI NG, ASHEVILLE ,
N.C. 28301.
SINCE THEY USE THE ABOVE FOPMAT, THE CARD DECK OBTAINED FROM
THEM CA'J OHECTLY 8E FED INTO THE PROGP AM. HOWEVER , IT SHOULD
FIRST BE CHECKED TO CORRECT MISTAKES AND FILL IN INCOMPLETE
INFORMATION.
THE ANALYSIS IS CARRIED CUT ON SEASONAL BASIS. THE LENGTH OF
A SEASON IS AT LEAST 1 MONTH AND AT MOST 12 MONTHS. A SEASON
CAN EXTErjO OVER TWO CALENDAR YEARS. THE CARD DECK SHOULD
CONTAIN ONLY THE DATA FOP THE SEASONS AND THE YEARS TO BE
ANALYZED. THE DATA SHOULD BE PLACED IN TIME ORDER. ANY SEASON
HAVING MISSING CATA SHOULD BE COMPLETELY REMOVED FROM
THE DECK.
NUMBERING CF RAINSTORM PARAMETERS.
RAIN
NO.
I
2
3
4
5
6
7
8
9
10
11
12
13
STORM PARAMETERS
DESCRIPTION
TIME BETWEEN RA INSTORMS.HR
RAINSTORM DURATION, HR (TIME BETWEEN THIS RAINSTORM AND THE
PREVIOUS ONE IS GIVEN BY RAINSTORM PARAMETER NC . 1J
TOTAL DEPTH 3F RAINSTORM, MM OP IN.
AVERAGE INTENSITY OF RAI NSTOR M, MM/HR OR IN./HR
STAMDARD DEVIATION OF PAINSTGPM DEPTH, *M OR IN.
FIRST MOMENT ARM OF THE HYETOGRADH WITH RESPECT TO THE
BEGINNING TIME OF THE RAINSTORM, HR
SECOND MDMENT ARM OF THE HYETCGRAPH WITH RESPECT TO THE
BEGINNING TIME QF THE R A INST3P v, , HR .SO .
DIMENSION A FOR TRIANGULAR REPRESENTATION OF HV ETOGRAPH, HR
DIMENSION B OF TPIANGLE.HR
DIMENSION H OF TFUANGLE, MM OR IN.
NONOIKENSIONAL STANDARD DEVIATION OF RAINSTCPM DEPTH
FIRST MOMENT MM OP THE NCMDIMFNSIONAL HYETOGRAPH
SECOND MOMENT ARM OF THE NOND IMENSIGNAL HYETOGRAPH
175
-------
- 1* DIMENSION A FOR TRIANGULAR REPRESENTATION OF THE
C NONDIMENSIONAL HYETOGRAPH
C 15 DIMENSION B FOR TRIAMGUL4R REPRESENTATION OF
C NONDIMENSIGNAL HYETQGRAPH
C 16 DIKENSION H FG* TRIANGULAR REPRESENTATION OF
C NONDIMENSIQNAL HYETOGRAPH
C
C
C PROGRAM 1 ...
f
DIMENSION X(2500),D(50),P(50,16),;>A{500,16),PA1(500,16),A(500I
DIMENSION 6(500) ,FRE0(150), PCTU50),OEN(150),CDEN<150),XX(150)
DIMENSION MONTH( 12),KYR( 12), KMOl12).KNUMOY112),IT(2^),IP(16)
DI ME NS IC N U B0 (3 ) , ST 4T S ( 5), NOV ( 2 ), U 3D 2 ( 3, 2 ) , STt>. T1 ( 3, 20) , S TAT 2 ( 3, 20)
DIMENSION XINSZ(2) |FRE01120 ,20),PCT1(20,20),DEN1(20,20),XX1(21)
DIMENSION XX2(21)
C
C *** INPUT 1-1
READ(5,57) IDIMX.IDIMPA
57 FORMAT(214)
C ***
C
C *** INPUT 1-2
READ(5,1) IUNIT
READJ5,!) IPRSI
1 FORMATU1)
Z ***
C
C *** INPUT 1-3
READ(5,2) NUMMON
2 FORMAT(I 2)
READ(5,3) (MONTHd ), 1 = 1, NUMMON)
3 FORMATU2)
C ***
C
C *** INPUT 1-^
REAO(5,4) NUMYR.LASTYR
4 FORMAT(212)
C ***
C
C *** INPUT 1-5
READ(5,5) IPR1
5 FORMAT(Il)
C ***
C
WRITE(6,6»
6 FORMAT('1')
176
-------
WRITE(6,7) NUMY*
7 FOR'4AT(//////////43X,' FREQUENCY ANALYSIS OF HOURLY PRECIPITATION 0
1ATA«///////57X,I2,« VEARS OF RECORD1 /////63X , 'MONTHS •/)
WRIT E( 6,8 ) (MONTH! I ) , I=1,NUMMON)
8 FORMAT(65X,I2)
WRITE(6,o)
WRirE(6,60)
60 FOR1ATJ /I X,' RAINSTORM PARAMETERS « ENGLISH UNITS ) •/ )
WRITE(6,48J
48 F3RMAT(2X,'i.TIME BETWEEN RAINSTORMS , HP «/2X, '2. RA I NSTORM DURATION,
1HR'/2X, '3. TOTAL DEPTH OF RAINSTORM , I N1 /2X , «4. AVERAGE INTENSITY GF
2RAINSTOPM, IN/HRV2X, «5. STANDARD DEVIATION OF RAINSTORM DEPTH, !N'/2
3X, '6. FIRST MOMENT ARM OF THE HYETOGRAPH WITH RESPECT TO THE BEGINN
4ING TIME OF THE RA INSTORM ,HR« /2X ,' 7. SECOND MOMENT ARM OF THE HYETO
5GRAPH WITH RESPECT TO THE BEGINNING TIME OF THE RAINSTORM, HP .SQ. ' I
WRITE(6,49)
49 FORMATC2X, '8. DIMENSION A FOR TRIANGULAR REPRESENTATION OF HYETOGRA
1PH,HR'/2X,«9.D!MENSI 3N B OF TR IANGLE, HP • /1X, • 10.DI MENSI ON H OF TP I
2ANGLE, IN1 /1X, ' 11 .NONTI MENSICNAL STANDARD DEVIATION OF RAINSTORM DE
3PTH' /1X, ' 12. FIRST MOMENT ARM OF THE NCNDI MENSICN* L H YETCGRAPH ' / IX ,
4«13. SECOND MOMENT ARM OF THE NONDI MENS ICNAL HYETC GRAPH •/ IX, • 14. U 1M
5ENSION A FOR TRIANGULAR REPRESENTATION OF THE NOKDIMENS ICNAL HYETD
6GRAPH'/1X, ' 15.0IMENS ION B FOR TRIftNGULAR SEPRESENTATnN OF NONDIME
7NSIOMAL HYETOGRAPH'/IX , "16. DIMENSION H FOR TRIANGULAR REPKESENTATI
8CN DF NOMDIMENSIONAL HYETOGRAPH1///)
WRITEC6.61 )
61 FORMAT(////1X,'RAINSTCRM PARAMETERS ( SI UNITS )•/)
38 FORMAT (2X ,' 1. TIME BETWEEN RAINSTORMS , HR '/2X ,' 2. RA I NSTCRM DURATION,
1HR'/2X, '3. TOTAL DEPTH OF RAINSTORM , MM' /2X ,' 4. AVER AGE INTENSITY CF
2RAINSTORM,MM/HR'/2X, '5. STANDARD DEVIATION CF RAINSTORM DEPTH, MM'/2
3X,«6. FIRST MOMENT ARM OF THE HYETOGRAPH WITH RESPECT TO THE BEGINS
4ING TIME OF THE RAINSTORM ,HR • /2X, ' 7. SECOND MOMENT ARM OF THE HYET3
5GRAf>H WITH RESPECT TO THE BEGINNING TIME OF THE P A INSTO"M,HR. SO. ' )
WRITE(6,39)
39 FORMATI2X, ' 8. DIMENSION A FOR TRIANGULAR REPRESENTATION OF HYETOGRA
1PH,HR'/2X, '9. DIMENSION B OF TR IANGLE , HR • /lx, « 10.DI MENSICN H CF TRI
2ANGLE,MM'/1X,' 11 .N3NID! MENSIONAL STANDARD DEVIATION OF RAINSTORM DE
3PTH1 /1X, • 12. FIRST MOMENT ARM OF THE NONDI MENSIONAL HYETCGRAPH' /1X,
4«13.SECCND MOMENT ARM OF THE NGNDI MENS IONAL HYETCGRAPH '/ IX, ' 14. D I M
5ENSION A FOR TRIANGULAR REPRESENTATION CF THE NONDIMENS IONAL HYET3
6GRA3H'/1X, ' 15. DIMENSION B FOR TRIANGULAR REPRESENTATION OF NGNDI ME
7NSIONAL HYETOGRAPH«/1X ,' 16. DIMENSION H FOK TRIANGULAR REPRESENTATI
SON OF NONDIMENSIONAL HYETOGRAPH' ///I
MC3=0
55 DO 9 I=1,IDIMX
177
-------
9 xm=o.o
c
11=0
12=0
13=0
10 12=11*12
IFU2.EQ.O) GO TO 20
IF(INXOY.EQ.l.AND.I MO.£0.MCNTKUNUMMCN)) GO TO 11
20 IFJIUNIT.EQ.2) GO TO 12
C
C *** INPUT 1-6
READ(5, 13) IYR,IMO.IDY,(X(I+12),1=1,241,1NXDY
13 FORMAT(6X,3I2,lX,i2F3.2/13X,12F3.2,29X,I2)
C ***
C
GO TO 14
C
C *** INPUT 1-6
12 READ(5,15) IYR, IMO, I DY, (X ( I «-I2 ), 1= 1, 24 ), INXDY
15 FOR1AT(6X,3I2,1X,12F3.1/13X,12F3.1 ,29X,I2)
C ***
C
14 I1=24-*(INXDY-IDYJ
IFII1.GT.O) G3 TO 10
GOrO(16,18f16,17,16,17,16,16,17,16,17fl6),IMO
16 N'JMDY=31
GO TO 19
17 NUMOY=30
GO TO 19
18 X1=IYR
X2=I YR/4
X3=Xl/4.
IF(X2.EO.X3) NUMDY=29
IFCX2.NE.X3) NUHDY=28
19 I1=(N'JMDY+1-IDY)*24
13=13*1
KYR< I 3) = I YR
KMOl I3»=IMO
KNUMDY(!3)=NUMDY
GO TO 10
C
11 IC=0
LC = 0
MC = 1
NC = 1
DO 59 1=1,12
IF(X(I).GT.O.» GO TO 21
IF(NC.NE.MC) GO TO 22
178
-------
IC=IO1
GO TO 59
21 IFCLC.GT.O) GO TO 23
P(MC,1)=IC
IC=1
23 LC=LC+1
0(LC)=X(I )
GO TO 59
22 P(M:,2)=tC
SUM1=0.
DO 24 J=1,LC
24 SUM1=SUMI + 0( J»
P(MCt3)=SUMl
P(MC,4)=P{MC,3)/P(MC,2)
SUM2=0.
SUM3=0.
SUM4=0.
DO 25 J=1,LC
XJ = J
SUM2=SUM?+D( J)*(XJ-0.5 )
SU»13 = SUM3*D( J»*( XJ-0.5)**2
25 SUM4=SUM4+(D( J)-PI HC,4) J**2
PCMC ,5)=SQRT(SUM't/P(MCf2J )
P(MCf6) = SU^2/P(MC,? J
P(MCt7)=SUM3/P(MCt3)+l./12.
P(MC,8)=3.*P( MC,6)-P(MCf2)
P(MC,9)=P(XC, 2)-P(MC,8)
PJ MC »10)=2.*P
-------
c
c
c
c
IFUPRl.EO.O) GO TO 27
19=0
WRITE<6,6)
DO 29 1=1,24
29 ITU )=I
IFCIUNIT.EO.l) GO TO 28
54 WRITE<6,30I
30 FOPMATJIX,'HOURLY PRECIPITATION DATA.'/)
WRITE16.31)
31 FORMATUX, 'VALUES ARE IN MILLIMETERS.'///)
WRITE(6,32>
32 FORMAT(67X,'HOUR ENDING1/)
WRITEt6,33) (ITU) ,1=1,24)
33 FORMAT(11X,24(3X,I2)/)
44 CONTINUE •
17=0
00 34 1=1,13
I4=KYR(I)
I5=KMO(I)
I6=
46 FORMAT (IX,12,IX,I 2,IX,I 2,3X,24F5.2)
34 CONTINUE
IFUUNIT.E0.1.AND.I9.EQ.1) GO TO 47
WRITE(6,6)
WRITE16,62»
62 FORMATCIX,'VALUES OF RAINSTORM PARAMETERS IN SI UNITS.'//)
GO TO 50
47 HRITE(6,6)
WRIT E(6,63)
63 FORMAT(IX,•VALUES OF RAINSTORM PARAMETERS IN ENGLISH UNITS.'//)
50 CONTINUE
WRITEI6.41) ( IP< I),1=1,16)
180
-------
41 FORMAT(4X,16(12,6X)/)
C
WRITE(6,42) t(P(I,J) ,J=1,16J,1=2,MCI)
42 FORMAT(16(1X,F7.2))
WRITE(6,56) MC2
56 F3RHAT(///1X,'NUMBER OF RAINSTORMS=»,14)
C
1FCIUNIT.E0.1.AND.I9.E0.1) GC TO 51
GO TO 27
C
28 WRIT£(6,30)
WRITE(6,43)
43 FORMATUX,'VALUES A3 E IN INCHES.1///)
WRITE(6,32)
HRITE(6,33) ( IT( I) ,1 = 1,24)
19=1
GO TO 44
51 IFIIPRSI.EO.O) GO TO 27
19=2
DO 52 1=1., 12
52 X(I)=X(I)*25.4
DO 53 I=2,WC1
P(I,3)=P(I,3)*25.4
P(I,4)=P(I,4)*25.4
PII,5)=PCI,5)*25.4
53 P(I,10)=P(I,10)*25.4
WRITE(6,6)
GO TO 54
C
27 IFUYR.LT.LASTYR) GO TO 55
C
C
C PROS RAM 2 ...
C
C *** INPUT 2-1
READ(5,80) IPR2
80 FORMAT (ID
C ***
C
IF(IPR2.EQ.O) GO TO 105
C
DO 83 J=l,16
00 83 1 = 1,HC3
83 PAH I,J» = PA(I,J)
C
C OPERATION 2-1 ...
181
-------
c
C *** 1NPJT 2-2
READ(5,81) IPR21
81 FO'.IATdl)
C ***
C
IFC IPR21.EQ.O) GO T3 90
C
110=1
19 = 0
K1=MC3
GO TO 93
C
C OPERATION 2-2 . . .
C
C *** INPUT 2-3
90 REA015.32I I PR22
82 FOR1ATII1)
C ***
C
IF(I PR22.EQ.O) GO TO 98
C
110=2
19=0
DO 95 J=l,16
DO 96 1=1 ,K1
A( I J=PA( I,J)
96 Bl I)=?A1 I, J )
CALL SORT(Ad) ,Kl,Btl ) )
DO 97 I=1,K1
97 PAlt I,J)=A(IJ
95 CONTINUE
GO TO 93
C
C OPERATION 2-3 ...
C
C *** INPUT 2-^
98 READ(5,85) IPR23
85 FOR>1AT(I1)
C ***
C
IF(I PR23.EO.OJ 30 TO 104
C
110=3
K1=MC3
111 = 0
182
-------
C *** INPUT 2-5
READ(5,99) NUHT
99 FORMATCI2)
C ***
C
103 19=0
IFCI11.GE.NUMT) GO TO 104
C
C *** INPUT 2-6
READ(5,94) NCPAR
94 FORMATU2)
C ***
C
D3 100 J=ltl6
DO 101 1=1,Kl
All)=PA(I,NCPAR)
101 B(I)=PA(I,J)
CALL SORT (A(H,K1,B(1I )
DO 102 1=1 ,K1
102 PAH I , J )=B< I )
100 CONfINUE '
111=111*1
GO TO 93
C
C OPERATION 2-4 ...
C
C *** INPUT 2-7
104 REAO<5,126) IPR24
126 FORMAT(Il)
C ***
C
IF(IPR24.EO.O) GO TO 105
C
110=4
111 = 0
C
C *** INPJT 2-8
READ(5,127) NUMT
127 FORMAT(12)
C ***
C
106 19=0
IF(Ill.GE.NUMT) GO TO 105
C
C *** INPJT 2-9
READJ5.107) NOPAR.XLL.XUL
107 FORMATU 2f8X,2F10.0)
C ***
183
-------
Kl=0
DO 108 1=1, MC3
COUNT=0.
IF(PA(I,NOPAR).GE.XLL.AND.PA< I ,NOPAR J .LT .XUL ) COUNT=1.
IFtCOUNT.EO.O.) GO TO 108
K1=K1+1
DO 109 J=l,16
109 PAK Kl,Ji=PA< I , J)
108 CONTINUE
DO 121 J=l,16
IF(J .EO.NOPAR J GO TO 121
DO 122 1=1, Kl
A(I) =PAK I ,NQPAR)
122 Ed )=PAU I,J)
CALL SCRT(A(1),K1,B(1) )
DO 123 1=1, Kl
123 PAK I,JI = B(I I
121 CONTINUE
DO 129 1=1, Kl
129 PAK I ,NrPAR)=A(I )
111=111+1
93 IFCIUNIT.EQ.l > 19=1
91 WRITE(6,6)
GO TQ( I10,lll,112r 113) ,110
110 WRITE(6,llH)
114 FORMATl IX,' VALUES OF RAINSTORM PARAMETERS.'///)
GO TO 120
111 WRITEC6.1151
115 FORMAT! IX,' VALUES OF RAINSTORM PARAMETERS SORTED IN ASCENDING ORDE
GO TO 120
112 WRITE(6,116) NOPAR
116 FORMAT! IX, 'VALUES OF RAINSTORM PARAMETER ',12,' SORTED IN ASCENOIN
1G ORDER AND CORRESPONDING VALUES OF CTHER PARAMETERS REARRANGED SI
2MULTANEOUSLY. •///)
GO TO 120
113 IFUUNIT. EG. 1. ANO.I9.EQ.il GO TO 118
IFUUNIT. EQ. 2) GO TO 125
IF«NOPAR.NE.3.AND.NOPAR.NE.'f .AND.NOPAR .NE. 5.AND.N&PAP. .NE . 10) GO TO
1125
XLL=XLL*25.4 *
XUL=XUL*25.
-------
3RRA1GED SIMULTANEOUSLY.'///)
GO TO 120
118 WRITE(6,119) NOP4R ,XLL,XUL
119 FDRHAT< IX, 'VALUES OF RAINSTORM PARAMETER ',12,' GREATER THAN OR EO
1UAL TO ',F7.2t' AND LESS THAN SF7.2,' (IN EN UNI TS) '/1X, 'SORTED I
2N ASCENDING ORDER ArlD CORRESPONDING VALUES OF OTHER PARAMETERS REA
SR^A-gGED SIMULTA'JEDUSLY.'///)
120 CONTINUE
IF(I UNIT. EO. LAND. 19. E0.1J GO TO 87
HRITE(6,130)
130 FORMAT(1X,'( SI UNITS )••///)
GO TO 88
87 WRITE(6,131 )
131 FORMATI1X.M ENGLISH UNITS )•///)
88 CONTINUE
WRITE(6,41) ( IP( I), 1 = 1, 16)
WRITE (6, 42) ( (PAKI.J) ,J=1,16),I=1,K1)
WRITE16.56) Kl
IFHUNIT.EO.l .AND.I9.EQ.1) GO TO 89
84 GO TO (90, 98, 103, 106), 110
89 IFUPRSI.EQ.O) GO TO 84
19=2
DO 92 1=1, Kl
PAH I,3) = PA1( 1,31*25.4
PAUI,4) = PAll 1,41*25.4
PAH I,5) = PA1(I,5)*25.4
92 PAK I ,10)=PA1( I, 10 >*25.4
GO TO 91
105 CONTINUE
c
C PRD3RAM 3 ...
C
C *** INPJT 3-1
READ(5,150) IPR3
150 FORM AT (11)
C ***
M
IF(IPR3.EO.O) G3 TO 151
C
DO 156 I=1,IDIMPA
156 A(I )=0.0
C
C CPESATION 3-1 ...
C
185
-------
C *** INPUT 3-2
REAO(5,152) IPR31
152 FORMATU1 )
C ***
C
I.FCIPR31.E0.01 GO TO 153
C
110=1
111=0
C
C *** INPUT 3-3
REA5(5f154) NUMT
154 FORM AT (12)
C ***
C
183 19=0
IF(Ill.GE.NUMT) GO TO 153
C
C *** INPJT 3-4
READ(5,155) NOPAR , U30 ( 1) ,UBO( 3) ,UBO( 2 )
155 FORHATJ I2,8X,3F10.0)
C ***
C
00 157 1*1, MC3
157 A(I)=1.0
CALL TAB1 (PAt A,NOPAR,UBOt FREQ.PCT, STATS t ID IMP A , 16)
111=111*1
GO TO 177
C
C OPERATION 3-2 ...
C
C *** INPUT 3-5
153 FEAD<5,184) IPR32
184 FORMAT(Il)
C ***
C
IF(IPR32.EQ.O) GO TO 151
r
110=2
111 = 0
C
C *** INPUT 3-6
REAO(5,185) NUMT
185 FORMATU2)
C ***
C
186 19=0
186
-------
IFdll.GE.NUMT) GO TO 151
C
C *** INPUT 3-7
READ<5,187) NOP4R8 ,NOPAR,XLL,XUL,UBO ( 1 ) ,UBO(3 ) ,UBO(2 )
187 FORMATC2I2.6X, 5F 10.0)
C ***
C
Kl = 0
00 188 I=ltMC3
IFIPAd .NOPAR6J . GE.XLL.AND.PAdtNOPARBKLT.XUL ) A( I) = 1.0
IFtACI ).EO.O.) GO TO 188
188 CONTINUE
CALL TABl(PA,A,NOPAR,U80iFREQfPCTf STATSt IDIMPA.16)
111=111*1
C
177 XnTSZ = tUBO(3)-UBC(l) J /( UBO ( 2 >-2.)
IUB02=UBO<2)
SUM=0.
00 158 I=1,IU302
PCT( I)=0.01*PCTtI)
DENC I)=PCT(I)/XINTSZ
SUM=SUM*PCT ( I )
158
160 IFdUNIT.EO. 1) 19=1
161 WRITE(6,6)
IFdUNIT.EO. LAND. 19. EO.l) GO TO 169
WRITE(6,170)
170 FORMATC61X,' ( SI UNITS )•//)
GO TO 171
169 WP,ITE(6,172 )
172 FOR1AT(58X,' ( ENGLISH UNITS )•//)
171 HRITEI6.159) NOPAR
159 FORMATdX.'FREQJENCV ANALYSIS FOR RAINSTORM PARAMETER «,I2/)
IF(IlO.EQ.l) GO TO 162
IFdUNIT.E0.1.AfMD.I9.E0.1J GO TO 164
IFdUNIT.EQ.2) GO TO 164
IFlviDPARI.NE.a.iND.NOPAP.B.NE.^.AND.NOPARB.NE.S.AND.NOPARB.NE.lO)
1GO TO 164
XLL=XLL*25.4
XUL=XUL*25.A
164 WRIT E16.165) NOP AR B, XLL.XUL
165 F3R^AT(1X,'ONLY VALUES CORRESPONDING TO VALUES OF RAINSTORM PARAME
ITE^ SI2,1 GREATER THAN CR EQUAL TO tfF7.2t' AND LESS THAN «,F7.2(
2' CCNS IDERED.V )
162 CONTINUE
HRITE16.167) UBO (1 » , UBC ( 3) ( UBO( 2) , XI NT SZ
187
-------
167 FORMAT (//IXt 'LOWER L IM IT = » t F7.2 , 5X, 'UPPER LI MI T = ' , F7. 2, 5X,« NUMBER
10F ;LASS INTERVALS=',F<».Of 5X, ' INTERVAL SIZ E= ' , F7. 2// )
WRirE(6,168) STATS <*******<<»*:*******
C
C
C
C
C
C
***
200
***
PROGRAM
INPJT 4
READ(5,
F3R^AT(
<» .
200)
11)
1PM
188
-------
IFCIPR4.EQ.O) GO TO 201
C
DO 202 I=1,IDIMPA
202 A(I)=0.0
111=0
C
C *** INPUT 4-2
READ(5,203) NUMT
203 FORM AT(12)
C ***
C
204 19=0
IF(Ill.GE.NUMT) GO TO 201
C
C *** INPUT 4-3
REAO(5,205) NOV (1),UB02(1,1),UB02( 3,1)
REA3(5,205) NOV12),UB02(1,2),UB02<3,2)
205 FORMAT(I2,8X,2F10.0)
C ***
C
UB02t2,l>=20.
UB02(2,2)=20.
C
00 206 1=1,MC3
206 A(I)=1.0
CALL TAB2(PA,A,NOV,UB02,FREQ1,PCT1,STAT1.STAT2,1 01 MPA, 16)
111=111+1
XINSZtl)=
189
-------
F3RMATUX,'RAINSTORM PARAMETER',I2fIX,•LOWER LIMIT*•,F7.2,5X,'UPP
1ER LIMIT=',F7.2,5X,'NUMBER OF CLASS INTERVALS=«,F4.0,5X,•INTERVAL
2SIZE=«.F7.2//I
213 CONTINUE
DO 216 1=1,21
IF(I.E0.1.CR.I.EQ.21) GO T3 217
XXK i )=uao2 (i,i)+( I-Z)*XINSZ(II
GO TO 216
217 XXK I )=100000.
216 CONTINUE
DO 218 1=1,21
IFU.EQ.1.CP..I.EQ.21) GO TO 219
XX2(I)=U302(1,2)+(I-2)*XINSZ(2)
GO TO 218
219 XX2II)=1COOOO.
218 CONTINUE
DO 215 K=l,4
L=(K-1)*5
HRITEC6.220)
220 FORM AT(///64X,1FREQUENCIES'//)
WRITE(6,221) NOV(2)
221 FOR>1AT(5X, 'RAINSTORM ',<*$*, 'RAINSTORM PARAMETER ',I2/>
WRITE(6,222) NOV(1), (XX2(I*L ) ,XX2(I*1+L),1 = 1,51
222 F3R^AT(^X,•PARAMETER ' , 12,7X,5IF7.2t2X,F7.2,3X)/)
WRITE (6,223) (XXI ( I ), XXK 1+1), FRE3 11 I , 1*L J ,FREQ1( I ,2*L) ,FRE01( I ,3*
1L),FREQ1U ,^f+L) , FREQ1(I,5*L), 1 = 1,20)
223 FOR1AT(1X,F7.2,3X,F7.2,9X,F10.A,9X.F10.4,9X,F10.^,9X,F10.A,9X,F10.
14)
215 CONTINUE
DO 224 K=l,4
L=(K-1)*5
WRITE(6,225)
225 FORMATI///59X,'RELATIVE FREQUENCIES'//)
HRITE(6,221) NOVI2)
WRITE (6,222) NOV(l), ( XX2( I «-L) ,XX2( I*1 + L) , 1=1, 5 )
HRir£«6,223) (XXK I ), XX 1( 1 + 1) , PCTKI,1+L), PCT1(I,2 + LJ, PCT1(I,3*
ID, PCT1(I,4 + L), PCT1(I,5 + L),I = 1,20)
224 CONTINUE
DO 226 K=l,4
L=(K-1)*5
WRITEI6.227)
227 FORMAT(///59X,1PROBABILITY DENSITIES'//)
WRITE(6,221) NOV12)
WRITE(6,222) N3V( 1),(XX2(I+LI,XX2(I + l + L),1 = 1*5)
WRITE(6,223) (XXI (I) ,XX 1 ( 1 + 1) , DENKI,1*L), D6NKI,2+L), DEN1(I,3 +
ID, DEN1(I,4+L>, OEN1(I,5 + D ,1=1,201
226 CONTINUE
WRITE(6,56) MC3
190
-------
IFUUNIT.E0.1.AND.I9.EQ.1) GO TO 228
229 GO TO 230
228 IFUPRSI.EQ.O) GO TO 229
19=2
IF(NOV<1).NE.3.AND.NOV(1)
1GO TO 231
U802( 1,1)=UB02«1,1)*25.4
UB32(3t 1)=UB02<3,1)*25.4
XINSZU)=XINSZ(l)*2:i.4
231 IFMOV(2).NE.3.AND.NOV(2)
1GO TO 232
UB02
-------
c***********************************************************************
c
C SUBROUTINE SORT
r
w
C IDENTIFICATION
C SORTS A REAL ARRAY A AND REARRANGES SIMULTANEOUSLY
C THE CORRESPONDING ELEMENTS OF AN ASSOCIATED REAL ARRAY 8.
C
C PURPOSE
C TO SDRT JJ ELEMENTS OF A REAL ARRAY A (BEGINNING AT All) AS
C SPECIFIED BY THE USER I IN ASCENDING ORDER.IN ADDITION,
C THE CORRESPONDING JJ ELEMENTS OF AN ASSOCIATED REAL ARRAY B
C (BEGINNING AT BID AS SPECIFIED 8Y THE USER) ARE REARRANGED
C SIMULTANEOUSLY. SORT ALLOWS SORTING UP TO 2**22-l ELEMENTS.
C
C USAGE
C CALL SORTIA(I), JJ,B(LM
C
c DES:RIPTION OF PARAMETERS
C ft(I) . -ELEMENT OF ARRAY A AT WHICH SORTING IS TO BEGIN.
C (INPUT-OUTPUT)
C JJ -NUMBER OF ELEMENTS OF A,BEGINNING AT A(I),TO BE SORTED.
C (INPUT)
C B(L) -THE JJ ELEMENTS OF B,BEGINNING AT B(L)tARE REARRANGED
C SIMULTANE3USLY.
C (INPUT-OUTPUT)
C
C REMARKS
C IF B IS AN INTEGER ARRAY THEN DELETE STATEMENT '1EAL NT,NTT'
C AND ADD A NEW STATEMENT 'INTEGER B1.
C
C***********************************************************************
C
SUBROUTINE SORT(A,JJ,B)
REAL NT,NTT
DIMENSION IU(21),IL(21),A(JJ),B(JJ)
M=l
11=1
1 = 11
J = JJ
5 IF(I.GE.J) GO TO 70
10 K=I
IJ = ( J+IJ/2
T=A(IJ)
IF(4(I).LE.T) GO TO 20
NT-B(IJ) ' ,
A(IJ)=AII)
B(IJ)=B(I)
192
-------
A(I)=T
BCI)*NT
T=M IJ)
20 L = J
IFU(J).GE.T»GO TO 40
NT=BU J)
A(IJ) = A(J)
B
-------
90 I
IFII .EO.JI GC TO 70
T=A( I + H
IF(ACI).LE.T) GO TO 90
NT=3(I*1)
K=I
100 A(K«-1)=A(K)
B(K«- 1J =B(KJ
K=K-1
IFIT.LT.AIKM GO TO 100
GO TO 90
END
194
-------
c
c
c
c
c
ft
V
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
4*
c
c
SUBROUTINE TAB1
PURPOSE
TO TABULATE FOR A GIVEN VARIABLE IN AN OBSERVATION MATRIX,
THE FREQUENCIES (NUMBER OF CBSERVATICNS) AMD PERCENT
FREQUENCIES 3VER GIVEN CLASS INTERVALS. IN ADDITION, FOP THE
SAME VARIABLE, TOTAL, MEAN, STANDARD DEVIATION, MINIMUM, AND
MAXIMUM ARE CALCULATED.
USAGE
CALL TABKA,S,NOVAR,UBO,FF.EQ,PCT,STATS,NO,NV)
DESCRIPTICN OF PAR/SMETERS
\ - INPUT MATRIX OF OBSERVATIONS, NO BY NV
S - INP'JT VECTOR SPECIFYING OBSERVATIONS TO BE CONSIDERED.
CNLY THOSE OBSERVATIONS WITH A CORRESPONDING NOM-ZERQ
S40 - NUMBER OF DBS ERVATICNS.
NV - NUMBER OF VARIABLES.
REMARKS
INTERVAL SIZE IS COMPUTED AS FOLLOWS
(uBO(3)-uBcm )/(uecm-2.)
A COUNT IS CLASSIFIED INTO A PARTICULAR INTERVAL IF THE VALUE
IS .GE. THE LOWER LIMIT OF THAT INTERVAL BUT .LT. THE UPPER
LIMIT CF THE SAMS INTERVAL.
THE DIVISOR FOR STANDARD DEVIAflON IS ONE LESS THAN THE NUMBER
195
-------
C OF OBSERVATIONS USED.
C IF S IS A MULL VECTOR, THEN TOTAL, PEAN, AND STANDARD
C DEVIATION = 0, MIN=1.E75 AND M4X=-1.E75
C SUBROUTINE TA81 IS IN IBM SYSTEM/360 SCIENTIFIC SUBROUTINE
C PACKAGE VERSION III.
C
[***<
SUBROUTINE TASK A, S, NOVAR.UBO.FREQ ,PCT , STATS ,NO,NVi
DIMENSION ACU.Sd ), U80( 1) , FREQC 1 ), PCT ( 1 ), STATS (1 )
DIMENSION WB013J
DO 5 1=1,3
5 HBO(U=UBOm
C
C CALCULATE MIN AND MAX
C
VMIN=1.0E75
VMAX=-1.0E75
IJ=NO*(NOVAR-l)
DO 30 J=1,NG
IJ = IJ*1
IFtS(J)) 10,30,10
10 IFtA(IJ)-VMIN) 15,20,20
15 VMIN=AUJ)
20 IFU(IJ)-VMAX) 30,30,25
25 VMAX=A(IJ)
30 CONTINUE
STATS(5»=VMAX
C DETERMINE LIMITS
C '
IF(UBO«1J-UBO(3) ) 40,35,^0
35 U3CI1)=VMIN
UBOl 3)=VMAX
40 INN«UBO(2)
C
c :LEAR OUTPUT AREAS
c
DO %5 1=1, INN
FREQU )=0.0
45 PCTl I)=0.0
DO 50 1=1,3
50 STATS(I)=0.0
C
C ' CALCULATE INTERVAL SIZE
C
(UBO(3)-UBO( II J/(UBO(2)-2.0)
196
-------
c
C TEST SUBSET VECTOR
C
SCNT=0.0
IJ=NO*(NOVAR-1)
00 75 J=1,NO
IJ=1J+1
IFISIJM 55,75,55
55 SCNT=SCNT+1.0
C
C DEVELOP TOTAL AND FREQUENCIES
C
STATS (1)=ST ATS C1)+A( IJ)
STATS(3)=STATS13)+A(IJ)*A(IJ)
TEM?=UBO( ll-SINT
INTX=INN-1
DO 60 1=1 ,INTX
IF«Vt IJI-TEMP) 70,60,60
60 CONTINUE
IFU( IJ)-TEMP) 75,65,65
65 FREQ( INN) =FREQ( I NN) +1 .0
GO TO 75
70 FRE3U I = FREO(I)*1.0
75 CONTINUE
IF (SCNTJ79.105, 79
C
C CALCULATE PERCENT FREQUENCIES
C
79 DO 80 1=1 , INN
80 PCTU)=FREQm*100.0/SCNT
C
C CALCULATE MEAN AND STANDARD DEVIATION
C
IF(SCNT-l.O) 85,85,90
85 STATS(2I=STATS(1)
STATS<3»=0.0
GO TO 95
90 STATSC2) =STATS(1 I/SCNT
STATS(3)=SOPT(ABS( CSTATS1 3>-STATS( 1 )*STATS( 1 ) /SCNT)/ (SCNT-1 .0) » >
95 DO 100 1=1,3
100 UBO( I)=W30tI)
105 RETJRN
END
197
-------
C *»*************************************»**********»•*#*»•****************
C
c
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
0»
w
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
w
c
c
c
c
SUBROUTINE TAB2
PURPOSE
TO PERFORM A TWO-WAY CLASSIFICATION FOF, TWO GIVEN VARIABLES IN
AN OBSERVATION MATRIX, OF FREOJENCIES (NUMBER OF OBSERVATIONS),
PERCENT FREQUENCIES, AND SOME STATISTICS OVER GIVEN CLASS
INTERVALS.
USA3E
CALL TAB21A,S,NOV,U30,FREQ,PCT,STAT1,STAT2,NO,NV)
DESCRIPTION OF PARAMETERS
4 - INPUT MATRIX OF OBSERVATIONS, NO BY NV
S - INPUT VECTOR SPEC IFlYING OBSERVATIONS TO BE CONSIDERED.
ONLY THOSE OBSERVATIONS WITH A CCRRESPONDING NON-ZERO
S(I) ARE CONSIDERED. VECTOR LENGTH IS NO
NOV - INPUT VECTOR OF LENGTH 2
NOV(1) = THE NUMBER OF THE FIRST VARIABLE
TO BE CROSS-TABULATED.
NOV(2)= THE NUMBER OF THE SECOND VARIABLE
TO BE CROSS-TABULATED.
JBO - INPUT MATRIX OF LENGTH 3 BY 2
UBOll,J»= LOWER LIMIT OF THE FIRST CLASS INTERVAL
FOR THE J TH VARIABLE, J=l,2
JBO(2,J)= NUMBER OF CLASS INTERVALS FDR THE
J TH VARIABLE, J=l,2
UBO(3,J)= UPPER LIMIT CF THE LAST CLASS INTERVAL
FOR THE J TH VARIABLE, J=l,2
IF UBO( 1,Jl=UBO(3, J), THE PROGRAM USES THE MINIMUM AND
THE MAXIMUM VALUES OF THE VARIABLE J AS UBOd.J) AND
UBO(3,J), RESPECTIVELY.
UBD(2,J) MUST INCLUDE FOR EACH VARIABLE TWO CLASS
INTERVALS FOR THE VALUES UNDER AND ABOVE LIMITS.
FREO - OUTPUT MATRIX OF TWO-WAY CLASSIFICATION OF FREQUENCIES.
ORDER 3F MATRIX IS INT1 BY INT2, WHERE INT1=UBG(2,1)
AND INT2=UBO(2,2»
?CT - OUTPUT MATRIX OF TV.C-WAY CLASSIFICATION OF PERCENT
FREQUENCIES. SAME ORDER AS FRE'J
STAT1 - OUTPUT MATRIX SUMMARIZING TOTALS, MEANS, AND STANDARD
DEVIATIONS FOR EACH CLASS INTERVAL OF VARIABLE 1
ORDER OF MATRIX IS 3 BY INT1
STAT2 - SAME AS STAT1 BUT FOR VARIABLE 2
ORDER CF M4TRIX IS 3 BY INT2
NO - NUMBER OF OBSERVATIONS.
NV - NUMBER OF VARIABLES.
REMARKS
198
-------
C INTERVAL SIZE F(H EACH VARIABLE IS COMPUTED AS FOLLOWS
C CUBOO, J»-UBO(l.J)l/(UBO(2,J)-2.)
C FO* EACH VARIABLE, A COUNT IS CLASSIFIED INTO A PARTICULAR
C INTERVAL IF THE VALUE IS . GE. THE LOWER LIMIT OF THAT INTERVAL
C BUT .LT. THE UPPER LIMIT CF THE SAME INTERVAL.
C THE DIVISOR FOR STANDARD DEVIATION IS ONE LESS THAN THE NUMBER
C OF OBSERVATIONS USED.
C IF S IS A NULL VECTOR, OUTPUT AREAS ARE SET TO ZERO.
C SUBROUTINE TAB2 IS IN IBM SYSTEM/360 SCIENTIFIC SUBROUTINE
C PACKAGE VERSION III.
C
C ********* ********* ******************************************************
C
SUBROUTINE TAB2(A,S,NOV,UBO,FREO,PCT,STAT1,STAT2,NO,NV)
DIMENSION Alll,S(l),NOV(2),UBO(3t2l,FREQ(ll,PCTIlJ .STATlll),
1STAT2(2),SINT(2)
DIMENSION HBO(3,2)
DO 5 1=1,3
DO 5 J=l,2
5 WBO( It J)=UBOtI,Jl
C
C DETERMINE LIMITS
C
DO *0 1=1,2
IFCUBGd.I )-UBC(3,I)> 40, 10, 40
10 VMIN=l.OE75
VMAX=-1.0E75
I J = NO*(NOVU )-l)
DO 35 J=1,NO
IFCS(J)) 15,35,15
15 IF(A(IJ)-VMINI 20,25,25
20 VMIN=A(IJ)
25 IF(AdJ)-VMAX) 35,35,30
30 VMAX=A( IJ)
35 CONTINUE
UBO< 1, I ) = VMIN
UBO(3,I>=VMAX
AO CONTINUE
C
C CALCULATE INTERVAL SIZE
C
45 DO 50 1=1,2
50 SINTCI ) = A3S((UBO(3,I )-UBOCl,I) )/ (UBO ( 2, I )- 2.0) )
C
C CLEAR OUTPUT AREAS ,
C
INT1=UBO(2,1J
199
-------
INT2=UBOt2,2)
INTT*INT1*INT2
DO 55 I=1,INTT
55 PCT(I)=0.0
INTY=3*INT1
DO 60 I=1,INTY
60 STAT1(I)=0.0
INTZ-3*INT2
DO 65 I=1,INTZ
65 STAT2II)=0.0
C
C TEST SUBSET VECTOR
C
SCNT=0.0
INTY=INT1-1
INTX = INT2-1
IJ=NO*JNOV<1»-1)
IJX=NO*(NQV(2)-1 »
DO 95 J=l tNO
IJX=IJX+1
IFtS(J)) 70,95,70
70 SCNT=SCNT*1.0
C
C 3ALCULATE FREQUENCIES
C
TEMPI =UBO( 1,1 >-SINT(l)
DO 75 IY=1,INTY
TEMP1 = TEMP1«-SINT (1)
IFIA1UJ-TEMP1) 80,75,75
75 CONTINUE
IY=INT1
80 IYY=3*(IY-1)*1
STATK IYYJ=STAT1(IYY)+A(IJ»
IYY=IYY*1
STATK IYY)=STAT1 (I YY) + 1.0
IYY=IYY+1
STAri(IYY» = STATl ( I YY ) +AU J)*A( I J)
TEMP2=UBO(1,2)-SINT(2)
DO 85 IX=1, INTX
TEMP2=TEMP2+SINT(2)
IF(A ( IJXJ-TEM"2) 90,85,85
85 CO NT IfJUE
I X=I NT2
90 IJF=INT1*( IX-1)*IY
FREQ(IJF)=FREQ(IJF)*1.0
IX=3*( IX-D + 1
200
-------
STAT2IIX)«STAT2(IX)+A(IJX)
IX-IX+1
STAT2(IX)=STAT2(IX)+1.0
1X=IX+1
STAT2(IX)=STAT2UX) + A(IJX)*A(IJX)
95 CONTINUE
IF (SCNT)98,151,98
C
C CALCULATE PERCENT FREQUENCIES
C
98 DO 100 I=1,INTT
100 PCTU>«=FREQ(I|*100.0/SCNT
C
C CALCULATE TOTALSt MEANS, STANDARD DEVIATIONS
C
IXY=-1
DO 120 I=1,INT1
IXY=IXY+3
ISD=IXY + 1
TEMP1-STAT1UXY)
SUM=STA,TIUXY-1)
IF(TEMPX-l.O) 120,105,110
105 STAT1(ISDI=0.0
GD TO 115
110 STAT 1 ( ISO ) =SCRT ( ABS t (ST ATI (ISD)-SUH*SUH/TEMP1)/(TEMP 1-1.0)) I
115 STAT1(IXY)=SUM/TEMP1
120 CONTINUE
IXX=-1
DO 140 I=1,INT2
lXX^IXX+3
ISO=IXX+1
TEMP2=STAT2(IXX)
SUM=STAT2«IXX-1)
IFITEMP2-1.0) 140,125,130
125 STAT2CISD)=0.0
GO TO 135
130 STAT2CISD)=SORT(ABS(CSTAT2tISD)-SUM*SUM/TEMP2)/
-------
APPENDIX B
LISTING OF COMPUTER PROGRAM FOR THE ILLINOIS
SURFACE RUNOFF MODEL
The Illinois Surface Runoff Model is programmed in Fortran IV
language for computer solutions. The input to the computer program is
the drainage basin characteristics and the rainfall hyetographs. The
output is the catchment hydrographs and pollutographs which serve as the
input to the sewer system.
The computer program allows the consideration of a maximum
number of 100 gutters at a time. Along each gutter, the subcatchments
can be approximated by as many as 10 rectangular strips. As many as
five different zones of rainfall can be considered for the entire basin.
Quality computations for two different pollutants are performed at a
time. The computations can be proceeded for as many as 100 time steps.
The storage requirement for the computer program in its present form is
400K. It can be modified to consider larger basins by changing the
arrays in DIMENSION statements if more storage is available.
A listing of the computer program of the Illinois Surface
Runoff Model is given below.
202
-------
MAIN PROGRAM FOR ILLINOIS SURFACE RUNOFF MODEL
MAIN PROGRAM PEKF^PMS THE ROUTING COMPUTATIONS FOR QUANTITY AND
QUALITY OF GUTTER FLOW
COMMnN/Zl/TEMP.lRYN»RYNST,RYNENu»DAYRYN,TM»RAlN
COMMHN/72/FFIN.F INS« C I NF . OR A I N » OPS I
COMMnN/Z3/XNU»Cl»FK»SO,ULL«NL/v»YCJlN,DT»Ql)VERL
COMMnN/Z4/OCACH.TfcA»GRtiTYP.K»GKoL,GO.T»VP»YP,-OP.FLOLAT,OPR
CtlMMnN/ZS/SG.HAFCUF.YlNG.cOTiH.SlNTH.COSTh.TANTH.vUpsT.OUPST
CUMMnN/Z20/CUNIi
*CICl,CIC2
COHMON/2B/G.L1 »TYh
COMMON/ Z°/ J TOTAL. INLET, LI Ml T.TYME NO, OTM,NSW,SWRTYP»BYSFLO
COMMON/ Z10/L2»L3»NOTM.NGTR
COMMON/21 t/Gl
COMMON/Z1 2/CJ1 «CJ2
DIMENSION CJK 100) .CJ2C 100)
RYNST(5).|iYNEND(5>)»TM(5.100)»RAIN(5,luO).nAYRYN(5)
NUKTAC5)
OCACH(100,10U),T8A(100),r,H8TYP(100)»W(loO),GR8L(100)
r,TRTYK{100),B(loO)'TlMEC100)fQUP(100,100).OPR(100),NGTR(
DIMENSION
DIMENSION
DIMENSION
DIMENSION
*100)
DIMENSION
.DIMENSION
DIMENSION
DIMENSION
DIMENSION
DIMENSION
DIMENSION
DIMENSION
DIMENSION
DIMENSION
DIMENSION
DIMENSION
SG(100)»RAFCUF(100),YING( 100)
Kl(100)'K^(100).K3(100).K<4(100),K5(100).K6flOO)
nRl(lUO),UR2(100),OH3(100),OK'i(100).OK5flon)»OR6(100)
00( II ).UM( 1 1 ). rn( 1 1 ) ' YN( 11 )
Nn(10t«plO}.FFINAL(100,10),FlNSML(inO,10).CINFLT(100.10)
OKClOO»±0),OS(loO,10)»OL(100.10)»YINO(100,10)
Ic,CINFPC 100)
INTEGER UNIT
FOLLOWING ARE THE STAT£rtt.NT FuNcTIONN DESCRIBING rROSS-SEc T I f)N AL
PROPERTIES AND THE FRICTlQN SLUpE EVALUATION FOR Thf! TYPES OF GUTTER
Al(E)=0.5*E**2/COTTH
T1
-------
T3(E)*II
R3(E)«U*E/(U+2.*E>
A3P(E)«U
T3P(E)»0.0
R3P(E)«U**2/(U»2.*E)**2
A3PP(E)=0.
R3PP(E)*il**2*t"4./(U*2.
R3PPPCF. ) = 24.*U**2/
A4(E)«ll*E
A«P(E)=U
R«P(E)«(H**2+U*riHHH)/(U*HHHH+E)**2
T«P(E)*0.0
A1PP(E)=0.n
R«PP(E)=0.0
R6PPP(De6.*(U**2*U*HHHH)/(U*HHHH + E
Gl(E)=CDF*(.67*Al(E)*RlPlE)/Hl(E)**.33+AlP(E)*Rl(r)«*.67)
G2(E)«rOF*(.67*A2(E)*H2P(E)/«2(E)**.33+A2P(E)*R2(r)**.67)
• 1(E)*R1P(E)**2/'«HE)**1.33+.67*A1(E)*R1PPCE)/R1(E)*«.33)
G?P(E)=CnF*(A2HP(L)*R2(t)**.67*1.33*A2H{E)*R2P(E)/R2{E>**.33-.22*A
*?(E)*R?P(E)**2/H2lE)**1.33+.67*A2(E)*R2PP(E)/R2(El**.33)
G3P(E)«cnr*(AJPP(L)*R3lE)**.67*1.33*A3P(E)*R3P(E)/R3tE>**.33".22*A
*3(E)*R3P(E)**2/K3(E)**1.33+.b7*A3(E)*R3PP(E)/P3(E)**.33)
GftP(E)BCOF*(A«PP{E)*Ra(L)**.67+1.33*A«P(E)*R«P(E)/R«(E)**.33-<22*A
=-.
*E)*«1.33+2.*AlP(EJ*rtIPP(El/RJ(EJ«*.33-.67*Ai(E)*RlP{E)*RlPPCE)/RH
*E)**1.33*.296*A1(L)*R1P(E)**3/R1(E)**2.3«*.67*A1(F)*R1PPP(E)/R1(F)
G2PP(E)sCOF*(2.»R^PtE)*A/pP(t)/t?2CE)**.33-.67*A2P(E)*R2P(E'**2/R?(
*E)**1.33+2.*A2PtE)*R2PP(t)/Hi:(E)**.33-.67*Az{E)*R9P(F)*R2PP(E)/R?(
G3PP(E)«COF*(2.*R3p(E)*A3pP(E)/R3(E)**.33-.67*A3P{E)*R3P(EJ**2/R3(
*E)**1.33+2.*A3P(E)*K3PP(E)/R3(E)**.33-.67*A3(E)*R-»P(F)*R3PP(E)/R3(
*E)**1.33*.?96*AJ(t)*R3PCE)**3/K3(E)**2.3<»+.67*A3(r)*R3PPP(E)/R3(E)
G«PP(E)«COF*(2.*R«PCE)*A<»pP{E)/R1(E)**.33-.67*A4P(E)*RlP(E>**2/R«(
***.33)
01(E)«CnF*Al(E)*Rl(E)**0.67
03(E)=COF*A3(E)*R3(E)**0.ft7
K( L)«H*i( E.)**u«67
r-.23«uArt*Tl(E)/OT+.5*Gl(E)*YD/DX)/(.2S*Tl(E)/DT+.5*Gl(
•-(BET-.2'j*f'AH*T2(E)/DT^.5*G2(E)*YO/OX)/(.25*T2(E)/OT+.5*G )(
*E>/OX + A1_FA)
F3(E)=F-(OFT-.2t>*&AM*T3(E)/DT+.5*G3(E)*YD/DX)/(.2S*T3(E)/nT+.5*G3(
*E)/DX+4LFA)
F«(E)«F-(BET-.2b*l-AM«T4(Ei/OT+.5*G4(E)*YD/DX)/(.2S*Ta(E)/DT*.5*G«(
*E)/OX+4LFA)
= .?5*GAM«T3P(E)/DT-.50*G3P(E)*YI)/DX
PAYA{E)«.25*GAM«T«P(E)/UT-.50*(i4P(E)*YO/DX
204
-------
«(E)/DX)
PA2(E)"(.25*T2P(E)/DT+.5*G2P(E)/DX)*{BET-.25*GAM*T2{E)/OT*.5*YD*G2
«(E)/DX)
PA3(E)«(.25*T3P(E)/DT*.5*G3P(E)/DX)*(BET-.25*GAM*T3{E)/DT*.5*YD*G3
*(E)/PX)
PA«CE)*(.25*T1P(E)/DT+.5*G«PlE>/DX>*{BET-.25*GAM*Tfl(E)/DT+.5*YD*G«
*(E)/OX)
PAYDl(De0.25*TJ(E;/DT*0.5*Gl(E)/DX + ALFA
.
F1P(E)»1.0+PAY1CEJ/PAYU1(E)+HA1(E)/PAYD1(E>**2
1.0+PAY3(E)/PAYu3lr)*PA3(L)/HAYU3(E)**2
1.0 + PAY«(E)/PAYUt(E)»PA'»(E)/PAYD«(E)**2
PAY1P(E)«0.5*YD*G1PP(E)/UX*C1.)
PAY2P(E) = 0.5*YD*Gi:PP(E)/DX*("l
PAY3P(F)=0.5*YD*GJPP(E)/OX*(-1I
PAYflP(F)=0*5*YD*G«PP(E)/DX*Cl
,25*GAM*Tl(E)/DT+.5*Yn*r,l(E)XDX)*(.2
*5*TlP{E)/DT+.5*lilP(E)/DX)*(-.2b*GAM*TlP(E)/DT+.b*YD*r,lPtE)/DX)
*5*T2P(E)/BT+.5*u2P(E>/OX)*(-.25«GAH*T2P(E)/OT*.5*vD*r.2P(E}/DX)
PA3P(E)=(*.5*G3KPtE)/DX)*(BET-.25*GAM*T3(E)/nT+.5*Yn*R3(E>/DX)*(.2
*5*T3P{E)/DT*.5*b3P(E)/OX)*(-.2b*GA«*T3P(E)/DT+.5*YD*r,3P{E)''OX)
.
PAYD1PCE)=.25*T1P(E)/01*.S*G1P(E)/DX
PAYD2P(E}=.25*T«PtE)/DT+.5*G2P(L)/Ox
rlPP(E)=(PAYlP(t)*PAY01(E)-PAYOlP(E)*PAYl(E))/PAYnl(F)**2+(PAlP(E)
**PAYDl(E)**2-2.*P«YUl(E)*PArDlP(E)*PAltE))/PAY01(F)**«
F2PP(E) = (PAY2P{L)*PAYPi!lE)-PAYUi?P(E)*PAY2(E))/PAYn2(F)**2*tPA2P(F. )
«*PAYD2(E)**2-2.*PAYlJ?(E)*PAYU2P(E)*PA2(E))/'PAYD2(F)«*«
F3PPfE)=(PAY3P(E)*PAYD3(E)-PAYD3P(E)*PAY3(E))/PAYn3(F)**2+(PA3P(E)
**PAY03(E)**2-2.*PAYU3{E)*PAY03H(E)*FA3(E})/PAY03(E)**i
F«PP(E) = (PAYlP(E)*PAYD'«(E)-PAYO
-------
6001 CONTINUE
101 FDRMAT(«I5.6F10.0)
C
C FOLLOWING IS RAINFALL DATA
C
C
C RYNST TS TIME AT WHICH RAIN STARTS
C RYNEND IS TIME AT *HlCH RAlN S'OPS
C OAYRYN IS TOTAL DAILY KAlMFALL
C NOKTA IS THE NUHBLR OF POINTS 10 DESCRIBE A MYETOr.RAPH
C TM AND RAIN RESPECTIVELY IS THE TIME AND RA INF ALL I NTENSI TV
DO 201 IRYN=I .NUZUNE
READ<5,202)RYNSr(iRYN),RYNENO(lRYN).DAYRYN(lRYN).NOKTA{IRYN)
RYNST(IRYM)=RYNST(IRYN)*60.
RYNENDCIRYN)=RYHE*D( IRYN)*60.
IF(UNIT.EQ.2)DAfRYN(lRYN)=DAYRYN(IRYN)*0. 03937
NOK=NDKTA( IRYN)
RE*D(5.203)(TM(jKYN.JJ)f KAlN(JRYN.JJ),Jj«l.NOK)
203 FORMATf 16F5.0)
DO fl91 NNOK=1 »NjK
. i )GO TU 491
NNOK) = RAlN(IRYN,NNOK)*0. 03937
CONTINUE
201 CONTINUE
C
C FOLLOWING IS GUTTER. I NLET AND &UBCATCHMENT DATA
1 = 1
C
C NGTR IS THE GUTTER NUMBER
C N IS THE NUMBER UF COMPUTATIONAL GRID POINTS FOR r.UTTER ROUTING
c IRNZDN is THE P-AIN-/ONE NUMBER THE GUTTER BELONGS TO
c GTRTYp.GL.B.s&.rt.ANo T&A ARE RESPECTIVELY THE TYPE-LFNGTH.^mTH.
C SLOPE. PEPTH AND THE ANGLE BETWEEN THE VERTICAL AND PLANE Of GUTTER
C. GRPTYP.K.GRflL.AHD OPR ARE RESPECTIVELY THE TYPE. WIDTH LENGTH AND
C OPtNING RATIO OF GRATE I*LET
C PL IS WIDTH OF STREET PAVEMENT
C YlNG IS INITIAL H^TER DEPTH lM GUTTER
C RAFCOF IS MANNING'S FRICTION FACTOR FOR GUTTER
C
36 READ(5.102)NGTR{IJ,N(I).lPNZOM(I),r.TPTYP(I),GL(l).R(T).SG(I),
*H(I).TBA(I).GReTYP(I).w{I).GKeL(I).UPR(I).PL(I).YlNG{I).RAFCOF(I)
IF(UNIT.E0.1)GU TO 600i?
GL(I)=RL(I)*3.2o
IF(GTRTYP(I).EQ.O.O)GLU J=GL(I)*3.2d
YINGCI)=YING( I )*3.28
R( I)sBC I)*3.2»
H( I)EH( I )*3.28
W( I )=W( I )*3.2»
GR6L(I )«GRPL( I )*3.20
PLU)«PLCI)*3t2o
8002 CONTINUE
102 FOR.".AT(3I5.13F5.0)
C
C FlNSG AND FFING ARE THE INITIAL AND FINAL INFILTRATION CAPACITY
C DF GUTTER SURFACE
C CIUFG TS CONSTANT OF OECAy OF INFILTRATION FOR GUTTER SURFACE
C FlNSP.FFINP AND CINFP ARE THOSE FOK STREET PAVEMFNT
C CGTR1 AND CGTK2 AKE INITIAL CONCENTRATIONS OF 1ST ANn 2ND POLLUTANTS
C IN GUTTEP
206
-------
c
READ(5.«02)FINS&C1).FF1NG(1).CINFG(1>,FINSP(I).FFINP( I).
• CGTR1CM.CGTR2C1)
«02 FDKMAT(Brb.O)
IFCUNIT.EG.l )GO TU 8003
FIfcSGU) = FINS(>U>*0. 03937
FFINGU)»FFING(U*0. 03937
FINSP( T)rFINSP( I) "0.0 39 3 7
FFINP(I)»FFINP(J)"0.0393/
8003 CONTINUE
C
C Kl.K2.K2.K3.K4.K5.K6 ARE IDENTIFICATION NUMBERS OF SIX IMMEDIATELY
C UPSTREAM INLETS
C nRl.OR2.OR3.OR1.OK5.OH6 ARE CAKKY-pVER DISTRIBUTION FACTORS
C
READ(5.502)Kl(Ij.K2(I)»K3(I).K.Kb(I).K6(I>. OR1 < I > .OR2( I )
502 FORMAT(6I5.6Fb.u)
NNN=N{ I )
DO 1 J&2.NNN
C
C nL.OS.OK- RESPECTIVELY ARE THE LENGTH SLOPE AND SURFACE ROUGHNESS
C OF SUBCATCNMENT STRIP
C FlNSHL.FFlNAL.CINFLT ARE THE INFlLTHATJQN PARAMETERS FOR sUflCATCHH
c YIKO is INITIAL DLPTH OF WATER IN SURCMCHMENT
C COl AND CO? ARE INITIAL CONCENTRATIONS OF THE 1ST AND 2ND
C POLLUTANTS IN SuBCftTCHMENT
c NO is NUMBER OF COMPUTATIONAL GRID POINTS ALONG SUBCATCHMENT STRIP
c
READ(5.103)OL(I. J).OS(l.J).OK(I.J).FJNSHL(l.J).FFlNAL(l.J)'CINFLT(
»I.J),YINn(I.J).C01(I.J).Co2(l,J).NOCl,J)
103 FORHAT(9F5.0.I5J
IF(UNIT.EQ.1)GO TU
nL
BYSFLOf JllNO = BYSFLO( JUNC )* 3 . 28*3. 26*3.28
8005 CONTINUE
202 FOKMATt 3F10.0. Ib)
703 FORMATCI5.4F5.0.1WI5)
702 CONTINUE
207
-------
DO 900 jK«lt7
00 113 I«ltITOTAL
GO TOC901»902«903»90«.V05.906.907).JK
901 IF(K1(M.EO.O)Gu '0 90b
GO TO 113
902 iFGO TO 9oe
GO TO
903 TF(K3(i).EO.O.ANO«K?(i).NE.O)GO TO 9oe
GO TO 1
90* IFtKfld
GO Tn 1
905 IF(K5(T
GO Tn i:
906 IF(K6d ).EO.O. AND.KSC I).NF..O)GO TO 908
GO TO 113
907 IF(K6d).NE.O)GU 10 906
GO TO 113
908 CONTINUE
2 TYM=0.0
LI * 0
.CO.O. AND^K3( I ).N£.0)GO TO 908
TO
L3M
GTfc=GTRTYP( I )
IFCGTR.EO.O. )00 TU 500
IF
-------
rr iFGO TO
FINS«FTHSPCI )
C1NF«CINFPCI)
IRYN«IPNZON( I>
CALL PASNlN
QL1«ORSI*PL< I)
FFIN«FFING( I >
FlNS«FIKSG< I)
CINF»CINFG(I>
IRYN»IRNZON< I>
CALL RASNIN
OL3=ORSI
GO TO 8
r IF(OK(I.jn.EO.UK(I.Jl-l).AND.Us(I.Jl).EO.OS(I»Jl-l>.AND.nUI»Jl).
*EO.OLU.J1-1))GU '0 9
GO TO F
9 IF(FFINAL(I» JD.E».FFINAL(I»J1-1).AND.FINSHL(I»J1).EQ.FINSHL(1'J1
*-l).ANn.ClNFLTCJ. Jl).EQ.CINFLT(l»Jl-l»
1. Jl)
Jl )
Jl )
IRYNalRNZONCI)
CDNINl=CPl( I. JU
COMN2eC02( 1 » JU
CALL OVLFLD
OL2=OOVERL
12 QLAT=OL1+OL2+QL3*U
FLOLAT=OLAT
YA*YO(J1-1)
YD»YN{J1-1)
GAHeYD-YA-YB
MTT = 0
YO=YD
lF(Jl.E0.2)YO=YU+tDX*XN*OLAT/SORT(100.*SS))**(3./BO
DELYO = YO-YING( I)
IF(GTR.EQ. 1 .OJGO TO 19
IF< YO.GT.HHHH)GU 10 70
RE1eOLAT-.25*T3(Yl>)*GAM/OT*.5*YD*G3(YD)/DX
ALFA=.25*T3(YD)/DT*.5*G3tYD)/OX
GO TO ?3 '
70 ALFA«.?5*Tfl( YD)/OT+.5*G
-------
ALFA«.?5*T1(YD)/DT+.5*G1(YD)/DX
BETeOLAT-.25*Tl(YD)«GAM/DT».5*YU*GKYD)/DX
GO TO 69
60 ALFA=.25*T?(YO)/OT+.5*G2(YD)/DX
69 IF( YO.GT.GOGD 10 31
CONVER«Fl*«2)
30 IF(CONVER.LT. 1 .0)00 TO 22
32 IF(MTT.NF.n)GU TO 35
YO«=YO-nELYO/20.
IF( YO.LT.YING(I))YO=YINGCI)
MT=MT+1
TF(MT.LT.20)GO TO 23
35 IF
FUNCP=r3P(YO)
GO TO 40
38 IF( YO.GT.GDJGO TO 41
FUNC=F1( YO)
rUNCP=F!P( YO)
GO TO 40
01 FUNC=F?(YO)
FU'JCP = r?P( YO)
GO TO 40
39 FUNC=F4(YO)
FUNCP=FflP( YO)
40 Yl "YO-FUMC/FUNCP
IF(ABS(Yl-YO).Lt.0.00001)GO To 42
37 YO»Y1
42 YC»Y1
«5 YN(J1)=YC
AYSE=AYSF*QL1
FATMA=FATMA+QL2
BOKl=BnKl*COVE«l«UL2
IF( Jl ,NE.NNN)GO TU 5
YP«YC
IF(GTR.E0.2.)GO TO 301
IF( YP.GT.GD)GO TO 302
OP»Ot(YC)
AsAl(YC)
VP=OP/A
T=T1(YC)
GO TO 303
302 OP=02(YC)
A«A2(YC)
VP«OP/A
T*T2(YC)
GO TO 30?
210
-------
301 JM YP.GT .HHHH)GU '0 30«
OP»03( YC)
A*A3CYC)
VP«OP/A
T«T3(YM
CO TO 303
OP=Ofl(YC)
VP=OP/A
303
C2LATspnK2/TOTA|_
CJLATeRCIK 1/TDTAL
ZlMeTOTAL/(NNN-l )
IFIL1.E0.1)CIB1=CGTR1(I>
IFCLl.f0.1)CIB2»C6TR2(I>
C1C1=(C1LAT*CIB1*A/DT*CUPST1*OP/GL(I))/(A/DT+QP/GL(I)+ZIM)
CIB1=CJC1
CIB2=CIC?
3303 CALL DPHBQU
5 CONTINUE "
DO 55 J2=1.NNN
55 YO(J2)=YN( J2J
2678 IF(TYM.LT.TYMLNu)GO TO 500
113 CONTINUE
900 CONTINDE
CALL SHRJNT
STOP
f ND
SUBROUTINE UPSBO
C THIS SUBROUTINE COMPUTES uPSTRtAH BOUNDARY CONDITION AT EACH
C TIME LEVEL
COMMON/Z5/SG.RAFCOF.YING.COTTH'SINTH.COSTH.TANTH.YUPST.QUPST
COMMON /26/Kl»K2.K3tK«»K6»ORl»OH2»UR3»OR<»»UR5»OR6.K5
COHMnN/27/I»GlRTYP.B.TlHE,OUP,^(iTRl.CGTR2.CUPSTl.cUpST2,CCCl.CCC2.
•C1C1.CTC2
COMHON/Zfl/G«Ll • TYM
DIMENSION Sf.(100)»RAFCOF(100).YING(100)
DIMENSION KHlOO)»K2(100).K3tlOO).xa(100)»K5(100).K6(100)
DIMENSION DC1(6).UC2(6)
DIMENSION nRl(100).OR2(100),CR3(100),OR1(100)»OR5(100)»OR6(100)
DIMENSION UPINF(6).CGTR1(100).CGTR2(100)
OIHENSIPN CCCl(10U.100).CCC2(100.100)
DIMENSION GTRTYP(100),H(loO),TlME(100).OUP(100,100).nPR(100).NGTR(
*100)
DO 1 MM=1.6
TF(MM.E0.2)MIN=K2(
IF(MH.E0.3)MIN=K3(
IF(MM.F0.4)MIN=Kl(
IF(MH.F0.6>MINsn6( I )
IF(MIN.FO.O)UPI.NF(MM)=0.0
IF(MIN.EC.O)DCl(Mh)=CGTRl{I)
IF(M1H.EO.O)DC2(MM)=CGTR2(I)
IF(MIN.EO.O)GO TO 1
MALe?
69 IF(1 YH-TIME(MAL) )1 12p 1 11»1 10
211
-------
110 MAL=MAL+1
GO TO f,9
111 UPlNF(MM)sQUP(HjN.MAL)
DC1(MM)=CCC1 (MIN.HAU)
nC2(HH)sCCC?(rt!.'*.MAL)
GO TO 1
112 |IPlNF(MM) = OUP(MlN'MAL-l) +
*TIME(MAL-l))*(TrM-TIME(MAL-l))
RATln=(TYM-TIME(MAL-l)>/(TlML(MAL)-TlMECMAL-l)>
DC2(MM)sCCC2(MlH»MAL-l)+(cCC2(MlN.MAL)-CCC2(MlN,HflL-l))*HftTlO
1 CONTINUE
DISCH = M.O)*(UR1(I)*UPINF(1)»OR2(I)*UPINF(2)+CIR3{I)*IIPINF(3)*ORO(I
*)*UPINrC«) + fJKb( I)*UPINF(5)+OK6(I)*l)PINF(6))
IFCDISCH.EO.O. )CUPSTl=CGlRl(I)
IFCDISCH.EQ.O.KUPST2 = CGlR2(l)
IFCDISCH.EQ.0.0 JGU TO 200
DlSJ=npl(I)*UPlHF(l)»DCl(l)+l.lR2(l)*uPINF(2)*nCH?i«-OR3(I)*UPINF(3)
*(6)*DC1 (6)
DIS2aOPl(I)*UPIiGR8TYP.w>GRdL.GO»T»vP»YP.OP»F|OLAT»OPR
COMMON/Z7/I.GTRTYP,d»TlME,QuP.CGTRl.CGTR2»CuPSTl,rUPST2.CCCl.CCC2i
*C1C1,CIC?
COMMON/Z10/L2.L3.NOTM.NG1R
N' cr.TRH 1UO).CGTH?( 100)
N Crd ( 100. 100) .CCC21 100» 100)
N GTRTYP(100),B(loO)'TlME()00).QUP<100,109).nPROOO),NGTR(
*100)
DIMENSION QCACH(100.100)tT8A(100).RR8TYP(100).W(lnO),r,R8L(100)
TIHE(1)=0.
OUP( I . 1 )=O.C
IF(GR8TYP( I ) .LO.U. )GO TO 20
IF(GR8TYP(I).E0.1.0.0R.GR«TYP(I).E0.2.)COF=1.0
IF(r,P8TYPC I ) .LO. 3.0)COr = 0. «
IF(GR8TYPtI).EO.t.n)COF=0.8
lF(f.R8TYP(I).E0.5.)COF = 0.7
IF(GTRT YP( I ) .E.0.2. )GO TU \ t
DEH^WC T )/TAf.'( T«A( 1 ) )
lF(YP.r,T.DFP)l'0 TU 13
KLET*OPR(I)/SIN(TtiA(I))
YM^YP/?.
212
-------
GO TO 12
13 HI«*C
YM«YP-nEP/2.
GO TO 12
11 WL=w(I)»OPR( I )
YM»YP
12 HFLe3.*COF*h'L*(YM«-VP*VP/(2.*G))*M.5
IF(NFL.LE.CP)GO TU 11
HFLs0.6*WL*GH6L(I>*SQRT
-------
c
C THIS SUBROUTINE COMPUTES THE QUANTITY AND QUALITY OF FLOW IN
C SUBCATCHMENT STrtlPS
C
COMMON/Zl/TEMP.IRYN.RYNST.RYNENO.OAYRYN.TM.RAlN
CDHMON/Z3/XNU»Cl»FK»SO«OLL»NUv»rOINtDT»QOVERL
CDHMDN/ZB/r,,Ll • TYM
DIMENSION YOLD(n>.YNEn(ll),00LD(m.QNEN(ll)
DIMENSION PYNST(5).r
GL(E)=3.*ALF*E**2
GLP(F)=6.*ALF*E
FUNClLfE)=YD*GLP(E)/(GL(E)+GAM/HAT )*(•!.)
rUNCPLfE)=(YD*GLPK(E)*((iL(E)+GAH/«AT)-YO*GUP(E)**?)/(GL(E)+GAM/RAT
1 )**2*(-l . )
PAY{E)=RAT*GLP(E)*(HET+YD*RAT*GL(E))
PAYDP(E)s2.*«AT*GLP(E)*(bAM + RAT
rLP(E)=1.0+fUNClL(E)*PAY(E)/KAYD(E)
=FUKiCPL(E) + (PAYP(E)*PAYU(t.)-PAYOP(E)*PAY(E>)/(PAYD(E})**2
E)=ALF*F**3
..
GTPP(E)=SQRT(b.*G«SU/E**3)*(-0«75*ALOGlO(2.*E/FK)-O.B5)
FT(E)=F"(BFT+KAT*TO«GT(E))/(GAM+HAT*GT(E))
.
TUNCPT(E)=( YD*GTPH(E)*(GT( E)+GAM/RAT )-YD*GTp(E )**?)/( GT(E)+GAM/RAT
1)**2*(-1 . )
PAT(E)=RAT*GTP(E)*(8ET+YD*RAT*GT(E))
PATP(E)=RAT*(BE1+RA1*YD*GT(E))*GTPP(E)+YO*(RAT*GTP(E))**2
/(PATD(E))**2
PATOP(F) = 2.*«AT»&TPCE)*(CjArt + KAT*GT(E))
FTP(E)=1.0*FUNClT/PATu(E>
FTPP{E)=FUMCPT(L)+(PATP(E)*PATU(E)-PATOP(E)*PAT{E))
OT(E)=SOPTf3.*G*SU*E**3)*(2.»ALUGlO(2.*E/FK)+1.7«)
REYCR(E)=Cl*(2.*ALUblOl2.*E/FK)*1.7q)**2
IF(OLL.EO.O.O)QJVLHL=0.0
=.
IF(OLL.EO.O.O)QJVLHL=0.0
IF(OLL.EO.O.O)Gu '036
CAUL RASNIN
OOL=ORSI
ALF=B.*G*SO/(C1*XNU)
RA7=OT/OOX
IF(Ll.NE.nGn Tu 2
FLOw=l.
DO 1 lnv=l .NOV
1 YOLDC lnv>=YOlN
2 no 3 KOV=I .NOV
IF(KOV.FO. 1 )GO TO 4
YOO»YNFW( KOV-1
YOB=YnLD(KOV)
IF(FLOV.'.F0.2.0.ANt).YOB.GT.FK.ANO.YOD.GT.FK)Gn TO 3 1
BET = 2.*oOL*DT-(lfOLi-YOA-YOci)+RAT*Yon*GL(YOU)
GAH=1 . + RAT«GL( YUO)
MOT = 0
MOTT = 0
214
-------
YO=«ONFW(KOV-1)4UOL*DOX)/ALF)**0.33
Dl Y0=( YO-YOIN)/«:0.
D?Y() = YO/10.
YYO=YO
CDNVFR = FL(YO)*Fl_PH'(YO)/(FLP] YO
MOTT=MnTT+1
IF(MnTT.|_E.20)GU <0 7
IF(MOT.EO.0)YO=YYO
MOT=MOT+1
IF(MOT.LE.75)GO TO 7
GO TO 31
5 HO 9 Mj =1 ,50
YO=ftRS( YO)
Yl=YO-FL( YO)/FLH( YO)
IF( ABS( YO-Yl) .Lt. 0.00001 )r,0 TO 10
9 YO=Y1
10 YQC=Y1
IF(Yf)C.L'T.YOIN)YOC =
LD NFW(KOV)/XNU
JF(REYNO.LT.REYtK( YOO.OR. YOU.LT.FK)GO TO 3
lF(REYNO.LT.RtYCH( YOO.UR. YUti.LT.FKJGO TO 3
FLOwr?.
31 BET = ?,n*QO|_*OT-(YUD-YOA-YnB) + RAT*YOD*GT
GAM=1 ,+RAT*GT( YUO)
MAT = 0
MATT = 0
YO=YOB400L*DT
01YO=( YO-YOIN)/20.
D2YO=YO/10.
YYO=YO
37 CONVER=FT(YO)*FTPP(YO)/(FTP(YO)**2)
CONVFR=ABS(CONVER)
IFCCOMVER.LT.l. JGO TO 35
IFCMAT.NE.OJGO TO 36
YO=YO-ni YO
IF(YO.LT.FK)YO=FK
HATT=MATT+1
IF(KATT.LE.203GJ lQ 37
36 IF(MAT.EO.O)YO=TYO
YO=YO+n?YO
HAT=MAT+1
IF(HAT,LE. 75)60 TO 37
WRITF(f.P)
8 FORMATC2X . 'OVERLAND ITERATION FAILS1)
STOP
35 00 39 M?=1.50
YO=ABSfYO)
Y!=YO-rTC YO/FTPC Yn)
IF( ABSC YO-Y1 ) .LL. 0.00001 )QO TO <<0
39 YO=Y1
«0 YOC=Y1
IF(YnC.LT.YOIN)rOC=YOIN
YNEWCKPV)=YOC
GO TO 3
215
-------
« YNEW(l)sYniU
QNE«(1)«0.0
3 CONTINUE
If (LI .F.0.1 )CONH1=CUNIN1
QOVERLrQNEwCNUV )
*/OLL+QOL )
COVrR2c(CONB2*rnE«(NOV)/DT+CONlN2*QOVERL/ULL>/(YNrK.OUP(100,100)
DIMENSION CGTftli 100) .CGTn?( 100)
DIMENSION CCC1 ( 100, 100),CCC2( 100. inO)
DIMENSION GL(100).W<10U>.NGTft(100),(iR8L(100).OPR(lOo)
IF(L1.NE.1)GO TU 1
AREA=W(I)*RR6LCj}*UPR(I)
XK=-(0.80*AREA*sQ«T(2.*G))/2.
XC1=GL( I)/OT
XI1=0.
XH1=0.
X01=0.
1 XI2=OP
XC2s.5*(XIl*XI2*XHl)+XCl*XHl
DISC=XK**2+«.0*^C1*XC2
XHO = (4-XK + SQRT(OISC))/C2.*XC1)
IFCXHO.LT.O.O)XHO=0.
XH2=XHO**2
X02=-2.*XK*XHO
IFCL2.NE.O)GO TU 2
L2=MDTH
L3=L3+1
OCACH{NGTR( I ) .L3) = Xii2
IFCX02.EO.O.)CCC1(NGTR(I),L3)=0.
IF(X(32.EO.O.)CCC2(NGTR(I),L3)=0.
1F( XQ2.EO.O..)Gri TU 2
XQ1=XQ?
XH1=XH?
XI1=XI2
RETURN
SUBROUTINE S«RJNT
216
-------
c
C THIS SUBROUTINE COMPUTES THE QUANTITY AND QUALITY OF DlREcT
C TO SEWER NODES
C
COMMnN/Za/QCACH.TbA»GR6TYP.ri»GR6L.GO«T»VP.Yp,UP.FLOLAT.OPR
COHMON/Z7/IiGT«IYP.B.TlME,OUP.talRl,CGTR2»CUPSTl.cUPST2»CCCl,CCC2.
*CIC1.CIC2
COMMON/ZB/G.L1•TYM
CfiMMnN/Z9/jTUTAL.I'''
DIMENSION CCrl(100,100).CCC2(100.100)
DIMENSION GTPTYtJ('00).B(loO).TlME(100).QUP(100,100).OPR(100),NGTR(
*100)
IHD=TYMEND
IDTM=DTM
ISTP=IHD+IOTM
DO 1 JJ=1,JTOTAL
,
200 FORMAT(2X> '*************•***************************************• )
WKlTE(6.eO)NSW( JJ)
80 FORMAT f 1 X . TLOW IMTU SErit-R JUNC T 1 ON= ', I 5 , i ********************* • )
WRITE(6,?00)
WRITE(ft.81 )
81 FORMATf///2X» 'TIME1, 3X. 'DISCHARGE'. 3X. '1ST POLLUTANT ', 3X .' 2ND POLL
• U T A N T • )
WRITEC A.P6)
86 FORMAT( IX.'C SEC )'»6X.'(CFS)'»6X,'CONCENTRAT ION', 3x» 'CONCENTRATION'
*)
TIME(l)sO.
ODIS( 1 )=RYSFLO( JJ)
COB1 ( 1 ) = r Jl ( JJ)
COB2( 1 )=CJ?( JJ>
DO 2 Nto = ?. LIMIT
ITYM(NH)=TIME(NN)
6 DIS=0.
CIS1=0.
CIS2=0.
DO 3 MM=I , 10
K=lNLET( JJ.MM)
IF(K.EO.O)GO TO 15
CIS1=CIS1+OCACH(K.NN)*CCC1(K»NN)
CIS2=CIS?+OCACH(K.NN)*CCC2(K»NN)
3 DlSsDI S + OCACHC K.NN )
15 ODIS(NN)=DIS+BY6FLO( JJ)
) = CC]Sl+uYSFLn
-------
22 KRITE(7.?5)(lTYM(MH),OOIS(MM),MMal,NL)
25 FORMATtI5»lX.F5.2»lX.I5.1X'F5.<:.lX»I5tlX»F5.2.lX.i5,lX.F5.2,lX»I5«
*lX»F5.?.lX,l5tlX.K5.2)
IF(SWRTVP(JJ).Eg.2.0}GO TO 3C
31 FORMflTC ' JENO' )
GO TO 100
30 WR1TEC7.32)
32 FDKMATC 'FEND* )
100 DO 82 MAMAsl.LIrtIT
WRlTE(A,83)TIME(M«MA).aUls(MAMA).COBHMAMA),COB2(MAHA)
83 FORMAT(F7.2.2X.F6.2»9X.F6.2.9x»F6.2)
82 CONTINUE
1 CONTINUE
PETURN
END
218
-------
APPENDIX C
LISTING OF COMPUTER PROGRAM FOR THE ILLINOIS SEWER
SYSTEM WATER QUALITY MODEL
The Illinois Sewer System Water Quality Model is programmed in
Fortran IV language. The input to the computer program includes the
depth and discharge hydrographs at the entrance and exit of each sewer
and the volume of water at each storage junction at given times as provided
by the output of the ISS model. In addition, the direct inflow hydrographs
and pollutographs to the sewer junctions as obtained from surface runoff
computations and the sewer system layout are also input to the program for
runoff quality routing in the sewer system. The output from the computer
program are the pollutographs at the sewer system outlet, at the storage
junctions and at the entrance and exit of all sewers.
The computer program allows the consideration of as many as 100
sewers at a time. The runoff quality routing is performed for two different
pollutants. The computations can be proceeded for as many as 100 time
steps. The storage requirement for the computer program in its present
form is 300K. It can be modified to consider larger sewer systems by
changing the arrays in DIMENSION statements if more storage is available.
A listing of the computer program for the Illinois Sewer System
Water Quality Model is given below.
219
-------
c
C SEWER SYSTEM WATER QUALITY MODEL
C
DIMENSION Q01llOO>.Q02<100>.UIl(100).Ol2C10o)
DIMENSION OOUTH100.100).QOul2C100.100).OINl(100.lOOJ.OIN2tlOO»100
o
DIMENSION TYM(100),QU(100).HU(100).DD(100)»Hn(100)
DIMENSION OOU(luO).UQDC100).HHU(100),HHO(100)»AAU(100>»AAn(100)
DIMENSION ON(10U)»HN(100).NOUE(100)
DIMENSION AREA(100),C1IN(100),C2INUOO>
DIMENSION NOOEUP(*00).NODOrtN(100),PLENGT(100)»D(loO).clU(lOO)
DIMENSION C1D(100).CN1(100).CN2(100)
DIMENSION C2U(luO>.C2DUOO).Y(100).TM(100)
C
C FOLLOWING ARE ThE GENERAL INPUT PARAMETERS
C
C NPIPE. TOTAL NUMBER OF SEHERS IN THE SYSTEM
C NSTJUN. TOTAL NU.tBLR OF REsER VOl R-T YPE JUNCTIONS IN THE SYSTEM
C NPYP08. NUMBER OF POINTS USED TO DESCRIBE EACH OF THp DISCHARGE AND
C STAGE HYURUGRAPHS AT THE ENTRANCE AND EXIT Or EACH SEWER
C NJNCD8. NUMBER OF POINTS USED TO DESCRIBE DIRECT INFLOW HYD«OGRAPHS
C AT EACH UF THE SEKER NODES
C NROOT, IDENTIFICATION NUMBER UF THE SEwER NODE AT THE SYSTEM OUTLET
C DT.THE TIME INTERVAL AT WHICH GRDlNATES OF THE DIRECT INFLOW
C HYDROGRAPHS AT THE SE«ER NUDES AHE PROVIDED
C
READ(5.1)NPIPE»NSTjUN,NPYpDtt»NJNCDB»NROOT.DT
1 FORMAT(5I5.F5.0J
C
C FOLLOWING ARE TnE INPUT PARAMETERS FOR EACH SEWER
C NODEUP(I). IDENTIFICATION NO OF JUNCTION NODE UPSTpEAM OF sEwER I
C NDDUWN(I).IDENTlFlCATI&N NO UF JUNCTION NODE DOWNSTREAM OF SEWER I
C PLF.NGT(I). LENGTH OF SEwER I
C D(D. DIAMETER OF SEwER I
C ClUU ).. INITIAL CONCENTRATION OF THE FIRST POLLUTANT AT ENTRANCE OF I
C ClOCI ). INITIAL CONCENTRATION OF THE FIRST PULLUTAwT AT EXIT OF 1
C C?U( I ). INITIAL CONCENTRATION OF THE SECOND POLLUTANT ftT ENTRANCE OF I
C C2D(I), INITIAL CONcENThATtUN OF TnE SECOND POLLUTANT AT ExlT OF I
C TYM(J1).T1ME AT Jl'ST TIME LEVEL
c out ji). DISCHARGE AT ENTRANCE AT JI»ST TIME LEVEL
C HU(J1).FLOK DEP1H AT ENTKANCE AT Jl'ST TIME LEVEL
C OD(Jl). DISCHARGE AT EXIT AT Jl'ST TIME LEVEL
C HD(J1).FLOW DEPTH AT EXIT AT Jl'ST TIME LEVEL
C
DO 2 !»!• NPIPE
READ(5.3)NODEUP(I),NODUWN(l),pLENGT(I)»D(I).ClU(I>.C10f I J
IFtClDfI).EO.Clu(I))ClD(I)=1.05*ClU(I)
IF73)
73 rORM«T(!X.«HATEK OuALITY CONDITIONS AT SE«ER SYSTF.M OUTLET****** ' )
74 FORMAT ( IX. * ft***************************************************** )
HRITE(6.75)
75 FORMATt////"X. 'TIME'. 5X. 'DISCHARGE1. 3X. '1ST POLLUTANT «. 3X, ' 2ND POL
*LUTANT' )
220
-------
NRITE(6,76>
76 FORMAT(3X.•(SEC)'.6X,'(CFs)''6X,'CONCENTRATION',3X.'CONCENTRATION'
*)
GD TP 900
66 FORMATfF7.2»5X.FB.2.9X.F6.2.9x»F6.2)
230 WRITE(6.76X,'(CFs)'»6X.'CONCENTRATION', 3-^. t CONCENTRATION'
*p7x,«(CFS>'•ibXt'CONCENTRATION'.ax."CONCENTRATION•>
900 DO 6 LIN=1.1
GO T0(7.8,9>10),LiN
7 DO 11 Klel.NPYPUfl
11 Y(K1)=OU(K1)
GO TO 15
8 DO 12 Kl=l,NPYPu8
12 Y(Kl )-=HU(Kl )
GO TO 15
9 00 13 K1=1,NPYPU8
13 Y(K1)=OD(K1)
GO TO 15
10 DO 14 K1=1.NPYPD8
Ifl Y(K1)=HD(K1)
15 K=l
TE = OT
17 K=K*1
J2 = 2
69 IF(TE-TYM(J2))112»111.110
110 J2=J2+1
GO TO 69
111 YE=Y(J?)
GO TO 200
112 YE = Y(J2-1) + CY-Y(J2-1).)/(TYM(J2)-TYM< J2-n>*
AlU=.2^*nd)*0(I)*ATAN(ARr,U)
A2U=(.5*nd)-HHU(KK))*SuRT(0(I)*HHU(KK)-HHU(KK)*HHU(KK))
AAU(KK)=AlU"A2U
ARGD=SQRTCD(I)*HHO(KK)-HHD(KK)*HHD(KK))/(.5*D(I)-HHDfKK))
221
-------
A1D«.25*0(I)*D(I>«ATAN
A2D=(.5*DCI)-HHU(KK))*SOKT(DU>*HHD-HHD{KK)*HHD{KK)>
AAD(KK>BA1D-A2D
80 CONTINUE
DO an L«?»KKK
AA«AAU(L-1)
OB-ttOD(L-l )
AB»AAD(L-1)
ODDmQQtl(L)
AD=AAU(L )
OCBOQO(L)
ACBAADCL)
DO 31 LL=li2
GO TOCfll ,12) .LL
41 lf
B1«-OC/PLENGT(I)
B2=AD/nT-QDD/PLENGT(I)
C«AC*CR/DT
CC=(C*R2-E*B1)/(X1*B2-X2*B1)
CD=(Xl*F-C»X2)/(Xl*d2-X2*Bl)
If (CC.LT.O.)CC=0.
If
-------
533 FORHAT(2F10.2.2X.Flo.2.9X,FlO,2.15X,F10.2.2X.F10.?»9X.F10.2)
30 CONTINUE
2 CONTINUE
DO 50 JJel.NSTJuN
C
C FOLLOWING ARE THE INPUT PARAMETERS FOR EACH RESERyOl R-T YPE JUNCTION
C
C NODE(JJ). IDENTIFICATION NijHBER OF THE JUNCTION JJ
C AREA(JJ). CROSS SECTIONAL AREA OF THE JUNCTION JJ
C CllN(JJ), INITIAL CONCENTRATION OF THE FIRST POLLUTANT IN JUNCTION JJ
C C2lN(JJ). INITIAL CONCENTRATION OF THE SECOND POLLllTANT
C TYM(J3).TIME AT J3'RD TIME LEVEL
C ON(J3).RATE OF UlKECT INFLOw AT J3'RD TIME LEVEL
C CNKJ3), CONCENTRATION OF THE FlKST POLLUTANT OF DIRECT INFLOW
C CN2(J3»,CONCENTKATIUN OF THE SECOND POLLUTANT OF nlRF.CT INFLOW
C HNCJ3). WATER STAGE IN JUHCTIUN AT J3'RD TIME LEVEL
C
READC5.151>NODECJJ).AREA(jJ).ClIN(JJ).C2IN{JJ>
151 FDRMAT(I5.3F5.0)
RE ADC 5. 152)(TYM( J3).GN( J3).CNl( J3).CN2( J3),HN( J3)» J3«l'NjNCD°>
152 FORMAT(F6.0»4F5.0«F6.0»4F5.0»F6.0.0F5.0)
WRITEtf ,54)NOOE( JJ)
54 FORHATMXt 'CONDITIONS AT JUNCTION NU« •• I 5. '•**•**»******•***••**')
WRITEC6.74)
WR1TE(A.55)
55 FDRMAT(///2Xt'TlMt'.6X» 'STORAGE*. 4X('1ST POLLUTANT '. 3X. '2ND POLLUT
«ANT»)
WRITE(6.56)
56 FORMAT( IX. '( SEC J'»6X.'(FT3)''6X.'CONCENTRATION'.3X' 'CONCENTRATION'
O
IF(KKK.GT.NJNCDU)KKK=NJNCD8
DO 90 11*1. KKK
002(II)=0.
OI1(II)=0.
90 Ol?CII)=0.
DO 91 HMel.NPlPE
IF(NODE( JJ).NE.NODEUP(MM))GO TO 92
DO 93 L1=2.KKK
ODl(Ll)*001(Ll>*aflUTl(MH»Ll)
93 002(L1 )=002(Ll)*OUUT2(MM.Ll)
92 IF(NODECJJ).NE.NOOOHN(MM))GO TO 91
DO 94 L2»2.KKR
OI1(L2»«OI1(L2)-»QIN1(MM.L2)
94 QI2(L2)=OI2(L2)+01N2(MM,L2)
91 CONTINUE
DO 57 K=?»KKK
IF(K.E0.2)C10LD«C1IN(JJ)
IFON(K)*CN2(K)"002CKl
HNfW=HNCK)
HOLD=HN(K-1 )
CNEW1=(C10LD*HOLD*BET1«DT/AREA(JJ))/HNE»«
CNEH2=(C20LO*HnLO*BET2*OT/AREA(JJ))/HNEw
|F(CNEKJ ,LT.O.)CNErf!=0.
IF(CNEW?.LT.O. )CNEW2=0.
S«HNEW4AREA( JJ)
WRITE(6.66)TYH(K)>S'CNEW1.CNEW2
C10LD»CNEW1
C20LD«=CNEW2
223
-------
57 CONTINUE:
50 CONTINUE
STOP
f NO
224
------- |