-------
The work upon which this publication is based was performed
pursuant to Contract No. EHSD 71-22 with the National Air
Pollution Control Administration, Environmental Health Service,
Public Health Service, Department of Health, Education and
Welfare.
-------
CONTENTS
SECTION PAGE
I INTRODUCTION 1
A. Perspectives on Modeling Needs 1
B. Overview of Previous Contributions 1
C. Advances Made in the Present Work 5
II KINETIC IMPROVEMENTS IN THE PHOTOCHEMICAL DESCRIPTION 7
A. Additional Oxidation Chains 7
B. Influence of Initial Composition 14
III MATHEMATICAL AND METEOROLOGICAL IMPROVEMENTS 21
A. Application of Fade" Approximants to Smog
Chamber Calculations . 22
B. Application of Pade1 Approximants to the
Atmospheric Model Equation 26
C. Meteorological Reformulation of Model 30
IV TESTS OF THE MODEL AGAINST LOS ANGELES BASIN DATA 38
A. Objectives of the Tests 38
B. Review of the Previous Validation Study 38
C. Sensitivity Studies on 1969 Trajectories
with the Expanded Model 40
D. History Analysis of the 1030 El Monte
Trajectory 45
V CONCLUDING REMARKS 51
REFERENCES 57
-------
ii
-------
ILLUSTRATIONS
NO.
1
2
3
4
5
6
7
8
Schematic of Diffusion Model for Air Pollution Simulation
Expanded Kinetic Mechanism
Computed Curves Compared With Experimental Points
Hydrocarbon Decay With Two Reactivities
Ozone Buildup With Two Reactivities
Influence of NO :HC-Ratio on Oxidant Production
X
Estimated Trajectory of Air Mass Arriving at El Monte,
1030 hours, 29 Sept. 1969
Estimated Trajectory of Air Mass Arriving at El Monte,
1130 hours, 29 Sept. 1969
PAGE
4
9
12
15
16
17
33
34
9 Estimated Trajectory of Air Mass Arriving at El Monte,
1230 hours, 29 Sept. 1969 35
10 Estimated Trajectory of Air Mass Arriving at El Monte,
1330 hours, 29 Sept. 1969 36
11 High Reactivity Hydrocarbon Concentration at Huntington
Park on Type 2 Days, 1968. f = 1/4; r = 1/2 39
12 Ozone Concentration at Huntington Park on Type 2 Days,
1968. f = 1/4; r = 1/2 41
13 Reactive Hydrocarbon History Along the 1030 Trajectory.
f = 1/4; r = 1/2 (Ground level concentrations) 46
14 Nitric Oxide History Along the 1030 Trajectory, f = 1/4;
r = 1/2 (Ground level concentrations) 47
15 Nitrogen Dioxide History Along the 1030 Trajectory.
f = 1/4; r ° 1/2 (Ground level concentrations) 47
iii
-------
ILLUSTRATIONS (Cont.)
NO. PAGE
16 Oxides of Nitrogen (NO + N02> History Along the 1030
Trajectory. f = 1/4; r = 1/2 (Ground level concentrations) 48
17 Ozone History Along the 1030 Trajectory. f = 1/4; r = 1/2
(Ground level concentrations) 50
iv
-------
I. INTRODUCTION
A. PERSPECTIVES ON MODELING NEEDS
Control of emission sources and of future land use offer possibili-
ties for managing air quality. Under statutory mandate, regional imple-
mentation plans must be submitted to show how these controls will be
applied. Mathematical models that give estimates of air quality in terms
of emissions in a particular airshed should play a key role in formulating
the plans.
The work described in this report has resulted in a diffusion model
with photochemistry that is applicable to pollutants that undergo chemical
changes. Although our further development of an earlier model is still
confined to the Los Angeles basin data, the method is applicable to the
broad class of problems. Most of the pollutants covered by criteria
documents participate in some forms of chemical reaction on time scales
comparable with air residence times in urban areas. To account for this
we have retained both the meteorology effects of dispersion and advection
coupled with photochemical effects of primary pollutants producing oxidant
compounds.
B. OVERVIEW OF PREVIOUS CONTRIBUTIONS
234
Since literature reviews ' ' of air pollution modeling in general
are available, our remarks here are confined to the previous work that
was specifically directed toward reactive pollutants. The discussion
centers on the relationship between the details emphasized by the model
and the levels of precision in either input or output data. At the very
least, any model must include the physics and chemistry of the governing
phenomena, but complexities which far exceed reasonable expectations of
specific validation are probably unjustified. In abstract language this
says that we can fit any observed data by using enough adjustable param-
eters in the model, but that we leave the realm of physics if we do not
subject most of the parameters to independent experimental verification.
-------
The models described by Calvert and Wayne, et al., stress chemical
kinetic mechanism. In the former, the photochemical smog mechanism is
reduced to 17 steps and in the latter to 31 steps. In the larger mecha-
nism care is taken to describe plausible reaction for the C-H../NO /air
system, but it is ultimately applied to the atmospheric case where dozens
of smog-forming hydrocarbons are involved. To preserve a consistent level
of detail, the number of reactions must be increased by a large factor;
however, this is avoided by adjusting the rate constants to obtain observed
results. Homogeneously stirred air parcels following wind trajectories
are assumed in order to retain first-order, ordinary differential equations.
Another way of restricting the mathematical complexity to first-
order, ordinary differential equations is given in Friedlander and
Seinfeld's paper. This method reduces the chemical system to seven
functional steps to be consistent with the level of observational detail
A
usually available for polluted atmospheres. Unlike the models discussed
above, this one accounts for finite rate atmospheric diffusion. Retention
of the mathematical simplicity is bought at the price of imposing the
Lagrangian similarity hypothesis. This means that the vertical profiles
of concentration in the polluted layer belong to simple parametric fami-
lies. Careful choice of the profile form functions may well allow this
method to be a fast and powerful computational tool.
8 9
Three-dimensional, time-dependent methods ' have been proposed
recently, but results for reactive atmospheres have not been reported at
this writing. Simplified chemistry must be employed in each of these
approaches because of the emphasis on details of advection and diffusion.
The body of meteorological data for most air basins falls far short of
the input requirements for any transport formulation of this complexity.
Here again it is difficult to avoid the problem of allowing too many
*
Even if more detail were included, there are no measured data on rate
constants for most of the additional reactions. Consequently, the model
is provided with additional degrees of freedom "wherever the kinetic
data are open to estimation".
-------
unspecified parameters to obscure the physically based portions of the
calculation. One of the new methods employs an Eulerian coordinate frame
which introduces the chronic problem of artificial (or numerical)
diffusion. Small time and space steps are the only means suggested
heretofore to remedy this problem. The other method uses a discretized
representation of species concentrations in a Lagrangian frame. For
multicomponent systems, the latter type of approach has previously been
limited because of the large quantity of rapid access memory storage
needed for each species. Both of these problems have been recognized
for some time with three-dimensional, time-dependent calculations. More
powerful computers and novel algorithms may well hold their solutions.
The method developed at General Research Corporation is more modest
in its chemical detail than the models first described above, and it is
far simpler in its transport formulation than the three-dimensional,
time-dependent approaches. An attempt has been made to keep the mathe-
matical descriptions near the level of detail and precision of observa-
tional data. This has resulted in a lumped-parameter, chemical kinetic
mechanism and a Lagrangian air parcel with finite-rate vertical diffusion.
The chief advantages of this over the approach of Ref. 7 are complete
flexibility in specifying vertical diffusion and in accounting for
distributed emission sources.
Our approach avoids the artificial numerical diffusion problem by
the choice of natural (or so-called "intrinsic") coordinates. We align
the coordinate system with the natural airflow at each instant so that
there is practically no advection across any grid line. This permits us
to follow air masses, to include vertical dispersion, and to account
for chemical changes simultaneously. Figure 1 schematically illustrates
these points. The large memory requirements do not appear in our approach
because we trace the history point-by-point and output the results sequen-
tially. Balancing the needs for realistic chemical and meteorological
descriptions against efficiency of calculation, the GRC model has achieved
-------
SUNLIGHT IS GIVEN AS
A FUNCTION OF TIME
00
c^
-------
It is useful to digress here to examine the relative magnitude of
the horizontal diffusion terms which are neglected. Carbon monoxide will
be chosen because it undergoes negligible atmospheric reaction over time
scales of interest. We shall form a ratio taken from the CO-balance on
the moving air parcel. The ratio of CO-loss horizontally to CO-gain
vertically gives a measure of how serious our assumption is. The air
2
parcel is taken to be 9 km in horizontal area because this approximates
the spatial resolution of the source grid. The vertical dimension is
taken to be 100 meters, typical of a high pollution day. In order to
force a high value for the ratio, the time of 1000 hours PST is selected
because the numerator of the ratio is high (the air having just been
loaded up by morning traffic) and the denominator is low (peak traffic
has passed). Horizontal loss is the turbulent diffusion coefficient
(-V300 m /min) times the average crosswind gradient (^3.8 x 10 mass
f 9
fractions per meter) times the area of all four sides (1.2 x 10 m ).
The horizontal gradient is maximized by using a concentration decrement
equal to the highest 10 A.M. reading in Huntington Park in 1969, which
was 19 ppm. The characteristic distance over which the concentration
drops from 19 ppm to zero is taken as the order of 30 miles. This gives
3
a loss of about 0.1 m /min. Using a grid square in the central basin
(containing both freeway and surface street contributions), we get a
3 2
gain of about 20 m /min averaged over the 9-km area. Consequently, the
error committed is of the order of 0.5% in the CO-mass balance. We
consider this to be a high estimate for the reasons outlined above.
C. ADVANCES MADE IN THE PRESENT WORK
Following the philosophy expressed above, we have advanced our
techniques to the stage where they can be adapted to planning and imple-
mentation of air pollution abatement measures. As described in the next
section, the photochemical description has been extended to account for
varying ratios of oxides of nitrogen to hydrocarbon. Some hypotheses of
carbon monoxide interactions have been tested. These kinetic improvements
add several reaction steps to our previously adopted mechanism. With no
-------
changes in numerical methods, the costs in computing time would have
grown roughly in direct proportion.
Numerical integration improvements have more than offset the
increases in chemical complexity. Another section describes these along
with refinements in the meteorological realism of the model. The moving
element shown in Fig. 1 replaces our previous formulation which was
limited either to a horizontally uniform slab or a two-dimensional
Eulerian coordinate system. The increased computing efficiency applies
to both the homogeneous (smog chamber) and the diffusion (atmospheric)
runs because we generalized it to apply to partial differential equations,
In the final sections of the paper, validation tests are presented
for Los Angeles basin data and the significance of the results are
discussed in terms of future applications of the technique.
-------
II. KINETIC IMPROVEMENTS IN THE PHOTOCHEMICAL DESCRIPTION
A. ADDITIONAL OXIDATION CHAINS
It has been known for some time that hydrocarbon consumption rates
observed in chamber experiments cannot be explained by ozone and oxygen
atom attack alone. This is discussed in the review by Altshuller and
Bufalini where it is noted that Schuck and Doyle termed the disparity
an "excess rate." Our previous modeling work treated this by increasing
both the 0_ and 0-atom rates of reaction with hydrocarbon.
To use known rate constants to a maximum extent in the present work,
we have added a hydrocarbon oxidation chain to reflect the attack of RO
12 13
upon the hydrocarbon. Because of its suspected dominance, ' the
hydroxyl radical (OH) was assumed to be the only RO reacting with the
hydrocarbon.
Following the inorganic cycle formed by
hv + NO,, -»• NO + 0
d>*
NO + 03 -> N02 + 02 (2)
we have allowed each hydrocarbon oxidation to generate multiple R0_
radicals as expressed by the b-factors in the steps
0 + HC -v b1(R02) (3)
The two reaction steps bracketed by (1) are combined by assuming 0-atom
stationarity internal to the logic of the computer program. Hence
reaction (1) nets the system hv + N02 -»• NO + 0».
-------
OH + HC -* b2(R02) (4)
03 + HC -»• b3(R02) (5)
Subsequent conversion of NO to NO,, was assumed to occur via
R02 + NO •* N02 + d(OH) (6)
where d expresses that fraction of the conversion responsible for
returning hydroxyl to the system.
Note that reactions (4) and (6) form something of a closed loop
and that a stability requirement like b_d < 1 may be needed to prevent
R02~runaway. This is not a precise requirement because of the variety
of external factors that influence the flow rate of free radicals around
the loop. Some of the R0_ is formed by the other two reactions, and
there are radical removal processes in the chain termination steps we
discuss in the following paragraph.
Both R0_ and OH are held in check by removal steps that
terminate the chains. We continue to use some lumped parameter reactions
such as
R02 + N02 -> c(PAN) (7)
and some elementary reactions like
OH + NO + M -> HONO + M (8)
OH + N02 + M -»• HN03 + M ' (9)
This rounds out the mechanism we have called Mark XII. It is diagrammed
on Fig. 2.
PAN is used to denote peroxyacetylnitrate.
-------
he
i
i ^
NO * 1 *
1
-
*
1
j
^-t— , '
HC ' 1— ^
1
1 J
'
>2+M
;
, i
J
3
r
^
r~
i
i ^~~
NO2
|
^U
V
V
P AKI
^
y
(C
^
,
l
! I
^ i
i i
• HI
r
[
R(
(HC
|
3
V
t
i
1
|
2
,)
1
k
Figure 2. Expanded Kinetic Mechanism
-------
In summary, the main changes from the previous scheme are:
addition of radical species RO (treated as OH), the provision for a
multiplicity of hydrocarbons (not shown in equations), and the elimination
of ozone production from an RO- reaction. The multiplicity of hydro-
carbons is done by adding parallel reactions to (3), (4), and (5). The
omission of the ozone production reaction is justified on the basis of
14
its suspected endothermicity.
For our initial validation tests we tried the rate constants com-
14
piled and estimated by Westberg and Cohen. Stoichiometric coefficients
b.. , b~ , and d were determined by hand-calculation analysis of some of
the propylene data represented by Ref. 15. We found that bd = 1 where
b is a composite branching factor, and b = 2 in order to explain the
decay rate of NO compared with that of propylene.
Explosions of R02~concentration characterized some early choices
of bd-combinations because of the positive feedback loop involving radical
reactions shown in Fig. 2. Since these occurred in the first few seconds
of real time, they were difficult to detect. Our usual printing interval
of 1 minute of real time gave computer outputs showing sudden drops of
NO and HC (to some nonzero levels) during the first minute followed
by the dynamic equilibrium between NO, N0» , and 0_ with a slow decay
£• J
of hydrocarbon. The underlying pathology was discovered by rerunning
with special diagnostic outputs.
Nominal rate constants were adjusted to reconcile computed results
with measured values. For example, the (OH + HC)-rate was cut to about
one-third the estimated value and the (NO + R0»)-rate was increased
eightfold. Reaction rates which have been reported individually in the
literature were held at their nominal values including the (0 + HC)-rate
which we had to increase many-fold in previous work. Also, the reten-
tion of reaction (5) was still necessary to describe the continued
C0H,-decay after the near disappearance of NO.
J o
10
-------
In order to reproduce the late-time behavior of N0_ and propylene,
we included chain termination reactions in addition to those indicated in
Fig. 2. In their more recent review, Altshuller and Bufalini point out
that HONO might be formed by the reaction with water vapor
NO -4- N02 + H20 -> 2HONO (10)
which is likely to proceed in the two steps
NO + N02 -» N203 (11)
N2°3 + H2° ~" 2HONO
12 13
A possible source of OH-radical ' is the photodissociation of HONO
hv + HONO ^ OH + NO (13)
Figure 3 shows the computed concentration history using this scheme
to describe the experimental observations of Altshuller, et al., and
Table 1 summarizes the reaction rate constants that were used. The column
in Table 1 labeled "nominal" gives the values we started from on the basis
of tabulations in Ref. 14. Following an empirical procedure, our calcula-
tions effectively combined H-0-concentrations with the rate constant for
reaction (10) as indicated; therefore no comparable nominal value is shown.
Higher levels of precision than that shown in Fig. 3 were not among our
objectives, since a detailed investigation of propylene is not justified
for describing the potpourri of hydrocarbons in the urban atmosphere.
Orders of magnitude for conversion times and levels are sufficiently good
to proceed from here. The objective that we have fulfilled is to repre-
sent the so-called "excess rate" of hydrocarbon oxidation with added
reactions instead of artificially increasing the 0-atom rate with
hydrocarbon.
11
-------
2.5 ppm -
2.0 -
Z
o
o
Z
o
A PROPYLENE
• NITRIC OXIDE
> ruirN is ur ,.
• NITROGEN DIOXIDE j ALTSHULLER, et ol "
T OZONE
60 80 100 120 140 160 min
TIME
Figure 3. Computed Curves Compared With Experimental Points
12
-------
TABLE 1
RATE COEFFICIENTS FOR MODEL
HYDROCARBON/NITRIC OXIDE MECHANISM
(Stoichiometry imbalances may occur
because of lumped parameter assumptions.)
Reaction
Model Values
From Validation
Nominal Values for
Propylene System-^
hv + NO- -> NO + 0
0.4 min"1
1.32 x 10"5
40 ppm min
6100 ppm min
80 ppm min
R02 + NO -»• N02 + 0.5 OH 1500 ppm'-'-min"1
_i _i
03 + NO
0 + HC
OH + HC ^ 2RO,
2
R02 + N02 •> PAN
OH + NO -»• HN02*
OH + N02 •*• HN03*
0_ + HC -»• RO^
NO
2HNO,
**
, -1 . -1
6 ppm min
10 ppm min
30 ppm min
0.0125 ppm" min"
0.01 ppm min
hv + HNO ->- NO + OH 0.001 min"1
0.4 min
"1
"5
1.32 x 10
22-44 ppm min
6100 ppm min
244 ppm min
122 ppm min
122 ppm min
99 ppm min
300 ppm min
0.00927 - 0.0125
Rate constant lumps third body concentration
**
Water vapor lumped into rate coefficient
13
-------
B. INFLUENCE OF INITIAL COMPOSITION
1. Synergism Between Hydrocarbons
The fact that free radicals produced in the photooxidation of one
hydrocarbon can accelerate the attack on a coexistent hydrocarbon has
been illustrated with experimental findings by Altshuller and Bufalini.
Because of the limited understanding of detailed interactive mechanisms,
we limited our multiple hydrocarbon investigation to a simple parametric
study. Using a truncated version of the mechanism in Table 1, we added
a parallel series of oxidation steps. This scheme is shown in the flow
chart of Fig. 2 with dashed lines along the bottom portion of the diagram.
In the parametric study, the same initial values of total hydro-
carbon and nitric oxide as those in Fig. 3 were employed. The hydro-
carbon, however, was partitioned into two hypothetical compounds one
having triple the (OH + HC) rate constant and the other having one-third
the (OH + HC) rate constant as propylene. The compound with the tripled
rate is called species "B" and the one with the decreased rate is called
species "A". All other reactions in the scheme are the same in all
respects as their counterparts in the (C-H, 4- NO )-system.
JO X
Figure 4 shows the hydrocarbon decay curves from the parametric
study. On the same coordinates we plot curves for model calculations of
pure hydrocarbon "A", pure hydrocarbon "B", and pure propylene for com-
parative purposes with the half-and-half mixture of hypothetical compounds,
The most interesting aspect of Fig. 4 is that the mixture decays faster
than propylene even though it contains equal quantities of "A" and "B"
which bracket propylene symmetrically. The coupling between the two
hydrocarbons occurs via the R0_ conversion of NO to N0» and the
production of OH . Evidently the combination of the rapid incubation
(afforded by the reactive hydrocarbon) with the twofold chain-branching
is sufficient to increase the decay rate over that of pure propylene.
14
-------
SPECIES 'A1 HAS AN (OH + HCJ-RATE ONE-THIRD THAT OF C3H6
SPECIES 'B1 HAS AN (OH + HCJ-RATE THREE TIMES THAT OF C3H6
ppm
Z
o
t-
<
at
i-
Z
ui
Z
O
y
u
2.0
1.5
1.0
0.5
PURE
'A-
CVJ
I
I
I
I
I
I
20
40
60
80
TIME
100 120 140 160 min
Figure 4. Hydrocarbon Decay With Two Reactivities
Figure 5 shows corresponding curves for ozone buildups predicted
in the computer experiment. Again the combination seems to exceed the
pure propylene in its rapidity to produce ozone. This occurs because
of the enhanced conversion rate of nitric oxide to nitrogen dioxide
stimulated by the early abundance of RCL generated by the high-reactivity
fraction.
These findings illustrate a plausible mechanism for the interaction
of multiple hydrocarbons. In future applications of our photochemical/
diffusion model this type of representation is likely to be more realistic
than a super-detailed single hydrocarbon mechanism with many adjustable
constants. Also, other radical interactions like R09 + HC and RCO, + HC
16
might provide plausible coupling paths in addition to the hydroxyl
attack on the hydrocarbon. Indeed, other forms of coupling may favor the
slower reacting hydrocarbon in a similar numerical experiment with a
half-and-half mixture.
15
-------
HYDROCARBON 'A1 HAS AN (OH + HC)-RATE ONE-THIRD
THAT OF C3H6
HYDROCARBON 'B' HAS AN (OH + HCj-RATE THREE TIMES
THAT OF C3H6
ppm
I.Oi-
Z
o
0.5 -
«J
z
o
u
CO
O
0 20 40 60 80 100 120 140 160 min
TIME
Figure 5. Ozone Buildup With Two Reactivities
2. Modeling the Influence of NO /HC-ratio on Oxidant Production
Our main purpose in including the chain termination details for
nitrogen oxides in the previous subsection was to improve the modeling
of initial mixture effects. It will be recalled that our earlier report
issued a caveat restricting the seven-step mechanism to a narrow range of
NO /EC ratios. Later, Seinfeld pointed out that the slope of the
curve of peak oxidant versus NO /HC-ratio did not even have the correct
sign. If modeling is to be used to evaluate hypothetical abatement
strategies, it is clear that this deficiency must be remedied.
*
Figure 6 shows the results from the twelve-step mechanism given
in Table 1. Experimental results for pure and mixed hydrocarbons are
shown on this plot of peak oxidant versus NO /HC-ratio. The model
X
mechanism displays a monotonically declining trend of oxidant production
as the oxides of nitrogen increase. Both of the pure hydrocarbons
Stationary state treatment of 0-atoms reduces this to eleven.
16
-------
z
q
t—
<
ec
i—
Z
o
u
»-
z
<
o
x
o
2 ppm C3H6
REF. 18
3 ppm n-C4H1Q
2 ppm HC
MODEL
RESULTS**
1 TO 2 ppm NMHC*
AUTO EXHAUST REF. 19
LOS ANGELES
INVENTORY
RATIO
AUTO EXHAUST\
2 ppm NMHC* \
REF. 20
CO
to
CXI
I
'j;
0.4 -
0.2 -
\
\
0 0.5 1.0 1.5
MOLE RATIO NOX TO HYDROCARBON
*NMHC:NON-METHANE HYDROCARBON
**THIS CURVE IS LOWERED CONSIDERABLY BY OUR
SUBSEQUENT ADOPTION OF ONE-HALF THE
PROPYLENE RATES FOR ATMOSPHERIC CALCULATIONS
Figure 6. Influence of NO :HC-Ratio on Oxidant Production
17
-------
measured by Altshuller, et al. , show peaks at widely separated values
of the mixture ratio. It seems plausible then that the typical mixture
of atmospheric pollutant hydrocarbons should show a decrease of oxidant
with increasing NO /HC-ratio. This is supported by the fact that the
relative abundance of hydrocarbon species decreases with increasing
reactivity. Now since the low reactivity species peak at lower
NO /HC-ratios the combined effect should give a negative slope.
X
The two plots on Fig. 6 of oxidant production from diluted and
19 20
irradiated automobile exhaust ' both have negative slopes comparable
to that of the model. As stated earlier, our purpose in using chamber
data for propylene to validate the mechanism is to capture the main
features of the experiments that apply to the atmosphere. If, on the
other hand, propylene were the subject of our study, we would proceed
to expand on the mechanism to get the proper curve shape shown in Fig. 6.
Having subjected the mechanism to these laboratory tests, we use
it in its present form for atmospheric validation studies described in
a later section.
3. Examination of Carbon Monoxide Effects
13 21 22
Recent work in photochemistry ' ' of smog has suggested that
carbon monoxide may play a role in accelerating the photooxidation of
hydrocarbon/nitric oxide mixtures. The mechanism suggested is
OH + CO ->• CO. + H
H + 02 + M -> H02 + M (15)
H02 + NO ->• OH + N02 (16)
(17)
18
-------
Reaction (15) is postulated to contribute H02 rapidly as it responds
to the presence of H-atom, and reaction (16) enhances the NO-to-NO.
conversion already occurring in reaction (6). Likewise, the hydroxyl
radical supplied by reaction (16) goes on to react with the hydrocarbon
in reaction (4), producing R02~radicals in a branched chain already
described. Reaction (17) is a chain termination step suggested for H02.
The key feature we wish to examine in this alternate path is the
competition of CO with hydrocarbon for the hydroxyl radicals. The
rate constant for OH + C-H, suggested in Ref. 14 is approximately equal
to that suggested in Ref. 21 for OH + CO . We tried a wide range of
(OH + CO)-rates holding the basic system at the values in Table 1. Since
reaction (15) is fast, reactions (14) and (15) were added to give the
overall reaction
02 + OH + CO -»• C02 + H02 (18)
The rate constant for reaction (16) was taken to be the same as its
RO..-analogue, reaction (6), and the rate constant for reaction (17) was
-1 -1
assigned the value of 1.2 ppm min based on the estimate in Ref. 14.
The species 02 is carried implicitly in the calculation, and the rate
constant for reaction (18) is varied parametrically for this study.
Beginning with (CO + OH)-values at the level suggested in Ref. 21,
(k..£ - 200 ppm min ) we found very rapid NO-to-N09 conversion and HC
disappearance with 100 ppm added CO . Only when k,0 was lowered to
-1 -1
1 or 2 ppm min did we obtain only small effects. Since much experi-
mental verification is yet to be done, only a limited effort was devoted
to this aspect of -the study. Table 2 summarizes the main findings. It
indicates that the half disappearance time of NO is the strongest
effect of k-R on the system, and the half-disappearance time of hydro-
carbon is least affected. Both the peaking time for N02 and the half-
rise time of ozone are moderately sensitive to the rate constant.
19
-------
TABLE 2
INFLUENCE OF 100 ppm CARBON MONOXIDE ON THE SYSTEM
SHOWN IN FIG. 3 WITH VARIOUS VALUES OF THE RATIO
OF (OH + HC)/(OH + CO) REACTION RATE COEFFICIENTS
k4/k!8 tl/2 ~ N0(mi
•k
0.4 <1
4 6
40 25
80 29
System in
Fig. 3 36
(No CO)
n) t , - N02(min) t.
4
13
43
52
66
L/2 ' HC(m±n)
40
48
73
81
92
t1/2 - 03(min;
10
22
62
73
83
This is with k..g set at the estimated value from Ref. 21.
In any event, whether we use the estimated rate of (OH + CO)-
reaction or even only preserve its ratio to that of the (OH + HC)-reaction
12
suggested in Ref. 14, we get very large effects. Other work suggests
a much larger (OH + HC)-rate constant for propylene. In view of our
findings and of the wide disagreements on some of the rates, a very real
need exists for definitive experimentation to achieve a deeper under-
standing of the kinetics.
20
-------
III. MATHEMATICAL AND METEOROLOGICAL IMPROVEMENTS
Our previous study brought out some of the limitations which hamper
the original time-dependent atmospheric model. This model is based on a
23
slab geometry and is solved by a Crank-Nicolson technique. The omission
of advection limited the calculation to places and times characterized
by stagnant conditions (hence our choice of the October 1968 episode for
Huntington Park validation data). Admittedly, this is a case that the
plume models (described in Refs. 2, 3, and 4) are incapable of handling,
but applications demand that cases with wind transport be treated. Our
original response to this need took the form of an Eulerian governing
equation for species concentration:
It
- £ - * (» If ) + •
where c = concentration
t = time
x = downwind distance
u = wind speed
z = height above ground
D = diffusion coefficient
R = chemical reaction rate
Solutions of Eq. 19 are costly because of the need to suppress artificial
*
diffusion by refining mesh and step size. Smog chamber analyses use
just the first and last terms in Eq. 19 so that they depend on ordinary
differential equations. These are solutions which describe the time-
dependent behavior of a homogeneous gas mixture. We used standard Runge-
Kutta techniques to solve them in the original code.
Artificial diffusion is an error from numerical differencing of the
equation. It produces effects that resemble physical diffusion in
the x-direction.
21
-------
The combined needs for computing efficiency and meteorological
realism have been met by improvements described in this section. The
first two subsections describe the incorporation of a classical approxi-
mation method into our numerical integration. The implementation of
these computing techniques strengthen the study by allowing numerous
modeling trials on both laboratory and atmospheric cases. The final
subsection outlines the change from Eulerian to semi-Lagrangian coordi-
nates in the meteorological formulation. By following air masses we use
a natural coordinate system thereby eliminating the artificial diffusion
errors.
A. APPLICATION OF FADE APPROXIMANTS TO SMOG CHAMBER CALCULATIONS
The chemical model previously described is implemented in rate
equations which must be integrated numerically. Because of the widely
varying reaction rates, the equations possess a characteristic known as
"stiffness" which makes their integration difficult. In the interest of
economy and accuracy, special techniques must be employed in obtaining a
24 25
solution. Below we describe a method which utilizes Pade approximants. '
For these ordinary differential equations, the procedure has yielded
significant speed improvements over the Runge-Kutta method previously
used.
For this class of problem, the implicit approach is the basic
advantage of the new method, thus guaranteeing stability. Consequently
accuracy controls step size. This provides the principal gain in
computing time.
The rudiments of the integration method are twofold. One involves
a linearization of the nonlinear chemical term at every step; the other
an approximation of the exponential term that appears as a result of the
linearization operation. Thus, consider the following system of equations
that describe the changes in species concentration due to chemical reaction:
22
-------
(20)
where c = vector of species concentrations
R(c) = vector of chemical reaction rate functions
c =
V
C2
•
cs
R(c) =
\M~
R2(c)
•
R8(c)
s = number of chemical species
To linearize the rate vector R(c) , we expand it about c in a Taylor
series. The subscript "o" denotes evaluation at some time t = t . Thus
R(c) = R(CQ) + J(c - CQ) + 0(Ac)'
(21)
where
-J — 1 O Q
. I J — _L J ^ J • • * > a
is the matrix of order (s * s) of first derivatives (the Jacobian
matrix) of the rate function. Substituting Eq. 21 into Eq. 20 yields
dc
dt
= Jc + B
(22)
where B = R(c ) - Jc . It should be noted that at time t = t the
o o o
quantities c , R(c ) , and J are known. Hence, the matrix differential
23
-------
equation (22) is a linear equation with constant coefficients and a
constant forcing function. Its formal solution is, therefore,
c(t) = eJAt [CQ + J~VJ - J^B (23)
where J denotes the inverse (or reciprocal) of the matrix J.
An approximate solution of Eq. 22 is obtained from Eq. 23 by
suitably approximating the matrix exponential e . This is accom-
plished by means of the Pad£ approximants of the exponential function.
These Fade approximants are rational functions of the form
JAt = P(JAt) o -1
6 Q(JAt) ^ V
where P(JAt) and Q(JAt) are matrix polynomials of degree p and q ,
respectively. They are defined by
(p + q - k)!p! xk f .
(25a)
(p+q)!k!(p -k)!
k=0
q
n - (P + Q ~ k)!q! f ,k f ,
Q ~ p + q)!k!(q - k) ! ("JAt) (25b)
k=0
Thus, substituting Eq. 24 into Eq. 23 we obtain
c(t) = Qc^ + JB - JB (26)
The choices p = q = 1 and p = q = 2 are especially useful. The
former yields a third-order integration formula, the latter a fifth-order
formula. As an example, for p = q = 1 we obtain from Eq. 25:
P = I + 1/2 JAt and Q = I - 1/2 JAt (27)
24
-------
where I is the identity matrix. Substituting Eq. 27 into Eq. 26 yields,
after the appropriate manipulations,
c(t) = c + [I - 1/2 JAt]~1AtR(c ) (28)
o o
Equation 28 has a form that is especially well-suited for digital computa-
tion. Note that the process requires at every time step: (1) the evalua-
tion of R(c) and J , and (2) the inversion of the matrix Q . The
latter is by far the more computationally costly of the two requirements.
In our case, the number of species determines the size of the matrix to
be inverted. Since this number is not large in our model, the matrix
inversion poses no special problems. Fortunately the number of species
is usually smaller than the number of reactions.
Table 3 shows examples of the gains obtained using the new computa-
tional scheme. The typical real-time/computer-time ratio was increased
from 36/1 to 180/1. Perhaps more significant is the fact that the Fade"
method has allowed us to obtain acceptable solutions in situations where
Runge-Kutta either failed to converge or produced spurious solutions.
One such instance is the integration of full differential equations for
the free-radical species: the Fade" method successfully computed solutions
without algebraic substitution of stationary-state assumptions whereas
Runge-Kutta failed to produce any solution.
As mentioned previously, the step size must still be controlled
in order to preserve accuracy. Referring to Eq. 21 we can see that the
difference Ac = c - c determines the order of the truncation error of
o
the linearization operation. The step size must be constrained so that
Ac will remain within some bound specified by the user. Our codes con-
tain a variable-step-size feature which ensures that Ac will stay within
the prescribed bounds.
25
-------
TABLE 3
CDC 6400 CENTRAL PROCESSOR TIMES TO RUN
C.H, + NO + hv CHAMBER EXPERIMENT SIMULATION
3 b
Type of Run
Runge-Kutta
Fade Approximant
Rapidly changing free radical
concentration
Typical run for 150 minutes
real time
Run which diverged for
Runge Kutta
2400 seconds
250 secon'ds
267 seconds
/diverged at\
180 minutes,
\real time /
122 seconds
50 seconds
50 seconds
/ran to completion,\
I 150 minutes real 1
\time /
It should also be noted that the nature of our chemical model is
such that it leads to a very good first-order approximation as indicated
in Eq. 21. This occurs because the full expansion of our particular rate
function, R(c) , contains only three terms. Hence in the first-order
approximation only one term has been dropped from the expansion.
B. APPLICATION OF FADE APPROXIMANTS TO THE ATMOSPHERIC MODEL EQUATION
We now describe the use of Pade techniques to solve the diffusion
equation with chemical reactions. The equation in question is referred
to a moving air mass. Use of the Lagrangian approach eliminates the
second term in Eq. 19 which is reduced to
(29)
where D is the vertical diffusion coefficient and c and R(c) have
been defined previously. Applying the "method of lines" to Eq. 29 the
partial differential equation is transformed into an ordinary differential
equation. Each line represents a different elevation and the resulting
ordinary differential equation for a line describes the concentration
26
-------
changes for a fluid mass moving along the line. The equation for the
i_th line is derived from Eq. 29 to give
dt
- 2c
(Az)'
where
R. = R(c. ) , i
1, 2,
,M,
(30)
where M is the number of mesh
points from the ground up to the mixing depth. The symbol c. defines
a vector of s species at the ith point in space. Similarly, R.
defines a vector of s rate functions at the ^th mesh point. Thus,
"R,
ci=
2i
8i
21
R
.
si
In order to solve Eq. 30 by the Fade" method, R. is linearized
as in Eq. 21 to yield
Ri " Rio + Ji(ci -
(31)
where the zero subscript again denotes evaluation at time t = t , and
J. is the Jacobian matrix of the rate function at the ith spatial mesh
point. Substituting Eq. 31 into Eq. 30 yields
dt
Xtci-l-2ci+ci+l]
+R
io
(32)
where X = D/(Az) . Writing Eq. 32 more compactly we obtain, after
collecting terms,
27
-------
= (XA + J)C + R - JC
dt oo
where H=XA+J,B=R -JC , and
o o
A =
-21 I 0
I -21 I
0 I -21 I
J =
R =
o
0
20
(33)
0
0
0
I -21
RMO
Equation 33 now has the same form as Eq. 22 and its formal solution is
thus analogous to Eq. 23. Hence
C(t) =
[co + H-V] - if
(34)
Approximating the exponential by
HAt _ I + 1/2 HAt £ -1
6 I - 1/2 HAt Q ~ ^
28
-------
and substituting in Eq. 34 and simplifying we obtain
where
C = C + QG At (35)
o o
G = XAC + R
o o o
It should be noted that any boundary conditions are implicitly contained
in Eq. 35.
Equation 35, which is the analog of Eq. 28, may appear to be simple,
but the inversion of Q poses some problems arising from its dimensions.
H is an (M * M) matrix each one of whose entries is an (s x s) matrix.
Hence Q is an (sM x sM) matrix, and its storage as well as 'its inversion
may be difficult. For example, if we have 10 species (s = 10) and 10
space grid points (M = 10) then Q is (100 x 100) and requires rather
long computing times as well as large amounts of core storage. However,
we take advantage of the structure of Q to develop an algorithm that
reduces the problem to the inversion of M matrices of size (s * s) .
It is noted that since J is a diagonal matrix and A is tri-
diagonal, then H is tridiagonal. Moreover, the off-diagonal entries
of H are identity matrices. Forming Q we have
Q = I - 1/2 HAt
and it is apparent that Q is tridiagonal, its off-diagonal entries
being identity matrices also. Thus only the M diagonal entries of Q
2
need be stored, thereby reducing memory requirements to Ms cells
compared to the maximum of M s . Writing Eq. 35 in the form,
QC = QC + G At
o o
29
-------
we note that the tridiagonal property of Q allows us to solve for C
without explicitly inverting Q . This simplification is afforded by an
23
extended form of Gaussian elimination which requires the inversion of
M matrices of size (s x s). This is preferable to inverting the full
(sM x sM) matrix.
Note that we can also take advantage of the fact that the off-
diagonal elements of Q are identity matrices to reduce the amount of
computation.
Because series are truncated in these step approximations, some
remarks about the resultant errors are necessary. The order of the
3
truncation error of the solution is 0(At) . This occurs because the
approximation of the exponential is
HAt s I + 1/2 HAt
e I - 1/2 HAt
Also, as in the solution of the chemical rate equations, a linearization
2
error of order 0(AC) appears in the approximate solution. Thus, the
same precautions that were taken to preserve accuracy in solving the
ordinary chemical rate equations are required in the solution of the
atmospheric model equation.
The new scheme has yielded a better than fourfold increase in the
real-time/computer-time ratio, i.e., from 7/1 to 30/1.
C. METEOROLOGICAL REFORMULATION OF MODEL
Long computer runs and unacceptable error propagation hampered the
operation of the earlier time-dependent advection and kinetics (TADKIN)
code. Its application was limited to an analysis of carbon monoxide to
examine the influence of advection added to that of vertical diffusion.
At the inception of the present work, these severe limitations dictated
a reformulation of the advective description before any atmospheric runs
30
-------
were attempted. Combined with the mathematical advances just described,
the semi-Lagrangian embodiment of the diffusion with kinetics (DIFKIN)
code extends the application of photochemical diffusion modeling to a
great many more cases than previously thought practicable.
Following a center-of-mass fixed coordinate system tied to an air
mass, we make use of intrinsic coordinates to avoid artificial diffusion
in the horizontal direction. Physical diffusion, therefore, is distinct
and identifiable because the moving control volume can be allowed to
undergo mass exchange with a neighboring air mass in a prescribed fashion.
The question of horizontal spatial resolution is answered by a selection
of source grid size and the vertical resolution is set by the choice of
the interval size in the z-direction.
Mixing depth or vertical extent of the pollution layer can be
handled in one of two ways. In the first, the vertical diffusion coeffi-
cient is assigned spatially and temporally dependent values according to
the combined effects of shear and buoyancy as turbulent energy sources.
The vertical extent of the grid is chosen far above the conventional
mixing depth so that the pollutants are automatically confined vertically
to a degree which depends on the diffusion coefficient profile. In the
second way of specifying the vertical behavior, however, the vertical grid
either continuously adjusts to the input values of mixing depth or assumes
an average value for the inversion base altitude. We have favored the
use of this second approach because it limits the field of computation
to the polluted region and thereby conserves computer time. A zero-flux
boundary condition is imposed at the top of the net. Vertical diffusion
coefficients employed in the earlier work are still used since other reports
2fi ?7 9R
.in the literature appear to corroborate those values. ' '
The advection patterns are obtained by taking a weighted average
of wind-station readings surrounding the point in question. This is truly
a receptor-oriented modeling approach since the trajectories are constructed
31
-------
in hourly upwind segments from the measuring station. Reciprocal distance
weighting is employed because of the dominant plane flow pattern which is
A
governed by combinations of sources, sinks, or vortices. In most cases
the three nearest neighbor stations are included in the stepwise upwind
tracing or air movement, but occasionally conditions of proximity suggest
including only the nearest one or two stations where coastal crossings
are involved in the analysis.
Figures 7, 8, 9, and 10 show examples of the paths of air masses
estimated in the manner just described. The times of arrival in El Monte
are used to identify each trajectory in any subsequent references in this
report. The date of September 29, 1969, is chosen because of the rich
variety of data that we have available for that day. Note that the
earlier morning meandering patterns give way to the dominant onshore
flows for all trajectories after 0800 to 0900 hours Pacific Standard
Time. The meteorological formulation that we have adopted takes the time
and location information from these trajectories to establish the initial
conditions and the boundary conditions. Initial conditions are specified
as vertical profiles of concentration and boundary conditions as time
histories of surface-based pollutant emissions.
Using this approach with data from the two monitoring sites, we
are able to match the computing requirements to the data base. Little
justification exists at this time to apply complicated grid methods to
these validation studies because we only have measurements in great
detail from Commerce and El Monte. Just as overly-elaborate chemical
models cannot be adequately tested in the atmosphere, neither can highly
complex diffusional and atmospheric formulations be checked out with our
present store of field observations. Certainly, the day will come when
frequent high-resolution readings are made and some of the more detailed
theories can be tested.
See Ref. 1 for a discussion of characteristic dimensions.
32
-------
AZUSA
METROPOLITAN
LOS ANGELES
8 10 mi
Figure 7. Estimated Trajectory of Air Mass Arriving at El Monte,
1030 hrs, 29 Sept. 1969
33
-------
1030^0530
T n~inj|flrii|ifhiiiii I iu
COMMERCE! I *
METROPOLITAN
LOS ANGELES
024 6 8 10 mi
Figure 8. Estimated Trajectory of Air Mass Arriving at El Monte,
1130 hours, 29 Sept. 1969
34
-------
METROPOLITAN
LOS ANGELES
8 10 mi
Figure 9. Estimated Trajectory of Air Mass Arriving at El Monte,
1230 hours, 29 Sept. 1969
35
-------
METROPOLITAN
LOS ANGELES
8 10 mi
Figure 10. Estimated Trajectory of Air Mass Arriving at El Monte,
1330 hours, 29 Sept. 1969
36
-------
The next section summarizes our results in testing the improvements
we attempted in the chemical, mathematical, and meteorological aspects
of the problem. Indeed, many times it appeared that the earlier simple
concepts gave better results, but working with better data and more
demanding theoretical descriptions than before, we uncovered some sig-
nificant areas for further investigation.
37
-------
IV. TESTS OF THE MODEL AGAINST LOS ANGELES BASIN DATA
A. OBJECTIVES OF THE TESTS
With the development of a more complete chemical model than the
previous one, the mathematical and meteorological improvements offset
the potential computing slowdowns due to added reactions. Having
described these advances, we turn now to the task of subjecting the
model to tests that are more severe than the previous ones. El Monte
is the station of interest in the present work. It is closer to the
edge of the urbanized areas than Huntington Park or Commerce, and is,
therefore, likely to show larger random fluctuations of pollutant con-
centrations. In view of the variability of wind directions, this is
expected because the air over El Monte may have originated from either
high or low pollution areas (see Figs. 7-10). Hence the accuracy of
the initial conditions for the test calculation may be as important as
the accuracy of the flux boundary conditions.
A primary purpose of this part of the work is to isolate the quality
of the model description from possible limitations in the data base. Such
limitations can arise either from an insufficient amount or from inaccura-
cies. We selected September 29, 1969, because it was in an episode of
29
high oxidant levels, and (not coincidentally) the field programs col-
lected a broad sample of airborne as well as ground station measurements.
Owing to the prior experience of Scott Research Laboratories, all data
were of higher quality than those obtained in 1968. If we can localize
the shortcomings of the model (or any that occur in the data) we then
possess a guide for action in seeking remedies.
B. REVIEW OF THE PREVIOUS VALIDATION STUDY
A slab treatment of the central Los Angeles basin is the chief
distinguishing feature of the previous study. Unlike box models, it
accounted for vertical variations of concentration; in fact, a compari-
son with reactive box modeling demonstrated the need for considering
38
-------
vertical diffusion. The choice of a centrally located station,
Huntington Park, minimized advection effects because the pollution
history in surrounding air was not unlike that at the station. Only
when considerable ventilation occurs will this break down. To provide
added assurances of testing our hypotheses under the assumed conditions,
we averaged over a series of relatively stagnant days occurring late in
October of 1968.
Figure 11 reviews the results of using our original kinetic model
with a conventional Crank-Nicolson diffusion algorithm. This particular
plot is for reactive hydrocarbon and it shows dashed lines for the box
model results. The symbols f and r refer to the fractions assumed
for the nitric oxide flux and the propylene oxidation rate. Experimental
results of C9H»/NO -ratios suggested as much as a fourfold deficit in
£» £ X
the NO -balance; therefore, we set f = 1/4 . The aggregate of the hydro-
X
carbon reactivity was one-third to one-half that of C-H/-, so that r = 1/2
pphm
100
z
o
u
z
o
u
10
1
0600
1
1
1
1
1
0800 1000 1200
TIME OF DAY
MEAN OF OBSERVED DATA30 ^
iSTANDARD DEVIATION 1
FIVE STATION DIFFUSION MODEL
09QO i HOMOGENEOUS MIXING MODEL
WITH INITIAL CONDITIONS AT
0600 ) TWO DIFFERENT TIMES
i
1400PST
Figure 11. High Reactivity Hydrocarbon Concentration at Huntington
Park on Type 2 Days, 1968. f = 1/4; r = 1/2
39
-------
*
was employed. Note that the diffusion model is preferable to the box
model, but that advection is inadequately represented in the mid to late
morning hours.
Figure 12 shows the oxidant buildup at Huntington Park for the late
October 1968 smog episode. Again it is seen that the finite difference
diffusion model is superior to the box model. The solid curve represent-
ing a solution using a five-station vertical mesh gives nearly the correct
peak height and time.
C. SENSITIVITY STUDIES ON 1969 TRAJECTORIES WITH THE EXPANDED MODEL
Based on the semi-Lagrangian formulation of the photochemical/
diffusion model, the computed end-point composition of the air masses
depends on initial conditions, flux from the ground along the trajec-
tories, and on reaction rates. For our tests, we concentrate on El Monte
data because much of the polluted air there comes from somewhere else.
This is believed to be a more severe test of the model than that at
Huntington Park. The initial conditions are based on measurements insofar
as possible. The principal initial values for the 1030 trajectory are
as follows for 0730 hours (PST).
c = 51 pphm (6 ppm C total HC)
r = 24 q
CN02 *'-
= 3.5 pphm
31
These were interpolated from the Azusa station data from the Los Angeles
County Air Pollution Control District. Reference to Fig. 7 shows the
basis of this choice. The reactive hydrocarbon value is 'derived from the
assumptions of 34 percent reactive fraction and an average carbon number
of 4 for the reactive family of compounds. These are seasonal averages
Propylene was used for the baseline validations of smog chamber results.
40
-------
pphm
lOOr
z
o
u
z
o
10
Q.I
0600
,/,/ MEAN OF OBSERVED DATA3<>
//// ± STANDARD DEVIATION
FIVE STATION DIFFUSION MODEL
i
0900
0600
HOMOGENEOUS MIXING
MODEL WITH INITIAL
CONDITIONS AT TWO
DIFFERENT TIMES
i
0800 1000 1200
TIME OF DAY
1400 PST
Figure 12. Ozone Concentration at Huntington Park on Type 2 Days,
1968. f = 1/4; r = 1/2
41
-------
29
over the 1969 measurement months at El Monte. The analysis that
32
generated the values from gas chromatographic data is discussed elsewhere.
Similarly the 0730 initial values for the 1130 trajectory are obtained
29
from the Commerce data log. They are:
cu., = 94.5 pphm (9.7 ppm C)
rlL
CNO = 43'9 Pphm
CN00 - 17'<
,m
The hydrocarbon conversion factors here are 41 percent reactive fraction
and an average carbon number of 4.2.
Since the 1230 trajectory begins near the coast, we employ station
76 (La Cienega Boulevard) measurements for 0630 hours (PST) initial
values of nitrogen oxides. No hydrocarbon data are available from that
station and an HC/NO -ratio of 2 was chosen. Although this is not
representative of the inventory ratio, it is representative of many
morning air samples. Such an observation further reinforces the suspicion
32
of an NO -deficiency on high oxidant days. Thus, the values are:
X
CRC = 54 pphm
CNO = 18 pphm
CN02 = 9 Pphm
The location of the 1330 trajectory origin is near no station and the
data indicate that the initial pollution cannot be neglected. Conse-
quently, results were not obtained for the 1330 trajectory.
42
-------
For boundary conditions we followed nearly the same procedure as
before. The fluxes of reactive hydrocarbon and nitric oxide were
obtained by superimposing a source grid on Figs. 7 through 10. The sta-
tionary sources were distributed as in the previous work, but updated
information was utilized for the vehicular fluxes as they are disposed
in space and time. We are especially grateful to Dr. Philip M. Roth,
Director of Environmental Studies, Systems Applications, Inc., for pro-
viding these data that his organization developed under the sponsorship.
of the Air Pollution Control Office. Our own emission factors were
applied using averages weighted by the vehicle age distribution as
33
described elsewhere. Ultraviolet data from Ref. 29 were employed to
get photodissociation rates for N0_ by means of calibration functions
32
developed in the data analysis study. Inversion base altitudes were
averaged from the airborne temperature measurements provided by Scott
29
Research Laboratories.
To examine parametrically the influence of inventory levels, we
altered the rate of emissions calculated for ground level. Halving the
nominal NO-fluxes is suggested by the departures of the CO/NO -ratios
32 x
and the C-H.-/NO -ratios obtained for high oxidant days. Subsequent
£t & X
halving of both reactive hydrocarbon and nitric oxide fluxes brings us
down to the adjustment represented by f = 1/4 in Figs. 11 and 12 but
preserves the HC/NO -ratio. In both cases the full propylene oxidation
X
rates were employed.
Table 4 summarizes the results for trajectory end points. Advancing
down the table we note an ever-growing sensitivity of the results to the
flux adjustments. This occurs because the increasing wind speed lengthens
the trajectories and there is larger relative exposure of the air mass
to high pollutant fluxes.
43
-------
TABLE 4
RESULTS OF MODEL SENSITIVITY STUDY
USING 1969 EL MONTE DATA
Concentrations at El Monte in pphm
El Monte
Trajectory
Arrival
Time (PST)
1030
1130
1230
NO
Measured
Values
3.4
3.2
1.5
Fu" *BC
0.5
1.5
1.1
™£
0.3
1.5
1.0
HC
Measured
Values
43
64
38
Fun'!"0
34.3
87.4
67.5
(1/4)*NO
24
70.1
48.1
NO
2
Measured
Values
16.4
21.3
10.3
(1/2) »NO
Ful1 *HC
20.5
58.5
49.2
(1/4>«NO
<1/2>«HC
14.5
48.6
34.6
Q
3
Measured
Values
13.4
23.6
24.5
(1/2)*HO
Fu" «HC
32.6
28.4
30.9
(1/4)*NO
(1/2)»HC
35.3
24.2
24.8
NOTE: Nominal values of emission fluxes are denoted by
halved the fluxes derived from Inventories.
Hence, (1/2)$ means that we
In the 1030 case, both hydrocarbon and nitric oxide are low, but
N0? and ozone are high. Changes in the fluxes have marked effects on
everything but the ozone. Comparing these end-point compositions with
the initial values cited above, we note a strong dependence on initial
conditions accuracy for this short time scale. Fortunately the Azusa
data are available as an aid in this respect.
For the 1130 results, the ozone and nitric oxide come in closer
to the data but
N02 is extremely high. This is symptomatic of the
large departures from the radiation-supported quasi-equilibrium that
32
we noted in analyzing the El Monte data. Although the model does
not assume quasi-equilibrium among the three reactions
hv + NO,
NO + 0
NO
the differential equation solutions nearly always approximate it rather well.
-------
Again, with the 1230 trajectory the ozone and nitric oxide levels
match observations reasonably well, but hydrocarbon and NCL are both
high. Despite the passage of this air from the seashore and then over
extensive regions of the central basin, we note that the initial values
of the primary pollutants still influence the final concentrations in
a dominant manner.
D. HISTORY ANALYSIS OF THE 1030 EL MONTE TRAJECTORY
Because of the relative completeness of initial conditions that we
can relate to the Azusa station, we have chosen the 1030 trajectory to
discuss in some detail. Examining Table 4 we note an overabundance of
ozone at 1030 and a correspondingly rapid completion of NO -> NO- con-
32 1
version. Reactivity analyses and previous modeling studies suggest
reduction of the oxidation rate constants. To achieve some level of
comparative assessment with the previous work we assign one fourth the
nominal NO flux and one half the oxidation rate; hence, f = 1/4 and
r = 1/2 describe the conditions as before. This time, however, we
preserve the HC/NO -ratio as in the entries in Table 4 and reduce hydro-
carbon fluxes by a factor of 2. This means that the difference in end-
point concentrations between this case and the l/4<(>M , l/2 entries
rJU liC*
is solely due to the rate constant reduction. This differs from the
earlier work in which hydrocarbon fluxes were not reduced.
Figure 13 shows the reactive hydrocarbon history starting at Azusa
at 0730 and ending at El Monte at 1030 on September 29, 1969. Despite
the sharp reduction in hydrocarbon fluxes, the calculated curve stays
above the Azusa levels until 1000 hours when it begins to slope-off in
the observed manner. The model output clearly bears a closer relationship
to the Azusa measurements than it does to the El Monte observations.
This problem is typical of the difficulties we encountered in attempting
to model conditions at the eastern portion of the basin at El Monte.
45
-------
60 pphm. -
Z
O
*-
<
et
»—
Z
iu
Z
O
(J
9 10
TIME OF DAY
11
12 PST
Figure 13. Reactive Hydrocarbon History Along the 1030 Trajectory.
f = 1/4; r = 1/2 (Ground level concentrations)
On Fig. 14 we note a sharp drop from the interpolated NO-level
between 0730 and 0740 hours. This is reflective of the previously noted
failure of the data to approach quasi-equilibrium between NO, N0?, 0» ,
and sunlight intensity under high oxidant conditions. The NO-conversion
seems to proceed at roughly the observed rate after the transient is
absorbed in the system; however, the level ends up closer to Azusa values
than to El Monte values. Note that if we were to employ 0830 El Monte
concentrations as initial values, we would have a lower HC/NO-ratio and
could expect still slower nitric oxide conversion rates. Thus, both
nitric oxide and hydrocarbon decay more like the Azusa data than the
El Monte data.
The NO- behavior on Fig. 15 exhibits more nearly what one would
expect than do either the reactive HC or the NO . Proceeding from
the end of its transient adjustment to observed sunlight intensity,
46
-------
10 pphm
8
Z
o
u
z
O
u
I
9 10
TIME OF DAY
11
J
12PST
Figure 14.
Nitric Oxide History Along the 1030 Trajectory.
r = 1/2 (Ground level concentrations)
40 pphm -
1/4;
z
o
I—
<
Of
t—
z
111
z
o
u
20
10
e
EL MONTE
I
I
I
8
11
12PST
Figure 15.
9 10
TIME Of DAY
Nitrogen Dioxide History Along the 1030 Trajectory,
f = 1/4; r = 1/2 (Ground level concentrations)
47
-------
the air mass gradually undergoes N0--transition from Azusa levels to
El Monte levels as it meanders about in the northeastern area of the
basin. Nitric oxide conversion supported by increasing sunlight inten-
sities drives the NCL upward. The ground-based sources continue to
feed in NO as it reacts and diffuses upward.
Summing over NO and NO- , we get excellent total balance
behavior for these oxides of nitrogen. Figure 16 indicates that not
only slopes, but also magnitudes are represented rather well in the
gradual change of composition moving from initial to final conditions
over the three-hour period. Apparently the dominant effect of the good
behavior of the
the data.
NO- brings the balance into favorable agreement with
It will be recalled that this was a key test in the original
choice of f = 1/4 after conducting our original modeling study.
40 pphm
20
z
o
>-
<
at
>—
Z
UJ
Z
O
u
10
8
6
8
11
12PST
9 10
TIME OF DAY
Figure 16. Oxides of Nitrogen (NO + NO-) History Along the 1030
Trajectory, f = 1/4; r = 1/2 (Ground level concentrations)
-------
Finally, the history of ozone concentration may be seen on Fig. 17.
The measurements at both Azusa and El Monte show a remarkable degree of
smoothness and regularity compared with the jagged curves of the species
displayed in the previous figures. Again, due to the rapid transient
response of the 0_/NO/NO_-system to the sunlight ultraviolet intensity,
there is a sharp change in the first ten minutes of the model curve. As
in the case of N0_ , however, the ozone undergoes a smooth transition
from its origin to its destination. Considering the usual sensitivity
of ozone to competing rate determining factors, this degree of realism
is gratifying.
If we applied this same set of assumptions to the other trajectories,
we would get too little ozone at the end points. Table 4 indicates fair
agreement with f = 1/4 and r = 1 as shown in one of the columns.
The hydrocarbon agreement for these other cases would be even worse
than it stands if r were decreased to 1/2. Because of low confidence
in the initial values that apparently dominate, the other trajectories
would likely be less than satisfactory as tests of the model.
49
-------
40 pphm
20
Z
o
at
t-
Z
iu
Z
O
10
8
I
9 10
TIME OF DAY
11
J
12 PST
Figure 17. Ozone History Along the 1030 Trajectory, f = 1/4; r = 1/2
(Ground level concentrations)
50
-------
V. CONCLUDING REMARKS
An overview of these further developments of the GRC photochemical/
diffusion model holds answers to many questions and points the way to
future research needs. Before detailing the accomplishments and recom-
mending new directions, we might make the general observation that apply-
ing and adapting the modeling tools that we possess now appears preferable
to initiating effort on far more detailed chemical dispersion formulations.
We have maintained a balanced posture in which the model advances and the
improvements in measurements have kept pace with one another.
The chemical descriptions offered here hold the main features of
the NO -chain termination that controls composition effects on secondary
pollutant production. In concert with this, the addition of hydroxyl
radical chemistry is a first step in explaining the "excess rate" of
hydrocarbon oxidation (over that attributable to 0-atom attack). Probably
one of the most uncertain aspects of the system is our inclusion of HONO
formation (from NO + N02) and opposed by photodissociation. If this is
an important feature of the free radical balance, we really need to know
the photodissociation absorption coefficients for HONO as well as some
of the radical-molecule collisional reaction rates. In his task force
34
report for Project Clean Air, H. S. Johnston stresses the need to con-
sider nitrous and nitric acids in the photochemical mechanism. This
suggests a very real need to monitor these compounds in the laboratory
experiments and in the atmosphere. Their involvement in aerosols and
surface reactions may well hold the key to the apparent anomalies in
the nitrogen balance. The identification of chemo-mutagenic effects of
HONO is yet another reason that it should be investigated. The brief
examination of hydrocarbon synergism shows at least one way that the
radical chains are cross-linked to influence combined reaction rates.
Experimental data on mixed hydrocarbons will help sort out the most
likely mechanisms. Likewise, our small study on possible CO-interference
indicates that suggested values of these rate constants indicate a very
strong influence on the propylene/nitric oxide system at concentration
51
-------
levels as high as 100 ppm. It may be that atmospheric hydrocarbons that
are not involved in the OH-cycle are far less affected. In any event
atmospheric levels of CO are much lower than 100 ppm except in the
vicinity of traffic arteries.
Reacting to the need for extended chemical representations, we have
applied Fade" approximant techniques to the numerical integration algorithms,
Speed increases of at least four or fivefold are obtained in the computa-
tion comparing the Fade" method for the more complex chemical system with
previous methods for the simpler system. For equal-sized chemical
matrices the improvements would most likely be even better. In the
interest of making these advances available to the community at large we
show the mathematical steps in some detail.
Besides the mathematical improvements, the atmospheric model has
been adapted to a semi-Lagrangian formulation. By following selected
air masses we avoid commitments of large quantities of memory and the
incursions of artificial diffusion errors. Most important, we do not
end up with stacks of computer printout that relate to regions where
there are no measurements. It should be added that predictive calcula-
tions will become ever more useful, but our present levels of resources
and sophistication demand that effort be concentrated on validation.
Only in this way can the confidence be built that is needed for applying
modeling techniques to implementation planning.
Validation tests of the model for 1969 Los Angeles basin data
turned up several useful findings. For time intervals within a diurnal
scale, we must use extreme care in selecting initial conditions. A way
to overcome this need is to compute for several days of real time on a
continuous basis. Increasing the time scale makes sense because episodes
tend to cycle through several days. The cost in computing will go up
more than proportionately, however, because as time scales increase so
must the spatial domain. Perhaps 100 kilometers of urban/rural influence
52
-------
zone might be needed in lateral expansion. Correspondingly, a larger
vertical field will have to be included to account for slow transport
within the inversion layer, if there is one. Models that require arti-
ficial image flows above the marine layer or that assume artificially
inflated "background" levels at the edge of the net will no longer be
useful. What used to be boundaries will now be included in the
computational field.
The sensitivity tests in the validation studies confirmed some of
32
the suspicions that we expressed earlier; namely, that the quasi-
equilibrium relationship between ozone and the oxides of nitrogen does
not seem to be recovered in the data. The largest departures are for
the highest ozone levels. Attempting to represent the physical setting
in a consistent manner, we find it difficult to use the measured ultra-
violet intensity to account for the observed ozone buildup and NO-
conversion simultaneously. The inconsistency even appears in the initial
behavior of a modeling run as a transient "induction" process that
rapidly adjusts concentrations to satisfy the rate equations. Unlike
some models, ours does not impose the quasi-equilibrium, but rather
solves for time history of species with 0-atom, RCL-radicals, and OH
in the stationary state.
Because of the extreme dependence on initial conditions, our history
analysis concentrates on an air mass with relatively well-defined concen-
trations at the beginning and the end of its travel. Giving it the
initial values we see the concentrations unfold as the air parcel moves
through the computed simulation procedure. In view of the sensitivities
discovered, the transition of oxidant species 0_ and N0_ proceeds
better than one might expect. The previously adopted biases on the
NO-flux and the propylene oxidation rates were confirmed in this run
having very different conditions than those in Huntington Park represented
by 1968 data.
53
-------
Much more work remains to be done. We have so many uncertainties
in both the data and the model assumptions, that systematic tests must
be devised to isolate the individual influences of the various parameters.
Our very limited sensitivity study illuminates some of these very serious
questions that need to be answered. The anomalous pattern of NO-flux
followed from our original work, through the data analysis project, and
into this work. Either some very serious deficiency exists in the trans-
port formulation or there are some rapid loss mechanisms for oxides of
nitrogen that reduce the apparent emission strengths. In view of the
extensive efforts directed toward NO -controls high priority must be
X
given to measuring the terms of the nitrogen balance equation in an urban
area. Nitric acid vapor, nitrous acid vapors, and particulate nitrates
could hold the key to this question. The modeling is not precise enough
to pin down the deficit, but the CO/NO and C?H /NO analyses of the
oo " £• £• fi
Los Angeles data showed significant deviations on high oxidant days.
Since the same types of questions have appeared in laboratory experiments,
it seems that a greater latitude of nitrogen-bearing compounds should be
sought in future atmospheric studies.
The plans announced for an APCO regional study should hold many of
the solutions to the problems identified here. Conversely, application
of the model should provide valuable feedback to such programs. Clarifi-
cations in the chemical mechanism should take the form of better values
of rate constants for the controlling reactions. Perhaps the very selec-
tion of which reactions to include might be altered for the air contaminant
hydrocarbon mix. Mathematical techniques seem well enough advanced to
match the simpler formulations of photochemical modeling schemes to the
computing machinery available. The more complex schemes proposed recently
will continue to shed light on questions of dispersion parameters and
inventory accuracy.
-------
Continued use of the model we have developed should generalize it
to a point where it can be integrated in the process of implementation
planning. The job ahead is a difficult one and only a judicious mixture
of experimental and theoretical effort will lead to a useful tool that
is consistent with the requirements facing us now.
55
-------
56
-------
REFERENCES
1. A. Q. Eschenroeder and J. R. Martinez, A Modeling Study to
Characterize^Photochemical Atmospheric Reactions to the Los
Angeles Basin Area, General Research Corporation CR-1-152,
November 1969.
2. R. C. Wanta, "Meteorology and Air Pollution," Air Pollution, Vol. 1,
2nd ed., A.C. Stern, Ed., Academic Press, New York, 1968, Ch. 7.
3. H. Moses, Mathematical Urban Air Pollution Models, Paper 69-31
62nd Annual Meeting of the Air Pollution Control Association, New
York, N.Y., June 1969.
4. J. H. Seinfeld, Mathematical Models of Air Quality Control Regions,
Symposium on the Development of Air Quality Standards, Los Angeles,
Calif., October 1969.
5. S. Calvert, "A Simulation Model for Photochemical Smog," California
Air Environment, Vol. 1, No. 3, Jul-Sep 1969, pp. 1-3.
6. L. Wayne, R. Danchick, M. Weisburd, A. Kokin, and A. Stein, Modeling
Photochemical Smog on a Computer for Decision-Making, Paper 70-18,
Air Pollution Control Association 63rd Annual Meeting, St. Louis,
Mo., June 14-19, 1970.
7. S. K. Friedlander and J. H. Seinfeld, "A Dynamic Model of Photo-
chemical Smog," Environmental Science and Technology, Vol. 3,
No. 11, November 1969, pp. 1175-1182.
8. R. C. Sklarew, A New Approach: The Grid Model of Urban Air
Pollution. Air Pollution Control Association Paper 70-79, 63rd
Annual Meeting, St. Louis, Missouri, June 14-18, 1970.
9. P. Roth, Development of a Simulation Model for Estimating Photo-
chemical Pollutants, Presentation at an Information Meeting of
Atmospheric Diffusion Modeling Studies by NAPCA Contractors,
Raleigh, N.C., November 16-17, 1970.
10. A. P. Altshuller and J. J. Bufalini, "Photochemical Aspects of Air
Pollution: A Review," Photochemistry and Photobiology, Vol. 4,
1965, pp. 97-146.
57
-------
REFERENCES (Cont.)
11. E. A. Schuck and G. J. Doyle, Photooxidation of Hydrocarbons in
Mixtures Containing Oxides of Nitrogen and Sulfur Dioxide, Air
Pollution Foundation Report No. 29, San Marino, Calif., October 1959,
12. D. H. Stedman, E. D. Morris, Jr., E. E. Daby, H. Niki, and B.
Weinstock, The Role of OH Radicals in Photochemical Smog, American
Chemical Society Division of Water Air and Waste Chemistry, Chicago,
Illinois, September 13-18, 1970.
13. J. R. Holmes, A. D. Sanchez, and A. H. Bockian, Atmospheric Photo-
chemistry: Some Factors Affecting the Conversion of^ NO to NO^,
Pacific Conference on Chemistry and Spectroscopy, San Francisco,
October 6-9, 1970.
14. K. Westberg and N. Cohen, The Chemical Kinetics of Photochemical
Smog as Analyzed by Computer, AIAA Third Fluid and Plasma Dynamics
Conference Paper No. 70-753, Los Angeles, June 29-July 1, 1970.
15. A. P. Altshuller, S. L. Kopczynski, W. A. Lonneman, T. L. Becker,
and R. Slater, "Chemical Aspects of the Photooxidation of the
Propylene-Nitrogen Oxide System," Environmental Science and
Technology. Vol. 1, No. 11, Nov. 1967, pp. 899-914.
16. A. P. Altshuller, and J. J. Bufalini, "Photochemical Aspects of
Air Pollution: A Review," Environmental Science^and Technology,
Vol. 5, No. 5, January 1971, pp. 39-64.
17. J. H. Seinfeld, private communication, letter dated June 2, 1970.
18. A. P. Altshuller, S. L. Kopczynski, D. Wilson, W. Lonneman, and
F. D. Sutterfield, "Photochemical Reactivities of n-Butane and
Other Paraffinic Hydrocarbons," Journal of the Air Pollution
Control Assoc.. Vol. 19, No. 10, October 1969, pp. 787-790.
19. M. W. Korth, A. H. Rose, Jr., and R. C. Stahman, "Effects of
Hydrocarbon to Oxides of Nitrogen Ratios on Irradiated Auto
Exhaust, Part 1," Journal of the Air Pollution Control Assoc.,
Vol. 14, No. 5, May 1964, pp. 168-175.
20. W. G. Agnew, "Automotive Air Pollution Research," Proc. Roy. Soc.,
Series A, Vol. 307, 1968, pp. 153-181. (Data reproduced from APCA
paper by C. S. Tuesday, B. A. D'Alleva, J. M. Heuss, and G. J.
Nebel, June 1967).
58
-------
REFERENCES (Cont.)
21. J. Heicklen, K. Westberg, and N. Cohen, The Conversion of NO to NO^
in Polluted Atmospheres, Pennsylvania State University Center for
Air Environment Studies Publication 115-69, July 1969.
22. K. Westberg, N. Cohen, and K. W. Wilson, "Carbon Monoxide: Its
Role in Photochemical Smog Formation," Science, Vol. 171, March 12,
1971, pp. 1013-1015.
23. W. F. Ames, Nonlinear Partial Differential Equations in Engineering,
N.Y. Academic Press, 1965, pp. 341-342.
24. R. S. Varga, "On Higher Order Stable Implicit Methods for Solving
Parabolic Partial Differential Equations," J. Math and Physics,
Vol. 40, 1961, pp. 220-231.
25. D. Magnus and H. Schecter, Analysis and Application o^f the Fade"
Approximation for the Integration of Chemical Kinetic Equations,
Technical Report 642, General Applied Science Laboratories, Inc.,
1967.
26. S. S. Wu, "A Study of Heat Transfer Coefficients in the Lowest
400 meters of the Atmosphere," Journal of Geophysical Research,
Vol. 70, No. 8, April 15, 1965, pp. 1801-1807.
27. C. R. Hosier, "Vertical Diffusivity from Radon Profiles," Journal
of Geophysical Research. Vol. 74, No. 28, December 20, 1969,
pp. 7018-7026.
28. G. De Zorzi, "Applicazione della Equazione di Diffusione allo
Studio della Bassa Atmosfera," Rivista di Meteorologica Aeronautica,
Vol. 30, No. 1, March 1970, pp. 3-28.
29. Atmospheric Reaction Studies in the Los Angeles Basin, Vols. I-IV,
Scott Research Laboratories, February 1970.
30. Final Report on Phase I, Atmospheric Reaction Studies in the
Los Angeles Basin, Vols. I and II, Scott Research Laboratories,
June 30, 1969.
31. Monitoring Data from Los Angeles County Air Pollution Control
District, 1969.
59
-------
REFERENCES (Cont.)
32. A. Q. Eschenroeder and J. R. Martinez, Analysis of Los Angeles
Atmospheric Reaction Data from_19j68 and 1969, General Research
Corporation CR-1-170, July 1970.
33. A. Q. Eschenroeder, Some Preliminary Air Pollution Estimates for
the Santa Barbara Region 1970-1990, General Research Corporation
IMR-1278, March 1970.
34. H. Johnston, "Reaction in the Atmosphere," Project Clean Air Task
Force Assessments, Vol. 4, Task Force No. 7, Section 3, University
of California, September 1, 1970, p. 3-1.
60
------- |