United States Office of Air Quality EPA-450/4-86-003
Environmental Protection Planning and Standards January 1986
Agency Research Triangle Park NC 27711
Air
Continued Analysis
And Derivation Of
A Method To
Model Pit Retention
-------
-------
EPA-450/4-86-003
CONTINUED ANALYSIS AND
DERIVATION OF A METHOD TO MODEL
PIT RETENTION
By
K.D. Winges
C.F.Cole
TRC Environmental Consultants, Inc.
Englewood, Colorado 80112
Contract No. 68-02-3886
EPA Project Officers:
J.L Dicke
J.S. Touma
U.S. ENVIRONMENTAL PROTECTION AGENCY
Office of Air and Radiation
Office of Air Quality Planning and Standards
Research Triangle Park, N.C. 27711
January 1986
-------
DISCLAIMER
This report has been reviewed by the Office of Air Quality Planning and Standards, U.S. Environmental
Protection Agency, and approved for publication as received from TRC Environmental Consultants, Inc.
Approval does not signify that the contents necessarily reflect the views and policies of the U.S. Environmental
Protection Agency, nor does mention of trade names for commercial products constitute endorsement or
recommendation for use. Copies of this report are available from the National Technical Information Service,
5285 Port Royal Road, Springfield, Virginia 22161.
-------
TABLE OF CONTENTS
SECTION TITLE PAGE
1.0 SUMMARY AND PURPOSE 1
2.0 BACKGROUND 5
2.1 SUMMARY OF PREVIOUS DATA ANALYSIS 6
2.2 OVERVIEW OF CURRENT DATA ANALYSIS 9
3.0 SIGMA-THETA DATA ANALYSIS 11
3.1 SIGMA-THETA AUDIT 11
3.2 DATA AVERAGING 13
3.3 METEOROLOGICAL DATA ANALYSIS 15
4.0 ANALYSIS OF ESCAPE FRACTION EQUATIONS 23
4.1 DERIVATION OF THE CANDIDATE ESCAPE FRACTION EQUATIONS. . . 24
4.1.1 THE ORIGINAL WINCES EQUATION 25
4.1.2 ALTERNATIVE 1 CONSTANT-K USING LINEAR MODEL 29
4.1.3 ALTERNATIVE 2 CONSTANT-K USING A MORE DETAILED MODEL. 32
4.1.4 ALTERNATIVE 3 VARIABLE-K USING LINEAR MODEL 35
4.1.5 ALTERNATIVE 4 VARIABLE-K USING THE MORE
DETAILED MODEL 38
4.2 EVALUATION OF CANDIDATE EQUATIONS WITH EXPERIMENTAL DATA . 39
4.2.1 EVALUATION DATA 39
4.2.2 COMPARISON OF THE CANDIDATE EQUATIONS WITH THE
EVALUATION DATA -. 44
5.0 ESCAPE FRACTION ALGORITHM FOR ISC 51
6.0 CONCLUSIONS AND RECOMMENDATIONS 53
6.1 CONCLUSIONS 53
6.2 RECOMMENDATIONS FOR FUTURE WORK 56
REFERENCES 61
APPENDIX A - AIR SCIENCES AUDIT REPORT A-l
APPENDIX B - HOURLY METEOROLOGICAL DATA BASE B-l
APPENDIX C - FORTRAN LISTING OF MODIFIED ISCST PROGRAM . . . . C-l
APPENDIX D - TEST RUNS AND SAMPLE INPUT FILE D-l
iii
-------
-------
1.0 SUMMARY AND PURPOSE
This report continues the analysis of pit retention meteorology and
predictive escape fraction equations begun in EPA's "Dispersion of Airborne
Particulates in Surface Coal Mines" (EPA, 1985). The purpose of this work,
which is described in this report, was three-fold:
Examine the existing meteorological and smoke release data
base to determine the relationship between in-pit and
out-of-pit sigma-theta and alphabetic stability class in
order to identify trends and other systematic behavior.
Incorporate other physical or meteorological parameters
(particularly wind speed) into the original Winges escape
fraction equation. Refinements to the basic equation are
to be tested against the existing field data.
Prepare and document a computer algorithm to predict escape
fraction for use in the ISC model.
The analysis of the meteorological data in and out of the pit yields an
important finding: the sigma-theta (standard deviation of horizontal wind
direction) inside the pit is almost always greater than the sigma-theta value
measured simultaneously outside the pit. This indicates that the horizontal
turbulence in the pit is greater than outside, and it is suspected that the
enhanced in-pit sigma-thetas are induced by mechanical turbulence as air
passes over, and in the wake of, the mine pit wall. The degree to which the
in-pit sigma-theta exceeds that out-of-pit^ ' inci
is not related to Pasquill-Gifford stability class.
in-pit sigma-theta exceeds that out-of-pit increases with wind speed, but
Both the in-pit and out-of-pit sigma-thetas appear to provide a
reasonably good measure of alphabetic stability class, when computed over a
one-hour time period. The alphabetic stability classes measured in and out of
the mine pits are identical to, or only one class removed from, the
Pasquill-Gifford stability class for roughly 80% of the data base hours.
1. As measured by the ratio of out-of-pit sigma theta divided by in-pit
sigma-theta.
-1-
-------
In an effort to incorporate other physical and meteorological parameters
(especially wind speed) into the original Winges escape fraction equation,
four alternative modifications to the Winges equation were derived. The
alternative escape fraction equations differ in simplifying assumptions and in
complexity:
ALTERNATIVE 1: CONSTANT-K LINEAR MODEL. The derivation of this
equation assumes a constant value of eddy diffusivity with
pit depth, and assumes that eddy diffusivity varies linearly
with wind speed.
ALTERNATIVE 2: CONSTANT-K DETAILED MODEL. Like the previous
derivation, the alternative 2 escape fraction equation assumes
that vertical diffusivity is constant with pit depth. However
the influence of both wind speed and stability class on dif-
fusivity is taken into account by introducing the Monin-Obukhov
length as a measure of stability.
ALTERNATIVE 3: VARIABLE-K LINEAR MODEL. The derivation of this
equation recognizes that eddy diffusivity is not constant with
pit depth.
ALTERNATIVE 4: VARIABLE-K DETAILED MODEL. The most complicated
of the four alternatives, this derivation uses variable eddy
diffusivity with pit depth, and incorporates Monin-Obukhov
length as a measure of stability. An involved numerical solu-
tion is required to compute escape fraction with this alter-
native.
The four alternative escape fraction equations were evaluated by
comparing values of escape fraction computed from the alternative equations
with values of escape fraction inferred from the smoke release data. In
general, the alternative equations predicted smaller escape fractions than did
the original Winges equation. Furthermore, all of the alternative equations
exhibit a much greater change in escape fraction with wind speed than does the
original Winges equation, and the increase in predicted escape fractions with
wind speed matches the trend observed in the smoke release data. In this
sense, the introduction of wind speed into the Winges equation is successful.
-2-
-------
However, the overall conclusion drawn from examining all of the
alternative equations' predicted escape fractions is that they do not perform
as well as would be liked. The correlation coefficients between predicted
escape fractions and those inferred from the smoke release data are never
greater than 0.39, and attempts at optimizing the agreement by introducing
linear coefficients into the alternative escape fraction equations show very
little improvement. Discrepancies between analytically predicted escape
fractions and those inferred from the smoke release data are attributed to two
factors. First, it must be remembered that the smoke release data do not
provide a direct measure of escape fraction, and it is possible that some
differences in measured and predicted escape fractions are due to
misinterpretation of the smoke data. Second, the original Winges equation,
and all of the alternative equations, assume that dust is removed from the
mine pits by dispersion rather than by convection. This suggests that the
Winges equations may be better predictors of escape fraction during stable
conditions than during unstable or neutral conditions. A re-examination (and
possibly re-interpretation) of the smoke release data gathered during stable
conditions may be warranted, particularly since it is the stable atmospheres
that induce peak concentrations downwind of surface mines.
Each of the four alternative escape fraction equations was coded into a
FORTRAN algorithm, and tested in the ISC model with input data from a
hypothetical surface coal mine. Run times for the four different algorithms
were recorded during the tests. As expected, the equations using the more
detailed analysis technique required more computer processing time. The two
techniques based on the linear model (Alternatives 1 and 3) required
approximately the same processing time as the original version of ISCST.
Alternative 2 (Constant-K, detailed model) increased the run time by roughly a
factor of 1.5, while Alternative 4 (Variable-K, detailed model) increased the
run time by roughly a factor of 5.
-------
-------
2.0 BACKGROUND
Pit retention is the term used to describe the tendency for particulate
matter released inside a surface mine pit to remain inside the pit. The pit
retention phenomenon is important because most air quality models that are
used to simulate particulate dispersion from surface mines treat these
emissions as if they occurred at grade level, and ignore the possibility that
a portion of the particulate matter may be trapped inside the pit, or that the
characteristics of the dust plume may be altered by the presence of the pit.
Two years ago the U.S. EPA's Office of Air Quality Planning and
Standards initiated a comprehensive study of the pit retention phenomenon
(EPA, 1985). This investigation began with a data collection field study at
four Western surface coal mines. Meteorological parameters were measured
simultaneously in and out of the mine pits for a total duration of
approximately 300 hours. In addition, a smoke release program was conducted
to provide data concerning air motion within the pits. At each of the four
mines, smoke generators at the bottoms of the pits were used to release
discrete 10-second puffs of diesel fuel smoke. An observer positioned at the
top of the pit filmed each smoke release on a video cassette recorder (VCR).
Roughly 800 such smoke release experiments were conducted at the four mines,
and the VCR observations were synchronized with the in-pit and out-of-pit
meteorological measurements.
These field data were later reduced and interpreted in order to
investigate relationships between meteorological variables and the behavior of
the smoke puffs. For each smoke release experiment, the time from initial
smoke release until the smoke puff exited the pit, or until the smoke puff was
no longer visible, was determined by viewing the VCR tape. This time was used
to define a discrete smoke release "episode". All of the data determined by
analyzing the VCR tapes, were organized into episodes. Meteorological data
(wind speeds, wind directions, temperatures, etc.) were averaged over the
episode duration for analysis, along with subjectively determined variables
-5-
-------
(characteristic flow pattern and location of plume exit), and elapsed time
duration of the smoke release episode. This information formed the data base
for subsequent analysis.
2.1 SUMMARY OF PREVIOUS DATA ANALYSIS
Several different kinds of analyses were made with the data base, as
discussed in "Dispersion of Airborne Particulates in Surface Coal Mines"
(EPA, 1985). A comparison of winds in and out of the pits during smoke
releases showed that in-pit wind speeds are, on the average, 25% less than the
out-of-pit wind speeds, and wind speeds both in and out of the pit were
positively correlated. Wind direction in and out of the pits, however, was
not correlated, so that a knowledge of wind direction at the top of the pit
(ie., at grade level) cannot predict wind direction within the pit.
The smoke puff observations by themselves did not provide a quantitative
measure of particulate pit retention. Consequently, a part of the data
analysis was devoted to inferring escape fraction from the smoke puff
observations by using two independent methods one based on a simple
settling model, the other based on the source depletion particle deposition
model. Both methods relied on assumed particle size distributions: one for
particles smaller than 30 microns aerodynamic diameter (called the universal
distribution), and one for particles up to 130 microns aerodynamic diameter
(called the EDS distribution). It was found that the value of escape fraction
inferred from both the settling and the deposition models is greater for
unstable and neutral atmospheric conditions, as shown in Table 2.1 This
suggests that stable atmospheres may suppress vertical motion causing
particulate matter to be retained in the mine pits. In a similar manner, the
1. A quantitative measure of pit retention is expressed by the escape
fraction, e , which is equal to the total mass of particulate that escapes from
the pit, divided by the mass of particulate emitted within the pit.
-6-
-------
TABLE 2.1
ESCAPE FRACTION SHOWN BY STABILITY
PARTICLE SIZE
DISTRIBUTION
UNIVERSAL
EDS
STABILITY^)
UNSTABLE
NEUTRAL
STABLE
UNSTABLE
NEUTRAL
STABLE
SETTLING
MODEL
1.00
1.00
1.00
0.81
0.90
0.70
DEPOSITION
MODEL
0.93
0.81
0.58
0.59
0.36
0.21
WINCES
EQUATION
0.99
0.92
0.58
0.90
0.59
0.20
1. "A" stability class used for unstable; "D" used for neutral; "F" used
for stable.
values of escape fraction determined by the settling and deposition models
were grouped by National Weather Service wind speed class, as shown in Table
2.2. This analysis indicates that the escape fraction increases with
increasing wind speed, as would be expected higher wind speeds tend to
remove more particulate matter from the pits.
TABLE 2.2
ESCAPE FRACTION BY WIND SPEED
DISTRIBUTION
UNIVERSAL
EDS
WIND SPEED
(CLASS)
1
2
3
4
5
1
2
3
4
5
EXIT VELOCITY
(SETTLING)
1.00
1.00
1.00
1.00
1.00
0.75
0.85
0.96
0.96
0.99
SOURCE DEPLETION
(DEPOSITION)
0.78
0.84
0.86
0.88
0.88
0.35
0.46
0.43
0.43
0.43
WINCES
EQUATION
0.90
0.91
0.95
0.95
0.96
0.70
0.70
0.73
0.69
0.76
-7-
-------
Two analytical expressions which predict escape fraction from
meteorological and mine pit parameters were tested. The Winges equation
(Winges, 1981), which expresses escape fraction as a function of pit depth,
vertical diffusivity, and deposition velocity, was found to be superior:
(-6"
where e is the escape fraction
V, is the larger of deposition or settling velocity, m/s
d 2
K is vertical diffusivity, m /sec
z
H is pit depth, m.
The Winges equation was applied independently to each of the smoke
release episodes, and the average values of predicted escape fraction were
grouped by Pasquill-Gifford stability class and by wind speed. These
predicted escape fractions are shown in Tables 2.1 and 2.2, where they are
compared with the escape fractions inferred from the measured field data.
Reasonably good agreement is indicated between the escape fractions inferred
from the settling and deposition models and those predicted by the Winges
equation when the data are grouped by stability class. The values of the
Winges escape fraction decrease as the atmosphere becomes more stable, just as
the measured values do.
When the data are grouped according to wind speed, as in Table 2.2, the
agreement between escape fraction inferred from the field data, and escape
fraction predicted by the Winges equation, is not especially good. One reason
for this may be that the Winges equation does not include wind speed
explicitly in estimating escape fraction. This suggests that the performance
of the Winges equation may be improved by incorporating wind speed into the
equation.
-8-
-------
2.2 OVERVIEW OF CURRENT DATA ANALYSIS
The findings from the previous analyses (EPA, 1985) suggest two kinds of
follow-on investigations. First, the moderate success of the Winges equation
in predicting escape fractions inferred from the field data leads to a
question of whether the Winges equation can be improved. In particular, can
the agreement between predicted and inferred escape fractions be improved by
introducing new variables (eg., wind speed), or by modifying the equation to
take into account more accurate representations of dispersion. These
questions are explored in Chapter 4 of this report.
The second follow-on investigation concerns the meteorological data
collected in and out of the pits. EPA's analysis in January 1985 looked at
meteorological conditions that were coincident with smoke puff releases, and
were averaged over a time period equal to the episode duration of the smoke
puff release. This meant that the values of sigma-theta measured in and out
of the pits were converted to alphabetic stability class over time periods
equal to the smoke puff episodes, which were generally between 30 seconds to
ten minutes in duration. The equivalent alphabetic stability class for these
short sampling times was predominantly "D", and as a consequence, further
analyses of sigma-theta stability class were not performed. In this present
report the values of sigma-theta stability class are recomputed over fifteen
minute and one-hour time intervals, as described in the "Guideline on Air
Quality Models (Revised)" (EPA, 1984). The details and findings of this
investigation are discussed in Chapter 3.
-------
-------
3.0 SIGMA-THETA DATA ANALYSIS
The use of sigma-theta as a measure of atmospheric stability is
especially attractive in the analysis of the field data because this is the
only turbulence parameter that was measured independently and simultaneously
inside and outside the pit. ' In addition, the recent "Guideline on Air
Quality Models (Revised)" (EPA, 1984) recommends the use of sigma-theta as an
acceptable measure of stability class, and provides a uniform method to
convert short-term values of sigma-theta to one-hour stability classes.
Because a majority of the alphabetic stability classes computed previously
were Category "D", there was some question about the accuracy of the field
data itself. A quality assurance audit of the instrumentation and the
software used to measure sigma-theta was made.
3.1 SIGMA-THETA AUDIT
Air Sciences, Inc., the company responsible for collecting the field
data, was asked to perform an audit of the wind direction and wind speed
instrumentation and the software in the data logger that were used to collect
the 1983 field data. Their findings are included in Appendix A of this
report. In summary, there were two separate causes of error discovered in the
collection and calculation of sigma-theta values:
POLLING FREQUENCY. The data logger used to interrogate the wind
direction sensor was programmed to poll once every 10 seconds, and
then compute a one-minute standard deviation from these six samples.
This is a low polling frequency. The effect of the low polling
frequency would be to introduce random errors in computed one-minute
values of sigma-theta. That is, some values of sigma-theta would be
artificially too big, and some values would be too small, but over
1. Sigma w, the standard deviation of vertical wind speed, was only measured
out of the pit (EPA, 1985).
-11-
-------
a large number of computed sigma-thetas the random errors would cancel
one another. The effect of this error would tend to diminish with the
number of computed sigma-theta values, and it would be expected that
over a full one-hour time interval the error inherent in individual
one-minute sigma-theta values would cancel out. Consequently, the
polling frequency error is not important in one-hour sigma-theta values.
COMPUTATION OF STANDARD DEVIATION. In computing sigma-theta, the
software used in the data logger employed an equation for sample
variance , as opposed to population variance. The difference in
variance computed with the two equations is insignificant for a large
number of samples, but when the number of samples, n, is small, the
difference can be significant (Mendenhall, 1968):
"...It can be shown that for small samples (n small) the sample
variance tends to underestimate [sigma squared], and that the
formula
n - I
provides better estimates."
The error introduced in this manner is systematic, and can be corrected very
easily by multiplying the individual one-minute sigma-theta values by 1.1,
which is derived as follows:
1/2
sigma = (6/5) sigma
e corrected ' 6 one-minute
= l.ls igma
0 one-minute
where sigma corrected = corrected value of sigma-theta
sigma . = value of sigma-theta from field data
one-minute 6
(6/5) = ratio of n and n-1
1. Variance is standard deviation squared,
-12-
-------
The conclusion from the instrumentation and software audit is that sigma-theta
values computed from the field data will be accurate if 1) the averaging time
for computation of overall sigma-theta is increased so that random errors will
cancel, and 2) individual one-minute sigma-theta values are multiplied by 1.1
before processing.
3.2 DATA AVERAGING
When the instrument and data logger software audit was completed, the
field data were averaged into discrete, consecutive, one-hour time intervals.
The one-hour averaging time was chosen since this is the standard time
interval used for dispersion model input, and because the "Guideline on Air
Quality Models (Revised)" (EPA, 1984) relates alphabetic stability classes to
one-hour sigma-theta values. The field data base was that submitted at the
conclusion of the field data gathering effort (Hittman and Air Sciences,
1983), except that spurious, illegible characters introduced into the data
base during transcription from cassette to magnetic tape had been removed.
The period of record for the field data base is shown in Table 3.1, by mine.
TABLE 3.1
FIELD DATA PERIOD OF RECORD IN 1983
MINE PERIOD OF RECORD
YAMPA JUNE 28 (1000 - 1400 HRS)
JUNE 29 (0800 - 1400 HRS)
JUNE 30 (0800) - JULY 2 (0700)
CABALLO JULY 11 (1200)- JULY 16 (0200)
SPRING CREEK JULY 19 (0700) - JULY 22 (1000)
ROSEBUD AUGUST 1 (1100) - AUGUST 5 (1000)
One-hour averages of wind speed (in and out of the pit) and wind
direction (in and out of the pit) were computed by scalar averaging and unit
vector averaging, respectively.
-13-
-------
The one-hour averages of sigma-theta in and out of the pit were computed
as follows. First, individual sigma-theta values were multiplied by 1.1 to
correct for the error in the data logger software. Next, fifteen consecutive
one-minute values of sigma-theta were summed and averaged by the
root-mean-squared (rms) method for each quarter hour within a fixed one-hour
time period. Finally, the four consecutive 15 minute averages of sigma-theta
were combined into an hourly sigma-theta value with the equation
2 2 2 2 1/2
+ sigma.,. + sigma ,. + sigma.
sigma. _ + sigma.. _ + sigma. _ + sigma. _ ^
( ; )
where sigma.. is one-hour average sigma theta
sigma..- is fifteen minute sigma-theta
This procedure has been recommended to compute average sigma-theta values in
order to minimize wind direction meander effects (EPA, 1984).
Hourly daytime stability classes were computed using sigma-theta
classifications and wind speed criteria shown in Table 9-2 and hourly
nightime stabilities were computed using Table 9-3 criteria (EPA, 1984). Both
of these methodologies eliminate unrealistic occurrences of stable and
unstable conditions that would not occur with the Pasquill-Gifford stability
typing scheme. For each hour of data, two independent sigma-theta stabilities
(one in the pit, one outside the pit) were calculated.
The Pasquill-Gifford (P-G) stability class was determined from cloud
cover and ceiling recorded in the field observer's logs, combined with average
out-of-pit wind speed during each hour. The procedures used to compute P-G
stability class were those used by the National Climatic Data Center (NCDC) in
deriving STAR distributions.
The final parameter computed was the ratio of the hourly average
sigma-theta measured put of the pit, divided by the sigma-theta measured in
-14-
-------
the pit. This ratio of sigma-theta values is a dimensionless variable that
indicates whether turbulence (as measured by sigma-theta) is greater in or out
of the pit. Values of sigma-theta ratio smaller than 1.0 indicate greater
turbulence in the pit than out of the pit.
The processed hourly data averages are included in Appendix B of this
report. Each data record (one horizontal line) shows a one-hour average of
the meteorological parameters. Data fields filled with 999s indicate missing
or invalid data.
The one-hour averaged data base shown in Appendix B is different from
the data base used previously (EPA, 1985) to examine escape fractions in this
respect: the data base in Appendix B presents one-hour averages of
meteorological variables, as opposed to averages computed over the smoke puff
release episode time.
3,3 METEOROLOGICAL DATA ANALYSIS
It is seen from the one-hour meteorological data in Appendix B that the
value of sigma-theta outside the pit is almost always smaller than sigma-theta
inside the pit. This is indicated by the value of the parameter RAT (ratio of
the hourly sigma-theta out of the pit, divided by the hourly sigma-theta in
the pit) which is almost always less than 1.0. In fact, of the 247 valid
hourly observations during which both sigma-theta (out) and sigma-theta (in)
were present, only 21 of the observations^ ' indicate that the
sigma-theta ratio is greater than 1.0. This indicates that the horizontal
1. The first measured value of sigma-theta ratios on Julian day 181, at time
8:59, was discarded from the data set. This is the first reading of a new
measurement run, and it appears to be erroneous. The subsequent value of
sigma-theta in the pit at time 9:59 was flagged as incorrect by Air Sciences,
Inc.
-------
wind direction fluctuation inside the pit is greater than outside the pit,
which is as expected. The in-pit wind sensor responds to the mechanical
turbulence caused by airflow over the edge of the pit, and by the wake created
downwind of the pit walls. The out-of-pit sensor is not subject to these wake
effects since it is located above the mechanically induced pit turbulence
region. It should be remembered that sigma-theta does not necessarily measure
vertical mixing, and it would be a mistake to conclude that pollutants inside
the pit would be more thoroughly dispersed vertically than those outside the
pit.
To examine the relationship of sigma-thetas in and out of the mine pits,
all data were segregated into groups defined by values of sigma-theta ratio
less than and greater than 1.0. The average wind speeds in and out of the
mine pits, segregated by sigma-theta ratios, are shown in Table 3.2.
TABLE 3.2
AVERAGE WIND SPEED (kts.) OUT AND IN PIT
GROUPED BY SIGMA-THETA RATIOS
WIND SPEED
LOCATION
OUT
IN
SIGMA-THETA
RATIO < 1.0
7.99
5.71
SIGMA-THETA
RATIO >1.0
4.62
4.12
Table 3.2 shows that when the ratio of sigma-thetas is less than 1.0 (ie.,
when sigma-theta out of the pit is less than in the pit) that wind speeds are
appreciably larger than when the sigma-theta ratio is greater than 1.0. The
relationship of sigma-theta ratios with wind speed can be seen by listing
average values of sigma-theta ratio with wind speed categories, in Table 3.3.
It is evident that the ratio of in-pit and out-of-pit sigma-thetas depends
strongly on wind speed. Values of the ratio decrease with increasing wind
speed.
-16-
-------
TABLE 3.3
COMPARISON OF WIND SPEED AND
SIGMA-THETA RATIO
OUT OF PIT SIGMA-THETA RATIO
WIND SPEED (kts.)
0-3 0.80
4-6 0.66
7-10 0.60
11 - 16 0.52
17 - 21 0.35
GREATER THAN 21 0.40
The relationship in Table 3.3 could be caused by either 1) increased
turbulence in the pit with greater wind speeds, or 2) decreased turbulence out
of the pit with greater wind speeds. Examining the sigma-theta values in and
out of the pit as a function of wind speeds, suggests that the second
explanation is correct.
TABLE 3.4
SIGMA-THETA OUT AND IN PITS
AS A FUNCTION OF WIND SPEED
OUT OF PIT SIGMA-THETA OUT SIGMA-THETA IN
WIND SPEED (kts.) (deg) (deg)
0 -
4 -
7 -
11 -
17 -
GREATER
3
6
10
16
21
THAN 21
14.0
13.9
13.4
10.2
8.3
7.2
19.5
21.7
23.7
20.2
24.8
17.8
The values of sigma-theta out of the pit decrease with increasing wind speed,
as seen in Table 3.4. This means that the horizontal fluctuations of wind
direction decrease with greater wind speed out of the pit, as would be
expected, since higher wind speeds tend to increase wind direction
persistence. Inside the pit, however, mechanical turbulence of the pit itself
-17-
-------
dominates the flow, and horizontal wind directions are more nearly constant
regardless of wind speed.
Table 3.5 shows the number of occurrences of alphabetic stability class
determined by the Pasquill-Gifford method versus those indicated by
sigma-theta in the pit. If the two methods agreed perfectly, then all of the
values in Table 3.5 would lie along a line drawn from the upper left of the
Table to the lower right. When grouped by general stability category (ie,
unstable, neutral, and stable) the agreement is fairly good. For any given
value of Pasquill-Gifford stability class, the sigma-theta stability class
in the pit tends to be slightly more unstable. Similarly, the comparison
of Pasquill-Gifford stability with out-of-pit sigma-theta stability (shown in
Table 3.6) exhibits similar general agreement, with the sigma-theta method
showing more neutral stability ("D" class) than the Pasquill-Gifford method.
Finally, a comparison of sigma-theta stabilities in and out of the pits shown
in Table 3.7 indicates that stabilities inside the pit are, by and large, more
diverse than stabilities outside the pit. While there are 136 occurrences of
neutral ("D") stability measured out of the pit, there are only 55 measured in
the pit.
TABLE 3.5
NUMBER OF OCCURRENCES OF STABILITY CLASS
DETERMINED BY PASQUILL-GIFFORD AND
SIGMA-THETA MEASURED IN-PIT
PASQUILL-GIFFORD
A
B
C
D
E
F
A
2
31
30
3
0
0
SIGMA-THETA
B
0
9
21
10
0
0
C
0
5
8
2
0
0
IN-PIT
D
0
4
7
24
10
11
E
0
0
0
1
6
19
F
0
0
0
0
10
38
-13-
-------
TABLE 3.6
NUMBER OF OCCURRENCES OF STABILITY CLASS
DETERMINED BY PASQUILL-GIFFORD AND
SIGMA-THETA MEASURED OUT-OF-PIT
PASQUILL-GIFFORD
A
B
C
D
E
F
A
0
13
1
0
0
0
SIGMA-THETA OUT-OF-PIT
B
4
19
7
0
0
0
C
1
22
26
0
0
0
D E
1 0
11 0
36 0
42 0
24 7
41 32
F
0
0
0
0
0
18
TABLE 3.7
NUMBER OF OCCURRENCES OF STABILITY CLASS
DETERMINED BY SIGMA-THETA OUT-OF-PIT
AND SIGMA-THETA MEASURED IN-PIT
SIGMA THETA
A
B
C
D
E
F
IN-PIT
A
12
0
0
1
0
0
B
14
2
1
0
0
0
SIGMA-THETA
C
20
13
6
1
0
0
OUT-OF-PIT
D E
18 0
24 0
8 0
45 7
17 5
24 17
F
0
0
0
1
4
7
-19-
-------
In general, however, the agreement between all three stability typing
schemes (Pasquill-Gifford, sigma-theta in the pit, and sigma-theta out of the
pit) is reasonably good. Table 3.8 shows the number of hours in which the
various stability classes differ by 0, 1, 2, 3, 4, or 5 categories. Table 3.8
shows that the Pasquill-Gifford method and the sigma-theta in-pit yield the
same alphabetic category 87 hours out of a possible 251 hours, and they differ
by one category 106 hours. From this, it can be seen that stability class
determined by the Pasquill-Gifford method and by measuring sigma-theta in the
pit are within one stability category for (87 + 106/251 = .77) 77% of the
valid data hours. Similarly, the P-G and the sigma-theta (out-of-pit)
stabilities agree within one stability class 82% of the hours, and sigma-theta
stabilities in and out of the pit agree within one stability class 64% of the
time.
TABLE 3.8
DIFFERENCES IN STABILITY CLASSES
CLASSES
DIFFER BY
0
1
2
3
4
5
P-G
& SIGMA-THETA
IN-PIT
87
106
55
3
0
0
P-G
& SIGMA-THETA
OUT OF PIT
112
138
54
1
0
0
SIGMA-THETA
IN & OUT
OF PIT
77
82
69
19
0
0
TOTALS
251
305
247
-20-
-------
This good agreement between sigma-theta and Pasquill-Gifford stabilites
seemingly contradicts the poor agreement between sigma-theta and P-G stability
detected in the previous examination of smoke release episode stabilities
(EPA, 1985) in which the majority of the sigma-theta stability classes were
found to be "D" class. The most likely reason for poor agreement between
sigma-theta stabilities and P-G stabilities in the previous investigation is
that the data sampling times in the smoke release data were limited to the
episode duration, which varied from one minute to, at most, 20 minutes. Over
these short time periods the horizontal wind direction fluctuations are small.
-21-
-------
4.0 ANALYSIS OF ESCAPE FRACTION EQUATIONS
In the previous report (EPA, 1985), TRC evaluated two equations for
computing the escape fraction. The evaluation data were the inferred escape
fractions from the video-tape interpretations. The comparison of the data
with the available formula indicated that the equation developed by Winges
offered promise. This section details the efforts to extend the original
Winges formula.
The original Winges equation was based on a theoretical analysis of
diffusion of particles from an open depression in the ground. The derivation
of the equation will be presented later in this document, but a general
discussion of the overall technique and assumptions is pertinent here. The
diffusion of particles from a sub-surface depression can be treated as a
steady-state process, such that three phenomena are in constant mass-balance:
the emission of dust in the pit, the deposition of dust on the surfaces of the
pit and the flux of dust out of the pit.
The mass-balance approach was augmented by the key assumption that the
transport of material within the pit could be completely characterized by the
diffusion process that is, that the mean properties of the wind do not
result in any transport of the material out of the pit, rather only the random
motions of the wind are responsible for dust loss to the atmosphere. This
assumption may be paraphrased as saying that there is no vertical wind within
the pit, but there is vertical turbulent diffusion. It is important to
emphasize that the concern here is only with vertical motion of the air since
the emissions are assumed to occur within a cut from a flat surface and
vertical motion is necessary for the escape of particles.
The key parameters are those that concern the vertical diffusion of
particles from the pit. In the original Winges equation a simple gradient
transfer approach was taken in which the vertical flux was assumed to be
proportional to the gradient of concentration of dust within each vertical
layer in the pit and the proportionality constant, called the eddy
-23-
-------
diffusivity, was assumed to be a constant for all heights. This approach,
which will be called the Constant-K approach, yielded a simple equation for
computation of the escape fraction. For implementation of the above equation,
a value of the eddy diffusivity was taken from the literature. Evidence from
Draxler (1979) suggests that different eddy diffusivities should be used for
different stabilities.
In addition to the experimental evidence offered by TRC, there is ample
evidence from the scientific literature that the eddy diffusivity, and hence
the escape fraction, is related to the wind speed (Draxler, 1979). The
purpose of the current investigation is to determine if wind speed could be
incorporated into the previous equation. The logical place to incorporate the
wind speed into the earlier formula is through the characterization of
diffusion (the eddy diffusivity). Horizontal wind speed influences the
vertical diffusion near a surface because a considerable portion of the
turbulence near the surface results from frictional shearing caused by the
wind as it passes over the surface. If the wind exerts more force on the
surface (as a result of greater wind speeds), then it can be expected to
create more turbulence.
The current report addresses two general avenues for incorporation of
the effect of wind speed on pit retention and within each of these avenues
there are options. All of these techniques are presented rather than
presenting a single method in the interest of completeness. Later in this
chapter, all of the newly-developed techniques will be compared with the
experimental data. In addition to a comparison with the experimental data as
interpreted in the earlier study, this report also offers a new interpretation
of the video tape data. Without detailed experimental data it is not possible
to determine if the more complex techniques result in greater accuracy.
4.1 DERIVATION OF THE CANDIDATE ESCAPE FRACTION EQUATIONS
As stated earlier, the original Winges equation was based on the
assumption that the eddy diffusivity is a constant throughout the pit. The
assumption of a constant eddy diffusivity is based in large part on the lack
-24-
-------
of understanding of what the dispersion behavior inside a pit really is. It
is likely that the flow and turbulence patterns inside a pit are highly
complex and not easily represented. Most research on turbulence
characteristics are for flow over uniform flat surfaces or simple geometric
shapes. Even the few studies which have been performed on shapes similar to a
mine pit would not be expected to generalize to all orientations or
configurations. Thus, the simple assumption of a horizontally well mixed
volume of air with a single value of the eddy diffusivity was used in the
original Winges equation because the actual behavior of the eddy diffusivity
in the pit is unknown.
4.1.1 THE ORIGINAL WINGES EQUATION
The derivation of the original Winges equation is not in the open
literature and may not be available to some readers. The fraction of material
which escapes the pit may be represented by the following equation:
e- (1)
where: e = escape fraction (dimensionless)
F = flux of material from the pit (g/sec-m2)
E = emission rate within the pit (g/sec-m2)
By a simple mass balance argument, the dust emitted in the pit must
either be deposited on the internal surfaces of the pit or transported as a
flux out of the pit. Mathematically this is represented by:
E = F + D (2)
2
where: D = deposition in the pit (g/sec-m )
The original Winges equation attempted to treat a very simplified
dispersion scenario, and a number of assumptions were made to simplify the
mathematical solution. These include:
-25-
-------
1. All emissions occur at the bottom of the pit.
2. The only mechanism for transport of material out of the pit is
turbulent diffusion. This assumption, discussed earlier, means that
vertical wind speeds will be ignored.
3. The vertical flux of material is constant with height. This must
occur if the flow is in steady-state, otherwise concentrations would
be building-up inside the pit.
4. The turbulence within the pit is constant throughout the pit. This
is the constant eddy diffusivity assumption.
5. Deposition occurs at the bottom of the pit and is proportional to
the concentration at the bottom of the pit. The assumption of
deposition being proportional to concentration at the ground is well
supported in the literature (see for example, Chamberlain and
Chadwick, 1953). The proportionality constant has the units of a
velocity and is termed the "deposition velocity".
6. Concentrations directly above the pit, resulting from pit emissions,
fall to zero at some height above the pit. This condition is
necessary as a boundary condition for the differential equations to
be solved. It is a reasonable assumption, since emissions that are
mixed to the top of the pit would be carried away by the prevailing
wind, so that the wind would provide a constant supply of "clean"
air at the top of the pit. The original Winges equation used the
assumption that concentrations fall to zero at the top of the pit,
because it turns out that this results in the greatest percentage of
material being lost and thus may be viewed as a conservative
assumption. This assumption is generalized here to simply say that
concentrations must fall to zero at some height above the bottom of
the pit, H, and that height may be specified by the user. To be
conservative, the user may select a value of H equal to the depth of
the pit so that the zero height is the top of the pit and thereby
maximize the escape of emissions.
The gradient transfer approach for dealing with turbulent diffusion is
to model the turbulent behavior using equations that match laminar flow. In
laminar flow the flux of material across any surface resulting from diffusion
is proportional to the concentration gradient between the two bodies of fluid
on either side of the surface. The proportionality constant is called
diffusivity. The turbulent motions called eddys result in far more transfer
of material than the laminar diffusion process. However, it is still the
gradient in concentration between two bodies of fluid that results in transfer
of material, since the eddy motions result in exchange of fluid across the
-26-
-------
boundary. Thus, the concept evolved of assuming the diffusion to be
proportional to the concentration gradient, but here a much larger
proportionality constant, called the "eddy diffusivity" was used
(Bird et al., 1960).
There is a large difference between a laminar diffusivity and an eddy
diffusivity. The laminar diffusivity is a function of the physical properties
of the fluid, such as its viscosity and temperature. The eddy diffusivity is
a property of the flow, and for a given fluid and temperature can vary widely
depending on the energy of motion of the fluid and the shearing forces and
other phenomena. For these purposes here, it is assumed the vertical motion
of particles emitted in the pit can be represented by a gradient transfer
equation:
F = -K*
where: K = eddy diffusivity (n^/sec)
X = concentration (g/m3)
z = vertical dimension (m)
The general approach in the Constant-K Model is to assume that the eddy
diffusivity, K, is a constant with respect to height of the pit. This
constant assumption allows easy integration of equation (3) as follows:
X = ~ K Z + C
where: C = constant of integration
It is necessary to evaluate the constant of integration with a boundary
condition, and for this purpose, we use the assumption that concentration
falls to zero at some height above the surface, H. This is accomplished as
follows:
0 = - | H + C (5)
C - H
where: H = height at which X = 0 (m)
-27-
(6)
-------
Now, equation (4) becomes:
X = |(H - z) (7)
The above equation can be used to evaluate the term "D" in equation (2) with
one additional assumption. The deposition at the surface must be proportional
to the concentration at the surface. Mathematically this can be represented
as (Chamberlain, 1953):
where: x = concentration at the surface (g/m3)
Z0
u, = deposition velocity (m/sec)
Equation (7) allows one to compute the concentration at the surface as follows
X - (H - zo) (9)
where: ZQ = some small height, usually called the roughness
height (further detail provided later) (m)
Reforming equation (1) and substituting from above as follows:
e -- - - (12)
1 + % - ZQ)
-23-
-------
Since the roughness height is usually very small when compared to H, it is
possible to ignore the roughness height and express the equation as follows:
(13)
The above equation is the one used previously (EPA, 1985) in the evaluation of
the alternate pit retention formulae and referred to as the original Winges
equation.
4.1.2 ALTERNATIVE 1 CONSTANT-K USING LINEAR MODEL
The simplest method of incorporating the wind speed into the above
formula is to keep the assumption of a constant eddy diffusivity and calculate
the value of the eddy diffusivity to be used as a function of the wind speed.
This report investigates two general methods for computation of the eddy
diffusivity as a function of wind speed. The first of these is based on an
assumed linear relationship between wind speed and eddy diffusivity. The
second, which will be presented later, involves a more detailed approach for
characterizing the eddy diffusivity. The linear assumption results from a
number of other assumptions about the relationship between turbulence and wind
speed and the derivation of this relationship is presented in the following
paragraphs.
The wind speed will generally be measured outside of the pit at some
reference height. It is well known that the wind speed in the lower layers of
the atmospheric boundary layer increases with height above the surface
(Turner, 1970). The shape of the wind speed profile, as it is called, is
reflective of the momentum balance of the flow at the surface. In one
simplified analysis, the wind speed profile is characterized by two
parameters, and these are usually expressed in a logarithmic equation known as
the logarithmic profile (Monin and Yaglom, 1971).
-29-
-------
u = ln<) (14)
where: u = wind speed (m/sec)
UA= friction velocity (m/sec)
k = von Karman constant, usually
assumed to be 0.35
The two parameters, friction velocity and roughness height, will be used
extensively in the analysis throughout this document, and need further
explanation. The wind moving over the surface of the earth creates a shear
stress at the surface. This shear stress, when divided ty the density of the
air to reduce it to its kinematic properties, has the units of a velocity
squared, and when the square root is taken the result is called the friction
velocity. The friction velocity, then, may be thought of as a measure of the
shear stress exerted by the wind on the surface of the earth. The surface
roughness height is a measure of the surface protrusions which create drag on
the wind as it passes. The greater the surface area offered by these
protrusions, the greater the drag, and the more gradual the increase of the
wind speed with height.
The shear stress at the surface is a way of expressing the transfer of
momentum by turbulent motion to the surface of the earth (Monin and
Yaglom, 1971). The process of the transfer of momentum in turbulent flows is
very similar to the process of the transfer of particles, thus it is useful to
examine the momentum transfer process as reflected in the wind speed profile
to see what it says about the particulate diffusion process. As with the
diffusion of particles, a gradient transfer representation can also be used
for the transfer of momentum. In the momentum transfer case, the gradient is
of the wind speed rather than the concentration of particles The
proportionality constant here is customarily called the kinematic eddy
viscosity instead of the eddy diffusivity as used for the diffusion of
particles earlier.
-30-
-------
However, the mechanism for momentum transport is exactly the same as the
mechanism for turbulent transfer of gasses and particles, and consequantly,
researchers have used the eddy viscosity as a measure of the eddy
diffusivity. The mathematical characterization of this process is as follows
(Monin and Yaglom, 1971):
U2=-Kfi (15)
where: K = kinematic eddy viscosity (n^/sec)
Differentiate the wind speed profile (equation 14), to develop the following:
Sz k z
Substituting and reforming, obtain the following:
Kvu*
(16)
(17)
(18)
As stated earlier, the wind speed will be measured at some reference height,
and consequently, one can compute the resulting eddy viscosity at the same
reference height by reforming the logarithmic profile to solve for the
friction velocity and substituting the resulting equation into the above
equation for the eddy viscosity. This is shown in the next few equations:
Ğ* - ? d9)
where: z ,. = reference height of wind speed
measurement (m)
,,
K -- US - z (20)
v z f ref v ^'
20
-31-
-------
Inserting the solution for the eddy vicosity in place of the eddy diffusivity
in equation (13) for the escape fraction yields a solution:
u, In C-) H
(21)
i
1
uk^z ,
ref
4.1.3 ALTERNATIVE 2 CONSTANT-K USING A MORE DETAILED MODEL
A major shortcoming of the previous approach is that, while it
incorporates wind speed in the equation, it has lost the capability to include
stability. The original Winges equation allowed the user to select eddy
diffusivities based on stability, if desired. The equation in Alternative 1
has used a simplified measure of the turbulence in the atmosphere to
substitute for the eddy diffusivity. A problem arises because the logarithmic
profile, while a reasonable approximation to the wind speed profile in uniform
flow over a flat plate, ignores the effect of the temperature structure of the
atmosphere in enhancing or inhibiting vertical mixing. Temperature structure
can have a significant effect on the vertical mixing of both mass and
momentum, and the atmospheric stability is an often-used concept to
characterize this influence.
There is an alternative approach to the one presented in section 4.1.2.
It involves considerably more detail and will be presented but not derived
here. It is fundamentally different than the previous approach in that it is
an empirical approach rather than a theoretical approach. It uses a parameter
called the Monin-Obukhov length to characterize the stability aspects of the
flow. The Monin-Obukhov length characterization of temperature structure
influences on dispersion is viewed as an improvement over the previous
stability classification scheme by the meteorological community.
The eddy diffusivity is computed using the following formula
(Draxler, 1979):
K = ^ (22)
where: $ = normalized temperature profile
h
-32-
-------
The normalized temperature profile may be computed by one of two formulas and
uses the Monin-Obukhov length L. In fact, it is necessary to compute the
Monin-Obukhov length first because the choice of formulas to use for the
normalized temperature profile is made with the quantity Z/L (also used as a
measure of the stability). If the stability is unstable (Z/L < 0) then the
following formula is used for the normalized temperature profile:
If the stability is stable or neutral (Z/L S 0) then the following formula
is used for the normalized temperature profile:
, = 0.74 + 5J- (24)
n L
The computation of the Monin-Obukhov length is complicated. First, one must
compute the Bulk Richardson Number, B, using the following equation:
where: g = gravitational acceleration
(9.81 m/sec2)
T = ambient temperature ( K)
A9 = potential temperature gradient
Then the Richardson Number itself, Ri, is calculated from the Bulk Richardson
using the following equation:
B = _^i (26)
~ m o
>
where: ^ and are defined below
m m
-33-
-------
For stable and neutral conditions:
*m (1 - 5R1)
(27)
While during unstable conditions:
*-= (29)
z0
+ 2(tan~1? - tan'^o)} (30)
5 = (1 - 15R1)1* (31)
Co Ğ (1 - 15Ri~^-r (32)
Zref
It will be noted that equation (26) cannot be solved directly for the
Richardson Number. In fact, solution of the equation is a tedious numerical
process. A computer algorithm for the solution of this complicated set of
equations is included in Appendix C. Alternative numerical solution
techniques may be an improvement and should be investigated. Once the
Richardson Number has been computed, the Monin-Obukhov length is computed by
one of two formulas. If the Richardson Number is less than zero, the
following formula applies:
- Ri (33)
-34-
-------
If the Richardson Number is greater than or equal to zero, the following
formula applies:
2 - Ri (34)
L (1 - 5R1)
The friction velocity is also computed using an..- empirical equation of the
form:
ku
u* =
z
Once the eddy diffusivity is computed using the above analysis, it is
again assumed to be a constant within the pit and escape fraction is computed
using equation (13), which has been repeated here for convenience.
E = 7^~
K
4.1.4 ALTERNATIVE 3 VARIABLE-K USING LINEAR MODEL
The previous two sections presented alternate methods of extending the
earlier equation to include wind speed, while at the same time maintaining the
assumption that the eddy diffusivity is a constant throughout the pit. It
will be noted that both of the alternatives presented thus far require the
input of height in the computation of the eddy diffusivity. Section 4.1.2
used the reference height of the wind speed monitor as the height to input
when computing the eddy diffusivity. Since the equations imply that eddy
diffusivity is a function of height, it seems logical to investigate the
implications on the escape fraction if the eddy diffusivity is input into the
current analysis as a function of height.
-35-
-------
As with the Constant-K methods, the Variable-K methods will investigate
two separate options for characterization of the eddy diffusivity: the linear
model and the more complex model, using the Monin-Obukhov length.
Equation (20) shows how the eddy viscosity can be calculated using the
linear model. The eddy viscosity is then assumed to be equivalent to the eddy
diffusivity based on the similarity of mass and momentum transfer processess
(Bird, et al, 1960). The height used to compute the eddy diffusivity appears
in two places in equation (20). To generalize the equation for application at
all heights, the z - is replaced in one of the occurrences by z. Equation
(20) then becomes:
K, = z (36)
It should be noted that z . was not replaced by z in the logarithmic
term because the friction velocity is still a constant established by a single
measurement of u at a reference height (using equation 19 which was then
substituted into equation 18 to produce equation 36). The variable
relationship for the eddy diffusivity in equation (36) is then substituted
into equation (3) to develop a new relationship for the vertical flux:
It is important to note that the vertical flux is still a constant, as
required by the steady-state assumption. Therefore, the integration of
equation (37) is still possible to develop a relationship for calculating
concentration as a function of height. The details of the integration are as
follows:
-36-
-------
1^ (38)
Fin (
8z (39)
Since X|Z=H = 0
H
The deposition at the surface is computed using equations (8) and (41) as
follows:
Finally, the escape fraction is computed from equations (11) and (42) as
follows:
(43)
The above equation can be used much as equation (21) is used.
-37-
-------
4.1.5 ALTERNATIVE 4 VARIABLE-K USING THE MORE DETAILED MODEL
Similar to the Constant-K models, equation (43) fails to allow the
escape fraction to be computed as a function of stability. It is possible to
overcome this limitation by using the more complex technique for computing the
eddy diffusivity as a function of the Monin-Obukhov length. The more complex
technique also reveals eddy diffusivity to be a function of height. As with
the Linear Model, for the more detailed approach in Section 4.1.3, this report
recommended using the reference height of the wind speed monitor to compute
the constant eddy diffusivity to be used in the escape fraction computation.
It is possible to generalize this process by allowing the eddy diffusivity to
vary with height in the computation of the escape fraction. The following
equation illustrates the generalization of equation (3):
F = - K(z) (44)
Equation (44) can be integrated over the range of z from the roughness height
to H as follows:
9X - - 3z (45)
3z
Since the concentration is zero at z=H, the following relationship is
developed for the concentration at the roughness height:
-38-
-------
Substituting equation (48) in equation (8) and the result into equation (11)
yields the following expression for the escape fraction:
e £ (49)
1 + u, r (~-r) 3z
d z0 vK(z)'
The integral of the inverse of the eddy diffusivity can be evaluated
numerically, by dividing the vertical extent of the pit (from the roughness
height to H) into a series of finite elements and computing the eddy
diffusivity at each height using the procedures outlines in equations (22)
through (35). The process is not as complicated as it appears. The
Richardson Number, Ri, the Monin-Obukhov length, L, and the friction velocity,
u^ need only be computed once with the height of the wind speed monitor,
z f, being used for z in all places in equations (25) through (35). Only
when computing the eddy diffusivity itself and the normalized temperature
profile in equations (22) through (24) should the actual height in the pit be
used. Once the eddy diffusivities are computed at each height, the inverse of
each is taken, and multiplied by the depth of each finite vertical element.
Finally, the resulting values are summed to calculate the integral in equation
(49).
4.2 EVALUATION OF CANDIDATE EQUATIONS WITH EXPERIMENTAL DATA
4.2.1 EVALUATION DATA
The only data available for the evaluation of the theoretical escape
fraction equations presented in the previous section are the video tape
recordings of the smoke releases documented in the earlier report. The major
problem with using the smoke release data to evaluate the escape fraction for
-39-
-------
mining dust is that the particle size distribution for the smoke particles and
the mining dust are very different. In fact, the density and size of the
smoke particles are sufficiently small to behave virtually like a gas, for
which no pit retention would be expected.
The video tapes do, however, give some information on the residence
times of the smoke in the pit. From these residence times it is possible to
infer some information about the pit retention behavior of actual mining
dust. In the earlier work (EPA, 1985) the escape fraction was inferred from
the measurement data using two separate techniques. The first involved the
computation of an escape velocity by dividing the vertical depth of the smoke
release by the residence time. For any fugitive dust source particles of all
different sizes would be released. Each particle would have a different
gravitational settling velocity and a deposition velocity. The particles were
grouped into classes dependant on size and a characteristic gravitational
settling velocity and deposition velocity were assumed for each class. The
escape velocity was then compared to the larger of the two characteristic
velocities (gravitational settling or deposition) for each class. If the
escape velocity was larger, all particles in the size class were assumed to
escape. If the escape velocity was smaller, all particles in the size class
were assumed to be retained.
Two separate particle size distributions were evaluated with the above
technique, and an overall escape fraction was computed for each distribution.
One particle size distribution came from the PEDCo and TRC 1982 study of coal
mines (PEDCo and TRC, 1982). The second size distribution came from a similar
study conducted by TRC for the mining industry (Shearer, et al, 1981). The
PEDCo/TRC study of size distributions considered only particulate matter
smaller than 30 microns, while the TRC mining industry study looked at
particles as large as 130 microns.
The second technique for computation of the escape fraction developed in
the earlier study involved using an additional theoretical expression. It
computed the escape fraction from the source depletion equation developed by
Van der Hoven (1968) and the meteorological data collected for each smoke
release.
-40-
-------
The escape fraction was calculated using each of the above techniques
and each of the two particle size distributions for all of the roughly 800
smoke releases. The results of this analysis have been reported in the
earlier study.
One observation in the course of the earlier study was that escape
fraction computed using the first of the two techniques above (using the
escape velocity) revealed little or no pit retention for virtually all cases
analyzed. This conclusion disagrees with that of the source depletion
analysis. Since the video tape data do not measure escape fraction directly,
but rather require the user to infer the escape fraction from the measure of
the escape velocity, one effort which was undertaken in the current study was
to re-evaluate the data extracted from the video tapes to determine if there
were other interpretations which could be used to infer the escape fraction.
The end result is that an alternate interpretation of the data was
developed. The escape fractions computed by this new technique provided a
different yardstick by which the various escape fraction equations could be
evaluated. The derivation of the new technique will be presented in the
following paragraphs.
It is necessary to compute how a fugitive dust particle would behave if
exposed to the conditions that the smoke puff was exposed to. It is assumed
that the residence time in the pit would be unaffected by the change from a
smoke puff to a fugitive dust puff, but during the residence time, many of the
fugitive dust particles would be deposited on the surface and walls of the
pit. Those particles which do not deposit during the residence time are
assumed to escape. For a puff of fugitive dust, the rate of deposition is
constantly changing as is the concentration within the puff. However, if a
mathematical characterization of the rate of deposition over time can be
established, the total deposition during the puff's residence in the pit can
be computed by integrating the deposition rate over time from the release time
to the exit time. This analysis is performed with four equations. Generally,
the techniques used to calculate the escape fraction here are the same as
equation (1). However, in equation (1) the concern was for a continuous
release.
-------
Here the concern is for an instantaneous release. Consequently, the flux
term in equation (1) has been replaced by a term representing the amount of
material which escapes after a certain residence time in the pit, R.
E - DEPO(t)I
e (5Q)
where: E = amount emitted in pit
DEPO(t) = amount deposited on the surface from
the release time through t
R = residence time (evaluate DEPO(t) at t=R)
The function DEPO(t) is the total deposition in the pit from the time of
release of the smoke puff until some time, t, later (but not later than the
residence time, R). The term DEPO is evaluated at t=R in the equation above
to determine the total deposition that occurs from the time of release until
the puff exits the pit. The function DEPO(t) is defined as follows:
DEPO(t) - /Q D(t)WL dt (51)
where: D(t) = the average deposition rate over
the area of the pit in g/m2-sec at
any time t
W = the width of the pit
L = the length of the pit
Note that the deposition rate, D(t), is a different function than
DEPO(t). D(t) is the instantaneous value of the deposition rate at any point
in time. As discussed earlier for equation (8), the deposition is assumed to
be proportional to the concentration at the surface for uniformly sized
particles and in the absence of any change in the meteorological conditions or
surface conditions. The concentration is also continuously changing variable
in the puff analysis, so a mathematical representation of this
proportionality, similar to equation (8), is as follows:
-42-
-------
D(t) = X(t)ud (52)
where: x(t) = average concentration in the pit
at any time t
u = deposition velocity (m/sec)
Finally, define the concentration as a function of time by simply
dividing the remaining suspended emissions (amount emitted minus amount
deposited from release time until some later time, t) by the dimensions of the
pit. This assumes the emissions are well mixed throughout the pit. It is
represented as:
x(t)
^ }
HWL
where: H = depth of the pit
The above system of four equations can be solved by first substituting
equation (53) into equation (52) as follows:
u E - u DEPO(t)
D(t) = -^- - (54)
u E u DEPO(t)
D(t) - HWL - - - (55)
Now, equation (51) is substituted into equation (55) and the result is:
D(t) = - D(t)WL dt
D(t) -'-^) dt
D(t) + "if /j; D(t) dt = ^=- (58)
-43-
-------
Equation. (58) is an integral expression which is solved by the following
expression for D(t):
_
D(t) = e" H (59)
Now use equations (51) and (59) to evaluate the function DEPO(t):
ud
u E - -q-t
DEPO(t) = rQ -g- e dt (60)
d -
DEPO(t) = -|- rQ e U dt (61)
Ud
uHE - -g-t
DEPO(t) --£-(- -S_ e H + ~ ) (62)
d d
DEPO(t) ğ -E e + E (63)
Finally, evaluate the escape fraction by substituting equation (63) into
equation (50) as follows:
->
e . EJJL... * (64)
-^(R)
e = e H (65)
Using equation (65), which will be called the residence time analysis
technique, the escape fraction was computed for each of the roughly 800 smoke
releases and for each of the two particle size distributions. This
established an additional measurement interpretation for evaluation of the
theoretical escape fractions.
4.2.2 COMPARISON OF THE CANDIDATE EQUATIONS WITH THE EVALUATION DATA
There are three different evaluation data sets: one based on the escape
velocity, one based on the source depletion equation, and one based on the
-44-
-------
residence time analysis. For each of these evaluation data sets, the escape
fraction has been computed for two different particle size distributions: the
Universal Size Distribution from the PEDCo/TRC Study of 1982 and the Emission
Factor Development Study (EDS) size distribution.
For each of the smoke releases a Pasquill-Gifford stability class
determined in the original study (EPA, 1985) was used here. Most of the
parameters needed by the candidate escape fraction equations were available
from the original data. It was necessary to specify some of the additional
parameters needed by the candidate equations presented earlier. Table 4.1
illustrates the values of the various parameters assumed in this analysis.
TABLE 4.1
PARAMETERS USED IN THE ESCAPE
FRACTION COMPUTATIONS
PARAMETER VALUE
Potential Temperature Gradient (°K/m)
Stability A
Stability B
Stability C
Stability D
Stability E
Stability F
Reference Height for Wind Monitor
Surface Roughness Height
-0.010
-0.007
-0.005
0.000
0.020
0.035
10 m
0.03 m
As depicted in Table 4.2, all of the candidate escape fraction
equations exhibit smaller escape fractions for stable conditions than for
unstable and neutral conditions, as would be expected. Alternative 2, based
on equations (22)-(35) and using equation (13) to compute the escape fraction,
demonstrates somewhat better agreement with escape fractions inferred from the
source depletion model than do the other alternatives.
-45-
-------
TABLE 4.2
ESCAPE FRACTIONS BY STABILITY CLASSa
Equation
Equation
Universal Size Distribution
Stability Class
EDS Size Distribution
Stability Class
Evaluation Data:
Escape Velocity
Source Depletion
Residence Time
Theoretical Formulae:
Original Winges
Alternative 1
Alternative 2
Alternative 3
Alternative 4
1.00
0.93
0.94
0.99
0.40
0.85
0.22
0.78
1.00
0.88
0.94
0.98
0.45
0.75
0.28
0.70
1.00
0.86
0.96
0.96
0.58
0.70
0.42
0.66
1.00
0.81
0.95
0.92
0.65
0.70
0.48
0.65
1.00
0.58
0.91
0.58
0.35
0.11
0.24
0.11
_£_
Evaluation Data:
Escape Velocity
Source Depletion
Residence Time
Theoretical Formulae:
Original Winges
Alternative 1
Alternative 2
Alternative 3
Alternative 4
0.81
0.59
0.59
0.90
0.23
0.47
0.10
0.38
0.85
0.46
0.63
0.84
0.25
0.34
0.13
0.28
0.93
0.43
0.70
0.73
0.37
0.30
0.22
0.26
0.90
0.36
0.68
0.59
0.44
0.30
0.27
0.25
0.70
0.21
0.51
0.20
0.17
0.03
0.11
0.02
a Escape fractions were not computed for P-G stability E due to the
infrequent occurrence of this stability class (EPA, 1985).
-46-
-------
A similar comparison, stratified by wind speed instead of stability is
shown in Table 4.3. All of the candidate escape fraction equations show a
much greater change in escape fraction with wind speed than does the original
Winges equation. The increase in predicted escape fraction with wind speed
matches the trend observed in the evaluation data. In this sense, the
introduction of wind speed into the escape fraction computation is successful.
However, the overall conclusion made from examining all of the candidate
equations' predicted escape fractions, stratified by both wind speed and
stability class, is that none of the candidate escape fraction equations match
the evaluation data very closely.
The reasons for the discrepancies are not known; however, it is likely
that a number of effects contribute to the error in prediction. Included
among these effects are the assumption that there is no vertical component to
the wind. It is possible that the vertical component of the wind is
responsible for considerably more direct transport of the smoke puff out of
the pit than turbulent diffusion. Another source of error is the
interpretation of the measurement data. Although the computation of the
escape fraction in the evaluation data sets is based on the measurement of
residence time for a smoke plume, it may not be possible to infer one from the
other. Only by actual measurement of particulate release data could such a
quantification be made.
Additional attempts were made to examine the degree of agreement of the
candidate equations with the evaluation data using linear regression. It is
not possible to perform a linear regression for one of the evaluation data
sets (the escape velocity techniques with the Universal Size Distribution)
because this technique yielded a value of 1.0 for the escape fraction in every
one of the puff releases experiments. Linear regressions were performed,
however, for the residence time evaluation data set and for the escape
velocity evaluation data set using the EDS particle size distribution. The
results of the linear regression using the escape velocity evaluation data set
(with the EDS particle size distribution) are depicted in Table 4.4. As the
table shows, the observed and predicted comparisons in all cases revealed a
-47-
-------
TABLE 4.3
ESCAPE FRACTIONS BY WIND SPEED CLASS*
Universal Size Distribution
Equation
Wind Speed Category
Evaluation Data:
Escape Velocity
Source Depletion
Residence Time
Theoretical Formulae:
Original Winges
Alternative 1
Alternative 2
Alternative 3
Alternative 4
1.00
1.00
0.78
0.92
0.90
0.33
0.65
0.20
0.61
1.00
1.00
0.84
0.94
0.91
0.49
0.61
0.32
0.56
1.00
1.00
0.86
0.97
0.95
0.59
0.68
0.43
0.64
1.00
1.00
0.88
0.97
0.95
0.72
0.78
0.54
0.72
1.00
0.88
0.97
0.96
0.80
0.85
0.61
,0.77
EDS size Distribution
Equation
Wind Speed Category
Evaluation Data:
Escape Velocity
Source Depletion
Residence Time
Theoretical Formulae:
Original Winges
Alternative 1
Alternative 2
Alternative 3
Alternative 4
a Wind speed categories
defined as follows:
0.75 0.83 0.96
0.35 0.46 0.43
0.54 0.62 0.73
0.70 0.70 0.73
0.16 0.27 0.36
0.31 0.25 0.27
0.08 0.15 0.22
0.27 0.20 0.24
0.96
0.43
0.73
0.69
0.53
0.36
0.31
0.30
are those used by the National Climatic
Category 1: 0-3 knots
2: 4-6
3: 7 -10
4: 11 -16
5: 17 -21
6: above 21
0.99
0.43
0.76
0.76
0.66
0.45
0.38
0.34
Data Center,
-48-
-------
TABLE 4.4
LINEAR REGRESSION RESULTS FOR THE ESCAPE VELOCITY
EVALUATION DATA SET AND THE EDS PARTICLE SIZE DISTRIBUTION
Original Winges Equation
Alternative 1
Alternative 2
Alternative 3
Alternative 4
a
0.73
0.65
0.76
0.59
0.67
Regression Parameters
b
0.17
0.53
0.29
1.02
0.68
r2
0.03
0.23
0.04
0.39
0.11
very low correlation. This implies that it will be extremely difficult to
improve the prediction accuracy by adjustment of the theoretical formulae with
arbitrary constants.
Similarly, linear regressions were performed with the residence time
evaluation data set and each of the candidate equations for both the Universal
and the EDS size distributions. The results are depicted in Table 4.5. As
with the earlier table, the agreement between measured and predicted is not
encouraging.
-49-
-------
TABLE 4.5
LINEAR REGRESSION RESULTS FOR THE
RESIDENCE TIME EVALUATION DATA SET
Universal Size Distribution
Regression Parameters
a b r2
Original Winges Equation
Alternative 1
Alternative 2
Alternative 3
Alternative 4
0.86
0.90
0.92
0.90
0.91
0.10
0.09
0.04
1.15
0.06
EDS Size Distribution
Regression Parameters
a b
Original Winges Equation
Alternative 1
Alternative 2
Alternative 3
Alternative 4
0.57
0.52
0.63
0.48
0.59
0.12
0.41
0.07
0.94
0.26
0.11
0.22
0.07
0.34
0.11
r2
0.03
0.22
0.00
0.39
0.03
-50-
-------
5.0 ESCAPE FRACTION ALGORITHM FOR ISC
The previous sections have detailed the development of four separate
equations for computing the escape fraction as a function of commonly
measurable parameters. This chapter discusses the adaptation of these
equations into one of the standard air pollution models, the Industrial Source
Complex Model (ISC). Subroutines are developed for each of the four
techniques in the previous section for incorporation into the ISC Short-Term
Model (ISCST). In addition, it was necessary to make certain changes to the
main section of the program and two existing subroutines, INCHK and MODEL.
Appendix C contains a listing of the complete program as modified for
this purpose. The version of the ISCST Model that is shown in Appendix C is
identical to the version currently available from the National Technical
Information Service (NTIS) in the UNAMAP series, with a few changes. These
changes are listed as follows:
1. The version of ISCST in Appendix C has been adapted to run on an
IBM-PC Computer. The changes necessary to accomplish this were very
minor. OPEN statements were added and the character strings were
explicitly declared. Also, all quotation marks were changed to
apostrophes.
2. A new subroutine called ESCAPE for computation of the escape
fraction was added. In Appendix C there are four separate versions
of this subroutine corresponding to the four separate techniques
developed for the escape fraction computation in the previous
chapter.
3. The addition of the subroutine ESCAPE required several new user
inputs which were added to the ISCST main program and the
subroutines INCHK and MODEL through the addition of a new COMMON
block called DEPO. The new parementers were deposition velocity
array (a separate deposition velocity for each source and each
particle size group within each source - very similar to the way
gravitational settling velocities are included in the original
code), the reference height, ZREF, the surface roughness height, ZO,
and a pit depth for each source (incorporated in one of the
previously unused storage spaces in the SOURCE array), Changes were
made to the ISCST program to allow for user input of these
variables. ZREF and ZO were added to the end of the card group 2,
card number 2. The pit depth is read for each source at the end of
card group 6, card number 1 (after the building height) and the
deposition velocities are read as a new card appearing after card
group 6, card number 4.
-51-
-------
4. The call to the new subroutine was added to the MODEL subroutine at
two places: in the loop over particle size classes for the
concentration calculation and in the loop over particle size classes
in the deposition calculation. The subroutine ESCAPE returns the
values of the escape fraction, ESCP, which is then used to reduce
the vertical distribution function, V.
Using each separate version of the new subroutine, ESCAPE, four separate
versions of a compiled and linked ISCST Model were made and tested with a
sample data set to verify that they were operational. Appendix D presents the
sample outputs for each of these test data sets. The input file is also shown
in Appendix D. The same input file is used for all four versions of the model.
Run times for the four different versions of the model were recorded
during the tests. While the absolute values of the runtime is not of interest
here, the relative times are significant. As might be expected, the equations
using the more detailed analysis technique required more computer processing
time. The two techniques based on the linear model (alternatives 1 and 3)
required approximately the same processing time as the original version of
ISCST. Alternative 2 (Constant-K, detailed model) increased the run time by
roughly a factor of 1.5, while Alternative 4 (Variable-K, detailed model)
increased the run time by roughly a factor of 5.
-52-
-------
6.0 CONCLUSIONS AND RECOMMENDATIONS
6.1 CONCLUSIONS
A number of conclusions can be formulated based on the foregoing
analysis. Some of these conclusions have been stated in the previous report
(EPA, 1985) and will be restated here for completeness. Other conclusions
presented here have not been previously stated and will be discussed in more
detail.
It is clear from the analysis of stabilities discussed in chapter 3.0
that the standard deviation of the wind direction measured inside the pit is
always larger than that measured outside the pit. This suggests that the
horizontal turbulence is greater inside the pit than outside, and a reasonable
explanation for this would be that sensors inside the pit respond to
mechanical turbulence caused by airflow over, and in the wake of, the pit
wall. Outside the pit this mechanically induced turbulence is absent. It
must be remembered that the measurement of sigma-theta in the pit says nothing
about vertical turbulence inside or outside the pit.
Another observation concerning the stabilities calculated from the
standard deviation of the wind direction is that they agree well with the
stabilities estimated from the Pasquill-Gifford method which uses cloud cover
and ceiling height. The stabilities computed in and out of the pits are
either identical, or only one class different from, the P-G stability for
about 80% of the data hours.
Four new equations were developed for the computation of the escape
fraction as a function of commonly measured parameters. When compared to
escape fractions inferred from the smoke release data, none of these new
equations were seen to provide accurate predictions of the escape fraction
over the full range of stability classes and wind speeds. There are many
reasons for this descrepancy between measured and predicted values, but it is
important to reiterate that the smoke release experiments did not measure all
of the important quantities which define the escape fraction in and out of the
pit.
-53-
-------
Another possible reason for the discrepancy between measured and
predicted escape fractions is that none of the techniques used to calculate
the escape fraction considered the vertical motion of the wind (called
convection). It is clear from the smoke release video tapes that in many of
the experiments, the smoke was moved from the pit by convection of the air
rather than by dispersion. In all four of the escape fraction analysis
techniques developed here, an essential assumption is that the only mechanism
for transfer of material from the pit to the external air is by dispersion
convection is ignored.
The question of the influence of stability is very important. The smoke
release data imply that during unstable and neutral conditions a large
percentage of the dust emitted in the pit escapes. The various theoretical
escape fraction equations disagree with the inferred escape fractions from the
smoke releases for unstable and neutral conditions. It is likely that for
these conditions, the escape fractions inferred from the measurement data are
more accurate than the theoretically calculated escape fractions, because the
vertical motion of the air may be quite significant for the unstable and
neutral conditions and the theoretical formulae do not consider such motion
while the residence time extracted from the video tapes is influenced by the
vertical winds. Although analysis of vertical wind speeds might be
instructive in these instances, characterization of the chaotic flow in the
pit would be extremely difficult.
For stable cases, however, the situation is quite different. Here again
the smoke release data infer that a large percentage of the emitted dust
escapes the pit, but all the theoretical formulae conclude that only a small
fraction of the emitted dust escapes. For stable conditions one would expect
very little vertical motion of the air, thus the primary mechanism for
vertical transport should be dispersion, and in fact that is precisely the
assumption made in the development of the candidate equations. While
confidence increases in the candidate equations for stable conditions,
confidence in the experimental data, and in particular in the ability to infer
the escape fraction from the smoke release video tapes, decreases for the
stable cases. The reason for this is that in both the interpretations of the
smoke release data (the escape velocity evaluation data set and the residence
time evaluation data set) the escape fraction is computed from the ratio of
-54-
-------
the pit depth to the residence time a quantity called the escape velocity.
Vertical motion is implicit in both evaluations, when in fact for stable
conditions, there may be no vertical motion at all, and the residence time of
the puff in the pit may be as long as the stable conditions persist. Material
will leave the pit, but the mechanism is by dispersion, not by a vertical
escape velocity. The centerpoint of the puff (the point of maximum
concentration) remains in the pit. The smoke release video tapes did not
allow for such long residence times, and in fact the residence time in many
cases where the puff appeared to disperse in the pit without any vertical
motion was arbitrarily defined based on the time when the camera was turned
off, or when the puff was no longer visible.
Focus on the stable cases here is appropriate because they are the most
important cases to consider. Computer modeling studies done for permitting of
surface mining operations typically predict the peak concentrations under
stable, low-wind-speed conditions. The inferred escape fractions from the
smoke release data imply that a large percentage of the dust escapes the pit
during these conditions, while the four candidate equations predict low escape
percentages. Since the ability to infer the escape fraction from the smoke
release data is the least reliable for the stable conditions, and since the
assumptions made to develop the theoretical formulae are most representative
of the stable conditions, we conclude that the theoretical formulae are likely
to be more correct for these stable conditions than the escape fractions
inferred from the smoke releases.
Selecting between the four theoretical formulae for calculation of the
escape fraction is not an easy task. None of the equations work particularly
well for unstable and neutral conditions, and for the most important
conditions, the stable conditions, the evaluation data are suspect and do not
provide reliable selection criteria. In all the techniques the escape
fraction is defined by the amount of mixing in the pit which allows emissions
at the surface of the pit to be mixed upward to where the external flow of
wind can carry them away. The amount of mixing is characterized by the eddy
diffusivity. For the two models based on the linear model, the amount of
mixing is determined from the shearing of the wind speed profile caused by the
-55-
-------
surface drag of the earth a reasonable assumption for small scale
dispersion over a uniform flat plat in the open boundary layer. The other two
techniques, called the more detailed models, allow consideration of the
stability of the atmosphere as it affects the vertical mixing in the pit.
It is our conclusion, therefore, that one of the more detailed
techniques (Alternatives 2 and 4) should offer improved prediction accuracy
for the escape fraction for stable conditions. It is also evident from the
data evaluation that for particles smaller than 30 microns (the universal
particle size distribution) there is little difference between Alternative 2
and Alternative 4 (Constant-K vs Variable-K). We conclude that Alternative 2
is the most reasonable method to use for the escape fraction computation for
stable conditions for several reasons. There is no basis on which to
postulate a relationship for the eddy diffusivity with height within a mining
pit and the assumption of a constant value is the simplest assumption which
can be made. Also, the Variable-K method (Alternative 4) significantly
increases the computation time when used in the ISC Model, without providing
significantly different values than the Constant-K method (Alternative 2).
6.2 RECOMMENDATIONS FOR FUTURE WORK
The attempt to develop a characterization for the escape fraction for
mining pits is made difficult by the complexity of the dispersion scenario and
the difficulty in collecting meaningful data. In the course of our work, we
have identified a number of shortcomings in the data and the analysis
techniques and we have considered numerous methods for overcoming these
shortcomings. However, we believe that our recommendations should be guided
by the practical applications of the analysis techniques that are to be
developed. Consequently, we will not attempt to provide recommendations for
the resolution of every uncertainty we have identified in the course of our
study.
While it would be interesting to determine the vertical motions inside
mining pits that result in the escape of dust via vertical winds during
neutral and unstable conditions, it is very likely that such motions are too
-56-
-------
complicated to be predicted and simulated in an operational air quality
model. The ultimate use of the air quality model is for permitting purposes,
and the most important consideration is the conditions which produce the peak
impacts. Invariably it is stable low-wind speed conditions that result in the
peak concentrations. We will therefore not recommend such studies as wind
tunnel investigations of the flow in and around mining pits, because such
studies cannot yet adequately reproduce all atmospheric stability classes. We
would also not recommend further equation development to take vertical wind
velocities into account, since velocities are likely to be highly dependent on
site specific conditions, which would not be known prior to a mine's
construction, and ultimately it is not such vertical wind conditions which
lead to peak impacts as predicted by the air quality model.
Since we have developed several equations for prediction of the escape
fraction, which are hoped to work best for the stable conditions of most
concern, our primary recommendation concerns the need to develop a data base
for the escape fraction during stable conditions which can be used to evaluate
the theoretical formulas. One of the fundamental problems with the smoke
release data was that they were collected entirely during daylight hours, when
stable conditions rarely exist. The best technique for measuring the escape
velocity would be the use of dual tracer experiments where a tracer gas is
emitted simultaneously with a tracer particle, such as zinc-sulfide. By
measuring the relative concentrations of the two tracers downwind, the amount
deposited can be determined and the escape fraction readily determined. There
are problems with such techniques, because re-entrainment of the tracer
particles from previous experiments could contaminate later experiments, thus
restricting the number of experiments which could be run at a single mining
operation. The dual tracers could be run at night, however, when the stable
conditions are most persistent. We have not attempted to estimate costs for
such a program but it is assumed that such a program would be a relatively
high cost option.
-57-
-------
An alternative program which would be less costly, and would not suffer
the problems associated with the re-entrainment of tracer particles is to
perform a series of experiments with a single tracer gas such as
sulfur-hexafluoride. The single tracer experiments could still be run at
night when the stable conditions persist, and would provide a direct measure
of the ground-level concentration from pollutants emitted at the surface of
the pit. While such tracer gases would not undergo deposition, the knowledge
of ground-level concentrations at a grid of points throughout the pit as a
function of time after release would allow us to fully quantify the dispersion
process in the pit. With some assumptions concerning the deposition
velocities, the deposition rates of particles on the surface of the pit could
be inferred (with much greater accuracy than the smoke releases) and the
escape fraction determined. The disadvantage to this technique is that it
would require measurement of the tracer concentration at a large number of
locations in the pit and it would still involve assumptions concerning the
deposition velocities, which would be directly measured in the dual-tracer
experiments described earlier. The single tracer program would be
significantly cheaper than the dual tracer program, because the sampling
equipment for tracer gases is typically low-cost gas syringes which can be
analyzed in a remote laboratory with a gas chromatograph.
Another option for establishing a data base for evaluation of the escape
fraction equations during stable conditions is to re-evaluate the video-tape
recordings. Although the bulk of the experiments were in unstable or neutral
conditions, there were roughly 60 experiments during stable conditions. At
the time the original viewing of these tapes was performed, the viewers did
not know the stability. If given the opportunity to examine these tapes
again, the limited number of tapes and the knowledge of the stability class
might allow the evaluator to more accurately determine the residence time in
the pit. The particular items desired would be the trajectory of the puff and
the amount of surface contact experienced by the puff. Also those cases where
the puff stayed in the pit and dispersed will be evaluated using a much longer
residence time than previously used. It is possible that by reviewing the
tapes a more accurate representation of the escape fraction for the stable
cases can be established. If a new data base for evaluation of the escape
fraction equations can be developed, the equations can be evaluated with the
same technique used in the current study.
-58-
-------
A final option for future work would be a very simple investigation to
determine a "ballpark" estimate of the magnitude of pit retention. Using
existing hi-vol data and meteorological data already collected in the vicinity
of surface mines, a comparison can be made of actual measured concentrations
just downwind of a pit (C ,), and modeled concentrations determined
measured
from the ISCST model (C , - ,), which idealizes the terrain as flat and
modeled
unaffected by the presence of the pit. Emission rates could be estimated from
AP-42, Supplement 14 fugitive dust factors, and a representative background
concentration (perhaps from an upwind hi-vol) would be subtracted from the
measured concentrations. Any departure in the value of
(C ,/C , , .) from 1.0 would be due to errors in the emission
measured modeled
factors, or to errors in the model. If a long time period is considered
perhaps by examining annual average concentrations then random errors in
the model and emission factors will cancel out. Difference in the value of
(Cmeasured/Cmodeled) from unity would be due to systematic errors, such as
pit retention or plume perturbation caused by the pit. In the absence of
systematic errors in the emission factors or in idealizing the dust plume, the
ratio of (cmeasured/cmo(ieled) would be just equal to the escape fraction
for the particle size distribution collected by the hi-vols. This approach
would be a "first-cut" at estimating the magnitude of pit retention. Of
course, it would offer no insight into the physical mechanisms that control
dispersion from the pit, but it would provide an evaluation of the performance
of the emission factors and the ISC dispersion model. A study of this sort,
using existing data, would cost from $10,000 to $20,000.
-59-
-------
-------
REFERENCES
Air Sciences, 1985, letter from R. G. Steen, Principal, Air Sciences, Inc., to
C. F. Cole, TRC Environmental Consultants, Inc., April 17, 1985.
Bird, R. B., W. E. Stewart, and E. N. Lightfoot, 1960, Transport Phenomena.
J. Wiley & Sons, New York.
Chamberlain, A. C. and R. C. Chadwick, 1953, "Deposition of Airborne
Radioiodine Vapor," Nucleonics 2, 22-25.
Draxler, R. R., 1979, "Estimating Vertical Diffusion from Routine
Meteorological Tower Measurements," Atmospheric Environment. Vol 13, pp
1559-1564.
EPA, 1984, "Guideline on Air Quality Models (Revised): Draft," U.S. EPA,
OAQPS, Research Triangle Park, NC, November 1984.
EPA, 1985, "Dispersion of Airborne Particulates in Surface Coal Mines
Data Analysis," EPA-450/4-85-001, prepared for U.S. EPA, OAQPS, Research
Triangle Park, NC, January 1985 (NTIS No. PB 85-185411).
Fabrick, A. J., 1982, "Technical Note: Calculation of the Effective Emissions
from Mine Pit Operations by Incorporating Particulate Deposition in the
Evacuated Pit," MEF Environmental, Del Mar, CA, 1982.
Hittman and Air Sciences, 1983, "Studies Related to Retention of Airborne
Particulates in Coal Mine Pits Data Collection Phase," prepared for
U.S. EPA, IERL, Cincinnati, Ohio, contract #68-03-3037, August 1983.
Mendenhall, W., 1968, Introduction to Probability and Statistics. 3rd ed.,
Duxbury Press, Belmont, CA., January 1968.
Monin, A. S. and A. M. Yaglom, 1971, Statistical Fluid Mechanics, The MIT
Press, Cambridge, MA.
PEDCo & TRC, 1982, "Characterization of PM-10 and TSP Air Quality Around
Western Surface Coal Mines," prepared for EPA, Air Management Technology
Branch, contract #68-02-3512, June 1982.
Shearer, D. L., et al, 1981, "Coal Mining Emission Factor Development Study,"
prepared by TRC Environmental Consultants, Inc., 0908-D10-15, Englewood,
CO, July 1981.
Turner, D. B., 1970, "Workbook of Atmospheric Dispersion Estimates," U.S. EPA,
OAQPS, AP-26, Research Triangle Park, NC, 1970.
Van der Hoven, I., 1968, "Deposition of Particles and Gases," appearing in
Meteorology and Atomic Energy 1968, ed. D. H. Slade, Technical
Information Center, U.S. DOE, TID-24190.
Winges, K. D., 1981, "Description of the ERTEC Mining Air Quality (EMAQ)
Model, "ERTEC Northwest, Inc., Seattle, WA.
- 61 -
-------
APPENDIX A
AIR SCIENCES AUDIT REPORT
A-l
-------
AIR SCIENCES INC.
12687 West Cedar Drive
Lakewood, Colorado 80228
303/988-2960
April 17, 1985
Project No. 5-2
TRC Environmental Consultants, Inc.
7002 South Revere Parkway
Englewood, CO 80112
Attn: Mr. Cliff Cole
Senior Project Manager
RE: Sigma theta audit for TRC project 2990-V82
Gentlemen:
Air Sciences has performed an audit of the wind direction
standard deviation signal generated for EPA in the summer of 1983.
The audit included an electronic evaluation of the speed and
direction circuits used in generating the input for the sigma
calculation, and a checking of the software aspects of the
calculation. The audit confirms the data as approximately correct
as presented in the August 1983 report to EPA.
The direction deviation was calculated on-site by a Campbell
Scientific CR-21 data logger taking instantaneous wind speed and
direction data from the meteorological sensors. Samples of speed
and direction data were taken every ten seconds and from these one
minute averages were calculated. Thus, six instantaneous values
make up each minute calculation. Note that wind speed and direction
vector average was also calculated by the data logger as one-minute
averages and any electronic error arising from the sensors, signal
conditioning or logger input programming would also affect the wind
speed and direction vector averages which were provided in the 1983
report to EPA.
There are several points in the signal conditioning and
calculation process where error could arise and a list of them is
given below.
1 calibration error in direction sensor
2 excess friction in speed sensor
3 calibration error in electronic conditioning for direction
4 calibration error in electronic conditioning for speed
5 improper matching of the output of conditioning cards to
input of logging unit
A-3
-------
6 incorrect algorithm built into the digital processing unit
7 incorrect field programming of inputs
8 incorrect field programming of outputs
9 error in transferring data from field tapes to archive tape
Each of these nine items has been investigated as a part of the
audit.
1-4 The speed sensor, direction sensor, signal conditioning
circuitry, and data logger in the pit were identical to those out of
the pit. All calibrations and alterations made to the in-pit
sensors and circuits were made also to the out-of-pit sensors. The
sensors and signal conditioning circuits from the out-of-pit station
have not been used since the 1983 study. These components were
recalibrated as a part of this audit and the calibrations compared
with their 1983 calibrations. The in-pit sensors are not available,
but because of the similarity of the in-pit and out-of-pit systems,
we consider a thorough study of the out-of-pit sensors sufficient to
demonstrate the condition of the in-pit system also. The
calibration documentation is attached. The comparison shows that
the components are still in calibration. The checking of the
friction of the speed sensor is not documented on the forms. It
was checked by touch of an experienced technician and no excessive
friction was detected.
5 The speed and direction conditioning card outputs were
checked and found to be in the range of 0 to 1 VDC as was earlier
assumed. The conditioning cards had been altered in 1983 to produce
an instantaneous signal output rather than an averaged signal. This
alteration was checked and found to be correct. The logger units
were programmed to accept 0 to 1 VDC inputs as designated by the
input program No. 1 (as shown on the logger documentation form,
input programming section). Logger program No. 1 scales the data to
engineering units by a linear equation. That equation requires a
slope and an intercept. Note from the programming list that after
the input program number the slope and intercept are given. Slope
for speed is 50 mph/VDC and for direction is 540 degrees/VDC. The
intercepts are both zero by default since they were not programmed
in. Thus, the logger was receiving the data in the proper units and
performing the proper scaling.
6 The wind direction standard deviation calculational
algorithm is attached. It is a numerical procedure for estimating
direction standard deviation with listed error of less than 1
percent for deviations less than 40 degrees, which is well within
the precision of the measurement. Because direction is a circular
function rather than a scalar, the exact mathematical formulation is
lengthy and the algorithm in the data logger is only an
AIR SCIENCES INC.
A-4
-------
approximation. It is based on the assumption that there is no
correlation of speed deviation with direction deviation (page B-9).
This assumption has been experimentally verified under certain
conditions as stated with the algorithm explanation. It is possible
that with 10-second scans making up the rather short one-minute
averages in the EPA program that this assumption may lead to some
error, but we suspect that with several minutes of data averaged
together random error of the type we are addressing will cancel.
The standard deviation routine calculates a variance by
dividing by (N) rather than (N-l). This introduces an error
(underestimate) of about 10 percent in the EPA application where the
standard deviation was composed of only six values.
7-8 The programming of inputs and outputs has been documented
in the report to EPA. These steps have been verified with the
programming manual and found to be correct. Whether these steps
were followed in the field cannot be checked, but since the speed,
direction and deviation data appears to be consistent among sites we
assume that no mistakes were made in the field.
9 Data were collected in the field on cassette tapes and
transferred to other magnetic media in the office. It is
conceivable that in this transferral process a column of field data
could have been truncated. The data is logged in engineering units
in the field and this truncation would not have affected the
location of the decimal point. The data transferral program (a
trivial program) cannot be located and rechecked, but the data has
been studied and truncation error is not apparent. Only truncation
of the left column would be of concern to us and if the left-most
column were truncated it could only be the hundreds column. Most
sigma data is in the range of 0 to 40 degrees, well below the
hundred level.
These steps complete the Air Sciences audit of the sigma theta
calculation. We will be happy to answer questions that should arise
from this audit.
Sincerely^
Rodger G. Steen
Principal
AIR SCIENCES INC.
-------
APPENDIX B
HOURLY METEOROLOGICAL DATA BASE
B-l
-------
-------
In this Appendix, hourly averaged values of the parameters defined in
Table B2.1 are shown for each hour of the meteorological data base.
TABLE B2.1
DEFINITION OF VARIABLES
NAME
MEANING
DAY
NTIME
WDOUT
WSOUT
SGOUT
ISGOUT
WDIN
WSIN
SGIN
ISGIN
RAT
IPG
IWS
JULIAN DAY ON WHICH DATA WERE RECORDED
TIME AT END OF HOURLY AVERAGE (HHMM)
OUT-OF-PIT AVERAGE WIND DIRECTION
OUT-OF-PIT AVERAGE WIND SPEED (mph)
OUT-OF-PIT AVERAGE SIGMA-THETA (deg)
OUT-OF-PIT SIGMA-THETA STABILITY CLASS
IN-PIT AVERAGE WIND DIRECTION
IN-PIT AVERAGE WIND SPEED
IN-PIT AVERAGE SIGMA-THETA (deg)
IN-PIT SIGMA-THETA STABILITY CLASS
SGOUT/SGIN
PASQUILL-GIFFORD STABILITY CLASS
WIND SPEED CLASS
B-3
-------
DAY
NTIME WDOLJT WSOUT SGQUT ISGOUT WDIN WSIN SGIN
(DEB) (MPH) (DEB) (DEB) (MPH) (DEG)
ISGIN RAT
IRQ
179 10
179 11
179 12
179 13
ISO B
180 9
180 10
180 11
ISO 12
ISO 13
131 B
isi 9
131 10
181 11
131 12
181 13
1S1 14
181 15
1S1 16
181 17
181 IS
181 19
181 20
181 21
181 22
1S1 23
182 0
182 1
1S2 2
182 3
182 4
182 5
182 6
132 7
182 B
182 9
182 10
182 11
182 12
182 13
1S2 14
182 15
182 16
182 17
182 18
182 19
182 20
182 21
182 22
182 23
183 0
183 1
1B3 2
183 3
183 4
183 5
183 6
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
1059
1159
1259
1359
859
959
1059
1159
1259
1359
859
959
1059
1159
1259
1359
1459
1559
1659
1759
1859
1959
2059
2159
2259
2359
59
159
259
359
459
559
659
759
B59
959
1059
1159
1259
1359
1459
1559
1659
1759
1B59
1959
2059
2159
2259
2359
59
159
259
359
459
559
659
164.3
15.6
11.7
349.9
15.3
351.2
354.4
359.0
329.2
320.8
257.3
288.5
268.6
273.8
273.4
272.3
280. 1
275.7
2B0.6
266.5
253.5
*-tTT **l
£.4.0 . ji.
196.9
199.7
203.5
203.2
209.3
206.0
96.7
139.0
154.7
239.4
218.6
320 . 9
277.6
224.6
341.S
286.5
257.6
248.8
231.0
212.8
261. 1
247.6
245.9
220.3
1B2.0
139.5
191.6
203.6
200.9
190.3
200.4
221.9
214.4
216.7
223.7
5. 1
6.2
6.3
7.2
3. 3
3.2
3.5
4. 1
6.4
10. a
4.3
9.0
15.0
16.0
17.3
13.7
14.3
17.0
15.6
11.6
11.6
6.9
6.0
8.8
12.5
14.5
15.4
14.0
3.3
2.9
4.0
3.7
5.0
2.1
2.6
4.3
5.5
10. 1
17.2
17.3
13.7
13. 1
12.3
8. S
8.6
10.7
5.9
4.9
6.0
9.4
11.4
9.3
9.7
13.9
15.4
15. 1
13.9
2O. 91
16. 12
12. 12
11.01
16.59
15.34
16.58
13.99
16.44
12.52
27.37
21.32
11.33
10. 15
9.80
1 1 . 50
11.49
9.50
8.94
10.58
9.42
7.20
1O.99
5.83
4. SB
5. 12
5.58
5.77
21.12
16.77
16.92
22.76
23. 33
23.66
30.88
27.09
23. S6
13,19
10.39
11.25
11.94
14. 10
12. 12
12.73
9. BO
9.53
11.67
11. 7B
17.95
14.22
7.19
15. 10
13.24
8.66
8.00
7.39
9. 14
2
3
4
4
3
3
3
3
3
3
1
2
4
4
4
4
4
'4
4
4
4
5
4
5
4
4
4
4
6
5
5
6
6
1
1
1
1
3
4
4
4
3
4
3
4
4
4
4
5
4
4
4
4
4
4
4
4
257.9
48.8
57.2
51.2
37.5
67.7
43.7
34.6
74.7
71.5
330. 6
313.3
228.9
226.6
224.6
234.3
226. 1
2^T1 O
-i*. O
215.3
228.2
238.4
253.8
253.7
247.2
270.7
274.8
271 = 7
270.7
312.0
176.2
175.3
285.2
265. B
303 . 8
334.5
10.7
123.9
224.7
234.9
241.3
262.6
276.8
237.3
24 1 . 0
241.5
263.7
2BO.B
179.9
242.5
265.0
28O.8
254.7
244.0
246.4
262.9
257.5
247.7
4. 1
7.5
6.6
7.2
3.2
3.5
4. 1
4.2
5.9
a. a
i- f-t
O Ğ ji.
32.31
18.52
26.03
32.43
29.85
33.93
13. 15
25.28
43.26
32.97
2.91
3.6999.90
11.2
10.3
13.2
10.9
11.5
12.5
12.6
9.1
8.6
5.7
4.6
4.9
8.3
10.4
12.3
10.8
5.5
1.7
2.0
a. 3
4.7
4.4
4.2
4.5
6.7
a. s
12. a
13.3
11.3
10. a
9.7
7.2
7.9
3.4
4.0
2.2
3.3
6.6
10. 1
5.7
4.6
9.4
10.7
10.6
9. 1
29.07
19.77
20.31
22.68
21.67
19.51
19.66
19.11
15.43
15.54
21.51
21.76
13.51
10. 13
9.44
11.32
25.20
29.78
31.74
18.41
30.37
26.62
31.20
29.35
36 . 09
27.37
16.02
15.54
19.44
21.68
19.48
20 . 25
14.25
15.35
23.14
25.74
26.43
21.62
11.63
30. 46
33. 18
15.50
15.67
12. 6O
46.77
1
"1
*~
1
1
1
1
2
1
1
1
4
9
1
2
2
1
2
2
2
^
4Ğ
3
4
6
6
4
4
4
4
6
6
6
4
6
1
1
1
1
1
3
j>
2
2
2
2
3
4
6
6
6
5
4
6
6
4
4
4
4
.65
.87
.47
.34
.56
.41
.91
.55
.33
. 33
9.42
9.99
.39
.51
.48
.51
.53
.49
.45
.55
.61
.46
.51
.27
.36
.51
.59
.49
.84
.56
.53
1.24
.77
.89
.99
.92
.66
.48
.65
.72
.61
.65
.62
.63
.69
.62
.50
.46
.68
.66
.62
.50
.40
.56
.51
.59
.20
2
2
1
2
2
2
*"^
jL
f^
fi
j-\
*Lğ
3
?
4
4
jj
3
3
4
4
3
4
6
6
5
4
4
4
4
6
6
6
6
6
2
2
i
2
o
3
3
3
3
3
-3*
3
5
6
6
6
5
5
5
5
4
4
4
4
B-4
------- An error occurred while trying to OCR this image.
-------
DAY
195 0
195 1
195 2
195 3
195 4
195 5
195 6
195 7
195 3
195 9
195 11
195 12
195 13
195 14
195 15
195 16
195 17
195 18
195 19
195 20
195 22
195 23
196 O
196 1
196 2
196 3
196 4
196 5
196 6
196 7
196 8
196 9
196 10
196 11
196 12
196 13
196 14
196 15
196 16
196 16
196 19
196 20
196 21
196 23
197 0
197 1
200 7
200 10
200 11
200 12
2OO 13
2OO 14
200 15
200 16
200 17
200 19
200 22
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
NT I ME
59
159
259
359
459
559
659
759
359
959
1159
1259
1359
1459
1559
1659
1759
1359
1959
2059
2259
2359
59
159
259
359
459
559
659
759
859
959
1059
1159
1259
1359
1459
1559
1659
1859
1959
2059
2159
2359
59
159
759
1059
1159
1259
1359
1459
1559
1659
1759
1959
2259
WDOUT
(DEG)
158.4
310.5
317.7
305. 7
319.3
325.0
30 1 . 3
146.9
151.4
157.0
157.2
160.9
147.3
132.9
191.8
197.8
207. 1
198.8
229.5
209.8
226.3
332.2
332.2
333.3
345.6
331.3
335. 6
343.6
341.7
335.6
340 . 4
357.0
8.1
346.0
332. 1
343.8
349 . 0
8.6
323.4
79.8
121.3
140.5
138.6
315.9
319.4
347.2
116.6
116.0
106.5
153.3
131.7
144.1
130.1
108.8
118.2
212.0
T-l f O
WSOUT SGDUT
(MPH)
-------
AY NTIME
WDQUT WSOUT
(DEB) (MPH)
201
201
20 1
201
201
201
201
201
201
201
201
201
201
201
201
201
201
20 1
201
202
202
202
202
202
202
202
202
202
202
202
202
202
202
202
202
202
202
202
202
203
203
203
203
203
203
/-if-sT
4.UO
1 59
2 59
3 59
4 59
6 59
7 59
a 59
9 59
10 59
11 59
12 59
13 59
14 59
15 59
16 59
17 59
18 59
19 59
20 59
0 59
1 59
2 59
4 59
5 59
6 59
7 59
8 59
9 59
10 59
11 59
12 59
13 59
14 59
15 59
16 59
17 59
19 59
21 59
23 59
1 59
3 59
5 59
6 59
7 59
8 59
9 59
159
259
359
459
659
759
859
959
1059
1159
1259
1359
1459
1559
1659
1759
1859
1959
2059
59
159
259
459
559
659
759
859
959
1059
1159
1259
1359
1459
1559
1659
1759
1959
2159
2359
159
359
559
659
759
B59
959
281.9
280.3
323.0
358.0
279.3
253.5
106. 1
153.8
66.2
11.9
9O.4
23.8
102.7
101.3
107.2
116.8
112.6
134.4
266.3
306 . 9
265.2
290.6
277.5
284.0
278.0
60.5
94.6
102.2
107.2
97.4
107.3
105.1
103.0
127.8
95.0
106.8
83.7
32.3
270.8
281.9
271.1
269.8
282.5
353. 1
307.7
299.5
6.6
6.2
3.8
3.9
4.2
3.2
2.3
4.6
3.0
2.3
2.8
3.6
8.6
14.0
10.7
6.0
6.5
7.5
15.8
7.4
4.7
5. 1
6.5
6.1
5. 1
3.5
3.8
6.0
5.5
6.0
7.6
10.3
9.7
10.5
a. i
7.0
9.2
3.1
4.6
3.4
6.2
4.8
5.5
2.0
7.0
6.4
SGQUT ISGOUT WDIN WSIN SSIN ISBIN
(DES)
9.22
7.95
18.78
12.32
9.66
8.46
1 1 . 8O
14.26
9.33
1 7 . 20
13.10
20.71
8.04
6.27
7.66
10.21
9.26
5.93
9.96
10.88
20 . 58
13.37
7.96
9.81
8.98
10.45
11.75
15. 11
23.64
20.79
20.52
15. 16
18.87
13.35
16.32
18.73
7.23
24. 15
7.57
12.79
3. 87
5.87
8.76
1 4 . 6O
14.66
11.71
4
4
6
4
4
4
4
3
4
3
3
~\
4
4
4
4
4
5
4
4
6
5
4
4
4
4
4
3
1
2
2
3
2
3
3
2
5
6
5
5
5
5
4
-r
-'
3
4
(DEB) U1PH> (DES)
274.8 4.4 17.80
273. 1 3.8 22.16
119.1 3.9 15.71
142.8 3.0 16.92
290.9 3.4 16.25
322.0 2.2. 14.01
113.7 2.1 12.84
99.7 4.1 16.43
106.4 3.0 11.64
142.9 1.8 24.21
122.6 2.8999.90
999.0999.0999.90
106.8 8.7999. 90
109.9 11.3 9.69
109.6 9.0 6.94
103.6 4.9 11.57
108.0 6.4 8.69
111.9 4.9 10.62
286.1 11.6 18.83
293.8 4.2 10.46
266.3 3.3 21.56
281.2 5.3 8.13
274.7 6.6 7.87
279.5 5.2 17.18
97.0 2.6 15.63
108.7 4.3 5.33
108.1 3.6 11.43
106.7 5.9 13.79
136.4 5.6999.90
270.0 6.6999.90
999.0999.0999.90
999.0999.0999.90
999.0999.0999.90
999.0999.0999.90
999.0999.0999.90
999.0999.0999.90
999.0999.0999.90
999.0999.0999.90
999.0999.0999.90
999.0999.0999.90
999. 0999. 0999. 9O
266.4 3.O 19.78
288.3 4.2 11.67
121.3 1.9 19.63
276.8 6.0 13.81
231.7 6.4 18.20
6
6
5
tr
vJ
5
3
.T
-'
4
i
9
9
9
4
4
4
4
4
4
4
6
4
4
5
5
4
4
J"
9
9
9
9 .
9
9
9
9
9
9
9
9
9
6
4
<-l
i.
_>
2
RA
.
.
i.
.
*
,
.
9.
9.
9.
m
1.
.
1.
ğ
m
1.
m
i.
l.
.
l.
1.
1.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
.
.
B
1.
i
T
52
36
20
73
59
60
91
87
BO
71
99
99
99
65
10
aa
07
.56
53
04
95
64
01
57
57
96
03
10
99
99
99
99
99
99
99
99
99
99
99
99
99
30
75
74
06
64
IPG IW3
6
6
6
6
6
*-l
jL.
ji.
J-l
T!
O
1
1
2
4
-i>
r?
--'
5
4
6
6
6
6
6
6
2
2
2
2
2
2
2
2
;r
.*.
*-ğ
^.
5
6
6
6
6
6
6
2
T;
.
^.
2
1
1
*7
i
i
^_
i
i
i
i
"^
4
^'
_
/-I
*-l
4
*-ğ
.£.
-l
j-
1
2
2
.hi
i
*"?
.^
R-7
-------
DAY
NT I ME
213
213
213
213
213
213
213
213
213
213
213
213
214
214
214
214
214
214
214
214
214
214
214
214
214
214
214
214
214
214
214
214
214
214
214
214
215
215
215
215
215
215
215
215
215
215
215
215
215
215
215
215
215
215
215
215
215
215
215
215
11
12
14
15
16
17
18
19
20
21
22
23
0
1
2
3
4
5
6
7
8
9
1O
11
12
13
14
15
16
17
18
19
20
21
22
23
0
1
*?
3
4
5
6
7
3
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
39
59
59
59
59
59
59
59
59
59
59
59
59
59
1159
1259
1459
1559
1659
1759
1859
1959
2059
2159
2259
2359
59
159
259
359
459
559
659
759
859
959
1 059
1159
1259
1359
1459
1559
1659
1759
1859
1959
2059
2159
2259
2359
59
159
259
359
459
559
659
759
359
959
1059
1159
1259
1359
1459
1559
1659
1759
1859
1959
2059
2159
2259
2359
WDOUT WSO'JT
(DEG) (MPH)
999.0999.0
999.0999.0
246.3
153.0
ISO. 2
151.7
249.9
227.9
186.1
213.4
234. 1
215.4
216.0
228.6
232.0
231. 1
227.4
229.9
240. 1
244.4
246.6
270.8
264.0
297.8
283.3
271.3
286.9
290.9
280.7
199. 1
163.8
164.2
174.2
195.4
182.7
203.4
308.0
274.0
227.8
207 . 9
202.4
223.8
230 . 5
276. 1
266.7
237 . 7
266.6
259.9
246.2
251.9
228.2
187. 1
204.5
173.3
145.4
154.6
156. 1
164.9
214.7
112.3
13.7
15.0
12.9
11.4
8.7
7.9
7.2
10.8
11.6
12.6
13.7
11.7
11.2
10.6
9.3
10.6
9.0
10.4
10. 1
11.2
10.6
11.2
13.6
13.3
12.1
11.3
8.6
9. 1
19.0
17.6
14.3
11. 1
10.6
7. 1
4.4
6.2
3.6
3.5
5.7
4.8
3.7
6. 1
7.7
11.3
B.B
a. i
7.5
6.8
6.4
7.4
5.5
6.4
7.0
8.7
11.7
9.2
6.4
4.8
SGOUT
(DEG)
999 . 9O
999.
IS.
a.
8.
9.
IS.
14.
13.
17.
13.
15.
14.
15.
13.
13.
15.
11.
14.
15.
17.
16.
17.
12.
15.
IS.
15.
13.
14.
11.
11.
8.
14.
19.
14.
90
02
46
59
09
14
72
68
63
41
15
45
34
30
75
01
70
92
38
48
91
83
70
08
73
76
33
33
50
29
39
54
40
53
15.21
8.
1O.
13.
22.
23.
26.
22.
17.
18.
15.
18.
20.
23.
20.
21.
24.
25.
IS.
12.
9.
6.
8.
17.
11.
67
20
76
58
12
72
28
67
36
51
25
06
24
09
91
15
60
33
57
26
19
64
70
30
ISBQUT WDIN WSIN
(DEB) (MPH)
9 352.8 13.2
9
4
4
4
4
2
4
4
4
4
4
4
4
4
4
4
4
4
3
3
3
2
3
4
2
3
3
3
4
4
4
4
4
4
4
4
4
5
6
6
6
6
2
2
C'
2
2
1
2
2
1
1
2
3
4
4
4
5
4
345
341
105
114
115
47
7
95
346
345
337
355
1
358
357
357
1
359
349
353
357
351
316
335
344
331
320
332
94
111
114
108
S3
112
74
305
315
257
2O4
126
12S
123
354
4
354
353
356
346
359
90
SO
96
131
129
120
117
299
159
.7
. 1
.4
.5
.5
.5
. 3
.3
.5
.8
.7
.5
.6
.4
.4
.2
.9
. 1
. 3
. 0
.8
.6
.7
.5
.9
.4
.4
.7
.0
.9
.4
.4
.7
.7
.9
.9
.7
. 1
.9
.7
.9
.2
.7
.8
.4
.6
. 1
. 1
.1
.2
.2
.0
.4
.2
.0
.7
.0
. 1
*">
11.6
12.9
8.8
9.6
9. 1
7.3
tJ O
3.2
4.S
7.8
S. 7
7.8
8.9
7.4
7.6
7.3
7. 1
6.2
8-r
m -'
8.7
8.4
8.5
9.4
10.7
11.0
9.7
9.2
7. 1
5.2
11.7
9.5
8.2
4.9
6.7
4.2
3.7
4.O
2.0
2.9
3.9
4.4
2.9
4.6
6.8
10.0
7.0
6.5
6.4
5.7
5.5
5.2
4.7
4.8
5.3
5.4
5.9
5. 1
3.7
3.4
SO IN
(DEG)
24.69
21.
24.
20.
10.
10.
23.
27.
38.
29.
15.
17.
28.
21.
24.
23.
20.
21.
21.
18.
23.
29.
25.
14.
22.
23.
25.
20.
22.
31.
IS.
IS.
30.
41.
35.
37.
13.
13.
21.
19.
35.
28.
29.
17.
24.
21.
33.
25.
29.
46
56
47
74
94
72
45
09
65
60
97
42
30
71
23
03
40
57
16
39
35
96
47
11
12
66
16
14
57
41
67
31
SI
44
12
44
5O
06
49
26
48
66
45
77
76
23
56
42
25.87
30.
32.
29.
23.
16.
21.
21.
27.
27.
21.
34
21
61
94
14
30
84
08
73
57
ISGIN RAT
1 9.99'
*"ğ
-1.
1
2
4
4
1
6
6
6
4
4
5
4
5
5
4
4
5
2
1
1
1
3
2
1
1
2
2
1
2
4
4
6
5
6
5
5
6
6
6
6
6
?
-~t
1
2
1
1
1
1
1
1
1
1
3
6
5
6
6
6
9.99
.73
.41
.80
.33
.76
.54
.36
.59
.86
.84
.51
.70
. 54
.59
.75
.55
.69
.35
.75
.58
.69
.88
.63
.31
.61
.66
.65
.36
.61
.48
.48
.46
.41
.41
.65
.76
.65
1. 16
.66
.94
.75
1.01
.74
.71
.55
.78
.79
.78
.72
.75
.86
.77
.78
.43
.28
.32
.64
.55
IPS
9
9
3
4
6
5
5
4
4
5
.5
5
5
3
3
3
3
3
4
4
4
5
5
6
6
6
6
6
6
6
6
3
2
2
2
2
2
2
2
3
5
5
5
6
6
B-8
-------
DAY
NT I ME
WDQUT WSQUT SSQUT ISBQUT WDIN WSIN SBIN ISBIN RAT
(DEB) (MPH) (DES) -\
JL
2
1
1
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
6
7;
o
f^
B
B
.
.
B
.
.
.
.
.
.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
9.
n
m
m
m
59
37
62
53
60
59
56
s4
97
98
79
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
43
87
73
88
6
6
6
6
6
6
6
t
2
rj
*-
1
1
^_
3
T;
4
4
5
6
6
6
5
5
6
6
6
6
6
6
jfl
O
2
1
1
1
1
1
1
^
1
^,
.-}
*^k
*"1
-i.
2.
T
_<
3
4
4
4
.^i
2
1
1
3
T
H
2
1
1
1
1
1
*i
3
B-9
-------
APPENDIX C
FORTRAN LISTING OF MODIFIED
ISCST PROGRAM
C-l
-------
------- OCR error (C:\Conversion\JobRoot\0000068W\tiff\2000N80A.tif): Saving image to "C:\Conversion\JobRoot\0000068W\tiff\2000N80A.T$F.T$F" failed.
-------
40
50
60
15 = 14
16 = 15
17 = 15
IF (NXWYPT
16 = 15 +
17 = 16 +
18 = 17
IFUSWC4)
19 = IS
IFUSWU7)
I1O = 19
111 = 110
112 = 110
IF(ISW<1G>
I1O = 19 +
111 = I1O
112 = 111
113 = 112
+ NYPNTS
.ED. O)
NXWYPT
NXWYPT
GOTO 5O
NE. O) IB
17 + NPNTS
NE. O) 19 = IB
.LE. 0> BOTO 60
15O*NAVG*NGRQPS
5O*NAVG*NGROPS
NAVG*NGROPS
215*NSOURC - 1
4*NN
: DETERMINE IF CALCULATED STORAGE ALLOCATION EXCEEDS LIMIT.
IF(I13 .LE. LIMIT) BOTO 7O
WRITEUO,9004) 113, LIMIT
STOP
: CALL INPUT ROUTINE.
7O CALL INCHK ,QF(I1) ,QF(I2) ,QF
------- OCR error (C:\Conversion\JobRoot\0000068W\tiff\2000N80C.tif): Unspecified error
-------
3O
50
6O
70
BO
90
iOO
110
12O
130
140
ISO
16O
170
ISO
19O
200
LINE = IOO
DO 310 I = 1,NSDURC
ITYPE = SOURCE(1,1)
IWAK = ITYPE/8192
ITYPE = ITYPE - (ITYPE/16)*16
IF(ITYPE-l) 4O ,11O,14O
HB = SOURCE(11,I)
HW = SOURCE(12,I)
IF(HB .LE. 0.0 .AND. HW .LE. 0.0) GOTO 19O
H = HB
IF(HW -LT. HB) H = HW
DO 5O J = 1,36
SOURCE(Bl+J,I) = (1.2*H/SAS1BZ(J))*#**SQ(J)
SOURCE(Bl+J,I) =
(SIGZO/SASIGZ(J))**
(l./SBSI6Z(J))
.LT.
.LT.
O.O) SOURCE(75+J,I)
O.O) SOURCE(Bi+J,I)
21O
220
23O
GOTO 16O
XO = SOURCE(9,I)
DO ISO J = 1,36
SOURCE(Bl+J,I) = .OO1*XO
NO VIRTUAL DISTANCFS CAN BE LESS THAN ZERO.
DO 17O J = 1,6
IF(SOURCE(75+J,I)
DO 1BO J = 1,36
IF(SOURCE(Bl+J,I)
Al = 99.99
IF(ITYPE-l) 200
XOP =O.O
H = HB
IF(HW .LT.
Al = 3.*H
IF(A1 .LT.
GOTO 230
XOP = 2.15*SIGYO
GOTO 23O
XOP = .5641896*SOURCE(9,I)
NSO = SOURCE(2,I)
XS = SOURCE(4,I)
YS = SOURCE(5,I)
IF(NXPNTS .EQ. O
POLAR = .FALSE.
IF(ISW(2) .EQ. 2
0.0
O.O
HB) H
,210 ,220
HW
99.99) Al = 99.99
.OR. NYPNTS .EQ. O) GOTO 27O
.OR. ISW(2) .ED. 4) POLAR
.TRUE.
DO
YR
II
DO
XR
240
260 J = 1,NYPNTS
= GRIDY(J)
= YR
26O K = 1,NXPNTS
= GRIDX(K)
IF(.NOT.POLAR) GOTO
YR = XR*COSNUM(I1)
XR = XR*SINNUM(I1>
A3 = YR - YS
XR = XR - XS
240
S0300600
SO3OO610
SO3OO620
S0300630
S03OO640
SO3O065O
S0300660
SO30O670
S03006BO
S0300690
S03OO7OO
SO3OO71O
S0300720
S0300730
SO30O74O
S0300750
SO3OO76O
SO3OO770
SO3OO7BO
SO3OO79O
SO3OO8OO
S03O0810
SO3OO82O
S0300830
SO30O84O
S0300850
S03OO86O
SO3OO87O
SO3OOBBO
SO3OO89O
SO3OO90O
SO3O091O
S03OO920
SO3OO93O
SO30O94O
SO3OO95O
S0300960
SO30O970
S0300980
S030O99O
S0301000
SO3O1O1O
50301O2O
SO3O1O30
50301O4O
SO3O1O5O
SO301O6O
S03O1O7O
SO301O8O
SO301O9O
S0301100
S0301110
S0301120
SO3O113O
S030114O
S0301150
50301160
50301170
SO3011BO
S03O119O
SO3O120O
50301210
S030122O
SO3O123O
S0301240
5O3O1250
SO3O126O
50301270
SO3O128O
C-6
-------
25O
260
270
280
290
300
310
C
C
320
C
C***
C
A2 <= BQRT(XR*XR + A3*A3) - XOP
IF(A2 . GE. Al) GOTO 26O
IF(LINE .LT. 57) GOTO 25O
WRITE(IO,9O11)
WRITE(IO,9OO5) TITLE
WRITE(I0,9002) CDNDEP
LIME = 16
WRITE(10,9003) NSO,GRIDX(K),GRIDY(J),A2
LINE = LINE + 1
CONTINUE
IF(NXWYPT .ED. O) GOTO 310
POLAR = .FALSE.
IF(ISW(3) .ED. 2) POLAR = .TRUE.
DO 30O J = 1,NXWYPT
YR = YDIS(J)
XR = XDIS(J)
IF(.NOT.POLAR) GOTO 280
II = YR
YR = XR*COSNUM(I1)
XR = XR*SINNUM(I1)
YR = YR - YS
XR = XR - XS
A2 = BQRT(XR*XR + YR*YR) - XOP
IF(A2 .GE. AD GOTO 3OO
IF(LINE .LT. 57) GOTO 29O
WRITE(IO,9O05> TITLE
WRITE(IO,9OO2) CONDEP
LINE = 16
WRITE(10,9003) NSO,XDIS(J),YDIS(J),A2
LIME = LINE + 1
CONTINUE
CONTINUE
INITIALIZE NUMBER DAYS, HOURS & HOURS
INDEX.
NTDAY - O
IFUSWU9) .GT. 1) GOTO 320
NHQURS = 24
IHM = 1
IFCISW(20> -BT. 0) IHM = 2
PER DAY. SET MIXING HEIGHT
BEGIN LOOP OVER DAYS OF METEOROLOGICAL DATA.
DO 169O IDY = 1,NDAYS
WRITE(*,*> ' STARTED DAY NO.',IDY
IF(ISW(19) .ED. 1) GOTO 38O
C INPUT A DAY OF CARD MET DATA.
DO 370 I = 1,NHOURS
READ(IMET,90O4) JDAY,AFV(I),AWS(I),HLH(I,1),TEMP(I),DTHDZ(I>,
1 ISTAB(I),P(I),DECAY(I)
IF(ISTABd) .GT. 6) ISTAB(I) = 6
AFVR(I) = AFV(I)
IF(JDAY .LT. 1) JDAY = 1
IF(I.EQ.1) JDY=ODAY
IF(ISW(21) .ED. 3 .AND. ISW(22) .EQ. 3) GOTO 35O
C COMPUTE WIND SPEED CATEGORY IN ORDER TO LOAD DEFAULT VALUE FOR
C P OR DTHDZ.
1ST = ISTAB(I)
DO 33O J = 1,5
ISP = J
IF(UCATSCJ) .GE. AWS(I)) GOTO 34O
33O CONTINUE
ISP = 6
340 IF(ISW(21) .NE. 3) P(I) = PDEF(ISP,1ST)
IF(ISW(22) .NE. 3) DTHDZ(I) = DTHDEF(ISP,1ST)
35O IF(ISW(6) .NE. 2) GOTO 37O
IF(I .GT. 1) GOTO 360
WRITE(ID,9OO1) ODAY
S03O1290
S03O13OO
SO30131O
SO30132O
SO3O1330
S030134O
S0301350
50301360
S03O1370
SO3O13SO
S03O139O
5030140O
S0301410
S03O142O
SO3O143O
S030144O
S030145O
S030146O
S03O147O
503014BO
S0301490
S03O15OO
S0301510
S0301520
S03O153O
503O154O
SO30155O
SO3O156O
SO3O157O
S03O15BO
S03O159O
S0301600
S0301610
SO301620
S030163O
S0301640
S03O165O
S030166O
SO3O167O
S03016SO
SO30169O
S03O170O
S0301710
S030172O
SO3O173O
50301740
50301750
S03O176O
S030177O
SO3017BO
S03O179O
SO3O1795
S03O180O
S0301B1O
50301820
S0301B3O
SO301B4O
SO301B5O
SO30186O
SO3O187O
S0301B8O
SO3O1B9O
S03019OO
SO30191O
SO3O192O
50301940
C-7
-------
WRITE(ID,9005) TITLE S03O1950
WRITE(ID,90O7) JDAY SO30196O
WRITE(IO,9OO6) SO3O197O
36O WRITE(IO,90O8) I,AFV(I),AWS(I>,HLH(I,1>,TEMP(I),DTHDZ(I), SO3O19QO
1 ISTABU) ,P(I) ,DECAY(I) SO3O199O
37O CONTINUE . SO3O2OOO
LINE = O SO3O201O
6DTD 49O SO3O2O2O
INPUT PRE-PRDCESSED MET DATA. SO30203O
3BO IFCIDAYUDY) .ST. O) GOTO 410 S030204O
II = IDY + 1 SO3O2O5O
IFdDAYdl) .BT. O) GOTO 39O SO3O2O6O
READUMET) ISTAB SO3O2O7O
GOTO 1690 SO3O2OBO
390 READ(IHET) JYR,IMO,DAY,ISTAB SO3O209O
LSTAB = ISTAB(1) SO3021OO
IF = HLH(J,I) SO3O224O
DO 430 I = 1,48,2 SO3O2250
J = .5*1 + 1 SO302260
430 HLH(J,1) = RLHU) SO3O227O
DO 44O I = 2,49,2 S030228O
J = .5*1 S0302290
44O HLH(J,2) = RLH(I) . SO3023OO
IF(IDY . EQ. 1) LSTAB = ISTAB(l) S03O2310
IF(LSTAB .GT. 6) LSTAB = 6 S0302320
DO NOT ALLOW STABILITY TO VARY RAPIDLY & ADJUST FOR URBAN MODES. S03O233O
DO 46O I = 1,24 SO3O234O
IFdSTAB(I) .GT. 6) ISTAB ISTAB(I) = 4 SO30244O
460 LSTAB = ISTAB(I) SO30245O
IF(ISW(6) .NE. 2) GOTO 4BO SO30246O
WRITE(IO,90O1) IDY SO3O247O
WRITE(IO,9OO5) TITLE S03O248O
WRITE(I0,9007) IDY S03O2490
WRITE(I0,9009) SO3O25OO
DO 470 I = 1,24 S0302510
470 WRITE(IO,9O1O) I,AFV(I),AFVR(I),AWS(I),HLH(I,IHM),TEMP(I), SO30252O
1 MSTAB
-------
490
50O
C
c***
C
c
c
510
52O
530
540
550
560
C
C***
C
570
CONTINUE
IMO = 12
CONTINUE
ISEA = ISEAS(IMO)
BEGIN LOOP OVER MET DATA FOR EACH HOUR.
DO 167O IHR = 1,NHOURS
1ST = ISTABUHR)
IF URBAN MODE 2, ADJUST STABILITY FOR CALCULATION OF SIGY & SIGZ.
ISTUM2 = 1ST
IF(ISW(20) . EQ. 2) ISTUM2 = 1ST - 1
IFUSTUM2 .LT. 1) ISTUM2 = 1
UBAR = AWS(IHR)
FV = AFV(IHR)
FVR = AFVR(IHR)
HM = HLHUHR,IHn>
SET MIXING HEIGHT TO 1OOOO.O SO THAT ONLY FIRST TERM OF VERTICAL
EQUATION IS COMPUTED (RURAL MODE, E & F STABILITIES ONLY).
IF .EQ. 2) GOTO 530
PP = PDEF(ISP,IST)
DTH = DTHDEF(ISP,IST)
DECAY
-------
TS = SOURCE(8,IS) SO3O3310
NSO = SOURCE(2,IS) SO3O332O
IWAK = ITYPE/B192 S03O3330
QFLS = ITYPE/512 - (ITYPE/B192)*16 SO3O334O
NVS = ITYPE/16 - *24 + IHR S030355O
630 QTK = SOURCE(I2+119,ID SO3O356O
640 QTK = SOURCE(3,IS)*TK*QTK S03O357O
CALCULATE EFFECTIVE WIND SPEED. SO30358O
UBARS = UBAR SO3O359O
IF(PP) 670,670,650 SO3036OO
650 IF(HS) 670,670,660 SO30361O
NOTE7. ZR IS IN RECIPROCAL FORM. S03O362O
660 Al = HS S0303630
IFCHS .LT. 10.0) Al = AMIN1(1O.O,1./ZR) SO3O364O
UBARS = UBAR*(A1*ZR)**PP S03O365O
670 UBARI = 1./UBARS . S03O366O
BEBIN PLUME RISE CALCULATIONS FOR STACK-TYPE SOURCES. ' SO3O367O
IF(ITYPE-l) 6BO,840,840 S0303680
6BO WAKE = .FALSE. SO30369O
IF(VS) 690,690,700 SO303700
CHECK FOR DOWNWASH STACK HEIBHT ADJUSTMENT, VS - O. S03O371O
690 IF(ISW<25) .EQ. 2) HS = HS -3.*D S03O372O
IF EXIT VELOCITY, VS, EQUALS O THEN DHA = 0. S03O373O
DHAWAK = HS SO303740
IF(HS .LT. 2.5*HB .AND. HS .LT. HB-H.5*HW) WAKE = .TRUE. SO3O375O
BOTO 840 S0303760
7OO VSD = V9*D S03O377O
CHECK FOR DOWNWASH STACK HEIBUT ADJUSTMENT, VS > O. SO3O378O
IF(ISW(25) .EQ. 2 .AND. VS .LT. 1.5*UBARS) HS = HS + (VS*UBARI SO3O379O
1 -1.5)#(D+D) SO3O38OO
BAMJI = l./(.33333333+UBARS/VS) S03O381O
BAMJI = BAMJI*SAMJI SO3O3B2O
IF(DTH .LE. O.O) BOTO 71O SO3O383O
S = 9.B*DTH/TA SO303S4O
SI = l./S S0303850
S3 = SQRT(S) S0303B60
SSI = l./SS SO3O3B7O
710 IF(TS-TA) 730,730,720 SO3O3B80
IF SOURCE TEMPERATURE = O, SET EQUAL TO AMBIENT AIR TEMP. SO3O3B9O
72O IF(TS) 730,730,740 SO3O390O
730 FM = VSD*VSD#.25 SO303910
F = 0.0 S0303920
FZERO = .TRUE. SO3O3930
BOTO 77O SO3O394O
740 TOT = TA/TS SO3O395O
FM = TOT*VSD*VSD*.25 SO3O396O
F = 2.45*VSD*D*(i.-TOT) SO3O397O
IF(F .BT. 55.0) BOTO 750 SO30398O
C-10
-------
FC = .O727*VSD**1.3333333 S03O399O
GOTO 76O SO3O400O
750 FC = .0141*VSD**1.6666667 SO30401O
760 IF(F .BT. FC) GOTO 77O SO30402O
FZERO = .TRUE. S030403O
F = O.O SO30404O
770 IF(HB .LE. O.O) GOTO BOO SO3O4O5O
IFCDTH .GT. O.O) GOTO 7BO BO3O4O6O
C TEST FOR WAKE EFFECTS-CALCULATE XPLUME. SO3O4O7O
DHA = 3.*FM*(HB+HB)*GAMJI*UBARI*UBARI S03O4OBO
DHAWAK = DHA**.33333333 SO30409O
GOTO 79O SO3041OO
7BO DHA = 3.*FM*GAMJI*UBARI*SSI SO3O4MO
IF(i.57O796327*UBARS*SSI.GT.HB-i-HB) DHA = DHA*SIN(SS*(HB+HB)*UBARI)S0304115
DHAWAK = DHA**.33333333 SO30412O
DHA1 = 3.*VSD*UBARI SO30413O
IF(DHAWAK .GT. DHA1) DHAWAK = DHA1 SO3O414O
79O DHAWAK = HS + DHAWAK SO30415O
IF(DHAWAK.LT.2.5*HB .AND. DHAWAK.LT.HB+1.5*HW) WAKE = .TRUE. 5O3O416O
BOO IF
-------
NEXTR = 2
95O 1=1+1
IF (I .LE. NXPNTS) GOTO 96O
NEXTR = 1
GOTO 92O
96O IJ = IJ + 1
XR = GRIDX(I)
GOTO 990
97O 1=1+1
IF (I .GT. NXWYPT) GOTO 140O
YR = YDIS(I)
IF (.NOT. POLAR) GOTO 98O
IYR = YR
YRS = SINNUM(IYR)
YRC = COSNUMUYR)
98O IJ = NXPNTS*NYPNTS + I
XR = XDIS(I)
99O CONTINUE
IF (POLAR) GOTO 1OOO
XR1 = XR -XS
YR1 = YR - YS
GOTO 1O10
1000 XR1 = XR*YRS - XS
YR1 = XR*YRC - YS
C CHECK IF TERRAIN ELEVATION IS LOWER THAN STACK HEIGHT.
1010 IF(ISW(4).NE. l.OR.HS+ZS-GRIDZCIJ) . GT. O. O.OR. ITYPE. EQ. 2)
IF(LINE .EQ. O) WRITE(IO,9O11)
WRITE(IO,9012) NSO,XR,YR
STOP
C CALCULATE DOWNWIND DISTANCE, XBAR.
1020 XBAR => -(XR1*FVRSIN + YR1*FVRCOS)
IF (XBAR .LE. O.O) GOTO 92O
IF
-------
1O6O
107O
1O8O
*
1O9O
1O95
11OO
DHA = DMA**.33333333
BDTO 1O90
DHA = 3.*FM*SAMJI*UBARI*SSI
IF(XBAR .BE. XPLUME) GOTO 1OBO
XP1 = SS*XBAR*UBARI
DHA = DHA*SIN(XP1)
DHA = DHA**.33333333
DHAi = 3.*VSD*UBARI
IF(DHA . GT. DHAI) DHA = DHAI
EFFECTIVE PLUME HEIGHT.
H = HS + DHA
ADJUST H DUE TO TERRAIN
IF(ISW(4) .NE.l.QR.ISW(l).NE.1.OR.NVS.NE.O.OR.ITYPE.EQ.2)
Al = ZS - GRIDZUJ)
IF(A1 .GT. O.O) GOTO 1100
H = H + Al
CONTINUE
IF SO3O596O
1200 S0305970
XBARK + XY S03059BO
A2 =
SIGY
GOTO
XBARY
TH = .O17453293*(SC(ISTUM2)-SD(ISTUM2)*ALOG(XBARY))
SIGY = 465.11628*XBARY*TAN(TH)
1200 SIGYI = l./SIGY
IFdTYPE .ED. 2) GOTO 121O
Al = .5*(YBAR*SIGYI)**2
IF(A1 .GT. 5O.O) BOTO 92O
S03O599O
SO3O6OOO
S0306O1O
S0306020
S03O6O30
SO3O6O40
-------
121O IF(SeZDQN) GOTO 122O SO3O6O5O
CALL SIBMAZ(XBARZ,SIGZ,BBAR,ISTUM2,IXDIST,1,SASIBZ,SBSIBZ,DUMMY) SO3O6O6O
IFCSIQZ .BT. 5OOO. .AND. NVS .EQ. O) SIBZ = 5OOO. SO3O6O7O
122O SIBZI = l./SIBZ SO3O6O8O
C CALCULATE DECAY TERM. SO3O6O9O
XBARU = XBAR*UBARI S03O61OO
DECAYT = l.O SO3O611O
IF(DECAY(IHR) .BT. 0.0) DECAYT = EXP(-DECAY(IHR)*XBARU) SO30612O
C CHECK CONCENTRATION-DEPOSITION SWITCH. SO3O613O
IFUSW(l) .EQ. 2) SOTO 132O SO30614O
C CONCENTRATION EQUATION. S03O615O
C CHECK FOR PARTICIPATES WITH SETTLINB VELOCITIES. SO30616O
IF(NVS .ST. 0) BOTO 1260 SO30617O
IF(SISZ*HMI .LT. 1.6) BOTO 1240 SO3O618O
C CALCULATE "BOX-MODEL" CONCENTRATION SO30619O
IFdTYPE .EQ. 2) BOTO 123O SO3O620O
CHI = QTK*UBARI*SIBYI#HMI*EXP(-A1)*DECAYT#.39894228 S0306210
BDTO 139O SO3O622O
1230 A3 = .7071O678*SIEYI SO3O6230
A4 = (XOP+YBAR)*A3 SO3O624O
AS = -(XOP-YBAR)*A3 SO3O625O
A3 = ERFX(A4,A5) S030626O
CHI = QTK*XO*HMI*UBARI*A3*.5*DECAYT SO3O627O
BOTO 139O SO30628O
C CALCULATE VERTICAL TERM FDR ALL SOURCE TYPES W/0 PARTICLE S0306290
C BETTLINB VELOCITIES. SO3O63OO
124O V = O.O S03O631O
A2 = O.O SO306320
1250 VL * V SO306330
A2 * A2 + 2.0 SO30634O
HMA2 - A2*HM SO30635O
A3 ğ
-------
A3 = (Hri-HHM-H+XBARUV)*SIBZI S030671O
A5 = -.5*A3*A3 SO30672O
IF(A5 .GT. -5O.) SUM1 = EXP(A5) SO30673O
C THE FOLLOWINB LINE OF CODE ALTERED TO COMPUTE ESCAPE FRACTION
IF(GAMMA . LE. O.O) SOTO 12BO SO3O674O
A4 = (HM+HM+H-XBARUV)*SIGZI SO3O675O
A5 = -.5*A4*A4 SO3O6760
IF(A5 .BT. -50.) SUM1 = SUM1 + EXP(A5>*SAMMA SO30677O
CALL VERT(H,HM,XBARUV,SISZI,GAMMA,A2,SUM1) SO3O67BO
C THE FOLLQWINB LINES OF CODE ALTERED/ADDED TO COMPUTE ESCAPE FRACTION
12BO CALL ESCAPE(ZREF,ZO,TA,IST,UBAR,UD(IS,K),SOURCE(14,IS),ESCP) TRC OO3
V = V + .5*PHI*(SUM+SUM1)*ESCP SO3O679O
129O CONTINUE SO3O68OO
C CALCULATE CONCENTRATON FOR ALL SOURCE TYPES WITH VERTICAL TERM. S0306B10
13OO IFUTYPE .EQ. 2) BOTO 131O SO3O682O
CHI = QrK*UBARI*SIGYI*SISZI*EXP(-Al)*V*DECAYT*. 31830989 SO306B3O
BOTO 1390 SO306B4O
1310 A3 = .7O7iO678*SIEYI S0306B5O
A4 = (XOP+YBAR)*A3 S03O6B6O
A5 = - S0306930
WRITE(10,9013) NSO 50306940
STOP S0306950
C CALL SIBMAZ TO COMPUTE AVERABE EFFECTIVE DOWNWIND DISTANCE, BBAR. SO306960
1330 CALL SIBMAZ(XBARZ,SIBZ,BBAR,ISTUM2,IXDIST,2,SASIBZ,SBSIBZ,DUMMY) S030697O
V = O.O SO3O698O
DO 1370 K = 1,NVS S030699O
JP70 = K + 55 SO3O7OOO
BAMMA = SOURCE(JP7O,IS) S0307010
JP70 = K + 15 S0307020
PHI = SOURCE(JP7O,IS) S0307O3O
JP70 = K 4- 35 SO307040
XBARUV = XBARU*SOURCE(JP70,IS) SO307O50
A5 = (l.-BBAR)*XBARUV S03O7060
BAM1 = 1.0 S03O707O
BAM2 = BAMMA SO30708O
A2 = O.O S0307O9O
SUM = O.O SO3O71OO
1340 SUML = SUM SO3O711O
A2 = A2 + 2. SO3O712O
HMA2 = A2*HM SO3O713O
A3 = (HMA2-H+XBARUV)*SIGZI SO3O714O
A6 = O.O SO3O715O
A3 = -.5*A3*A3 S03O716O
IF (A3 .BT. -50.) A6 = EXP(A3)*BAM1*(BBAR*(HMA2-H)-A5) SO30717O
IF(SAMMA .BT. O.O) BOTO 135O S03O71BO
SUM = A6 S03O719O
BOTO 136O SO3O72OO
1350 A4 = (HMA2+H-XBARUV)*SIBZI S03O721O
A7 <= 0.0 50307220
A4 = -,5*A4*A4 SO307230
IF(A4 .ST. -50.) A7 = EXP(A4)*BAM2*(BBAR*(HMA2+H)+A5) S030724O
SUM = SUM + A6 + A7 SO3O7250
IF(ABS(SUM-SUML) .LT. l.E-B) BOTO 136O S03O726O
BAM1 = BAM2 S03O727O
BAM2 = BAM2*BAMMA SO3O72BO
BOTO 134O SO3O729O
136O A3 = (H-XBARUV)*SISZI SO3O73OO
A7 = -.5*A3*A3 5O3O731O
A3 = 0.0 50307320
C-15
-------
IF(A7 .BT. -5O.) A3
(BBAR*H
A5)*EXP(A7)
S0307330
THE FOLLOWING LINES OF CODE ALTERED/ ADDED TO COMPUTE ESCAPE FRACTION
CALL ESCAPE(ZREF,ZO,TA,IST,UBAR,UDUS,K> , SOURCE ( 14, IS) ,ESCP) TRC O04
V = V + ( 1. -GAMMA )#PHI*( A3 + SUM)*ESCP SO30734O
1370 CONTINUE
C FINISH DEPOSITION CALCULATIONS.
IFUTYPE .EQ. 2) BOTO 138O
CHI = QTK*SIGYI*SIGZI/XBAR*EXP(-A1)*DECAYT*V#. 15915494
SO TO 139O
138O CHI = QTK*XO*SIGZI/XBAR*DECAYT*V*ERFX((XOP+YBAR)*SIGYI*.7O71O67B
1 ,-(XOP-YBAR)*SIGYI*.7O71O67B)ğ. 39894228
C STORE CONCENTRATION OR DEPOSITION INTO CALC ARRAY. GO GET
C NEXT RECEPTOR.
1390 CALC(IJ) = CHI
UP = IJ + NPNTS
CALC (UP) = CALC (UP) + CHI
GOTO 92O
1400 CONTINUE
IF
-------
IB = 1 S03O7960
IFCNBROUP .LE. 0) GOTO 154O S0307970
153O IMS = NSOGRP(IG) SO3O79BO
ITO = NSUM + NS - 1 S0307990
C SO3O8OOO
C BEGIN LOOP OVER ALL TIME PERIODS FOR THIS HOUR. S030801O
C SO3OBO2O
1540 IAVB = O SO3O803O
DO 164O I = 1,8 SO3OBO4O
C FOR DAILY TABLES COMPUTE AVERAGES, WRITE TO TAPE & PRINT. S03OBO5O
IFdSWd+6) .NE. 1) GOTO 164O SO30BO6O
IAVG = IAVG + 1 SO30B070
IF(.NOT.IFLAG(I)) GOTO 164O SO3OBOSO
II = NPNT5*((IG-1>*NAVG + IAVG - 1) SO30B09O
IF
-------
DO 1660 I = 1,NPNTS
IP6 = I + NPNTS
166O CHIAN(I) = CHIANU) + CALC(IP6)
: END HOURLY LOOP.
167O CONTINUE
: CLEAR DAILY AVERAGES ARRAY BEFORE GOING TO NEXT DAY.
NPNTS2 = NAVG*NPNTS
IF(NGROUP .GT. O) NPNTS2 = NPNTS2*NGROUP
DO 16BO I = 1,NPNTS2
16BO CHIAV(I) =O.O
NTDAY = NTDAY + 1
169O CONTINUE
: END OF MET DATA.
NDAYS = NTDAY
NSUM = 1
IB *> 1
IFCNGROUP .LE. O) GOTO 171O
170O NS = NSOGRP(IG)
ITD = NSUM + NS - 1
: PRINT "N"-DAY TABLE
171O IF(ISM(15) -NE. 1) GOTO 173O
NHTOT = NTDAY*24
IF(ISW(19) .NE. 1) NHTOT = NDAYS*NHOURS
HTOT = 1./FLOAT(NHTOT)
IF(ISWU) .EQ. 2) HTOT = l.O
II = (IG-1)*NPNTS + 1
12 = II + NPNTS - 1
DO 172O I = 11,12
1720 CHIAN(I) = CHIAN(I)*HTOT
CALL DYOUT(BRIDX,GRIDY,XDIS,YDIS,CHIAN(I1),75,IDY,IHR,1,NSUM,ITQ,
1 16)
IF(ISW(5) .EQ. 1) WRITEUTAP) NHOURS, NTDAY, NGROUP, (CHIAN (I) ,
1 1=11,12)
C
C
BEGIN LOOP OVER TIME PERIODS.
173O IAVG = 0
DO 1750 I = 1,8
IF(ISW(I+6> .NE.
IAVB = IAVG + 1
: PRINT HIGHEST &
IF(ISW(17) .NE.
1) GOTO 1750
SECOND HIGHEST CONCENTRATION(DEPOSITION) TABLES.
1) GOTO 174O
IP6 = 4*NPNTS*((IG-1)*NAVG + IAVG - 1) + 1
CALL DYOUT(GRIDX,GRIDY,XDIS,YDIS,CHIMAX(IP6),KAVG(I),IDY,IHR,2,
1 NSUM,ITO,IG)
IP6 = IP6 + NPNTS + NPNTS
CALL DYOUT(GRIDX,GRIDY,XDIS,YDIS,CHIMAX(IP6),KAVG(I),IDY,IHR,3,
1 MSUM,ITO,IG)
: PRINT MAXIMUM 5O
1740 IF(ISM(IB) .NE. 1) GOTO 175O
IP6 = (IG-1)*NAVG + IAVG
CALL MAXOT(CHI5O(1,IP6),GRIDX,GRIDY,XDIS,YDIS,IPNT(1,IP6),
1 ICOUNT(IP6),KAVG(I),NSUM,ITO,IG)
175O CONTINUE
IG = IG + 1
IF(IG .GT. NGROUP) GOTO 1760
NSUH = NSUM + NS
1760
1770
9001
9002
1) GOTO 1770
GOTO 17OO
IF(ISW(5) .NE
ENDFILE ITAP
ENDFILE ITAP
RETURN
FORMATC1 ',121X,9HMET. DATA/122X,3HDAY,I4>
FORMAT(31X,69H* SOURCE-RECEPTOR COMBINATIONS LESS THAN 1OO METERS
1DR THREE BUILDING/34X,25HHEI6HTS IN DISTANCE. NO ,6A4,
2 16H IS CALCULATED *///46X,25H- - RECEPTOR LOCATION - -/SIX,
3 1HX,BX,10HY (METERS),1OX,8HDISTANCE/31X,6HSOURCE,1IX,
4 23HOR RANGE DR DIRECTION,9X,7HBETWEEN/31X,6HNUMBER,1IX,
5 21HCMETERS) (DEGREES),1IX,BH(METERS)/3OX,3O(2H- )/)
SO3O863O
S030B640
SO3O8650
S030S66O
SO3OS67O
SO3OB6aO
SO3OB69O
SO3OB7OO
SO3OB71O
S030B720
SO3OB73O
SO30874O
SO3OB75O
S030B760
S03O877O
SO3O87BO
S0308790
SO3OBBOO
SO30BB1O
SO30BS2O
S030BS30
S030S840
S0308B5O
SO3O886O
S030BB70
S03OBB8O
SO3O8B9O
SO3O89OO
SO3OB91O
50308920
SO3O893O
S03OS940
S0308950
SO30S96O
SO3OS97O
SO30898O
S030S990
SO3O9OOO
SO3O9O1O
S03O9O2O
S0309030
SO3O9O4O
S0309050
S0309060
SO3O9O7O
SO3O9O8O
SO3O909O
SO3O91OO
SO3O911O
SO3O912O
SO3O913O
SO3O914O
S03O915O
S0309160
SO3O917O
S0309180
SO3O919O
S03092OO
S0309210
SO3O922O
SO3O923O
S0309240
SO3O925O
50309260
SO3O927O
SO3O928O
SO3O929O
SO3O93OO
C-18
-------
9003 FORMAT(31X,I5,BX,2F13.1,7X,F1O.2> BO3O931O
9004 FORMAT(IB,5FB.O,I8,2F8.0> SO30932O
90O5 FORMAT(32X,4H*** ,15A4,4H ***//) SO3O9330
90O6 FORMAT(//6SX,1OHPOT. TEMP./29X,4HFLOW,7X,15HWIND MIXINB,13X, SO3O934O
1 BHBRADIENT,17X,16HWIND DECAY/2BX,16HVECTOR SPEED,5X, B03O935O
264HHEIBHT TEMP. (DEB. K STABILITY PROFILE COEFFICIENSO3O936O
3T/20X,92HHOUR S030939O
9OOB FORMAT(21X,12,Fll.i,F10.2,Fll.l,F9.1,F12.4,I9,F13.4,E15.6> SO3O94OO
9OO9 FORMAT(//47X,6HRANDOM/3BX,2(4HFLOW,6X),16H WIND MIXING,15X, SO30941O
1 19HIMPUT ADJUSTED/37X,2(6HVECTOR,4X),27H SPEED HEIGHT S030942O
2 TEMP. ,2(3X,9HSTABILITY)/29X,6HHOUR ,2UOH (DEGREES) ) ,3X, SO3O943O
3 3OH(MPS> (METERS) (DEG. K) ,2(BHCATEGORY,4X>/27X,4O(2H -)/)S030944O
9010 FORMAT(30X,12,Fll.1,F10.1,F1O.2,F11.1,F9.1,19,I12) SO3O945O
9011 FORMATC11) S0309460
9012 FORMAT(1OX,46H*** ERROR *** PHYSICAL STACK HEIGHT OF SOURCE,IS/ SO30947O
1 10X,52HIS LOWER THAN THE TERRAIN ELEVATION FOR THE RECEPTOR/1OX, S03094SO
1 12HLQCATED AT (,F 9.1,1H,,F 9.1,19H). RUN TERMINATED.) S03O949O
9013 FORMAT(1OX,25H***ERROR*** SOURCE NUMBER,16,41H HAS NO BRAVITATIONASO3O95OO
1L SETTLING CATEGORIES,/1OX,52HWITH WHICH TO CALCULATE DEPOSITION. SO30951O
2 RUN TERMINATED.) SO3O952O
END S0309530
C-19
-------
SUBROUTINE INCHK(CALC,CHIAV,CHIAN,GRIDX,BRIDY,XDIS,YDIS,BRIDZ, S020001O
1 CHIMAX,CHI50,IPNT,ICOUNT,SOURCE) SO20002O
C SUBROUTINE INCHK (VERSION BO339), PART OF ISCST.
C THIS ROUTINE READS THE REST OF THE INPUT VARIABLES AND PROVIDES S020003O
C DEFAULT VALUES IF REQUIRED. ALSO TABLES LISTING THE INPUT VARI- SO20004O
C ABLES ARE CONTROLLED BY THIS ROUTINE. SO2O005O
C . SO20OO6O
C THE FOLLOWING LINES OF CODE ALTERED TO RUN ON IBM-PC
CHARACTER*! ATHRUF
CHARACTER*4 TITLE,METER,SEASON,IBLANK,IQUN,ICHIUN,CONDEP
LOGICAL DONE SO200O7O
INTEGER WAKE,QFLG,QFLGS SO2O008O
COMMON /LOGIX/ ISW<4O),NSOURC,NXPNTS,NYPNTS,NXWYPT,NGROUP, S02OOO9O
1 NSOSRP(15O) ,IDSOR < 2OO),IPERD,NPNTS,NAVS,NHOURS,NDAYS,NTDAY,LINE, SO2OO1OO
2 IO,IN,TITLE(15),IQUN<3),ICHIUN(7),CONDEP<6>,LIMIT,MIMIT S020011O
COMMON /MET/ IDAY(366),ISTAB(24),AWS(24),TEMP(24),AFV(24), SO2OO12O
1 AFVR(24),HLH(24,2>,P(24),DTHDZ(24),DECAY(24),PDEF(6,6), SO2OO13O
2 DTHDEF(6,6),GAM1I,GAM2I,ZR,DDECAY,IMET,ITAP,TKTUCATS(5) SO20014O
C THE FOLLOWING LINE OF CODE ADDED TO COMPUTE ESCAPE FRACTION
COMMON/DEPO/UD(2OO,2O),ZREF,ZO TRC OO5
DIMENSION GRIDX(l),GRIDY(1),XDIS(1),YDIS(1),GRIDZ(1),SOURCE(215,1)SO20O15O
DIMENSION METER(2) ,SEASON(2,4),ATHRUF(6),UCTDEF(5) SO2O016O
EQUIVALENCE (ISW(23),QFLGS) SO2O017O
DATA ATHRUF / 'A','B','C','D',FE','F' / SO2OOIBO
DATA UCTDEF / 1.54,3.09,5.14,8.23,1O.B / S020019O
DATA METER /'(MET1,'ERS)'/ SO20O2OO
DATA SEASON /'WINT','ER ,'SPRI','NG 't'SUMM','ER ','AUTU', S020021O
1 'MN ' / S020022O
DATA IBLANK /' '/ SO2OO23O
C CHECK "ISW" AND SET DEFAULT VALUES. SO2OO24O
C DEFAULT TO CONCENTRATION ON RECTANGULAR GRID & DISCRETE POINTS. S02O025O
1O IF(ISWU) .LE. 0) ISW(l) = i S02O026O
IF(ISW(2) .LE. 0) ISW(2) = 1 S020027O
IF(ISW(3) .LE. 0) ISW(3) = 1 SO2002BO
C DEFAULT CARD MET PARAMETERS. SO20O29O
IF(NDAYS .LE. O) NDAYS = 1 S02O030O
C DEFAULT TO PRE-PROCESSED MET DATA WITH RURAL OPTION. S020O31O
IF(ISW(19) .LE. O) ISW(19) = 1 S020032O
IF(ISW(19) .EQ. 2) ISW(20) = 0 SO2OO33O
C DEFAULT TO PROGRAM"S WIND PROFILE EXPONENT AND VERTICAL POTENTIAL SO2O034O
C TEMPERATURE GRADIENT VALUES. SO2OO35O
IF(ISW(21) .LT. 1) ISW(21) = 1 S020O36O
IF(ISW(22) .LT. 1) ISW(22) = 1 S02OO37O
C DEFAULT TO FINAL PLUME RISE FOR ALL RECEPTORS. S02OO3BO
IF(ISW(24) .LT. 1) ISW(24) = 1 SO2OO39O
C DEFAULT TO NO STACK DOWNWASH ADJUSTMENT. SO2O04OO
IF(ISW(25) .LT. 1) ISW(25) = 1 SO20041O
C READ GRID THEN DISCRETE POINTS SO20O42O
IF(NXPNTS .EQ. O .OR. NYPNTS .EQ. O) GOTO 70 S02O043O
IF(ISW(2) .NE. 3) READ(IN,9O20) (GRIDX(I>,1=1,NXPNTS> SO2OO44O
IF(ISW(2) .LT. 3) READ(IN,9O2O) (GRIDY(I>,1=1,NYPNTS) SO2OO45O
IF(ISW(2) .NE. 3) GOTO 3O SO2OO46O
C GENERATE GRID, THEN READ DISCRETE POINTS. S020O47O
READ(IN,9O20) GRIDX(1),DX SO2OO4BO
12 = NXPNTS - 1 SO20O49O
DO 2O I = 1,12 SO2OO50O
II = I + 1 S020051O
20 GRIDX(ID = GRIDX(I) + DX SO2O052O
3O IF(ISW(2) .LT. 3) GOTO 5O SO20O53O
READ(IN,9O2O) GRIDY(1),DY S020054O
12 = NYPNTS - 1 S0200550
DO 40 I = 1,12 S0200560
II = I + 1 S0200570
C-21
-------
4O
50
60
70
80
9O
1OO
110
12O
121
125
C
C
130
14O
150
160
170
180
19O
200
t- DY
, ISW(2) .NE.
VALUES.
4) GOTO 7O
,OR. BRIDYd) .BT. 36O.O) BRIDYd)
. GT. 36O.O) YDIS(I)
36O.O
110
READ
NO OF SOURCE
GRIDYU1) = BRIDY(I)
IF(ISW(2) .NE. 2 .AND
SET DEFAULT DIRECTION
DO 60 I = 1,NYPNTS
IFCGRIDYd) .LE. O.O
IF(NXWYPT .EQ. O) BOTO 9O
READ(IN,9020) (XDIS(I),I=1,NXWYPT)
READ(IN,9O2O) (YDISd),I=1,NXWYPT)
IF(ISW<3) .NE. 2) GOTO 9O
SET DEFAULT DIRECTION VALUES.
DO BO I = 1,NXWYPT
IF(YDISd) .LE. O.O -OR. YDIS(I)
CHECK FOR TERRAIN HEIGHTS
IF(ISW<4) .NE. 1) GOTO 125
IF(NXPNTS .EQ. O .OR. NYPNTS .EQ. O) GOTO
READ TERRAIN FOR GRID AND DISCRETE REC"S;
DO 10O J = 1,NYPNTS
II = .NE
DO 17O J = 1,6
READ
-------
240
IFUCHIUNU) .NE.
CONTINUE
. EQ. 2) GOTO 250
= '(MIC*
= 'ROBR'
= 'AMS/'
= 'CUBI'
ME'
IBLANK) BOTO 26O
25O
IF(ISWd)
ICHIUN(l)
ICHIUN(2)
ICHIUN(3)
ICHIUN(4)
ICHIUN(5) = 'C
ICHIUNC6) = 'TER)'
ICHIUN(7) = '
GOTO 260
ICHIUN(l) = '(BRA'
ICHIUN(2) = 'MS/S'
ICHIUN(3) = 'DUAR*
ICHIUN(4) = 'E ME'
ICHIUN(S) = 'TER '
ICHIUN(6> = '
ICHIUN<7) = '
READ "DAY" ARRAY & MET IDENTIFICATION.
260 IF(ISM(19) .NE. 1) BOTO 27O
READ(IN,9022) ISS, ISY, IUS, IUY
NDAYS = 365
IF(MOD(ISY,4) .EQ. 0) NDAYS = 366
READ(IMET) ISSI,ISYI,IUSI,IUYI
IFUSS.EQ. ISSI.AND. ISY.EQ. ISYI.AND. IUS.EQ. IUSI.AND. IUY.EQ. IUYI)
1 GOTO 2SO
WRITE(10,9O25) ISS,ISSI,ISY,ISYI,IUS,IUSI,IUY,IUYI
STOP
FOR CARD MET DATA SET RURAL-URBAN SWITCH TO RURAL.
27O ISV-K20) = O
2BO IF(MSOURC .ST. O) GOTO 29O
WRITE(IO,9O26)
STOP
C*
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
S0201270
S02012BO
S02O1290
S02O1300
S02O1310
S02O132O
S020133O
S0201340
S02O1350
SO2O136O
S0201370
S02013BO
S020139O
S02O14OO
S0201410
S0201420
S0201430
SO201440
S0201450
SO2O146O
S0201470
S02O14BO
SO2O1490
S0201500
S0201510
SO2O1520
S02O1530
S02O1540
S0201550
SO2O1560
S0201570
SO2015BO
SO201590
50201600
SO201610
S0201620
S0201630
290
300
READ SOURCE DATA.
MOST VARIABLES ARE READ DIRECTLY INTO THE "SOURCE" ARRAY WHICH
HAS 215 STORAGE LOCATIONS ALLOCATED PER SOURCE. STORAGE LOCATION S0201640
1 CONTAINS WAKE, QLFG, NVS & ITYPE PACKED INTO THE FIRST LOCATION.S0201650
STORAGE LOCATIONS 2-13 CONTAIN/. NSO, Q, X, Y, ZS, HS, TS OR SO2O1660
SIGZO, VS OR SIGYO OR XO, D, HB, BUILDING LENBTH, AND BUILDING S02O1670
WIDTH, RESPECTIVELY. STORAGE LOCATIONS 16-35 CONTAIN PHI, 36-55 S02O16BO
CONTAIN SETTLINS VELOCITIES AND 56-75 CONTAIN GAMMA. STORAGE S020169O
LOCATIONS 76-81 CONTAIN STABILITY-DEPENDENT LATERAL VIRTUAL S02O1700
DISTANCES AND LOCATIONS 82-117 CONTAIN STABILITY AND XBAR- SO2O171O
DEPENDENT VERTICAL "IRTUAL DISTANCES BOTH OF WHICH ARE COMPUTED SO2O1720
IN SUBROUTINE MODEL. STORABE LOCATIONS 120-215 CONTAIN G! SO2O173O
ADJUSTMENT FACTORS AS A FUNCTION OF EITHER TIME OF DAY-SEASONAL ORSO2O174O
STABILITY-WIND SPEED VARIATIONS. STORAGE LOCATIONS 14, 15, 118, &SO2O1750
119 ARE CURRENTLY NOT BEINS USED. SO2O1760
SO2O1770
II = 1 SO2O17BO
IFUI .GT. NSOURC) BOTO 32O S02O179O
THE FOLLOWING LINES OF CODE ALTERED TO COMPUTE ESCAPE FRACTION
READdN,9O27) NSO,ITYPE,WAKE,NVS,QFLB,(SOURCE(I,11) ,1=3,13),
1 PITDEP(II)
l^**#*****#*****#*#*#ğ***#****ğ#*###*ğ*^^^^^(.^^****^^^^^^^^^^^^^^^(.^^^^^|.^^*^^^^^^*^^^(.^
IF(NVS .LE. O) BOTO 31O
INPUT VARIABLES RELATED TO PARTICULATE SOURCES.
READ(IN,9020) (SOURCE(15+1,II),1=1,NVS)
READ(IN,9020) (SOURCE(35+1,II),1=1,NVS)
READ(IN,9020) (SOURCE(55+1,II),1=1,NVS)
READ(IN,9020) (UD(11,I),1 = 1,NVS)
31O CONTINUE
PACK SOURCE VARIABLES WAKE, QFLB, NVS & ITYPE INTO FIRST
ALSO STORE SOURCE NUMBER.
SOURCE(1,II) = ITYPE + NVS#16 + QFLB*512 + WAKE#B192
SO2O1BOO
TRC 006
Ğ********
S0201B10
S0201B20
S0201B30
502O1S40
SO2O1B5O
SO2O1S6O
LOCATION.5020187O
S0201BBO
SO2O1B9O
C-23
-------
SOURCE(2,II) = NSO S02019OO
II = II + 1 SO2O191O
GOTO 30O S02O192O
ENTER SOURCE EMISSION RATE SCALARS. 30201930
32O II = 1 SO2O194O
IF(QFLGS .LT. 1 .OR. QFLBS .ST. 5) GOTO 330 SO2O195O
DONE = .TRUE. S02O196O
DFLB = QFLBS SO2O197O
BOTO 350 SO2019BO
33O DONE = .FALSE. S02O199O
34O IF(II .ST. NSOURC) GOTO 43O SO202OOO
ITYPE = SOURCE(1,11) SO2O2O1O
QFLB = ITYPE/512 - (ITYPE/8192)*16 SO20202O
IF(DFLB .LT. 1 .OR. QFLB .BT. 5) GOTO 420 SO2O2O3O
35O J = 1 SO202O4O
1=4 SO2O2O5O
GOTO (400,360,370,300,390), QFLB SO2O2O6O
360 I = 12 S0202070
BOTO 400 SO2O2OSO
370 I = 24 S0202090
BOTO 400 S0202100
3BO J = 6 50202110
1=6 SO2O212O
BOTO 40O S02O213O
390 J = 4 S0202140
I = 24 SO20215O
4OO DO 410 II = 1,J S020216O
IFR = ,I=1,NGROUP) S02O245O
13 = 0 S0202460
DO 46O I = l.NBROUP SO20247O
460 13 = 13 + NSOGRP(I) S02024BO
WRITE(IO,9O58) (IDSOR(I),1=1,13) S020249O
PRINT UPPER BOUND OF FIRST 5 WIND SPEED CATEGORIES. SO20250O
470 LINE = LINE + 6 SO20251O
WRITE(IO,9OO1) (UCATS(I),1=1,5) SO20252O
IF(ISW(19).EQ.2.AND.ISW(6).EQ.2) GOTO 530 SO20253O
IF(ISW(21) .EQ. 3) GOTO 5OO BO2O254O
PRINT WIND PROFILE EXPONENTS. SO2O255O
LINE = LINE +12 . . S02O256O
IF(LINE .LT. 57) GOTO 480 SO2O257O
C-24
-------
48O
49O
50O
510
52O
530
54O
LINE = 15
WRITE (11,11 = 1,6)
DO 490 I = 1,6
WRITE(IO,9017) ATHRUF(I),(PDEF(J,I),J=1,6)
IF(ISW(22) .EQ. 3) BOTQ 53O
PRINT VERTICAL POTENTIAL TEMPERATURE GRADIENTS.
LINE = LINE + 12
IF(LINE .LT. 57) GOTO 51O
LINE = 15
WRITE(IO,9O29) TITLE
WRITE(IO,9O6O)
WRITE(IO,9016) (11,11=1,6)
DO 52O I = 1,6
WRITE(IO,9017> ATHRUF(I) , (DTHDEF(J,I) ,J=1,6>
PRINT RECEPTOR INFO.
NYPNTS
550
56O
570
580
59O
6OO
WRITE(IO,9038>
WRITE(IO,9O39)
WRITE(IO,9O41)
WRITE(IO,9O42)
IF(NXPNTS .EQ. 0 .O.I. NYPNTS .EQ. O> GOTO 550
LINE = LINE + 20
IF(LINE .LT. 57) GOTO 54O
LINE = 6
WRITE(10,9029) TITLE
IF(ISW(2) .EQ. 1 .OR. ISW(2) .EQ. 3)
IF(ISW(2) .EQ. 2 .OR. ISW(2) .EQ. 4)
WRITE(IO,9040) (GRIDX(I),I=1,NXPNTS)
IF(ISW(2) .EQ. 1 .OR. ISW(2) .EQ. 3)
IF(ISW(2) .EQ. 2 .OR. ISW(2) .EQ. 4)
WRITE(IO,9O4O) (GRIDY(I),1=1,NYPNTS)
IF(NXWYPT .EQ. 0) GOTO 57O
LINE = LINE + 5 + NXWYPT/5
IF(LINE .LT. 57) GOTO 560
LINE = 6
WRITE(IO,9029) TITLE
IF(ISW(3) .EQ. 1) WRITE(IO,9O43)
IF(ISW(3) .EQ. 2) WRITE(IO,9044)
WRITE(I0,9045) (XDIS(I),YDIS(I),I=1,NXWYPT)
PRINT TERRAIN HEIGHTS.
IF(ISW(4) .NE. 1) GOTO 5BO
CONDEP(3) = 'HGT '
CALL DYOUT(GRIDX,GRIDY,XDIS,YDIS,GRIDZ,99,IDY,IHR,1,O,O,O)
PRINT OUT SOURCE INFO.
CONTINUE
LINE = 100
13 = O
DO 60O I = 1,NSOURC
IF(LINE .LE. 56) GOTO 59O
WRITE(10,9029) (TITLE(J),J=1,15)
WRITE(IO,9046) ((IQUN(J),J=1,3),12=1,2),(METER(l),METER(2),J<
LINE = IB
CONTINUE
ITYPE = SOURCE(1,1)
GET WAKE OPTION, SOURCE NO., NVS & TYPE FROM FIRST WORD.
NSO = SOURCE(2,I)
WAKE = ITYPE/8192
QFLG = ITYPE/512 - (ITYPE/B192)*16
NVS = ITYPE/16 - (ITYPE/512)*32
ITYPE = ITYPE - (ITYPE/16)*16
IF(NVS .GT. 0) 13 = 1
WRITE(IO,9047) NSO,ITYPE,WAKE,NVS,(SOURCE(J,I),J=3,13)
LINE = LINE + 1
CONTINUE
IF(13 .NE. 1) GOTO 63O
PRINT OUT PARTICLE CATEGORY INFORMATION.
LINE = 1OO
DO 620 I = 1,NSOURC
IF(LINE .LT. 43) GOTO 610
WRITE(10,9029) (TITLE(J),J=1,15)
WRITEUO.9049)
SO2025BO
S02O259O
S02026OO
S0202610
S0202620
S0202630
S02O264O
SO2O265O
S020266O
SO2O267O
50202680
SO20269O
SO2027OO
SO20271O
S02O272O
S0202730
S020274O
S0202750
S0202760
SO2O277O
S02027BO
S020279O
S0202BOO
SO2O2S1O
S0202B20
SO202B30
S0202B4O
SO2O2B5O
SO2O2S6O
SO202B7O
S0202B8O
SO2O2B9O
SO2029OO
50202910
S020292O
SO2O293O
S0202940
S0202950
S020296O
S0202970
S02029BO
S0202990
SO203OOO
S020301O
S02O3O2O
SO2O303O
S0203O4O
=1,1O)SO2O3O5O
SO2O3O6O
S0203O7O
SO203OBO
S0203O9O
S0203100
SO2O311O
S020312O
S02O313O
SO2O314O
S020315O
S0203160
9O20317O
S02031BO
S0203190
SO20320O
S0203210
S020322O
SO20323O
S020324O
SO2O325O
C-25
-------
LINE = 1O SO2O326O
610 CONTINUE SO2O327O
ITYPE = SOURCE(1,1) SO2O328O
NSO = SOURCE<2,I) SO2O3290
NVS = ITYPE/16 - (ITYPE/512)*32 SO20330O
IF(NVS .LE. O) BQTD 620 SO2O331O
WRITE(I0,9050) NSO S02O3320
12 = 15 * NVS SO2O3330
WRITE(IO,9051) ,J=16,I2) SO2O334O
12 = 35 + NVS SO2O335O
WRITE,J=56,12) SO2O33BO
LINE = LINE + 14 SO2O339O
620 CONTINUE SO2O340O
PRINT SOURCE EMISSION RATE SCALARS. SO2O341O
63O 1=1 SO2O342O
IF(QFLQS .LT. 1 .OR. QFLGS .BT. 5) GOTO 64O SO2O343O
DONE = .TRUE. S02O3440
DFLG = QFLGS S0203450
LINE = 100 SO20346O
GOTO 670 S0203470
64O DONE = .FALSE. SO2O34BO
J = 1 SO20349O
650 IF(J .GT. 5) GOTO B2O S02035OO
LINE = 10O SO20351O
1=1 . SO2O352O
660 IF(I .GT. NSOURC) GOTO BIO SO20353O
ITYPE = SOURCE(1,1) SO2O354O
QFLG = ITYPE/512 - (ITYPE/B192)*16 SO2O355O
IF(DFLG .NE. J) BOTD BOO SO20356O
NSQ = SOURCE(2,I) S0203570
670 GOTO (680,700,720,740,770), QFLG SO2O35BO
680 IF(LINE .LT. 54) GOTO 69O SO203590
WRITE(I0,9029) TITLE SO2036OO
WRITE(10,9002) S02O361O
IF(DONE) WRITE(ID,9003) SO20362O
WRITE(IO,9004> {(SEASON(I1,12),I1=1,2),12=1,4) SO20363O
LINE = 14 SO2O364O
690 IF(.NOT.DONE) WRITE(10,90O5) NSO SO2O365O
WRITE(IO,9OO6) (SOURCE(I1,I),I1=120,123) SO20366O
IF(DONE) GOTO 82O SO203670
LINE = LINE + 3 SO2O36BO
GOTO BOO SO2O369O
7OO IF(LINE .LT. 54) GOTO 71O SO2O37OO
WRITE(ID,9O29) TITLE SO2O371O
WRITE(IO,9OO7) SO2O372O
IF(DONE) WRITE(I0,9003) SO2O373O
WRITE(IO,900S) S0203740
WRITE(IO,9013) S0203750
LINE = 14 SO20376O
71O IF(.NOT. DONE) WRITE(IO,9O09) NSQ SO20377O
WRITE(IO,9O1O) (SOURCE(II,I),11=120,131) S02037BO
IF(DONE) GOTO B2O SO20379O
LINE = LINE + 3 S0203SOO
GOTO 8OO SO203B1O
720 IF(LINE .LT. 50) GOTO 73O S02O382O
WRITE(IO,9O29) TITLE SO203B3O
WRITE(IO,9011) SO2O384O
IF(DONE) WRITE(IO,9O03) SO2O3B5O
WRITE(10,9012) SO203B6O
WRITE(10,9013) SO203B7O
LINE = 14 S0203BBO
730 IF(.NOT.DONE) WRITE(IO,9OO9) NSO SO203B9O
WRITE(IO,9O14) (II,SOURCE(119-H1,I),11=1,24) SO2O39OO
IF(DONE) GOTO 82O SO20391O
LINE = LINE + 7 SO20392O
GOTO 800 S0203930
C-26
-------
740 IF(LINE .LT. 49) GOTO 750 SO20394O
WRITE(I0,9029) TITLE SO203950
WRITE*6 + 12O SO204O3O
ITO = IFR + 5 SO204O4O
760 WRITEdO,9017) ATHRUF (11) , (SOURCE (12,1) , I2=IFR, ITO) S0204O5O
IF(DONE) GOTO 82O SO204O6O
LINE = LINE + 8 SO204O70
GOTO BOO SO204080
770 IF(LINE .LT. 37) GOTO 7BO SO2O4090
WRITE(IO,9O29) TITLE SO2041OO
WRITEdO,9O1B) S020411O
IF(DONE) WRITE(IO,9OO3) SO20412O
WRITEdO,9012) S0204130
WRITEdO,9013) S0204140
LINE = 14 S0204150
7BO IF(.NOT.DONE) WRITE(10,9OO9) NSO SO20416O
DO 79O II = 1,4 S02O417O
IFR = dl-l)*24 + 119 SO2041BO
WRITEdO,9O19) SEASON(1,I1) ,SEASON (2,11) S02O419O
790 WRITE(IO,9O14) (12,SOURCE(I2-HFR,I),12=1,24) S02042OO
IF(DONE) GOTO B2O SO204210
LINE = LINE + 22 S020422O
BOO 1=1+1 S02O4230
GOTO 66O SO2O424O
BIO J = J + 1 S02O425O
GOTO 65O SO20426O
C STORE RECIPROCAL SQUARED OF BETA1, BETA2 AS GAM1I, GAM2I AND STORES020427O
C RECIPROCAL OF ZR. S02O42BO
B20 CONTINUE S02O429O
GAM1I = 1./(BETA1*BETA1) SO2O43OO
GAM2I = l./(BETA2*BETA2) S020431O
ZR = l./ZR SO2O432O
C COMPUTE EFFECTIVE BUILDING WIDTH FOR ALL SOURCES & STORE IN SO2O433O
C LOCATION 12 OF "SOURCE" ARRAY. BUILDING LENGTH & WIDTH WILL NO S020434O
C LONGER BE NEEDED. ALSO, RELOCATE AREA SOURCE COORDINATES FROM S020435O
C THE SOUTHWEST CORNER TO THE CENTER OF THE AREA SOURCE. S02O436O
DO 83O I = 1,NSOURC S02O437O
C 2/SQRT(3.14159265) = 1.12B3792 S02043BO
SOURCE(12,I) = 1.12B3792*SQRT(SOURCE(12,1)*SOURCE(13,1)) S020439O
ITYPE = SOURCE(1,1) SO2O44OO
IF(ITYPE-(ITYPE/16>*16 .NE. 2) GOTO B30 SO2O441O
Al = .5*SOURCE(9,I) SO2O442O
SOURCE(4,I) = SOURCE(4,I) + Al . SO2O443O
SOURCE(5,I) = SOURCE(5,I) + Al SO2O444O
B30 CONTINUE S020445O
C SET HEADING. SO2O446O
IFdSW(l) .ED. 1) GOTO B4O SO2O447O
CONDEP(l) = ' TO* SO2O44BO
CONDEP(2) = 'TAL ' SO2O449O
CONDEP(3) = 'DEPO' SO2O45OO
CONDEP(4) = 'SITI' SO2O451O
CONDEP(5) = 'ON ' S020452O
CONDEP(6) = ' ' S020453O
GOTO B50 S020454O
B4O CONDEP(l) = 'AVER' S02O4550
CONDEP(2) = 'AGE ' SO20456O
CONDEP(3) = 'CONC* S02O4570
CONDEPU) = 'ENTR' SO2O45BO
CONDEP(5) = 'ATID' S0204590
CONDEP(6) = 'N ' SO2046OO
B50 CONTINUE SO20461O
RETURN SO2O462O
C-27
-------
9OO1 FORMAT(/34X.64H*** UPPER BOUND OF FIRST THROUGH FIFTH WIND SPEED CS02O4630
1ATEBORIES ***/6OX,12H(METERS/SEC)//46X,5(F7.2,1H,)) S02O464O
9O02 FORMAT(39X,54H* SOURCE EMISSION RATE SCALARS WHICH VARY SEASQNALLYSO2O465O
1 *//> S02O466O
90O3 FORMAT(56X,19H* FOR ALL SOURCES *//) S02O467O
9OO4 FORMAT(4OX,4(2A4,7X)/20X,4O(2H- )/) S02046SO
9OO5 FORMAT(/20X,12HSOURCE NO. = ,I6> S020469O
9O06 FORMAT(3SX,4(E1O.5,5X)) S02O47OO
9007 FORMAT(41X.51H* SOURCE EMISSION RATE SCALARS WHICH VARY MONTHLY * S02O4710
1 //) S0204720
9008 FORMAT(7X,51HJANUARY FEBRUARY MARCH APRIL MAY , S02O4730
1 5BHJUNE JULY AUGUST SEPTEMBER OCTOBER NOVEMBER , S02O4740
2 8HDECEMBER/) S02O475O
90O9 FORMAT SO2O4B7O
9018 FORMAT(32X,68H* SOURCE EMISSION RATE SCALARS WHICH VARY SEASONALLYS02O4BBO
1 AND DIURNALLY *//) S02O489O
9019 FORMAT(59X,9HSEASON = ,2A4) SO2O49OO
9020 FORMAT(BF10.0) SO2O491O
9021 FDRMAT(EB.O,3F8.O,3A4,7A4,212> S0204920
9O22 FORMAT(BOI1) SO20493O
9023 FORMAT(2014) SO20494O
9024 FORMAT(1316) SO2O495O
9025 FORMAT<'1',1OX,63H*ğ* ERROR *** MET DATA REQUESTED DOES NOT MATCHSO2O496O
1 MET DATA READ./10X ,28H"REQUESTED/READ" VALUES ARE7./1OX, SO2O497O
2 21HSURFACE STATION NO. =,16,1H/,I6,23H YEAR OF SURFACE DATA =,I6,S02049BO
3 1H/,I6/10X,23HUPPER AIR STATION NO. =>, 16,1H/, 16, S0204990
4 25H YEAR OF UPPER AIR DATA =,16,1H/,I6/10X,15HRUN TERMINATED.) SO2O5OOO
9026 FORMAT('i',10X,73H**ğ ERROR *** NUMBER OF SOURCES TO BE READ EQUASO2O5O1O
1LS ZERO. RUN TERMINATED.) SO2O5O2O
C THE FOLLOWINS LINE OF CODE ALTERED TO COMPUTE ESCAPE FRACTION
9027 FORMAT(I5,2I1,I2,I1,EB.O,2F7.0,9F6.0) SO205O30
9O2B FORMAT(9X,II,5F1O.O) SO2O5O4O
9O29 FORMAT('1"//32X,4H*** ,15A4,4H ***//) SO2O5O5O
9O3O FORMAT(1BX,4OHCALCULATE (CONCENTRATION*!,DEPOSITION=2),29X, SO2O5O6O
1 BHISW(l) =,I4/1BX,55HRECEPTOR GRID SYSTEM (RECTANGULAR=1 OR 3, POSO2O5O7O
2LAR=2 OR 4),14X,8HISW(2) =,I4/ SO2O5OBO
3 18X,48HDISCRETE RECEPTOR SYSTEM (RECTANGULAR=1,POLAR=2),21X, SO205O90
4 BHISW<3) =,I4/,1BX,4OHTERRAIN ELEVATIONS ARE READ (YES=1,NO=0), SO2O51OO
5 29X,8HISW(4) =,I4/,1BX, SO2O511O
6 45HCALCULATIONS ARE WRITTEN TO TAPE (YES=1,NO=O),24X,BHISW(5) =, SO2O512O
7 I4/1BX, S0205130
8 4SHLIST ALL INPUT DATA (NO=O,YES=1,MET DATA ALSO=2),21X, S020514O
9 BHISW(6) =,I4//18X,39HCOMPUTE AVERAGE CONCENTRATION (OR TOTAL, SO20515O
O 12H DEPOSITION)/18X,32HWITH THE FOLLOWING TIME PERIODS7./2OX, SO2O516O
1 19HHOURLY (YES=1,NO=0),4BX,8HISW(7) =,I4/2OX, SO20517O
2 19H2-HOUR (YES=1,NO=O),4BX,8HISW(8) =,I4/2OX, SO2O51BO
3 19H3-HOUR
-------
3 25HDAILY TABLES (YES=1,NO=O),41X,9HISW(16) =,I4/20X, S02O52BO
4 44HHIBHEST & SECOND HIGHEST TABLES (YES=1,NO=O),22X,9HISW(17) =, S02O5290
5 I4/2OX,3OHMAXIMUM 5O TABLES (YES=1,NO=O),36X,9HISW(IS) =,I4/18X, SO2053OO
6 57HMETEQRDLOGICAL DATA INPUT METHOD (PRE-PROCESSED=1,CARD=2),1IX,30205310
7 9HISWU9) =,I4/1BX,5SHRURAL-URBAN OPTION (RURAL=O,URBAN MODE 1=1,30205320
BURBAN MODE 2=2),1OX,9HISW(20) =,I4/1BX,57HWIND PROFILE EXPONENT VASO20533O
9LUES (DEFAULTS=1,USER ENTERS=2,3),1IX,9HISW(21) =,I4/1BX, SO20534O
O 64HVERTICAL POT. TEMP. GRADIENT VALUES (DEFAULTS=1,USER ENTERS=2,S02O535O
13),4X,9HISW(22) =,I4/1BX,49HSCALE EMISSION RATES FOR ALL SOURCES (SO2O5360
2NO=O,YES>O),19X,9HISW(23) =,I4/1BX,53HPROGRAM CALCULATES FINAL PLUSO20537O
3ME RISE ONLY (YES=1,NO=2),15X,9HISW(24) =,I4/1BX, SO2O53BO
4 59HPROGRAM ADJUSTS ALL STACK HEIGHTS FOR DOWNWASH (YES=2,NO=1), 50205390
5 9X,9HISW(25) =,14) S02O54OO
9032 FORMAT(/1BX,23HNUMBER OF INPUT SOURCES,46X,BHN5OURC =,I4/1BX, S020541O
1 40HNUMBER OF SOURCE GROUPS (=O,ALL SOURCES),29X,BHNGROUP =,I4/1BXS02O542O
2,53HTIME PERIOD INTERVAL TO BE PRINTED (=O,ALL INTERVALS),17X, SO20543O
3 7HIPERD =TI4/18X,31HNUMBER OF X (RANGE) GRID VALUES.38X,8HNXPNTS S020544O
4=,I4/18X,31HNUMBER OF Y (THETA) GRID VALUES,38X,BHNYPNTS =,I4/1SX,S02O545O
5 28HNUMBER OF DISCRETE RECEPTORST41X,BHNXWYPT =,14) SO20546O
9033 FORMAT(18X,46HNUMBER OF HOURS PER DAY IN METEOROLOGICAL DATA,23X, S02O547O
1 BHNHOURS =,I4/iBX,37HNUMBER OF DAYS OF METEOROLOGICAL DATA,32X, SO2054BO
2 8H NDAYS =,14) S020549O
9O34 FORMAT(1BX,44HSOURCE EMISSION RATE UNITS CONVERSION FACTOR,27X, S02O550O
1 6H TK =,E1O.5/18X,47HENTRAINMENT COEFFICIENT FOR UNSTABLE ATMOSPS020551O
2HERE,22X,BH BETA1 =,F5.3/1BX,45HENTRAINMENT COEFFICIENT FOR STABLESO20552O
3 ATMOSPHERE,24X,8H BETA2 =,F5.3/,18X, SO2O553O
4 52HHEIGHT ABOVE GROUND AT WHICH WIND SPEED WAS MEASURED,18X, S0205540
5 7H 2R =,F7.2,BH METERS/1BX, S02O555O
6 42HLOGICAL UNIT NUMBER OF METEOROLOGICAL DATA,29X,6HIMET =,14) SO2O556O
9035 FORMAT(18X,52HDECAY COEFFICIENT FOR PHYSICAL OR CHEMICAL DEPLETIONSO20557O
1 ,18X,7HDECAY =,E12.6/1BX,19HSURFACE STATION NO., SO2O55BO
3 53X,5HISS =,I6/1BX,20HYEAR OF SURFACE DATA,52X,5HISY =,I3/18X, SO20559O
4 21HUPPER AIR STATION NO.,51X,5HIUS =,I6/1BX, S02O56OO
5 22HYEAR OF UPPER AIR DATA,SOX,5HIUY =,13) S02O561O
9036 FORMAT(1BX,39HLOGICAL UNIT OF CALCULATION "SAVE" TAPE,3OX, SO2O562O
1 8H ITAP =,I4> 30205630
9037 FORMAT(44X,43H*ğ* METEOROLOGICAL DAYS TO BE PROCESSED ***/ S020564O
1 63X,6H(IF=1)//B(11X,5(1012,2X)/>) SO2O565O
9O3B FORMAT(//42X,4BH*ğ* X-COORDINATES OF RECTANGULAR GRID SYSTEM ***/ S02O566O
1 62X,8H(METERS)/) SO20567O
9039 FORMAT(//47X,35H*#* RANGES OF POLAR GRID SYSTEM ***/62X, SO2O56BO
1 BH(METERS)/) S020569O
9O40 FORMAT(1OO<5X,1O(F1O.1,1H,)/)) SO2O57OO
9O41 FORMAT(//42X,48Hğ** Y-COORDINATES OF RECTANGULAR GRID SYSTEM *** SO2O571O
1 /62X,8H(METERS)/) SO2O572O
9O42 FORMAT(//45X,42H**ğ RADIAL ANGLES OF POLAR GRID SYSTEM ***/ SO20573O
1 /62X.9H(DEGREES)/) SO2O574O
9O43 FORMAT; (M/SEC)5,12X,3(5HBLDG.,4X)/1OX, SO2O5B2O
3 20HY A NUMBER TYPE=2,25X,4HBASE,12X,53HVERT.DIM HORZ.DIM DIAMSO205B3O
4ETER HEIGHT LENGTH WIDTH/3X,19HSOURCE P K PART. ,3A4.5X, SO205B4O
5 1HX,BX,43HY ELEV. HEIGHT TYPE=1 TYPE=1,2 ,4<6HTYPE=O,S020585O
6 3X)/3X,31HNUMBER E E CATS. *PER METER**2,2(5(1X,2A4),IX)/ SO20586O
7 63(2H -)/) SO205B7O
9047 FQRMAT(IB,I3,I2,I5,3X,E11.5,2F10.1,FB.1,2F9.2,IX,5F9.2) S02O58BO
904B FORMAT(IB,I3,I2,I5,3X,2A4,A3,1X,2F9.1,3F9.2,IX,5F9.2) SO205B9O
9O49 FORMAT(SOX,31H*#* SOURCE PARTICULATE DATA **#//) SO2O590O
9050 FORMAT(/10X,19H*** SOURCE NUMBER =,I6,4H ***) SO2O591O
9O51 FORMAT(/10X,15HMASS FRACTION =/2(1OX,1O(F7.5,1H,)/)) SO20592O
9052 FORMAT(/1OX,31HSETTLING VELOCITY(METERS/SEC) =/2(1OX,10(F7.4,1H,) S02O593O
1 /)) SO2O594O
9053 FORMAT S0205960
C-29
-------
9054 FORMAT(34X,27H* SEASONAL SOURCE STRENGTHS,3A4, SO2O597O
1 26HFQR EACH HOUR OF THE DAY *//20H *** SOURCE NUMBER =,I5,4H ***)SO2O598O
9055 FORMAT
-------
ESCAPE FRACTION SUBROUTINE
ALTERNATIVE 1
CONSTANT-K, LINEAR MODEL
C-31
-------
-------
SUBROLJTINE ESCAPE (ZREF,ZO,TA,IS,U,UD,H,EBCP)
SUBROUTINE ESCAPE, ALTERNATIVE-! CONSTANT-K LINEAR MODEL
A5=ALOG(ZREF/ZO>
A6=1./.123/ZREF
A1=UD*A5*H/U*A6
ESCP=1./<1.+A1>
RETURN
END
C-33
-------
SUBRQUTINE ESCAPE < ZREF,2O,TA,IS,U,UD,H,ESCP)
DIMENSION DTDZ(6)
DATA DTDZ/-.O1,-.OO7,-.005,0.,.02,.O35/
A5=ALOS(ZREF/ZO>
A6=l./.123/ZREF
B=9.B1*ZREF*ZREF*DTDZ(IS)/TA/U/U
A1=AI_OS(ZREF/ZO>
IFCB.LT.O.) BOTD 1
XH=0.05
X=O.15
FH=XH/ ( (A1+.33333)/I. 33333) **2-B
1O IF #*. 25
ZETA=U.O-15*X)**.25
ZETAO=(1.0-15*X*ZO/ZREF)**.25
ARB1=ALOB<(ZETA-1.0)*(ZETA-H.O)/((ZETA+1.O)*(ZETAO-1.O)))
ARG2=2.O*(ATAN < ZETA)-ATAN(ZETAO))
F=X/((Al-PS)/PH)**2-B
IF(ABS
-------
ESCAPE FRACTION SUBROUTINE
ALTERNATIVE 3
VARIABLE-K, LINEAR MODEL
C-37
-------
SUBROUTINE ESCAPE(ZREF,ZO,TA,IS,U,UD,H,ESCP)
A5=ALOB(ZREF/ZO)
A6=l./.123/ZREF
A1=UD*A5*ALOG(H/ZO)/U/.123
ESCP=1./(1.+A1>
RETURN
END
C-39
-------
ESCAPE FRACTION SUBROUTINE
ALTERNATIVE 4
VARIABLE-K, DETAILED MODEL
C-41
-------
SUBROUTI ME ESCAPE < 2REF, ZO, TA, IS, U, UD, H, ESCP)
DIMENSION DTDZ(6)
DATA DTDZ/-.O1,-.OO7,-.OO5,O.,.02,.O35/
A5=ALOB
A6=l./.123/ZREF
DZ=H/1O.
TEDDY=0.
DO 3 1=1,10
Z=DZ*I-DZ/2
CALL KCAL(ZREF,Z,ZO,TA,DTDZ
RETURN
END
SUBROUTINE KCAL < ZREF,Z,ZO,TA,DTDZ,U,UD,H,EDDY)
B=9.Si*ZREF*ZREF*DTDZ/TA/U/U
A1=ALOB(ZREF/ZO)
IFtB.LT.O.) BOTO 1
XH=0.05
X=0.15
FH=XH/((A1+.33333)/I.33333)**2-B
1O IF
FS=-5*X*PH
IF(ABS(F)-LT.O.OOO1) BOTO 1OO
IF(X.ED.XH) BOTO 100
SL=(F-FH)/(X-XH)
BINT=F-SL*X
XNEW=-BINT/SL
IF(ABS(XNEW-X).LT.0.0001) BOTO 10O
XH=X
FH=F
X=XNEW
BOTO 1O
1 XH=O.
FH=-B
X=-.05
20 IFCX.BE..O6667) X=.O6666
IF(X.E0.0.) X=0.0001
PH=1.O/(1.0-15*X)**.25
ZETA= <1.0-15*X)**. 25
ZETAO=(1.O-15*X*ZO/ZREF> **. 25
ARB1=ALOS((ZETA-1.O)*(ZETA-H.O)/((ZETA+1.O)*(ZETAO-1.O)))
ARB2=2.O*(ATAN(ZETA)-ATAN < ZETAO))
F=X/((Al-PS)/PH)**2-B
IF(ABS(F).LT..OOO1) BOTO 1OO
SL=tF-FH)/(X-XH)
BINT=F-SL*X
XNEW=-BINT/SL
IF(ABS(XNEW-X).LT.0.0001) BOTO 1OO
XH=X
FH=F
X=XNEW
BOTO 20
100 IF(X.LT.O-) ZOL=X
IF(X.BE.O.) ZOL=X/(1.0-5*X)
USTAR=0.35*U/(Al-PS)
IF(ZOL.LT.O. ) PHH=O.74/U.-9*ZOL)**.5
IF(ZOL.BE.O.) PHH=.74+5*ZOL
EDDY=.35*USTAR*Z/PHH
RETURN
END
C-43
-------
APPENDIX D
TEST RUNS AND SAMPLE INPUT FILE
D-l
-------
-------
n
o
SAMPLE INPUT FILE
o
n
o o \n -a -a -jj -a in 't v> n r> M * * -t * * * in * "> n o o o o o o o o o o o o o o o o o o o o o o o o
O 0~ O- N O O Cf- 0- N O O O O O '-' O O O O O O O O O O O O O O O O O O O
^ OQ i- in o- 6694 (
t
O
o
,~.
CN
CD
o'
,3
O
00
O'
d n
o
o
o
'H d
O 1'"'
o o
n
CO
CN
.--I
N
CN
i-^i
0-
-0
-0
,-.
o
d
Ğr
<
o
o
d
o in in iii o o in in in in o o o o o
in n CN n CN in in * in in in oo o o- o~
CN CN
00 O
d d
MinQJa>oDin*
-------
M
>
_l
s
a
Q
a. *
N. CL
CC D
I O
1 CC
t CD
M Ul
*
*
*
1-
(J
Ul
n
a
cc
a.
*
CD 1
a a * a
CC M
u cc cn cc
Ğ O UI CD
E U. U
~ cc cc
* a o
Z M O 1-
o ui a.
w CC UJ
t- 3 _J U
~ CC h-
Z 3 U.
a cc
U CD * O
Z U.
Ul >-.
CD a *
^ Z
CC Ul
UJ
> *
^
cc
n
o
I
Ğ*
M
>
>
i
01 in
w CC 1
X Ul
E 1
1
o o M in N
o o in N w
O O 03 M 111
O O N -0 00
O O 0 M M
-<ĞĞ*
o * in M o
O O * N CO
O O CO -0 O
O O 0ğ M *
O O N M O
' Ğ
-------
Q
Q. #
j-: x
> tC
_1 I O
i i
*
a *
.
1
1-4
^[
a
if
1
1
1
1
1
1
1
1
1
1
1
1
O 1
O 1
0
0 1
o-
1
1
1
o
1
o
O 1
in
t^ i
i
i
**
i
t- tn
o
i
s: o
a o i
r n
*" * r^ \
X
r
i
*
i
O 1
0 1
o
111 1
1
1
1
1
X X
1
tn 5
w CC 1
x Ul
Z 1
o o in <* o-
o o r-i N *
o o w N in
O 0 -0 * CO
o o o Ğ -*
-< * M
O bl -< O N
o o * * o
O O N 0- O-
O O O M tN
o o CM * r-j
^ 10 tN
X X X X X
o o o o o
o o o o o
0 O O O O
Q O O O O
M -H O ff> CD
Ğğ4 i-Ğ Ğ-l
SAMPLE OUTPUT FILE
ALTERNATIVE 2
Constant-K, Detailed Model
D-5
f>
ft
-------
Q
0. *
X >. 0.
> IT 13
_l I O
M I <£
1*0
Q PI CD
*
*
^_
CJ
UJ
a
a:
a.
a.
UJ
cc
0
u.
UJ
in
u
t-
m
UI
t-
*
£
~
UI
t-
UJ
E
U
hH
IB
3 *
U
in
E
(C> *
(D *
E 1
1
O o Ch - -<
O O M O -0
O O CC CD CD
o o in o i*>
o o h- in CD
M IN
O M N M CO
o o - * -*
o o o ft n
O O Q *-> "<
o o o- o- o
CM CN
*X "N, ^ X, "X
o o o o o
o o o o o
o o o o o
o o o o o
M -. O 0- tD
*-l ^-1 ^-Ğ
SAMPLE OUTPUT FILE
ALTERNATIVE 3
Variable-K, Linear Model
D-6
Ik *
-------
Q
Q. *
x -x a.
> CC U
Q CM CD
*
*
^
CJ
UJ
1-3
a
tr
0-
*
CD *
a:
o
i
i
C^
J-
_J
1^
Q
41
O
O
o
o
0-
o
o"
o
111
^
I-
-O
-0
hi
cn
_j
-------
-------
TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing/
1. REPORT NO.
EPA-450/4-86-003
I. RECIPIENT'S ACCESSION NO.
4. TITLE AND SUBTITLE
Continued Analysis and Derivation of a Method to
Model Pit Retention
5. REPORT DATE
January 1986
6. PERFORMING ORGANIZATION COOS.
7. AUTHOR(S)
K. D. Winges
C. F. Cole
8. PERFORMING ORGANIZATION REPORT NO.
9. PERFORMING ORGANIZATION NAME AND ADDRESS
TRC Environmental Consultants, Inc.
7002 South Revere Parkway, Suite 60
Englewood CO 80112
10. PROGRAM ELEMENT NO.
11. CONTRACT/GRANT NO.
68-02-3886
12. SPONSORING AGENCY NAME AND ADDRESS
13. TYPE OF REPORT AND PERIOD COVERED
Monitoring and Data Analysis Division
Office of Air Quality Planning and Standards
U. S. Environmental Protection Agency
Tn'anlp Park. NP. 77711 _
14. SPONSORING AGENCY CODE
EPA/200/04
15. SUPPLEMENTARY NOTES
ng
NO
Project Officer: J. S. Touma
16. ABSTRACT
This report summarizes the results of a continuing effort to better understand the
dispersion and transport of particulate matter released within surface coal mines.
The report examines the relationship between critical meteorological parameters in
an effort to refine an existing model algorithm to determine escape fraction. Methods
to incorporate calculating particulate matter escape fraction into a regulatory
air quality model are proposed and FORTRAN program listings of four alternatives
are included.
17.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.IDENTIFIERS/OPEN ENDED TERMS C. COSATI Field/Group
Air Pollution
Coal Mining Emissions
Particulates - Escape Fraction
Meteorology
18. DISTRIBUTION STATEMENT
19. SECURITY CLASS (ThisReport)
Unclassified
21. NO. Of PAGES
133
Unlimited
20.
SECURITY CLASS (This page)
unclassified
22. PRICE
EPA Form 2220-1 (Rtv. 4-77) PREVIOUS EDITION is OBSOLETE
-------
INSTRUCTIONS
1. REPORT NUMBER
Insert the EPA report number as it appears on the cover of the publication.
2. LEAVE BLANK
3. RECIPIENTS ACCESSION NUMBER
Reserved for use by each report recipient.
TITLE AND SUBTITLE
"itle should indicate clearly and briefly the subject coverage of the report, and be displayed prominently. Set subtitle, if used, in smaller
ye or otherwise subordinate it to mam title. When a report is prepared in more than one volume, repeat the primary title, add volume
mber and include subtitle for the specific title.
5. REPORT DATE
Each report shall carry a date indicating at least month and year. Indicate the basis on which it was selected (e.g., date of issue, date of
:py -oval, date of preparation, etc.).
6. PERFORMING ORGANIZATION CODE
Leave blank.
7. AUTHOR(S)
Give name(s) in conventional order (John R. Doe, J. Robert Doe, etc.). List author's affiliation if it differs from the performing organi
zation.
8. PERFORMING ORGANIZATION REPORT NUMBER
Insert if performing organization wishes to assign this number.
9. PERFORMING ORGANIZATION NAME AND ADDRESS
Give name, street, city, state, and ZIP code. List no more than two levels of an organizational hirearchy.
10. PROGRAM ELEMENT NUMBER
Use the program element number under which the report was prepared. Subordinate numbers may be included in parentheses.
11. CONTRACT/GRANT NUMBER
Insert contract or grant number under which report was prepared.
12. SPONSORING AGENCY NAME AND ADDRESS
Include ZIP code.
13. TYPE OF REPORT AND PERIOD COVERED
Indicate interim final, etc., and if applicable, dates covered.
14. SPONSORING AGENCY CODE
Insert appropriate code.
15. SUPPLEMENTARY NOTES
Enter information not included elsewhere but useful, such as: Prepared in cooperation with, Translation of, Presented'at conference of,
To be published in, Supersedes, Supplements, etc.
16. ABSTRACT
Include a brief (200 words or less) factual summary of the most significant information contained in the report. If the report contains a
significant bibliography or literature survey, mention it here.
17. KEY WORDS AND DOCUMENT ANALYSIS
(a) DESCRIPTORS - Select from the Thesaurus of Engineering and Scientific Terms the proper authorized terms that identify the major
concept of the research and are sufficiently specific and precise to be used as index entries for cataloging.
(b) IDENTIFIERS AND OPEN-ENDED TERMS - Use identifiers for project names, code names, equipment designators, etc. Use open-
ended terms written in descriptor form for those subjects for which no descriptor exists.
(c) COSATI FIELD GROUP - Field and group assignments are to be taken from the 1965 COS ATI Subject Category List. Since the ma
jority of documents are multidisciplinary in nature, the Primary Field/Group assignment(s) will be specific discipline, area of human
endeavor, or type of physical object. The application! s) will be cross-referenced with secondary Field/Group assignments that will folio
the primary posting(s).
18. DISTRIBUTION STATEMENT
Denote releasability to the public or limitation for reasons other than security for example "Release Unlimited." Cite any availability tc
the public, with address and price.
19.&20. SECURITY CLASSIFICATION
DO NOT submit classified reports to the National Technical Information service.
21. NUMBER OF PAGES
Insert the total number of pages, including this one and unnumbered pages, but exclude distribution list, if any.
22. PRICE
Insert the price set by the National Technical Information Service or the Government Printing Office, if known.
EPA Form 2220-1 (Rev. 4-77) (Reverse)
------- |