&EPA
United States
Environmental Protection
Agency
Environmental Monitoring
and Support Laboratory
P.O. Box 15027
Las Vegas NV89114
EPA-600/4-78-030
June 1978
Research and Development
Environmental
Monitoring Series
Optimum Meteorological
and Air Pollution
Sampling Network
Selection in Cities
Volume I:
Theory and Design
for St. Louis
-------
RESEARCH REPORTING SERIES
Research reports of the Office of Research and Development, U.S. Environmental
Protection Agency, have been grouped into nine series. These nine broad categories
were established to facilitate further development and application of environmental
technology. Elimination of traditional grouping was consciously planned to foster
technology transfer and a maximum interface in related fields. The nine series are:
1. Environmental Health Effects Research
2. Environmental Protection Technology
3. Ecological Research
4. Environmental Monitoring
5. Socioeconomic Environmental Studies
6. Scientific and Technical Assessment Reports (STAR)
7. Interagency Energy-Environment Research and Development
8. "Special" Reports
9. Miscellaneous Reports
This report has been assigned to the ENVIRONMENTAL MONITORING series/This series
describes research conducted to develop new or improved methods and instrumentation
for the identification and quantification of environmental pollutants at the lowest
conceivably significant concentrations. It also includes studies to determine the ambient
concentrations of pollutants in the environment and/or the variance of pollutants as a
function of time or meteorological factors.
-------
EPA-600/4-78-030
June 1978
OPTIMUM METEOROLOGICAL AND AIR POLLUTION SAMPLING NETWORK SELECTION
IN CITIES
Volume I: Theory and Design for St. Louis
Fred M. Vukovich, Walter D. Bach, Jr.,
and C. Andrew Clayton
Research Triangle Institute
P. 0. Box 12094
Research Triangle Park
North Carolina 27709
Contract No. 68-03-2187
Project Officer
James L. McElroy
Monitoring Systems Research and Development Division
Environmental Monitoring and Support Laboratory
Las Vegas, Nevada 89114
ENVIRONMENTAL MONITORING AND SUPPORT LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
LAS VEGAS, NEVADA 89114
-------
DISCLAIMER
This report has been reviewed by the Environmental Monitoring and
Support Laboratory-Las Vegas, U.S. Environmental Protection Agency, and
approved for publication. Approval does not signify that the contents
necessarily reflect the views and policies of the U.S. Environmental
Protection Agency, nor does mention of trade names or commercial products
constitute endorsement or recommendation for use.
ii
-------
FOREWORD
Protection of the environment requires effective regulatory actions
which are based on sound technical and scientific information. This infor-
mation must include the quantitative description and linking of pollutant
sources, transport mechanisms, interactions, and resulting effects on man
and his environment. Because of the complexities involved, assessment of
specific pollutants in the environment requires a total systems approach
which transcends the media of air, water, and land. The Environmental
Monitoring and Support Laboratory-Las Vegas contributes to the formation and
enhancement of a sound monitoring data base for exposure assessment through
programs designed to:
• develop and optimize systems and strategies for moni-
toring pollutants and their impact on the environment
• demonstrate new monitoring systems and technologies by
applying them to fulfill special monitoring needs of
the Agency's operating programs
This report discusses the theoretical bases for a method of designing
meteorological and air quality monitoring networks and the application of
the method to the metropolitan St. Louis area. Regional or local agencies
may find this method useful in planning or adjusting their aerometric moni-
toring networks. The Monitoring Systems Design and Analysis Staff at the
EMSL-LV may be contacted for further information on the subject.
^
George B(. Morgan
l Director
Environmental Monitoring and Support Laboratory
Las Vegas
iii
-------
PREFACE
This document deals with the theoretical aspects of a methodology for
optimum meteorological and air quality monitoring network selection and the
application of the methodology to the metropolitan St. Louis area. Subse-
quent reports will be concerned with verification of the methodology for
St. Louis with regard to the airflow and with regard to the air quality.
James L. McElroy
Project Officer
Environmental Monitoring and Support Laboratory
Las Vegas
iv
-------
ACKNOWLEDGEMENTS
This report was prepared by the Research Triangle Institute (RTI),
Research Triangle Park, North Carolina, under Contract No. 68-03-2187
for the U.S. Environmental Protection Agency (EPA). A major portion
of the funding of Phase I of this research project came from the
National Science Foundation, Research Applied to National Needs. The
support of these agencies is acknowledged. The Project Officer for
EPA was Dr. James L. McElroy. Many individuals from RTI participated
in this project. They are too numerous to cite here. However, cer-
tain individuals stand out and their efforts shall be acknowledged.
Mr. J. W. Dunn was responsible for developing the computer algorithm
for the three-dimensional primitive equation model. Mr. Bobby
Crissman developed all the major parameters for the primitive equation
model and helped in running the field program to verify the results
of the entire modeling effort. Mr. Clifford Decker was responsible
for the management of the field program.
Special acknowledgements are to be cited for the cooperation of
various individuals and organizations allowing us to set up field sta-
tions on their property in St. Louis, Missouri. These are Sister
Rosita Hyland of Incarnate Word Academy; Mr. E. M. Krings of the Chan-
cery Office in St. Louis who allowed us to set up the site of Kenrick
Seminary; Mr. Paul Mydlar of the St. Louis City Air Pollution Control
Agency; and Mr. James Spanos, District Engineer of the East Side Levy
and Sanitary District. Special thanks to Mr. Ashwin Gajjar of the
St. Louis County Air Pollution Control Agency for not only allowing
us to place our instruments in some of his air monitoring stations
but also for providing us air pollution and meteorological data.
We would also like to acknowledge the cooperation of Mr. Robert
Browning of the EPA, Research Triangle Park, for providing us with
Regional Air Pollution Study (RAPS) data which will be used in a
major part of our data analysis to follow.
-------
TABLE OF CONTENTS
Foreword ill
Preface iv
Acknowledgements v
Table of Contents vii
List of Figures viii
List of Tables ix
Glossary of Terms
Primitive Equation Model X
Regression Model and Site Selection Algorithm xii
Objective Analysis Model xiv
Summary xvi
1.0 Introduction ' 1
2.0 Theory 4
2.1 Applications of Model 6
3.0 Model Development 9
3.1 Primitive Equation Model 9
3.1.1 Grid, Finite Differencing and Boundary Conditions 11
3.1.2 Boundary Layer Parameterization 14
3.1.3 Initial Conditions and Integration 20
3.2 Determination of a Wind Monitoring Network 20
3.2.1 General 20
3.2.2 Development of Statistical Model Form 23
3.2.3 Selection of Candidate Sites 31
3.2.4 Design Selection Criteria 34
3.2.5 Design Selection Algorithm 39
3.2.6 Application of the Algorithm 41
3.2.7 Comparison of Alternative Wind-Monitoring
Networks 44
3.3 Objective Variational Analysis Model 59
3.3.1 Description of the Model 62
3.3.2 Analysis Technique 67
3.3.3 Numerical Solutions 69
3.3.4 Weighting Factors 71
3.3.5 Test Results 72
4.0 Test and Evaluation 84
4.1 Operational Optimum Network 84
4.2 Field Programs 87
4.3 Analysis of Data 89
Appendix A - A Theoretical Study of the St. Louis Heat Island:
The Wind and Temperature Distribution 91
Appendix B - Area Sources with Objective Variational Analysis Model 116
Appendix C - Sensitivity Analysis and Testing of the Variational
Objective Analysis Model 121
References 133
vii
-------
LIST OF FIGURES
Number Page
1 Flow diagram for developmental concept 7
2 Model grid points in the metropolitan St. Louis area 12
3 Distribution of H(x,y) (m) field for the St. Louis, Missouri
region 15
4 Distribution of the 66(x,y) (°K) field for the average heat
island in St. Louis, Missouri 16
5 Design selection strategy 21
6 Thirty-point designs from square grid and diamond grid 43
7 Variance function for design AO 48
8 Variance function for design Al 49
9 Variance function for design A2 50
10 Variance function for design A3 51
11 Stations in the A2 network (the OSN) for St. Louis, Missouri 53
12 Variance function for the OSN* 58
13 'Concentration isopleths (yg/m3) at ground level from emissions
(10 gm/s) normally distributed within the shaded area, at
neutral stability and a constant wind, U, of 3 m/s 73
14 Horizontal distribution of £( ), S (— • — •), and y-1
( ) along y = 5 km for Case I 74
15 Concentration distribution (yg/m3) with x as determined using
R~2 interpolation of randomly selected observations (*) from the
distribution of Figure 12 78
16 Concentration distribution (yg/m3) with x for Case VI showing
i|>* (— • —), ty(x x),i(i from Case II ( )* along y = 5 km 79
17 Analysis error variance versus iteration in Case VI ( ) and
Case VII ( ). _2 81
18 Concentration isopleths (yg/m3) determined from an R inter-
polation of perturbed observations at randomly selected data
points (*) of the distribution of Figure 25. 82
19 Concentration distributions with x for Case VII showing ip*
(— • — •)> ?(' — ') and iK ) along y = 5 km 83
20 Station array in the RAPS network and in the city/county
network in St. Louis, Missouri 85
21 Stations in the OSN* for St. Louis, Missouri 86
22 Location, name and number of stations maintained by RTI 88
viii
-------
.LIST OF TABLES
Number
1 Terms Included in 49 Candidate Models 28
2 Means of R2 and S2 Over 48 Cases for 49 Models 29
3 Correlations Between Simulated and Predicted Values for
Selected Models 30
4 Four 19-Point Networks Selected for Evaluation 45
5 Average Correlation and Sums of Squares of Deviations for
Each Network 45
6 Symbol Key to Variance Function Plots 46
7 Summary of Correlations Between Observed and Predicted
Values Over 441 Points of Innermost Grid—By Wind Shear
and Design 56
8 Specification of Design Points for The OSN and OSN* 57
ix
-------
GLOSSARY OF TERMS
Primitive Equation Model
x ~ horizontal coordinate in a Cartesian coordinate system which is positive
in the direction of the initial flow field
y ~ horizontal coordinate perpendicular to the x axis in the Cartesian
coordinate system which is positive to the left of the initial flow field
z ~ vertical coordinate
u ~ x-component of the horizontal wind field
v ~ y-eomponent of the horizontal wind field
w ~ vertical velocity
p ~ pressure
0 ~ potential temperature
po ~ 100 centibar
R ~ gas constant for dry air
f ~ Coriolis parameter
p ~ density
g ~- the acceleration of gravity
AX ~ x-component of the exchange coefficient for heat and momentum
Ay ~- y-component of the exchange coefficient for heat and momentum
A ~ vertical component of the exchange coefficient for the x-component of
xz
momentum
A ~ vertical component of the exchange coefficient for the y-component of
yz
the momentum
A.. ~ vertical component of the exchange coefficient for heat
( )' ~ perturbation quantities
( ) ~ initial values
H(x,y) ~ the height of the natural and manmade topography above the river height
at a grid point
z* ~ the height of any grid point above H
A9 ~ the temperature associated with the heat island at any time
66 ~ the total amplitude of the heat island
Y ~ rate constant for the production of the heat island
-------
UA ~ x-component of the friction velocity
v^ ~ y-coraponent of the friction velocity
T* ~ the scale temperature
U* ~ J2~_2
Vu* + VA
z ~ roughness length
k ~ von Karman constant
K = 0.287
a = 18
£ ~ mixing length for momentum
X, ~ mixing length for heat
n.
Ri ~ Richardson number
ai = 18
ct = 11
xi
-------
GLOSSARY OF TERMS
Regression Model and Site Selection Algorithm
x~horizontal coordinate in a Cartesian coordinate system; increasing from
west to east
y~horizontal coordinate in a Cartesian coordinate system; increasing from
south to north
H(x,y)~ height (above Mississippi River) of the natural and manmade topography
at the point (x,y)
P.(x,y)~a polynomial term (identified by the index j) of the form x y where a
and Y take on non-negative integral values
Z, (x,y)~the kth wind speed component (k=l for WE and k=2 for SN) at the point
U,y)
B,.~the estimated regression coefficient associated with Pj when component
Z, is used as a dependent variable
Q.(x,y,H)~a polynomial term (identified by the index j) of the form x y [H(x,y)]
Cv.~the estimated regression coefficient associated with Q. when component
KJ J
Z^ is used as a dependent variable
Zjc(x,y)~the predicted value for wind component k at the point (x,y)
s ., — the residual variance for wind component Z^ from fitting a particular
model (index m) to a particular set of data (index i)
9
R ., •— the proportion of total variation in wind component Z accounted for
by a particular model (index m), when applied to a particular set of
data (index i)
b.(k)~the jth estimated regression coefficient in a particular (13-term)
model which relates wind component Z, to geographic location and
elevation
X^~ the matrix of independent variables for a particular design (network
of stations.) D
bv ~the vector of b (k) values obtained from a design D
" (D)
Z (x,y)~the predicted value of the kth wind component at the point (x,y)
k
obtained from applying a particular model to data from a design D
a 2~the (unknown) error variance associated with component Z,
v_.(k,x,y)~ the variance of Z, ^ '(x,y)
xii
-------
s\
s (x,y) = vr)(k,x,y)/a , , the value of the variance function associated with
D K
design D, evaluated at the point (x,y)
w(x,y)~ the value of a weighting function at the point (x,y)
R~ a set of 521 grid points contained in the region {(x,y): |x| _<_ 16,
|y| <_ 16}
T = • 2-i w(x)y)s^(x,y), the value of the design selection criterion
for design D
A_~a class of r-point designs
D „ ~a particular element (identified by the index 1) of A
Af f
xiii
-------
GLOSSARY OF TERMS
Objective Analysis Model
x ~ horizontal coordinate in a Cartesian coordinate system; increasing from
west to east
y ~ horizontal coordinate in a Cartesian coordinate system; increasing from
south to north
z ~ vertical coordinate in a Cartesian coordinate system increasing with
distance from earth's surface
u — x-component of the horizontal wind vector
v ~ y-component of the horizontal wind vector
w ~ vertical velocity
fy ~ pollutant concentration
Q ~ pollutant emission rate
a ~ standard deviation of the pollutant concentration in the horizontal
7
crosswind plane from a point-source emission
0" ~ standard deviation of the pollutant concentration in the vertical plane
Z
from a point-source emission
U ~ mean wind speed
4> ~ an arbitrary, unspecified variable
ft ~ the domain of the integration; relaxation coefficient
K-- ~ eddy diffusivity in the horizontal plane
K ~ eddy diffusivity in the vertical plane
Z
Z ~ maximum altitude of the surface layer
S ~ local standard deviation of pollutant concentration in the vertical
Z
plane
L ~ characteristic length of the pollution system
E ~ characteristic value of o
z
q ,q ~ standard deviation of the area emission distribution
x y
R ~ turbulent Reynolds number (= K/UL)
(-) ~ characteristic magnitude of a variable
z
dz
f
- = /
JQ
(*) *- analytic solution to the pollution distribution
(^) ~ "observed" quantity
xiv
-------
(ij; ) ~ "background" value
B
i,l~ indices of increments in the x direction
j,ra ~ indices of increments in the y direction
v ~ iteration number
-------
SUMMARY
This report describes the theory and the models which have been
developed to establish an optimum meteorological and air pollution sampling
network in urban areas. The theory and models are generalized so that they
can be applied to any urban area. However, the network developed is specific
for a given urban region. The network yields both meteorological and air
pollution data so that not only may the distribution of a pollutant be
obtained over the domain of the network, but that the distribution may be
predicted over short time intervals.
The optimum network provides a basis whereby the distribution of any
given pollutant may be obtained at a given time over the domain of the
network. The data may be used to determine long-term statistics or short-
term analysis. Long-term statistics may be determined from the background
or maximum pollutant concentrations. The short-term analysis will provide a
means whereby regions may be determined where a pollutant concentration
exceeds the standards.
Three specific models are required in order to determine the optimum
network and acquire the air pollution distribution over the domain of the
network. These are: a three-dimensional hydrodynamic model; a statistical
model; and an objective variational analysis model.
The basis of the network is the wind field in the urban area rather than
the air pollution distribution. The wind field was chosen because it
provided a solution with longer-term stability than the air pollution distri-
bution. The addition of two or three major sources in an urban area can
markedly changed the statistics of the air pollution distribution in a short
period of time. Any network based on today's air pollution distribution may
become inadequate with the addition of major sources.
The primitive hydrodynamic model is used to simulate the wind field over
a variety of cases in which the major factors influencing the flow (essen-
tially the urban heat island, surface roughness variability, topography,
etc.) are included in the model. With proper definition of the urban forcing
functions, a reasonable approximation to the statistical nature of the flow
xvi
-------
can be obtained. By statistical nature of the flow, we mean that the
specification of the exact magnitude of the flow is not necessary if its
order of variability is reasonably well represented. This further asserts
that the magnitude of the difference between minima and maxima in the flow
may not be exact, but the location of minima and maxima be reasonably well
represented. More specifically, the simulation data were used to develop the
form of a statistical model which could be used to describe the wind distri-
bution over the domain of the network. The magnitude of the flow and wind
speed maxima and minima are determined after the parameters of the model are
estimated. In the applied sense, observed wind speeds and directions from
stations in the network would be used to estimate parameters (i.e., actual
data will be used to determine the magnitude of the flow and of the wind
speed maxima and minima).
These simulated data were used to determine the form of a regression
model which approximates the various wind fields generated by the hydro-
dynamic model. The regression model form was then used, along with a set of
potential network sites and a criterion for judging alternative networks
(e.g., minimizing the average variance of predicted values), to derive an
"optimum" network for sampling winds. The method used to develop the network
involved the successive elimination of candidate sites—until a reasonably
sized network was achieved.
The air pollution distribution is obtained through an objective varia-
tional analysis model. The model simultaneously minimizes the error variance
by comparing the observed pollution concentration field with the derived
pollution field and the error variance of the constraint equation. In order
for the variance to be a minimum, it is shown that, a second-order Euler-
Lagrange equation must be satisfied under specific boundary conditions. The
model requires as input parameters: the wind field derived from the optimum
sampling network, the observed pollutant concentration data from the optimum
sampling network, and the source emissions.
The models were applied to develop an optimum network for St. Louis,
Missouri. St. Louis was chosen as a test case because the Environmental
Protection Agency's (EPA) Regional Air Pollution Study (RAPS) provided
xvii
-------
corollary data to evaluate the optimum network. After the optimum network was
determined, existing stations were included as part of the optimum network
which coincided or nearly coincided with those in the RAPS or city/county
network. Of the nineteen (19) stations in the optimum network, sixteen (16)
coincided or nearly coincided with existing stations.
In the summer of 1975 and winter of 1976, field programs were conducted
concurrent with a RAPS intensive study period to collect data in order to
test and evaluate the optimum network. Besides the sixteen (16) stations that
overlapped with the city/county and RAPS network, three (3) stations had to be
set by the Research Triangle Institute. At each of these stations, wind speed
and direction, carbon, monoxide, methane, and total hydrocarbon data were
collected. Carbon monoxide was considered the principal pollutant to be used
in the evaluation of the network. The results of this evaluation will be
presented in a future report.
rvili
-------
1.0 INTRODUCTION
The establishment of the air quality standards has brought about the
need to develop air quality sampling networks in urban areas. Such networks
have two primary purposes. The first is to provide data for analysis of air
pollution over the urban region to determine if and where a particular pollu-
tant concentration exceeds the standard and to determine the effectiveness
of long-term control strategies. The former information is required to
implement control measures to reduce the concentration of the particular pol-
lutant to an acceptable level. The second purpose is to provide a data base
for both the short- and long-term prediction of the concentration of a par-
ticular pollutant over the urban region. This information provides a base
for the development of long-term control strategies. Diagnosis and predic-
tion of air pollution is particularly necessary now that the energy crisis
has made it probable that some industries and power plants will burn high
sulfur coal.
Sweeney (1969) attempted to establish objective guidelines for deter-
mining the number and location of sampling stations in an optimum sampling
network (OSN) for urban areas purely by statistical means. He concluded that
the present state-of-the-art of statistical theory does not provide a basis
for a general solution to the problem. The basis for this conclusion lies in
the fact that the air pollution distribution and meteorological conditions
are dynamic phenomena that vary in space and time in an urban region and vary
from one urban region to the next. However, Sweeney indicated that a solu-
tion specific for a particular urban region can be obtained if a representa-
tive and statistically significant data set is available for the urban
region. This data set should be made up of observations of n number of pollu-
tants and m number of meteorological parameters. For most urban regions, n is
greater than or equal to three, and m is at least one (the wind distribu-
tion) but can be larger, depending on the accuracy required.
In order to acquire this data set, Sweeney suggested that a high reso-
lution sampling network be set up in a given urban region for the collection
of data on pollution concentrations and meteorological conditions over a
period of years. Though this is scientifically appealing, it is economically
unfeasible. The Regional Air Pollution Study (RAPS) being sponsored by the
-------
Environmental Protection Agency (EPA) in the metropolitan St. Louis area will
cost millions of dollars and does not yield the high resolution data required
by Sweeney. Sweeney required a resolution of the order of 2 to 4 km. To
obtain data at those resolutions would require from 64 to 256 sampling
stations. It is unreasonable to ask each urban region to support such
programs.
In order to ease the economic pressure, Sweeney suggested that the urban
region should be divided into sections. A high resolution network should be
set up in one section and data collected over a number of years. When
sufficient data are obtained in the first section, the sampling network would
be moved to another, and data collection would commence and continue until
sufficient data were obtained. This iterative process would continue until
the data set was complete.
There is, however, a dilemma that arises when one determines an OSN
based solely, or in part, on air pollution data. Most, if not all, urban
regions are constantly changing. New industries are moving in, old industries
are moving out; that is, new point sources move in and old point sources move
out. It is important to note that the location of the old and new point
sources are very often different. Furthermore, over the years, rural roads
which often have insignificant traffic, become major line sources with the
establishment of housing projects and accompanying commercial areas near or
along the road. Other kinds of urban expansion results in similar changes.
Urban expansion can produce a major change in the air pollution distribution
in a 5- to 10-year period. Any OSN based on today's representative air
pollution data sets may soon become obsolete. This effect would prevent
Sweeney's suggested technique for obtaining a representative data set from
being practical. In the intervening time between collecting data in the first
segment of the urban region and collection of data in the last segment, major
changes in source configurations may have occurred, making the segmented data
sets totally independent of each other.
It is obvious for economic reasons that the basis of the OSN must be very
stable over a relatively long period of time. A technique to establish an OSN
on a more stable basis than is provided by the air pollution distribution is
described in this report. The models developed have been completed and an OSN
-------
for St. Louis, Missouri has been created. St. Louis, Missouri was chosen as
a test case because the RAPS program will provide corollary data to evaluate
the OSN.
-------
2.0 THEORY
The conservation equation for any pollutant, i|/, defined for a point in
space away from the earth's surface, yields the exact nature of the variables
necessary to completely specify the distribution of the pollutant, i.e.,
= - V • 7* + V • (KVifO + ty L, k $ + S . ,
n n
n
The first two terms on the right yield the meteorological influence on the
distribution of ty. The first term is the advective term where V is the wind
vector, and the second term is the diffusion term where K is the exchange
coefficient. The third term is the chemical reaction terra where k is the
reaction constant, i[i* is the pollutant with which reactions take place in the
case of second-order reactions, and n is the index of the summation indicating
that more than one, second-order reaction may take place. The last term is
the source (S) term. In this development, it was assumed that the pollutant
E*
k i|> =0.
For a given urban region, the source parameters [source strength, source
location, source type (line, point, etc.)] are known or can be determined.
Any changes in the source parameters due to urban expansion or industrial
relocation can also be specified. In the case of an inert gas, therefore, the
unknowns are the meteorological parameters, i.e., the wind field and the
distribution of K. Since a reasonable estimate of the diffusion coefficients
can be obtained through wind observations and estimates of the stability, it
»
is only necessary to establish the wind field.
The wind field should be relatively unaffected by urban expansion and
industrial relocation in comparison to the effect on the pollution distri-
bution. The changes in surface roughness that occur when a rural area is
converted into a residential area will be important, but will not completely
alter the nature of the flow. For example, winds from the north over a
formerly rural area will have a slight westerly component over the rougher
residential area, but the major wind component will be from the north; that
is, something as significant as the wind turning and becoming southerly due
to the roughness will not occur.
-------
Urban expansion increases the areal dimensions of the urban heat island,
which is probably the major forcing function affecting the mesoscale flow in
an urban region. However, this -expansion should be an effect manifested over
many years (20 to 30 years) and would not alter the flow significantly if the
location and heat content of the center of the heat island did not change
markedly. The expansion of the suburban residential area would tend to
flatten the temperature gradient in the suburbs, which would reduce accel-
erations due to pressure forces rather than intensify them. The establishment
of another heat island center would be an important factor which would affect
the flow, but the creation of such a center with an affective areal extent in
most cases would also take many years (again 20 to 30 years). This would
indicate that the urban wind field should be stable in the mesoscale over a
relatively long period of time compared to changes that could take place in
the air pollution distribution due to urban expansion.
If the parameters in the conservation equation are specified over an
urban region, then it is possible to determine the distribution of a parti-
cular pollutant using an Objective Variational Analysis Model (OVAM)
(Sasaki, 1970). This technique provides the pollution distribution by making
proper use of limited observed data and the governing equation. In terms of
the air pollution distribution, the governing equation (the conservation
equation), in effect, allows the sources to be treated similarly to observa-
tions of the particular pollutant and thus, a higher order variability of the
particular pollutant can be ascertained than would otherwise be obtained using
limited observed data. Therefore, an OSN for the wind field would suffice to
determine the distribution of a particular pollutant if: 1) air pollution
data were also available at the stations in this network, 2) a proper source
inventory were available, and 3) the above data (1 and 2) including the wind
field, were incorporated in an OVAM. These data may also be used to predict
the air pollution distribution over a short period of time.
However, we are still left with the prospect of establishing a represen-
tative data set for the wind field in order to determine an optimum network.
A numerical simulation of the wind field in an urban region under the action
of a multitude of synoptic meteorological conditions circumvents the economic
problem of high resolution data collection over long periods of time.
-------
With proper definition of the urban forcing function, a reasonable approxi-
mation to the statistical nature of the flow can be obtained. By statistical
nature of the flow, we mean that specification of the exact magnitude of the
flow is not necessary if the order of variability is reasonably well repre-
sented. This further asserts that the magnitude of the difference between
minima and maxima be reasonably represented. The simulation data would be
specifically used to develop the form of a statistical model which could be
used to determine the flow over the domain of the network. Actual magnitudes
for the wind speed components and between wind speed maxima and minima are
determined after the parameters of the statistical model have been estab-
lished. In the applied sense, parameters of the statistical model would be
determined using observed data from stations in the network. Since only the
form of the statistical model is developed using the simulation data, it is
not necessary that the simulation data precisely describe the exact magnitude
of the flow. Only the order of variability must be well represented. These
simulated data can then be used to obtain an approximation to the OSN by
integrating them with a model based on regression and sampling theory.
2.1 Applications of Model
Figure 1 is a flow diagram for the entire development concept. The
primitive equation model is used to establish a representative data set for
the wind field in a given urban area. In the case of St. Louis, Missouri,
the data set was composed of forty-eight (48) case studies. The data are used
to determine the form of a regression model that adequately characterizes the
wind field for the given urban region. The resultant model form is then used
in the site selection algorithm.to yield the OSN. The theoretical development
phase of the modeling ends here; that is, in the figure the theoretical
development phase is the region above the lower dashed line.
The applied phase is the region below the upper dashed line. For this
phase, it is assumed that the OSN for the wind field has been set up in the
urban region. Wind and air pollution data are obtained at the stations in
the OSN. The wind data are used to develop the wind distribution through
the same regression model determined in the theoretical development phas«.
The wind distribution, the observed air pollution data, and the emissions
-------
__ ___
Regression
Model
\
f
Site Selection
Algorithm
— —
_ _—
Op t imun
Sampling
Network
Air
FolluClon
Data from
SCations In
Network
Kind Data
from
Stations in
Network
Regression Model
Emissions
Inventory
Objective Variational
Analysis Model
Air
Pollution
Distribution
Theoretical Development
Applied Phase
Figure 1. Flow diagram for developmental concept.
-------
inventory are input data for the OVAM, which subsequently yields the air
pollution distribution.
-------
3.0 MODEL DEVELOPMENT
The principal purpose of this phase of the research project was to
develop the models necessary 1) to yield an approximation to the statistical
nature of the wind field, 2) to determine an OSN for the wind field, and
3) to establish the air pollution distribution using the wind distribution,
air pollution data from the OSN and the source parameters.
3.1 Primitive Equation Model
A three-dimensional primitive equation model was developed to simulate
the flow in an urban region. The model included those factors (essentially
the urban heat island, surface roughness variability, and topography—
including effective building heights) which would significantly change the
synoptic-scale flow. Solutions were obtained for forty-eight (48) synoptic
cases. These included variations in wind direction (8 points on the compass),
three wind speed distributions (weak, medium and strong wind speeds with cor-
responding weak, medium and strong vertical wind shear in the boundary layer),
and two categories of stability (stable and adiabatic boundary layer). These
solutions were specific for the city of St. Louis, Missouri in that the
topography and the parameters needed to generate the urban heat island were
obtained from data specific for that city. The city of St. Louis is located
in rather flat terrain and is not affected by large bodies of water. In order
to obtain solutions for other urban regions, it is only necessary to supply
the model with these parameters for the desired urban regions.
The forty-eight (48) cases mentioned above on which the optimum network
was based were made up of surface wind speed categories which range from
approximately 0.5 m/s to 5.0 m/s. There was no apparent increase in the order
of variability when the wind speed was further increased. Analysis also
showed that increasing the stability or decreasing the stability beyond the
range depicted had no appreciable affect on the order of variability. Since
the form of the statistical model which would be used to describe the flow
over the domain of the network using the data obtained from the network,
requires that the order of variability is well represented in the data set
used to determine the model form, it was concluded that the wind speed
-------
categories and stability categories were sufficient to yield a reasonable
estimate of the model form.
The basic equations governing the dynamics and thermodynamics for dry
convection in the primitive equation model are listed below:
3u 3u 3u JL" _ _ M _.-
(1)
3x x 3x 3y y 3y 3z xz 3z
3v 3v 3v 3v _ R6 /P_\K 3p f
3t U 3x V 3y W 3z p pQ' 3y ~ U
(2)
3x x 3x 3y y 3y 3z yz 3z '
(4)
P K
_ 12*
3z R9 vp ' (5)
10
-------
The model is dry. Long- and short-wave radiation is also not included.
In the above equations, the horizontal diffusion coefficients for momentum and
heat were considered identical. The local rate of density change was set to
zero (gT= 0)« This assumption yielded equation (4) for the vertical velocity.
The most important approximation used in the model is the hydrostatic
assumption through which the pressure field was derived. This assumption is
justified since we are dealing with a relatively shallow atmosphere; i.e.,
the top of the atmosphere is 4.0 km in the model. It is expected that the
pressure perturbations produced by the forcing functions in the urban areas
should be less than 10~^ g which suggests that the hydrostatic approximation
ia a reasonable approximation.
Preliminary computations of the vertical velocity using equation (4)
indicated that the contribution of the second term on the righthand side is
negligible. The simplified version of (4) was used in this study; i.e.,
2 2
8 w , 0 3p 3w . 3 p 3 , ,3u 3v. ,
P7T+ 2-3t 3l+W72 = - 3¥ {p(^+37)} (6)
dZ dZ
Perturbation principles were applied to the above equation. In the per-
turbation analysis, variables designated by the overbar (p, for instance) are
used to represent the synoptic state which was assumed to be both in geostro-
phic and hydrostatic equilibrium. Primed variables are used to represent the
perturbations produced by the available forcing functions in the city area.
3.1.1 Grid, Finite Differencing and Boundary Conditions
The horizontal dimensions of the grid are 144 km by 144 km, which more
adequately covers the urban and suburban regions of St. Louis. A nested grid
was employed in which the central region (21 x 21 grid points) centered about
the city had a grid spacing of 1 km. Outside the central region, the grid
spacing increased in both x and y according to the expression 2 km, where
i = 1, 2, 3, 4 until the grid spacing equalled 16 km. Then the grid spacing
was held constant for three grid points. Figure 2 shows a graphical descrip-
tion of the grid. The three equally spaced grid points 16 km apart are not
included in the figure.
11
-------
o
• ••••••••**• •/•
•=16.
.=&-
16
(EAST)
Figure 2. Model Grid Points in the metropolitan St. Louis area
12
-------
In regions where grid space was constant, center difference formulas were
used for first- and second-order space derivatives. In the region of changing
grid spacing, upstream differencing was employed to give first-order space
derivatives, and second-order, non-centered derivatives were used to determine
the diffusion terms. Since upstream differencing is dissipative, the hori-
zontal diffusion terms were adjusted in the region of unequal grid spacing so
that overdampening did not occur.
At the horizontal boundaries, mass, momentum, and heat were allowed to
flow in or out except at the upstream boundary; i.e.,
at x = - 72 km,
u' - v' = 8' = p1 = 0
at x = + 72 km,
3x 3x 3x 3x
at y = _+ 72 km,
3y ~ 3y ~ 3y ~ 3y
There were eight vertical levels in the model. The vertical levels
were: z = 0, 100 m, 300 m, 600 m, 1000 m, 1500 m, 2500 m and 4000 m.
Heights are with respect to the z = 0 plane which will be defined in
Section 3.1.2. Upstream differencing was used exclusively in the vertical.
The upper boundary conditions we-re
w1 = v1 = u1 - 0' = p1
and
u = u, p = p, and 6=9
13
z = 4 km
-------
The lower boundary conditions will be discussed in detail in Section 3.1.2.
3.1.2 Boundary Layer Parameterization
The primary forcing functions in urban areas are the surface heat flux
and local terrain effects. In order to describe the effect of terrain, the
topographic distribution was parameterized for the model. The average height
of the local terrain plus the average building heights were computed in areas
about each grid point. The building heights were determined from Sanborn
Maps, courtesy of the St. Louis City Planning Office. The size of the area
depended on the grid-spacing, i.e., for a 1-km grid-spacing, the area about
the grid point was 1 km , and for a 2-km grid-spacing, it was 4 km , etc. The
minimum averaged terrain height was determined from the set of average values
and this value was subtracted from each terrain height. Than the building
heights were added to yield average deviations from the reference height. The
z = 0 plane is defined where H(x,y) = 0.
At z = H(x,y) and at the side walls of an obstacle that may be higher
than the first grid point above z = 0 where the prediction equations are
integrated, the boundary conditions were
u = v = w = 0.
If z* - H(x,y) < 15 n, where z* is the height of any grid point above H,
the prediction equations become unstable for the time increment used.
Therefore, values of H(x,y) which did not satisfy this criteria were reduced
to H(x,y) = z* - 15. For St. Louis, all values of H(x,y) were less than 100 m
which is the height of the first level in the model; about six H(x,y) values
exceeded the stability criteria and were reduced (Figure 3).
In order to characterize the effects of an urban heat island and/or large
bodies of water (sea or lake breeze), a field of potential temperature
departures, 60(x,y) at z = H(x,y), were used. The
-------
Figure 3. Distribution of H(x,y) (m) field for the St. Louis,
Missouri region. See text for definition of H(x,y)
15
-------
—I
r!
\
Figure 4. Distribution of the 60(x,y) (°K) field for the average
heat island in St. Louis, Missouri. See text for
definition of 69(x,y).
16
-------
A6(x,y) = 69(x,y) [1-exp (-
where
Y is the rate constant and for this study was such that A9(x,y)«69(x,y)
in three hours
The boundary condition at z = H(x,y) on the potential temperature was
9(x,y,H) = ~9"(y,0) + A9(x,y)
If an obstacle was higher than one or more grid levels, then at z = H(x,y) and
I
at the side walls of the obstacle (which was not the case in St. Louis):
6(x,y,H) = ¥(y,H) + A9(x,y)
The initial synoptic field represented by the overbar is determined as if the
terrain were not present.
The heat and momentum flux at z = H(x,y) were determined through simul-
taneous solution of the boundary layer profile equations for the friction
velocities, UA and v^, and the scale temperature, T*. These parameters are
proportional to the momentum and heat fluxes, respectively, at the surface.
The profile equations used were:
u =
U
k IJ-"V ' 9 u*
,1/2
iy , ijj. nig ^z—n j J
' 8 U*
O
1/2
» Zl ri rz"Hv -»- aT*fcCg(z-H)}J
= 1. l-"l\_ ' ' Q TI*
V = k~ liTl(~) + 8U* ] ' and
17
-------
where U* = -J u£ + v£ , and a = 18. The above equations are adaptations of
those used by Estoque and Bhumralkar (1969). The roughness length, ZQ , was
allowed to vary in the following manner. In regions other than the built-up
sections of the city
ZQ = 0.0015 H(x,y)
Using the above expression and for the given values of H(x,y) in the St. Louis
area other than in the built-up sections, 1 cm _<_ ZQ ^_ 13 cm. In the built-up
sections,
ZQ = .02 H(x,y)
which yielded values as large as 1 meter. In the model only one major
built-up section was designated. This region was bounded by Forest Park to
the west, the Mississippi River to the east, approximately St. Louis Avenue
to the north, and Russell Avenue to the south. This technique is an adapta-
tion of a technique suggested by Lettau (1969).
The formulas for the vertical eddy diffusion coefficients were
xz
du
dz
A
A
yz
dv
dz
AH,
9
^ p te
H
dV
dz
where
~ + v . The mixing lengths, £, were given by
18
-------
A - k(z-H) (l-o^ Ri)
1/4
Ri)
lM
and
k(z-H) (l+ct2Ri)
-1/2
it
(1+a Ri)
-1/2
where ctj = 18, 02 = 11,. and Ri is the Richardson number. The above equations
were derived from expressions given by Rossby and Montgomery (1935) and
Holzran (1943). The vertical profile of the vertical diffusion coefficients
was constrained to have a second-order variability with height, and the
coefficients were set at zero above 800 m.
The horizontal diffusion coefficients were assumed to be proportional to
the square root of mean grid spacing; i.e.,
A = b
x o
VAx +
-^
Ax,
VAy +
—V
where the subscripts, u and d define the upstream and downstream grid-spacing
of the grid point, respectively. The constant b depended on the initial
boundary layer stability. Though AX and A increase as Ax or Ay increase,
19
-------
the ratio A /Ax or A /Ay decreases markedly, which reduces the horizontal
x y
diffusion in the unequal grid spacing regions. This was done to compensate
for the dampening produced by the upstream differencing in that region.
3.1.3 Initial Conditions and Integration
The initial wind field was assumed to have a wind direction parallel
to the x axis and a speed which varied in z alone. An initial potential
temperature distribution was given at x = -72 km and y = -72 km. Since the
synoptic state (initial conditions) was in geostrophic and hydrostatic
equilibrium, the necessary initial pressure and potential temperature field
were derived through integration of the geostrophic and hydrostatic
equations.
In order to study air flow in the city for various wind directions, the
H(x,y) and S6(x,y) fields which define the city area were rotated to
effectively yield a new angle for the wind. In order to accomplish this, two
H(x,y) and 66(x,y) fields had to be developed. One field was used for the
northwest, southwest, southeast and northeast wind directions. The other for
north, west, south and east wind directions.
For the time integration, forward time differencing was combined with
the leap-frog method to control the divergence of the solution at even and
odd time steps. Forty-eight (48) different simulations were performed for
St. Louis; these included variations of synoptic wind direction and wind
speed and of boundary layer stability. Exact specifications of the initial
stability and wind field and other model parameters used in the simulation is
presented in Appendix A. The Appendix presents a paper published in the
Journal of Applied Meteorology.
3.2 Determination of a Wind Monitoring Network
3.2.1 General
The basic 'strategy for selecting a network of wind monitoring stations is
illustrated is Figure 5. The objectives of this procedure are twofold:
1) to determine a statistical model form which could be used
to characterize the horizontal wind field over the St. Louis
area for a variety of meteorological conditions, and
2) to determine a network of sites for monitoring winds.
20
-------
CLASS OF
STATISTICAL
MODELS
SITE SELECTION
CRITERION
SITE SELECTION ALGORITHM
Figure 5. Design selection strategy.
21
-------
Once such a network is selected and set up, actual wind data would be
collected at each site in order to estimate the parameters of the statistical
model. This fitted model would produce an estimated wind field over the
region which would be combined with the emissions source inventory (in the
OVAM) to produce an estimated pollutant concentration field.
As indicated in Figure 5, the three basic ingredients required for
selecting a wind-monitoring network are
(a) a statistical model form (which relates winds to geographic
location),
(b) a set of M candidate sites (from which a subset of n is to
be selected as a wind monitoring network). and
(c) an objective criterion for comparing alternative networks.
Items (a), (b) and (c) are discussed, respectively, in Subsections 3.2.2,
3.2.3, and 3.2.4.
Given the set of M candidate sites, there are
C)
M!
n! (M-n)!
possible networks which should be compared via the chosen criterion. In order
to make any reasonable claim of optimality, it is clear that M should be much
larger than n (i.e., the class of possible networks is large). On the other
hand, the evaluation all possible networks in such a situation is
computationally unfeasible. (For instance, if M is only 50, there are
47,129,212,243,960 possible networks containing 20 stations which can be
selected. Ideally, of course, M would be much larger). Hence, it was
necessary to develop an algorithm which evaluates only some of the possible
networks. This algorithm is described in Subsection 3.2.5. Results of its
application are given in Section 3.2.6. In general, application of this
algorithm would yield a set of sites for monitoring winds which would permit
good wind field estimation. Because of several economic and practical
constraints, it was necessary to make several modifications to the strategy
depicted in Figure 5. First, it was necessary to apply the algorithm to two
different sets of candidate sites. Secondly, it was necessary to combine the
two resulting networks into a single network.
22
-------
Finally, for validation purposes, it was necessary to make some modification
to this network in order to take advantage of data available from existing
stations. These modifications to the general strategy are described in
Section 3.2.7 along with the resulting network.
3.2.2 Development of Statistical Model Form
In general, the development of the statistical model form involved two
basic steps:
(1) specification of a class of models which, for any given set of
meteorological conditions (covered by the simulated data, since
these are the available data), is assumed to contain a model
which will permit an adequate approximation of the wind field,
(2) determination of a single model form from within this class
which will provide a reasonable approximation of the wind
fields under all such conditions.
Selection of a statistical model form was based on data generated by the
primitive equation model. These data consisted of quasi-steady-state values
for the west-east (WE) and south-north (SN) wind speed components at various
geographical points (grid points). Forty-eight (48) sets of initial conditions
were used to create 48 data sets as indicated by the cells in the table below:
Wind Initial Prevailing Wind Direction
Shear Wind Stability
Code Speed Condition* 0° 45° 90° 135° 180° 225° 270° 315°
1 1.5 1
2 3.0 1
3 6.0 1
4 1.5 2
5 3.0 2
6 6.0 2
*1 = stable boundary layer (nighttime)
2 = neutral boundary layer (daytime)
23
-------
For each cell in the 0°, 90°, 180°, and 270° columns of the table, wind
speed data were available at each of the points shown in Figure 2. For the
other prevailing wind directions, the data were available for the points
obtained by rotating the grid points in Figure 2 by 45°. These sets of points
will be referred to as square grids and diamond grids, respectively. Primary
interest was in the central portions of the areas defined by these grids.
Hence, the model fitting was based only on the data at the 521 grid points
within:
(a) the square: {(x,y) : - 16 £ x _< + 16, - 16 _< y _< + 16} , or
(b) the diamond: {(x,y) : |x| + |y| <_ 16 /2 }
Two classes of statistical models were considered. The first of these
permitted models of the form
Zk(x,y)
where
Z, (x,y) = the kth wind speed component (k=l for WE and k=2 for SN)
at the point (x,y)
B, . = coefficients to be estimated
ay
P.(x,y) = a term of form x y , where a and j can take on non-negative
integral values.
This first class of models consisted of all terms P.(x,y) such that
0 <_ cc + Y £ 9
and
0 <_ a <_ 6
0 j< Y J; 6.
These 43 terms are represented by the products of row and column headings in
the following array:
24
-------
1
y
y2
y3
y4
y5
y6
1 X X2 X3 v* *5 V6
A. A X X A X
•••••••
•••••••
• ••»«••
• ••»•••
• •*•••
* • • • •
• • • •
Q (x,y, H(x,y)J
•J J
(2)
Since there are 43 such terms, this class contains 2 ^ subsets of terms (i.e.,
2^3 candidate model forms).
The second class of model forms differs from the first in that it permits
models involving the topographic evaluation. In particular, this class
contains models of the following form:
Zk(x,y) -'
where the Z, are as in (1), the C, . are parameters to be estimated, H(x,y)
denotes the elevation (above a base plane) at the point (x,y) , and the Q_.
represent terms of the form
H6 xa yY
The exponents are restricted to be non-negative integers in the ranges
0 <_ 6 <_ 2
0 <_ a <_ 6
0 — Y — &
In addition, the following restrictions were employed:
then a = 6 = 0
25
-------
Thus, the second class of model forms is represented by the following
arrays of terms:
f —C\ 1 V V^ V-* V^ V^ ^^O
1
y
y
6=1
i
y
y2
y3
y4
6=2
1
y
•
1 X XX X
• • • • •
• • • •
• • •
• •
•
1 X X2 X3
• • • •
• • •
There are, thus, 2 ^ potential model forms in this second class, some of which
are also in the first class.
In order to determine a "good" model for a particular simulated case and a
particular wind component, stepwise regression was applied twice—once for each
of the two classes of terms. This procedure allowed a determination of which
of the B, . and C. . could reasonably be assigned a zero value (i.e., determine
KJ KJ
which terms in (1) and (2) should be included). In all, 192 stepwise regres-
sions were performed (2 components x 48 sets of conditions x 2 classes of
terms). The Statistical Analysis System (SAS) (Service, 1972) was used to
perform these analyses.
26
-------
As might be expected, different models appeared "best" for the different
situations and components. In order to arrive at a single model form (i.e., a
single set of terms) for each component, the results of the 192 stepwise
regression analyses were tabulated in terms of how frequently the various
terms were picked by the procedure as being of significance. The lower order
terms, as might be expected, were judged significant in virtually all of the
cases whereas, the higher order terms occurred only rarely. This analysis
led to a set of 49 potential model forms which were judged to warrant further
consideration (see Table 1).
Each of these 49 models was fitted to the wind component data for each
set of conditions. For each model form and each case, the following summary
statistics were computed in order to determine the adequacy of the least
squares fit:
1) si, = residual mean square for the kth component (k = 1 for WE,
k = 2 for SN when the mth model (m = 1, 2, ..., 49) is
i
fitted to case i (i = 1, 2, ... 48)).
2) R~, = the proportion of the total variation in component k
accounted for when model m is used for the ith case.
2 _ 2 2
~ "ilm Si2m
3) s
977
Thus, for the mth model form, 48 values of sf-i , Si2m* si-m5
R? _ were determined. The means of these 48 values are shown in Table 2.
J. » Ul
Table 3 shows, for selected cases and models, the correlations between the
simulated data and the predicted values.
Selection of a model form from among these 49 candidates was, to a large
degree, subjective. One would prefer one of the simpler (i.e., smaller)
models in that fewer observations would generally be necessary for estimating
the model parameters (i.e., smaller networks would be possible). On the other
hand, one must be assured that all major variations in the wind field—within
the region of interest—can be characterized by such a model. Results such as
those shown in Tables 2 and 3 indicated that the 13-term model might satis-
factorily achieve these two conflicting goals. This was verified by visual
comparisons of several simulated wind fields and the corresponding wind fields
predicted via this model (Model #2 of Table 1).
27
-------
TABLE 1. TERMS INCLUDED IN 49 CANDIDATE MODELS —
—' An x indicates that the term, or set of terms, In the column heading is
included in the model.
4
* 7
28
-------
TABLE 2. MEANS OF R AND S OVER 48 CASES FOR 49 MODELS
No.
Terms
7
U
14
14
15
15
'~ts"~
16
16
16
16
16
16
17
17
17
17
17
17
17
18
18
18
18
18
18
18
18
19
_J9
19
19
19
19
20
"20-
20
20
20
20
21
21
21
21
22
22
23
Model
No.
1
2
3
4
5
6
7
8
9
10
11
12
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
"39
40
41
42
43
44
45
46
47
48
49
2
.506690
,612684
.610855
,619309
.621956
.626621
.-.628498
.632617
.631938
.628740
.637016
.632402
.633196
.639205
.639641
.649766
.641896
.641776
S647259
A649483
.660507
.643447
.659748
.654901
.644181
j.648129
.652812
.654352
.667160
.667651
.655282
.654789
.662321
.663166
.673714
.676539
.671257
^65888 7
.670341
,668447
.6791-20
.682313
.683092
.674997
.690203
.687850
.686697
.693491
.702166
Ri2m
_i505472
1.6 0 1 4 7 0
.602028
'.'608337
.610881
1620353
i627943
'.625636
.627727
I630877
'.619352
1625615
.61L836
1632162
.619730
^629442
.632942
^643474
.644662
'.651902
.643192
.644561
.646289
'.648201
.645746
.649922
'.648361
;655630
,649763
'.658291
,654599
^,655356
.655660
J659408
U68854
.663122
"".659430
.663070
1662787
1672308
.673161
"".678870
"".671031
.681627
1677992
.686474
1692708
2
3ilm
_.ft32?J70
^0250040
.0243681
.0247029
.'0237868
*fl245! 17
..0241404
•J123.47L9
.0234377
.0226128
.0242622
.0238663
.(J231702
.0228322
.0219225
.0236405
.0236537
.0229790
.0225609
.0215696
.0227069
.0213931
.0216514
.0225244
..0233922
.C227316
.6222756
.0212449
.0209566
.0219^46
.0229084
.C220989
.0220731
.0210597
.0206323
.0208271
.0217872
.0207275
.0218383
.0207915
.0203177
.0200279
.02J3450
.0200606
.0201613
.0198961
.0199033
.0193341
812m
.0273231
.0226008
.0228007
,0222233
,0224006
.0317743
,0208269
tP2}4662
,0211935
,0207780
,0220495
.0214110
.0206676
.0211175
.0203724
,0216140
,0209432
.0202812
,0206673
.0107121
.0206260
.0200606
.0204004
.0200157
,0200003
.0198670
.0204065
,0195608
.0202740
.0195745
.0195914
.0194937
.0199057
,0191196
,0199041
.0189649
,0192599
,0192769
.0192319
,0187527
,0196717
.0188329
.0187873
.0183824
.0191709
.0184982
,0184696
.0181627
.0177598
4,
1506081
j607077
^606441
jtlll"
^628220
>2"9227
'629633
,629808
.628184
J'629009
6326 16
^635693
C637419
>42625
.645960
7650693
7651649
1644004
^653019
1644963
(649025
.650?8fc
1654991
*65«472
.654940
7655072"
.659100
*672696
,667190
j659l59
.666706
^670939
.670953
^677310
.676126
7676934
7680617
1'684739
7682344
7689983
.697437
2
.'06G2<)01
.'0476048
.0471866
^0461874
.0462860
,0449672
^0449381
.04(44566
.0442157
.0446622
.'0456732
.0445339
.'0442877
.'0432046
.0435366
,0445837
.0439349
.0436463
.042273C
.0422156
.0427675
.'0417935
.'0416672
,0425247
.0432592
.043138C
.0418363
.0415186
.0405312-
.0415C6C
.0424021
.0420046
.04H927
J0409636
.0395972
.0400871
.0410641
.0399594
.0405910
.'04QA632
.'0391506
.0368152
.0397273
.0392316
.0386595
."0383659
.'038066C
.0370939
29
-------
TABLE 3. CORRELATIONS BETWEEN SIMULATED AND PREDICTED VALUES
FOR SELECTED MODELS!/
Wind
Shear
Code Angle
1 0
90
180
270
2 0
90
180
270
3 0
90
180
270
4 0
90
180
270
5 0
90
180
270
6 0
90
180
270
Model Number
#1
WE
.69
.67
.47
.67
.68
.74
.47
.56
.62
.75
.51
.55
.79
.76
.90
.60
.89
.75
.65
.83
.80
.78
.62
.76
(7)
SN
.59
.53
.61
.58
.56
.60
.75
.66
.75
.56
.75
.57
.81
.74
.74
.70
.86
.72
.76
.41
.89
.76
.78
.88
#2
WE
.80
.78
.58
.81
.77
.85
.56
.73
.72
.78
.56
.65
.84
.80
.92
.68
.90
.78
.76
.86
.85
.80
.69
.80
(No. of Terms Shown in
(13)
SN
.73
.69
.72
.67
.69
.73
.82
.73
.78
.64
.77
.64
.84
.83
.78
.76
.88
.80
.79
.65
.91
.78
.82
.89
#14
WE
.82
.82
.67
.83
.79
.87
.64
.75
.73
.79
.64
.68
.84
.81
.92
.68
.91
.78
.82
.86
.86
.81
.70
.80
(16)
SN
.76
.77
.77
.75
.72
.75
.83
.76
.79
.68
.78
.67
.86
.83
.80
.77
.89
.82
.80
.66
.92
.79
.83
.89
#36
WE
.87
.82
.71
.86
.84
.89
.69
.80
.77
.80
.68
.74
.85
.84
.92
.71
.92
.81
.82
.88
.87
.82
.72
.82
Parenthesis')
(19)
SN
.80
.80
.77
.77
.78
.78
.83
.80
.84
.71
.80
.70
.87
.84
.82
.77
.91
.83
.82
.69
.93
.79
.86
.90
#48
WE
.87
.82
.71
.88
.85
.89
.70
.81
.77
.80
.70
.75
.87
.84
.93
.72
.92
.81
.85
.89
.88
.82
.73
.82
(22)
SN
.81
.80
.77
.77
.79
.79
.83
.80
..84
.71
.80
.71
.87
.84
.83
.78
.91
.84
.83
.69
.94
.79
.86
.90
— Estimation of the model parameters was based on the simulated data at
the 521 points of the square grid [+16, +16], The correlations
were computed using the 521 pairs of simulated and predicted wind
components.
30
-------
Thus, the form of the prediction model chosen for input to the site
selection algorithm was
Zk(x,y) = bQ(k) + b1(k)x + b2(k)y + b (k)H(x,y)
(3)
b4(k)x2 + b5(k)y2 + b6(k)xy + b?(k)x3
bg(k)y3 + b9(k)x2y + b1Q(k)xy2
bn(k)x4 + b12(k)y4,
where k = 1 for the WE component and k = 2 for the SN component. It is
assumed that the b^(k) would be estimated, for a given case, by ordinary
(least squares) regression.
3.2.3 Selection of Candidate Sites
A second required input for the network selection algorithm, in addition
to the model form, is a set of candidate sites or candidate design points.
The choice of a large number of such points is obviously preferable, in that
a large number of potential networks covering a wide range of possibilities
should be included in the evaluation.
The set of points selected should encompass, and extend somewhat beyond,
the geographic area for which good estimation of wind fields is desired. On
the other hand, all of the points should fall within an area over which the
model can be expected to hold. This is particularly important, since the
design selection criteria (see Section 3.2.4) are aimed at minimization of
variance as opposed to bias. Such criteria not only assume the correctness
of the model but also tend to select some design points at the extremities
of the allowable region.
J It should also be noted that many practical restrictions can, in general,
be incorporated into the algorithm via the specification of candidate sites.
For instance, one can force certain points to be in the design (e.g.,
31
-------
previously existing stations); similarly, certain areas not amenable to the
establishment of wind/pollution monitoring can be taken into account by simply
not including any candidate points within such areas. In the context of this
study, a natural restriction on the set of candidate sites was introduced by
the fact that the model includes the elevation H(x,y). This dictated that one
know the elevation of any candidate site. For this reason, any candidate site
was required to be coincident with a point from either the diamond or square
grid.
As previously indicated, the model fitting was judged on the basis of two
regions (actually! on a discrete set of points within each of two regions).
These were (the 521 grid points in)
{ (x,y) : |x| l 16 and |y| <_ 16 }
and (4)
{ (x,y) : |x| + |y| <_ 16 /!"} .
As indicated above, the set of potential design points should, therefore, fall
within a subregion of these regions. This subregion should also cover the
urban area. One option for the set of candidate sites would be the union of
the grid points within the following two sets:
{ (x,y) : |x| <_ 10 and |y| £ 10 } (5)
{ (x,y) : |x| + |y| <. 10 } . (6)
It should be noted that each of these regions is contained in the intersection
of the two regions given by (4). There were 441 points in each of these
regions; since there was little exact overlap of the grid points in these two
regions, a large number of points are contained in this union. This number of
points not only would be computationally impractical in the design selection
algorithm, but also would include many essentially redundant points because
32
-------
of near-overlap. A second approach would, therefore, be to select a subset of
this union—for example, removing points which nearly overlapped other points.
Both of these approaches, however, involve a problem which steins from
the fact that the set of candidate sites would be made up of points from both
the square and the diamond grids. Hence, a design selected via the algorithm
would typically involve points which originated from the two grids.
Consequently, evaluation of such a design with respect to the simulated data
would be difficult since simulated data would be unavailable for some of the
design points. For instance, suppose that a 20-point design were obtained and
that 10 of the points originated from the square grid and the remaining 10,
from the diamond grid. Suppose one wished to fit model (3) to the
simulated data for a particular case (say, a north wind case) using these
20 sites. Since simulated data for a north wind case would be available only
at the 10 points originating from the square grid, one would not be able to
perform this evaluation (or would need to use some type of interpolation to
produce data at the remaining 10 points in order to do it).
In order to avoid these types of problem, it was determined that:
1) the design selection algorithm would be applied twice using
two different sets of candidate sites (namely, grid points
in the region (5) and those in region (6)),
2) the two resultant networks would be combined into a single
design by clustering together points which were close
together, and
3) the points of the single design would be forced to coin-
cide with points of one of the original grids (so that
evaluation of half of the 48 simulated cases could be
carried out).
It should be noted that modifications are currently being made to the
simulation program which will permit it to generate wind data over a single
grid for all cases. This will eliminate many of the somewhat subjective or
arbitrary decisions which were required to circumvent problems arising from
the use of two grids.
33
-------
3.2.4 Design Selection Criteria
In addition to the model form (3) and the above-described sets of
candidate sites from which potential sampling networks (i.e., designs) can be
chosen, the design selection algorithm requires a statistical criterion for
objectively evaluating various designs. The objective of the algorithm is,
thus, to optimize the given criterion by choice of design. Such design
optimality criteria are generally concerned with one (or both) of two sources
of error: bias and variance. In the present case, two sources of bias are
actually involved: (1) the failure of (3) to represent the simulated
data; and (2) the failure of the simulation model to portray reality. The
model selection strategy was directed toward choosing a model form for which
the first of these bias components is sufficiently small. Most practical
design selection strategies are aimed toward the variance source of error
(i.e., variation about the fitted model) since assessing bias requires know-
ledge of the correct model. Some protection against large bias is necessary
and is achieved indirectly by restriction of the region of interest to a
smaller area that that covered by the initial grid; unless the selection of
design points are constrained in such a manner, variance criteria will
typically yield unrealistic designs in that some points will be selected at
the extremities of the allowable region. The model (3) , or any similar
model, will obviously not be appropriate over too large a region.
In order to discuss design selection criteria, it is convenient to
rewrite the fitted model (3) in vector notation:
Zk(x,y) = _x' _bfc (7)
where x/ = (1, x, y, H, x2, y2, xy, x3, y3, x2y. xy2, x\ y4) ,
bk' = (b0(k), bjtk), .... b12(k)), and
k = 1 for the WE component, k = 2 for the SN component.
The vector b1. is the least squares estimate of the corresponding parameter
^"~ 1C
vector
(60(k)' 8l(k)' ••" Bl2(k))
34
-------
i.e., the underlying model for the observed wind component data is assumed
to be
Z(x,y) =
(9)
where the deviations e, from the surface Jl'jLi are assumed to be uncorr elated
K ^~ ic
. • 7
and randomly distributed with zero mean, and with variance a£, which does not
K.
depend on the point (x,y) .
A design of size n consists of a set of n points (x,y). The
coordinates, elevation, and observed value, respectively, for the ith such
point are denoted as (x.,y.), H. = H(x.,y.) and Z, . = A, (x.,y.). Corres-
ponding to each such point is the vector:
=
Xy
ii'
Xiyi' Xiyi'
The following notation is made:
(D) -
"
k2
where the subscript (or superscript) D is used to indicate that a matrix (or
vector) depends on the design Oi.e., the particular set of points involved).
The estimates b.
(D)
are determined as:
(D)
r1 x' z (D)
; XD -5-k
35
-------
At an arbitrary point (x,y), then, the predicted value for the kth wind
speed component, based on the particular design, D, is
z«V7) = ,• bk .
The variance of Zfe (x,y). denoted by vDCk,x,yl is given by
Var [Zk(D) (x,y)] E vD(k,x,y) = x1 (X^)~1 x a^ .
For the given model form and a particular design, the variance of the
predicted value at the point (x,y) is, thus, proportional to the variance
function s~(•) evaluated at the point (x,y):
sD(x,y) =jc' (Xjj Xjj)'1 x ;
that is, vD(k,x,y) = a^ sD(x,y).
Various functions of STN(*) can be used for evaluating designs of size n.
For example, one may regard D as "better than" design D* if
(a) sD(x0,y0) <_ sD*(x0,y0) for a particular point (xQ,yo) ,
(b) max sD(x,y) <_ max s,^(x,y) where R is a specified continuous
region of the x-y plane O£ a specified set of discrete points,
or
(c) average sD(x,y) _< average s ^(x,y).
The criterion chosen was similar to (c) above. Note that the criterion (c)
depends on
36
-------
(i) a class of n-point designs (containing at least two designs),
(ii) a set R (either continuous or discrete)— , and
(iii) the model form.
In the present case, however, the availability of the simulated data
provided some additional information. If certain areas within the St. Louis
region consistently tended to experience lower (or higher) wind speeds than
other areas, regardless of the initial set of conditions, then the fitted
surfaces should reflect this pattern. Unless one or more design points fall
within such an area, however, then it would be unlikely that the fitted
surfaces would behave in this manner. Because such patterns were expected,
it appeared that a criterion \diich incorporated this additional information
should be developed for use in the selection of designs. Accordingly, the
following criterion was chosen for evaluating designs of size n:
Design D is better than design D* if
average w(x,y) sD(x,y) ^_ average w(x,y) sDA(x,y) . (10)
(x,y)eR (x,y)eR
The weight w(x,y) attached to a point (x,y) was chosen to reflect patterns
which were consistent over varying conditions, i.e., if, at a particular point
(x ,y ), the wind speed was consistently higher and/or lower than was typical
for most points, then the weight w(xo,yo) would be large relative to the
weights for other points. Intuitively, the form of (10) suggests that
a "good" design will attempt to make Sjj(x,y) small when w(x,y) is large in
order to yield a small average value. The construction of the weights and
design selection criterion is described below for the square grid. A
completely analagous development holds for the diamond grid.
Twenty-four cases were available (6 wind shear codes x 4 prevailing wind
directions) which provided values of the wind speed components for each of the
521 grid points contained in the region {(x,y) : |x| <_ 16, |y| <_ 16} . This
set of points is denoted by R. Let Z, .(x,y) denote the observed (simulated)
I/
The use of a continuous region R, while somewhat more appealing, was
regarded as infeasible because the elevation H was not available on a conti-
nuous basis (but rather, only at the grid points).
37
-------
value for the kth component (k = 1 for WE, k = 2 for SN) at the point (x,y)
for the jth case (j = 1,2,3 ..., 24). The wind speeds were then calculated;
(x,y)eR
The mean and sum of squares were then computed for each case.
Zi = 5ll Z (x,y)
J R J
ssi = Z (Z.(x,y) - Z.)2
J R J J
j = 1,2,...,24.
The individual deviations from the mean, (Z.(x,y) - Z.), were then squared
and standardized by SS^:
(Z,(x,y) - ZJ
_ T ,2 (x,y)eR
Yx'y) = ss
j = 1,2,...,24.
Finally, the weights were determined by averaging the w (x,y) over conditions:
, 24
w(x,y) = -TT 2^ w (x,y) (x,y)eR.
J'=l J
The criterion for selecting designs of size n is, thus, minimization of T-Q
with respect to designs, where
521
R
(ID
38
-------
Thus, the function of the weights is essentially to reduce, in a probabilistic
sense, the size of the class of models of the form (3). Obviously, particular
models of this form can yield wind component surfaces (over the x-y plane)
with several hills, valleys, ridges, etc. The complete class of models of
this form would allow such patterns to occur in any part of the region. Many
of these particular models (i.e., specific j^ vectors in (9))may never occur—
e.g., regardless of wind conditions, the surfaces may tend to be flat over
subareas of the region of interest. Obviously, several design points in such
a subregion would be unnecessary. On the otherhand, in subregions where
hills, valleys, or ridges are common, the presence of several design points is
essential in order to characterize the surfaces.
The above-described strategy was repeated for the diamond-grid situation.
Criterion (11) with w(x,y) = 1 (i.e., criterion (c), as previously described)
was also examined.
3.2.5 Design Selection Algorithm
As previously indicated, the potential designs points (i.e., candidate
sampling sites) were restricted to a smaller region than the initial grid in
order to achieve protection against bias. In particular, for the square-grid
situation, the set of candidate sites was restricted to the set of grid points
contained in the region {(x,y) : |x| _< 10, |y| <_ 10}. This set contains 441
points (see Figure 2).
In the process of site selection, if one wished to consider all possible
designs containing 19 stations, one would need to evaluate the criterion 9.737
x 10 ^ times in order to obtain the best design. For designs of other sizes,
another enormous number of criterion evaluations would be required. Since
s-(*)i and most other design bptimality criteria depend on the inverse of the
matrix XI XD, a tremendous number of matrix inversions would be required,
which would result in insufferably long computer runs. Hence, there was a need
to develop an algorithm which would yield a good design without evaluating all
possible designs of a given size. It was also desirable to program this
algorithm so that it would be adaptable to various choices of models, criteria
and potential points.JV
\J This flexibility is achieved by requiring the user to specify the model and
criterion (in user-supplied subroutines) and the set of candidate sites
(as input to the program).
39
-------
A "backward elimination" algorithm satisfying these 'objectives was
developed. Basically, the algorithm works in the following manner. Consider
M candidate sites (numbered 1 through M) and suppose a design of size n is
desired. Let Aj^ be the class of designs (Dj", D2~ ..., D"} con-
taining M-l points, where DJM~^ denotes the design of size M-l obtained from
the set of M points by omitting the ith point. Let r> ' be the value of the
criterion for this design._/ The programmed algorithm evaluates the criterion
for all M designs in the class, A.,_i » and determines that the minimum value
is achieved for the design. D, (M~1' , i.e., r.01""1^ = min {r, (M~1) ,
/vr "I \ . f*jr <* \ fc *^ -^*
T« , ... FM }. the kth point is then removed from, consideration
as a possible design point. Hence, the set of candidate sites now consists of
M-l points and there are (M-l) possible designs of size M-2 in the class
A., 2« The .designs in this class are obtained by omitting the ith point
(assuming the set of points has been renumbered from 1 to M-l). In the same
manner as above, the program evaluates the criterion for each design.
f\r p\
The minimum value F, , is determined, and the set of points in D,
becomes the set of candidate points for the next iteration; this process
continues until n points remain. It should be emphasized that such a tech-
nique cannot guarantee that the "best" design of a given size is found (where
"best" is defined in terms of the particular optimality criterion) . Although
the resultant design should be reasonably good, the term "optimal sampling
network" should not be interpreted in a rigorous mathematical sense.
The backward elimination procedure requires
(M-n)(M+n+l)/2
evaluations of the criterion. For instance, for M = 441 and n = 19, the
design selection criterion would be evaluated 97,271 times. Explicit inver-
sion of matrices is largely avoided by the use of recurrence relationships
which express the inverse of a matrix in terms of previously computed matrices,
jL/It is assumed in this discussion that optimization is achieved by minimi-
zing, as opposed to maximizing, the particular criterion.
40
-------
3.2.6 Application of the Algorithm
The design selection algorithm was first applied to the square grid using
the weighted average variance criterion (11) where the average was taken
over the set of 521 grid points in {(x,y) : | x| £_ 16, | y| <_ 16}. The
13-term model (3) was employed and the initial set of candidate sites was
restricted to the set of 441 grid points in the square {(x,y) : | x| <_ 10,
| y| _<_ 10). The resulting designs of sizes 13, 14, ..., 30 are specified below:
13-Point Design
X
-10
-10
-10
-7
-7
-7
0
0
6
7
10
10
10
-10
0
7
-7
7
10
-10
10
7
-10
-7
0
10
Number of
Points
14
15
16
17
18
19
20
21
22
X
7
0
-7
7
0
-10
0
7
10
Z
-7
6
0
0
-1
10
-7
7
7
Number of
Points
23
24
25
26
27
28
29
30
X
10
7
-7
-7
6
5
-10
-7
-10
10
6
-10
-7
1
-7
-6
The points in the 14-point design consist of those in the 13-point design
plus the point labeled "14". The 15-point design contains the points in the
14-point design plus point number "15", etc.
The algorithm was then applied to the diamond grid. The average was
taken over the rotated set of 521 grid points in the region {(x,y) : | x[ + | y|
_< 16 /2 } and the initial set of candidate sites consisted of the grid points
in {(x,y) : |x| + |y| <^ 10 /2~ } . The resulting design points were as follows:
41
-------
13-Point Design
X
-14.14
-4.95
0.00
-7.78
1.41
-9.19
2.12
7.07
7.07
0.71
0.00
8.49
14.14
2.
0.00
-9.19
-14.14
-0.71
-8.49
4.95
-2.12
-7.07
0.00
7.78
14.14
5.66
0.00
Number of
Points
14
15
16
17
18
19
20
21
22
X
-9.19
-4.95
5.66
7.78
2.12
0.00
9.90
2.12
-7.78
I.
-4.95
9.19
8.49
0.71
-0.71
7.07
-4.24
-7.78
0.71
i
Number of
Points
23
24
25
26
27
28
29
30
X
-2.12
0.71
1.41
-0.71
0.71
-13.44
-9.90
7.07
L
12.02
13.44
8.49
-2.12
-9.19
-0.71
2.83
-1.41
The two 30-point designs are shown in Figure 6. It should be noted that
many of the selected sites occur on the boundaries of the allowable region. As
previously indicated, this is typical of variance-minimization criteria. Con-
sequently, the choice of the shape and size of the region of interest is very
important since it will typically have significant influence on the resulting
design. It was for this reason that the design points were constrained to be
in the innermost grids, and it is assumed that interest in estimating wind
fields is confined to this region. It should also be noted that since, for a
given value of y, the model form is a quartic in x (apart from the H term), one
requires a minimum of five distinct levels (values) of x in a design in order
to be able to estimate the model parameters. Projection of the design points
arising from the square grid, for example, into the x axis shows a clear
pattern of five essentially distinct levels (at -10, -7, 0, 7 and 10 km). A
similar pattern holds for the y coordinate. That such patterns persist even
with the inclusion of the H term in the model form and the inclusion of weights
42
-------
KEY
Square grid design
Diamond grid design
+10
-10
+10
Figure 6. Thirty-point designs from square grid and diamond grid.
43
-------
in the design selection criterion suggests that neither of these was over-
whelmingly important in affecting the site selections. With regard to the
influence of the weights in (11), this was substantiated by applying the
i.,
algorithm with an unweighted criterion (i.e., w(x,y) =» 1 in (11)). The
chief influence of the weights, as might be expected, was a slight shift of
some points near the urban St. Louis area into or closer to this urban area.
It should be emphasized that the influence of topography and/or the
value of a weighted criterion could be much greater for other areas of the
country—for instance, in areas of more complex terrain. In such areas, one
would also expect that a more complex model form would be required (assuming
the same size region of interest) to characterize the more complex wind-
component surfaces that would generally occur in such an area.
3.2,7 Comparison of Alternative Wind-Monitoring Networks
To arrive at a single design consisting of less than twenty stations, as
dictated by economic considerations, the clusters of points evident in Figure 6.
were reduced to single points. In particular, four 19-point designs were
chosen for futher evaluation. In order to have simulation data available for
evaluating these networks, the points of these designs were forced to be
points contained in the full square grid.-i/ The four designs, designated as
AO, Al, A2 and A3,are indicated in Table 4. The designs AO, Al and A2 were the
result of subjective clusterings of points. Note that AO and Al differ by only
one point and that A2 tends to concentrate more design points in the urbanized
area. Design A3 was determined by applying an objective statistical technique
(hierarchial fusion) for clustering the sixty points shown in Figure 6.
BMD^statistical programs package contains the clustering routine employed. The
The 19 points listed under A3 in Table 4 represent the subset of square grid
points closest to the centroids of the resulting nineteen clusters determined
by the procedure.
For each of the four designs, the 13-term model was applied to WE and SN
wind speed components for the 24 sets of simulation data associated with the
I/The choice of the square grid, as opposed to the diamond grid, was arbitrary.
The rationale for forcing design points to coincide with points contained in a
single grid was described in Section 3.2.3.
—/Biomedieal Computer Programs, Health Science Computing Facility, Department
of Biomathematics, School of Medicine, University of California, Los Angeles,
March, 1971, (Supported by NIH, Special Research Resources Grant RR-3).
44
-------
TABLE 4. FOUR 19-POINT NETWORKS SELECTED FOR EVALUATION
AO
X
-16
16
0
0
-10
-10
-7
-7
-7
0
0
1
0
7
7
5
8
10
10
y
0
0
-16
16
-10
10
0
-7
7
-10
-7
-2
6
-7
0
6
6
-4
10
Al
X
-16
16
0
0
-10
-10
-7
-7
-7
0
0
1
0
7
7
5
10
10
10
y
0
0
-16
16
-10
10
0
-7
7
-10
-7
-2
6
-7
0
6
-10
-4
10
A2
X
-16
16
0
0
-10
-10
-7
-7
-7
0
6
1
1
7
5
5
10
10
10
y
0
0
-16
16
-10
10
0
_-j
7
-10
-2
-2
8
-7
1
6
-10
-4
10
A3
X
-16
16
0
0
-10
-10
-8
-8
-8
1
-6
1
0
7
7
-6
10
10
8
y
0
0
-16
16
-10
10
1
-6
6
-8
-10
-1
8
-8
0
10
-10
-6
8
TABLE 5. AVERAGE CORRELATION AND SUMS OF SQUARES
OF DEVIATIONS FOR EACH NETWORK
Network
AO
Al
A2
A3
Average
Correlation
Coefficient
.703
.705
.684
.711
Average
Sums of
Squares of
Deviations
13.7
13.2
15.1
14.2
45
-------
TABLE 6. SYMBOL KEY TO VARIANCE FUNCTION 128
— Given in Figures 7, 8, 9, 10 and 11.
46
-------
square grid (6 wind shear codes times 4 prevailing wind directions: 0°, 90°,
180°, 270°). Average correlations and sums of squares of deviations between
"observed" and predicted values (over the innermost grid) are given in Table 5,
The differences between the average correlation coefficients and between the
sum of squares of deviations were too small to make these results conclusive.
Comparison of the designs by use of the variance function sn(x,y),
however, revealed a clear preference for design A2. Plots of these functions
are given in Figures 7 through 10. Table 6 yields the key to the plots.
Figure 7 indicates that AO is likely to provide poor prediction in the south-
east corner of the innermost grid; design A3 suffers a similar deficiency at
the northeast corner (see Figure 10). Actual plots of wind vectors using these
networks showed that, indeed, very large errors were found in the respective
regions. For this reason, designs AO and A3 were dropped from further consi-
deration. The regions of large errors apparently did not produce smaller
correlations and larger sums of squares of deviations because these designs
permitted sufficiently good estimates to be obtained elsewhere on the grid.
A comparison of the variance function plots for designs Al and A2 indi-
cated a preference for A2; in particular, A2 appeared to provide for better
prediction in the urbanized area. The following frequency distributions of
the variance functions s^(x,y) and s^2(x»y) > for instance, can be derived
from Figures 8 and 9.
Range
of SD(«)
<_ 1/2
1/2 to 1
1 to 2
Region:
Al
219
205
17
(-10 < x < 10\
1 10 £y <_ 10/
"I
: A2
224
211
6
Region:
Al
80
54
2
/ 0 < x < 7\
\-8 <_ y <_ 8 /
A2
100
36
0
47
-------
._.v*Bmcr .rcNciicu
2o,oooo«ooo
00
12.90000090
o.oooooooo
'
It
"
-1,03000000
-.12,00000000
•20,00000000
•
7 d
6 4
2
« 1
i
2 1
1
I I
1
i 1
1
4 1
2
5 1
6 5
J
1 t
I .1
1 1
... -I..O
1 1
1 1
1 1
1 1
2 1
2 2
2 2
1 2
1 |
1 1
1
i
2... .
4
t
i
i
Q
0
t
1
1
1
1
2
1
1
0
1
0
1
1
4
-
1 .
1 1
J - -1
1 1
t 1
0 0
• 0 0
1 0
1 0
1 0
•1 0
1 1
1 0
2 0
| |
1 0
0 0
0 0
Q. 0 0
1 I 0
^
2
1
1 2
..2 „. i
1 2
l__l._
1 1
a o
0 0
0 fl
0 0
1 0
t ft '
1 0
1 1
0 0
0 1
ft )
..J ....
2
. 1
2 . Z
If P
2 '2 1
2 1 J
1 1 1
1 1 0
000
000
0 0 1
0 1 1
o o i
0 0 1
000
0 C 0
1 1 i
1. .! 1
.- -
1
.......
2
2 2 (
.-. 2.._.2._.1...._
1 1 1
_.J..... 1 1. .
1 1 1
» ° »
000
o o i
~ «'" 1 i
t i I
I 1 i
0 1 1
0 i 0
00 1
-o" r "i —
tut
._.!:...«...]
i
-••
? •
2
2 ;
1 i
1
9
1 (
I (
0 (
1
1
I
i
1
0
(I
1
I. .(
2
3
J
> 1
I_L_
i
1
0
> >
) 0
0""
0
—o
0
0
0
0
o
1
1
! .?. ..
*
2
i i
i. »
0 0
> 9
0 • 0
0 0
b o
i i
I t
"0. 0
0 0
0 0
0 1
o i
i |
i 2
2 2
3 3
4
*
!.
S « !•
23 5 I
! . . " ~1
02 '5 !
0
6 . i . . „
0 '-.
0| 2
1
1 _ — .. ..
li i !
1 ~.
0 1
0 1 J
1
12 1 - -
? y
2 i n 5
« i
4 « 5
b 0 7 .
!
T « '
I
•l^j&OOOOOOO -9T50000000 «3lSOOOOCOO 2.^0000000 e.'SOOOOOCO
Figure 7. Variance function for design AO.
14.SOOOGOGC
-------
_VAR1»NCE FilWCriOM
" "" "~
VO
24.00000000
12.10400040-
4.00000000
y
11.10000000
MZ.'tOOVOOOO
•
-20.1C1UOOIIO
.... . ._
6 S. .
S 1
2
4 1
J .
,
2 1
1
1 1
1
2 1
t
« . 1
2
S 1
6 S
-
.1.
1 1
.L. J
1
0
t
1
I
. ...I
2 t
2 2
2 t
1 2
.1. 1
1 I
1 1
1 0
1 1
1 1
1 1
2
..I. .
2__ L
1 1 1
.0-0- - t
o o t
00 0
00 0
1 0 0
1 0
1 0
t 0
1 0
- 1 . .1
t 0
2 0
2 t
1 0
0 0
1 0
4 00
0 00
1 0 0
1 1 0
1 1
—
i
2
J
1
1
1
I
0
0 0
0 0
0 0
0 9
0 0
0 0
00
» 0
o Q
1 0
0 0
6 6
0 0
0 I ~
4 0
0 1
1
8
. ....
t
2
2
1
t
I
0
0
0
0
0
0
0
0
0
0
0
0
i
0
1
I)
-
2
2
1
I
0
0
0
0
A
0
t
0
4
0
0
0
0
A
1
I
1
2 ..
1
t
1
0
0
0
0
0
0
a
t
0
o -
0
1
1
I
t
-
1
~
I'.'.-l 1 ."."
1 t 1
1 1 1
1 1 1
1 1 I
0 0 1
0 1 0
o I a
Q t L
0 4 1
040
0 1
0 1
0 't
0 I
0 1
046
0 4 t
0 1 1
1 0 1
1 1 0
1
• • -•
- •-
2 . . _
J I
1 2
J 1
1 1
1 1
1 1
I 0 0
000
0 1 0
1 1 0
1 1 0
10 0
i i 6
1 1 0
2 1 0
1 4 0
ooo
0 0 0
040
iTot
000
0 0 1
1 2
• -
»
— -•
/I
t 1 I
t t 1
0 '1 I
ooo
ooo
0 00
000
o •>
4 A
0 ~t'~
0 I
r ii —
090
000
ooo
040
•o .50000001 oj.'scoooooo ^.snoooooo it.Sbuounuo
Figure 8. Variance function for design Al.
-------
.. LtbLNJl-A S. 1. UBS .,. U.J_<; .LlUi. ,-LlL.
: 29.19900000
i
- 12.10300000
4.40000009
••1.10000090
.f 12.10000000
' '" " •21.'400oOO«b
t>. s
5 J 2
2 1
|
4 1 1
V
1 0
1
2 1 t
1
1 i
1
1 1 » . a ._
2 2
1 1
2 1 1
1 0
0
1
2 1
5 } 2
A 5
4
1
0
0
0
1
1
t
1
1
t
1
2
I
1
0
1
0
0
0
1
. 1
4
t
0
0
A
0
1
0
0
0
0
1
1
.VAH1
PL
0
g
1
j
0
1
1
1 .
1
1 _
1
1
I-.-
1
. 1
1
0
1
1
ANCE riiNCfi
or ..OP-*— va..
2
1 1 1
. l._4 1 —
l)M
t
1 1
1 I 1
1 .1 .L
III 211
0 t i 111
00| 110
000
000
000
1 0 0
. 1 . 0.. .0 .
1 0 0
0 1 0
. l._l 0 _..
o i o
0. 1 I....
0 0 1
0 0 t
001
0 1 I
1
2
o n o
000
000
000
0 0.. 0
040
000
a. e.._i
o n |
o.o.o
000
1 1 1
1 1 1
1 1 1
l-.l. J
L )
t
1
. t
1
1
0
0
o
0
0
0
0
A
0
..... 0
0
1
1
1
I 1
1 1
1 1
1 u
4 0
0 0
Q i
1 0
0 0
4 1
.1. 1
0
0.
0
1
1
1
1
2
t 1
III
1 1 1
1 1 1
0 1 1
0 1 t
000
000
u 0 0
• —
1
_L
0
ft
0
0
1
0
1
1
t
i
- 1 . . ...S 6
1
i s !j
- 1
III t
11 1
01 t X
& I)
« 0 1
0 0 J
00 I 2
0 1
411 j
C 0 0 0 .0 I : i
o o o o o i i i 1
000
1 0 0
000
MOO
i o o
1 0 0
I 0 0
0
0
_-»
0
o
0
0
0
t_l_l I
»
'I 1
0 0 1
.«..»._ 0 ;
0000 2 •
000 !
0001 1
0001 4 il
1 J 1 J II
2 j » s :
1
'156
• IS.&IIOOOOOO t9. 5000040
0
-i.'sooooooo
2
,'iOOOOOOO
0
j
,50000000 14.50000000
Figure 9. Variance function for design A2.
-------
VARIANCE FIINP. riuH_
'
20.40000000 t
12,10000000
1.100000JJO
.1
-'t. 10000000
-ll.'IOOOOOOO
••29.10146000
_
. ....<• .5
_ . ...» i
2
1 1
1
2 1
.._... 1
I 1
1
1 1
I
5 I
1
6 4
7 o
2
...I «
1 1
A 0
1 A
1 0
1 1
1 1
. .1 1
2 1
2 2
2 2
1 2
1 1
1 0
1 1
• ft- i
0 ft
t 0
1 1
1 1
2
4 _
2
0
0
0
0
0
0
0
1
2
1
2
1
2
2
0
0
1
1
0
0
0
t
4
I
000
o .0 L_
0 1
0 0
0 ft
6 0
0 1
1 0
1 0
I • o
2 0
1 0
2 0
1 |
0 0
0 0
0 0
I q
i i
1 0
1 0 0
t_
3
I
I i
1 .2
2 2
1
0
1
1
0
0
0 0
0 0
0 0
1 0
1 0
1 0
1 1
0" 0
0 0
0 1
0 1
0 t
(
2
L
2 2
._ 1.... I
2 2
2 2
1 1
1 1 .
0 1
1 o
I o
o a
ft 0
o i
I 2
d I
i i
0 0
0 " 0
1 1
" 1 i
1 I
1 1
-
™-^— *^—
1
t
1
1
1
1
1
1
1
1
0 0
o "o
0 1
i 1
1 1
2 1
2 2
1
. „ . {..
1
1
1
1 1
1
•
1
.
- •
1
t
.1 2
I 2
1 1
2 2
2 2
2 J
2 2
2 2
~Z If
0 1
"t I
1 I
1 1
1
..
-J - -
1 2
222 2
1 2 2
U 1 2
1 1 I -
10 »
000
1 1 1
2 2 1
2 2 I "
211 1
2 2' » 1
220 1
I 2 0 0
2 0' 0 " 0
U 0 0 0
r o o~- o
ft 0 0 0
t a o o
0000
1122
1 1
1
«.
1 4
2 2 J.
2 2
2 2
1 "2
1 1
1 1
1
1
I
1
' j
1
I
0 0 1
0 J 0
00" i>
000
40 0" '
1) 0 i
f. 2 1
1 i
'I
.
~
" 1
ft ?
.
S 6 !
4 .
i i
2
1
1 1 1
1 i
i 2 ; '
1 |
1
~r Hi
s s
i
5 6
,„...
•li.iOOOOOUO
t9.'50004000
•4.50000000
2.50001)000
8,5000001)0
H.SOOOO^OO
Figure 10. Variance function for design A3.
-------
The preference for A2 was further demonstrated by comparisons between pre-
dicted and simulated wind fields. The predicted wind fields were generated,
for several cases, by estimating the model parameters from simulated data at
the design points of these two networks. Hence, A2 was chosen as the best
19-point design. All of the evaluation criteria (average correlations, sums
of squares of deviations, variance function plots, and direct examination of
predicted versus observed wind fields) indicated that wind field prediction
outside the region [_+_ 10, _+_ 10] was highly unreliable (for any of the designs).
Figure 11 is a plot of network A2 which will, from this point on, be referred
to as the OSN (Optimum Sampling Network).
Implementation of this 19-point design would typically proceed on a
sequential basis; consequently, it may be important to know
1) the order in which stations should be established, and
2) the minimum number of stations which must furnish wind data (i.e.,
at what point in the sequential implementation can one begin to
utilize the data to estimate wind fields?)_/
Also, even after full implementation of the network, not all 19 stations would
always be able to furnish wind data. Hence, it was also important to have some
idea of the impact of missing data (non-reporting stations) on the wind field
estimation. The results below, though not specifically directed toward this
latter problem, do provide some insight into how many stations must be func-
tioning at a particular point in time in order for the network data to be
useful.
To obtain designs containing fewer than 19 points, the design
selection algorithm was applied once again using the same model and criterion
as before (i.e., Model (3) and Criterion (11). The initial set of
candidate sites for this run, however, were the 19 points of the OSN (i.e.,
M=19). The resultant designs can be specified as follows:
^/Estimation of the parameters of model (3), of course, requires a minimum
of 13 stations; the questions here are concerned with how many additional
stations are needed for reliable estimation and which stations should these
be.
52
-------
\
Figure 11. Stations in the A2 network (the OSN) for St. Louis,
Missouri. See text for discussion of numbered
stations•
53
-------
Design
No. Points
Specification
B
C
D
E
F
G
18
17
16
15
14
13
Remove (7,-7) from OSN
Remove (-7,7) from Design B
Remove (5,1) from Design C
Remove (5,-2) from Design D
Remove (-7,-7) from Design E
Remove (1,8) from Design F
The distributions for the variance functions for the above designs over 441
points in the innermost grid are characterized in the table below:
Design
G
F
E
D
C
B
OSN
No. of
Sampling
Stations
13
14
15
16
17
18
19
Frequency
-------
4
.m-4 ^SSDi:Jlm+SS1W- (12)
Values of this quantity are shown below:
Design
G
F
E
D
C
B
OSN
No.
Points
13
14
15
16
17
18
19
1
39.9
32.7
27.5
27.8
26.2
25.6
25.3
2
81.8
58.5
50.4
49.7
50.6
48.9
48.5
Wind Shear
3
113.9
85.2
63.2
63.3
65.0
63.5
64.7
Code
4
6.4
4.6
3.1
3.2
3.4
3.4
3.3
5
19.9
12.6
8.7
9.9
10.1
10.0
9.4
6
76.1
43.8
26.5
29.3
30.6
31.4
30.2
The above tables, and the results shown in Table 7 (which presents a summary
of the correlations between observed and predicted values for each design) are
quite consistent. They indicate that the 13- and!4-point designs are substan-
tially inferior to the larger networks. The 15-point design, however,
does not appear to be unduly bad with respect to the 19-point OSN.
Because it was economically unfeasible for RTI to establish a 19-station
field program for validating the results, it was necessary to make a further
modification of the basic strategy (Figure 5). This modification, which is
described further in Section 4.0, involved some shifting of the points in the
OSN (in order to be coincident with existing RAPS stations or city/county
stations). The resultant design points, referred to as the.operational OSN
and denoted by OSN*, are shown in Table 8 along with the corresponding points
of the OSN. The variance function plot for the OSN* is given in Figure 12.
The distribution of the variance function for the OSN (Figure 9) and the OSN*,
over two regions of interest, are given below:
55
-------
TABLE 7. SUMMARY OF CORRELATIONS BETWEEN OBSERVED AND PREDICTED VALUES
OVER 441 POINTS OF INNERMOST GRID—BY WIND SHEAR AND DESIGN
Design
G
F
E
D
C
B
OSN
G
F
E
D
C
B
OSN
G
F
E
D
C
B
OSN
G
F
E
D
C
B
OSN
G
F
E
D
C
B
OSN
G
F
E
D
C
B
OSN
No.
Points
13
14
15
16
17
18
19
13
14
15
16
17
18
19
13
14
15
16
17
18
19
13
14
15
16
17
18
19
13
14
15
16
17
18
19
13
14
15
16
17
18
19
Correlations for 4 angles x 2 components
Minimum Mean Maximum
.02
.44
.47
.48
.48
.50
.49
.24
.34
.29
.28
.29
.46
.46
.16
.31
.30
.29
.29
.38
.37
.24
.50
.53
.48
.48
.54
.58
.37
.44
.55
.58
.59
.59
.60
.54
,55
.59
.57
.57
.57
.59
.50
.63
.65
.65
.66
.67
.67
.46
.58
.61
.60
.60
.62
.62
.48
.51
.54
.53
.52
.54
.53
.65
.69
.72
.71
.71
.70
.72
.63
.67
.71
.70
.70
.70
.72
.66
.69
.73
.72
.71
.71
.72
.72
.79
.79
.79
.79
.79
.79
.78
.82
.83
.83
.83
.83
.83
.67
.68
.71
.70
.71
.71
.70
.90
.88
.89
.89
.89
.88
.88
.80
.81
.83
.79
.78
.80
.80
.82
.85
.88
.86
.85
.85
.85
56
-------
Region
r-10 ^ x _< 10) ( 0 _< x <_ 1\
: (-10 <_ y <_ 10/ Region: )-8 _< y <_ 8/
Range
of SD (.)
< 1/2
1/2 to 1
1 to 2
2 to 4
OSN
224
211
6
0
OSN*
223
205
11
2
OSN
100
36
0
0
OSN*
108
28
0
0
The above results (derived from Figures 9 and 12) clearly indicate that the
performance of the OS;N* relative to the OSN is substantially poorer over the
western portions of the innermost grid; over the urban area, however, the
two networks compare favorably.
TABLE 8. SPECIFICATION OF DESIGN POINTS FOR THE OSN AND THE OSN*
Number
OSN
OSN*
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
x
-16
16
0
0
-10
-10
- 7
0
1
5
10
10
10
1
- 7
6
5
- 7
7
1
0
0
-16
16
-10
10
0
-10
- 2
6
-10.'
- 4
10
8
- 7
- 2
1
7
- 7
- x
-20
18
5
0
- 9
-16
- 7
0
1
5
10
9
10
0
- 6
4
6
- 4
9
1
- 3
0
-16
16
- 8
8
2
-10
- 1
7
-10
- 2
12
10
- 6
- 3
1
8
- 6
City/ County 6
EPA 109
EPA 118
City/ County 8
EPA 119
EPA 120
City/ County 9
City/ County 2
EPA 106
EPA 102
EPA 104
EPA 108
EPA 113
EPA 105
EPA 101
EPA 110
57
-------
VARIANCE FUNCTION
10
10
Ul
00
• 10
.11
•id *
it..
•20
2
1
2
2
2
2
2
2
1
1
it
2
2
2
2
1
1
1
1
I
>
2
2
3
4
t
\
\
?
S i
1
1 1
1 1
1 1
0 0
0 0
A
1 1
0 A
0 0
0 I)
1 0
1 1
1 1
> 2 1
1 3
1 1
1 1
0 1
0 0
o n
0 0
i
0
1
0
i 1
1 1
0 0
0 0
0 0
0 0
0 0
1 1
1 1
2
ItlOOlOll 10 1
lOnoiiiii >o i
1 t A 0 0 0 0 1 1 00 I
001)000001 00 1
DOA000001 01 1
QOflOOOOOOOOl 1
OOA000000011 1
1A*AAAAA»QA01 I
QnQOUVUM/uvi i
JlrtOOOOOOOOO 1
OOAOOOOOOOOO 1
oonooooooooo i
OOA00000100001
oonooooooooooo
OOA 1001 1000000
OOAlOOOOOOOOOO
OOA| 1110000000
OOnOI 100000000
00111111000000
01111111000001
11(11111000011
1 1 1 1 1 1 1 1 I 1 1 1 1 >
2 1 11 222
• 15
• 10
10
15
20
Figure 12. Variance Function" for the OSN*.
-------
3.3 Objective Variational Analysis Model
The air pollution distribution is obtained using an Objective Variational
Analysis Model (OVAM). The model simultaneously minimizes the integrated error
variance of the observations from the analysis and the weighted error variance
of the analysis from a steady-state solution. In order for the integrated
variance to be a minimum, a second-order, Euler-Lagrange equation must be
satisfied under specific boundary conditions. The Euler-Lagrange equation is
solved by relaxation techniques. The model requires the wind field, the
observed pollution concentrations, and the source emissions. This model is
generalized so that these parameters can be specified for any urban area.
From an analysis point of view, the distribution of pollution concentra-
tions is primarily dependent upon 1) the distribution of pollutant sources,
2) the mean and turbulent properties of the airflow, 3) measurement locations,
and 4) the measured concentrations. The first two of these ingredients are
clearly evident in the continuous point source equation:
iKx,y,z) = —~ exp <- 2-(?~~)2 + (g^f (13)
2ira a U ( y z )
y z * '
(See Glossary for definition of terms). Equation (13) is a steady-state
solution to the mass continuity equation for the pollutant, for constant wind
and stability conditions. The distribution of the pollutant concentration
can be estimated when the wind, the stability, and the emissions are known.
A class of steady-state Gaussian dispersion models such as the Air Quality
Display Model (AQDM) (TRW Systems Group, 1969) or the Climatological
Dispersion Model (CDM) (Busse and Zimmerman, 1973) estimate the concentrations
at given locations under specified conditions. Even more simplified models
such as proposed by Gifford and Hanna (1973) are used to estimate concen-
trations downwind of uniform emission areas. In general, the estimated value
will disagree with a measured concentration for the same location because of
inappropriate assumptions, poor knowledge of input variables, or both.
Within an urban area, the pattern of pollutant concentrations is not
always apparent from a set of concentration measurements at fixed locations.
59
-------
The presence of sources, the limited number of sampling locations, and the
variable wind velocity mask the true distribution from those obtained by
conventional subjective or objective analysis.
In conventional objective analysis, the estimated value of a variable
at a given location is extrapolated/interpolated from observed values of the
variable at specified locations. Cressman (1959) and Barnes (1964) used
different distance-dependent weight functions to interpolated within a given
radius of an observation. Endlich and Mancuso (1968) used an elliptical
weight function oriented along the wind direction and proportional to the
wind speed. These techniques smooth the distributions and are incapable of
showing maxima or minima where there are no observations. In an urban area
with nonuniform sources, it is necessary to have an analysis which accounts
for variable distributions of sources and, possibly, locally high or low
concentrations where there are no observations.
Sasaki (1958, 1970a, b, c) proposed an analysis technique, first as an
initialization procedure for a numerical forecasting scheme, that uses the
observed variables at given locations and a model describing the behavior of
those variables to produce analysis which is consistent with both the obser-
vations and the model. The technique is based on minimizing the weighted
error variance of 1) the variable, , and 2) the
model equations, E(), over the domain of the integration. The 'total
weighted error variance, EV, is given by
EV
• / oU-402 + aE2 ((()) dfl •
The terms a and a are the weighting functions and are inversely propor-
tional to the typical variance of their respective terms. This formulation is
called the "weak constraint" by Sasaki since it is E2() which is being mini-
mized rather than E(4>). The distinction is significant since E(<|>) is normally
j Q
zero. For example, in some applications the adiabatic equation TT = 0 = E(<|>)
might be an appropriate constraint. In the weak constraint formalism, 8, the
potential temperature is not absolutely conserved over the domain, but its
variance while moving a parcel of air is conserved in accordance with the
observations.
60
-------
The distribution of , which minimizes the integral, EV, is given by the
solution to the Euler-Lagrange equations (see Lanczos, 1970). Those equations
specify the partial differential equation of which must be satisfied on
the interior of the domain, and specify the boundary conditions which must be
met. In most instances, the solution can be obtained by numerical integration.
Sasaki and Lewis (1970) used the technique to analyze three-dimensional
mesoscale distributions of wind components, temperature*and moisture asso-
ciated with severe storms. Stephens (1970) obtained synoptic scale distri-
butions of pressure height fields consistent with the balance equation. All
these users have shown that the solution technique acts as a low-pass filter
for real or for induced (simulated) errors in the observation field.
Wilkins (1971, 1972) applied the variational objective analysis princi-
ples to urban air pollution analyses. He described a method of adjusting a
pattern of concentration isopleths using the classical diffusion equation as
a constraint, whenever discrepancies in temporal continuity occur between two
successive analyses. Errors arising from misanalysis or from a variable wind
field were also considered. Using realistic estimates, Wilkins showed that the
analysis error of the sequences of an analysis could be reduced to one-tenth
the initial error. He also outlined the potential application of the method-
ology as a part of the continuous urban air monitoring surveillance program.
The approach used in this research project reverts to more basic forms.
The differential form of the mass continuity equation for a trace material is
used as a constraint instead of an integrated form (Eq. 13). The concentration
observations are used where and when they are available. Wilkins1 approach
could be used to adjust the resulting analyses.
The purpose of the objective analysis is to obtain an estimate of the
concentrations where observations are missing. This analysis approach
requires that the analysis fit the observed concentrations and the existing
emission, transport, and diffusion characteristics throughout the urban area.
The "best analysis" is defined as the distribution which minimizes the
weighted sum of the error variance a) between the observed and analyzed .
concentrations, and b) the departure of the results from a steady state.
61
-------
3.3.1 Description of the Model
Consider a volume having dimensions X (west-east) , Y (south-north) and
Z (ground-up) which encompasses a flat urban area. The observations of the
pollution concentration, i|i, are known for various locations. The ground-level
emissions density, Q(x,y) as well as the velocity components, u, v, w (in the
x, y and z directions, respectively) are known throughout the volume. The
equation of mass conservation of a passive (inert) gas in the atmosphere is
written:
<»>
where
q
ip is the concentration (ug/m ) ,
2
K^ is the eddy dif fusivity in the horizontal (m /s) , and
K is the eddy diffusivity in the vertical (m /s) .
Z
Since observations are normally made only near the ground, little is
actually known about the vertical distribution of the pollutant or about the
local vertical velocity, w, from measurements. Therefore, seeking a way to
express the concentration within a given volume in terms of its ground-^level
concentration, the following assumptions are made:
1) iKx,y,z) = 0(x,y,0) exp [ - -| (z/Sz)2]
2) w w 0 in the layer between z = 0 and z = Z
3) Z is small
4) u and v do not vary substantially between their measurement
height and Z
5) All emissions are from ground-based (z = 0) sources. (This is a
simplifying assumption, not a model requirement).
Integrating Eq. (14) from z = 0 to Z, gives
Z
9t 3x 3y
Z
~.. _ >, ~ , \LJ)
62
-------
where
and
L
I
\|> dz
2*
'
3y
By Assumption 2,
0
By Assumption 5,
K
z 3z
.00
0
Q(x,y)
Given the functional form of if%
/\
.
x,y) • / exp - j
•/Q
dz
or
(x,y) er
= Y
and
K
z = Z
= - K(Z)
i* u
Substituting Eqs. (16), (17), (18) and (19) into (15) gives
(16)
(17)
(18)
(19)
Vu2 i|> + K (Z) 3 • * - Q = 0. <20)
H -O z o
63
-------
Since reported concentrations are averaged over a given time increment
or are an instantaneous value during the period, the time variability of the
concentration over that interval cannot be accounted for. If the emissions
data, Q, and the wind data are appropriate for that time interval, a nearly
steady state (i.e.—La-0) may be assumed for the time interval. With that
3t
assumption, the emissions must be locally balanced by advection and
diffusion out of the volume. (The assumption of a steady state is not a
requirement of the analysis model. If only one set of observations is
available, it is the best assumption. If data for more than one time are
available, then the time variability can be included).
A dimensional analysis of the terms of the continuity equation shows the
relative contribution of each term and will also aid in computations by
reducing truncation errors. Defining characteristic magnitudes of each
quantity given
x,y = L • (x',y')
Q - Q • Qf
K = K • K' , and
u,v = U (u',v')
where
fy is a characteristic magnitude of fy , (gm/m^),
L is a characteristic length of the urban pollution system (m) ,
Q is a characteristic surface emission magnitude (gm/m^ s),
U is a characteristic wind speed (m/s) ,
and K is the characteristic magnitude of K^ or K (nr/s) ;
Equation (20) becomes
>:> + * B • *' - Q Q' - 0 . (21)
64
-------
Dividing Eq. (21) by —J-L, dropping the prime (') from the non-dimensional
Li
variables, gives
where
R = K/UL
n =
5 = KBL/yU.
The advection terms, u'-r-^, are on the order of unity and the other terms
oX
are scaled by the non-dimensional coefficients.
In the atmosphere, K, may take on a wide range of values, depending upon
stability conditions. In this study, the conditions are classified by the
Pasquill stability categories (Turner, 1969) and the value of the vertical
diffusion coefficient, a , for that category. The two terms are related
Z
through the Fickian solution to the steady-state diffusion equation which
gives
a 2 = 2Kt
z
Letting Z be the characteristic value of a and T = L/U be the characteristic
Z
travel time gives
K = I2 U/2L . (23)
Choosing Z~ 0" (x = 1 km), the minimum value of a was encountered. That value
Z Z
is related to the characteristic values of S and y> the coefficient of fy
and n, the coefficient of Q.
The magnitude of the coefficients determines the relevance of the
respective terms as compared to the advective term. The simple dimensional
terms like L, , U and Q can be defined. The terms £ and n are, however,
functions of the standard deviation of the local vertical profile of the
concentrations (S ) which is a function of the local source. Sz is
65
-------
dependent upon but not identical to a , the vertical standard deviation of a
z
point source of material. S reflects the upwind emissions, transport,
z
diffusion, and "background" so that
*(x,y,Z)/*(x,y,0) = exp I - f(Z/Sz)2J
or
sz = z/
At any location, the ty and i|> , will be dependent upon the vertical
diffusion of the atmosphere and upwind sources. I tends to maximize the
coefficients of $ and Q. This definition permits K to adjust with changes
in stability classification.
The 3 term can be written as
- \ (Z/Sz)2J .
BOO = T Z/S J exp - f (Z/Sy\ . (24)
2
The maximum value of 3 occurs for (Z/S ) 'v 2.0. If 3 is evaluated there
(^0.75/Z), the term in Eq. 21 representing vertical flux through the top of
the volume will be maximized, and other characteristic values can be evaluated.
Choosing Z = /2~ S gives a characteristic value, y> f°r Y»
Z
Y = 0.68 /ir/2 S ,
Z
or
Y = 0.68 /ir/4 Z
The characteristic value of £, ?, is given by
I ' K L
When (Z/S )2 ^ 2.0, C « 0.62 Z2/Z2. If Z2 ^ a 2 ^ S 2, then
Z Z 2
T«0. 31
The characteristic value of n» n, is given by
"n = "Q L/? UY
66
-------
Choosing
Q = 10-8 gm/m2 s
L = 4 km
* » 10-6 g/m3
U = 3 m/s
Y = 0.6 Z
and Z = 102m
gives n = 0.22
With those assumed values, and Eq. 22, R is given by
R = E2/2L2
« 1.56 x 10-4 .
The horizontal diffusion term, therefore, is negligible with respect to the
other terms in the analysis as Wilkins also showed.
3.3.2 Analysis Technique
The approach to the analysis techniques taken here was first implemented
by Sasaki (1957). Following the hypothesis that the best analysis is one
which minimizes the error variance of the observations and the continuity
equation through the domain (X by Y), the error variance integral, EV, is
formed, i.e.,
EV
/ / [a('JH') + oA2Jdydx (26)
(26.)
where ij; is the non-dimensional value of the observations,
<|) is the non-dimensional value of the analysis,
a is a weight function inversely proportional to the
variance of ^,
and a is a weight function inversely proportional to the
variance of A.
67
-------
The distribution of i|> which minimizes the integral is given as the
solution to the Euler-Lagrange equation (See Stephens, 1970). For the
integral in Eq. (26) that equation is
a OMO + aUA - - ) = 0 (27)
subject to the boundary conditions that ^ is invariant on the boundary
(6ip = 0) or the normal gradient of i|> is invariant ( = a(n |^ - nSQ) - a ^ = H(x) (30)
which is a one-dimensional Helmholtz equation of the form
- y2* - H(x) .
3x
The forcing function, H, is governed by the weighted observations, the
emissions, and the emissions gradient. The magnitude of a and a also influence
the solution. If a/a is large (-10), Eq. (28) reduces to fy = ty. If a/a is
small (~ 0.1), the observations have far less weight and the distribution of
sources governs the solution, e.g.,
68
-------
-n-neQ - (3D
3x 3x
This is exactly the analysis objective—use the observations where they are
known and use the emissions and meteorology where the observations are missing
to determine the distribution of fy over the area.
Since (a£ + a)/a > 0, Eq. (30) is an elliptic partial differential
equation with constant coefficients and is amenable to solution by numerical
relaxation. Similar numerical techniques should be valid for integrating
Eq. (26) to obtain a solution when the coefficients are variable.
3.3.3 Numerical Solutions
A finite difference grid (I x J) with uniform grid spacing in the x and
y directions is assumed to encompass the urban area. At each grid point
(i,j) the wind components (u,v) are known; the emission for the area bounded
by unit increases of i and j is assigned to that point and is assumed constant
over the unit area; the measured (observed) and estimated concentrations are
also assigned. The measurements do not necessarily occur at the location,
but we assume they have been interpolated to the grid point.
A first guess to the concentration distribution is made at z = 0 and at
z = Z, assuming that a steady-state condition with a constant wind, the
given emissions, and a representative atmospheric stability exists. A
background concentration, \1> , is added to give fy and 0.999 i|)_ is added at
B o o
z = Z to give fy , so that S will be defined at locations upwind of any
Z z
source. Although a CDM estimate at the two altitudes might be used, a more
rigorous integration has been chosen than the virtual point source approach
used in the CDM (see Appendix B). The estimates of fy at z = 0, and Z are
analyses used to estimate S , y, n and £ for each grid point.
Z
With these data, a solution to the analysis equation is obtained using
the point-by-point Richardson relaxation procedure (Thompson, 1961). Schemes
with more rapid convergence might be used, but may present some problems when
the coefficients are variable. If the Richardson procedure will not give a
solution, neither will any of the others.
69
-------
of
At the v th iteration, A. . is computed with existing grid-point estimates
\\ -J
, (fy . ) , using the finite difference analog of Eq. (26a), i.e.
where r = L/2Ax
and n=fO
4.0
is the mean value of the area emission adjacent to the point (i,j). At the
•
south and west edges of the grid, Q is estimated by linear extrapolation
outward from interior grid points. Along the boundaries, the normal gradient
of ip is computed using the current values of ty at the boundary and the first
interior row or column of the grid. It is advantageous to compute A since the
value of the integrand EV of Eq. (26) is evaluated at each iteration.
v
The residual, R^. , at each interior grid is given by the finite differ-
ence analog of the analysis equation, (Eq. 27), i.e.,
(vA)M-l] J
+ (*M - * l.J> '
The improved estimate of ijj , tjjv, , is given by
i» J ij
where
Successively applied, the recursion formulas lead to convergence of the
solution in the constant coefficient case. The solution is obtained when
,
' *
«
at every point in the interior of the grid.
< 0.01
(32)
70
-------
This criteria is quite stringent, indicating that the solution is known
to one percent or less at every point in the area. Realistically, even in
the best conditions, errors of concentration are 5- to 10-tImes larger.
In some instances, v = 90 also served as a cutoff value. The trend of the
solution is well established by that iteration and convergence with more
iterations is assured.
3.3.4 Weighting Factors
The weighting factors, like the other terms of the analysis equation are
non-dimensional. They are inversely proportioned to the variance of ty and -2i-
about the domain of interest. Their relative value (i.e., the ratio a/a) is
more important to the resulting analysis than is the absolute magnitude of
either factor.
The non-dimensional value of A at any point increases in proportion to
the ratio L/U (r\, C, and the advection derivative depend upon L/U) . The non-
dimensional variance of A is, therefore, proportional to (L/U) . The
dimensional form of T^- may be written as (ij>/T)A. Therefore, the variance
_ _
of dimensional form $$- is proportional to (ipL/UT)2. The variance of i|> is
— o t ^
proportional to ij> , so that the ratio of a/a is given by
a/a « (L/UT)2 .
This ratio was taken as unity in developing Eq. (23) and is so taken
here. That choice means that the time scale of variability and of sampling
is proportional to the transport time. Typical values of L (=4 km) and
U (=3 m/s) gives T ^ 20 minutes. If the sampling rate is greater than the
transport rate, then T should be proportional to the sampling rate. If the
analysis is being used to update the pollution distribution based on recent
measurements, the choice of T should depend upon the update frequency.
In practice, however, a/a should be large where there are observations
since there is no substitute for a good observation; yet, it should not be
so large as to force an observation which is inconsistent with other obser-
vations or with the dynamical considerations into the analysis. Where there
are no observations, a/a should be small, thereby relying on the dynamics
of the system to determine the concentrations.
71
-------
3.3.5 Test Results
There are many assumptions and interpretations left to the user which
have a differing influence on the solution to the analysis equation. Either
boundary condition is acceptable and will give a minimum,but not necessarily
equal, error variance for the integral. The solution will be different in
each case. The choice of space scales, L and Z, may be important to the
solution as they affect the characteristic numbers TI and £. The magnitude of
the weighting terms a and o, while left to the user, should be representative
of the situation being modeled. The distribution of sources affects the
"measured" concentrations, the error integral, and the magnitude of the
characteristic numbers. The atmospheric stability assumption also governs
the magnitude and distribution and n and £. The choices are interdependent.
A series of test cases were developed to verify that the solution can
be obtained to determine the effects of different assumptions or parameter
choices upon the solutions, and finally, to examine the applicability of the
model in a simulated urban area where there are a limited number of obser-
vations which contain some error. In these tests, ^*, is the concentration of
the steady-state numerical solution obtained for i|i(x,y,0), ty will be the
"observed" concentrations throughout the area which may be wholly or partially
dependent upon iji*. and if> will denote the concentration given by the solution
of the Euler-Lagrange equation. The urban area modeled is 9 km (south-north)
by 20 km (west-east). .The center of the area source is at x = 7.5 km,
y = 4.5 km. The basic source configuration (Appendix B) was chosen so that
the effects of gradients in the emissions patterns can be investigated.
Emissions do not occur out of the grid area.
The distribution of concentrations for a basic source configuration in a
neutrally stable atmosphere (the Pasquill-Turner D stability) and a background
concentration ip = i|i is shown in Figure 13. Since the distribution is symme-
trical about the downwind centerline (i.e., y = 4.5 km), cross sections along
the y = 5 km axis will characterize the solutions which are obtained. The
distributions of S^y-*, and £ along y = 5 km determined from that distribution
when L = 50 m are given in Figure 14. Upwind of any source contribution to
concentrations at z = Z, background relationships prevail and S_ = 22.35 Z.
Z
72
-------
U)
0
Figure 13.
10
X, km
15
Concentration isopleths (yg/m ) at ground level from emissions (10 gm/s) normally
distributed within the shaded area, at neutral stability and a constant wind, U,
of 3 m/s.
20
-------
0.2,—r
20'0
0,
ex.
-,27
-25
-23
-2, 3,
-19
- 17
Figure 14. Horizontal distribution of £ ( ), S (— • — •),
—1 z
and Y ( ) along y = 5km .for Case I.
74
-------
The details and results of the sensitivity analyses for the effects of
a/a, L, source size, boundary conditions, Z, and observational errors are
given in Appendix C. Table C-l of Appendix C summarizes the input data, A
brief description of each case study follows:
Case I: a/a When a/a = 10, the analytic solution was returned. As a/a
decreases, the influence of the observations decrease so the
difference between the observed and the solution increases.
The maximum value of the analysis is displaced slightly
downwind from the analytic solution, so the analysis
predicts lower concentrations upwind of the source and
predicts higher concentrations downwind of the source.
Total mass is approximately conserved.
Case II: L The effects of the constraint equation were studied by
choosing a/a = 0.01. Under these conditions, the
analysis reduces to a one-dimensional Poisson equation
whose solution depends upon the choice of L. When
L * 1 km, the solution was achieved in fewer iterations,
but differed only slightly with the solutions when
L = 4 or 10 km.
Case III: (a) Source Configuration:
The total amount of emissions was kept constant but was
distributed over larger areas by increasing the standard
deviation, qx> of the Gaussian-shaped emission pattern
from 1.25 to 2.50 to 5.0 km. As the emissions became
i
more uniform, i.e., qx increased, the concentration
estimates of the analysis equation (a/a = 0.01) approached
the analytical solution for that emission pattern.
Case III: (b) Boundary Conditions:
The boundary conditions were changed from invariant
(4>= constant) to constant flux conditions across the
boundaries, and the same source-dependent experiments
were rerun. The greatest error and error variance
75
-------
occurred in the q = 1.25 case when a negative concen-
tration was analyzed at the upwind boundary. That was
caused by the lack of emissions there. At the downwind
boundary,maintaining the gradient required that the
concentrations there were overestimated by 1.67 times
the observed. As the emissions became more dispersed
(q = 2.5),the analysis approached the analysis with
constant boundary conditions and the analytic solution.
For q = 5.0>the analyses with different boundary
conditions were indistinguishable.
Case IV: Effects of Z
The effects of choosing different values of Z, the
height of the urban volume upon the analysis equation
(a/a = 0.01) were investigated using Z = 100 m, 200 m,
and 31.5 m (= Sz) and comparing the resulting analyses
with solutions when Z = 50 m (*v» /2~ S_). As Z increases,
2
the analysis progressively and substantially under-
estimates the maximum concentration and displaces it
further downwind. When Z = 200 m, the analysis under-
estimated the true concentration at every point. The
smaller value of Z showed good agreement upwind of the
source area but consistently overestimated the concen-
tration downwind of the peak concentration. The choice
of Z^/2 S gives the best overall agreement, minimum
Z
error variance, of all of the values used.
Case V: Observational Errors
Random errors were drawn from a normal distribution
with zero mean and 1(1 standard deviation and were
D
added to the analytic solution at each point forming
the "observed" distribution. The analysis, with
a/a = 1.0, reduced the magnitude of the perturbations
and the analysis approached the analysis of an
unperturbed field. Thus, the OVAM filters the high
frequency perturbation.
76
-------
The results of the sensitivity analyses and tests clearly show that the
OVAM works well with a well-defined initial distribution of concentration,
i.e., when there are good estimates at every point. However, in practice,
only a few observations are known in the urban area. To test the
technique more realistically, observations (^ = <(;*) were selected and
retained at 20 randomly selected grid points. The weight ratio, a/a, at those
points was assigned a value of 1.0. "New observations" were obtained at the
remaining boundary and interior grid points by interpolation from the
retained observations using a 1/R^ weight function, e.g.,
where
R2 = (i-£)2 + (j-m)2 (# 0)
fy ( i,m) are the retained values of fy and the summation is over all retained
variables. For these "new" observations a /a =0.1.
In this manner, the solution for iji is initially estimated, based solely
on the observations. Where there are observations, greater weight is accorded
than where there are interpolated values. The relaxation coefficient, fl,
changes accordingly at each point of the interior. An analytical test to
guarantee a convergent solution with a variable, fi , was not performed. The
experience of Cases I-V suggested that the relaxation would converge to a
solution even if a -»• 0.
The initial distribution of the "observed" concentrations, i|> ; and the
location of the retained observation is shown in Figure 15. The analytic
solution, tp*, for the same conditions is shown in Figure 13. Obviously, ty is
a poor representation of ij;*, especially in the areas near the source. This
clearly points out the potential problems of using a straightforward inter-
polation scheme to obtain concentration estimates in an urban area.
The solution concentration along y = 5 is shown with the analytic
solution, the initial distribution, and the results obtained in Case II in
Figure 16. The solution is an obvious improvement to the observations in
77
-------
I
M
>-
oo
2 '
i i
10
X.km
15
20
3 -2
Figure 15. Concentration distribution (yg/m ) with x as determined using R interpolation
of randomly selected observations (*) from the distribution of Figure 12.
-------
Figure 16. Concentration distribution Cug/m ) with x for
Case VI showing <|>* ( ) , iji(x x),
i|i from Case II (••••)* along y = 5 km. Asterisks
(*) indicate locations of observations on the
y axis.
79
-------
the domain. In this case, 3 of the 15 retained observations lie
along the y = 5 km axis, at x = 3, 4 and 14 km. At 14 km, the solution
(8.006 yg/m^) is nearly identical to the observed value (7.928 yg/m^) and off
o
by about 0.5 yg/mj at the other two locations. In general, this solution
behaves much the same as those cases when a/a = 0.01. The different boundary
values and the observation at x = 14 km seem to anchor the solution, letting
the constraint equation fill in the remainder of the distribution.
The expression p. . - (fy. . - ij;*. .)/i|»*. . was used to measure the error of
^•jJ 1-3 3-J !J
the analysis relative to the analytic solution. The initial error average
showed a mean of 0.265 and a standard deviation of 0.743. The solution had a
mean error of 0.169 and standard deviation of 0.450. The weighted mean
square error of the analysis, i.e., the numerical value of the integral, EV,
decreased to 10.8 percent of its original value, as shown in Figure 17.
A different set of observations were picked from fy*, and the procedures
used were used to get new observations. Those data were then perturbed (as
in Case V) as a final simulation of a case using field data to give the
distribution of Figure 18. The observations, an approximate solution, and
the analytic solution are shown in Figure 19.
80
-------
4x I04
-------
oo
to
Figure 18.
X,km
3 ~2
Concentration isopleths (yg/m ) determined from an R interpolation of perturbed
observations at randomly selected data points (*) of the distribution of Figure 25.
-------
0)
o
o
o
I I I I I I I I ! I I I I
Figure 19. Concentration distributions with x for Case VII
showing i|)*( ), * ( ) and $ ( ) along
y = 5 km. Asterisks (*) indicate locations of
observations on the y = 5 km axis.
83
-------
4.0 TEST AND EVALUATION
4.1 Operational Optimum Network
It was economically unfeasible for RTI to set up a 19-station
network in St. Louis. To reduce the cost of the field program required to
validate these results, it was decided to use stations in either the RAPS or
city/county air pollution network (see Figure 20) when there was overlap or
near overlap. This tractable version of the OSN is designated the OSN*, and
the station plot is given in Figure 21. Once again, the numbers next to the
stations are used to designate stations which are to be deleted if less than
optimum subclasses are desired. The analysis of the OSN* was discussed in
Section 3.2.7.
It should be noted that one of the most important aspects of the entire
technique was the selection of the regression model to describe the airflow.
The 13 term model was selected because 1) the 13-term model was
considered the smallest model which would yield reasonably good results,
2) most urban areas could not suffer the economic burden that would be
required to establish an OSN using a regression model with many more terms
than 13 and 3) models with as many as 10 more terms did not yield an
extremely significant improvement over the 13-term model.
The sampling network selection technique can be also applied using
models with a smaller number of terms in order to establish an initial
network with a small number of stations. This might be done because a certain
city may only be able to initially afford R number of stations, but the
minimum regression model which gives adequate results may have S number of
terms where S < R. So, for economic reasons, a regression model with less
terms than S must be used. This approach is very limited because a point will
be reached rather quickly when, no matter how many stations are added at some
later time, the cost of the stations will over-balance the improvement; i.e.,
a base error exists due to the missing terms in the model which cannot be
reduced by improving the estimate of the coefficients in the model (adding,
stations) but only by adding the required terms. Unfortunately, if the terms
are added as the network is enlarged, the OSN obtained from the inadequate
model may not be optimum for the improved model.
84
-------
v RAPS Station and Number
* City-County Station and
Figure 20. Station array in the RAPS network and in the City/County
Network in St. Louis, Missouri.
85
-------
0-Overlap or Near
Overlap with (•—
RAPS Stations »
(J)-Over lap or Near ^
Overlap wich City/
County Scations
^-So Overlap
\
\
Figure 21. Stations in the OSN* for St. Louis, Missouri.
86
-------
4.2 Field Programs
The major emphasis of the second phase of this research project involved
the preparation and execution of a summer and winter field program in
St. Louis. These field programs were held during the period when EPA was
performing an intensive study in St. Louis (July and August 1975 and February
and March 1976). During this time, there was a concerted effort to maintain
a high level of performance of the RAPS stations.
As indicated in Figure 19, three stations did not overlap with either
a RAPS station or a city/county air pollution station. It was necessary for
RTI to establish its own stations in these locations. Sampling stations were
located on the grounds of Incarnate Word Academy in northwest St. Louis
county; on the grounds of Kenrick Seminary in southwest St. Louis county; and
on the grounds of the East Side Sanitary District's South Pumping Station in
East St. Louis, Illinois. Figure 22 shows the location of these stations
(the black dots).
The facility on the grounds of Incarnate Word Academy consisted of the
RTI Mobile Ambient Air Monitoring Laboratory which is a 31-foot motorized
laboratory equipped for ambient air monitoring. A Beckman 6800 was used to
measure THC, CH^ and CO at this and all other locations maintained by RTI.
The Beckman 6800 was chosen to conform with the instrumentation in the RAMS
stations and, thus, avoids a possible bias due to different instruments.
Outputs of the respective instruments were recorded in analog and digital form
in the mobile laboratory. An onboard mini-computer provides for real time
data processing and display of selected parameters in terms of their 5-
minute values as well as computed hourly averages. The facility contains a
10-meter crank-up tower for the wind system. Wind speed and direction was
recorded in analog (strip chart recorder) and digital form. The facility was
located at the southwest end of the baseball field on the campus of Incarnate
Word Academy. There were no obstructions to the wind flow to the north and
east for more than a mile. The school complex was located to the west
approximately 100 m from the van. To the south was a housing project.
The facilities at both the South Pumping Station and on the grounds of
Kenrick Seminary were similar. These consisted of an 8-foot by 7-foot
metal house and a 10-m tower. The metal house contained facilities for both
87
-------
204
®
WIEDMASIl
^-Overlap or Near
Overlap with
RAPS Stations
-Overlap or Near ><-x"J_— .-•^""
Overlap with City/
County Stations
-No Overlap/RTI Station
Figure 22. Location, name and number of stations maintained
by RTI.
88
-------
heating and cooling, and housed the Beckman 6800 and the strip chart recorders,
The air intake for the Beckman 6800 and the wind system were placed at the
10-m level of the tower.
The facility at Kenrick Seminary was located in a wooded portion of the
campus, approximately 100 m southwest of the seminary complex and about 100 m
north of the heat plant for the seminary. The South Pumping Station is
located in the rural portion of Cohokia, Illinois. The RTI station was
located near the top of a flood control levee about 100 m southeast of the
pumping station, the only obstruction to the flow in the general area.
Through the cooperation of the St. Louis City/County Air Pollution
Agencies, Beckman 6800's were also placed in each of the city/county stations
which overlap with the OSN. As in the case of Kenrick Seminary and the South
Pumping Station, the data were recorded on strip chart records. Wind data
were collected at these stations by systems already on site and maintained by
the city/county agencies. The exposures at these stations were generally good,
Dynamic calibration techniques were used to calibrate each Beckman 6800
air quality chromatograph at specific intervals. The Beckman analyzers at
each of the seven stations manned by RTI were calibrated every two days. The
calibration consisted of a dynamic zero and an upscale calibration at
80 percent of the instrument's range. These instruments were also calibrated
by EPA. The calibration data were used to update the transfer equations for
converting strip chart readings to physical units.
A dynamic calibration procedure was also used by EPA to calibrate their
Beckman 6800's. The EPA Beckman 6800"s were calibrated daily to obtain
updated transfer equations. Their daily calibration consisted of a dynamic
zero and an upscale calibration of the instrument's range. Multi-point
calibrations were accomplished weekly.
4.3 Analysis of Data
To date, all data from the winter and summer field programs have been
reduced. The analysis described below is in progress. Results of these
analyses will be presented in later reports.
89
-------
1) The data from the field programs will be used to estimate
wind speeds and directions over the domain of the network
by fitting the regression model to data from the OSN*.
2) These predicted winds will be statistically compared to
actual data from stations not in the OSN*.
3) Emissions inventory obtained for St. Louis, the wind
distribution obtained using the data from the OSN*,
and the air pollution data from the OSN* will be
combined with the objective variational analysis model
to yield the air pollution distribution over the domain
of the network.
4) Air pollution concentrations established using the
objective variational analysis model will be compared
with observed data at locations where air pollution
data are available but are not stations of the OSN*.
90
-------
APPENDIX A
A THEORETICAL STUDY OF THE ST. LOUIS HEAT ISLAND:
THE WIND AND TEMPERATURE DISTRIBUTION
Reprinted from the Journal of Applied Meteorology,
Vol. 15, p. 417-440, by permission of the American
Meteorological Society.
91
-------
Reprinted from JOURNAL OF APPLIED METEOROLOGY, Vol. 15, No. 5, May 1976
American Meteorological Society
Printed in U. S. A.
A Theoretical Study of the St. Louis Heat Island:
The Wind and Temperature Distribution
FRED M. VUKOVICH, J. W. DUNN III AND BOBBY W. CRISSMAN
Research Triangle Institute, Research Triangle Park, N. C. 27709
(Manuscript received 21 July 1975, in revised form 26 January 1976)
ABSTRACT
A three-dimensional primitive equation model was used to study the St. Louis heat island. In this paper,
the influence of synoptic wind speed and wind direction on the heat island is presented. With respect to
the synoptic wind speed, it was found that the temperature and wind distribution associated with the St.
Louis heat island changed markedly as the wind speed increased. When the synoptic wind speed was small,
the intensity of the heat island was independent of the wind direction. However, for large syrioptic wind
speeds, the intensity of the heat island changed, and the change was dependent on the wind direction. These
changes were due to the influence of the local topography.
1. Introduction
It has been shown that an urban complex acts as
a heat reservoir producing the so-called urban heat
island (Landsberg, 1956; Chandler, 1964, 1965, 1967;
Woollum, 1964; Lowry, 1967; Mitchell, 1961; Clarke,
1969). Some rather generalized theoretical studies have
shown that the urban heat island should create a
solenoidal circulation (Delage and Taylor, 1970; Vu-
kovich, 1971, 1973, 1975; Bornstein, 1972); i.e., the
differential heating between the urban region and the
surrounding suburbs and rural region sets up a mass and
thermal contrast which results in the conversion of
available potential energy into kinetic energy through
convective circulation. These results are generalized
because they do not apply to any specific urban heat
island. Furthermore, the results treat a two-dimen-
sional heat island, and are therefore limited in
application.
This paper gives some of the results of a three-
dimensional numerical simulation of the heat island
circulation in St. Louis, Mo. The numerical model
incorporates the basic forcing functions characterizing
St. Louis; i.e., the topography, the building heights,
the roughness variations, and the temperature per-
turbations describing the heat island. St. Louis was
chosen for study because the verification of the various
models developed in the research program in which
these numerical simulations were performed, will be
performed using data from the Regional Air Pollution
Study (RAPS) which has begun in St. Louis.
2. The model
The basic equations governing the dynamics and
thermodynamics for dry convection in the model are
as follows:
d« du du du
—\-u—\-v—H0—=
dt dx By dz
dv di> dv dv
dt dx dy dz
AV—+AZI— (1)
36 dd
dt dx
—+AV,— (2)
dO dd d*6
v \-v>—=A z \-A „
dy dz dx2 dy2
dp dw
p -- f-2
dz2 dz dz
32p drd(pu) d(pv)~\
dz2 dzl dx dy J
dp_= _PS/fr\'
dz ~RO\.p /
(3)
(4)
(5)
The symbols used are the conventional ones used
in meteorology. In the above equations, the horizontal
diffusion coefficients for momentum and heat were
considered identical. The local rate of density change
was set to zero (6V/6V=0), an assumption which
yielded Eq. (4) for the vertical velocity. The most
important approximation used in the model is the
hydrostatic assumption through which the pressure
field was derived. The use of. this assumption was
justified since we dealt with a relatively shallow
92
-------
418
JOURNAL OF APPLIED METEOROLOGY
atmosphere, i.e., the top of the atmosphere is 4.0 km
in the model. Furthermore, the vertical accelerations
produced by the forcing functions in the urban area
are less than 10~5 g which suggests that the hydrostatic
approximation is reasonable in any case.
Perturbation principles were applied to the above
equation. In the perturbation analysis, variables des-
ignated by the overbar (p, for instance) are used
to represent the synoptic state which is assumed to
be both in geostrophic and hydrostatic equilibrium.
Primed variables are used to represent the perturba-
tions produced by the available forcing functions in
the city area.
The horizontal dimensions of the grid on which
the model is structured are 144 km by 144 km which
more than adequately covers the urban and suburban
regions of St. Louis. A nested grid is employed in
which the central region (20X20 grid points centered
over the city) has a grid spacing of 1 km. Outside
the central region, the grid spacing increases in both
* and y according to the expression, 2* [km]], where
*=1, 2,3, 4 until the grid spacing equals 16 km. At
this point, the grid spacing is held constant for three
more grid points. Fig. 1 shows a portion of the
total grid.
In regions where grid spacing is constant, centered
difference formulas are used for nrst- and second-order
space derivatives. In the region of changing grid
spacing, upstream differencing is employed to give
first-order space derivatives, and non-centered dif-
ferencing is used to determined second-order space
derivatives.
At the horizontal boundaries, mass, momentum and
heat are allowed to flow out except at the upstream
boundary; i.e.,
at x= — 72 km (upstream boundary)
at x=+72 km (downstream boundary)
du dv de dp
dx dx dx dx
at y=±72 km (lateral boundaries)
du dv BB dp
— =_=_= — =o.
dy dy dy dy
\
FIG. 1. Horizontal grid. The grid spacing in the center, uniformly spaced area is
1 km. For the first outer set of grid points, the grid spacing is 2 km; for the second
outer set, 4 km; and for the third outer set, 8 1cm. The hatched region indicates
the area of increased roughness in the city of St. Louis.
93
-------
MAY 1976
F . M . V U K 0 V I C H , J . \V . DUNN AND B . \V C R I S S M A N
419
There are eight vertical levels including the lower
boundary and the upper boundary. The vertical levels
are s=0, 100, 300, 600, 1000, 1500, 2500 and 4000 m.
These are based on a definition of the 3=0 plane for
the local area. This will be discussed later. Upstream
differencing is used exclusively in the vertical. The
upper boundary conditions are:
«=tZ, p=p, 0=0
4km.
The lower boundary conditions will be discussed in
detail later.
The primary forcing functions in urban areas are
heating and local terrain effects. In order to describe
the effect of terrain, the topography is parameterized
for the model. The average height of the local terrain
and the average building heights were computed in
areas about each grid point. The building heights
were determined from Sanborn Maps, courtesy of the
St. Louis City Planning Office. The size of the area
depends on the grid spacing; i.e., for a 1 km grid
spacing, the area about the grid point is 1 km2, and
for a 2 km grid spacing, it is 4 km2, etc. The minimum
averaged terrain height was determined from the set
of average values, and this value was subtracted from
each terrain height before the average building height
was added. The result, after adding in the average
building height, is the reference height H(x,y) for the
area around each grid square. The 2=0 plane is de-
fined where H(x,y)=Q which was located over the
Mississippi River in St. Louis. Fig. 2 gives the dis-
tribution of H(x,y) used in the model. At z=H(x,y)
and at the side walls of an obstacle that may be
higher than the first grid point above 2=0 where the
prediction equations are integrated, the following
boundary conditions are used for the wind:
If z*— H(x,y)<15 m, where z* is the height of the
first grid point above H, the prediction equations
became unstable for the time increment used in
integration. Therefore, values of 3(x,y) which did not
satisfy this criteria were reduced to H(x,y)=z*—lS.
For St. Louis, all values of H(x,y) are less than 100 m,-
which is the height of the first level above the z=0
plane in the model. About six H(x,y) values exceed
the stability criteria and were reduced.
In order to characterize the effects of an urban
heat island, a field of potential temperature depar-
tures, 80(x,y), defined at z=H(x,y), were used. For
St. Louis, this field was obtained from data presented
in a Stanford University report (1953) and from
available surface weather data for the region. The
resultant field of 60(x,y) determined from these data
are average values and are given in Fig. 2. The nega-
tive values of 50(x,y) to the west, north and south
are due to radiational processes in Forest Park, in
the combined areas of Bellefontaine Cemetary, Calvary
Cemetery and O'Fallon Park, and in Tower Grove
Park. The urban heat island represented by the 5B(x,y)
values was allowed to develop gradually with time
by letting
where 7 is the rate constant and for this study is
such that &ff(x,y)»M(x,y) in 3 h.
The boundary condition for the potential tempera-
ture distribution at z=E(x,y) is
If an obstacle was higher than one or more grid levels,
then at z=H(x,y) and at the side walls of the obstacle
The initial synoptic field represented by the overbar
was determined as if the terrain were not present.
The heat and momentum fluxes at z=H(x,y) were
determined through simultaneous solution of the
boundary layer profile equations for the friction
velocities M* and v+ and the scale temperature T*.
These parameters are proportional to the momentum
and heat fluxes, respectively, at the surface. The
profile equations used were
eu*
*f I (Z~H\
u=— ml ]
kl \ z0 /
9U*
where U*= («*2+v*2)» and was allowed to vary in
the following manner. In regions outside the city
z0=0.0015 H(x,y).
Using the above expressions and for the given values
of H(x,y) in the St. Louis area other than, in ;the
built-up sections, 1 cm 4 z<> ^13 cm. In the buiit-up
sections :
which yielded values of s0 as large as 1 m. This tech-
nique is an adaptation of one suggested by Lettau
(1969). In the model only one major built-up section
was designated. This region was bounded by Forest
Park to the west, the Mississippi River to the east,
approximately St. Louis Avenue to the north, and
Russell Avenue to the south (see shaded area, Fig. 1).
94
-------
420
JOURNALOF APPLIED METEOROLOGY
95
-------
MAY 1976 F. M. VUKOVICH, J. W. DUNN AND B. W. CRISSMAN 421
The formulas for the vertical eddy diffusion coeffi- 4 _
dents are
Azl=l
Av,=l
du
dz
dv
dz
dV
dz
where 72=MZ-|-5Z. The mixing lengths are given by
dz
liHl dO
where ai=18, a2=H> and Ri is the Richardson
number.
Since the magnitude of diffusion brought about by
subgrid-scale motion is dependent, among other things,
on eddy size, and since larger grid spacings incorporate
larger subgrid-scale eddies, the horizontal diffusion
coefficients are directly proportional to the square
root of the local mean grid spacing. The constant of
proportionality depends mainly on the boundary layer
stability.
The initial wind field is assumed to have a wind
direction parallel to the x axis and a speed which
varies in z alone. The initial potential temperature
distribution is given at x = — 72 km and y=—72 km.
Since the synoptic state (initial conditions) was as-
sumed in geostrophic and hydrostatic equilibrium, the
necessary initial pressure and potential temperature
fields were derived through integration of the geo-
strophic and hydrostatic equations.
The initial conditions are represented by two cate-
gories of stability and five different wind profiles.
Fig. 3 shows the vertical distribution of potential
temperature. Profile 1 describes an atmosphere with
a shallow, near-isothermal layer from the surface to
the 100 m level. Above the 100 m level, the atmo-
sphere is stable with a lapse rate equivalent to the
Standard Atmosphere lapse rate. Profile 2 describes
an atmosphere in which the first 600 m has an adia-
3 -
2 -
1 -
i i i i i i i i
286 290 294 298 302
Potential Temperature (°K)
4-,
0 2 4 68
Wind Speed (m s'1)
FIG. 3. The vertical distributions of the initial potential tem-
perature (upper) and the initial wind speed (lower) used to inte-
grate the equations of motion.
batic lapse rate. Above 600 m, the lapse rate is again
equivalent to the Standard Atmosphere lapse rate.
Profile 1 may be considered an evening sounding and
Profile 2 a noon sounding.
The initial wind profiles used are also shown in
Fig. 3. "Since the urban heat island normally develops
in the late afternoon and is most intense in the evening
when the boundary layer is stable, the results to be
presented here will employ the potential temperature
Profile 1 combined with wind Profiles 1-4. In a future
paper, the influence of an adiabatic boundary layer
FIG. 2. Upper: The distribution of topography heights (m) above the z=0 plane
in the St. Louis region. The values of topography include both natural topography
and area-averaged building heights.
Lower: The distribution of surface temperature (°C) used to generate the urban heat
island in St. Louis.
96
-------
422
JOURNAL OF APPLIED METEOROLOGY
VOLCME U
Fio. 4. Upper: The horizontal distribution of the quasi-steady-state, perturbation
potential temperature (K) at the 100 m level for the northwind case. The initial condi-
tions were .P(l,l).
Lower: The vertical cross section of potential temperature (K) along the line AA'.
97
-------
MAY 1976 F . M . V U K 0 V I C H , J . \V . DUNN AND B . W . C R I S S M A N
423
FIG. S. As in Fig. 4 except for initial conditions P(2,l)>
98
-------
424
JOURNAL OF APPLIED METEOROLOGY
VOLUME I:
-!US-
1 I I I I I I I I I I I I F 1 I i I
-11 -It I 11
FIG. 6. As in Fig. 4 except for initial conditions P(3,l).
99
-------
MAY 1976 F. M. VUKOVICH, J. W. DUNN AND B. \V. CRISS.MAN
425
! ! mTTTTTTTTTTT m I m !
-II -I 0 I
FIG. 7. As in Fig. 4 except for initial conditions P(4,l).
100
-------
426
JOURNAL OF APPLIED METEOROLOGY
VOLUME 15
101
-------
MAY 1976
F . M . V U K 0 V I C H , J . W. DUNN AND B . W . C R I S S M A N
427
will be presented. In the text, the notation P(2, 1)
will be used to specify the initial wind and tempera-
ture profiles, i.e., wind Profile 2 and potential tem-
perature Profile 1 in this example.
In order to study flow in the city for various wind
directions, the H(x,y) and 58(x,y) fields which define
the city area were rotated to effectively yield a new
angle for the wind. Eight wind directions were used,
representing an eight-point compass rose. For the time
integration, forward time differencing was combined
with the leapfrog method to control the divergence
of the solution at even and odd time steps. For the
purposes of this paper, the quasi-steady-state solution
will be presented. In future papers, the time evolution
of the heat island circulation will be presented. Steady
state was achieved in S h or more. This varies some
from case to case.
3. The thermal distribution
The following describes the results of integrating
the equations of motion using the initial conditions
P(l,l) through P(4,l); i.e., the boundary layer was
stable from case to case, and only the wind speed
and shear changed. The initial (geostrophic) wind
direction is from the north hi the initial cases. How-
ever, later in this section, the influence of wind
direction will be discussed. The horizontal diffusion
coefficients are of the order 10s m2 s"1 in each case.
Fig. 4 shows the perturbation potential temperature
distribution at the 100 m level associated with the
St. Louis heat island for initial conditions P(l,l). The
temperature perturbation is concentrated in the city
limits since there is little thermal advection. The
central temperature perturbation is 0.92 K. Fig. 4
also shows the vertical distribution of potential tem-
perature along the line AA'. Potential temperature
data were acquired at grid points nearest to the line.
The heat island is depicted as a region having an
unstable boundary layer surrounded by a stable
boundary layer in the suburbs and rural regions. The
top of the heat island is found over the center of
the city. In this case, it has been determined from
detailed analysis that the top of the heat island is
slightly higher than the 300 m level. The higher
penetration of the heating hi the center of the city
is thought to be a combined effect of increased rough-
ness and the larger SO values.
For initial conditions P(2,l), the thermal plume at
100 m (Fig. 5) is oriented north-northwest and south-
southeast, and extends about 10 km southeast of the
center of St. Louis, due to the increased thermal
advection. It is believed, however, that the actual
downstream extent of the heat plume in this case
and in the preceding case should be greater since
advection is hampered by the dampening produced
both by using upstream differencing and by expanding
the grid. The central perturbation temperature in this
case is 0.85 K, a decrease of 0.07 K from the previous
case. The vertical cross section (Fig. 5, lower) along
AA' shows that the layer between the 100 and 300 m
levels is less stable downstream of the heat island
than upstream. This is due to the thermal advection
and is identical to the effect described by Clarke
(1969) in the Cincinnati boundary layer. Strong sta-
bility developed near the surface immediately down-
stream of the heat island due to the upper level heat
transport. The maximum vertical extent of the heat
island (~300 m) was located in the central regions
of the city where there was increased roughness and
large values of SO.
As the initial wind speed and shear is increased
further, the heat plume from the St. Louis island is
displaced further downstream [Tig. 6 gives the results
for P (3,1)3- The displacement is about 15 km. The
maximum positive temperature perturbation is 0.83 K,
which is about a 0.02 K decrease from the previous
case |~P(2,1)]. The vertical cross section along AA'
again demonstrates that the stability in the layer
between 100 and 300 m is stronger upwind of the
heat island than downwind. The strong near-surface
stability found in the previous analysis downwind of
the heat island and produced by the upper level heat
transport is not present in this analysis. This may be
a result of the orientation of the cross section used.
The top of the heat island is at about the 250 m level.
In the case when initial conditions P(4,l) are used,
the heat plume is displaced 22 km downwind of the
center of St. Louis (Fig. 7). The maximum positive
temperature perturbation in this case is 0.80 K. The
general decline of the temperature perturbation from
case to case is directly related to the increased wind
speed and shear (Vukovich, 1975). The vertical cross
section indicates weaker stability downstream of the
heat island between 100 and 300 m as before. A strong
near-surface stable layer, produced by the upper level
advection of heat, extends 10 km downstream from
the outer fringes of the city. The top of the heat
island is at approximately the 200 m level. It is sug-
gested that the decrease in the vertical extent of the
heat island from case to case is due to the increase
in wind speed. In general, the top of the heat island
found using this model agrees well with that de-
termined using the data given in the Stanford Uni-
versity report (1953).
Fio. 8. Upper: The horizontal distribution of the quasi-steady-state horizontal wind
vector at the 100 m level for the north wind case. The initial conditions were P(l,l).
Lower: The horizontal distribution of the quasi-steady-state wind vector at the 300-m
level for the north wind case. The initial conditions were P(l,l).
102
-------
428
JOURNAL OF APPLIED METEOROLOGY
VOLUME 1:
FIG. 9. Upper: The horizontal distribution of the vertical velocity (on s~l) at the
300 m level for the north wind case. The initial conditions were P(l,l). Positive values
indicate upward vertical motions.
Lower: The vertical cross section of vertical velocity (cm s~l) along line AA'. Positive
values indicate upward vertical motions.
103
-------
MAY 1976
F . M . V U K 0 V I C H . J . \V . DUNN AND B . W . C R I S S M A N
429
Changes in wind direction have little effect on the
development of the heat island when initial conditions
P(l,l) or P(2,l) are used, i.e., the heat island is con-
fined mainly to the city limits with only slight down-
wind displacement. When the wind speeds increase,
as in initial conditions P(3,l) and P(4,l), the down-
wind transport of the heat plume and other factors
describing the thermal distribution are generally the
same for the other wind directions with the exception
of (geostrophic) winds from the east and southeast.
In those cases, the maximum positive temperature
perturbation was smaller (~0.60K), and the heat
plume was displaced downstream ~6 km less than
for the other cases. The downwind region between
the 100 and 300 m levels of relatively less stable air
than the upwind region did not extend as far down-
stream in this case as in the other cases. The down-
wind stability in that layer became equivalent to the
upwind stability about 12 km from the heat island
center. The top of the heat island was at the 250 m
level. The differences were attributed to the combined
effects of increased friction and adiabatic cooling
produced by flow over the hills west of the city. The
topography also affected the wind field. This will be
discussed in the next section.
4. Wind distribution
The cases demonstrating the wind distribution asso-
ciated with the St. Louis heat island circulation will
show the results of integrating the equations of mo-
tion for the same initial conditions used for the
thermal distribution when just the effect of wind
shear and speed are of interest. Fig. 8 shows the
distribution of the horizontal wind velocity at the
100 m level for the initial conditions P(l,l). Various
streamlines are drawn in the figure to help in the
analysis of the flow. Due to the weak initial wind
field, a rather pronounced zone of convergence and
strong cyclonic curvature in the flow are established
in the central region of the city. The strongest winds
are found upwind (north) of the city's central region
and are a result of relatively large pressure gradient
accelerations. The weakest winds are found southeast
of the city's center and, for the most part, are a result
of pressure gradient accelerations opposite to the
direction of the synoptic flow. The 300 m level is
characterized by divergence (Fig. 8, lower). The
strongest winds at this level are found downwind
(south) of the city whereas the weakest winds are
north of the city. The horizontal distribution of the
vertical velocity at the 300 m level (Fig. 9) shows
1
FIG. lOa. As in Fig. 8 (upper) except for initial conditions P(2,l).
104
-------
430
JOURNAL OF APPLIED METEOROLOGY
VOLUME 15
i r i r r r i i i
\ nnrni
r mtfrrm
r r-f r r r r r/// T \r in 1
r i i !
FIG. lOb. As in Fig. 8 (lower) except for initial conditions P(2,l).
a concentrated region of upward vertical motion over
the city's central region. The maximum upward ve-
locity is 19 cm s~'. The vertical cross section of the
vertical velocity along line AA' (also shown in Fig. 9)
suggests that the upper boundary of the St. Louis
heat island circulation in this case is at about the
1 km level (where the vertical velocity is zero).
The horizontal wind field at the 100 m level for
initial conditions P(2,l)vgiven in Fig. lOa shows
a zone of convergence in the central regions of the
city, but the strong cyclonic curvature in the flow,
obvious in the previous case, is not prevalent hi this
case. The strongest winds are again found upwind
(north) of the city, but the weakest winds are located
in the center of the city. The relative difference in
wind speed upstream and downstream is again due
to pressure gradient accelerations. The pressure gra-
dient acceleration opposes the increased factional
effect in the urban boundary layer due to the thermal
instability, and apparently dominates the flow. At
the 300 m level (Fig. lOb) divergence prevails, but
there is no obvious difference between upstream and
downstream wind speeds. The initial wind speeds in
this particular case are similar to those observed by
Ackennan (1972) in a case study of flow over St.
Louis at night; however, the wind direction was from
105
-------
MAY 1976 F. M. VUKOVICH, J. W. DUNN AND B. \V . C R I S S M A N
431
now <•>
FIG. 11. As in Fig. 9 except for initial conditions P(2,l).
106
-------
432
JOURNAL OF APPLIED METEOROLOGY
VOLUME IS
the east-northeast in her case. She observed con-
vergence over the city at the 100 m level and diver-
gence at the 350 m level.
The vertical velocity distribution at 300 m (Fig. 11)
shows a concentration of upward vertical motion in
the central regions of the city, with some slight down-
wind displacement. In this case, the maximum upward
speed is 16 cm s"1. The suburban and rural regions
west and south of the city are characterized by weak
downward motion (~ —3.0 cm s~l) at 300 m, as is the
case in East St. Louis and Belleville, 111. The vertical
cross section along AA' suggests that the top of the
urban heat island circulation is found at about the
900 m level. Upward vertical motions, for the most
part, characterize the boundary layer below the 300 m
level north and south of the city.
As the wind speed and shear increase further, as in
the case of initial conditions P(3,l), the inflow at
100 m, present in the previous cases, is replaced by
a line of convergence which is displaced downstream
of the city's center (Fig. 12a). The strongest hori-
zontal wind speeds are found immediately upwind of
the city's center. The weakest wind speeds are found
in the southwest corner of the figure, and are a result
of increased frictional effects due to the bluffs in that
region. The wind speeds immediately downwind
(southeast) of the city's center are weaker than those
immediately north, east or west of the city. At 300 m
(Fig. 12b), the horizontal divergence is still obvious,
but upstream and downstream wind speed differences
are not immediately discernable.
The vertical velocity distribution over the city
(Fig. 13) shows some downstream displacement, but
for the most part the region of upward vertical mo-
tion is concentrated over the city. The maximum
updraft velocity is 12 cm s"1. Minor centers of down-
ward motion are located immediately southwest of
the center of upward motion and downstream of the
region of upward motion over the city. According to
the vertical cross section along AA', the top of the
heat island circulation is at about the 1400 m level.
Regions of relatively strong downward motion are-
developed aloft due to the imposition of the strong
initial wind speeds and wind shear.
When initial conditions P(4,l) are employed, the
horizontal wind field at the 100 m level (Fig. 14a)
indicates that the line of convergence is displaced
farther downstream than in the previous case [P(3,l)3-
The horizontal velocities immediately upstream of the
central regions of the city are still somewhat larger
than those immediately downstream. At 300 m (Fig.
14b), the divergence found at this level in the pre-
N
\ \
N
\ \
A
FIG. 12a. As in Fig. 8 (upper) except for initial conditions P(3,l).
107
-------
MAY 1976
F . M . V U K 0 V I C H , J . W DUNN AND B . W C R I S S M A N
433
M 11 I
FIG. 12b. As in Fig. 8 (lower) except for initial conditions P(3,l).
vious cases appears to be located on the east side of
the city.
The vertical velocity distribution (Fig. 15) is char-
acterized by an elongated region of upward vertical
motion extending from the center of the city down-
stream about 20 km and containing three centers of
maximum updraft. The largest updraft velocity is
located in the center farthest downstream of the city
and is 11 cm s~l. The increased vertical velocity
downstream of the city may be a result of the inter-
action of flow over a barrier, since it appears that
the downstream center is located upwind of the bluffs
southeast of St. Louis, and of the heat island circula-
tion. The marked downstream displacement of both
the zone of convergence and the updraft vertical
velocity distribution in the case of strong initial
winds was predicted by Vukovich (1971); however,
his model predicted downward motion over the center
of the city which does not occur in this case. Rela-
tively strong downdrafts are found immediately south-
west of the elongated region of updraft, and weaker
downdrafts are found to the northeast. The pattern
of an elongated region of upward motion with down-
ward motion on either side almost coincides with the
pattern detected by Angell et al. (1973) over Oklahoma
City in the evening with strong horizontal winds.
During the day they found the strongest updraft
velocities immediately upwind of the city with a
secondary maximum downwind. They attributed the
upward motion during the day to the barrier effect
108
-------
434
JOURNAL OF APPLIED METEOROLOGY
FIG. 13. As in Fig. 9 except for initial conditions P(3,l).
-------
MAY 1976
F . M . V U K 0 V I C H , J . W . DUNN AND B . \V . C R I S S M A N
435
of the city rather than to a thermal effect of the
city; but during the evening, they indicated the pos-
sible influence of the urban heat island.
The vertical cross section along AA' (Fig. 15, lower)
suggests that the top of the circulation induced by
the heat island is 1500 m over the city, but somewhat
less downstream. The height of the top of the per-
turbation has increased as the initial wind speed and
wind shear increased. Vukovich (1975) has shown
that such an increase is to be expected when the
characteristic Richardson number decreases; this is
brought about in this case by increasing the initial
wind shear. Furthermore, the magnitude of the up-
draft velocity generally decreased, which may also
be attributed to the increase in the wind shear (Vu-
kovich, 1975).
Changes in wind direction did not significantly
affect the momentum perturbation produced by the
St. Louis heat island except when initial condition
P(4,l) was employed. The effects of changes of wind
direction are best shown using the vertical velocity
distribution. The most intense upward vertical veloci-
ties, among all cases in which initial conditions P(4,l)
were used, were found when the initial wind was
from the northeast (Fig. 16). The tongue of upward
vertical motion extended to a point 30 km down-
stream of the city's central region. Two updraft
centers existed, with the downstream center, located
about 17 km from the city's central regions, having
the largest updraft velocity (14 cm s"1). A secondary
maximum was found immediately upstream of the
center of the city (~6.0 cm s~'). Weak updrafts were
also found over the city to the north and northwest.
Weak downdrafts existed immediately east and west
of the tongue of upward motion. In this particular
case, it is suggested that the line of convergence
produced by the urban heat island was intensified by
topographically induced convergence as air was chan-
neled southward over the Mississippi River between
the bluffs.
The weakest upward motion occurred when the
initial wind was from the west (Fig. 16). As in the
previous cases, the tongue of updraft velocities ex-
tended about 30 km downstream from the city's
center. Two relatively weak updraft centers (3.0 cm
s~x) were apparent—one located over the center of
the city and the other about 8 km from the city's
center. Downdraft centers were found to the west,
northwest, southeast and southwest of the city's
center. The rapid decrease of surface elevation east
\ ' \ \
A
FIG. 14a. As in Fig. 8 (upper) except for initial conditions P(4,l).
110
-------
436
JOURNAL OF APPLIED METEOROLOGY
VOLUME 13
FIG. 14b. As in Fig. 8 (lower) except for initial conditions P(4,l).
of St. Louis due to the river bottom would produce
vertical expansion, horizontal divergence, and down-
ward vertical motions in the boundary layer flow.
This would act in contrast to the vertical contraction,
horizontal convergence, and upward vertical motion
produced by the urban heat island, resulting in a
reduced heat island effect downstream. Though this
effect occurred in both the north and northeast wind
cases, it appears that in this case the location of the
downstream updraft center is immediately over the
region of rapid, surface elevation decrease. This did
not occur in the other cases. Furthermore, in the
west wind case, general downdrafts would be created
over the city due to the general decrease in surface
elevation from the west county region to the Mis-
sissippi River, which could account for the reduction
of the overall intensity of the heat island in this case.
The most complex pattern of vertical velocity was
found when the initial wind was from the east (Fig. 16).
The tongue of updraft velocities present in all the
previous cases is not well defined in this case. Two
updraft centers were evident, but these appeared to
be separate entities. The maximum updraft speed in
the center over the city's center was about 8 cm s~l,
and that in the downstream center was about 7 cm s~l.
The downdraft center immediately west-northwest of
the city's center, is located over a low topographical
region, suggesting either a topographically induced
111
-------
MAY 1976 F. M . VUKOVICH, J. W. DUNN AND B. W. CRISS.MAN
437
~\ I iTTTTTTi ifi ii 11 i i
FIG. 15. As in Fig. 9 except for initial conditions P(4,l).
112
-------
438
JOURNAL OF APPLIED .METEOROLOGY
VOLUME 15
113
-------
MAY 1976
F . M . V U K 0 V I C H , J . W . DUNN AND B . W . C R I S S -M A N
439
FIG. 16. Upper: The vertical velocity (cm s"1) distribution at the 300 m level for
the northeast (upper), west (middle) and east (lower) wind case. The initial conditions
were P(4,l). Positive values indicate upward vertical motions.
downdraft or an intensification of a compensating
downdraft associated with the heat island perturba-
tion by the topography. The downstream updraft
center may be associated with the heat island per-
turbation or with a topographical perturbation induced
by flow over the bluffs to the southwest or with a
combination of both effects. The third explanation
appears to be the most plausible. The narrow region
of upward motion oriented northeast-southwest, found
immediately northwest of the downdraft center, is
probably a product of the flow over the bluffs found
to the west.
5. Conclusions
The influence of synoptic wind speed and direction
on the St. Louis heat island have been studied. It is
found that, as the wind speed increases, the intensity
of the heat island circulation decreases, and the heat
plume associated with the heat island extends further
downwind. Thus, the heat island circulation changes
from a system characterized by cyclonic rotation to
one characterized simply by a convergence line. In the
case of wind profile 4, the heat plume and the con-
vergence line associated with the urban heat island
extend farther downstream, which is to be expected.
However, a tongue of positive vertical motion also
extends downstream from the heat island center and
contains a number of centers of positive vertical
motion. The center containing the largest positive
upward motion in the tongue is generally the one
farther downstream from the heat island center. It is
suggested that this may be a result of topographical
effects. In the future, the interaction of topography
and the heat island will be examined much more
thoroughly.
It is found that when the wind speeds are small,
the intensity of the St. Louis heat island is inde-
pendent of wind direction. However, under strong
wind speeds, the intensity of the heat island changes
markedly from one wind direction to the next. The
most dramatic changes occur when the wind direction
is from the west, east and southeast, in comparison
to those cases when the wind direction is from the
northeast or north. In all cases, the changes are at-
tributed to the influence of local topography on the
urban heat island. These effects do not appear until
the initial wind speed near the surface reaches a
critical value which apparently is greater than or
equal to 3 m s""1.
114
-------
440
JOURNAL OF APPLIED METEOROLOGY
VOLUME 15
Acknowledgment. This research was supported by
NSF/RANN Grant GI-34345.
REFERENCES
Ackerman, B., 1972: Winds in the Ekman layer over St. Louis.
Preprints Conf. Urban Environment, Philadelphia, Amer
Meteor. Soc., 22-28.
Angell, J. K., W. H. Hoecker, C. R. Dickson and D. H. Pack,
1973: Urban influence on a strong daytime air flow as de-
termined from tetroon flights. /. Appl. Meteor., 12, 924-936.
Bomstein, R. D., 1972: Two-dimensional non-steady numerical
simulation of nightime flow of a stable planetary boundary
layer over a rough warm city. Proc. Conf. Urban Environment,
Philadelphia, Amer. Meteor. Soc., 89-94.
Chandler, T. J., 1964: City growth and urban climates. Weather,
19, 170-171.
, 1965: The Climate of London. Hutchinson and Co., 250 pp.
, 1967: London's heat island. Biometeorology, Part n. Proc.
Third Intern. Biometeorological Congress, Pau, France, 1-7
September 1963, Oxford, Pergamon Press, 589-597.
Clarice, J. F., 1969: Nocturnal urban boundary layer over Cin-
cinnati, Ohio. Man. Wea. Ra., 97, 582-589.
Delage, Y., and T. A. Taylor, 1970: Numerical studies of heat
island,circulations. Boundary-Layer Meteor., 1, 201-226.
Landsberg, H. E., 1956: The climate of towns. Proe. Intern.
Symp. Man's Role in Changing the Face of the Earth. The
University of Chicago Press, 584-606.
Lettau, H., 1969: Note on aerodynamic roughness-parameter
estimation on the basis of roughness-element description.
/. Appl. Meteor., 8, 822-832.
Lowry, W. P., 1967: The climate of cities. Sci. Amer., 217, Xo. 2,
15-32.
Mitchell, J. M., Jr., 1961: The thermal climate of cities. Air over
Cities, U. S. Public Health Service Rep. A62-5, Cincinnati,
Ohio, 131-145.
Stanford University Aerosol Laboratory and the Ralph M.
Parson's Company, 1953: Behavior of aerosols within cities.
Joint Quarterly Reports, Nos. 3, 4, 5 and 6.
Vukovich, F. M., 1971: Theoretical analysis of the effect of mean
wind and stability on a heat island circulation characteristic
of an urban complex. Man. Wea. Ra., 99, 919-926.
, 1973: A study of the atmospheric response due to a diurnal
heating function characteristic of an urban complex, if on.
Wea. Ra., 101, 467-474.
, 1975: A study of the effect of wind shear on a heat island
circulation characteristic of an urban complex, if on. Wea.
Ra., 103, 27-33.
Woollura, C. A., 1964: Notes from a study of the micrometeorology
of the Washington, D. C. area for the winter and spring sea-
sons. Weatherwise, 17, 263-271.
115
-------
APPENDIX B
AREA SOURCES WITH OBJECTIVE VARIATIONAL ANALYSIS MODEL
116
-------
AREA SOURCES WITH OBJECTIVE VARIATIONAL ANALYSIS MODEL
During the development of the model, an estimate of the area emissions
for the test grid was needed. Rather than use an existing emissions
inventory, a binormal, Gaussian distribution centered at a prescribed x ,y
with a total emission, Q. and standard deviation q and q was chosen
Thus, at every point (x,y) the emission was given by:
In this formulation the amount of material emitted, the location of sources
and the shape of the distribution can be specified with great versatility.
For example, a line source in the y direction can be approximated by
choosing q > 10 Ax, and q < Ax/2; a point .source by q < Ax/2, q < Ax/2;
y x y x
or large nearly uniform area sources (q ,q > 5Ax) .
y x
RELATIVE DISPERSION OVERLAY
Rather than adopt the approach that area sources can be treated as
virtual point sources, such as in the models of the EPA's UNAMAP Library,
the concentration at a receptor downwind of an area source was computed as
the integrated concentration arising from a set of smaller areas of
uniform emissions.
Referring to Figure B-l, the concentration at x,y due to an emitting
area centered at x ,y is the integral of the subareas of the larger area.
s*. **
A representative area is located at (x,y). Assuming a steady-state
Gaussian point source at that point, the concentration at (x,y), dKx,y),
is given by
.
2irUa (x-x) a (x-x)
z y
exp
[/ ~ \*n
i/lo—)]
\a (x-x)/ J
Letting x = x + v, y = y + u and integrating gives
u
Ax/ 2
-Ax/ 2
\s
Q(xQ,yo)
ir U a (x-x -v)
dv
117
-------
00
t
Y1 A
y
\j
Vo
-— — — -
•9
..I-*
* T
X X,
x-
Figure B-l. Coordinates^for determining the contribution of emissions from a
subarea (x,y) within the finite difference grid (x ,y ) upon a
receptor at a general grid point (x,y). ° °
-------
Evaluating the interior (bracketed) integral, B, gives,
B = /2iT erf (p(v))
when
i (y-y -Ay/2)
P(v) =
/I
The complete integral is
Ax/2
"i
erf
* U J ° (x-x -v)
-Ax/2 z °
which is evaluated numerically for specific locations,
m
x = (XQ - Ax/2) + JIA x A = 1,2,3,...,!
y = (yQ - Ax/2) + mA x m = 1,2,...,5
For a given atmospheric stability condition, a i|i(&,m) can be deter-
mined. The total concentration at a given point in the integral of all
of the contributions from upwind sources which can be rapidly calculated
using the appropriate (&,m) array, emissions and wind.
When the actual wind is not oriented parallel to the x direction, the
same array can be oriented along the direction of the v component and that
contribution can be computed.----- If the wind component changes sign between the
source and the receptor, it is assumed that the source did not affect
that receptor.
ATMOSPHERIC STABILITY
In practice,for models based on the Gaussian distribution, atmospheric
stability is accounted for in the magnitude of a and az, the cross wind,
and vertical plume spread coefficients. These quantities are expressed,
in general, as functions of distance downwind from a source, x, by
0y = al + a2xb
az = Cl + C2xd '
119
-------
Martin's (1976) values of a.., a*, b, c.., c« and d for the six Pasquill
Turner stability categories can be used. If the relationship of a and
a to the eddy diffusion coefficient is used, then
z
and
a 2 = 2K t
y y
a 2 = 2K t
z z
where t = x/U is the travel time from source to receptor.
120
-------
APPENDIX C
SENSITIVITY ANALYSIS AND TESTING OF THE
VARIATIONAL OBJECTIVE ANALYSIS MODEL
121
-------
SENSITIVITY ANALYSIS AND TESTING OF THE
VARIATIONAL OBJECTIVE ANALYSIS MODEL
CASE I - THE ROLE OF a/a
The most favorable conditions for obtaining a solution to the analysis
in the idealized case occur when errorless observations are made at every
grid point, and they are known to satisfy an analytical distribution. In such
a case, the solution analysis should be the analytical solution. Departure,
therefrom, might indicate weaknesses in the formulation of the analysis
technique, the numerical techniques used, or the differences in the solution
of the constraint equation and the analytical concentration.
Accordingly, choosing 5/a = 10 everywhere should make 41 = i|» = ip*. Using
the input variables listed in Table C-l, this solution was confirmed.
Decreasing the ratio of 1.0, and then 0.1 progressively, increases the
influence of the constraint equation upon the solution. As shown in
Figure C-l, ij/ continually departs from y. The maximum value decreases and is
displaced further downstream, the upwind concentrations decrease, and the
downwind concentrations increase accordingly. The total mass is approximately
conserved, but it is redistributed by the analysis. These results show that
the constraint equation does influence the solution away from *f = 41*- That
solutions were obtained as the i|>/ifi* ratio changes suggest that a solution
can be obtained when the ratio is variable; this is verified later.
CASE II - THE ROLE OF THE CONSTRAINT EQUATION
With the assurance of a convergent solution, the role of the ocnstraint
equation was examined by letting a/a = 0.01; this ratio minimizes the
contribution of the observations to the solution and emphasizes the role of
the constraint equation. In effect, this case investigates the solution which
could be obtained if there were no observations and only emissions, wind and
stability were known. Test parameters are given in Table C-l.
The solution, ij;, along y = 5 km is shown in Figure C-2 with the analytic
concentration, i|>*. Clearly, i|> underestimates i|»* as much as 100 percent
(10.8 versus 5.2 @ x = 7) when x < 10 km and overestimates (10.5 versus 8.6
at x = 13 km) along the remainder of the axis. The peak value of ip is about
1 km downwind of the maximum of ij»*. The integrated deficit of (ip-ip*) is 1.08
122
-------
I—I—1—I—I—1—I—I—\—I—I
1—I—i—I—1—I
» I 1 I I I I 1—I—I—1—I—I—I—1
X, km
Figure C-l. Concentration distributions with x^for Case I,
showing i|»* (— • — •) and i|> when a/ a = 10
(— • — •)> a/a = 1.0 ( ) and ci'ci = 0.1
( v •
123
-------
18
16
14
•o 12
10
8
o
c
0) _
o 6
o
o
4
2
I r i i r
i r
i i i i i l
i i i i i i i
i i i i i
i i i i t
10
X ,km
15
20
Figure C-2. Concentration distributions with x for Case II showing
^*(— . — .) and ^ when 1 = 1 ]asit ( ) and L = 4 or
IQ km ( ) for Case II.
124
-------
TABLE C-l. PRINCIPAL INPUT VARIABLES AND CHARACTERISTIC NUMBERS
CASE
la
b
c
Ha
b
c
Ilia
b
IVa
b
c
V
VI
VII
L
(km)
1
1
4
10
1
1
4
4
4
Z "a/a q
(m) (km)
50 10 1.25
1
0.1
50 0.01 1.25
50 0.01 2.50
5.00
31.5 0.01 1.25
100
200
50 1.0 1.25
50 1.0 1.25
50 1.0 1.25
BC» y
CV 35.0
CV 35.0
CV/CG 35 . 0
CV/CG
CV 27.0
39.4
39.5
CV 35.0
CV 35.0
CV 35.0
n
28.5
28.5
114.0
285.0
28.5
37.1
25.4
25.3
114.0
114.0
114.0
£
max
.19
.19
.76
1.89
.16
.11
.30
.05
.01
.76
.76
.76
COMMENT
Test for Solutions
Basic Experiment
Cases II-IV
Solutions Dominated
by the Constraint
Equations
Observations Errors
Observations Removed
See text for oVa
Observations Removed
and Perturbed
o
>VBC - Boundary Conditions, CC - Constant Gradient, CV - Constant Value
Cases VI and VII are discussed in the text (Pages 77 to 84).
-------
times the integrated surplus, indicating that the total mass density of
pollutant was nearly conserved i-n this plane.
The reasons for the solution profile became apparent after examining
Eq. (31) which is appropriate for small a/a , i.e.,
_.
- 5 *- n 3- nCQ
ox
For 5 < 0.2 and the parameters of this test, the £*ty term is much smaller
than the remaining terms and the solution is dominatedby the curvature term
of fy. For the source configuration used, 3 = ~ ( — 2/Q» so the e 0, the curvature of ij; is positive and ty is
concave upward. In this case, (x-x + C <1 ) < 0 for x < 7.1 km. Beyond
o x
that point, the curvature remains negative for as long as 0 makes a contri-
bution; e.g., x-x > 3 q (x > 12.5 km). After that point, Q is very small,
" 3*
the curvature of i|> is near zero. As evident in the solution, the slope of y>
is nearly constant for x > 13 km.
When the curvature is a minimum, the function should be a maximum. The
curvature is an extreme when
which, with the Gaussian source distribution, constant r\ and C,occurs when
For small values of £ (~0.l) and small values ofq (<2), maxima and minima
2C
of ij; should occur when
126
-------
x - X0 « ± qx . (c-i)
The numerical solution of ty agrees with this result near the maximum concen-
tration. The disagreement about the location of the expected minimum concen-
of \l>, i.e., a minimum @ x~xo = ~ and n in the development.
Near x = 3.5, the gradients of 5 and n are large (Figure 14) and should be
considered.
The magnitude of £ can be increased by increasing L. However, this
proportionally increases the coefficient of all of the terms of the finite
difference equation so that the solution is generally unaffected. Test cases
with L = 4 km and 10 km showed only a slight difference in the numerical
solution (Figure C-2) .
The convergence criteria was met after 67 iterations in the
L = 1 km case. After 90 iterations, the solution criteria was not achieved
at 26 points when L = 4 km. At that point, the constraint field, A, was
approximately 10 times larger, residuals were 10 times larger, fl was 100 times
larger so the changes of ty were a tenth of those made in the L = 1 km case.
Thus, convergence was slower. The trend toward convergence left no doubt that
a solution would have been attained with further iterations. A different,
faster relaxation procedure would have helped reduce the number of iterations
needed to obtain a solution.
CASE III - SOURCE DISTRIBUTION/BOUNDARY CONDITIONS
The analysis of Case II is continued by examining the influence of the
source shape upon the distribution. The total emission was held constant but
the distribution parameters q and q were increased to 2.5 km and then to 5.0 km.
The analytic concentration \p* and the solution concentrations, 1(1, along the
y = 5 km axis when a/ a =0.01 are shown in Figure 29 for both constant value
(6i|» = 0) and constant gradient [ 6 (-) = 0] boundary conditions. The solution
9n
indicated by Sif> = 0 in Figure C-3a was discussed in Case II. Changing the
boundary conditions changes the solutions, especially in areas where the
emissions are very small. Upwind of the source, the effects of "background"
127
-------
Figure C-3.
Concentration distribution of iji* (— • —) and i|>
with constant boundary conditions ( ) and
constant flux boundary condition ( ) for
the Case III with a) qx = 1.25 km, b) q =
2.50 km and c) q = 5.00 km. X
128
-------
are lost to the analysis. Downwind the constraint equation overestimates the
concentration by a large amount. Between XQ - 2qx and XQ + 2qx, distribution
of emissions rather than the boundary conditions control the solution.
When the emissions are more uniformly distributed q = q = 2.5 kin,
x y
i|» (Figure C-3b) becomes more like i|>*, and the boundary conditions have less
effect upon the solution. The maximum ty is displaced about 1.5 km further
downwind of the maximum i|/* (~1.5 km) as Eq. (C-l) predicts, but the maximum
values differ only slightly. Further spreading of the emissions (q = q =
x y
5.0) further brings the analytic and constraint solutions closer together.
The boundary conditions showed no substantial effect on the solution since
emissions were effective at the upwind and downwind boundary of the grid.
Clearly, the constraint equation gives a better reproduction of the actual
distribution of concentration when the source is uniformly distributed,
regardless of the boundary conditions.
CASE IV - EFFECTS OF Z
In the development of the constraint equation, y(x) was nearly maximized
by choosing Z ~ /2~ Sz. The choice of Z affects the values of y» n and £ and
hence, the solution. Furthermore, the background concentration of ty at Z is
arbitrarily given as 0.999 times the background at z = 0 to give a finite
value of S at points upwind of the sources. Increasing Z effectively makes
the background even more uniform with height. The effects of assuming
different values of Z on the constraint equation solution (a/a = 0.01) were
investigated numerically. The experimental data are given in Table C-l.
Although the convergence criteria had not been met after 90 iterations, when
Z was increased, the distribution of ty along the y = 5 km axis shown in
Figure C-4 indicate the trend of the solution. As Z increases, the solution
progressively underestimates i|>*, and the maximum <|» is displaced progressively
further downwind.
A reduction of Z (z = S (Ax)) reproduces the magnitude and location of
the maximum i|>* quite well but overestimates the downwind concentrations. The
error variance, initially 36.6, is reduced to 13.5 in contrast to the Z = 50 m
case when the error variance begins at 46.0 and reduces to 8.3. By that
measure, the overall solution has not been improved using the smaller Z.
129
-------
18
14
ro
>> 12
10
£ 8
(U
o
I 6
4
i i i i i n?Fi;i i i i i i i i i r
• \*«
i T—t.^
I I I I I I till
10
15
20
x, km.
Figure C-4. Concentration distributions with x for Case IV,
showing ty* (— • — •), and ty when Z = 200 in
( ), Z = 100 m ( ), Z = 50 m (x
and Z = 31.5 m ( ).
130
-------
The solutions for variable Z clearly show that the choice of Z and
the background value of Z strongly influence the resulting concentrations.
This does not present a major problem when the impact of area sources upon
ground level air quality is considered. However, elevated sources will
require careful consideration and possibly some revision of the analysis
approach.
CASE V - OBSERVATIONAL ERRORS
A certain amount of error exists in each measurement of a pollutant
concentration. In an effort to simulate a more realistic situation and
the performance of the analysis technique in that situation, the observations,
ij;, were redefined so that at all interior points,
*i,j = *i,j + ^ (0'V
where E.. is a number drawn from a normal (Gaussian) population having a zero
mean and standard deviation of iL,. All other parameters were the same as
D
the Case I(a) of Table C-l.
The t|/*, ijj and i|) are shown in Figure C-5. The perturbations of ^ are
reduced by the analysis. This means that the analysis filters some high
frequency noise of the observations. The trend is toward the solution of
Case I(a). The experiment was performed for different sets of E.., with
similar results.
131
-------
I I I I I I I I I I I I I I I I
Figure C-5. Concentration distributions with x for Case V showing
$(x x), ij>*(— • — •), and i|» ( ) along
y = 5 km.
132
-------
REFERENCES
Barnes, S. L. , 1964: "A Technique for Maximizing Details in Numerical Weather
Map Analysis", J. Appl. Meteor., _3, 396-409.
Busse, A. D. and J. R. Zimmerman, 1973: "User's Guide to the Climatological
Dispersion Model", EPA-R4-73-024, (NTIS Accession Number PB-227-346).
Cressman, G. P., 1959: "An Operational Objective Analysis Systems", Mon. Wea.
Rev., _8_7, 367-374.
Endlich, R. M. and R. L. Mancuso, 1968: "Objective Analysis of Environmental
Conditions Associated with Severe Thunderstorms and Tornadoes", Mon.
Wea. Rev., 96_, 342-350.
Gifford, F. A., Jr. and S. R. Hanna, 1973: "Modeling Urban Air Pollution",
Atmos. Environ., ]_•> 131-136.
Holzman, B., 1943: "The Influence of Stability on Evaporation", Ann. N. Y.
Acad. Sci., 44, 13-18.
Lanczos, C., 1970: The Variational Principles of Mechanics, Fourth Edition,
University of Toronto Press, Toronto.
Lettau, H., 1969: "Note on Aerodynamic Roughness-Parameter Estimation on the
Basis of Roughness-Element Description", J. Appl. Meteor., 8_, 822-832.
Martin, D. 0., 1976: "The Change of Concentration Standard Deviations with
Distance", J. Air Poll. Contr. Asso., 26, 145.
Rossby, C. G. and R. B. Montgomery, 1935: "The Layer of Frictional Influence
in Wind and Ocean Currents", Papers Phys. Oceanog. Meteorol., ^, No. 3.
Sasaki, Y., 1958: "An Objective Analysis Based on the Variational Method",
J. Meteor. Soc., Japan, 36, 77-88.
, 1970a: "Some Basic Formalisms in Numerical Variational Analysis",
Mon. Wea. Rev., 9£, 875-883.
, 1970b: "Numerical Variational Analysis Formulated Under the
Constraints as Determined by Long-Wave Equations and a Low-Pass Filter",
Mon. Wea. Rev., j>8., 884-898.
, 1970c: "Numerical Variational Analysis with Weak Constraint and
Application to Surface Analysis of Severe Storm Gusts", Mon. Wea. Rev.
^8., 899-910.
133
-------
REFERENCES (Cont.)
, and J. M. Lewis, 1970: "Numerical Variational Objective Analysis
of the Planetary Boundary Layer in Conjunction with Squall Line Formation",
J. Meteor. Soc. of Japan, 48, 381-398.
Stanford University Aerosol Laboratory and the Ralph M. Parson's Company, 1953:
"Behavior of Aerosols within Cities", J. Quar. Reports, Nos. 3, 4, 5 and 6.
Stephens, J. J., 1970: "Variational Initialization with the Balance Equation",
J. Appl. Meteor., 9_, 732-739.
Sweeney, H., 1969: "Development of Sampling Guidelines for Regional and Local
Planning", Final Report, NAPCA, Contract No. CPA-22-69-57.
Thompson, P. D., 1961: Numerical Weather Analysis and Prediction, The
MacMillan Company, New York.
TRW Systems Group, 1969: "Air Quality Displace Model, USDHEW, Public Health
Service, NAPCA, Washington, D. C.
Turner, D. B., 1969: "Workbook of Atmospheric Dispersion Estimates", U.S.D.
HEW, PHS, Publication No. 999-AP-26, Cincinnati, Ohio.
User's Guide to the Statistical Analysis System, 1972: Jolayne Service,
North Carolina State University.
Wilkins, E. M., 1971: "Variational Principle Applied to Numerical Objective
Analysis of Urban Air Pollution Distribution", J. Appl. Meteor., 10;
974-981.
, 1972: "Variationally Optimized Numerical Analysis Equations
for Urban Air Pollution Monitoring Networks", J. Appl. Meteor., 11;
1334-1341.
134
*U.S. GOVERNMENT PRINTING OFFICE: 1978 - 785-926/1219 Region No. 9-1
-------
TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing)
. REPORT NO.
EPA-600/4-78-030
3. RECIPIENT'S ACCESSION NO.
OPTIMUM METEOROLOGICAL AND AIR POLLU-
4. TITLE AND SUBTITLE
'ION SAMPLING NETWORK SELECTION IN CITIES
Volume I: Theory and Design For St. Louis
5. REPORT DATE
June 1978
6. PERFORMING ORGANIZATION CODE
. AUTHOR(S)
'red M. Vukovich, Walter D. Bach, Jr., and
C. Andrew Clayton
8. PERFORMING ORGANIZATION REPORT NO.
9. PERFORMING ORGANIZATION NAME AND ADDRESS
Research Triangle Institute
P.O. Box 12094
Research Triangle Park, North Carolina 27709
10. PROGRAM ELEMENT NO.
1HD620
11. CONTRACT/GRANT NO.
68-03-2187
12. SPONSORING AGENCY NAME AND ADDRESS
U.S. Environmental Protection Agency-Las Vegas, Nevada
Office of Research and Development
Environmental Monitoring and Support Laboratory
Las Vegas, NV 89114
13. TYPE OF REPORT AND PERIOD COVERED
14. SPONSORING AGENCY CODE
EPA/600/07
15. SUPPLEMENTARY NOTES
The report is the first in a series on this topic. For further information contact
James L. McElroy, Project Officer (702)736-2969. ext. 241 in Las Veeas, NV
16. ABSTRACT
A technique was developed to establish an optimum meteorological and air pol-
lution sampling network in urban areas. The basis of the network is the wind field in
the urban area rather than the air pollution distribution because it provided a solu-
tion with longer-term stability than the air pollution distribution.
Three specific models are required in order to determine the optimum network. These
are: a three-dimensional hydrodynamic model; a statistical model; and an objective
variational analysis model. The primitive equation model is used to simulate the wind
field for a variety of cases. These simulated data were used to determine the form of
a regression model which approximates the various wind fields. A regression model form
was then used, along with a set of potential network sites, and a criterion for judging
alternative networks to derive the sampling network for the winds. The method used to
develop the network involved the successive elimination of candidate sites until a
reasonably sized network was achieved. The air pollution distribution is obtained
through an objective variational analysis model. The model simultaneously minimizes
the error variance by comparing observed pollution concentrations with derived pollutioi
concentrations and the error variance of the constraint equation.
7.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.lDENTIFIERS/OPEN ENDED TERMS
c. COSATl Field/Group
Air quality
Monitoring
Air monitoring network
design
Ambient air monitoring
Carbon monoxide network
design
Carbon monoxide monitorir
network design method
t.
04B, 09B, 12A
18. DISTRIBUTION STATEMENT
RELEASE TO PUBLIC
19. SECURITY CLASS (This Report)
UNCLASSIFIED
. NO. OF PAGES
154
20. SECURITY CLASS (Thispage)
UNCLASSIFIED
22. PRICE
EPA Form 2220-1 (Rev. 4-77) PREVIOUS EDITION is OBSOLETE
------- |