EPA-600/3-83-035
May 1983
A REGIONAL SCALE (1000 KM) MODEL
OF PHOTOCHEMICAL AIR POLLUTION
Part 1. Theoretical Formulation
by
Robert G. Lamb
Meteorology and Assessment Division
Environmental Sciences Research Laboratory
Research Triangle Park, North Carolina 27711
ENVIRONMENTAL SCIENCES RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
RESEARCH TRIANGLE PARK, NORTH CAROLINA 27711
-------
EPA-600/3-83-035
May 1983
A REGIONAL SCALE (1000 KM) MODE
OF PHOTOCHEMICAL AIR POLLUTION
Part 1. Theoretical Formulation
by
Robert G. Lamb
Meteorology and Assessment Division
Environmental Sciences Research Laboratory
Research Triangle Park, North Carolina 27711
ENVIRONMENTAL SCIENCES RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
RESEARCH TRIANGLE PARK, NORTH CAROLINA 27711
-------
DISCLAIMER
This report has been reviewed by the Environmental Sciences Research
Laboratory, U.S. Environmental Protection Agency, and approved for publi-
cation. Mention of trade names or commercial products does not constitute
endorsement or recommendation for use.
AUTHORS' AFFILIATION
The author is on assignment with the U.S. Environmental Protection
Agency from the National Oceanic and Atmospheric Administration, U.S.
Department of Commerce.
ii
-------
PREFACE
The original goal of this work was to develop a specific model of
regional scale photochemical air pollution. However, as the work pro-
gressed and new developments and ideas continually emerged, the need was
seen for a general modeling framework within which the various physical
and chemical processes that play important roles could be treated in
modular form. This would permit ongoing incorporation into the model of
state-of-the-art techniques without the need to overhaul the model each
time. The basic framework we have developed for this purpose is based on
phenomenological concepts and is presented in Section 2. Most of the sub-
sequent sections develop specific modules for use in implementing a first
generation model. Therefore, the techniques presented in those sections
should not be viewed as an integral part of the model but rather as one of
many possible methods that could be utilized to perform specific functions.
Moreover, in a few instances the techniques presented require further deve-
lopment before they can be employed in an operational role. This is par-
ticularly true of the mixed layer submodel given in Section 4.1. In short,
the reader will find in this report rather exact specifications of the
framework of a regional scale model but only rough sketches of some of the
components that are required to make the model work. In later parts of
this report we will describe the components that we assemble to construct
an operational, first generation regional scale model of photochemical air
pollution.
iii
-------
As we approach the task of simulating air pollution chemistry and dis-
persion over multi day 1000 km scale domains, we move into a realm of prob-
lems that lies largely beyond the scope of empirical science. The classical
definition of turbulence and the empirical data that form the basis of
current short range diffusion models are insufficient to treat dispersion
and chemistry at long range. As a consequence we must turn to theoretical
science to provide the additional information we need and to formulate ex-
pressions for the quantities we wish to predict. With this shift from em-
piricism to theory there arises the somewhat philosophical question of
whether the limitations of theoretical science preclude the formulation of a
model that can provide the information that regulatory officials require.
One way of phrasing this question is this: given perfect emissions data,
perfect knowledge of the chemical kinetics, an exact solution to the turbu-
lence closure problem, and an infinitely large and fast computer, could one
predict within the limits of accuracy required the quantities needed in
regulatory studies? In Sections 6 and 7 we discuss this important question
in detail, and we lay the basis for an effort parallel to the model develop-
ment work that in principle can provide quantitative estimates of the uncer-
tainty, or "error-bounds", associated with the predictions of regional scale
models. In the author's view this task is crucial to the meaningful utili-
zation of model predictions in decision making processes.
R. G. Lamb
February 1982
IV
-------
ABSTRACT
A theoretical framework for a multi-day 1000-km scale simulation model
of photochemical oxidant is developed. It is structured in a highly modular
form so that eventually the model can be applied through straightforward
modifications to simulations of particulates, visibility and acid rain.
The model structure is based on phenomenological concepts and consists
of three and one-half layers. The interface surfaces separating the layers
are functions of both space and time that respond to variations in the
meteorological phenomena that each layer is intended to treat. Among the
physical and chemical processes affecting passage and distribution of photo-
chemical concentrations that the model is designed to handle are: horizontal
transport; photochemistry; nighttime wind shear and the nocturnal jet; cumulus
cloud effects; mesoscale vertical motion; mesoscale eddy effects; terrain
effects; subgrid scale chemistry processes; natural sources of hydrocarbons,
NO , and stratospheric ozone; and wet and dry removal processes, e.g., washout
A
and deposition.
The predictibility of pollutant concentrations at long range is considered
along with such related problems as the parameterization of "mesoscale" diffu-
sion and the design of model "validation" experiments. A basis is established
for estimating quantitatively the levels of uncertainty associated with disper-
sion model predictions.
This report focuses on theoretical aspects of the model and the question
of predictability. Results of the model's performance and quantitative assess-
ments of predictability will be presented in subsequent parts of this report.
-------
CONTENTS
Preface H1
Abstract v
Figures viii
Tables x
Acknowledgments xi
1. Introduction 1
2. Derivation of Prognostic Equations for
Layers 1, 2, and 3 18
3. The Interface Surfaces H. 26
The model's top surface H3 26
The subsidence inversion oase surface, H~ 26
The nocturnal jet top, H, 27
The surface layer top, H^ 32
4. The Interface Fluxes: F.. . . 35
Flux across the subsidence inversion base: F? 36
Equations governing w2, fe/A6, and z~ 42
Flux across the jet layer top: F.. 53
Flux across the surface layer top: F 62
5. Derivation of Diagnostic Equations for Layer 0 65
Mass balance equations in V - 5V 71
The flux component F in the case of no
cumulus clouds . 75
Modifications of F to include effects of
cumulus clouds 77
Equations for the cell-averaged and
fluctuation concentrations in layer 0 80
6. The Representation of Atmospheric Motion:
Winds and Turbulence 83
A game analogy of diffusion modeling 88
The statistical bases for air pollution modeling ... 99
The global family of atmospheric states TOO
General considerations 100
The constraints of physical laws 103
7. Mesoscale Diffusion; Concentration Predictability;
Model "Validation"; and Related Topics 114
Remarks 125
The character of diffusion at long range '25
Parameterization of dispersion ]32
Predictability of concentration at long range . . '33
Model validation exercises 136
Predicting impact of future emissions
distributions 137
vi
-------
8. Formulation of Miscellaneous Fields and
Model Parameters 139
Source emissions rates i, S 140
The plume volume fraction, c 141
Turbulence parameters w , x , etc on H
and in Layer 0. ...".." ? 142
The plume entrainment velocity, v 144
9. Computer Solution of the Governing Equations 146
Introduction 146
Transformation of the governing equation . ..... 148
Solution of the r equation, (9-5) ... 154
Solution of the y equation, (9-17) 180
Solution of the Y/a)n equation, (9-24) 189
10. Summary of Model Equations 199
References 207
Appendices
A. Transformation of the governing equations
to curvilinear coordinates 210
B. Criteria for validity of the steady-state
assumption in layer 0 217
C. Explicit forms of the b and g matrix elements
b , g that enter in the ? equation 221
vn
-------
FIGURES
Number Page
1-la. Schematic illustration of the "dynamic" laye/ structure
of the regional scale model and the daytime phenomena
each layer is-designed to treat 7
1-lb. Same as 1-la except nighttime phenomena 8
1-2. Idealized view of the fate of polluted air parcels
in a cumulus humilis topped mixed layer. A, parcel
receives emissions; B, it rises and cools; C, saturation
occurs producing cloud droplets; 0, the parcel descends
and dries; and E, it receives new emissions and/or
deposits pollutants on surfaces. From E the cycle
begins anew 11
1-3. Hypothetical case of an urban area in which 03 and
NO are segregated into 16 cells. If this distribution
is simulated in a regional model whose cell size is
the large outer square, subgrid chemistry phenomena
arise 14
4-1. Schematic vies of cumulus clouds transporting material
aloft from the mixed layer and acting to lower the
surface Hg 35
4-2. Projection of a unit horizontal square on the surface
H2 38
4-3. Probability density of vertical velocity w~ near the
mixed layer top (H2) derived from the numerical model
of Deardorff (1974) 45
4-4. Comparison of the true solution of Equation (4-27a)
with the approximation (4-27c) 47
4-5. Schematic illustration of typical temperature, wind
speed and mixing ratio profiles over (a) rural and
(b) urban areas at night 54
4-6. Intersection of the surface z = z,, defined by
Equation (3-5), with elevated terrain features.
A fraction a-p of the grid cell's horizontal area
A is intersected by terrain 57
vi i i
-------
4-7. Coincidence of surfaces H and H, with terrain
surface. Surface area of H lying on terrain is
approximated as (ay - ) that comprise the global family of fluid
flow velocity Fourier Transforms. Set H represents
functions that satisfy purely mathematical constraints
like Equation (6-15). Set A represents functions
that satisfy mass continuity relationships, like
Equation (6-18); and set ir denotes functions that
satisfy all other laws like momentum and entropy.
Intersection of A and ir defines the global family
n of u. 106
7-1. Mapping of the set n and W in flow field function
space into the setT and C in concentration function
space for a given source distribution S 117
7-2. Plot, of a hypothetical example of the function
(x0>t|xs,t') for given source location xg and
receptor site x . The function $ has unit value
inside the shaded regions and is zero everywhere
else 121
7-3. Trajectories of particles eminating from x in the
ff V
ensemble W associated with the single wind obser-
vation u(x ,t) = u , where x = x 127
7-4. Plumes from the C ensemble that mark the limits
of lateral motion 134
9-1. Illustration of the grid, coordinate system, and
other parameters on which the 2-D representation
(9-46) and (9-47) of the function r(n,0 is based. ... 163
IX
-------
9-2a. Results after 50 time steps from a simulation by
the biquadratic scheme of the advection of a conic
shaped cloud in a rotating flow. Angular speed, of
rotation. Angular speed of rotation o> = (40A)
and 1 = 0 ........................ 167
9-2b. Same as 9-2a, except results are after 100 time steps. . . 168
9-3a. Same as Figure 9-2a using the bicubic scheme ....... 169
9-3b. Same as 9-3a, except results are after 100 time steps. . . 170
9-4. Comparison of simulations by 3 differencing schemes
of the advection of an ellipsoidal cloud in a rotating
flow. Panels a-e display different cross -sections of
the cloud (indicated in the upper right corner of each
panel) after 1 complete rotation of the cloud, 100
time steps in the case of schemes Q and S, 150 steps
in the case of Z .......... ............ 173
9-5. Sample solutions of the y equation (9-78) ......... 187
9-6. Temporal behavior of the reaction system (9-24) for
initial species concentrations: (a]cNO = 100 ppb,
N0« = 500 ppb, 02 = 30 ppb, 0 = 10 ° ppb; and (b)
one-tenth the vafues used in (a). Arrows indicate
the time steps used by the model to arrive at
concentrations at t = 100 s ................ 195
A-l. The curvilinear coordinate system used in the model. ... 210
B-l. Two-layer system for analysis of the steady-state
mixing assumption ..................... 217
TABLES
Number Page
9-1. Computer time requirements for the Q, S and
Z schemes 179
-------
ACKNOWLEDGMENTS
The model development work reported here was initiated in 1977 by
Dr. Kenneth Denierjian through a contract with Systems Applications, Inc.
Without his support in the intervening years, this work would not have
been possible.
Dr. Robert Papetti reviewed the original derivation of the model
equations presented in Section 2 and provided a number of suggestions
which lead to a greatly improved set of notation and also to a more elegant
formulation of some of the basic mathematical relationships. I am indebted
to him for this help and for his interest in this work generally.
In the course of this project I have benefited by discussions on many
topics with my colleague, Dr. Francis Binkowski. I would also like to
acknowledge the capable assistance of Mr. Chris Mitchell in developing and
testing the new differencing schemes presented in Section 9.
Finally, I extend my sincere appreciation to Ms. Barbara Hinton for
her patience and stamina in typing and retyping this manuscript and its
innumerable equations. My thanks go also to Ms. Sheila Pinkney for her
noble efforts in typing the original draft.
xi
-------
SECTION 1
INTRODUCTION
The objective of this study is to develop a model that can guide the
formulation of regional emissions control strategies. In this task the model
will be called upon to estimate the impact of sources on concentrations in
remote regions, to determine the pollution burden that cities impose on
distant neighbors, and eventually to analyze quantitatively emissions impacts
on acid rain, visibility and fine particulates. We believe that in all these
roles the utility and credibility of the model will be determined primarily
by the extent to which it accounts for all the governing physical and chemical
processes. Accordingly, in this report we will formulate a generalized
model that in principle can treat all of the chemical and physical processes
that are known, or presently thought, to affect the concentrations of air
pollutants over several day/1000 km scale domains. Among these processes
are (not necessarily in order of importance):
1. Horizontal transport;
2. Photochemistry, including the very slow reactions;
3. Nighttime chemistry of the products and precursors of
photochemical reactions;
4. Nighttime wind shear, stability stratification, and turbulence
"episodes" associated with the nocturnal jet;
1
-------
5. Cumulus cloud effects—venting pollutants from the mixed
layer, perturbing photochemical reaction rates in their shadows,
providing sites for liquid phase reactions, influencing changes
in the mixed layer depth, perturbing horizontal flow;
6. Mesoscale vertical motion induced by terrain and horizontal
divergence of the large scale flow;
7. Mesoscale eddy effects on urban plume trajectories and growth
rates;
8. Terrain effects, on horizontal flows, removal, diffusion;
9. Subgrid scale chemistry processes—resulting from emissions
from sources smaller than the model's grid can resolve;
10. Natural sources of hydrocarbons, oxides of nitrogen (NO ) and
A
stratospheric ozone (03);
11. Wet and dry removal processes (e.g., washout and deposition).
Of the 11 processes listed above, only the first and last have been
treated in any detail in the regional scale models of air pollution developed
to date. In fact, a review of these models (see, for example, reviews by
Drake and Bass in Henderson ejt al_. 1980) reveals that virtually all Eulerian
type models are in essence simply expanded urban scale models. They account
for the physical processes that are active during daylight hours and within
10 km or so of a source, but they neglect both the processes that are impor-
tant beyond this distance and those that are active at night.
When this model development work was initiated some 3 years ago, an
attempt was made to derive from the observational evidence available at
-------
that time an estimate of the minimum vertical and horizontal resolutions
necessary to describe regional scale air pollution phenomena. The aim was
to arrive at the best compromise between the restrictions imposed upon the
model by computer time and memory limitations and the need to describe as
accurately as possible all of the governing processes cited earlier. We
reviewed the nitric oxide (NO), Oo and meteorological c -^ported in
Siple et, al_. (1977) by the.participants of the 1975 Northeast Oxidant Trans-
port Study and made the following observations:
o Ozone concentrations at an elevation of 2 km are consistently
between 40 and 60 ppb regardless of the concentrations below
this level. Exceptions occur when the base of the synoptic
subsidence inversion is above 2 km. On these occasions, 03
values at 2 km can reach 80 ppb. This suggests that the top
of the model need not be much higher than 2 km .(" 800 mb).
o After the passage of a cold front and the onset of high pressure
in the study area, dry, strong subsidence inversions based
between 1700 m and 2000 m frequently occur. Ozone values within
these inversions usually appear to be about 50 ppb near the
base and to increase upward. Nitric oxide concentrations within
the inversion are very lean and are approximately constant with
height (see, for example, Siple eta1_. (1977), Flight 12,
Spirals 1 through 3; Flight 6, Spirals 1 through 3). These ob-
servations and the extreme dryness of the inversion air indicate
that this air might be of stratospheric origin. If so, entrain-
ment of inversion layer air into the mixed layer may result in
a contribution of stratospheric 0^ to ground-level oxidant
-------
concentrations. This can be handled in the model by appropriate
material fluxes across the inversion base.
o Nitric oxide is usually uniformly distributed in the mixed layer
(the layer between the ground and the base of the subsidence
inversion). Exceptions consist chiefly of periods when maximum
NO occurs near the ground. Both characteristics can be accom-
modated in a model that represents the mixed layer with only 2
vertical levels.
o-- For the majority of the flights, there is strong evidence of a
positive correlation between 0, concentration and dew point
temperature. The latter is a measure of the concentration
of water vapor (see Flight 3, Spirals 1 through 6; Flight 8,
Spiral 1). This suggests that the bulk of the 03 observed
in the study area is of an anthropogenic origin; if it originated
in the stratosphere, there would be a negative correlation between
Og and dew point. These observations also suggest that estimates
of the rate of transport of 0^ and its precursors through the
top of the mixed layer by large-scale vertical motion and cumulus
clouds can be obtained from estimates of the transport of water
vapor in the upper atmosphere by these same processes. The latter
estimates can be obtained from surface and rawinsonde meteorolo-
gical data.
o When the Boston urban plume spreads out to sea during daylight
hours, 0^ concentrations in the plume are usually at a maximum
near the sea surface, and they decrease steadily upward (see
for example, Flight 2, Spirals 2, 3, 6, 7, 8; Flight 9, Spiral 1).
-------
A plausible explanation of this pattern is that as the warm air
moves over the colder water, the layer of air in contact with
the water is cooled and thereby stabilized. As this occurs,
turbulent motions that otherwise would mix pollutants vertically
are attenuated, and rich mixtures of NO and hydrocarbons remain
A
at low altitudes. With sufficient sunlight, large quantities of
0- can be produced in this air. The same type of stabilization
occurs at night over land as air in contact with the ground is
cooled by long wave radiation, but, in this case, the absence of
sunlight results in a net decrease in 03 concentration. These
phenomena play a crucial role in determining the ultimate species
concentrations in an aged air mass that has spent part of its
life over the sea or over the land at night. To account properly
for these processes, the model must possess the ability to simu-
late the development of these shallow, stable layers adjacent to
the ground and the rates of chemical reactions within these layers.
o Marked changes in wind speed and direction often occur near the
interface between surface stable layers and the air mass above
and at the interface between the daytime mixed layer and the subsi-
dence inversion. The model must have separate layers to treat
each of these regions.
The characteristics of the 03 distribution described above would
require at the very least a 3-level -model: one level assigned to the
surface layer, a second level to the remainder of the daytime mixed layer,
and a third layer atop the mixed layer. The top level would be used in
conjunction with the mixed layer to account for downward fluxes of strato-
-------
spheric 0, as well as upward fluxes of 03 -nd Its precursors into the sub-
sidence inversion layer above. Material that entered this top layer could
be transported by winds aloft to areas outside the modeling region; it
could reenter the mixed layer by subsidence or entrainment; it could enter
precipitating clouds and be rained out of the atmosphere; or it could
undergo chemical transformation. Representing the subsidence inversion,
where cumulus clouds often form under stagnant high pressure conditions, the
top level of the model could be instrumental in simulating the chemical sink
effect of heterogeneous (within cloud droplets) reactions among Og, its
precursors, and other natural and pollutant species. Including cloud effects
in the model would be especially important in simulating sulphur dioxide
(S02) and sulfates.
Having 3 layers in a model is insufficient in itself to simulate the
phenomena we have discussed above. For example, 3 layers of constant thick-
ness are incompatible with the spatial and temporal variability that the
radiation inversion and mixed layer thicknesses are known to have. What is
needed in the model is 3 "dynamic" layers that are free to expand and con-
tract locally in response to changes in the phenomena they are intended
to treat. The model we will develop here possesses this property. Figure
1-1 illustrates the vertical structure of this model and the physical pheno-
mena that each layer is intended to simulate. The surfaces that comprise
the interfaces of adjacent layers in our model are variable in both space
and time in order that each layer can keep track of the changes that occur
in the particular set of phenomena that layer is designed to describe
(summarized in Figure 1-1). A consequence of this structure is that the
volumes of the grid cells vary in both space and time. By contrast, in
-------
«- >
I/I
c
.o
*j
o
3
LL
I
. re
«•- O o>
0 c s = •£
5 .2 g S -
„ «- S E O g
0)
**
I
0)
i_
3 <0
4J C
O O>
3 E
J_ O
•<-> C
10 O)
jC
i. 0..
d)
>» QJ
«3 E
U
I
(O
C
>1
•o
= T3 -M
c re
O) fO i— 4J
(U
M- TJ O
O O 4->
c E-o
O OJ ,
•M (O
(O JZ
J= U
-------
,0
*J
u
LL
OJ
CO
01
T5 S
il
u •£
•o
a E *
3 U
c 5
re •" c =
»- •> «_ IB
02 O *•
re
t»
re
O
r
o
Q..
= c • '
-I I*
k. O C C t/i
o -a 5 jj re c
I"S ;2™ "I-a E
y~ o to to "jj j? «• M
™ ^J i» at (U u M
« c -
= 5 a-**
2. ! «A •
C ** ;S
3 O «
- - _ . " •§ °
1.2 I S .£ -a ^
l-o e o „ e g
g 2 JC Jj g g O)
0)
S5
o <3
1
o
u
n
o S
•o
01
X
;^
•a IB
6 -1
•^^N— •
/
«v
Radiation
Inversion,
Nocturnal Je
1
.n
„
•S §1
^^
o
c
(U
CL
a.
0)
u
X
a>
10
.a
t—i
t
-------
conventional models the grid network and cell volumes remain fixed and
surfaces such as the mixed* layer top move through the grid system. In
the following several paragraphs, we elaborate on some of the phenomena
cited in Figure 1-1 that our model will take into account.
During the day the highest layer shown in Figure 1-la represents the
synoptic scale subsidence inversion, which may or may not contain cumulus
clouds. Stratospheric 0^ is transported downward through this layer and
anthropogenic CL and its precursors can be carried into it by cumulus
clouds or penetrative convection, the depth of the penetration being deter-
mined mainly by the intensity of the temperature stratification in the in-
version and by the scale and intensity of convection in the mixed layer.
The base of this layer is normally 1 or 2 km above ground level. Below it
pollutants are kept well mixed vertically by turbulent convection. If the
winds are strong or the surface heat flux is weak, wind speed and direction
may vary appreciably within the first several hundred m above ground. There
is usually also a marked difference in the wind speed and direction between
the inversion layer and the mixed layer below. Over large lakes and along
sea coasts there is frequently a second inversion layer below that generated
by synoptic scale subsidence. This lower inversion is produced by sea or
lake breeze regimes and it restricts the vertical mixing of pollutants
emitted over the water and within several km inland from the water's edge.
Air drawn into young cumulus clouds originates primarily in the lower
portion of the mixed layer. Fresh emissions of NO and hydrocarbons can
^
transported by the vertical currents that feed these clouds from ground-
level to altitudes well above the top of the mixed layer in one steady, up-
ward motion. In the process little or no mixing with aged pollutants in
-------
the mixed layer occurs. At night, cumulus clouds usually evaporate and, when
they do, they leave behind products of liquid phase reactions that can be
transported hundreds of km before sunrise.
When the daytime mixed layer is moist enough for cumulus humilis clouds
to form at its top, pollutant emissions can undergr . -asi-cyclic series of
events that may have a great effect on the chemical rate of pollutants and
particularly the formation of aerosol (Figure 1-2). As shown in the figure,
the thermals, or vertical jets of warm air, of which cumulus humilis are a
part, have their roots near ground-level. Fresh emissions of pollutants
that enter a thermal can rise in 1. steady motion all the way to the top of
the mixed layer. Individual clouds mark the locations where air parcels
that compose thermals have become saturated and liquid water has formed.
In the saturated parcels, air pollution chemistry will be altered by the
sudden presence of liquid water and also by the accompanying-attenuation
of sunlight. The presence of cumulus humilis signifies that the upper layers
of the atmosphere are too stably stratified for thermals to penetrate higher,
and so in this instance the cloud tops bend over and the air of which they
are composed dives back down into the mixed layer (see Figure 1-2). As a
saturated parcel descends, it is warmed by compressional heating and even-
tually the cloud droplets within it evaporate leaving behind the residue
of aqueous phase reactions, probably in aerosol form. After some time the
parcel will arrive again at ground-level where new emissions may be injected
into it, and/or deposition on surfaces may occur; subsequently, the entire
process is repeated. The time required for 1 complete cycle would typically
be 30-50 min with maybe one-tenth of the time spent in the cloud stage.
The chemistry associated with the cycle we have just described could possibly
be studied in the laboratory by measuring the composition of initially moist,
10
-------
Geostrophic
wind
Cumulus
Humilis
/ / f// / /" 7
S S S / ' S / /
Figure L2. Idealized view of the fate of polluted air parcels in a
cumulus humilis topped mixed layer. A, parcel receives
emissions; B, it rises and cools; C, saturation occurs
producing cloud droplets; D, the parcel descends and
drys; and E, it receives new emissions and/or deposits
pollutants on surfaces. From E the cycle begins anew.
11
-------
polluted air as it is subjected to a series of expansions and compressions
in a cloud-chamber apparatus. The framework of the model that we develop
later will allow treatment of this process.
Dramatic changes occur in the mixed layer at night. With onset of
surface cooling following sunset, a stable layer of air forms near the
ground that quenches the vertical momentum fluxes that give rise to fric-
tional drag on the horizontal flow. With retardation forces eliminated,
the wind just above the stable layer accelerates giving rise to the pheno-
menon known as the nocturnal jet. Wind speeds in the core of the jet, which
usually lies between 300 and 500 m above ground, may be 10-15 m/s while
at the same time the air is nearly calm at the surface. Emissions from tall
stacks and from sources within the urban heat island enter the jet region
at night. There they react with aged pollutants from the previous day and
are transported considerable distances by the strong flow. The remnant of
the previous day's mixed layer above the jet is isolated from the influence
of fresh emissions and it moves at a slower speed than air below.
Sporadic episodes of turbulence in the shear layer beneath the nocturnal
jet is a mechanism by which 0^ and constituents of urban plumes are brought
to ground-level at night. There deposition on surfaces and reactions with
emissions of small, low-level sources occur. This sporadic mixing process
is perhaps the only mechanism by which the reservoir of aged pollutants aloft
can be depleted at night.
One point that we wish to emphasize here is that 1-layer regional scale
air pollution models are incapable of simulating the effects on pollutants
like 03 of the vertical segregation of aged and fresh emissions that occur
at night. Being cut-off from contact with the ground and fresh NO emissions,
A
12
-------
03 above the nighttime radiation inversion ~ free to travel great distances
before it is mixed vertically by convection the following day. The effect
of this nighttime segregation of pollutants is to extend greatly "the effective
residence times of species like 03 in the lower troposphere. Consequently,
a multi -layered model seems to be essential to simulate accurately the long
range transport of photochemical air pollutants.
The horizontal resolution of the model that we now plan to make opera-
tional is about 18 km. We want this resolution to be as high as possible
to mitigate the effects of subgrid scale concentration fluctuations. Since
apparently few modelers are aware of this phenomenon, it is perhaps worth-
while to digress here with an explanation of it.
Suppose we were simulating 03 concentrations within an urban region
covered by the 16 grid cells shown in Figure 1-3. Suppose further that at
some instant transport and diffusion ceased leaving NO and 03 unmixed and
partitioned in the manner shown. If NO and 03 are the only species present,
then clearly there will be no change in the concentrations in any of the
cells. In this case our model, using the rate equations
would correctly predict
d NO _ d°3
_
-at -
where the overbar denotes an average over any 1 of the 16 cells shown in
Figure 1-3.
Suppose now that this same concentration distribution were contained
within a single cell (denoted by the heavy outer square in Figure 1-3) of
13
-------
NO
NO
NO
0.
NO
NO
NO
NO
NO
Figure 1-3.
Hypothetical case of an urban area in which ozone and NO
are segregated into 16 cells. If this distribution is
simulated in a regional model whose cell size is the large
outer square, subgrid chemistry phenomena arise.
14
-------
a regional model. Denoting mean values wit;.
-------
the meteorological and chemical data required to formulate the model. A
second goal of NEROS was to gather the data required to perform comprehensive
test runs and validation exercises of the model.
A second problem is limitations of computer storage capacity. To simu-
late air quality over the Northeastern United States w the horizontal
4
resolution we desire, our model will have roughly 10 g.id points and it
will treat 23 (eventually more) chemical species. Thus, the concentration
variables alone will require 250K words of storage and this is just under
the working limit of 260K words of memory on EPA's Univac computer. To
accommodate a model of the anticipated size we will have to develop special
techniques that permit handling the modeling domain in piecewise fashion.
This problem is addressed in Section 9 of this report.
Finally, due to the large number of processes that we will attempt to
treat, our model will be rather complicated. In order to alleviate the prob-
lems that this might cause in operating the model and in making future re-
finements, we will structure it so that its central core consists solely of a
set of algorithms for solving the coupled set of generalized finite different
equations that describe processes in each of its layers. The modeling
functions of describing the mixed layer dynamics, topographic effects on
wind, chemistry, cloud fluxes, etc., will be handled by a set of special
processors that are external to the central model and which feed the model
key variables through a computer file. Within this framework the techniques
used to describe the various physical processes can be altered without over-
hauling the model itself. An additional advantage is that execution times
are greatly reduced when several runs of a given scenario are to be performed
in which only 1 or a few parameter values are altered. Details of the model
16
-------
system are given in Part II of this report.
In the sections that follow we develop in detail the theoretical bases
of the regional model. Operational aspects, refinements, validation exer-
cises and other details will be presented in subsequent volumes of this
report.
17
-------
SECTION 2
DERIVATION OF PROGNOSTIC EQUATIONS FOR LAYERS 1, 2 AND 3
Let c denote the concentration of any pollutant species we wish to
model. Its mass continuity equation is
If + !H ' (c!} + fe (cw) = s + R • w (2'1}
where
VM = the horizontal del operator
v * the horizontal component of the wind,
w = the vertical component of the wind,
S - the source strength function of the species,
R * the net rate of production of c by chemical reaction,
W = the rate of removal by rainout and washout.
Deposition is taken into account by the boundary conditions on (2-1).
In cases where c represents a photochemical pollutant, (2-1) is coupled
through the reaction term R to the equations governing other pollutant con-
centrations.
The surfaces separating the model layers will be expressed in the
implicit form
H^x.y.t) = z.(x,y,t) - z, i = 0,1,2,3 (2-2)
18
-------
The terrain surface will be denoted by ZT nd we will allow it to extend
up through HQ and H, under certain conditions, which we will define later.
Each of the interface surfaces H is prescribed later, in Section 3.
Let be an arbitrary function of (r,t). We define the mean value
4> in Layer j (= 0,1,2,3) to be
<4>(x,y,t)>.= rr
rr
X+AX
X-AX
y+Ay
y-Ay t
'ZjU'.y'.t)
*(x',y1fzi,t) dz'dy'dx1
(2-3)
where
V,(x,y,t) =
J
X+AX
X-AX
y+Ay
y-Ay J
z,(x',y',t)
J
dz'dy'dx1
(2-4)
In order to derive from Equation (2-1) the equation governing the mean
pollutant concentration . in each layer, it is convenient at this
J
point to introduce some of the properties of time -and space derivatives
of volume averaged quantities defined by (2-3). Consider first the
time derivative
From (2-3) we obtain
|fj -
ffj-1
l£dz'] dx'dy' -^
. av
J 3tJ
(2-5)
where
dx dy denotes integration over an area 4AxAy centered at
19
-------
From Equation (2-4), we obtain
(!fd - Mi - 1) dx'dy'
'3t
at"
(2-6)
Let us define the following surface average:
rx+Ax
4»[xl,ylfzj,t]dyldx1
X-AX
(2-7)
where
A = 4AxAy
(2-8)
This is simply the average of 4>» evaluated on the surface H., within the
J
rectangular area A centered at (x,y). We emphasize that (2-7) is not an
average over the projection of A on the surface H.. Using (2-7) and the
abbreviated notation
2. = 9Z./3t
we can write (2-6) in the form
IV,-
3tJ
Then, Equation (2-5) gives
It
• v.
J
If
Consider next the space derivatives of <>:
17
(2-9)
(2-10)
j 3 <* vj
20
-------
where
•y+Ay
y-Ay
(2-13)
We can transform Equation (2-12) into more familiar terms by comparing it with
f^z'dy'dx'
z.
j-l
This equation can be simplified starting with
/•y+Ay
3zi 3zi 1
[(x',y'iZ,,t) T-4 - *(x',y'tz. ,,t) -r4
J OA J-l OA
y-Ay
(2-14)
|frdz']dy'
(2-15)
Combining Equations (2-14) and (2-15), we obtain
- If jj • ?.
- Vx - 4x)] + £'* If^-1
j
I2-:) (2-16)
and substituting from Equation (2-12), we find that
-j-l
_ 3
—•- .
3X J
By analogy
—j-l
—i 1
3yJ""
fy
(2-18)
21
-------
Finally, we see from Equations (2-3) and (2-7) that
V.
J
j -f ( 4J -I*3'"1) (2-19)
J
We now have the relationships needed to derive the equations governing
.. Performing the < > averaging on Equation (2-1) and making use of
J
Equations (2-11) and (2- 17) -(2- 19) we obtain
' <*> + J 5T & W <^j * T7 J $
J J
J j j 1 4 1
- 7, <- IfJ + - & + 7, '- If^1 + « ^-»
J J J
. - j - j - 1
+ §• (cw - cw ) = . + . - . " (2-20)
j J J J
where v = (u,v). Recall that z = 3z/8t.
Collecting terms representing averages over each of the surfaces H. and
J
H. ,, we can express (2-20) in the more concise form
J "* X
t j j t vj + !H ' ^j + j
J |y & V0 = 7
where
F. = c SSd (2-22)
J
22
-------
v = (u,v) and H. is given by (2-2). Note that F. is the total flux of the
— J J
species represented by c across the surface H. within the area formed by
J
the projection of A, centered at (x,y), on H..
In modeling atmospheric processes over 1000 km scale regions, earth
curvative effects must be taken into account. This can -~a done either
by transforming the governing equations into one of the _c:ilinear pro-
jections of the earth, such as polar stereographic, Albers equal area,
etc., or by casting the equations in curvilinear coordinates, such as
latitude - longitude. In either case the equations acquire a form different
from the well-known Cartesian form we have developed above. In Appendix
A we give the details of the transformation of Equation (2-21) into a
curvilinear frame in which latitude, longitude and elevation are coordinates
and the basis vectors point north, east and vertically upward at every point
on the earth. We have chosen this frame because it is a natural one from
which transformations to any rectilinear system are easily performed.
Also, it is the frame in which worldwide meteorological data are reported.
In the chosen frame Equation (2-21) takes the form
. 3. . 3- 0
3 j j + "x -3TJ + u+ ~5rJ (2-
«x
where y and y, are the metric factors
A ' 4> a
a is the earth's radius at msl; X is longitude; 4 is latitude;
v is the north-south wind component; u is the east-west component;
23
-------
and
The product mean terms that enter in (2-24) can be simplified through
use of the continuity equation for the atmosphere. In studies of flows
that are shallow comparable to the scale height of t • -nosphere, the
continuity equation can be simplified to the form
v • v * 0 (2-26)
(see Haltiner 1971, page 54).
Applying the < > averaging operation to this equation we obtain through
a process like that described in Appendix A,
8£nV_.
" * " * "
X , X j , J
7. t"x u fl"-1 + ", » fj-1 + « fj-H " <2'27'
J
- f, I"x u f J + % v f1 + « &1 ' °
In this equation, as in (2-24), the terms under the overbars are evaluated
on the surface H. whose derivative is the product term.
J
Let us define the flux terms
= -
= - (2-28)
applicable to any of the 3 model layers. Note that under these de-
finitions, and are not necessarily zero. Making use of (2-27)
and (2-28) we can express (2-24) in the equivalent but more manageable form
24
-------
It
where
3H~
F. . = c.v • VH. - , v • VH. + c. ||j (2-30)
J » 1^ J~~J' I\--J JoL
and
n — 0 i v O
-u A oA u) o u) vZ
and c. and v denote concentration and velocity, respectively, evaluated on
surface H.. By definition [Equation (2-30)], F. , is the rate of transport
J J >K
of material c across the surface H., to or from Layer k, per unit horizontal
J
area (not unit surface area of H.) with H. held fixed in time. Equation (2-29)
J J
is the general form of the model equation with which we shall work in the
remainder of this report. In the sections that follow, we develop expressions
for each parameter that enters in this equation.
25
-------
SECTION 3
THE INTERFACE SURFACES H,
J
We summarized earlier in Figure 1-1 the various roles that each of the
models' 4 layers are intended to serve. In this section we will specify
functional forms for the layer interface surfaces H. (A,,t) that we believe
J
are compatible with those roles.
The Model's Top Surface, H,
There are two conflicting constraints on the surface H^. The first is
that ZgCx.y.t) must be high enough that no appreciable turbulent or cloud
fluxes of material cross it. The second requirement is that the depth of
Layer 3, i.e., Zg - i?* mus^ ^e sma^ enough that subgrid scale concentration
variations within the layer are negligible. To satisfy both these require-
ments, we will prescribe zJx.y.t) based on observations of convective cloud
top elevations and on the estimated level z~ of the "mixed" layer top.
Further details will be presented in Part 2 of this report.
The Subsidence Inversion Base Surface, Hg
This is one of the most important parameters in the model; because acting
as a barrier (a "leaky" barrier when cumulus clouds are present) to the upward
movement of pollutants, the height of Hg above ground determines the concen-
tration that fresh emissions produce once they become well mixed vertically.
26
-------
As a consequence, H~ has an1 effect both on the rates of slow second order
chemical reactions that occur among the constituents of polluted air as it
moves over long distances and on the rate of surface deposition. In view
of its importance, the surface H2 (x,,t) is treated rather rigorously in
the model. We will postpone further discussions of it until Section 4,
where we describe in detail the procedure we use to calculate H£.
The Nocturnal Jet Top, H,
We believe that an essential requirement for the model's ability to
simulate accurately the long range transport of CL is the ability to handle
properly the deposition of nighttime emissions of hydrocarbons and NO .
A
Preliminary results of aircraft studies performed as part of project NEROS
in August, 1979 support our original speculation that nighttime surface
emissions are confined primarily to the lower one-quarter of the old daytime
mixed layer, and that 0^ and aged precursors reside in the upper three-
quarters where they undergo virtually no mixing with the fresh emissions below.
It is this segregation of the old and fresh emissions at night that allows
0^ to travel great distances before being scavenged by NO or hydrocarbon.
Thus, the surface H^ should be set just above the level that nighttime emis-
sions reach. This level is probably controlled by the temperature profile
of the nocturnal boundary layer, by the wind speed and by the buoyancy flux
of the urban plume. One of the objectives of the NEROS field program is to
determine just how these factors control the nighttime urban plume rise and
spread.
Another role intended for surface HI is to represent the barrier to
upward diffusion imposed by the weak subsidence inversion associated with
27
-------
lake and sea breeze regimes. The base of this inversion is lowest over the
water, but it slopes upward inland until at some distance from the shore
it vanishes altogether. This feature of the boundary layer probably plays
an important part in the dispersion of emissions around the Great Lakes and
along the Atlantic Coast.
A third role of Layer 1 is to describe the daytime shear layer in
which mean transport Has a different direction and speed than in the upper
part of the mixed layer.
- To serve in all these capacities, the surface H, should be calculated
from the set of equations that describes the nocturnal inversion layer and
urban plume rise at night, the lake/sea breeze inversion by day and the shear
layer depth. At the present stage of the model development we attempt to
describe H, rather rigorously under stable conditions but only approximately
in neutral and unstable situations. In the former case we prescribe the
volume flux T\, across H,, i.e.,
dH,
-j^ = rij t stable conditions. (3-la)
Using computed values of the fluid velocity (u,, v,, w,) in the stable boun-
dary layer and the estimated value of nj, we solve (3-la) for H, during
stable conditions. Details are given in Part 2 of this report. During
neutral and unstable conditions the methods employed in the stable case do
not apply and we rely then on an ad hoc formulation of z, itself. The ob-
jective is to account for the shear layer depth and terrain effects mentioned
above but not the lake/sea breeze inversions. The same ad hoc method is
used to describe H« under all conditions. We describe the method below.
28
-------
First, we assume that at any point the vertical velocity w is composed
of 2 components: one, WT, induced by the motion of horizontal flow over
terrain features and the other, WD, arising from turbulence and horizontal
divergence of the large scale flow. That is,
w = WT + WD (3-lb)
If the vertical velocity component WD were zero and the horizontal wind
(u, v) were steady, H, would be a material surface if
?H1 * !HZ1 = WT1 (3- ic)
where VHI
!H E u* IT + -4
and where the subscript "1" denotes evaluation at level z,. We can express
W in the form . .
where the term in brackets is a measure of WTQ, i.e., vertical speed at the
ground surface z,., and A, is an unknown function (which we approximate below)
If we express z, in the form
where H. is the depth of the shear layer over smooth terrain, then we find
on substituting this expression and (3-le) into (3-lc) that in order for z,
to satisfy (3-lc), H, must satisfy
v^ + (1 - Aj^z-p = 0 (3-lg)
Thus, once we have an estimate of A^ we can use (3-lf, g) to define a sur-
face HJ that will satisfy (3-lc).
29
-------
To approximate A we first define 2 smoothed terrain functions
x + 2.5
x - 2.5
x + 25
x - 25
•y + 2.5
zT(x',y')dy'dxl
y - 2.5
•y + 25
zT(xl,y')dyldx1
y - 25
(3-2a)
(3-2b)
where (x, y) are Cartesian coordinates centered at (x, ), ZT is local topo-
graphy elevation (MSL), and distances are measured in kilometers. We next
define a rms topography
rx + 25 ,y + 25
1
" 2500
x - 25 J
^ ^
[zT(x',y') - zT(x',y')]2 dy'dx1 (3-3)
y - 25
Using ZT we define a terrain Froude number F by
(3-4)
where N2 = g/e gj is the Brunt-Vaisal a frequency and U is the horizontal
wind speed at an elevation IT above the level z,.. According to Godowitch
ejt aj_. (1979), nighttime inversions over rural areas around St. Louis were
observed to have intensities of between 10 and 20°K per kilometer. Thus,
we will adopt N = 0.33 s as a representative value for use in our modeling
studies of the Northeastern United States.
Hunt e_t al_. (1978) found in laboratory studies that stably stratified
fluid tends to move horizontally rather than over an isolated hill of
30
-------
moderate slope when the Froude number is less than 0.3. Based on this
finding we propose as an interim approximation
Ak = exp [- (2JM2 - (tfk/3sT)2] , k = 1 or 0; (3-5)
and F is given by Equation (3-4) with
(0.03 s , nighttime hours
N--< (3-6)
j 0.0 s , daytime
The last term in brackets in Equation (3-5) is intended to account for
the decrease in the amplitude of vertical velocity with height that is asso-
ciated with flow over obstacles. We expect that Hj will have a value of
about 300 m. Thus, terrain features smaller than about 50 m will have
virtually no influence on the slope of z,.
In summary, during stable conditions we prescribe
dH,
-g-r- = n, » stable conditions; (3-7a)
where n^ is a function given in Part 2 of this report; and under neutral
and unstable daytime hours we define
Zj = ZT + Hl (3-7b)
where H^ satisfies
7HHr + (1 - A^v^ = o (3-7c)
in the interior of the model domain and has the value of the shear layer
depth at the boundaries. With A. given by (3-5) we assume that z, satisfies
!HI • !HZI = wn (3-7d)
In this case we have
31
-------
dH, 32.
~dT s nl = ~TF " WD1 ' neutral and unstable (3-7e)
•
Thus, in stable conditions ru is prescribed and z. is computed, and con-
versely in neutral and unstable conditions.
All our assumptions regarding z, are speculative and are intended to
serve only as a first order estimate of the effects of topography on the
vertical transport of pollutants. More refined representations of the sur-
face H, that are consistent with the roles that Layer 1 serves in the model
can be incorporated as they become available without having to overhaul the
model (see Part 2 of this report). Further discussions of the surface z,
and its roles will be given in Section 4.
The Surface Layer Top, HQ
The 2 main roles of Layer 0 are to treat the subgrid scale concentration
fluctuation effects on chemical reactions and to parameterize surface removal
processes. Both phenomena are confined to a relatively shallow layer next
to the ground. Therefore, it is possible through a judicious design to for-
mulate the Layer 0 equations in a way that does not require additional com-
puter storage. In other words, the Layer 0 equations can be written in a
diagnostic form in which concentrations in Layer 0 are expressed in terms
of concurrent concentrations in Layer 1, surface source emission rates and
the like. The details will be given later. The only restriction on such a
formulation is that the total mass flux of material through the top surface
of any Layer 0 cell be much larger than that through its sidewalls. If L
denotes the horizontal cell size and HQ the depth of Layer 0, then we require
32
-------
I 2
Vertical Flux A °w _ - a Q\
Horizontal Flux ~ ( '
where a is the effective turbulence velocity on the top surface of Layer 0
w
and u is the mean flow speed within the cell. This condition is satisifed if
L/HQ » u/aw (3-9)
The velocity ratio on the right side of this expression is typically
of the order of 10 or less throughout the depth of convective mixed layers.
In stably stratified flows over flat terrain, it varies from a value of
about 10 a few meters above the ground to a value of several hundred near
the top of the nocturnal boundary layer (roughly 100 m).
Thus, using a grid cell L ^ 20 km we can let H be as large as 100 m
during the day (this will allow us to encompass the bulk of the subgrid
scale concentrations variations) but at night we will have to limit H to
about 30 m.
Since we will neglect horizontal motions in Layer 0, we must define
HQ so that there is no net transport of air across it. That is
VQ • VHQ = 0 (3-10)
To meet this condition, we will assume that the average divergent component
of the vertical velocity, i.e., WD [see (3-la)] is zero at z = ZQ, rather
than at z = ZT, and we will define ZQ so that it follows approximately the
terrain induced flow, that is,
"TO ' Vo % + Vo TT (3-U)
where wTn is the component of w at z induced by flow over terrain. We
'TO
33
-------
can satisfy this requirement and at the same time fulfill the constraints
on the depth of Layer 0 defined above by making ZQ similar to z^ in form,
namely
= zT(X,*) + HQ(t) (3-12)
where AQ is given by (3-5), HQ satisfies
in the interior of the model domain and on the boundaries has the value
«(t)
75 , t$R + 3 < t < tss
75 - (t-tss)-15, tss <
0 30+(t-tSR)-15, tSR< t
-------
SECTION 4
THE INTERFACE FLUXES: F. ,
J ,K.
In this section we derive explicit forms for the interface fluxes F. .
J , K
that enter in the governing equations (2-29) of each of the model's layers.
Flux Across the Model's Top Surface: £ 3
We concluded in the previous section that the elevation z3 of the top
surface H3 of the model must be a prescribed function. In this case (2-30)
gives (see also 2-2, 2-22 and 2-23)
f3,3 B C!f ' !H3 - 3 ! • !H3 * C
dH
where c- denotes concentration on the given surface H3- We will assume that
dHT
3 , 1f-gf <0
(4-2)
, otherwise
where c is the concentration of species c in the free atmosphere above the
oo
model domain. Combining (4-1) and (4-2), we obtain
E
3z3 dH3
3lt • 1
3,3 - -) — (4'3)
T 9Z
(c - _) —TT- j. 1 , otherwise
3' dt * 3 "af
35
-------
This approximation of F, - assumes imp'i citly that there is no turbulent flux
of material across FL. No error should be introduced by this assumption if,
as we intend, Hg is above the top of the mixed layer and above the top of any
convective clouds that might be present.
Flux Across the Subsidence Inversion Base: F2 2 and ^ 3
We have defined H2 as the base of the synoptic scale subsidence inver-
sion. During daylight hours this surface marks the top of the region known
as the mixed layer in which, due to the action of convection, moisture and
pollutant concentrations become approximately uniformingly distributed
vertically. Observations have shown that potential temperature, water vapor
mixing ratio, and often pollutant concentration change abruptly at the
mixed layer top. Based on this knowledge we shall assume, following Lilly
(1968), that all 3 of these variables have discontinuities on H2- This
approximation will be used below to derive the expressions both for the
pollutant fluxes across H2 and for the spatial and temporal variations in
H2 itself.
H:
i
Figure 4-1. Schematic view of cumulus clouds transporting
material aloft from the mixed layer and acting
to lower the surface H2.
36
-------
We consider the general situation depicted schematically in Figure 4-1
in which cumulus clouds transport heat, moisture, and pollutants through
the mixed layer top into the cloud layer above. The upward motion in the
clouds produces a compensatory downward velocity between clouds that acts
to lower the mixed layer top.
On H2 we assume that the concentration c2 at any point is given by
(4-4)
c , beneath cumulus clouds
C2 *•
2 + c'2 , between clouds
where CG is the concentration of species c in the updrafts feeding cumu-
lus clouds in the given grid cell. This variable is a function of space
and time and will be evaluated later in this section.
Let H2+ be a surface parallel to H2 that is just high enough that the
dry thermals that overshoot level Z2 between cumulus clouds cannot reach it.
On H2+ we assume that the local concentration c2+ is given by
cc , beneath cumulus clouds
(4-5)
3, between clouds
t
The latter expression is a result of our choosing H2+ high enough that
there are no fluctuations in concentration on this surface due to the tur-
bulent convection below.
Let a denote the fraction (0 <_ a < 1) of any given cell area A
c ~~ c
covered by cumulus clouds. Spatial and temporal variability of oc will be
allowed but not explicitly stated.
For mathematical convenience we will express F? « and E, - in the form
t-,L (. , J
(4-6)
37
-------
where
dH
and
L. - L.
dt = at I2• ,
and v2 is the velocity vector measured on H2. Usi
obtain
dH;
dt
(4-74)
(4-7b)
) and (4-7a) we
+ da ] (4-8)
1 - a.
where the first integral is over the portion of H2+ beneath cumulus clouds
and the last integration is over the remaining area within a grid cell.
Assuming that there is negligible variation in c from cloud to cloud
within any grid cell, and that spatial variations of 3 within any cell
are small, we can reduce (4-8) to
(4-9)
where ( ) represents a mean value over the portion of H2+ beneath clouds
and ( ) represents a mean value over the environment between clouds
within any grid cell.
Through a process similar to that leading to Equation (4-9), we obtain
for the flux F2 on H2:
(2 + c2)
1 - a^
dt
da]
(4-10)
•Vca? + (i --e) t
-------
Continuity of material flux across H2 and the assumption that the
distance separating H2 and H2+ is small require that
F2 = F24. (4-12)
Thus, equating (4-9) and (4-11) and noting that since H2 and H2+ are
parallel,
dH? _ dH2
dt dt
we obtain
(4-13)
(4-14)
By definitions (2-2) and (4-7b) the downward mass flux of air across
the moving surface H2 (i.e., measured by an observer at a fixed point on
surface H2) is
(4-15)
where p is the a'ir density on H2, A is the area of the projection of a
unit horizontal square on H2 (see Figure 4-2), and F is the total downward
mass flux through area A. It follows
unit horizontal
Figure 4-2. Projection of a unit horizontal square on the
surface H~.
39
-------
from (4-15) that p dH2/dt is the dowm j mass flux of air per unit
horizontal area, rather than unit area of H2. Therefore, the last term on
the right hand side of (4-14) is the mean downward mass flux relative to
H2 that includes the effects of mean vertical fluid motion, mean motion
of surface H2 and turbulent convection.
For notational convenience we will write
Substituting this into (4-14) we get
^«^^^ o
BI2 = VAC (4'17a)
where
AC = 3 - 2 (4-17b)
An important point to note here is that 3 in this last expression
is used as a measure of the concentration on surface H2+ [see (4-5)] which
penetrating convective elements do not reach. This fact will be used
later when we develop an approximation for f and the rate equation
c
governing H2.
The analyses that led up to (4-17a) are applicable also to other scalar
quantities, such as potential temperature 8. This fact and (4-17a) indi-
cate that
fc/Ac = fe/A8 (4-17c)
where
V^IF6 M-17d)
and
A6 = 3 - 2 (4-17e)
40
-------
As in the case of concentration, 3 i- .reated as a known quantity
and is used here as a measure of 9 at the level just above the altitude
of the highest reaching convective parcels from the mixed layer, or equi-
valently, the level where the turbulent heat flux becomes zero. Since
it is fQ/A6 that we will parameterize later, we shall take advantage of
o
(4-17c) and substitute it for f_/Ac from here on.
The surface H2 is smooth in that it is not deformed abruptly within
individual cumulus clouds. As a result (4-17a) is equivalent to
7UZ2 - W = + f./A6 (4-18)
„ rl 6 o
where w is the mean vertical velocity over the environment between cumulus
clouds.
By definition, the mean vertical fluid velocity w2 at any instant
on fL is
W2 = acwc + (1 - ac) we (4-19)
where w is the mean upward velocity in cumulus clouds in the given cell.
This expression and (4-18) yield
~= w2 - a_w
3z? , — „ _ _ ^ c c . ,
The flux F2 2 [see (4-6)] may now be expressed using (4-9), (4-12), (4-13),
(4-18) and assumptions invoked above concerning the variations of H2 within
the horizontal scale of a single grid:
2,2 c c 1 - a 3 c
r ac ,- - 3Z~2, (4-21a)
ET r- (wo - wj + fe/Ae - TT-]
41
-------
and by similar means we obtain from (4-6)
F2,3 = ac[-TTir + fe/*9^c - 3>+ 3 "if
Except for c all of the independent variables in (4-21) are inter-
related [see, for example, (4-20)]. Thus, to obtain meaningful estimates of
F2|2 . (4-22c)
where e2 is potential temperature on surface FL [cf. (4-4)]. A plausible
assumption regarding e2 is
62 =
2 , if WR >_ 0
(4-23)
2 + AS' , WR < 0
•^»
Since fluid that descends through ^ originates somewhere between levels
and ^2+* tne fluctuation A61 in (4-23) must lie in the range
42
-------
0 < A91
-------
to the mixed layer depth, about 2 km. Here we shall extend this assumption
to include averages over the environment between cumulus clouds represented
by the overbar and superscript "e" in (4-22a). The validity of this approxi^
mation requires that the fractional coverage OG of cumulus clouds be less
than about 0.5.
With this assumption, (4-22a) can be written
fe - -
wRp(wR)dwR + e2
wRp(wR)dwR]
(4-24)
where P(WR) is the probability density of vertical velocity WR relative to
H2, as defined by (4-22b). Combining (4-22c), (4-23) and (4-24) we get
fe = - Ae'
wRp(wR)dwR
(4-24a)
where A9 is given by (4-17e).
Figure 4-3 shows the probability density of the vertical velocity w2
at the inversion base H2 as given by the numerical model of Deardorff (1974)
It indicates that to good approximation, w2 has a Gaussian density, that is,
Pz(w2)
(4-25)
where w2 is the mean vertical velocity at H9, due to subsidence, and a is
fc W
the variance, which is proportional to the turbulent energy at H?. Since
the 2 last terms on the right side of (4-22b) are deterministic variables,
it can be shown that P(WR) is also Gaussian with mean
WR = we - a!2/3t - v2H-vHz2 (4-25a)
and variance a . the same as that of w2. The w\ term enters this
w e
44
-------
ALTITUDE = 1100 M ^ MIXED LAYER TOP
0.90 H
0.80 -
0.70 -
0.60 -
050 -
0.40 -
0.30 -
0.20 -
0.10 -J
0.00
1.00 -0.00 1.00 2.00 3.00 4.00 5.00
•4.00 -3.00 -2.00
W (M/SEC)
Figure 4-3.
Probability density of vertical velocity w« near the mixed
layer top (H~) derived from the numerical model of
Deardorff (1974) .
45
-------
expression because we are averaging over the cumulus free environment only,
as indicated by (4-22a).
The Gaussian form permits the integral in (4-24a) to be evaluated analy
tically. We get
- _ _ 2
Wp Wp a Wp
f, - - «• if ( I - erf (-B-) ) - -*- exp (- ,£,)] (4-26)
/ 2 w /2ir w
where erf denotes the error function [see (4-23)]. Combining (4-25a) and
(4-18) we obtain
- WR = fQ/A9
and upon introducing (4-26) we get after some rearrangement of terms
J [1 - ?/2 (1 + erf J)] = -I— exp (- J2) (4-27a)
~
where
WR
J * - — - (4-27b)
The solution J(T) of the transcendental equation (4-27a) in the interval
0 <_ T < 0.9 1s approximated closely by
0 < T < 0.9 (4-27c)
3737 - 0.86T - 1.46T
This relationship is plotted with the exact solution of (4-27a) in Figure 4-4.
T
Thus, from (4-27b), (4-25a) and (4-19) we obtain the mixed layer growth rate
equation _ _
" (4-28a)
at"
and
fa/A9 = f /AC = ^aw J(T) (4-28b)
y c *•
46
-------
J(T)
1.60-
1.40-
1.20-
1.00-
0.80-
0.60-
0.40-
0.20-
0.00
exact solution
3.37 - 0.86T - 1.46T2
0.00 0.20 0.40
i
0.60
0.80
i
1.00
i
1.20
i
1.40
I
1.60
Figure 4-4.
Comparison of the true solution of equation (4-27a) with
the approximation (4-27c).
47
-------
The vertical velocity fluctuations on H« are generated by convection
within the mixed layer below. Hence, it seems reasonable to assume, as is
often done, that a is proportional to the convective velocity scale w*:
a = aw* (4-29a)
where
w* = (g/e (w'e'Lh) f'
(4-29b)
Here 0 is the mean potential temperature in the mixed layer:
/X+AX
2
..e(A,*,t) 'TTT11
JX-AX
e(x' ,<}>' ,t) d'dx'dz
(4-30a)
A = 4AXA4»a cos* is the grid cell area;
h is the mean mixed layer depth in the cell,
X+AX
X-AX
22 - Z.
Z2(X' ,' ,
(4-30b)
(w'e1)Q is the kinematic heat flux at the ground; g is gravity and a is a
constant. Laboratory investigations of Willis and Deardorff (1974) and
numerical model results of Deardorff (1974) indicate that the value of a is
about 0.4.
The heat equation, which governs 0, is of the same form as (2-29) so
we have at once
I1
- 9
dH, _
' j. Q \»
T-jr. T O— V-
dt T .
(4-30c)
48
-------
where
^H'9'^" ^H^M ' ^HV (4'30d)
Here M denotes the mean horizontal wind in the mixed layer, e2 and v2
are the potential temperature and wind vector on H2, HT represents the
ground surface where e = 8j, and is defined by (2-28).
The last term on the righthand side of (4-30c) is identically zero
because the ground is a material surface. The third term on the right side
of this expression is just the kinematic heat flux at the ground, namely
"• 9T-dt "^To = fes {4'31)
Comparing the first 2 terms on the righthand side of (4-30c) with (4-6),
which is the definition of F2, we see that the sum of these terms is
equivalent to (4-21), but with 2 replaced by 0; 3 - 2 replaced by
A6 [see (4-17b) and (4-17e)]; and with cj. replaced by
ec - se (4-32)
where £ >_ 1 (assuming that the ground surface is warmer than the mean
temperature e of the mixed layer). Substituting these results into (4-30c)
we obtain
90
aT + H>M * ^U ® + <^U ® M * ^U^2 "*" n^u * 3, the
potential temperature at a level just above the highest penetrating thermal s
49
-------
(i.e., at the lowest level where e is free of spatial and temporal fluctuations)
and 0, the mixed layer mean temperature, we propose that
A6 = FAh (4-34)
where r is the temperature lapse rate above z2,
r = 89/32 , z > zj (4-35)
and Ah is the distance above z2 that thermals penetrate.
To estimate Ah we assume that at the level z2 + Ah of maximum overshoot,
gApAh = pw2
where p is the density of thermals at level z2, Ap is the difference between
their density and that of the environment at z2 + Ah, and w is the effective
upward speed of thermals at z2. Approximating p/Ap by 0/A6, we get
A. e w2
Ah -
or
2 2
(AS) = row /g • (4-36)
In the case of a horizontally homogeneous atmosphere (4-28a) reduces to
»h w"2 - cr w_
• .' + ^cHp) (4-37)
and (4-30c) yields [see (2-21), (4-6) and (4-21)]
a w2 - w
^| (h0) - ac0 1 . 0C + tf awJ(p) [0 + (1 - ac)A9] + fgs (4-38)
C
where
f>V3 (4-38a)
16 • [re/g] 6ou (4-38b)
W
Here 8 is the ratio of the velocity w at H2 of the deepest penetrating
thermals and a. Clearly 8 > 1. Also, in deriving (4-38) we assumed that
W ~~
50
-------
the parameter c in (4-32) is unity.
We have performed preliminary tests of the model equations (4-37,38)
using data from the Wangara experiments for validation. These tests have
shown that the accuracy of the model is so sensitive to the choice of the
parameter T that accurate results cannot be achieved with c " • approximations
of this parameter, as originally hoped. Therefore, based o,. cne' findings
of Artaz and Andre (1980) who tested a number of mixed layer entrapment
models against observed data, we will adopt the simple empirical formulation
At this point we have developed methods for determining the values of
all independent variables in the expression for the flux F2 across H2 (4-21)
except the material concentration c~c entering cumulus clouds. We consider
this variable next.
Cumulus clouds are initiated by convective thermals, or updrafts, in
the mixed layer which in turn are generated and sustained by the surface
heat flux. We noted in the Introduction that the roots of the updrafts that
feed young cumulus clouds extend down into a relatively shallow layer of air
next to the ground. There they can collect fresh emissions of NOX and hydro
carbons and can transport them within a period of a few min into the clouds
above. This phenomenon is potentially quite important to the photochemistry
of the mixed layer because it is a mechanism by which fresh emissions are
somewhat selectively removed. It is also important to the photochemistry
of the atmosphere above the mixed layer because pollutants are injected
into this layer only after experiencing the saturated environment within
the clouds themselves.
51
-------
The depth of the cumulus roots, and hence the quantities of fresh
pollutants that can be removed, diminishes as the clouds grow. This occurs
because the shadows of the growing clouds attenuate the surface heat flux
until at some point continued cloud growth can be sustained only by the
latent energy of condensation within the clouds themsf as. If the atmosphere
is conditionally unstable and there is a sufficient supply "of moisture in
the mixed layer, cumulus'clouds eventually become self sustaining. Once
this occurs and the surface heat flux has been diminished to very low levels,
the updrafts feeding each cloud probably drain most of their air from upper
levels of the mixed layer, and with it more aged pollutants than those
found near the ground.
We will attempt to parameterize these cloud effects in our model by
allowing cumulus clouds to draw a fraction ip of their air from Layer 0.
The remaining portion will be taken from Layer 2. Keep in mind that con-
centrations in Layers 1 and 2 tend toward the same values during well mixed
daytime conditions. The fraction $ will be made a function of time to
reflect the changes (discussed above) that occur in the life cycle of the
cloud updraft. The form of the cloud updraft concentration parameter c~c
follows at once:
f »
(l-4/)2 . , in F9 9 formula > (4-39a)
^ ȣ
(l-*)2 + *[$£' + (l-C)c] , in £> , formula. (4-39b)
The term in brackets in (4-39b) is the concentration in Layer 0. It will be
discussed in detail later. The differences in the expression for c~c for
Layers 3 and 2 as indicated in (4-39) are the result of transport of material
from Layer 0 directly into Layer 3 without any intermediate mixing in
Layers 1 and 2. Thus, the net loss of material from Layer 2 to Layer 3 is
52
-------
as specified by (4-39a); but the gain by Layer 3 is the loss of Layer 2 plus
the portion of the loss of Layer 0 that enters cumulus updrafts. Further
consideration of this partitioning of material will be necessary later
when we derive F_ and f ..
0,0 0,1
Flux across the jet layer top: F and F.
i»i l ,£
In Section 3, we defined H. to be just above the tops of nighttime
urban plumes, typically about 300 m above ground-level. In the present
version of the model, H, will be maintained at approximately this same level
during the day [see Equations (3-5) - (3-7) for the mathematical definition
of H-], but eventually we plan to modify this formulation to allow Layer 1
to handle additionally the marine layers associated with lake and sea breeze
regimes.
A discussion here of some of the characteristics of the nocturnal
boundary layer that fall within the bounds of Layer 1 will help provide the
rationale for the flux formulation F that follows.
1, K
Figure 4-5 compares typical temperature mixing ratio and wind speed
profiles measured over rural and urban areas of St. Louis near dawn on
summer days. The top of the surface radiation inversion is at about 300 m
above ground level (AGL) at both sites. However, over the city there is a
shallow mixed layer approximately 100 m deep. At both sites the mixing
ratio decreases with height up to a point near the top of the inversion
where it becomes approximately constant. Its value at this point is com-
parable to that observed in the mixed layer the previous day. This suggests
that nighttime mixing of air between ground-level, the source of moisture,
and the upper levels of the radiation inversion is weak, particularly over
53
-------
(a)
(b)
~300m
\
\
\
T,q,lvl
T,q,M
Temperature
-._-._ Wind Speed
•—•— Mixing Ratio
Fi gure 4-5.
Schematic illustration of typical temperature, wind speed,
and mixing ratio .profiles over (a) rural and (b) urban
areas at night .
54
-------
rural areas where the largest'vertical gradients in the mixing ratios are
observed.
We have mentioned earlier that vertical fluxes of heat momentum and
drier air are often observed in rural areas at night during brief, sporadic
episodes of turbulence that are believed to be caused by the growth and
subsequent breakdown of waves in the nocturnal jet's shear layer. The jet
is evident in both sets of soundings in Figure 4-5 as the wind speed maximum
that occurs at about the top of the surface inversion layer. It is not
known whether these sporadic downward excursions of air are confined to a
shallow layer adjacent to the ground or whether they extend well up above
the jet itself. Telford et.a_l_. (1976) speculated on the basis of aircraft
observations that nighttime exchanges of air take place within layers
300 - 400 m thick at all levels up to about 1000 m. This is not supported
by vertical 03 profiles measured over St. Louis near sunrise (Evans 1979).
These data reveal the rather common occurrence of an 0, maximum in the core
of the nocturnal jet. If deep mixing of air occurs as suggested by Telford
£t !]_., 0^ maxima such as these could not persist through the night. One
of the goals of NEROS will be to determine the vertical extent of nighttime
mixing processes. This information is essential to the development of re-
liable regional scale oxidant models. Because if frequent exchanges of air
do occur between ground-level and layers above the top of the nocturnal in-
version, the maximum possible transport distance of 0, will be considerably
smaller than if deep mixing does not take place. Whatever the nature of the
mixing happens to be, we should be able to treat it in the model through a
suitable formulation of the flux term F .
55
-------
A feature of Hj not discussed in Section 3 that will affect the flux
formulation is that with z, defined by (3-7), it is possible in regions
of rough terrain for hills to extend up through the surface z = Zj.
Figure 4-6 illustrates a hypothetical case. We will denote by ay,(x,,t)
the fraction of the horizontal area A of the grid cell centered at (A,)
that is formed by the intersection of terrain with the surface z = Zi at
time t. This possibility requires that the definition of the top surface
H1 of Layer 1 be modified as follows:
Zi(A, ZT(A,) (4-40a)
ZT(A,<|>)-Z , otherwise (4-40b)
^
where Zj is given by Equation (3-7).
Using notation like that employed earlier in formulating F- , [see Equa-
c., K
tions (4-6) and (4-7)], we can express F, . in the form
1 ,K
dH,
* (4-40c)
where GI and Vj denote values measured on surface HJ, and the overbar
represents, as before, an average over the area A of a grid cell. To
evaluate (4-40c), we will divide the overbar averages into 2 parts: one over
the area ay, A intersected by terrain and the other over the remaining area
(1 - on)A.
Over the former, the first term in (4-40c) is the deposition flux onto
the ground surfaces that extend into Layer 2. We will parameterize this
flux by 62, where 6 is the deposition velocity of species c. Over this
same area the last term in (4-40c) is identically zero; because by virture
of (4-40b), H! is coincident with the terrain surface where it extends into
56
-------
_-2l(\.0.t)
•ZT^)
oT1A
gure 4-6. Intersection of the surface z = zi, defined by Eq. (3-5),
with elevated terrain features. A fraction a-r, of the
grid cell's horizontal area A is intersected by terrain.
57
-------
Layer 2 and hence there it is a material surface. Thus, the component of
F. . over the area Oy,A is
aTle2
(4-41a)
, if k = 1 .
dH,
Fl,k>l-aT1 - U - °T1^C1 -dT *
where Vj = (uj, v^, Wj). By design
dHl .
"dt ' nl
where n, is a given function in stable conditions (see 3-7a) and
nl = *1 ' WD1
in neutral and unstable situations. Consequently, (4-41b) becomes
l-o.
Tl
, stable ;
' (4-41c)
kwDl ' neutral
and unstable .
To approximate the terms in this equation that apply to unstable condi-
tions, we must take into account the effects of cumulus clouds that transport
material across H, from Layer 0 to Layer 3. The zones occupied by these
roots on surface H, contribute nothing to air exchanges between Layers 1
and 2. Given that the cumulus updraft velocity is w7, that a fraction a
w C
of the cell area A is covered by cumulus clouds, and that a fraction $ of
the cumulus flux is from Layer 0, we see that the "excluded" area for air
exchanges between Layers 1 and 2 across H, contains a total volume flux
*orw A. (4-4Id)
58
-------
We can express the total area of the excluded zones in terms of the
frequency distribution p(w) of vertical velocity on HI (i.e., all components
of wi except that induced by terrain w,.,).
If we assume that the horizontal area A of a cell is large enough that
spatial averages of velocity over it are equivalent to ~nsemble averages,
it follows from the definition of the probability dens.. p,-w) that the
fraction of A covered by vertical velocity in the range w to w + dw is
p(w)dw.
Keep in mind that we are concerned at the moment with unstable cases and
with the component of F, . through the area (1 - cr-r,)A of H, (see 4-41c).
Thus, it follows from (4-41d) and the properties of p(w) discussed above that
in the region (1 - a-r,)A of H, that contributes to E . , vertical velocity
WQ, is in the range - « < w £w , where w is the value that satisfies
wm
In summary, we assume that all regions of H! in which the instanta-
neous vertical velocity is greater than wm do not contribute to F1>k be-
cause they are part of cumulus updrafts that are transporting material
directly from Layer 0 to Layer 3.
Looking now at the first two terms on the righthand side of (4-41c) and
keeping in mind that WQ. includes turbulent velocity fluctuations, we have
GI*I - C!WDI = - dwR1 (4-42a)
where
wR1 = wD1 - 2l (4-42b)
59
-------
is the vertical component of the instantaneous fluid velocity relative to
surface Hj. As a first order approximation we shall assume that the concen-
tration c1 at any point on Hj is equal to : if the relative speed WR,
at that point is positive and that Ci = 2 if the local speed WD, < 0.
Kl —
That is
Ti , if wni > 0
(4-42c)
2 , if WR, <_ 0
This and the assumption that spatial averages over (1
to ensemble averages lead to
= 2
P(wR1)dwR1
— 00
Averaging (4-42b) we obtain
w
Rl
wR1 p(wR1)dwR1 = w
With this we can reduce (4-42d) to
= (2 -
- a-r,)A are equivalent
'wm *
WR1
wR1 P(wR1)dwR1
(4-42e)
wR1 P(wR1)dwR1]
(4-42f)
For simplicity we will assume that the probability density of WD, is Gaussian
with mean wQ1 and variance c2w,. Note that
P(WRI) = PD(WRI + 2J
where pQ denotes the density of w^. (see 4-42b). In this case we find with
the aid of (4-41e) that the term in brackets in (4-42f) is
60
-------
s z
and that
V2i
"1m
w
W Z.
erf
w,
__i
/2T
w2
[1 + erf
)]
(4-42g)
(4-42h)
where wD1 is given by (4-42e) and where
'Rl
rx
erf(x)
e"t2dt
(Under this definition, erf («) = 1.)
(4-43)
Applying a similar analysis to the stable case, where WR. then plays
the role of -ru (compare 4-41c and 4-42a) we find that in general E k can
be expressed in the form
Fl,k = aT!6k
-------
w
01
given , neutr . and unstable j
Wlm • stable
(4-44c)
and 6 =lifm = n;=0 otherwise.
mn
The variance a , of the vertical velocity component WD on surface Hj
is caused only by turbulence. Kaimal e_t at_. (1976) found from field studies
that in the daytime mixed layer,
aw2 = 0.4w*2 (4-45)
throughout the upper 80 percent of the mixed layer. [Here w* is the con-
vective velocity scale defined by (4-29b)]. Surface H. in our model lies
in this region most of the time.
Values of o at the level of H, are not as well known under stable
nighttime conditions. Pending the outcome of the NEROS nighttime studies
of vertical mixing, we will assume that under stable conditions at the
ground, there is no turbulence at the altitude of H, (approximately 300 m).
In summary, we will assume on surface H,
0.6 w*, L < 0
(4-46)
0 , otherwise
where L is the local Obukhov length scale.
Flux Across the Surface Layer Top: F and F ,
O 9 O U 9 Ji
Several factors complicate the form of F . . One is that in regions of
O 9 K
rough terrain, not all of surface H is common to both Layer 0 and 1; part
of it lies on terrain (see Figure 4-7). Another complicating factor is
that part of the volume flux that leaves Layer 0 may enter cumulus clouds
wl
62
-------
rather than Layer 1. A third factor is that to parameterize the subgrid
scale chemistry, which is one of the main roles of Layer 0, we require
different flux forms for each pollutant species. The last 2 factors will
be treated in the next section where we formulate the Layer 0 equations in
detail. Here we will develop the gross form of the F . expression and
0, K
consider the terrain effects cited above.
The general form of F.-. is the same as for all the other layers [see
0 » K
Equation (2-30)]
By design of surface H , the last term in this expression is identically zero
[see Equation (3-10)]. The remaining term is the total downward mass flux
of species c, across HQ per unit horizontal area.
Referring to Figure 4-7, we will assume that an area (a, - cr,.) of the
horizontal projection , of surface HQ lies on terrain. (Note that due to the
natural slopes of terrain, within any grid cell ay is always greater than
Assuming as before that the deposition flux can be approximated in
terms of layer-averaged concentrations, we have
Fo,k = 6k6kl + U -
where FQ is the mass flux through the portion of H common to both Layers
0 and 1, 6 is the deposition velocity of species c, and 6mn is the Kronecker
delta defined earlier following (4-44). In the next section we formulate
Fo'
63
-------
H
Figure 4-7.
To ' aTl)A
as
64
-------
SECTION 5
DERIVATION OF DIAGNOSTIC EQUATIONS FOR LAYER 0
The aim of this section is to develop a scheme for parameterizing the
physical and chemical processes cited in the Introduction that occur within
about 100 m of the ground. Among these are the effects on chemical reac-
tion rates of the spatial segregation of pollutant concentrations that re-
sult from the injections of material from point and line sources; the
stagnation of surface emissions in the calm surface layer of air over rural
areas and in mountain valleys at night; and the deposition of pollutants on
vegetation and other surfaces.
A basic constraint on the design of this scheme will be that it has
a diagnostic form; that is one that does not require computer storage in
addition to that taken by the other 3 layers of the model. (Appendix B
describes the constraints on the thickness of Layer 0 that allow the use of
diagnostic, rather than prognostic, equations to treat this layer.) With
this type of formulation, conditions in Layer 0 are expressed in terms of
current conditions in Layer. 1 immediately above and local, independent para-
meters such as emission rates, deposition velocity and the like. In essence
Layer 0 is simply a sophisticated lower boundary condition of Layer 1.
Indeed, since horizontal transport is negligibly small within Layer 0, it
is only through the vertical material flux F that the effects of events
in this layer are propogated to distant places.
65
-------
An added feature of the scheme we develop is that it provides an estimate
of the magnitude of spatial variations in concentration in Layer 0. This
feature is particularly valuable in comparing cell-averaged concentrations
that the model predicts with values observed at fixed point monitoring sites.
For simplicity we will formulate the scheme in My for the case in
which cumulus clouds are not present. Once this ba-.c "orm is established
it will be a relatively straightforward task to generalize it to include
cumulus cloud root effects.
The flux component FQ is given by [cf. Equations (4-47) and (4-48)]
dHo.
Fo " co dT = -co
-------
Line Source
Line Source Plume
Volume = 6V
Figure 5-1, Schematic illustration of Layer 0 and the notation used
to describe processes within it .
67
-------
and
x_ = 1 - \+ - (5-2b)
Let w_ be the mean value of |w | in the regions of descent (i.e., where
w < 0), and let w+ denote the mean relative speed of ascending fluid. In
terms of p we have
)ro(w)dw (5-3a)
wpro(w)dw (5-3b)
In a later section we will prescribe forms for p from which x_, x+, w_,
and w+ can be evaluated.
The area marked SA in Figure 5-1 is formed by the intersection of the
line source plume with the surface H . Since the source is at ground level,
fluid within 6A must on the average be ascending. We assume then that
C, 0AXA<|> (5-5)
In words, it is assumed that on the surface H the plumes of all
sources in Layer 0 are confined to zones of ascending fluid. Once these
emissions reach the cells above the lowest layer, thorough mixing of all
pollutants is assumed to result so that in all fluid descending through
HQ concentrations are uniformly distributed.
Fluxes through the side walls of cells within Layer 0 are negligible,
as a consequence of our definition of H [see the discussion proceeding
68
-------
(3-8)], compared to those through the surface HQ and will be neglected.
Within each cell of Layer 0, we distinguish between 2 domains: 6v,
the volume occupied by all emissions of sources in that cell and V - 6v,
the region outside
A+AA
max{[z0(Al,(f.l,t) - ZT(A',<(>')],0} (5-6a)
A-AA
It is convenient to define also an effective cell depth
AZQ E V/AQ (5-6b)
and
(5-6c)
0
Within the restrictions of assumptions introduced above, the only source
of fluid in V - Sv is subsidence through the surface H . However, there
are 3 avenues by which fluid leaves this domain. One is by direct trans-
port back though the cell top, that is, through an area (1 - c)A+A.
A second path of escape is turbulent entrainment into the plume volume dv.
Finally, pollutant species can leave V - 6v by deposition on the ground
surface. Chemical reactions within V -
-------
as 0, and NO, chemical processes within 6v are expected to be rapid because
the concentrations of these pollutants are greatly disturbed from their
chemical equilibrium values by the fresh injections of NO . We believe that
/\
0.,, NO and N09 are the only species that are sufficiently reactive to be
*j £
affected by subgrid scale concentration variations (SSV), as it exists
within 6v. Accordingly, we shall assume in this study that only Og and NO
undergo significant conversion within 6v and that, moreover, the rate of
reaction between these species is infinite. All other pollutants will be
treated as pseudo-inert substances inasmuch as the time scales of their
reactions are much larger than the characteristic residence time of material
within Layer 0.
Our plan now is to write the mass balance equations for each of the
pollutant species in sv and V - Sv and, from these, to formulate the desired
expressions for the net material fluxes through H and the mean pollutant
concentrations in Layer 0. For these purposes, it is convenient to introduce
the following special notations:
*
M(xO = net downward mass flow rate (mass/time) of given
species x across H into V - i = average concentrations of species x in Layer 1
in the cell atop the given cell of Layer 0,
70
-------
x = average concentration of species x in V - <5v,
X1 = average concentration of species x in <5v.
Many of the terms listed above are illustrated in Figure 5 1.
A. Mass Balance Equations in V - 5v
Consider first the mass balance of CL. We have
M(03+) = M(03^) + M(03i) (5-7)
Straightforward analysis gives
•
- 03w+(l - c)A+A (5-8)
(5-9)
In the last equation, v is the entrainment velocity of fluid into 6v and
This is the fraction (0 <_ ? < 1) of the lowest cell occuppied by source
emissions. (Note that c t «v/V unless OT =0.)
o
For the deposition mass flow rate of 03 in V - 6v, we assume
•
MCOgi) = e0 A03(l - 5) (5-12)
where $n is the deposition velocity of Ov The factor (1 - c) in Equation
U3
(5-12) is intended to account for the fact that the domain V - 6v does
not necessarily cover the entire ground surface. We assume that a fraction
71
-------
C of it is covered by 6v.
Substituting Equations (5-8), (5-9) and (5-12) into Equation (5-7)
and making use of the relationship
vc = w+A+£ (5-13)
imposed by conservation of mass, we obtain
w_A_i
By analogy with this relationship, we obtain
w A i
NO = w x + (! . g)e (5-i5)
w A !
N02 - w A ; ~(l /g)B (5-16)
and in general
w A i
B. Mass Balance Equations in 6v
For 03, we have
• * • •
M(03-^) = M(0^*) + M(O^t) + M(OJt) (5-18)
The first term of this equation was given earlier:
*
M(03+T) = v6a03 - (5-19)
The deposition term is given by Equation (5-6)
•
M(OJt) = 6Q cA03 (5-20)
and the upward mass flow by
) = W+A+C03A (5-21)
•
The chemical removal term M(03+*) is somewhat more complicated. Under the
72
-------
assumption adopted previously that NO and 03 react infinitely rapidly in
i) at which it is entrained into 6v no matter
how large the NO emission rate might be. Let SNQ reprecnnt the NO emission
rate averaged over the lowest cell in units of moles/are ... Then from
the considerations cited above, we conclude
03v6a if
SNQA otherwise
(5-22)
Substituting Equations (5-19) through (5-22) into Equation (5-18) and simp-
lifying, we obtain
0 if SNQ > vc03
(5-23)
n
°3V '
Vv '
otherwise
The mass balance equation for NO in 5v is
= M(NO'+*) + M(NO'i) + M(NO'i)
(5-24)
All of the terms here can be written down immediately based on information
supplied above:
MfNG-^) = NOvSa , (5-25)
M(NO'-h) = NO'w+X+cA (5-26)
•
M(NO'+) = 6wnNO'cA > (5-27)
if
S..QA otherwise
(5-^8)
73
-------
Combining these expressions with Equation (5-24), we obtain
NO =
-^+v(NO-OJ
S -5 J jr f
if
vNO
NO
, otherwise
The mass balance equation for NCL is
M(N0->I)
A -i-
Note that we have neglected the loss of NOo due to photolysis.
procedures like those used above we obtain
(5-29)
(5-30)
Following
if
SN09 + SNO
, otherwise
(5-31)
One of the assumptions we made earlier is that all pollutant species
except 03 and NO react too slowly to undergo significant chemical trans-
formation during the brief time they are in the lowest layer (Layer 0).
Let x' represent the mean concentration of 1 of these species in 6v.
Its mass balance equation then has the form
0 (5-32)
+ SYA = M(x'±)
A.
which reduces after analyses like those performed earlier to
6 + v
X
(5-33)
74
-------
Recall that x can be obtained from Equation , j-17).
C. The flux component F in the case of no cumulus clouds
We can now formulate the flux component F that enters in the general
expression (4-48) for the total flux Fn , . Let F (0,) be the net down-
O, K 0 o
ward flux of 03 through HQ over a horizontal area A. That is,
• •
FQ(03) = £ [M(030 - M(03t)] (5-34)
which reduces with the aid of (5-8) and (5-21) to
F (OJ = w X <0,>, - 0?w,x, (1 - 5) - w.X.cO; (5-35)
UO "*""O1 *3'' '«j
Making use of Equations (5-14) and (5-23), we obtain finally
Fn(0j = w X <0^>, - O^w.x, (1 - 5)
O*5 ~™Ox «ji^"
* ^
0 if SNQ > v;03
(5-36)
03v2; - SNQv
otherwise
Q.
Vv
where
w x <0o>,
(5-37)
Note that Equations (5-36) and (5-37) -give the flux in terms of the predicted
mean concentration <03>, in the cell above Layer 0, the NO emissions rate
in Layer 0, and other quantities, all of which are specified.
With the NO flux F (NO) defined in the same manner as Equation (5-34),
we obtain through a process similar to that leading to Equation (5-37):
75
-------
F (NO) = !W_X_
where
- 03)
if
1 +
5 NO
t otherwise
(5-38)
NO =
'NO
(5-39)
and 03 is given by Equation (5-37)
The N02 flux is found to be
(1-
vc(NO
^2 if§NO>vc03
9 + Swn + Swn
2 N02 N0 . otherwise
3NO,
1 + —
V V
(5-40)
where
(5-41)
76
-------
and 03 is given by Equation (5-37).
Finally, for all pollutants x other than NO, 0,, and NOo, we have
F0(x) = iw.^_
"(1 - c)6y *
[d - c)ex *
w+X+s
W X
v£x + S
X
6
1 + _2L
•\ \
where
lW_X_
V+
(5-42)
(5-43)
D. Mod.ifications of F,, to include effects of cumulus clouds.
When cumulus clouds are present overhead, they pull a volume flux
from Layer 0.
Fcu = *°cvcA
Here
V = -
- w
1 -
- fe/A9
(5-44)
(5-45)
[see Equations (4-15), (4-20) and (4-20a)]. The upward volume flux through
the top of Layer 0 is
w+X+ (1 - aTo)A (5-46)
This follows from the definitions of w+ and x+. Obviously, the cumulus
flux F cannot exceed that given by (5-46), so we must limit the values
we assign to \i> by
JTo-
(5-47)
We must impose this restriction on ^ at all times and at all points in
space during model simulation.
77
-------
Figure 5-2 illustrates the partitioning of the cumulus cloud flux.
Of the upward volume flux through H a fraction a given by
a =
(5-48)
is destined for cumulus clouds. The remainder enters Layer 1. Thus, when
clouds are present, the flux component F , which is the net downward flux
into Layer 0 from Layer 1, has the form
•
F0(x) - £ (M(x+) -
- Mc(x+)] - CM(x't) - Mc(x't)]} (5-49)
where M represents mass flux into clouds. Modifying Equations (5-36} (5-38)
and (5-40) in this manner we obtain the final forms of the flux component
F for use in Equation (4-48):
<03>lW-X- ' °3W+X+
- a)
-(1 - a)
0, if SNQ > v?03
" SNOV , otherwise
where 03 is given by (5-37);
FQ(NO) » iW_X_ - NOw+X+ (1 - 5) (1 - a)
-(1 - a)
NQ
- 0)
, if S
NO
, otherwise
where NO is given by (5-39);
(5-50)
(5-51)
78
-------
\KVCA;
) - Total upward,
volume flux' '
= w+X+(1-ato)A
Figure 5-2.
Partitioning of surface flux into cumulus clouds and
Layer 1 .
79
-------
FQ(N02) = iW_X_ - N02w+X+ (1 - 0 (1 - a)
-(1 - a)
2 + °3> + SNO,
if S
NO >
(5-52)
where N02 is given by 5-41;
F0(x) = iw_x_ - xw+X+ (1 - 5) (1 - a)
(5-53)
- (1 - a) (v£X + S
where x is given by (5-43).
6 /v)
A
-1
E. Equations for the cell averaged and fluctuation concentrations in Layer 0.
The average concentration of species x over the volume V of a cell in
Layer 0 is
o = V
(5-54)
where V is given by (5-6a). The assumed uniformity of concentrations
within each of the two domains V - 5v and 6V of each, cell leads to
'- (1 - c)x + ex1
(5-55)
where
Using this equation and those given earlier for 03§ Oi, NO, and so on, we
obtain the following:
80
-------
if
°
1 -
otherwise
where CU is given by Equation (5-37);
(5-56)
Q = NO
1 -
+ -^
where NO is given by Equation (5-39);
PNO T v
0 otherwise
if
Q = N02
5NO,
+ v
(5-57)
If
N0
otherwise
and N02 is given by Equation (5-41);
iW_X_[v + (1 - c)B ]
6 + v
X
(5-58)
(5-59)
As we noted earlier, it is not necessary to evaluate the concentration
formulas [Equations (5-56) through (5-59] at every time step and in every grid
81
-------
cell. In the computer program these equations are coded as a subroutine
that is called only at those times and sites where ground-level concentration
estimates are needed. Concurrent evaluations of the root-mean-square con-
centration variation a within selected cells of Layer 0 are also performed.
A
In terms of the concentrations x and x1 defined earlier, a is given by
A
- X')2 (5-60)
A
This expression can be applied to any of the pollutant species. For example,
the rms ozone concentration can be obtained by substituting (5-37) for x
and (5-23) for x' in (5-60).
No other air pollution model predicts quantities like a , but information
/v
about mean concentration variations within grid cells, such as a provides,
X
is of value in estimating concentration extrema, in assessing the magnitude
of modeling errors from comparisons of predicted cell averaged concentrations
with values measured at fixed monitoring sites within those cells, and in
studies where subgrid scale concentrations fluctuations play an important
role.
82
-------
SECTION 6
THE REPRESENTATION OF ATMOSPHERIC MOTION: WINDS AND TURBULENCE
In classical studies of short-range dispersion, atmospheric motion is
assumed to consist of 2 rather distinct components: a large scale, uniform
(in space and time) flow (the so-called "mean" or "transport" wind); and
small scale spatial and temporal fluctuations in the flow known as turbulence.
The latter are generated by 2 physical processes: the transformation into
kinetic energy of potential energy acquired by the fluid when its lower
boundary is heated (so-called buoyancy generated turbulence); and small scale
motions associated with the instability or rolling-up of vortex sheets pro-
duced on rigid surfaces by viscous effects (so-called shear generated tur-
bulence). It has been found from observations over flat, uniform terrain that
if the "mean" wind is represented by a time average of the instantaneous
velocity over a period of about 30 min, that time averaged variances, covari-
ances, etc., of the corresponding turbulent velocity fluctuations acquire
universal forms when properly scaled with the speed of the "mean" wind, the
depth h of the fluid in which turbulence is confined, and properties at the
lower boundary, namely the surface roughness and heat flux.
It is perhaps these universal characteristics that have given rise to
the common notion that turbulence, and turbulent diffusion in particular, are
intrinsic properties of atmospheric flow. This premise underlies the common
practice in air pollution models of representing the mean or transport wind
83
-------
by a mathematical expression fit to hourly averaged wind observations and
parameterizing the effects, of small scale wind fluctuations by dispersion
coefficients or an eddy diffusivity whose values are taken from empirical
data. One of the purposes of this section is to emphasize that this is not
a sound method of representing atmospheric motion in regional dispersion
models. Our principal objective here will be to formulate a conceptual
and mathematical basis for representing fluid motion that not only is con-
sistent with the classical definitions of transport and turbulence asso-
ciated with boundary layer diffusion studies at short range but also is
applicable to simulations of dispersion over arbitrarily large scales.
Consider for a moment the universality of the structure noted above
of turbulence defined in the classical sense. When a fluid of given depth
is forced by some uniform, external agent to flow over an "infinite", flat,
uniform surface of given roughness and temperature, the states accessible
to the velocity, pressure and temperature of the fluid comprise only that
limited region of the phase space of these fields whose members satisfy all
the external constraints and physical laws jointly. It is in part the
limited domain of accessible states that is manifest in the tendency for
time-averaged quantities in idealized boundary flows to exhibit universal
forms; and in essence it is the quantification of -the properties of those
states that is the objective of boundary layer turbulence studies.
Now in regional scale fluids the physical laws involved are the same
but the external constraints are different. The driving force of the fluid
is not uniform or steady; the bounding surface is not infinite, flat, or
uniform; the effective depth of the flow is not fixed; etc. Under these
conditions the accessible states of the fluid motion comprise a much
84
-------
different region of phase space whose members .an not be characterized solely
by properties of the states accessible to the idealized boundary flow.
The existence of a broader set of possible flow fields in these conditions
is evident in a superficial way in the fact that a set of discrete, wind
observations collected over a regional scale area do not uniquely define
a flow field. Indeed, when one attempts to fit a continuous field to dis-
crete data using some objective analysis technique, one finds that each
technique produces a different field. The ability of any technique to
fit a single function to a given set of data is a result of the imposition
by that technique of additional constraints on the system. In most cases
these added constraints are artificial and have no relevance to the "physics"
of the system. One also finds that the flow field derived from a set of
data varies as the number of locations of the observation sites changes.
In short, the mean or transport flow fields commonly used in regional scale
diffusion models are not unique.
This brings into question the practice of parameterizing the disper-
sive effects of small scale velocity fluctuations using empirical expressions
that are completely independent of the data used to represent the "mean"
wind. If the flow field is decomposed into a "mean" and a "turbulent" part,
then for a given flow the character of the "turbulence" will change if one
alters the description of the "mean" flow, as, for example, by the incor-
poration of additional wind data into the analysis. Recent studies such as
that of Carras and Williams (1981) have measured the spread of point source
plumes out to distances of 1000 km from the source and expressed the spread
as a function of travel time or distance. Many long range transport models
utilize these data to estimate the width of a plume, and flow fields fit to
hourly averaged wind data in the model domain to estimate the center!ine.
85
-------
This practice is simply an extrapolation of the concepts employed in the
empirical Gaussian model of short-range dispersion, and as we shall show,
it is not supported by theory.
That a basic problem exists in the current methods of treating diffusion
in regional scale simulations is evident in the following paradox, of which
apparently few modelers are aware.
In implementing K-theory in long range diffusion models, Taylor's
(1922) statistical theory of turbulence is usually invoked to support the
premise that for large-scale diffusion the eddy diffusivity is proportional
to the variance of the turbulence velocity, u1* say, times the Lagrangian
integral time scale T, of the turbulence, thus,
K * TL IT2" (6-1)
Since T, is not an easily measureable quantity, its value is usually inferred
from diffusion data [see, for example, Draxler (1976)]. However, the inte-
i
gral time scale T, of a stochastic process is nonzero only if the energy
spectrum of the process has finite amplitude in the limit as the frequency
u •*• 0 (we show this in Section 7). But the zero frequency component of at-
mospheric motion is contained not in the "turbulence" but in the time aver-
aged wind u" which is used to describe the "mean" flow. Consequently, the
integral time scale of "turbulence" as it is defined implicitly in regional
models should always be zero! The existence of this paradox is an indicator
that basic features of existing diffusion models are not well defined.
The questions and problems that we have cited above have come to the
forefront as the scope of diffusion models has been expanded from the micro
scale to the regional scale, and beyond. Short-range diffusion studies have
been based heavily on empiricism and in this way theoretical questions of
86
-------
definition have been largely avoided. For example, the Gaussian plume for-
mula has been the mainstay of short-range dispersion analyses for many
years. Being an empirical relationship (based largely on diffusion data
collected during the Prairie Grass experiments), the plume formula provides
quite well defined concentration estimates (approximately 10 min averages
at fixed points in space) in terms of well specified inp - -ameters.
Unfortunately, there does not exist a comparable empirical basis for regional
scale modeling, although we might add that there is some weak evidence that
such a basis might actually exist, at least for chemically inert material.
Specifically, the success of the Gaussian model in predicting 10-30 min aver-
aged concentrations out to 1 km from a point source is apparently due to the
small scatter in 30 m averaged concentrations that one finds from one obser-
vation to another conducted under the same flow conditions. And this small
scatter is apparently attributable to a gap in the turbulence energy spectrum
(Van der Hoven, 1957) at periods of about 1 h. It turns out that there is
another gap in the spectrum at periods of the order of 90 days, and this
suggests that perhaps there exists a regional scale counterpart of the em-
pirical Gaussian model that would provide estimates of seasonally averaged
concentrations out to 1000 km from a source in terms of certain seasonally
averaged meteorological fields. Demonstrating the viability of such a
model would require a large scale dispersion experiment conducted over a
period of at least 1 year.
It is imperative now that unresolved theoretical aspects of the repre-
sentation of atmospheric motion in diffusion models be addressed. We attempt
that in this section. By this effort we hope to derive sound definitions
of transport and turbulence and at the same time lay the basis for
answering the following closely related questions which have great relevance
87
-------
to the utility of pollution models in regulatory decision-making roles:
(1) What aspects of pollutant concentrations are models capable of
predicting?
(2) Given our present state of knowledge, what are the theoretical
limits on the accuracy with which these quantities can be pre-
dicted (assuming perfect emissions data, chemical information,
etc.)?
(3) How does one assess the accuracy of a given model?
A Game Analogy of Diffusion Modeling
In order to help convey some of the basic ideas that we introduce later,
we present in this section a simple game analogy that possesses many fea-
tures in common with the diffusion modeling problem but which are less
abstract. Our game involves the throwing of 3 dice 3 times to produce a
3x3 matrix of, for example, numbers A. (We consider 3 throws of 3 dice,
rather than a single throw of 9 dice, as a way of retaining space-time fea-
tures). We assume for simplicity that the faces of each of the 3 dice are
marked with an integer from the set 1, 2, 3, 4, 5 and 6. In this case the
g
matrix A is a member of the family of 6 * 10,077,696 matrices illustrated
in Figure 6-1. We will call this the global family of matrices because it
contains all matrices that are consistent with the fundamental definition of
the system: 3 dice with faces marked 1, 2, 3, 4, 5 or 6 thrown 3 times.
Every realization of this system must be a member of the global family
shown in Figure 6-1, regardless of whether the dice involved are "loaded",
or whether some were accidentally stamped with the same number on more than
one face, or how the dice are thrown, etc. These details are characteristics
-------
n=1
= 10,077.696
Figure 6-L The global family of possible outcome matrices A
resulting from the throw of three dice three
times (see text for more details) .
89
-------
of the specific system (i.e., the actual cice used in a particular instance,
the throwing method, etc.). For reasons that will become clear later, to
maintain a proper analogy with diffusion modeling we cannot restrict specific
details of the system beyond the 3 basic requirements given earlier. In
other words, in our game the choice of the dice, the throwing method, etc.,
are optional, so long as there are only 3 dice thrown 3 times with faces
marked with 1 of the integers, 1, 2, 3, 4, 5 or 6.
Suppose that a set of dice has been selected and is thrown. We ask
is it possible to predict in advance the outcome matrix A? The determinism
that we associate with physical law leads us to believe that we could pre-
dict A if we knew all the physical details of the dice, their exact positions
at the beginning of the throw, etc. But like so many natural phenomena,
these details are so numerous and intricate that they are outside the scope
of our capabilities of quantitative observation. And without precise
knowledge of the initial state of a system, physical laws are powerless as
predictors.
This is one of the key points that we wish to emphasize in this section.
We believe that one of the debilitating practices of atmospheric diffusion
modeling as it is now applied to large scales is the employment of physical
laws in deterministic roles when the observational data are too crude to
warrant it. In the next subsection we argue that under these conditions, a
more natural role for physical law is that of delineating the class or
family of possible atmospheric states (such as the family of A matrices
defined above). We will postpone further discussion of this point until
later in this section.
90
-------
Lacking ability to predict the specific matrix A that will result
on any outcome of our dice game, we next pursue the question of whether
there is not some attribute of each A that is so prevalent in the family of
A's that nearly all members possess it, or at least a value near some one
value. In seeking such an attribute, we consider only those that are analo-
gous to properties of interest in the diffusion modeling counterpart that
we introduce later. Here the space-time character of our A matrix will be
of value. Let's assume that each column in A represents the numbers on each
of the 3 dice on 1 of the 3 roles and is roughly analogous to 3 spatial
points at a fixed time. Similarly, each row of A contains the 3 outcomes
of each of the 3 successive roles and we will regard it as analogous to
temporal variation at a fixed space point.
Within this frame of view the attribute "the sum of the 9 elements of
A" is the analogue Qf a space-time average of a single realization or out-
come of A and therefore it is a property of interest. Suppose we scan down
the family of An shown in Figure 6-1 and compute the sum Sn of the 9 elements
of each member, viz.,
s"'A A i* (6-2)
where a.. are the elements of An.
We find that the sums Sn lie in the range 9 < S < 54 for all n. Only
n ™~ n ^~
1 of the A has the sum 9, namely the member n=l; and only the last member
has the sum 54. Nine members have the sum 10 and 9 possess the sum 53.
But an inordinate 767,394 have the sum 31 and the same number have the sum
32. In fact, over one-half of the entire family have sums that are within
±10 percent of the family mean sum 31.5; and almost 90 percent have sums
91
-------
within ±20 percent of the'family mean. T .= full distribution of sums Sn
is plotted in Figure 6-2. (We add in passing that since the profile shown
in Figure 6-2 is simply a property of the set of numbers illustrated in
Figure 6-1, we can argue that the Gaussian shape of this profile is attri-
butable to the "geometry" of the number system. Perhaps this is why suffi-
cient conditions for the occurrence of the Gaussian distribution, as set
forth by the Central Limit Theorem, are so weak.)
The tight clustering of the sums Sp about the family mean value suggests
that even though we cannot predict which of the matrices A will occur on
any occasion, we can predict with respectable accuracy what the sum of its
elements will be, provided that there is nothing inherent in the dice, the
throwing method, etc., that causes preferential occurrence of those matrices
whose sums are farthest from the family mean value. Thus, we are back to
the problem of predictability or, more generally, causality. What we lack
is information about the "mechanism" that selects from the family of A
matrices the ones we actually observe when a given set of dice are rolled.
Philosophically, one might argue that there is no causal "mechanism",
that the outcomes of A that we observe are simply revelations of the unfolding
of time. However, this point of view does not aid our understanding or
conceptualization of the problem so we will postulate a "mechanism" that is
consistent with observation. Specifically, we imagine that for every set of
dice, throwing method, surface, etc., there exists a well balanced wheel and
pointer with the family numbers n shown in Figure 6-1 inscribed on the cir-
cumference of the wheel. Before each throw (i.e., 3 throws of the 3 dice),
the wheel is spun (by "mother nature") and the number n indicated by the
pointer when the wheel stops determines the matrix A that we observe on
that throw.
92
-------
Number of A with given S
770,000
760,000
750,000
740,000
30,000
20,000
10,000
-I
I 1 1 1 I »*T 1 1 1
ho
1
1 I 1 1 1 1 1 ] J
20
l t 1 1 l 1 l l i i i l l I l i i l i i l i rJtl l l
30 40 50
33 _
Figure 6-2,
Number of matrices in the global family with given
element sum S .
93
-------
Now there are many ways to arrange the numbers n on the wheel . The
simplest is to include all the numbers, n=l,2,... 10,077,696 and to allocate
equal angles of arc to each. This configuation is the one we would associate
with "fair" dice games in which any one of the possible outcomes (A) is as
"likely" to occur as any other. Another possible arrangement of the numbers
is to include them all but to allocate more space to sor. ' jers than to
others. This would represent a game in which the dice were "loaded". Yet
another possibility is to exclude some numbers n altogether. For example,
if due to some manufacturing flaw 1 of the 3 dice had had a "5" stamped on
the face that should be marked "6", our corresponding "causal wheel" would
not contain the numbers n of Figure 6-1 that denote matrices with 6's in
row 2, say; and those numbers that represent matrices with 5's in row 2
would have twice the space allocations as the others.
Our wheel concept is quite crude but we believe that it provides an
easy way to conceptualize the link between the purely mathematical aspects
of the problem, such as the family of possible states of the system (or
sample space as it is usually referred to in probability theory), and the
sequence of system states that we would actually observe with a particular
system (e.g., set of dice, throwing method, etc.)- Also, in the context
of the wheel the probability of occurrence of a given matrix An is simply
Prob(Ap) = •£ (6-3)
where x., is the arc length in radians allocated to member n.
n a
In statistical analyses one deals with ensembles of the function or
variable in question. We can produce an ensemble of A's by spinning the
wheel many times and listing the A-s that result. One of the important
94
-------
purposes that an ensemble serves is to provide information about which
members of the global family of states is included on the wheel and the
relative space allocated to each (i.e., their probabilities).
In some instances we do not need to know anything about the "causal
wheel". For example, if we are designing a system and are concerned about
the occurrence of certain states that violate some safety criteria, and if
we are able to define the global family of possible states as we have done
in Figure 6-1 for the dice game, then there is no need for information about
the wheel if the states of concern are not members of the family. In most
cases, however, the violating states are members of the global family and
the problem is then one of estimating how frequently they will occur during
operation of the system. (The ultimate goal of the present study is the
estimation of the frequency of violation of air quality standards.) In
these situations information about probabilities is essential.
In this connection 2 important points must be made. The first is
that theoretical knowledge is generally incapable of describing exactly
which family members are on the wheel and how they are apportioned space.
We can only infer this information from observations, namely from an en-
semble as defined earlier. The second point, which is somewhat contradic-
tory of the first, is that the exact collection and layout of numbers on
the wheel cannot be determined from observation.
To see this, suppose we had a perfectly balanced wheel with the numbers
n=l,2,... N=10,077,696 arranged at equal intervals around its rim. Since
there exists no force to cause the wheel to stop on each of the N numbers
once and only once in N spins, an ensemble of N numbers produced by N spins
might not contain some of the n. Obviously, we would be wrong if we inferred
95
-------
from this ensemble that the missing n were not present on the wheel. It is
this same freedom of the action of the wheel that makes it impossible to
predict the frequency with which given states of a system will occur.
For example, our intuition tells us that throwing 9 dice and having
them all come up 6's (the last member of the A family of Figure 6-1), is
a "rare" event. Yet if the dice are "fair", this event is just as likely
as any of the others. Thus, all of the outcomes are equally rare so a "rare"
event occurs on every outcome of the throw! It follows that there is no
reason for "rare" events to occur infrequently. In the language of hydro-
logists, the "hundred year flood" may occur several years in succession.
Consequently, all we can do is base decisions on the expected frequencies
of events, and then to hope that we're lucky!
The expected frequency of the all-6 matrix of Figure 6-1 in a fair
game is 1 in 10,077,696 games. If, in some hypothetical game, the player
t
won $10 on every throw that the all-6 matrix did not come up but was beheaded
the first time it did, one could "expect" to make a million dollars a year
at this game for 100 years before suffering the ultimate loss. But there
is nothing in the universe to prevent the all-6 matrix from coming up the
very first game! Similarly, we can "expect" the sum S of the dice to differ
from 31.5 by no more than ±20 percent in only 1 game in 10. But there is
nothing to prevent exceedance in 10 successive games. Our purpose in elabo-
rating on these basic concepts is partly to emphasize the limits on the know-
ledge that we can hope to acquire about the future behavior of natural systems,
We especially want to emphasize that the notion of determinism that under-
lies current thinking and practices in air pollution modeling and in the
expectations that regulatory officials have in the utility of model pre-
96
-------
dictions is unsupportable. Another purpose is to lay the bases for the
general approach to diffusion modeling that we present in the next subsection.
Before proceeding we need to introduce a few more concepts and then
to summarize all that we have presented in this section.
We first define a subset of the global family of states as the collection
of all members that have some given property in common. For example, we
can speak of the subfamily of A's that have 1's along their main diagonal
(there are 6 members of this subfamily); or the subfamily whose element
sum S = 18 (there are 22,825 members of this subfamily); etc.
The ensemble mean of a parameter is obtained by generating an infinite
ensemble of system states (e.g., A's) in the manner described earlier, and
averaging the parameter values collected from each member of the ensemble.
We denote ensemble averages by angle brackets <>.
Intuitively, one can see that the ensemble mean can be computed from
the global family by weighting the parameter value associated with each
member of the family by the probability of that member, as defined in
Equation (6-3). For example, the ensemble mean of the matrix element sum
Sn discussed earlier is
10,077,696
= E Prob(AjSn " (6-4)
n=l ~n n
(In the remainder of this report the ensemble average <> and the cell average
<>. can be distinguished by the presence or absence of the subscript.)
J
If S represents the sum of the elements of the matrix A that occurs on
a given throw, then as we have already discussed, S is not predictable but
is. Thus, it is of interest to know how far to "expect" S to fall from
97
-------
the ensemble mean value . Denoting the difference by
E = S - (6-5)
we see that the ensemble mean square value of e is
10,077,696
= I Prob(An)(S_ - )2 (6-6)
n=l - n
If 2 is a small fraction of , then within the ensemble the values of
S are clustered around and we can expect values near to occur fre-
quently. (The distribution of Sp shown in Figure 6-2 is an example.) This
measure of the range of likely values is one of the easiest to determine
and is used widely in statistical analyses.
In this section we have distinguished between 2 inherent aspects of
natural systems: the global family of possible states of the system, and
the "causal wheel" which we shall use here to conceptualize the agent that
causes particular states to occur on particular occasions when observing a
particular system. A precise description of the global family (or sample
space) can often be derived from the fundamental features of the system.
We did this earlier for the case of a system of dice. However, a quantitative
description of the causal wheel is generally not derivable and can only be
inferred from actual observations of the system's behavior. The property
of concern is embodied in the concept of probability, which we defined in
Equation (6-3) in terms of our wheel. Having defined the global family and
estimated the probability of its various members, we can then calculate en-
semble mean values as in Equations (6-4) and (6-6) above. If our description
of the global family and its member's probabilities are accurate, we can ex-
tract certain ensemble statistics from this information that will delineate
bounds within which given system's parameter values can be expected to fall
98
-------
with given frequency. These expected freqjencies represent averages over
an infinite period of observation and may differ greatly from those seen
during a finite period of observation, even when the information used to make
the expected frequency estimates is exact. The important point is that the
parameter value itself is not predictable. The best that can happen is
that the interval over which 90 percent of the values are expected to occur
is so narrow that, on average, 90 percent of the observed values differ
negligibly from the prediction.
The Statistical Bases for Air Pollution Modeling
The quantities of pollutants and their precursors emitted into the
atmosphere and the times and places where they are released are variables
that, for the most part, are controlled by man. In fact, it is this con-
trol that is to be exercised to achieve concentration levels that are in
compliance with air quality standards. But once material is airborne, its
fate is determined jointly by atmospheric motion, chemistry and other pro-
cesses that are neither controllable nor observable (quantitatively) in
their entirety. Like the dice discussed earlier, atmospheric motion obeys
known physical laws but the information needed about the state of the atmos-
phere to apply these laws in a predictive role is far beyond our capabilities
to acquire. Consequently, the joint effects of atmospheric motion, chemistry,
etc., on pollutant concentration must be treated in a statistical manner like
that employed earlier with the dice game. Our first step, then, is to define
the global family of possible atmospheric states.
99
-------
The Global Family of Atmospheric States
1. General considerations
Let P denote the spatial region in which pollutants are to be modeled
and let T represent the period of interest. Let u(x,t) represent the in-
stantaneous fluid velocity at any point (x,t) in spac 3. Then, by ex-
pressing atmospheric motion in the truncated form
A
u(x,t), if x e V and 0 <. t <_ T
u(x,t) -
0, otherwise
•
we can represent u by the complex Fourier transform
(6-7)
u(x,t) =
1
*n+l
(6-8)
where n is the spatial dimensionality of the modeling region (n=l,2, or 3);
(1 is the complex Fourier transform of u; k is the wave number vector; u is
the angular frequency; and the integrations are over all k and u space.
By definition of U the inverse transform of (6-8) is
«•
r rT
ut)
dtdx
(6-9)
V o
Consider for the moment the case where V is one dimensional. Equation (6-8)
reduces to
u(x,t) =
U(k,u)ei(kx * ut>d«dk
(6-10)
where now u, Ut and x and k are scalar variables. We can imagine the
integrals in (6-10) to be the limits of sums, namely
100
-------
BO OB
u(Xft) = — L- lim Z I i,ne1(knX + wmt) ^^ (6-11)
(2ir)2 AkAurK) n=-» m=-» mn
where k = nAk, to = mAu and U are elements of the matrix U. In words,
given a one dimensional domain V and time period T, we can represent any
A
velocity fields u(x,t) in that region by a matrix (l usir , ^tion (6-11)
and its 1-D counterpart, Equation (6-9). For each possible velocity
distribution u(k,t) there exists a matrix U, and conversely.
A
Note that the matrices U are analogous to the A used in the previous
section. .to represent the possible states of the system of dice. Our reason
for working with U rather than the velocity itself is that the former pro-
vides a much more concise way of representing continuous functions, especially
when they are to be approximated by a finite set of discrete data.
In order to illustrate important points, it is necessary to work with
at least 2-D fluids. Therefore, considering a two dimensional flow, we
adopt the notations
u(x.t) = [u(x.t), v(x,t)] (6-12a)
x - (x,y) (6-125)
))] (6-12c)
k - (kx, ky). (6-12d)
Then, from Equation (6-8) we have
. f f i(k-x + ut)
u(x,t) = — i— U(k,o>)e ~ - doidk (6-13a)
(27r)3 J J -
v(x.t) =
(27T)
3
i(k-x + tot)
V(k,o))e - - dudk (6-135)
101
-------
In the case of the dice game we were able to define the entire family
of possible states A using only the definition of the system. To some extent
we can do that with the fluid system. First, we note that since u is a real
function, the transform U, which is complex, must satisfy certain conditions.
One of these is found by taking the complex conjugate of (6-9). We get
t/*(k,u>)
i(k-x + wt)
u(x,t)e - - dtdx (6-14)
V o
where U* denotes the complex conjugate. Comparing (6-14) and (6-9) we see
that
) = t/(-k, -co) (6-15a)
Also £/*(-k, -co) = U(k,co) (6-15b)
Multiplying (6-15a) by (6-15b) we find that
:)l (6-16a)
-k, -)| (6-16b)
In words, the amplitudes of the complex Fourier transforms must be even
functions of k and u in order for the fluid velocity components to be real
functions. Thus, the global family of transforms U is comprised of only
complex functions that satisfy Equations (6-15) and (6-16).
Strictly speaking this condition and any others that arise from the
realness of u are the only conditions that we should impose on the definition
of the global family. The constraints imposed by physical law should be
looked upon as properties of our causal wheel inasmuch as they are the re-
102
-------
suits of observations of the natural world. (This is the view we held
earlier in defining the global family of outcomes A of the dice game.)
Indeed, fluid flows that violate 1 or more of the so-called physical laws
may be possible but have such small probabilities that they have never been
observed. This philosophical argument notwithstanding, we shall proceed
to include physical law in the definition of the global family of velocity
transforms U.
2. The constraints of physical laws
Physical laws have the effect of imposing certain relationships among
the elements of each member of the U family. We will illustrate this here
considering only the mass continuity law in its incompressible form, and
the momentum law, both in the context of a 2-D fluid. Later we will consider
a more general case.
The continuity law is
Substituting (6-13a) and (6-13b) into this equation and performing the differ-
entiation, we get
kxu(kx,ky,w) + kyl/(kx,ky,u) = 0 . (6-18)
Thus, any function U = (U, I/) that does not satisfy both (6-15) and (6-18)
is not a member of the global family. All other U, of which there are in-
finitely many, are eligible members.
The momentum law in 2 dimension, and somewhat simplified, can be expressed
in the form
103
-------
where n is the vorticity defined by
- 3v 3u (z yn\
n- ' al'ay (6~2 }
and < represents molecular viscosity. Denoting the transform of n by
, we find from (6-19)
) + i W(kl,o»l)[kYt/(k - k',u - u1) (6-21a)
/ / ~ x ~ ~
+ k i/(k - k'.u - u'Jldk'du1 = -<(k* +
y «• x y
and from (6-20) we obtain
N(k,u>) = i[kYl/(k,u>) - kvU(k,a>)] (6-21b)
_ A - J -
Any complex function (1 = (U,V) that does not satisfy all of the relation-
ships, (6-21), (6-16) and (6-15) is not a member of the global family of
fluid velocity transforms. In atmospheric studies we deal with a 3-d i-
mensional fluid that is heated and cooled and that contains water. In this
instance to describe the state of the system we need in place of the vector
(u, v) used in the preceeding example a 6-dimensional vector (u, v, w, p,
e, q), where (u, v, w) are the fluid velocity components, p is the fluid
density, e its temperature, and q the mass of water vapor per mass of air.
If we let (1 represent this 6-D vector, or its equivalent, in Fourier space,
then in the manner illustrated above we could derive from the conservation
laws of momentum, mass and energy the relationships that define the global
family of a. Let's call this family n. The details of its definition,
corresponding to Equation (6-18) and (6-21), are not important at this point
and we can think of it simply as the intersection in function space of those
subsets of U each of whose members satisfies one of the physical laws.
104
-------
This is illustrated schematically in Figure 6-3.
The set ft represents all possible fluid motions. Specifying the fluid
(i.e., its viscosity, heat capacity, etc.) and suitable initial and boundary
conditions effects the selection from a of a single function that describes
the exact state of the fluid system in space-time under the conditions pre-
scribed. This is the essence of determinism, but it is generally unachiev-
able in practice for the same reason that the outcome of the throw of a
set of dice is not deterministic: the initial and boundary conditions can-
not be determined with adequate resolution.
When we specify the initial and boundary conditions of a fluid, we do
so only in some macro sense. That is, we say that a surface is "smooth"
but on the microscale it is very rough. Similarly, we say that a fluid is
initially "at rest", but on the microscale the fluid is in a constant state
of agitation (this is evident in Brownian motion). Generally speaking,
in low Reynolds number flows the state of the system is insensitive to micro-
scale details in the initial and boundary conditions. In this class of
problems a macro-specification of the initial and boundary conditions deter-
mines the macro state of the system uniquely. The analogous situation
with the dice is having a set that is so "loaded", that despite differences
in the way the dice are held or tossed, every throw results in the same
outcome matrix A. In effect, the causal wheel has only 1 system state on
it.
The same is not true at high Reynolds numbers. It is found, for example,
that as the pressure gradient driving fluid through a pipe is steadily in-
creased, a "critical" Reynolds number is eventually reached at which the
flow field suddenly changes from the organized, Hagen-Poiseuille state
105
-------
Global family of
Figure 6-3.
Illustration of the set n of complex functions
that
comprise the global family of fluid flow velocity Fourier
Transforms. Sets represents functions that satisfy purely
mathematical constraints like Eq. (6-15). Set A represents
functions that satisfy mass continuity relationships, like
Eq. (6-18); and set ir denotes functions that satisfy all
other laws like momentum and entropy. Intersection of A
and TT defines the global family fl of (1.
106
-------
characteristic of the low Reynolds number regime into a chaotic state known
as turbulence. No matter how many times the experiment is repeated, the
state of the fluid after the onset of turbulence is different and it differs
from one pipe to another. In these cases the flow is so sensitive to small
scale details in the initial state of the fluid and the K">undaries that
specification of only macro properties of the initial ar, ..^..ndary conditions
is only adequate to delineate the subset of fl to which the fluid state U
belongs, but not the exact state (1 itself (the causal v/heel now contains a
subset of n rather than a single member). Let's call this set W. We can
illustrate the differences among the members of W and something of the nature
of the set W itself with the following analogous problem in free convection.
Imagine a closed air chamber with thermally insulated walls and top
and with a bottom made of a high resistance metal that becomes hot when an
electrical current is made to flow through it. And imagine that a velocity
probe is placed at some arbitrary, fixed location inside the chamber.
Suppose now that after waiting several days for the air inside the
chamber to reach a "state of rest", an electrical current were applied to the
bottom of the chamber and simultaneously recording of the signal from the
velocity probe were begun. We would expect that after some brief interval
the velocity would depart from its initial value of zero and would vary in
an erratic manner from that time on. Let's suppose that in this first ex-
periment the velocity is recorded for 1 h and after that the heating current
is switched off.
Suppose that after waiting many days for the heat injected into the
chamber in the first experiment to escape and for the air to return once
again to a "state of rest" the same experiment were repeated. Our intuition
107
-------
tells us that the initial portion of the velocity record made this second
time might coincide with that made in the first experiment, but in time the
2 records would diverge and would differ in a "random" way for all later
times.
Since it is axiomatically assumed that fluids of physical laws
«
(whether those laws are represented accurately by the mathematical equations
currently used is another matter) and that the solutions of the equations
that express those laws are unique functions of the initial and boundary
conditions, the only way the velocity records collected in our hypothetical
experiments could differ would be for the initial, or less likely, the boundary
conditions to be different in the 2 cases. (Actually, to the author's
knowledge an experiment of this type has never been conducted, so the con-
clusion that the velocity would differ from one outcome to the next is only
speculation. We look upon this simply as a "thought experiment" that illu-
strates the nature of an ensemble of flow fields.)
Our statement that the air in the chamber is initially "at rest" is a
macro-specification of the velocity distribution that does not uniquely de-
termine the precise mathematical state v(r,tQ). In fact, there exists an
innumerable set of functions v(r,tQ), (r e P, where V is the space inside
the chamber) each of which is consistent with the criteria, implicit or
otherwise, that define the macrostate "fluid at rest". Let's call this set
of functions I. Then under our premise each member of I maps through the
equations that represent physical laws into a set W of functions v(r,t),
•w *»
r e P» ^0 1 * < •• Thus, W represents the subset of the global family fl of
velocity fields that is delineated by the conditions that define our chamber
experiment and it is the members of the set W that our velocity probe
108
-------
samples. Since there is no known mechanism that would tend to cause any
member of I to occur more frequently than any of the other members of the
set, then each velocity field in W is equally likely to occur and hence W
is also the ensemble of flow fields.
As in the dice game, there is no way of determining which member of I
characterizes the air in the chamber at the beginning of any experiment and
therefore there is no way of predicting the outcome v(r,t). We can only
hope to determine properties of the ensemble W.
If the angle brackets <> denote an ensemble average and rQ represents
the location of the velocity probe in the chamber, then and
are the mean and mean square velocities, averaged over infinitely
many experiments, that we could expect to observe.
Since the bounds on the subset W are determined by the parameters that
define our experiment, namely the dimensions of the chamber, the rate at
which heat is injected into it, physical properties of air, and the roughness
of the chamber walls, all ensemble mean quantities must be functions of these
variables. The aim of turbulence theory is to predict these relationships.
Two important points should be noted here. The first is that the "ensemble"
of flow fields is determined by the conditions that define the problem. In
the present example these include: thermally insulated, smooth- walled air
chamber with bottom heated uniformly at a specified rate; fluid initially at
rest; etc. In typical atmospheric boundary layer studies the defining con-
ditions include: flat infinite surface of given roughness length ZQ and
heat flux H ; friction velocity u* (which is a surrogate for the pressure
gradient that drives the flow); etc.
109
-------
The second point 1s that in conventional turbulence studies the condi-
tions that define the ensemble are usually "external" parameters such as
those just listed rather than "internal" parameters that we define below.
Perhaps the purpose for this is that once the relationships between ensemble
averaged state variables and the external parameters have been determined,
those expressions can be used to predict averaged values in other similar,
but quantitatively different, situations.
In regional scale air pollution studies we are generally not interested
in prognostic capabilities of this sort. We consider the flow fields in
the region of interest to be fixed, in a climatological sense, and we attempt
to determine how air quality in that region would change if the source emis-
sions were altered in a prescribed way. Furthermore, external parameters
are not readily definable in these problems because there are no walls or
tops on the airshed of interest, the earth's surface is not uniform, etc.
In these cases the conditions that define the flow fields are measurements
of "internal" parameters such as wind velocity, temperature, pressure, and
so on. Therefore, these data must define the ensemble of flow fields. And
as we show in the next section, they define in turn the corresponding en-
semble of concentration fields associated with a given source distribution.
To define the ensemble in terms of "internal" parameters, let's assume
that in the space-time domain (D, T) of interest we have available observa-
tions of the horizontal wind components (u, v), temperature e, etc.,
sites XL x2,... XN. Suppose further that these data are in the form of
-» •» *n
moving time averages.-
110
-------
t+T
Tjfn(t) = 2} J u(xn,t')dt' , n=l,...N (6-22a)
t-T~
t+T
yt) = ^l e(xn,t')dt' , n=l,...N (6-22fa)
t-T
etc. Let's assume also that in addition to the N surface monitoring stations
there are M upper air soundings made at locations XN+,, XM+2""XN+M ^at
give averages of winds, temperature and dewpoint over small vertical intervals
AZ:
{. ^ l^i-l C.
-Z 1 f
\(z>v = AT rv v *'•
Z + AZ/2
(6-23a)
2 " ^/2 m = N+l,...N+M
z + AZ/2
i
z " Az/2 m = N+l....N+M
etc. Since the functions on the lefthand side of (6-22) and (6-23) are given
(from the measurements), this set of equations constitutes a system of in-
tegral equations in the unknown flow fields u(x,t), e(x,t), and so on. We
now define the ensemble W of velocity fields associated with the space (P, T)
and observations un(t), en(t), etc in this space as the set of functions
u(x,t), e(x,t),etc (x,t) e (V, T) that satisfies jointly the system of Equa-
tions (6-22, 23) and the following set of equations that represents physical laws:
dV .
.£.. -IVp - f x V + Fu (6-24a)
|f = - pg (6-24b)
ff- = - V • (py) (6-24c)
111
-------
(6-24d)
where
da_qs9
at - —
^ j
5 H
= 1. If ;j±- <0 and q > q
""*"
= 0, if jt >p or_ q < qs
In Equation (6-24), V = (u.v); V = i •£ + j •£
(6-24e)
p is pressure; p is density; 0 is temperature; f is the Coriolis parameter;
Fj, is the friction term which is expressible approximately in terms of V
and other parameters and generally is significant only in the lowest 2 km
of the atmosphere; q is the specific humidity and q the saturation specific
humidity, which is a function of 9 and p; i is the concentration of liquid
water; c and c. are the specific heats of air at constant pressure and
water, respectively; R and Ry are the gas constants for air and water vapor,
respectively; g is gravity; 0 is the rate of heat loss or gain by radiation
processes or precipitation; and L is the latent heat of vaporization, or
sublimation, whichever"applies.
•
In our earlier discussion we defined the set n (see Figure 6-3) to be
the set of (u, v, w, p, e, q) that satisfy physical laws, namely the set
of equations (6-24). Thus, by imposing the additional constraints (6-22, 23)
we reduce the set of states (u, v, w, p, e, q) accessible to the system to
a subset of n, namely the subset we call W. By contrast, the classical
approach is to specify certain initial and boundary conditions in place of
112
-------
(6-22), (6-23) that in principle reduce the set of accessible states to a
single member of n. But as we discussed earlier, in the domain of n asso-
ciated with "turbulent" flows, specifications of the initial and boundary
conditions that provide only a macro-description rather than point details
are only adequate to delineate a subset of systems states.
The observations u , etc., that enter in Equations (6-22) and (6-23)
reduce the number of possible flow states from the set n to the subset W
for exactly the same reason that knowledge of the outcome after each roll of
1 of the 3 dice in the earlier game would reduce the number of possible matrices
A from the global set of 6 matrices to a subset of 6 . These reductions
effect an improvement in the certainty with which we can describe events
that have already taken place. Indeed, in the case where all the dice are
observed or where flow variables are measured everywhere in space-time, the
subset W of possible states is reduced to a single member of the global set
n. However, observations obviously do not exist for events that have not
yet taken place, and therefore in attempting to predict future states of the
fluid or outcomes A of the dice, the set of possible states is the ensemble
set which may be identical to the global set n. As we pointed out earlier,
the purpose of observations is to provide us with a knowledge of the ensemble
set and the probabilities of its members. We return to this point in Section 7.
In Part II of this report we will show how (6-22), (6-23) and (6-24) are
used to derive an ensemble of cell-averaged flow fields and cell layer heights
from which the regional model can generate the associated ensemble of pollu-
tant concentrations distributions, and we will begin to answer the 3 basic
questions raised at the beginning of this section.
113
-------
SECTION 7
MESOSCALE DIFFUSION; CONCENTRATION PREDICTABILITY;
MODEL "VALIDATION"; AND RELATED TOPICS
To derive the equations that describe the fate of pollutants in each
of the model's 3 layers, we started from the mass continuity equation
(2-1 j which involves the instantaneous, point fluid velocity v. The final
form (2-29) of the model equations governing the cell-averaged species con-
centrations in each layer contain the subgrid "eddy" flux terms
(7-1)
. = . - .
j j j
= j " j3
where <>. denotes a cell average in layer j, j =0,1, 2, 3; u and v denote the
components of the horizontal fluid velocity u; and primed quantities repre-
*•
sent deviations of the local, point value from the cell averaged value <>..
In the previous section we showed that from a discrete set of observations
of meteorological variables in a domain (V, T), one can determine the velocity
field v = (u, v, w) in that domain to within only a set of functions W. To
each member of this set there exists a . and . for each distri-
— — — j j — — -
bution of species sources. Thus, our first task in this section will be to
express these subgrid eddy fluxes in terms of * and .. Following that
— j j
we will address the 3 questions raised at the beginning of Section 6,
which involve the relationship between the set of associated with the
114
-------
set W of meteorological fields and the concentrations that one could expect
to measure under the meteorological conditions implicit in W. Finally, we
will look at whether the concepts we presented in the previous section are
consistent with the concepts of turbulent diffusion that form the basis of
the classical theories of microscale dispersion.
For simplicity we will approximate the subgrid eddy flux terms in the
gradient transfer form
-II*. \t v
j
_£
9. of the
" J
members of the set W is much broader at regional scale travel distances than
the spread of a plume due to subgrid scale velocity fluctuations in a single
member of W. One of the tasks that will be performed as part of the NEROS
field studies will be to assess the accuracy of Equations (7-2) and (7-3).
115
-------
Figure 7-1 illustrates that for a given distribution S(x,t) of sources
of J material species, that is, S(x,t) « [Si(x,t),...t Sj(x,t)], the prin-
ciple of mass conservation, which is represented in the present model by the
system of Equation (2-29), effects the unique mapping of each member of the
global family n of fluid states into a set r of functions in concentration
space. And each member of the subset W of n maps into a subset C of r.
If there is no process under the observed flow conditions that define W [see
(6-22), (6-23)] that would tend to favor the occurrence of any particular
members of W over others, then we must assume that each member of W is
equally probable and, therefore, that each member of C is equally probable.
Since there are infinitely many functions in the set W, the probability of
occurrence of any 1 of them is infinitesimally small. This fact precludes
our predicting the precise concentration distribution that will result in a
given situation and therefore, as in the case of the dice game, we must exa-
mine the members of C to determine whether they possess some property (like
the sum of the elements a^ -n of the matrices An that make up the ensemble
of dice game outcomes) whose value can be expected to lie close to some
predictable value, regardless of the member of C that actually occurs.
Consider the time-averaged value of each member of c(x,t) of C at a
given point x_:
-° t+T
c"(x,t)=4 c(x ,t')df , ceC (7-4)
t-f
where x e V, T <_ t <_ T - T, and T is an arbitrary time interval. (In order
to avoid awkward notation, we will use c in this section to represent the
members . of C).
116
-------
Flow field phase space
•-V
Law of mass conservation
(Eq. (2-29) with ,S (x, t) given)
Concentration phase space
Figure 7-1. Mapping of the set n and W in flow field function space
into the sets r and C in concentration function space for
a given source distribution S.
117
-------
By virtue of the assumption that the members of C are equally probable,
we obtain by averaging (7-4) over the set C the ensemble time averaqed con-
centration
=
_ 1
t+T
dt'
t-T*
(7-5)
Keep in mind that is the mean value of all possible time averaged con-
centrations c~ that could occur for a given source distribution S under
the discrete set of observed meteorological variables in "(6-22), (6-23)
that define W.
Let
e(x0,t) = c(xQ,t) -
(7-6)
denote the difference between the time averaged value c of any member of the
ensemble and the (predictable) ensemble mean value . Squarinq Equation
(7-6) we get t+T
- -~ J ]c(xn,t')c(xn,t")dt'dt"
41 t-T
t+T
(7-7)
c(x0,t)')dt'
t-T
Averaging this expression over the ensemble C we obtain
t+T
dt'dt"
- 2
(7-8)
Note that = 0.
118
-------
The value of provides a measure of how closely the time averaged
concentration Fat x can be expected to lie to the predictable ensemble
mean value .
From the Tchebycheff inequality of statistics, we find that with = 0
the fraction A of the ensemble population whose values c" are within ±eQ of
is
X i 1 - ±&- (7-9)
In words, in a large number of predictions the average frequency with which
the observed time average concentration c"(x ,t) exceeds the predicted value
by an amount larger than some fixed value eQ decreases as
decreases. Thus, our next task is to relate and ,
the quantities needed to compute and , to the known ensemble properties
of W.
Lamb (1975) has shown that for materials that are chemically inert or
undergo only first order reactions (i.e., the function R in (2-29) is a
linear function), the ensemble mean concentration is given by
ff*
= p(x,t|x',tl)S(xI,t')dt1dxl (7-10)
o
and
t f
=
I J J p2(x,t; x,f Ix'-.t'-.x"'^"1) (7-11)
0 0
S(xll,t")S(x"1,t"1)dt"dtllldx"dx"1
where S is the source strength function of the material species whose con-
centration is c.
119
-------
In Equations (7-10) and (7-11) the functions p and p2 are Lagrangian
properties of the ensemble W, and hence it is these functions that link the
properties of W and those of C for given S.
The function p is called the single particle displacement probability
density and is defined for chemically inert material by
p(x,t|x',tl) = <t|x1,t')>/6v (7-12)
Here, Sv denotes the (small) sample volume centered at x, <> denotes aver-
aging over W; and * is defined for each v(x,t) in W by
1 , if v(x,t) is such that a particle released
4>(x,t|x',t') =-
at (x',t') is in 6v centered at x at time t
0 , otherwise (7-13)
In the case of a point source of strength S(t) located at XQ,
S(x,t) = S(t)6(x - x0)
we obtain on combining (7-10) and (7-12)
1 ft
= (6V)"1 S(t')<4»(x,t|xn,t')>dt1 (7-14)
JQ
It is instructive at this point to comment on the function and the
physical significance of the forms it takes.
With the source location x$ and receptor site XQ given, we can plot 4>
as a function of t and t'. Figure 7-2 illustrates a hypothetical case.
Referring to the figure, we point out that a particle that leaves x at time
•W »
t'i enters the sample volume 5v at x at time ti and exits at time ti + At.
The value of the time integral of , which plays a role in (7-14), is also
indicated.
120
-------
t=t'
t' t'+I
Figure 7-2.
Plot of a hypothetical example of the function
,t|x ,t')
for given source location x$ and receptor site x . The
function 4> has unit value inside the shaded regions and is
zero everywhere else.
121
-------
If x is one of the points at which the meteorological data that enter
in the definition of W are collected, we can expect the functions $ corre-
sponding to each member W to have the same general features. By contrast,
if x and XQ are far from a meteorological station, there will generally be
large variations in the form of $ from one member of * " ensemble to another.
(The reasons for this and its ramifications on model p* -cncbions are dis-
cussed under Remarks later in this section.)
If the plots of all 4> within W are superposed and averaged at each point
(t,t')t we obtain the function <4»(x0,t|xs»t')> needed in (7-14) to compute
. If the distribution of <> is aligned parallel to the line t = t1
and if its .width and amplitude are constant along its length, then 41 is a
statistically stationary function. In this case it is easy to see that
f*
< $(t,t')dt'> = (7-15)
where <&t> is the average time spent by a particle released at any time t1
in fiv (cf. Figure 7-2). If, under these conditions, the source strength is
also steady in time, say S(t) » SQ, we have the simple result
= s7 (7-16)
In essence, it is this simple expression that the Gaussian plume formula
approximates in an empirical way.
The assumption of stationarity is applied almost routinely in turbulent
diffusion studies partly because it simplifies the mathematics. But it is
not hard to see by visualizing the shapes of (t,t') that various flows would
produce that stationarity is generally not applicable in regional modeling
studies. In these cases calculation of must be based on (7-13) and
122
-------
(7-10) rather than equations like (7-16) that assume stationarity.
The function p2 that enters in (7-11) is called the two-particle dis-
placement probability density and is defined as follows:
P2(x,t;x,t'|x",t";x"1,t"1) = <*2(x,t;x,t' Ix''^'"^" ,t"' )>/(Sv)2 (7-17a)
Here, as in the definition of p (see 7-12), the angle brackets <> denote
averaging over W; 5v is the sample volume centered at x; and 2 is defined
for any v(x,t) in W by
" 1 , if v(x,t) is such that a
particle released at (x",t")
and one released at (x"',t"')
are found at (x,t) and (x,t'),
respectively:
(7-17b)
0, otherwise
»
For a single point source the autocorrelation is derivable
from the contracted form 2(x,t;x,t' |xs,t";xs,t"'). This function contains
4 independent variables for given x and xs and is much more difficult to
visualize than $.
To summarize briefly the points made thus far, we have argued in the
previous section that a discrete set of meteorological observations, as given
in (6-22), (6-23), in a given space-time domain (V, T) is adequate to specify
the flow field v(x,t) to within only a set of functions W. For a given source
distribution S, each member of the set W maps through the species mass con-
servation equation [e.g., (2-29)], onto a function element of a subset C in
concentration space. Since it is impossible to say which member of W actually
describes the state of the fluid in (V, T) under the discrete set of con-
ditions specified, it is impossible to predict the precise concentration dis-
tribution that would occur under those conditions for aiven S. As in the
123
-------
dice analogy, we can only hope that there exists some property, say g(c),
such that g(c) - for nearly all c e C. Here is the mean
value of g(c) over the set C (or the ensemble mean), and it is a quantity
that, in principle, can be calculated for any g(c). We showed that if g(c)
is the moving time average defined by (7-4), then a fraction A of the mem-
bers of C possess g(c) = c~ values that are within ±e of the set mean value
, where e is an arbitrary interval and \ is given by (7-9).
At this point we have answered the first question raised at the end of
the introduction to Section 6 and we have obtained partial answers to the
other 2 questions. Specifically, we have shown that for a given source
distribution S and a given set of meteorological observations, a model can
delineate the set of functions C to which the actual concentration distri-
bution c(x,t) resulting from those conditions belongs, but not the function
c itself. From the model generated function set C, one can compute the set,
or ensemble, mean values of c such as the first and second moments given
by (7-10) and (7-11) and from these one can specify the frequency with which
actual concentration can be expected to fall within a given interval of values
in a large number of observations. The width of this interval [see (7-9)] is
a measure of the "predictability" of the concentration and it is determined
by the properties of the set C. These in turn are determined by S and the
set of meteorological data that define W. Thus, "validating" a model is
tantamount to testing whether observed concentrations fall within the inter-
val specified with the frequency specified and, if not, whether the failure
can be attributed to sampling fluctuation. Note that from the standpoint of
124
-------
regulatory needs the utility of a model is measured partly by the width
of the interval in which a majority of observations can be expected to fall.
If the width of the interval is very large, the model may provide no more
information than one could gather simply by guessing the expected concen-
tration. An important task in this regard is the "inverse" problem of
specifying the density and type of meteorological data that are needed to
achieve a given level of concentration predictability. We will consider
these topics again later in this section.
The conventional approach to turbulent diffusion modeling is to formu-
late equations for and possibly higher moments that yield these quan-
tities directly. By contrast, the approach we are following here is to
generate members of the ensemble set C and then to compute the set mean
values , etc., by averaging the individual members. In the author's view,
this is the only approach that can produce accurate estimates of ensemble
mean values, because the closure problem that must be overcome to compute ,
, etc., directly is made virtually insurmountable by the combination of
the nonlinearity of the chemistry and the large time and space scales of
atmospheric motion. In Part II of this report we will describe our progress in
implementing this approach.
Remarks
1. The character of diffusion at long range.
In this subsection we would like to examine the classical definition of
turbulence employed in microscale dispersion analyses from the perspective of
the mathematical set definition we are proposing here. Our main aims are
to examine the cause of the paradox described in the introduction to Section
125
-------
6 that results when the classical definition is applied to long range diffu-
sion and to determine whether this paradox is reconciled by our definition
of "turbulence".
Consider the frequent situation in microscale studies where only a
single wind observation site exists, at XQ, and a single point source of
material is located nearby at xs. Under our definition, this observation
establishes through (6-22), with N=l, and (6-24) a set of functions W that
contains the field v(x,t) that describes the actual fluid velocity in the
domain (D, T) in which the observations were made. Suppose that the ob-
servation at x is a moving time average, as defined by (7-4).
If we were to superpose the trajectories originating at xs (= XQ) de-
derived from each flow field in W, we might obtain a pattern like that illu-
strated in Figure 7-3. The spread of these trajectories is what we regard
as "turbulent diffusion".
Since each member of W satisfies (6-22), all trajectories have the same
general direction, namely that of ZT(x ,t), near the source; but farther down-
stream there is increasing freedom in the direction that the flow there, and
hence the trajectories, can take and still satisfy (6-22). For example, on
occasions when the wind vector is 17 at x , the flow can turn cyclonically
downstream while on other occasions, with an identical flow at x , it may
turn anticyclonically. Without more information than simply the wind vector
at XQ, there is no way to say which direction the trajectories will take.
Indeed, it appears obvious from this perspective that the spread of trajec-
tories far from the source is affected by large scales of motion that one
normally associates with transport.
126
-------
igure 7-3. Trajectories of particles eminating from xs in the ensemble
W~ associated with the single wind observation u(x ,t) = u ,
where x£ - XQ . - "
127
-------
We can verify this important point by analyzing the problem mathemati-
cally. We will show that with "turbulence" defined in the classical way,
its spectrum is an implicit function of distance from the site at which
the observation used to define the "mean wind" is made. In the present
example we are assuming that a wind observation is made at the single location
x and it provides the moving time averaged velocity
rt+T
TT0(t) ^j u(x0,t')dt' (7-18)
t-T~
where u(x,t) represents the instantaneous velocity at (x,t). Expressing u
in terms of its complex Fourier transform defined by (6-8) we have
u(Xft) = —— k,u,)e udkdu (7-19)
(2*)* J J- -
where U is the complex transform at wave number k and frequency u. Combining
(7-18) and (7-19), we get
u0(t) - u(k.«) e ' - dkdo, (7-20)
For notational convenience we will assume that the factor (2ir) is absorbed
in U.
Since u is a real function, we must have
f f -i(k-x + ait)
u(x.t) = U*(k,«a)e ' - dkdu (7-21)
- - 11-- -
where U* is the complex conjugate of U. Thus the product of u measured at
2 points separated in space-time by (6, T) is
128
-------
u(x,t)u(x + 6, t + T)
lu(k,u)U*(k',a)1)
(7-22)
i(k-x + ut) - i[(k'Hx + 6) + a>'(t + T)]
Averaging this product over the ensemble W defined by the observation at x£
^
[and (6-22) and (6-24)], we have
=
- ik'-6 + i(u>-u>')t-iVT
dkdk'dujda)'
By definition if the set W of velocity fields is homogeneous and sta-
tionary, then.
= R(6,-r) (7-24)
That is, the set average of the velocity product depends on the separation
(S,T) of the velocity observations but not on the location (x,t) in (fl, T)
at which the pair of observations are made. Comparing (7-23) and (7-24) we
see that stationarity and homogeneity can exist only if
=
(7-25)
where 6( ) is the Dirac delta function. Combining (7-23, 24, and 25) we get
-i(k-6
(7-26)
Now under the classical definition, the turbulent velocity component at
an arbitrary point (x,t) is
u'(x.t) = u(x.t) - u0(t)
(7-27)
Thus,
u'(x,t)u'(x + 6, t+t) = Cu(x.t) - u0(t)][u(x + 6, t+r)
- U0(t+r)3
(7-28)
129
-------
Taking the ensemble average of this expression and making use of (7-25) we
obtain
=
sinojT cik'(x - XQ)
- - -1 IDT
- XQ)
(7-29)
, /
' oJT "^
Evaluating this expression at 5 = T = 0 we obtain the spectral representation
of the turbulent energy at (x,t):
- XQ)
(7-30)
This expression shows clearly that while the fluid velocity u is stationary
and homogeneous, "turbulence" defined by (7-27) is only stationary — its
energy is a function of distance x - x from the wind monitoring site.
It is easy to see in (7-30) that as the distance from x increases, pro-
gressively larger scales of the velocity u become embodied in u1.
If we evaluate (7-30) at the wind measurement site x , we obtain the
conventional description
*(k.«)[l -
(7-31)
which conceals the inhomogeneous character of u1
130
-------
It is the inhomoqeneity of u1 that is responsible for the dispersion
paradox described in Section 6. To show this we evaluate (7-29) for 5 = 0
and use the result in the definition of the Eulerian integral time scale
V
1
J
dT (7-32)
0
Upon making use of the relationship
1
2:r
e"1a)tdt = 6(u) (7-33)
we obtain
-r / \ _ 2ir
(7-34)
This expression shows that at the point x where the properties of turbulence
are measured, the integral time scale T£ is identically zero. And it is from
this result and the assumption that the Lagrangian time scale T, is related
to TE by
TL = CONSTANT • T£ (7-35)
that the aforementioned dispersion paradox arises.
Eq. (7-34) reveals that the observed spreading of a plume at long range
is due to increasingly larger scales of motion, scales that are considered
part of the "transport" or "mean" wind at x . Consequently, if one infers
values for T. from atmospheric data and then uses these values in a regional
model in which wind observations are made at multiple sites, the result will
be that certain scales of motion are treated twice, once in the transport
flow and once (implicitly) in the turbulent diffusivity that is based on TL>
We consider this problem next.
131
-------
2. Parameterization of dispersion
If we increase the number of flow observation sites in (V, T) from 1,
as in the example above, to N, the set of equations (6-22), (6-23) and (6-24)
that defines the set W of flow fields acquires a different form and hence a
new set W of flows applies. Let Wj and W2 represent the sets of flow fields
corresponding to a single flow observation, as in Figure 7-3, and to 2 ob-
servations, respectively. Suppose that in the case of 2 measurement
sites, one is located at the same site x as the monitor in the single
observation set, and that the other is located at xi downwind of x . Then
the sets of trajectories, associated with the sets W, and W2 will differ.
For example, those trajectories in W, that pass through xi but which have
a speed and/or direction at xi that is inconsistent with the observation at
xi are not members of W2- Thus, the addition of a second flow observation
has effectively altered the character of dispersion. This example illustrates
that dispersion is an artifact of our observations, rather than an intrinsic
property of the atmosphere, and it must be treated in this manner when
attempting to approximate it.
The proper way to estimate the functional form of an effective eddy
diffusivity for use in a diffusion equation model is the following.
First, calculate for the ensemble of concentration distributions
associated with the set W of flow fields. Second, define the transport
wind field to be the ensemble mean wind in W. Finally, using the
known and solve the diffusion equation "in reverse" (i.e., the
inverse problem for the diffusivity K(x,t)). Lamb and Durran (1978) illus-
trate this process.
132
-------
In problems involving nonlinear chemistry, the eddy diffusivity derived
by th'is procedure would be an implicit function of the source distribution
S used to generate the concentration ensemble C, and therefore the model
would not be applicable to problems with different source distributions.
For this reason, we propose to generate the ensemble properties of concen-
tration like from a set of functions approximating C rather than from
an ad hoc model of the ensemble mean .
3. Predictability of concentration at long range
We have emphasized in this section that for a given set of meteorological
conditions and a given source distribution S, one can determine the concen-
tration distribution c(x,t) that would result to within only a set of func-
tions C(ceC). Therefore, in order to utilize a model in decision making
processes, it is necessary to know how much variation exists among the mem-
bers of C. For if the variation is large and each member of C is equally
likely to occur, then the information provided by the model has little or no
practical value.
In Equation (7-6) we defined a parameter E that provides, through Equation
(7-8) and (7-9), an estimate of how closely the time-averaged members F of C
are clustered around the set mean value . Using this parameter, we will
illustrate briefly below that the scatter of c~(x,t) about increases
as the distance x - x from the source increases (assuming x = xrt, the
•v ~5 *,S 4.0
wind observation site).
Referring to Figure 7-4, let's assume that the 2 plumes shown represent
the extreme lateral deviation from the x axis of all plumes in the C ensemble.
By "plume", we mean the superposition of particle trajectories originating
133
-------
Fi gure 7-4.
Plumes from the C ensemble that mark the limits of
lateral motion. .
134
-------
continually at the source over some small time interval for each u(x,t) in
W. Thus, the plume envelopes in the figure represent the locus of travel
of the particles released during this small time period. Also, the plumes
shown assume that C is an ensemble in which a wind observation site is lo-
cated at the source but that no other observations are available.
For simplicity, we will assume that the concentration ,. .orrelation
has the form
= •< ~ (7-36)
where T is arbitrary and independent of t and x. Assumino that within C,
plumes are equally distributed over the lateral ranges L, and 1_2 shown in
Figure 7-4, and their width I is nearly the same within C for a fixed x, we
obtain by straightforward calculations based on the definition of ensemble
averaging
= 2
2 (7-37)
= 2 L2/£2.
Now from (7-36) and (7-8) (and assuming stationarity) we obtain
= Y (2 - y)( - 2) (7-38)
where, it will be recalled, e is the difference between an observed concen-
tration averaged over the moving interval T, or equivalently a member c~ of
C, and the corresponding predicted ensemble mean value . Making use of
(7-37) and assuming t\ - £2 = £, we aet
2 2
L2 - I
-------
Thus, for the same averaging time T, concentrations measured at x2 can be
expected to differ from the mean value at that site by a larger fractional
margin than at the point xi closer to the source. Judging from the Prairie
Grass data, on which the empirical Gaussian plume formula is based, we ex-
pect that for points xi within 1 km of a source, the scoter in 15 min aver-
t
aged concentrations c is relatively small if meteorolog -"i conditions are
steady over comparable time periods. However, according to (7-39) the
scatter in similar measurements made at a site X2 farther from the source
will be greater by an amount proportional to the square root of the ratio of
the plume envelope widths L measured at sites X! and x2. This result is con-
sistent with our earlier finding, Equation (7-34), that the local integral
time scale of turbulence increases with distance from the flow observation
site xfl.
If the flow were stationary, the scatter in c~ measurements at x2 would
be smaller for concentrations averaged over longer time intervals than 15
min. But atmospheric flows are not stationary over arbitrarily long periods.
The alternative to reducing the spread of c" about is to restrict the
spread L within the flow ensemble W. This would require a redefinition of
W using additional meteorological observations at sites in the vicinity of
the source. For example, it is easy to see that with an additional wind
observation near x2» the spread of trajectories in the new W ensemble would
be less in general than in the single station ensemble shown in Figure 7-4
because some trajectories in the latter would be eliminated by the constraint
of the observation at x2.
4. Model validation exercises
Our discussion in the preceding subsection illustrates the natural
deviation between observations and model predictions that is not attribu-
136
-------
table to model error. The proper way to establish whether the model assump-
tions and data are correct is to examine the frequency with which observed
values satisfy (7-9). If the observed frequencies differ from that pre-
dicted by (7-9) by an amount that is not likely to be due to sampling fluct-
uations, the modeling hypotheses must be judged to be in error. Obviously,
we can find an interval 2e centered at within which all observations
will lie. But as we pointed out earlier, if the size of the interval in which
a majority of the observations can be expected to fall exceeds some value,
the model is worthless.
5. Predicting impact of future emissions distributions
Currently the method of judging the effectiveness of proposed emissions
control strategies is to predict concentration levels that proposed sources
would produce under so-called "worst case" meteorological conditions. Within
the context of the approach we have developed here, this entails defining
the ensemble W for the set of "worst case" meteorological conditions and then
computing and for the proposed source distribution S(x,t). However,
since meteorology and the source distribution jointly affect concentrations,
the concept of "worst case" meteorology is not well defined. Furthermore,
the flow regime associated with a past pollution episode might have such a
low probability that its expected frequency of occurrence is, say, once in
a century. A more meaningful approach, therefore, would be to compute
and for the ensemble C associated with the flow ensemble W where W is
the union of the sets W^, i=l,2,...M, each of which is defined by Equations
(6-22), (6-23) and (6-24) from meteorological observations in (V, T.),
1=1,2,...N; and the T.. are a collection of time intervals selected from the
period in which climatological records are available for the modeling region
137
-------
V. By this definition W is an estimate of the flow ensemble set, referred
to at the end of Section 6» which forms the basis for predicting future
states of a system.
138
-------
SECTION 8
FORMULATION OF MISCELLANEOUS FIELDS AND MODEL PARAMETERS
Model parameters not yet considered will be discussed in this section
and possible techniques for deriving them will be described. Chemical
kinetics schemes will not be treated in this report because no particular
scheme i-s an integral part of this model. Rather, the chemical mechanism
is a peripheral component that can be interchanged as desired. A number of
kinetics schemes are currently used in various urban modeling studies and
in the course of our regional modeling work we plan to utilize several of
them in comparative studies. At the time of this writing we are performing
operational tests of the regional model computer code and for these purposes
we have implemented the Demerjfan-Schere (1979) mechanism that simulates 36
reactions among 23 chemical species.
Deposition Velocity, 6.
Following Wesley and Hicks (1977) we shall express the deposition velo-
city e of each species in terms of the deposition resistance rg:
,t;n) = 0.4u* [Zn(10/z) + 2.6 + 0.4u*r - J"1 (8-1)
where the surface roughness ZQ and friction velocity u* are functions of
(x,,t) and the resistance r is a function of space and time and the pollu
tant n. The function t|»c in (8-1) is an empirical function given by
* = exp[.598 + .39 JM (-10/L) - .09 [In (-10/L)]2] (8-2)
139
-------
and L is the Obukhov length at (A,4>,t).
The surface roughness z (A,) is prepared from land use data. The
resistances r are functions of both land use and time of day. The methods
we used to formulate both these variables are described in Part II.
Source Emission Rates i, S
The stationary source emissions inventory will be prepared in such a
manner that the information necessary to make plume rise calculations will
be available for the set of sources that together are responsible for 50
percent of the total emissions of a particular pollutant in the modeling
region. Given the stack parameters and a plume rise estimate for each hour
of the simulation, the emissions of this set of sources will be partitioned
between i and S, depending on whether the estimated effective plume height
is above or below H .
The emissions of all remaining stationary sources and those from mobile
sources will be included in S. It is important to note that i is a
volume source, that is, it has units of mass per volume per time, and it
represents uniform emission rate of pollutant from each point within a
single grid cell. By contrast, S is an area source — units of mass per
area per time. Emissions from the 10 or so largest point sources in the
modeling domain (all power plants) will be handled separately and treated
by a puff type model embedded within the regional scale model. Special
treatment of these sources is necessary because the injection of material
into the upper layers of the model in the form of slender, highly concentrated
plumes can cause enormous subgrid scale chemistry effects at altitudes where
there is currently no provision in the model to treat them. Moreover,
140
-------
these sources can cause large subgrid scale variations in ground-level con-
centrations that the grid model alone cannot resolve.
The Plume Volume Fraction, c
This parameter was introduced in the Layer 0 formulation as part of
the scheme to treat subgrid scale chemistry associated with low level
point and line sources. It was defined in (5-11) as the volume fraction
of any Layer 0 grid cell occuppied by plumes from point and line sources
within that cell.
As in (5-6) let AZ represent the depth of a Layer 0 cell whose hori-
zontal area is A. In the surface layer of the atmosphere we can assume to
good approximation that the relative rate of expansion of a point source
plume is proportional to the friction velocity u* and that the centroid of
the plume rises at about the same rate. Hence, the centroid of a surface
source plume reaches the top of Layer 0 in a time
T * AZQ/U* (8-3)
and it has a diameter a there of order u*x = AZ . The volume 6v of a single
point source plume in Layer 0 is thus
6v * (AzQ)3 (8-4)
If there are p point sources per unit area, then
e= sir <• ><4zo>2 <8-5'
0
This expression is not valid if many of the point source plumes overlap.
In fact, since c <_ 1 by definition, we see that the total number N of point
sources in any grid cell must satisfy
141
-------
in order for (8-5) to be a consistent approximation. In our case where
A * 400 km2 and AZQ < 100m, Equation (8-5) is consistent as long as there are
many fewer than about 4.10*1 point sources within any grid cell.
Suppose there is a total length L of nonoverlapping line source plumes
in a given cell. The volume fraction of these plumes, derived by the same
analyses used for the point sources, is
LAz
"A
C * -IT2- (8-6)
When both point and line sources are present we shall assume that
5 = p(AZ0)2 t-j-i (8-7)
where p and L are values derived from the source inventory. A much more
detailed method of approximating 5 is given in Lamb (1976).
Turbulence Parameters w_, A_, etc., on H and in Layer 0.
In Equation (5-2a) we introduced the probability density p (w) of ver-
ro
tical velocity relative to surface HQ. In this section we shall derive values
for the parameters x_, X+, w_ and w+ that stem from pp [see Equations (5-2)
o
and (5-3)] by assuming that the turbulent velocity fluctuations at the level
of surface HQ have a Gaussian density, viz
w 2
P(wQ) = ^—exp (--?-) (8-8)
0 0
Implicit in this expression is the earlier assumption that the mean vertical
motion w"D at the level of H is everywhere and at all times zero. Since
o
142
-------
the velocity of surface HQ is 2Qj it follows from (5-1) and (8-8) that the
density pp (w) is
° i (V w>2
Pr M = — ^- exp [- -2 - ] (8-9)
Thus, from (5-2a) and (8-9) we get
A+ = h\.l - erf (—2—)] (8-10)
where erf ( ) is defined by (4-43); and from (5-2b)
A_ = 1 - A+ (8-11)
Using (5-3a) and (8-8) and making use of the integral equivalence (4-42h),
we obtain
a
W" I *° * 2o r. - , *o ,
\ ~ y T* «^~ L i ^ sri ^ — ^
2°»o ^o
and from (5-3a)
^n *n ^n *„
w = £ exp (- -2-) - / [1 - erf (-^— )] (8-13)
The speed 2 of surface H is obtained directly from the definition
of HQ [see (2-2) and (3-12)]. The variance aw of turbulent velocity fluc-
o
tuations is estimated from the following expression proposed by Binkowski
(1979):
o 'm
where
e-^
5 " L
and u* is the friction velocity. All of the terms appearing in (8-14) will
143
-------
be described in detail in Part II.
The Plume Entrainment Velocity, v
We use as a measure of the plume entrainment rate the rms spread rate
of particle pairs in the surface layer. Using similarity theory, Batchelor
(1950) showed that in homogeneous turbulence the rate of expansion of a
cloud is
I*- cie2/3 -£02/3T2 (8-15)
where
c, = a constant of order unity,
e = the energy dissipation rate of the turbulence,
T = the travel time of the cloud,
1Q - the initial cloud size,
Z2" = the mean square diameter after a time T.
Deardorff and Peskin (1970) have shown that Equation (8-15) is also
valid for shear flows, such as that characteristic of the lowest layer in
our study. Defining
and making use of Equation (8-15), we obtain
v = Cl1/2(E£0)1/3 (8-17)
Combining this expression with the surface layer formulation of e given by
Wyngaard and Cote (1971), namely
,,3 7 ?n 3/2
c * g{l + 0.5|£|*/d) (8-18)
where k is the von Karman constant (= 0.4) and L is the Obukhov length, we
get
144
-------
(zk)
1/3
(8-19)
where c. is a constant of order unity. If we take the release height z of
the cloud to be z = t (which is not unreasonable for motor vehicle emissions)
and £ to be approximately 1 m, then the term in bracket' 'i Equation (8-19)
has a value near unity for values of L typical of urban i.gions, consequently,
v = au* (8-20)
where a is a constant that we assume for the time being to be 1.
Rainout and Washout Processes
The term W on the righthand side of the general model equation (2-1)
represents the scavenging rate of pollutant species c by precipitation pro-
cesses. We can express it in the form
W = Ac (8-21)
where c is the material concentration and A is the scavenging coefficient.
The latter is a function of many variables among which are rainfall rate,
gas solubility, and in the case of aerosols, size distribution. McMahon
and Denison (1979) have summarized the few measurements of A available to
date and attempts to relate A theoretically to some of the parameters cited
above. We win postpone further consideration of wet removal processes
until a later report; because in the near-term oxidant studies that are
of concern to use, the meteorological regimes of interest are generally not
associated with precipitation.
145
-------
SECTION 9
COMPUTER SOLUTION OF THE GOVERNING EQUATIONS
Introduction
The differential equations that we formulated in the previous sections
to describe concentrations in each of the model's 3 layers cannot be
solved in general without the use of a large computer. If a digital machine
is employed, each of the concentration and input parameter fields must be
represented in discrete value form. In this case it is also customary to
model the governing differential equations by a system of finite difference
equations. There are countless methods of generating finite difference
analogues of differential equations but many of them do not preserve essen-
tial properties of the differential equations that are important in simu-
lating physical phenomena.
For example, in modeling regional scale air pollution it is not possible
in practice to represent the concentration fields in discrete form with
sufficient spatial resolution to resolve all of the major variations in the
concentration distribution. In some finite difference representations of
the advection equation, this limited resolution gives rise to the well
known phenomenon of pseudo-diffusion, which is not present in the true
advection process.
146
-------
Another source of potential trouble in applying a conventional method
to our set of equations is the large vertical concentration gradients that
can arise. In the atmosphere, there are extended periods (usually the night-
time hours) when vertical fluxes of material are very small. During these
periods very large vertical gradients in concentration can arise, and these
lead in a simulation model to large surges of material from one layer to
another once the physical processes that produce interlayer exchanges begin.
With many finite difference expressions, large transients of this kind can
cause computational instability.
Finally, the multitude of chemical reactions that make up the photo-
chemical air pollution phenomenon possess an extremely wide range of char-
acteristic time scales. A common method of mitigating the problems that
this disparity of scales causes the computer simulation is to adopt the
steady-state approximation for the fastest of the reactions. The effect
of this approach is to "hardwire" a particular chemical kinetics parameteri-
zation scheme into the model. Consequently, it is virtually impossible
to examine alternative representations of the chemistry and hence the
utility of the model as a research tool is diminished.
In this section we seek to develop ways of generating approximate
solutions of the governing equations that avoid the problems just cited.
We hope to achieve this in 2~ways. First, we will express the governing
equations in an alternate form that permits the advection, vertical flux
and chemistry phenomena to be treated separately and in modular forms.
This will allow the computer version of the model to have sufficient flexi-
147
-------
bility that various numerical schemes can be utilized without having to
overhaul the entire code. Second, we will develop numerical techniques
for treating advection, vertical fluxes and chemistry that are well suited
to handling the specific problems that arise in regional scale modeling.
These schemes are developed later in this section following a modification
of the governing equations.
Transformation of the Governing Equations
The differential equations that describe concentrations within each
of the model's 3 layers that the general form
(see Section 10) where cn = n< Here L is the linear advection and diffu-
sion operator (which for illustration purposes only we write in the K- theory
form)
1 5 u 3x~ +
Sn is the source emissions in layer n; Rn represents all chemical reactions
in this layer; and Fn is a function involving concentrations in layers
m f n that describes material fluxes from one layer to another.
Our Intent here is to restructure (9-1) so that the advection, vertical
flux and chemistry processes are decoupled and treatable separately. This
suggests looking for a product solution
-------
Substituting this into (9-1), dropping the subscript temporarily for
notational convenience, collecting terms in y and r, and assuming that spatial
variations in KH are of a much larger scale than those in r|, we obtain
Suppose that we set the first group of terms in brackets equal to (-F + R)/r
and the second group of bracketed terms equal to S/Y. Then if r is the
solution of the equation
with initial conditions
r(r. t0) = c(r, t0) (9-6)
and boundary conditions identical to those on Eq. (9-1); and if y is the
solution of
with initial conditions
r(r. t0) - 1 (9-8)
we see that c = yr satisfies Eq. (9-1) and its initial and boundary con-
ditions for a period t < t < t + AT. The magnitude of dT is determined
by the rate at which spatial variations in y arise. Note that r has units
of concentration while y is dimensionless. Note also that by construction,
149
-------
Y is initially uniform in space. However, according to (9-7) horizontal
variations in Y are generated by inhomogeneities in the vertical fluxes F
of material, in the rate R of chemical processes, and in the source strengths
S. Once horizontal gradients in Y have been generated, the last two terms
on the left side of (9-4), which effectively couple Y *^ r. became finite
and c = Y? is no longer an exact solution of (9-1). - r are also
coupled through the chemistry term R, which we treat later.)
Working within this limitation on AT, we can use r, as given by (9-5)
and (9-6), and Y» given by (9-7) and (9-8), to obtain approximate solutions
of the differential equation (9-1) at discrete intervals AT forward in time.
The procedure is as follows. First, we use the initial concentration dis-
tribution c(r, tQ) in (9-6) to solve (9-5) for r at time t = tfl + AT. Note
that r is independent of emissions, chemistry and vertical material fluxes.
Concurrently, we solve (9-7) and (9-8) for Y at time t = tQ + AT. We then
form the approximate solution
c(r. t0 + AT) = Y(t0 + AT)r(r, tQ + AT) (9-9)
We now use this expression in (9-6) to solve for r at time tQ + 2AT. The
equations governing Y» (9-7) and (9-8), have the same form during this second
time interval except r(r, tQ + AT) is used for r in (9-7). This process can
be continued indefinitely in time provided that the size of the interval AT
is kept within the limits discussed above.
At this point we have succeeded in splitting the basic model equation
into 2 parts, one of which is independent of the emissions, chemistry and
vertical fluxes. This component is given by Equations (9-5) and (9-6) and
it can be treated by any suitable 2-dimensional numerical scheme without
regard to the complicating effects of chemistry.
150
-------
The second part of the solution y is affected by emissions, vertical
fluxes and chemistry and we wish now to split this component into 2 subparts
which account for these processes separately.
As before, we shall try a product solution of (9-7) and (9-8); but
first it is convenient to exploit the smallness of the horizontal gradients
in Y during each interval AT and write (9-7) in .approximate form
f£+ ay + F/r = R/r+ s/r (9-10)
where y now pertains to arbitrary, local regions of space in which F, R, and
F have prescribed values. A better approximation is obtained by dropping
only the second order derivatives and translating to a system moving with
the fluid. In this case 3y/3t in (9-10) becomes dy/dt.
It is necessary now to resume the use of subscripts to denote the
model layer. We then have from (9-10)
tr+ Vn + Kn-l.n ' W ' Vn * Vrn <»
where ^ m is the flux across surface Hn to or from Layer m, and h is the
local thickness of layer n. The flux Fn k can be written in the general form
F „ = I [b/ c ] + g'
n,k m=l n
Thus, let
< Fn-l,n - Fn,n)/hn ' Sn + ancn = bnlcl + bn2c2 + bn3c3
Combining this and (9-11) we get
i\
+ bnl + bn2Y2 + b3 = 9n/Fn + Rn/rn
151
-------
where
b* = b r /r (9-15)
nm nm nr n v '
Note that the term a , which appears in (9-11), is absorbed into the coeffi-
cient bnn.
We look now for a product solution of (9-14) of the form y = yny'
in which the component y^ represents the chemical processes and yn describes
the effects of vertical material fluxes and sources on yn- Substituting the
proposed solution into (9-14) and collecting terms we have
A
nyn
A
The chemistry-free component y is therefore the solution of
A
Y«
with initial condition
Yn(t0) - 1 (9-18)
and
C ' bnm
Similarly, the component y^ that handles the effects of chemistry satisfies
v
nYn
152
-------
with initial condition
TjfV - 1 (9-21)
We must solve (9-20) and (9-21) for each of the pollutant species present
in each layer. For species a, we have
r(a)nY(a)n
(9-22)
with Y'( \n = 1 for all a and n. For bimolecular reactions of the form
prevalent in photochemical air pollution, R/a)n can be written
I I
^mcm",, (9-23)
where I is the total number of species present and where < >n denotes a cell
average in Layer n, as used throughout the analyses in this report [(see
(2-3)]. We discussed in Section 5 the existence of subgrid scale variations
in the concentration fields and pointed out that they are probably most
pronounced near the ground where they are generated by the highly inhomo-
geneous field of sources. Based on this assumption and motivated by the
need to keep the mathematical structure of the model tractable, we developed
a scheme for parameterizing the subgrid scale concentration variation effects
in the lowest layer of the model, Layer 0, but we chose for simplicity to
neglect subgrid scale phenomena in the other layers (n = 1,2,3). In keeping
with this assumption we will assume that the product averaged term in (9-23)
can be approximated by
,
Combining this with the product solution form c = Tyy1 we obtain finally
153
-------
« I T T
(<*)" , r r |<* ..V'/M Y'/M (9-24)
3t 1=1 J-l
where
k* - k r(Dnr(J)nY(i)nY(j)n (g_25]
k alj " "aid r. . C. . ( '
with initial condition
Y'(a)n(t0) =1 , all a and n (9-26)
Note that (9-24) is a coupled system of I nonlinear ordinary differential
equations. Solving this system for y1 for each species, in each layer, and
at each grid point of the model domain we obtain 1 of the 3 factors
A
in the solution = Y*Yr' The second factor y is derived from the
solution of the system of 3 linear ordinary differential equations (9-17)
for each pollutant at each grid point. And the final component r is ob-
tained by numerical solution of Equation (9-5) for each pollutant and within
each layer over the entire modeling region.
Numerous techniques already exist for treating each of these 3 sets
of equations but we seek new ways of treating them that possibly are better
suited to our specific modeling needs. In the sections below we develop
these alternative methods of solving the r, Y» and Y' equations given above
C(9-5), (9-17), and (9-24), respectively].
Solution of the r Equation, (9-5)
In formulating the equation for r, (9-5), from equation (9-4), we
collected only those terms that describe horizontal transport and diffusion.
154
-------
Chemical processes, vertical'. dispersion and emissions were relegated
to the Y equation. As a result (9-5) is a linear, homogeneous partial
differential equation with a solution of the form
r(r.t) = rfr'.^Mr.tlr'.t^dr' (9-27)
where p is the Green's function of (9-5). In the present instance we have
P(x,y,t|x',y',t J - -*— exp[- -i- (x-x'-x)2 - -i- (y-y'-y)2] (9-28)
0 2ira2 2a2 2a2
where
a2 = 2KH(t - t0) (9-29a)
*
. x(t) = u[x' + x(t'),y' + ytt'Kt'jdt1 (9-29b)
Jt
y(t) = j v[x' + x(t'),y' + y(t'),t']dt'. (9-29c)
We wish to emphasize that (9-5) is a specific form of the general equation
|£ + LT = 0 (9-30)
where L is an operator describing the transport and diffusion processes.
The solution form (9-27) still applies except in this case p is the Green's
function of (9-30). For reasons discussed in the previous section, we use
K-theory in our model to describe subgrid scale turbulence only; ensemble
average material fluxes, associated with the flows in W, are estimated in
a direct manner that does not require the gradient transfer assumption.
Thus, for each member of the flow ensemble W, the effects on cell averaged
concentration ., j=l,2,3, of subgrid scale variations in the velocity and
j
concentration fields are approximated by the gradient transfer expression
(7-2) with the diffusivity KH given by (7-3).
155
-------
Keep in mind that Equation (9-5) and its initial condition (9-6) apply
only to discrete-time intervals. If t and t, denote the beginning and end
of one of these intervals, solution of (9-5) and (9-6) yields r(tj). Within
the same interval the y equation (9-7), (9-8) is solved concurrently to give
y(t,). The product c(t,) a r(t,)y(t,) then serves as ~'-° initial condition
(9-6) for r in the next time interval tj •»• t2 [cf. (9-. ,j. 'We will use the
general solution (9-27) as the basis of a scheme for deriving r at the re-
quired time intervals.
In the computer model we will represent r(r,tn) on a fixed grid net-
work. A grid system that moves with the wind would yield more accurate
descriptions of the time evolution of r but such networks pose problems
operationally, particularly at the boundaries of the region and when the
wind fields are horizontally divergent, as they are in the problems of
interest to us. Thus, in the present study we will utilize a fixed grid
network.
Let the coordinates of the lattice point (I,J) be given by
x = IAX
y s JAy
where (AX, Ay) is the grid mesh dimension; and let
tn = nAt
denote the elapsed time between the initial instant t and the end of the
n-th time interval. Since (9-27) is applicable to any interval for which
the initial value of r is known, we have
n
r(iAx, JAy, nAt) = J Jr(x',yl,tn.1)p(IAx, JAy, tjx',y'.t^Jdx'dy' (9-31)
156
-------
To evaluate this expression we must express r(x,y,tn-1) as a continuous,
rather than a discrete, function. Regardless of whether we use polynomials
or trigonometric functions to represent the discrete set of ^t^j) values,
the integrals that result in (9-31) with p given by (9-28) can be evaluated
analytically to yield algebraic expressions for r(tn) at ---h grid point.
We can illustrate this most easily by considering a one d. =n:-ional problem.
In this case the counterpart of (9-31) is
r
r(IAx,nAt) = r(x',(n-l)At)(— i-)exp [- — (x-x'-uAt)2]dx' (9-32)
2o2
We have assumed here that the time interval At is small enough that (9-29b)
and (9-29c) reduce to the form
x = uAt (9-33a)
y = vAt (9-33b)
This assumption is adopted here to keep the math simple. In our general
scheme, x and y are calculated from (9-29) and this gives the method an
"upstream" characteristic that ensures unconditional computational stability.
The exponential term in (9-32) has a maximum value at the point
x1 = x = x - uAt = IAx - uAt (9-34)
and it has values significantly different from zero over an interval
AX' - 4c centered at this point. These facts indicate that the continuous
functional representation we use for r[x' ,(n-l)At] should be most accurate
within the interval x - 2a _< x1 <_ x + 2a.
Let X|< = KAx denote the grid point (on which the r values are stored)
that is nearest to the point x given by (9-34), and let
157
-------
e = XK - x - (K - I)Ax + uAt (9-35)
denote the separation between points XK and x. By virtue of the definition
of x we have
lei 1 W2 (9-36)
Introducing the new coordinate
n = x1 - XK = x-1 - x - £ (9-37)
into (9-32) we have
f*
r(lAx, nAt) = r(n * L (n-l)At)(— ^-)exp [- (n * g)2]dn (9-38)
2o2
Suppose now that we represent r[x',(n-J)At] by a quadratic expansion
about the point XK; that is,
r(x', (n-l)At) • a + b(x' - XK) + c(x' - xK)2 (9-39)
f
In terms of the grid point values of r, the coefficients of the expansion
(9-39) are found to be
_ a ' rK,n-l
b *
c =
where
rK,n-l " r(KAx» (n-DAt) (9-41)
Substituting (9-39) into (9-38) we get
]dn (
158
-------
The integral here is one of the many t^ulated by Gradshteyn and
Ryzhik (1965). (Note that it is equivalent to the expression for the
moments of a normal random variable.) Thus, we are able to reduce (9-42)
to the simple, closed form
ri,n = a - £b + (*2 + £2)c. (9.43)
Making use of (9-40) we have
ri,n " rK,n-l * 1
(9-44)
where £ is given by (9-35) and a2 = 2KHAt.
An interesting aspect of (9-44) is that if we set K = I (and if a2 = 0),
then c = uAt and we have the well known Lax-Wendroff finite difference
approximation of the advection equation, which was derived originally by a
method different from the one we have used.
The Lax-Wendroff scheme is known to be computationally stable only if
|u| <, AX/ At. From the perspective of the method we have used to derive
(9-44), it is easy to see why this criterion must be satisfied. If |u|
exceeds Ax/ At and yet we use K - I, as is done in the Lax-Wendroff scheme,
we introduce a systematic error into the calculations by representing r in
the vicinity of X'=KAX by a quadratic expansion about the distance point X'=!AX.
It appears that the method we have introduced above to derive difference
schemes inherently guarantees stability, consistency, and certain conser-
vation properties. This ability combined with its ease of application to
multi -dimensional problems makes this procedure particularly useful in
model development. Numerous finite difference schemes are available for
159
-------
1-dimensional problems but many lose their stability, consistency and con-
servation properties when applied straightforwardly to several dimensions.
For example, Leith (1965) showed that if the Lax-Wendroff representation of
the space derivatives is applied with forward time differencing to the 2-D
advection equation
f£*«£+vf=0 (9-45)
the resulting finite difference approximation is unconditionally unstable.
The stability properties possessed by 1-dimensional schemes can be preserved
in several dimensions by using either the alternating direction or the time
splitting techniques. Both of these methods are very popular but they are
not well suited to our needs. The principal difficulty is that both methods
require that the dependent variables at each grid point be treated twice
(3 times in 3-D simulations) in order to advance these variables to the
next time level. Thus, if the available computer memory is not large
enough to accomodate the entire model domain (our model exceeds the memory
capacity of EPA's Univac 1182 by a sizeable margin) costly I/O operations
are necessary to implement these schemes. The most desirable difference
scheme is one that allows small portions of the model domain to be processed
at a time, particularly subregions small enough to fit within the high-speed
memory of the computer; and that require only 1 sweep of the entire domain
to advance the dependent variables 1 time level. Difference schemes gene-
rated by the above method possess these properties.
Thus far, we have developed a series of 3 schemes using this technique
each of which is applicable to 2-dimensional equations, specifically (9-5).
Each permits the dependent variables to be stepped forward in time in a
160
-------
spatially piecewise fashion. These schemes have differing orders of accuracy
but each is unconditionally stable and the computer memory requirements of
each is independent of the size of the modeling region or the number of
dependent variables. We illustrate below the derivation of the lowest order
2-D scheme and comparative results obtained from test studies conducted with
it and the next higher order schemes in the sequence.
Recall in the earlier demonstration of the method using a 1-dimensional
problem that we expanded the dependent variable r in a quadratic series
about the point x, given by Equation (9-34), where the Green's function p
[see Equations (9-31) and (9-32)] has its maximum value. To derive an
analogous 2-D scheme we proceed in the same way except in place of the
quadratic approximation (9-39) we use the 2-D, biquadratic expansion
I7- (9-46)'
where (n,c) are the space variables (x.y) transformed as in (9-37). The
coefficients Anm in (9-46) can be determined using the Lagrange interpo-
lating polynomial
3 3
n (n-n,) n U-O
3 3 j=l J k=l k
r(n,0 = r r (r- ) (9-47)
rn-1 n=l nm 3 3
n (nn"ni^ n -
j=l n J k-1
where
n2 - C2 = 0 (9-48)
n3/Ax = 53/Ay = 1
161
-------
and r^ = r(KAx + nn, LAy + cm. NAt). Here (KAx, LAy) denotes the grid
point nearest [(lAx - uAt), (JAy - vAt)] where the Green's function (9-28)
has its maximum value [see (9-34) and the ensuing discussion]. This and
other aspects of the representation (9-46) and (9-47) of r(n,c) are illu-
strated in Figure 9-1.
Using (9-46) and a procedure similar to that which led to (9-38) we
obtain
33 f" r , m. , n
r(iAx, JAy, (n+l)At) =^ I I A (^ (|-)
n=l m=l J— J-- ttx oy
'd^1 (9-49)
where
n1 « nAx c1 * SAy
<*0 - (I-K)Ax - uAt (9-50a)
60 = (J-L)Ay - vAt (9-50b)
and
202
(C'-BJ2
exp[ -- ^— 3
2a2
162
-------
Y
i
IT:
13
AY
trajectory segment over a
period At ending at (I,J)
at time NAt
rigure 9-1.
Illustration of the grid, coordinate system, and other
parameters on which the 2-D representation (9-46, 47}
of the function r(n,0 is based .
163
-------
The integrals in (9-49) can be evaluated analytically. We get
1
2
27TU
m = 0
m = 1
m = 2
J — CO J — CO
n = 0
1
ao
AX
(AX)2
n = 1 n = 2
Ay (Ay)2
a |3 a(82+a2)
00 00
AxAy Ax(Ay)2
(o o\ /o o\/o *^\
a 2 + a2)g (a 2 + a2)(g 2 + a2)
(Ax)2Ay (AXAy)2
Making use of these results in (9-49) we obtain
r(lAx,jAy,(n+l)At) = A + A10ao + A01go + AllaoBo
Ax Ay AXAy
(Ax)2Ay Ax(Ay)
(AX)2 "(Ay)2 (AXAy)2
The coefficients A,,m are found from (9-46) and (9-47) to be
nm
Aoo ' F/2
A1(J - (M-C)/4
A01 s E/2
A = (H-B)/4 (9-52)
164
-------
A21 = (B-2E + N)/4
A12 = (G-AJ/4
A2Q = (C-2F + M)/4
AQ2 = D/2
A22 = (A-2D + G)/4
where
A = T'U - ZT'l2 + r;3
B = rJ3 - ril
c = 2r;2
0 = r21 ' 2r22 + F23 (9-53)
F = r' - r'
t r23 r21
F - 2r^2
G = r31 * 2r32 + r33
H - r ' r'
H ' r33 r31
M = ZT'32
The other 2 difference schemes that we have developed are derived in
exactly the same way as the biquadratic scheme (9-51) above except they
are based on bicubic and biquintic expansions, respectively, of the r field.
Whereas (9-51) uses r values at 9 points to estimate r(I,J,n+l) [note that
there are 9 terms on the right side of (9-51)], the bicubic scheme employs
16 points and the biquintic scheme 36 points. The last scheme requires
more computer time to operate than either of the 2 lower order schemes,
but by streamlining the algebraic operations executed by the computer
165
-------
code, we have reduced the CPU time requirements of the biquintic scheme to
levels of the same order of magnitude as the time required by other differ-
ence schemes in current use. We will show this later.
Because they are based on odd ordered polynomials, both the bicubic
and biquintic schemes possess much different (and more desirable) properties
than the biquadratic (9-51). We illustrate several of these in Figures 9-2
and 9-3 where we show the results of a simulation of a conic-shaped r dis-
tribution advected by a flow in counter-clockwise, solid body rotation about
the center of the square model domain. Tha angular speed of fluid rotation
is to = (40At)"1 s-1 and the eddy diffusivity KH = 0. The initial r dis-
tribution is
r(x,y,o)
1 - R/4A , if R < 4A
0 , otherwise
(9-54)
where R = [(x - 2lAx)2 + (y - ISAy^J^ and we assume AX = Ay = A. The results
of our solution of Equation (9-5) in a 30A x 30A domain using the parameter
and initial valuc> given above are shown in Figures 9-2a,b and 9-3a,b. The
former displays the results after 50 and 100 steps using the biquadratic
scheme and Figure 9-3a,b exhibits the corresponding solution given by the
bicubic scheme. In both cases A = 1 = At.
An obvious deficiency of the biquadratic scheme seen in Figures 9-2a,b
is its generation of "ripples" of positive and negative r values that appear
as a "wake" behind the moving conical cloud. Actually, the amplitude of
the ripples is exaggerated in the Figure by the heavy contours and tic
marks drawn along the lines of r = 0. Within the regions of negative r,
the amplitudes are mostly of the order of -0.02; but there is a region of
r * -0.1 immediately behind the cloud that has deepened to r = -0.13 after
100 time steps.
166
-------
0
tt
UJ
a.
10
a
a: o
t- UO
UJ
z a.
O ul
u r-
10
O
3 IU
o z:
_i •-•
O H-
Figure 9-2a.
Results after 50 time steps from a simulation by the
biquadratic scheme of the advection of a conic shaped
cloud in a, rotating flow. Angular speed of rotation
to = (40A)"1 - n
and
0.
167
-------
Figure 9-2b. Same as 9-2a, except results are after 100 time steps
168
-------
bJ
>•
IU
U
UJ
-------
to
Ul
o
Ul
o
o
It!
t—
to
O H-
Fiqure 9-3b. Same as 9-3a, except results are after 100 time steps
170
-------
Although the bicubic scheme also gener--es negative r values, they are
of much smaller amplitude and are confined to a ring that completely sur-
rounds the moving cloud. This is apparent in Figures 9-3a,b. More impor-
tantly, the size of the domain of negative numbers and the maximum magni-
tude of the values do not increase with time. In contrast, the crests and
troughs of r values generated by the biquadratic scheme increase slowly
in magnitude with time and successively more and more crests and troughs
are formed. We should add that many tests of difference schemes reported
in the literature refer to these negative concentrations as "bad numbers"
and refrain from showing them. We include them in this discussion because
they can have serious effects in air pollution simulations that treat non-
linear chemical reactions.
Another superior feature of the bicubic scheme is its maintenance of
axial symmetry of the moving cloud. The biquadratic scheme shows a slight
tendency to elongate the cloud along its direction of motion.
The biquadratic excels in its ability to preserve the peak concen-
tration in the cloud. After 100 time steps the highest concentration in
the cloud simulated by the biquadratic method is rv * 0.8 whereas the bi-
cubic method gives r = 0.7. Both schemes preserve total mass to within
4
±1 part in 10 out to 100 time steps, and apparently indefinitely.
We have found that the accuracy with which the differencing scheme
preserves both mass and symmetry is quite sensitive to the accuracy with which
the point (x, y) is specified. Recall that this is the point where a par-
ticle located at a given grid point (I, J) at time t = nAt was located at
time t = (n-l)At. In the vast majority of advection-diffusion differencing
schemes, the point (x, y) is assumed implicitly to be the point of inter-
171
-------
section of the straight line of slope v/u that passes through (I, J) and the
j,
circle of radius (u2 + v2)2At centered at (I, J). This assumption is good
as long as At and/or the flow speed are small; but in schemes such as those
derived above that allow arbitrarily large values of At, the estimation of
(x, y) by a simple linear extrapolation can cause significant errors in mass
conservation and symmetry preservation. To circumvent these errors, we
convert the flow speed (u, v) at each grid point and time step into the corre-
sponding (x, y) array. This is done outside the model using an algorithm of
high-order accuracy to evaluate (9-29). The array (x, y) is then used as
input to the model, rather than (u, v).
Another difference between the 2-D schemes derived from the technique
introduced above and most of those in current use is the neglect in the
latter of cross product terms in the polynomial representation of the de-
pendent variable. Examples of these terms are those with coefficients A,,,
A21, Ajg and A22 fo the biquadratic expansion (9-46) of r(n, c). The errors
incurred by the neglect of these terms are not revealed by the popular test
of advection differencing schemes used in Figures 9-2 and 9-3 that simulates
clouds of axially symmetric shapes in a rotating flow, because the cross
product terms in question are nonzero only when the cloud is asymmetric.
We illustrate this below in simulations of the advection of a cloud of ellip-
soidal shape in a rotating flow.
Figure 9-4 shows the results of tests of 3 differencing schemes: the
biquintic scheme (Q) derived by the procedure illustrated above with (9-46)
replaced by the 36-term biquintic polynomial; the upstream cubic spline scheme
(S) reported by Mahrer and Pielke (1978); and the fourth-order flux corrector
scheme (Z) of Zalesak (1979). All 3 are explicit schemes. Each of the
172
-------
• i. (U
2 -a
<4-i— O M-
O «*- i— IO C •
(J •»"•
V) Ol "-> M
C OJ O
.C 10 ra 4->
O 4-> 0)
O C -i- -C
C «r- 4^ «4- O •*->
O) O O O
S- "O CD t-i C
a; 3 to j- -i-
^ O I QJ *
**- i— 10 c -o 10
•r- ,-^ •*-> j= i—I
t/1 Q. J_ i.
C -r- OJ ««- »
Or— «t- S- O 00
•i- i— M- OJ
•M OJ -i- C. C -0
> +J
E 03 QJ fl3 O'
to o a.4-> o vi
V) L QJ
«*- C -r- C E
o o -o •<- a) Q_»— O O
ai
s.
01
173
-------
o
o
«3-
4>
0)
S-
174
-------
L-,
•o
-------
o
u
I
Ot
0)
176
-------
Q)
•x
'X-.
0)
Q.
O
O
0)
=3
cn
177
-------
5 panels of Figure 9-4 shows a different cross-section of the cloud (indi-
cated in the upper right corner of each panel) after 1 complete revolution
about the axis of fluid rotation (also indicated in Figure 9-4). In the
tests of the Q and S schemes, which are unconditionally stable, 1 revolu-
tion was completed in 100 time steps. A slower rotat"' rate was required
with the Z scheme to maintain computational stability: ,_ -therefore in its
test 1 revolution marks the completion of 150 time steps.
The results presented in Figure 9-4 show clearly the superior ability of
the biquintic scheme to preserve symmetry, shape, phase, amplitude and other
important properties. The upstream cubic spline distorts the shape and orien-
tation of the simulated cloud but the errors are somewhat smaller than those
generated by the Zalesak scheme (Z). The errors generated by both of these
schemes raise the question of whether either could provide a viable basis
for simulating advection in regional models of photochemical pollutants, where
the chemical reaction rates are nonlinear functions of concentration.
The principal asset of the Zalesak scheme is that it preserves the posi-
tive definite property of concentration. In both the Q and S schemes, nega-
tive concentrations are produced. However, it is apparent from Figure 9-4
that to maintain positive amplitudes, the Z scheme induces sizeable, non-
symmetrical errors whose effects in some applications, such as simulating
photochemical pollutants, are larger than those incurred in schemes like the
Q and S schemes simply by setting negative amplitudes to zero. The maximum
negative values generated by the biquintic scheme are about 3 percent of
the peak amplitude of the cloud, depending on the width of the cloud rela-
tive to the grid dimensions, AX, Ay.
178
-------
The computer time requirements of each of the 3 schemes are summarized
in Table 9-1. The time shown for the biquintic scheme is based on the opti-
mized code, mentioned earlier; and it includes simulation of both advection
and diffusion. The S and Z schemes treat advection only.
TABLE 9-1. COMPUTER TIME REQUIREMENTS FOR THE Q. S AND Z SCHEMES
Scheme CPU time (UNIVAC 1182) per grid point per time step
Biquintic (Q) 8.8 • 10"4 s*
Upstream ,,
cubic spline (S) 7.0 • 10~* s
Zalesak fourth 4
order flux corrector (Z) 3.0 • 10 s
* Advection and diffusion combined. S and Z schemes treat advection only.
179
-------
Solution of the Yn Equations, (9-17)
If the time interval At that we use as the time step in solving the
r equation is subdivided into periods small enough that changes in y^, the
"chemistry variable", are negligible, then within each of these subintervals
we can treat Equation (9-17) which govern rn in each of the 3 layers
n=l,2,3 as a system of linear, ordinary differential equations. We should
point out here that with the grid cell sizes we anticipate, the time scale
At of the horizontal transport and diffusion process is much larger than
that of vertical mixing, which is implicit in the coefficients b** of the
nm
**
Yn equation; and that both time scales are generally larger than that of the
fastest chemical reaction processes described in the yn equations (9-24).
With the concentration split into the 3 phenomenological
A
components rn, yp and y^, these intrinsic differences in time scales can
be handled separately and, moreover, locally in the simulations. This
permits optimal efficiency in solving the governing equations because
those phenomena that are changing slowest can be treated with large time
steps, faster phenomena can be simulated using smaller time steps. And
in regions of the modeling domain where approximate equilibrium prevails
among the various phenomena, larger time steps can be employed than in
those areas where transient processes are taking place.
Rather than approach the solution of Equation (9-17) by finite difference
methods, we propose instead to exploit the pseudo-linear character of these
equations during the subintervals in At to obtain analytic solutions. Aside
from the obvious advantage of precision that the analytic solutions afford,
there is the added asset that time steps can be chosen as large as the
180
-------
chemistry variable y1 will allow. We pointed jt earlier that during
nighttime hours when vertical fluxes among the model layers are negligible,
large gradients in material concentration can arise. Once vertical mixing
begins following sunrise, large transient fluxes of material result as
concentration levels proceed toward the well-mixed state. These transients
could cause computational instability in numerical solutions of (9-17);
but by using the analytic solution, we can use time steps of arbitrary length
without adverse effects. We will demonstrate this feature later.
Detailed mathematical descriptions of each of the coefficients bnm that
**
enter into the definition (9-19) of the coefficients b of the j equations
(9-17) are given in Appendix C. Also provided there are expressions for
the inhomogeneous terms gn that appear in (9-17). For notational conven-
ience we will represent (9-17) here by
H
(9-55)
First, we look for solutions of the homogeneous equations, (9-55),
with Gj = G2 = G3 = 0, of the form
x = kxext , y = kyext , z = kzext (9-56)
Substituting these into the homogeneous equations, expanding and collecting
terms we obtain
181
-------
kx(x + Dj) + kyD2
kE + k(x + E) + kE
xl
kxFl + kyF2
z3
F3)kz
This system of equations has a solution only if
= 0
X
E
F
+ Dl
1
1
D
X +
F
2.
E2
2
0
E3
X +
F3
whose roots are the eigenvalues of (9-55). Here
• F3 + E2 + Dl
q =
r = D,(E9F- - E-F-) -
- E,F.)
(9-57)
(9-58)
Expanding this determinant we obtain the characteristic equation
A3 + px2 + qX + r = 0 (9-59)
(9-60)
(9-61)
(9-62)
Due to the physical nature of the system that (9-55) represents, all 3
of the roots of (9-59) must be real and negative.
The characteristic equation (9-59) can be solved rather easily. First
we use Newton's iterative method to find the root nearest zero. With this
technique the estimate of a given root after n iterations is
^ « x* , S=*_ (9-63)
where f(x) represents the left side of (9-59) and f'(x) = df/dx. Starting
182
-------
(9-63) at X* = 0 to find the root nearest 0 we find from (9-59) and (9-63)
x*l = x*0 " r/q = " r/q
(9-64)
With only 4 more iterations of (9-63) we obtain an estimate of 1 root of
(9-59) that is accurate to 8 significant figures. Let us call this result
Xj, i.e., x| = \i. Then the remaining 2 roots X2 and X3 are found by
synthetic division of (9-59) to be
C- (P
/(P
- 4q -
xj
(9-65)
and a similar expression for X3 except for a minus sign before the square
root term.
Now for each eigenvalue X there is a system of equations (9-57) for the
components of the eigenvector k. For example, the components of the eigen-
vector associated with \i satisfy the equations
= 0
(9-66)
We can choose one component of this system arbitrarily. If we pick
k = - 1 (9-67)
then
+ PI
(9-68)
F +
F2],
(9-69)
provided that D2 and (F3 + x:) are not zero
Through a similar process we obtain k , kv , k , etc. Thus, the general
x2 x3 y2
form of the solution of the homogeneous system (9-55) is
183
-------
X = VjXj + V2X2 + V3X3
y =
Z =
v2y2 + v3y3
V2Z2 + V3Z3
(9-70)
where x: = k eXl , = Xl
AJ
mined by the initial conditions on (x,y,z). Following ,~jj\a!\ (1962) we
= k el , etc; and the v 's a- constants deter-
n
can solve the inhomogeneous system (9-55) using (9-70) with the v 's treated
as time dependent variables. Substituting (9-70) into (9-55) and making use
of the fact that (9-70) is the solution of the homogeneous system (9-55)
we get _
dv2 dv3
*i
X2
dt
X3
dt
dt
dv3
dt"
dt
dt
zidT +
Solving this system first for dvj/dt we obtain
where
dv'
k - k. k ) + G3(k k - k k )
xy xy
23 32
and
Xl y2 Z3 2
k,. ] - k [k
"„ ]
(9-71)
(9-72)
(9-73)
(9-74)
184
-------
Hence, from (9-72)
dt =
_-lfie-*it , if A! f 0
\ K e
dt
, otherwise
(9-75)
Solving for v2 and va in a similar manner and subst * 3 the results into
(9-70) we obtain finally for the general solution of (9-55)
* • '.j BI"^'* + KiW (9-76'
with similar expressions for y and z except with k . and k . replacing k ..
jr I Z I XI
In this equation
] (9-77)
t , otherwise
and the C^ 's are constants that are evaluated by setting t = 0 in (9-76)
and solving the resulting set of equations with x, y, and z assigned their
initial (t = 0) values. We will not write the expressions here for these
constants nor will we discuss any of the operational aspects of determining
the k's that enter in (9-76) in special situations (such as the case of
double roots of the characteristic equation). We have developed a computer
code for evaluating (9-76) which we use in the regional model to solve the
A
Yn equations (9-17). Below we illustrate some results of exercises per-
formed with this code which, in essence, constitutes a 3^ layer, 1-dimen-
sional diffusion model.
Returning to our original notation, we can express Equation (9-17) in
the concise form
185
-------
dY
~ + BY - G (9-78)
A A * ^
where y = (YI> Y2» YS)» B is a 3x3 matrix whose elements are the coeffi-
*» *•
cients of (9-17), and G = (G^ G2, G3). Figure 9-5 shows solutions we
have derived for 3 sample cases. The first, depicted in Figure 9-5a, is
the case of no sources, G = 0, no material sinks, and layers of equal
^
thickness. There is vertical mixing between each of the 3 layers but
mixing across the interface between layers 2 and 3 is only one-half as
strong as across the surface between layers 1 and 2. Starting with initial
values YI = Ya = 0, Y2 = 10. the evolution of y(t) is as shown in the Figure.
After a time t = 15, layers 1 and 2 are well mixed, but the concentration
in layer 3. is smaller reflecting the weaker mixing between layers 2 and 3.
In the limit as t •*• », the material initially present in layer 2 is equally
A
partitioned among the 3 layers, as is evident in the y values shown for
time t = 500.
Incidentally, the curves shown in Figure 9-5(a), (b) and (c) were ob-
tained by recursive application of the algorithm we developed above for
solving (9-78). That is, we solve (9-78) for y at t = 1 using the initial
conditions. This result is used as the initial condition in a second appli-
A
cation of the algorithm in which we solve for y at time t = 2. It is not
*
difficult to show that this second value of y is equivalent to the solution
of (9-78) at time t = 2 using the given value of y at t = 0 as initial con-
dition. The use of a small time step like At = 1 would be essential in
solving (9-78) by finite difference methods under the conditions given in
A
Figure 9-5(a) to reproduce accurately the yp values in the interval 0 < t < 10
A
where changes are occurring most rapidly. In fact, since dy2/dt = 1.5
186
-------
'.1 -.1 0
B-I-.1 .15 -.05
0 -.05 .05
G4o°)
"\ol
•n
500
b.
10 -
10
0
i(o)- I o
1 o;
-.1 0
B- (-.1 .2 -.1
,0 -.1 .1
10
[ 0
\10
&• same as in b.
G= 0
~ I 0
—• «. _ Total Mass
500
Figure 9-5. Sample solutions of the y equation (9-78) .
187
-------
initially, a time step larger than about 7 in a finite difference solution
A
of (9-78) would produce negative values of Y and possibly computational
instability. Using the analytic solution one can obtain exact solutions
for arbitrary values of t with essentially a single iteration of the algo-
rithm. Of course, more computations are required than in a single time
step of a numerical solution technique, but the total computer time require-
ments are certainly comparable, and the flexibility and accuracy provided
by the analytic method are most valuable.
In the second sample calculation presented in Figure 9-5(b), a source
of unit strength is located in layer 1 and there is also material deposition
in this layer (deposition rate = O.lyi). Mixing of equal strength occurs
between each of the layers and initially no material is present. The
solution y(t) shown in the figure displays the expected behavior, with
A A A
YI > Y2 > Y3 due to the source in layer 1. The steady-state is reached
in this system when the concentration YI has become large enough that the
A
deposition rate O.IYI is equal to the source emission rate GI = 1. Our
A A A
algorithm gives the correct steady-state solution YI = Y2 s Y3 = 1°> as
indicated in the figure for t = 500.
The final calculation shown in Figure 9-5(c) is for the same case just
presented except instead of a material source, we begin with a puff of
concentration 10 initially in layer 3. The results show that YI rises
from a value of 0 initially to a peak value of about 1.4 at time t = 20
and then declines after that. The total mass of material present, indicated
by the dashed line in the figure, is almost constant until t = 5, but after
that it begins a rather rapid decline as material filters into layer 1
and is subsequently removed by deposition.
188
-------
Using the 3 case studies just described and others, we have validated
the accuracy of the algorithm we have developed above for use in solving
A
the Y equations, (9-17). This algorithm is currently operational in our
regional model.
Solution of the Y'/ \n Equations (9-24)
Since the y and y' equations (9-17) and (9-24) are coupled, they
must be solved simultaneously. (By contrast the r equation (9-5) is vir-
tually independent of both y and y1.) The method we developed above for
solving (9-17) gives y at the end of time intervals of arbitrary length
J- -Ir
provided that the parameters b and Gn that enter (9-17) are approximately
constant during each interval. Both of these parameters are functions of
y', [see (9-19) and (9-55)], so the time interval used to evaluate yp is
determined by the rate of change of y'.
The system of equations (9-24) that governs y/ \n [recall that the
subscript (a) refers to the pollutant species] is nonlinear. Despite its
relatively simple form this system is particularly difficult to handle
numerically; because being representative of the chemical reactions that
Q
occur among air pollutants, some eigenvalues of the system are 10 times
larger than others. We develop below a numerical technique that appears to
be able to handle this particular system. It is computationally stable; it
does not produce negative concentrations, as most numerical schemes do when
the constants k .. have greatly disparate values; and it provides easy control
Q t J
over the accuracy of the solutions.
To develop this method we first write (9-24) in the equivalent form
3y' I II
«£* Zk .y'y! + I Z k ..y'.y!. (9-79)
at .=1 oaiWi . . . a
j)
189
-------
where I represents the total number of chemical species present. For no-
tatlonal convenience we have dropped the layer reference subscript n from
Y/ x and use the subscript to refer to the pollutant species. The first
summation on the right side of (9-79) represents all reactions of the type
aai
a + i -" t + m (9-80)
These reactions destroy species a by transforming it into different
species JL and m. It follows then that k . is negative.
act i
The last term on the right side of (9-79), which by design does not
contain y'» represents all reactions of the type
i + j •*• a + n (9-81)
These reactions produce species a from other pollutant species and hence
k . . is positive. Having noted these aspects of the summation terms in
(9-79) let us now write this equation in the more abbreviated form
where
P
(9'84)
Consider now a time interval tQ <, t <_ ti for which the values yj(t )
are known for all i = 1,2,... I pollutants. And let PaQ and QQO be the
values obtained from (9-83) and (9-84), respectively, using the given values
of yj at tQ. If PO and QQ were constant in (tQ, tj), Equation (9-82) would
have the solution
190
-------
- Qa0/Pa0] exp [- Pafl(t - tfl)] (9-85)
Obviously the degree to which this expression approximates the general solu-
tion of (9-82) is determined by the rate at which P and Q are changing
in (t0, ti).
Suppose that from the set of solutions (9-85) for *
-------
smaller than, or equal to the initial value
Y;(t0) -= Y;O (9-88)
As we noted earlier, (9-85) is a good approximation of the solution of
(9-82) only for times t such that
«
AYJ = |YJ(t) - Y!(t0)| lXY!(t0) • (9-89)
where A is a small fraction. The "smallness" of x is determined by what
we define a "close" approximation of the solution of (9-82) to be.
The limit on t imposed by (9-89) can be found by substituting (9-85)
into (9-89) for yj(t). We get
AY- = |YJ(tQ) - ^ICl - exp(-P1o(t -t0))] i AYJ(t0) (9-90)
or
1 - exp(-P1o(t - t0)) < }*|u°).tl| 'F, (9-9D
If F.. >_ 1 then any value of t >_ tQ satisfies (9-90); but when Fj < 1, we
require
ti - t0i- ln(l - F^/P^ (9-92)
Thus, a sufficient condition for (9-85) to be a close approximation of the
solution of (9-82) for all I species in the interval t <_ t <_ tx is that
ti satisfy
ti - tfl + min { - tn(l - F^/P^} (9-93)
We should add here that in evaluating F^ using (9-91), XYJ(t ) is replaced
192
-------
by X5i if YJ(tQ) < cr
In the studies we have performed to date with this method of solving
(9-82) we have found that the time intervals At = t - t , given by
n n n-1 3 J
(9-93) are unnecessarily small when the same value of X is used for all
species. Our available evidence is that x can be made much larger for
pollutant species like oxygen atom that are present in extremely small
concentration than for the'more plentiful species like NO, 03, etc.
For example, we used the scheme above to simulate the 4-species
mechanism
N02 + NO + 0
0+02+M + °3+M (9-94)
k3
with rate constants
ki « (21/36) • 10"2 s'1
M02)(M) -^- 10*s~l
k3 = (4/9) • 10"1 ppo"1 s'1
o
Note that the effective rate constant k2(02)(M) is about 10 times larger than
We solved the rate equations for this system
) - k3((U(NO)
VI V ^ v
etc, using a value x = .1 for each species.' The inital conditions were
chosen arbitrarily as NO = 100 ppb, N02 = 500 ppb, 03 = 300 ppb, 0 = 10" ppb.
193
-------
The results are shown in Figure 9-6a. We repeated this simulation using
the same initial conditions but with x = 0.1 for NO, N02, and 03 and
X = 0.4 for oxygen atom, 0. The results obtained in this calculation for
t = 100 were within approximately 0.1 percent of the values obtained in the
first run (shown in Figure 9-6a); yet the computer time required was only
2/3 that required in the first case. The time difference was due to the
fact that 18 time steps were needed to reach t = 100 with x = 0.1 for all
4 species whereas only 12 steps were required when the X value assigned to
oxygen atom was relaxed to 0.4. In the former case the first 7 of the 18
time steps were smaller than 0.001 s. During these 7 steps oxygen atom con-
centration increased gradually (by 10 percent increases) from its initial
value of 10~ ppb to 4.06 • 10" ppb. During these 7 steps the concentrations
of Og, NO and N02 were virtually unchanged from their initial values. Once
the 0 concentration had reached 4.06 • 10 ppb it was within 10 percent of
its equilibrium value for the NO, N02 and 0, concentrations existing at that
time and thus infinitely large subsequent steps were permissible "with respect
to 0." Consequently, control of the time step size then shifted to NO, N02
and 0^ at this point and, due to their slower reaction rate, the remaining
time steps were larger than 1 s. The actual time intervals are indicated
by the arrows in Figure 9-6a. Also shown is the sum (NO + N02) which,
according to (9-93), should remain constant. The scheme maintains this
property to within a few percentage points depending on the value of x chosen.
The interesting feature of the simulation is that when the value of
X used for oxygen atom was increased from 0.1 to 0.4, with the other 3 x
values held the same, the model used only 1 sub-millisecond time step ini-
tially, rather than 7, yet in doing so it commited no error (the values of
0 at t = 1.027 are identical in the 2 cases).
194
-------
t(sec)
c
C0
1.C
NO2
i L
so
100
t(«c)
Figure 9-6.
Temporal behavior of the reaction system (9-94) for initial
species concentrations (a) NO = 100 ppb, N02 = 500 ppb,
0, = 30 ppb, 0 = 10 ppb; and (b) one tentn the values
used in (a). Arrows indicate the time steps used by the
model to arrive at concentrations at t = 100 sec.
195
-------
As a further check of the effect of letting x vary with the species,
we performed a second pair of simulations as above but with initial condi-
tions l/10th those used in the first set — NO = 10 ppb, N0 - 50 ppb,
03 = 30 ppb, 0 = 10 ppb. The results of this experiment with A = 0.1 for
all 4 species are depicted in Figure 9-6b. The values obtained with
X = 0.1 for NO, NO, and 0, and X = 0.4 for oxygen atom were within 1 percent
C. 0
of the values derived in th.e first example at t = 100 s. However, with this
set of initial concentrations, the simulation that used x = 0.4 for oxygen
atom took only 5/11 the computer time needed with x set to 0.1 for all
species. This reduction in machine time is considerably larger than that
realized in the first pair of experiments where different initial concen-
trations were used. The reason is that with smaller initial concentrations
of NO, NOo and 0,* the rates of change of these species are considerably
reduced (compare Figures 9-6a and b) and thus a larger fraction of the
iterations required to reach time t = 100 are dominated by the rapid varia-
tions in oxygen atoms. Incidentally, it is interesting to note in Figure
9-6 that due to the nonlinear character of the chemistry, the temporal be-
havior of all of the species is reversed when the initial concentrations are
all divided by 10.
Having found that the performance of our solution technique is much
more sensitive to the behavior of some species than others, we restructured
the method to achieve maximum execution efficiency. This work was performed
through applications of the scheme to the 23 species kinetic mechanism re-
ported by Demerjian and Schere (1979). As a result of this work, the tech-
nique presented above for solving the y1 Equations (9-24) was modified and
it now has the following form:
196
-------
1. Compute preliminary estimates P*Q, Q*Q of the decay and pro-
duction coefficients (9-83, 84) using the known values of
2. Compute the sum CnM of the concentrations of NO, N0« and 0^ at
Q.,
time t :
CQN = [NO]Q * [N02]0 + [03]Q
Note that [NO]Q = r(NQ) YfNO) , etc. (cf. 9-3)
3. Determine the time ti using (9-93) but with the function minO
I
replaced by min{>, where K denotes the subset of the I species
K
whose concentrations [i]Q at tQ satisfy
m0 1 CON/IOO
4. Using the values PJQ, QJQ from Step 1 and the y^, compute pre-
liminary estimates Y'* = r'(ti) using (9-85)
5. Using the y'* values from Step 4 and (9-83, 84), compute P* , Q*
ai ai a
6. Compute corrected, final estimates of P „, Q „ by
au (xo
P _ = P*
oO a
(These expressions are the result of empirical findings rather
than theoretical considerations.)
7. Compute the final estimates of v1 using the P and Q from
r 'aj 3 oo xao
Step 6 and the y1
8. Return to Step 1 and repeat the steps above to obtain at y' , etc.
02
197
-------
In Part III of this' report we will present detailed comparisons of the
performance of this scheme with that of the well-known Gear routine. We
will show that the solutions of (9-24) given by the technique above are
within about 1 percent of those given by the Gear method for each of the
23 species in the Demerjian-Schere mechanism over 24 h simulations.
The method developed above achieves this level of accuracy with about % the
computer time required by the Gear method. If the control parameter A in
the scheme is increased until the maximum differences between the generated
solutions and those given by Gear are of the order of 25 percent, the exe-
cution time is about 1/5 the Gear routine requirement. This difference in
machine times is due chiefly to the fact that the method above is in essence
a 2-step scheme whereas Gear is multi-stepped. In our regional model where
A
the Y and y' equations are solved alternately over short time intervals,
the need for multiple values of y' by the Gear method to produce solutions
within each subinterval results in a large computational "overhead".
198
-------
SECTION 10
SUMMARY OF MODEL EQUATIONS
For convenience we summarize below the basic forms of the governing
equations in each of the model's 4 layers.
Layer 1:
9i
—
3i , 9i
— L + ~
where
_ 1
A a cos*
= I
a = earth radius at MSL (in meters)
A = longitude
9 = latitude
A = a2 (AcfiAA) COS9
199
-------
,AX = latitude, longitude grid cell dimensions
» constants (A = - AX = %
= a2cos
X+AX/2
X-AX/2
/2
MAX (0,
i = all chemical, rainout and washout processes.
! = all emissions of c in Layer 1 (includes stacks and surface
sources above nighttime radiation inversion.
i, i = layer averaged horizontal wind components
(u = east-west component), (v = north-south component)
1m
6 = deposition velocity of species c
aT-(x,4),t) = fraction of surface Hj penetrated by terrain
w> - i^f + r c1 + erf (•
200
-------
w
Rl
w
Dl
w
m
WQ, - z,, neutral and unstable conditions:
~ HJ> ( = given inversion layer depth growth rate), stable
mean vertical velocity on H, (terrain induced component
excluded)
rms vertical turbulence on H, ( = 0 in stable cases)
threshold of cumulus "root" updraft veK.it/ on H,
» a » w (see Layer 2 equations)
C* w
rx
erf(x) = —
J 0
o,l
e"t2 dt
61
- *To>Fo
fraction of surface H in given grid cell penetrated by
terrain
Ozone:
Fo(03) = <03>1 W-X- " °3W+X+
- (1 - a)
, otherwise
Nitric Oxide:
FQ(NO) = ! w_X_ - NO w+X+(l - 0(1 - a)
- (1 - o)
- 0) , if
otherwise
201
-------
Nitrogen Dioxide:
FQ(N02) = 1w_A_ - N02 w
- (1 - a)
- a)
03)
NO
1 + @
N02/v
J2
1 + 6
, otherwise
N02/v
all other species, x-
P0(x) = lW.X_ - XW+X+(1 - 5)(1 - a) - (1 - a)(v?X
+ ax/v)
-1
W - w
a =
"9 "r
- \ _ oc - fe/Ae]
See Layer 0 equations for specification of 0^, NO, N02, x> w+, x+, etc.
Layer 2:
- 2 + 2
3-
2
y, TT 2
A d A
V2(x,4»,t) = a2cos<()
A+AX/2
X-AX/2
- MAX
202
-------
2
• subgrid scale fluxes of c (see Section 7)
2 = all chemical, rainout and washout processes in Layer 2
2 = source emissions in Layer 2 (if any)
2. 2 = Layer 2 averaged winds.
tfw - w 3z f
" ac
fl
(cc - 2) + 2 -^ - ^ C-
Fl,2 = a
ar Wlm' see La^eri; **D1, see (4-44c)
a s fractional area coverage of cumulus clouds
w = mean upward velocity in cumulus clouds
— = turbulent entrainment rate of inversion air into mixed layer
w2 = mean vertical velocity (mean of both terrain induced component
wT2 and divergent component wD2)
TT-
= z2 = local time rate of change of mixed layer top
<|> = fraction of cloud air from surface layer
0 1 * 1 1 and ty <_ w+X+(l - a
w - w
Vc - -
c, c1, c> w+, A+ = See Layer 0 equations
203
-------
Layer 3:
3t
3
3-EnV:
3j3_ , ,^ 3:
3 ~3X + y,t) = mixed layer top
> ,t) = model top
3
3
subgrid scale flux of c (see Section 8)
3 = all chemical (wet and dry) rainout and washout processes
3, 3 = layer averaged winds
- w
2,3
r3,3
- 3) +
32.
3 -^ , if dH3/dt £0 •
(co, - 3) dHj/dt + 3 -g| , otherwise
c^ = concentration of species c above z«
dH3/dt = given volume flux through model top surface
where c1 and c are defined with the Layer 0 variables.
204
-------
Layer 0:
<03>0 =
0
3 w+x+
o3- -
0, if SNQ >
°3V " SNO/^ , otherwise
Vv
0 = (1 - c) NO +
w_X_!
N02
.NO
^'"NO
v(NO - OJ ,f ; >
J t 'I iMn •>
NO1
1 ,,Mn
, otherwise
0 » (1 - c) N02 +
w_X_i
SNO
2 , otherwise
205
-------
Species x other than 03> NO, and
o = (1 - c)x + ex1
XV
Parameters:
X+ = *[! - erf (-2—)]
\
>Q - -± exp
.
0
[1 - erf (A-)]
= w- - Z0
v = u* (plume entrainment velocity)
5 = plume volume fraction
a = rms vertical turbulent velocity on H
wo
206
-------
REFERENCES
Artoz, M. A. and J. C. Andre, (1980): "Similarity Studies of Entrainment in
Convective Mixed Layers", Bound. Layer Meteor.. Vol. 19, pp 51-66.
Batchelor, G. K., (1950): "Application of the Similarity Theory of Turbulence
to Atmospheric Diffusion", Quart. J. Roy. Meteor. Soc., Vol. 76, pp 133-146.
Binkowski, F. S., (1979): "A Simple Semi-Empirical Theory for Turbulence
in the Atmospheric Surface Layer, Atmos. Environ., Vol. 13, pp 247-253.
Carras, J. N. and D. 0. Williams, (1981): "The Long-range Dispersion of a
Plume from an Isolated Point Source"., Atmos. Env., Vol. 15, pp 2205-2217.
Deardorff, J.. W., (1974): "Three Dimensional Study of the Height and Mean
Structure of a Heated Planetary Boundary Layer", Boundary Layer Meteor.,
Vol. 7, pp 81-106.
Deardorff, J.W., and R.L. Peskin, (1970), "Lagrangian Statistics from
Numerically Integrated Turbulent Shear Flow," Physics of Fluids,
Vol. 13, pp. 584-595.
Demerjian, K. L. and K. L. Schere, (1979): "Applications of a Photochemical
Box Model for Ozone Air Quality in Houston, Texas. Proceedings,
Ozone/Oxidants: Interactions with the Total Environment II, Houston,
TX, 14-17 Oct. 1979, APCA, Pittsburgh, Pa., pp. 329-352.
Draxler, R. R., (1976): "Determination of Atmospheric Diffusion Parameters",
Atmos. Env. 10:99-105.
Evans, R.B., (1979) "The Contribution of Ozone Aloft to Surface Ozone
Maxima", Doctoral Dissertation, University of North Carolina, School
of Public Health, Chapel Hill, NC. 301 pages + viii.
Godowitch, J. M., J. K. S. Ching, J. F. Clarke, (1979): "Dissipation of the
Nocturnal Inversion Layer at an Urban and Rural Site in St. Louis",
Preprint Volume, Proceedings of the Fourth Symposium on Turbulence,
Diffusion and Air Pollution, Reno, Nevada.
Gradshteyn, I. S. and I. M. Ryzhik, (1965): Table of Integrals Series and
Products, Academic Press, New York.
Haltiner, G. J., (1971); Numerical Weather Prediction, John Wiley and Sons,
New York, New York, p. 54.
207
-------
Harming, R. W., (1962): Numerical Methoos for Scientists and Engineers.
McGraw Hill Book Co., New York.
Hay, J. S. and F. Pasquill, (1959): "Diffusion from a Continuous Source in
Relation to the Spectrum and Scale of Turbulence", Advances in Geophysics,
Vol. 6, Academic Press.
Henderson, R.G., R.P. Pitter and J. Wisniewski (1980) "Research Guidelines
for Regional Modeling of Fine Particulates, Acid Deposition and
Visibility", Report of a Workshop held at Port Deposit, MD October 29-
November 1, 1979.
Hunt, J. C. R., W. H. Snyder and R. E. Lawson, (1979): "Flow Structure and
Turbulent Diffusion Around a Three-Dimensional Hill", EPA-600/4-78-041,
Environmental Protection Agency, Research Triangle Park, North Carolina.
Kaimal, J. C., J. C. Wyngaard, D. A. Haugen, 0. R. Cote, Y. Izumi, S. J.
.Caughey and C. J. Readings, (1976): "Turbulence Structure in the Con-
vective Boundary Layer", J. Atmos. Sci.. Vol. 33, pp 2152-2169.
Kaplan, W., (1962): "Ordinary Differential Equations", Addison-Wesley Publishing
Co., Reading, MA. 534 pages.
Lamb, R. G., (1976): "Research in Mesoscale Air Pollution Simulation Modeling,
Vol. Ill: Modeling of Microscale Phenomena. Final Report by Systems
Application, Inc. to EPA, Contract 68-02-1237.
Lamb, R. G., W. R. Shu, D. R. Durran, J. H. Seinfeld and L. •£. Reid, (1977):
"Continued Research in Mesoscale Air Pollution Simulation Modeling,
Vol. VI: Further Studies in the Modeling of Mesoscale Phenomena", Draft
Final Report to EPA Contract No. 68-02-2216.
Lamb, R. G., (1971): "Numerical Modeling of Urban Air Pollution", Ph.D disser-
tation, Department of Meteorology, UCLA, Los Angeles, CA. 289 pages + xvii.
Lamb, R. G., (1975): "The Calculation of Long Term Atmospheric Pollutant
Concentration Statistics Using the Concept of a Macro Turbulence. Pro-
ceedings of the IBM Seminar held in Venice. (Available from the author).
Lamb, R. G. and D. R. Durran, (1978): "Eddy Diffusivity Derived from a Numeri-
cal Model of the Convective Planetary Boundary Layer. II Nuovo Cimento
1C:1-17.
Lass, H., (1950): Vector and Tensor Analysis. McGraw Hill Book Co., New York,
New York, p 5T.
Leith, C. E., (1965): "Numerical Simulation of the Earth's Atmosphere," Methods
in Computational Physics, £: 1-28.
Lilly, D. K., (1968): "Models of Cloud Topped Mixed Layer Under a Strong
Inversion", Quart. J. of Royal Met. Soc., Vol. 94, pp 292-309.
208
-------
Lumley, J. L. and H. A. Panofsky, (1964): The Structure of Atmospheric
Turbulence", John Wiley and Sons, New York.
McMahon, R. A. and P. J. Denison, (1979): "Empirical Atmospheric Deposition
Parameters - A Survey", Atmos. Environ., Vol. 13, pp 571-585.
Mahrer, Y. and R. A. Pielke, (1978): "A Test of an Upstream Spline Interpo-
lation Technique for the Advection Terms in a Numerical Mesoscale Model".
Mon. Wea. Rev., Vol. 106, pp 818-830.
Siple, G. W., et al. (1977): "Air Quality Data from the Northeast Oxidant
Transport Study", EPA-600/4-77-020, Environmental Protection Agency,
Research Triangle Park, North Carolina.
Taylor, G. I., (1921): "Diffusion by Continuous Movement". Proc. London Math.
Soc., Ser. 2, 20:196-202.
Tel ford, J. W., A. Vaziri and P. B. Wagner, (1976): "Aircraft Observations
in the Planetary Boundary Under Stable Conditions", Boundary Layer
Meteorology, Vol. 10, pp 353-377.
Tennekes, H-. , (1973): "A Model of the Dynamics of the Inversion Above a
Convective Boundary Layer", J. Atmos. Sci.. Vol. 30, pp 558-567.
Thorpe, A. J. and T. H. Guymer, (1977): "The Nocturnal Jet", Quart. J. Royal
Met. Soc., Vol. 13, pp 633-653.
Van der Hoven, I., (1957): "Power Spectrum of Horizontal Wind Speed in the
Frequency Range from 0.0007 to 900 Cycles Per Hour", J. of Meteor.,
Vol. 14, pp 160-164.
Venkatram, A., (1979): "The Expected Deviation of Observed Concentration from
Predicted Ensemble Means". Atmos. Env. 13:1547-1550.
Wesley, M. L. and B. B. Hicks, (1977): "Some Factors that Affect the Deposition
Rates of Sulfur Dioxide and Similar Gases on Vegetation", J. Air Pollutant
Control Assoc.. Vol. 27, pp 1110-1116.
Willis, G. E. and J. W. Deardorff, (1974): "A Laboratory Model of the Unstable
Planetary Boundary Layer", J. Applied Meteor.. Vol. 31, pp 1297-1307.
Wyngaard, J. C. and 0. R. Cote, (1971): "Budgets of Turbulent Kinetic Energy
and Temperature Variance in the Atmospheric Surface Layer", J. Atmos.
Sci., Vol. 28, pp 190-201.
Zalesak, S., (1979): "Fully Multidimensional Flux-Corrected Transport Algorithms
for Fluids", J. Comp. Physics. 31:335-362.
209
-------
Appendix A.
Transformation of the Governing Equations to Curvilinear
Coordinates.
Figure A-l. The curvilinear coordinate system used in the model.
To transform the divergence operator (v«) into the coordinates of the
frame depicted in Figure A-l in which the basis vectors e, and e^ are every-
-A ,
where parallel to latitude and longitude circles, respectively, and in
which eR is vertically upward, we use the vector calculus relationship [see
Lass (1950)]
l < A> + 3 (W*> + a! ( W*)] (A-1]
H70T
A K
where F is an arbitrary vector with components (F,, F., FD) and h., h. and
«, A
*
(A-2)
(A-3)
210
-------
Thus, we find from Equations (A-l) and (A-3) that in the curvilinear frame,
divergence of a vector field is given by
3F, , 3F, 9F_ F, 2F_
div F »
where
__!
*COS 3X R 3ldx' )dR (A-7a)
where the cell volume V. is defined by
rX+AX
X-AX
R2cos4.'dRd4.'dx'
(A-7b)
Since the angular dimensions (AX,A<|>) of the cells we plan to use are small
with respect to unity, and since we are concerned with elevations z that are
a small fraction of the earth radius a, we can simplify (A-7) to
211
-------
j(A.*,t)
A+AA
A-AA
Zj(A,*,t)
F(A',,t) * a2cos
A+AA
A-AA
Z.
(A-9)
z.
Consider now the forms that cell averaged spatial derivatives acquire.
For notational convenience we will drop the layer designation subscripts.
Let's look first at ~~ . From Equation (A-8) we find using Leibniz1 theorem
where
1(5)
In a similar fashion we obtain
31
AA) - I(A - AA))] -^ff
22(5.*')
(A-10)
Integrate this with respect to 5
d£ • 1(5 + A5) - 1(5 - A5)
212
-------
Letting 5 = x we get
I(x + AX) - I(x - AX) = -— [F(x,*,z2)
where
n . 1
X+AX
X-AX
32.
(A-ll)
(A-12)
A =
fX+AX
X-AX
1 = a2cos(4A
3X
3V.
3X
Substituting this into (A-ll) we obtain
3
TT
3V _ A
In the governing equation, we have the term
1 3F
Rcos<{> 3X
From the analyses above we conclude that
*COS$ 3X j
Let's now look at the term
- .
otp
213
-------
From Equation (A-8) and Leibniz1 rule,
a2cos4> r-T / . ,\ , i •
—^—*• LI^(* + A; - I^(
azcos
-------
Substituting this relationship into (A-16, yields
We conclude that on averaging the governing equation the term 4 T^ becomes
K 3i
< — If. - 1. r % • c 4- • 3V. A / p £2.- 1 r 3Z •
Performing the <> averaging on the governing equations and using results
(A-14.1, (A-17) and ( 2-11) we get
3UC 1 3 3£nV A r 3Z, ,.r 11
3X acos ax acos^ 3X aVcos^ L 3X ax
3WC A l Z? Zi n
"3Z~ > = ? ^ " WC "
* 0
2cw n
. 0
Adding these terms and recalling the form (A-6) of the wind vector v, we obtain
the governing equations of layer- averaged species concentrations in curvilinear
coordinates:
215
-------
V7
dH.
c dt°
where
acos<{>
u =
* a
and
- z.
216
-------
Appendix B. Criteria for Validity of the Steady-State Assumption
in Layer 0.
We shall investigate here the error incurred by the assumption of
steady-state conditions in Layer 0 made in the formulation of the governing
equations of this layer in Section 5. Since the meteorological conditions
under which the steady-state assumption is most erroneous are those that
confine mixing to Layers 1 and 0, we will consider only these 2 layers
in the analyses below.
Ci
1
GO ;
hi
t
[ho
fff/f/'ilf/i//t/iffr
Figure B-l. Two-layer system for analysis of the steady-state
mixing assumption.
Figure B-l shows the two-layer system we will consider. The equations
used in the model to describe concentrations Cj and CQ in the lowest 2
layers have the form (for inert species)
3Cj/3t = - (Cj - cQ)aw/h1 (B-l)
- eco/ho + S/ho +
-------
where a is the rms turbulence veloci: on surface H which separates the
W 0
layers. Here we have neglected horizontal transport and sources in Layer 1
The surface deposition process is so slow that no significant error arises
in treating it with steady-state approximations, so we shall set the deposi
tion velocity 8 = 0 in (B-2). The solution of the resulting equation is
= Cl(o) + |t - 3- (S + Vo)(l - .-») (B-3)
J. i. W
.
C0(t) - cl(o) +
Sli h A
.
where
Ao = cl(0) ' co(o)
awt/h0
and where we have assumed
(B-4)
(B-5)
(B-6)
(B-7)
The corresponding steady-state (i.e., 8c /3t = 0) solutions are
clss(t)
St/h
coss(t) = CI(D)
(B-8)
(B-9)
Thus, the error committed in assuming the steady-state in Layer 0 is:
'o '
- c
1
l
1
218
-------
The steady-state solutions will have the largest error during early
morning and later afternoon hours in rural areas as a undergoes a rapid
change. Suppose that in the hours just before sunrise, a has some small
w
value OWQ and that within 1 time step At of the model simulation it jumps
to a much larger value a. From (B-8) and (B-9) we get
W
AQ =
and hence
h S
.
1w ' • n w i .».
W 1 WO
where X is
a At
- X = -£_ (B-13)
0
and At and h are given. We have already assumed that
hl »> ho
so if we multiply e, by h,, we get a measure of the total mass lost or gained
AM as a result of the steady-state assumption. Using a typical regional
scale model time step At = 600 s, h = 30 m (h, > 300 m) and a post-sunrise
value o.. * 0.3 m s, we see that A = 6 and that
W
AM = -h S/a (sunrise) (B-14)
W
After sunrise there is a mass loss from the system equal to the quantity
emitted by the sources in the time interval h./allrt required to traverse the
0 WO
depth of Layer 0 moving at a speed aWQ. If hQ = 30 and aWQ is only several
centimeters per second, AM can be quite sizable. However, the source emission
rate prior to sunrise is usually a minimum and, moreover, the occurrence of
a large increase in a following sunrise occurs only in rural areas where S
W
219
-------
is generally quite small. Nevertheless, a constraint on the nighttime depth
of Layer 0 should be
nn/a» « < At (nighttime) (B-15)
0 WO *.
In the late afternoon the reverse condition arises and a drops abruptly
w
in a short period. Assuming a /a « 1 we get
W WO
AM = SAt (sunset) (B-16)
That is, the system gains mass at sunset equal to the emissions of the
sources in 1 time step.
During the nighttime hours when a is increasing slowly,
h
AM ~ (T£) SAt
nl
We conclude that the steady-state approximation that we used in formulating
the Layer 0 equation is acceptable provided that
Vhi" l
and that during nighttime hours
V'wo ! At
where a is the rms vertical velocity fluctuation on surface H and At is the
model time step.
220
-------
Appendix C. Explicit Forms of the b and g Matrix Elements b , g
That Enter in the y Equation. nm "
The Y equation (9-17), that describes the vertical material fluxes,
source functions, surface deposition and the like contains the matrix b
and the vector g, both of which are defined in Equation (9-13). Here we
prescribe the exact forms of these two parameters.
Rewriting Equation (9-13) here for reference we have
^ [Fn-l,n - Fn,n] + ancn ' Sn = bnlcl + bn2c2 + bn3G3 " ^n (C'^
where h is the thickness of layer n, c is the average concentration in that
layer, Sn is the emission rate, and an = 3£nV /3t.
Consider, for example, the Layer 1 equations. From Section 10 we have
F_ . - F. , = (arn - aT1)BC, + (1 - vcQn tQ^i
U = 4 N0 °3 3 l (C-4)
I 0 , otherwise
and
w X
QX = w X +"(1 - OB x = °3' N0> N02 or x (C'5)
221
-------
Let
ew ' <"TO ' CTTI>SX - (1 - an)(wi + V (c'6)
ex = w_x_ - W+X+(1 - e)(l - a)Qx (C-7)
Then on substituting (C-2, 3 and 4) into (C-l) and noting that
al = Mhl
and that for ozone the emission rate S in any layer is identically zero,
we get
Ozone (03)
1 °3 °3 C1 - U0)^ - a)v2£Qn
bll ' h7 ^ew3 * fii + (1 - aTo)(efl3 -- ^— - 93)]
1 U3 + v
b!2 - 4 d - -Tl^lm
(C-8)
b13 = o
-(1 - U0)(l - aTo)(l - a)vSNO
91 hl^Q3 + vj
In a similar manner we find that for species NO
Nitric Oxide (NO)
" °
b!2 ' h
- 0
.HO/UO
222
-------
Nitrogen Dioxide (NCL)
Species y (excluding 0, NO and N0)
ATT Species
b21^
-
b!2 - KJ (1 - "Tl>»im
((MO)
b!3 = 0
N02 (1 - a,
g, = S, 2 + __!
+ f9
b23 ' - hpe (1 - ac)
92 =S2
A ' ac
b!2s
h n (C-U)
b!3 = 0
cX (1 - «TTo)(l - a) .
9i •s^+ H /J;«i— vs..
l-ac
223
-------
Ozone (03)
1 - U
32
(C-13)
33
93 - h C(v
u
3U3C»
Nitric Oxide (NO)
NO r gv
— v + 6
NO
1-5]
(C-14)
33
_! TH u rNO y""'"o /"NO
93 - h, L«3 3C- ' v + 6NO< ;
Nitrogen Dioxide (N02)
NO,
'31 h, L v + 6
c T^T +1 - 5]
NO
^ (i - *)«
J, ( - M * H U )
h3 3
(C-15)
224
-------
NO ' SNO~ SNf)
-r1 - *.<°3>i * (1 - V -r
Species x (excluding 03, NO and
"31
532 = TT
0
b33 = h
1 ^
93 = H7 A3U3C- "
>j
where
ri , if H3 > 0
0 , if H3 <_0
- w(
7c
« -
M = ac[ ^ . gc + fe/A9] (C-18)
H3 = dH3/dt (C
We can check the correctness of the coefficients t>nm that we have just
formulated by determining whether they satisfy certain integral constraints.
In particular, in the absence of emissions (S = Sn = 0), chemical reactions
(R = 0), and surface deposition (e = 0), the total mass of pollutant in a
vertical column extending through all 3% layers must be constant. That is,
225
-------
h2Y2 + h3Y3] = 0
(C-20)
or
(C-21)
Since this condition must hold for any combination of values of YI» Y2 and
Y3, we find from (9-14) that (C-21) holds in general only if the following
three conditions are satisfied:
hlbll + h2b21 + H3b31 - fii = °
72
= 0
(C-22)
'1U13
'23
= 0
Consider for example the elements b in the equations governing ozone
(03). From (C-8, 12, and 13) we obtain with the aid of (5-13)
hlbll * h2b21 + h3b31
J-3
(C-23)
On substituting the expressions for a (5-48), M (C-18), and Qn (C-5) we
U3
find that the righthand side of (C-23) is identically zero, and hence the
first condition in .(C-22) is satisfied. In like manner we find that the
expressions given above for b for all species satisfy the consistency
conditions (C-22).
226
-------
TECHNICAL REPC "DATA
(Please read Instructions on the ret:-ne before completing)
1. REPORT NO.
2.
3. RECIPIENT'S ACCESSION NO.
4 TITLE AND SUBTITLE
A REGIONAL SCALE (1000 km) MODEL OF PHOTOCHEMICAL AIR
POLLUTION
1. Theoretical Formulation
5. REPORT DATE
6. PERFORMING ORGANIZATION CODE
7 AUTHOR(S)
Robert G. Lamb
8. PERFORMING ORGANIZATION REPORT NO.
9 PERFORMING ORGANIZATION NAME AND ADDRESS
10. PROGRAM ELEMENT NO.
CDWA1A/02-1335 (FY-82)
Same as 12
11. CONTRACT/GRANT NO.
;. SPONSORING AGENCY NAME AND ADDRESS
Environmental Sciences Research Laboratory - RTF, NC
Office of Research and Development
Environmental Protection Agency
Research Triangle Park, North Carolina 27711
13. TYPE OF REPORT AND PERIOD COVERED
In-house
14. SPONSORING AGENCY CODE
EPA/600/09
1 >. SUPPLEMENTARY NOTES
1;}. ABSTRACT
A theoretical framework for a multi-day 1000 km scale simulation model of photochemical
oxidant is developed. It is structured in a highly modular form so that eventually
•:he model can be applied through straightforward modifications to simulations of
^articulates, visibility and acid rain.- The model structure is based on phenomeno-
iogical concepts and consists of 3 1/2 layers. The interface surfaces that separate
che layers are functions of both space and time that respond to variations in the
meteorological phenomena that each layer is intended to treat. Among the physical and
chemical processes that the model is designed to handle are: horizontal transport;
photochemistry and nighttime chemistry of the products and precursors of pollutant
reactions; nighttime wind shear, stability stratification and turbulence "episodes"
associated with the nocturnal jet; cumulus cloud effects - venting pollutants from
the mixed layer, perturbing photochemical reaction rates, etc; mesoscale vertical
motion induced by terrain and horizontal divergence of the large scale flow; subgrid
scale Chemistry processes — resulting from emissions from sources smaller than the
model s grid can resolve; natural sources of hydrocarbons, NOX and stratospheric ozone;
and others. Considerable attention is given to the question of the predictability of
pollutant concentrations at long range and to the related problem of parameterization
of mesoscale diffusion, the design of model "validation" experiments, and the like
7. KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
1
18. DISTRIBUTION STATEMENT
RELEASE TO PUBLIC
b.lDENTIFIERS/OPEN ENDED TERMS
19. SECURITY CLASS (This Report 1
UNCLASSIFIED
20. SECURITY CLASS (This page)
UNCLASSIFIED
c. COSATI Field/Group
21. NO. OF PAGES
239
22. PRICE
:PA Perm 2220-1 (R«». 4-77) PREVIOUS EDITION is OBSOLETE
227
------- |