3SR-844
VOLUME I
NOVEMBER 1971
SYSTEMS, SCIENCE AND SOFTWARE
P.O. BOX 1620, LA JOLLA, CALIFORNIA 92037, TELEPHONE (714) 453-0060
A PARTICLE-IN-CELL METHOD FOR NUMERICAL SOLUTION
OF THE ATMOSPHERIC DIFFUSION EQUATION,
AND APPLICATIONS TO AIR POLLUTION PROBLEMS
BY
R. C. SKLAREW, A. J. FABRICK AND J. E. PRAGER
FINAL REPORT
FOR THE
DIVISION OF METEOROLOGY
NATIONAL ENVIRONMENTAL RESEARCH CENTER
RESEARCH TRIANGLE PARK, NORTH CAROLINA 27711
UNDER
EPA CONTRACT NO. 68-02-0006
-------
1OOR71015
3SR-844
VOLUME I
NOVEMBER 1971
SYSTEMS, SCIENCE AND SOFTWARE
P.O. BOX 1620, LA JOLLA, CALIFORNIA 92037, TELEPHONE (714) 453-0060
A PARTICLE-IN-CELL METHOD FOR NUMERICAL SOLUTION
OF THE ATMOSPHERIC DIFFUSION EQUATION,
AND APPLICATIONS TO AIR POLLUTION PROBLEMS
BY
R. C. SKLAREW, A. J. FABRICK AND J. E. PRAGER
FINAL REPORT
FOR THE
DIVISION OF METEOROLOGY
NATIONAL ENVIRONMENTAL RESEARCH CENTER
RESEARCH TRIANGLE PARK, NORTH CAROLINA 27711
UNDER
EPA CONTRACT NO. 68-02-0006
-------
3SR-844
ACKNOWLEDGMENTS
The authors gratefully acknowledge the technical as-
sistance of S3 scientists Burton Freeman, Manuel Rotenberg
and Albert Petschek; the comments and advice of the project
technical monitor, Mr, Kenneth Calder, Division of Meteorol-
ogy, Environmental Protection Agency, and the many technical
discussions and required input data supplied by Mr, Alan
Eschenroeder, General Research Corporation, Mr. Raymond Dick
son, Air Resources Laboratory, National Oceanographic and
Atmospheric Administration, Mr. Robert Lamb, University of
California at Los Angeles, Mr. Philip Roth, Systems Applica-
tions, Inc., and Mr, Gene Start, Air Resources Laboratory,
National Oceanographic and Atmospheric Administration, This
study would not have been possible without the support and
encouragement of the S3 management.
-------
3SR-844
TABLE OF CONTENTS
Page
VOLUME I
1.
2.
3.
4.
INTRODUCTION
DESCRIPTION OF THE PICK METHOD
1.1 The Diffusion Equation .
1.2 The PICK Method
1,3 Simulation of Emissions and Chemical
Reactions
1.4 One -Dimensional Example of the PICK Method .
1,5 Accuracy of the PICK Method
1.6 Summary of the PICK Method
TESTING THE PICK METHOD
2.1 Description of the PICFIC Code
2,2 Comparison of Methods for Pure Advection . .
2.4 PICFIC Simulation
2 . 5 Summary
NEXUS SIMULATION OF CARBON MONOXIDE IN LOS ANGELES
3.1 Description of the NEXUS Model
3.2 Input Data
3.3 Results of NEXUS Simulation
3.4 Testing of NEXUS
3.5 Summary and Conclusions
NEXUS/P SIMULATION OF PHOTOCHEMICAL SMOG
4.1 Description of NEXUS/P
4.2 Photochemical Simulation
4.3 NEXUS/P Testing
4.4 Input Data
1
3
j
9
15
16
23
25
27
27
32
47
55
63
67
70
78
89
104
116
119
120
123
130
134
111
-------
3SR-844
TABLE OF CONTENTS, contd.
Page
4.5 Results of the NEXUS/P Simulation 138
4.6 Summary and Conclusions 149
5. SUMMARY AND CONCLUSIONS 157
REFERENCES 161
VOLUME II - PICFIC and NEXUS Computer Code Listings
VOLUME III - NEXUS/P Operations Manual
IV
-------
3SR-844
LIST OF ILLUSTRATIONS
Numb e r Title Page
1 Calculation of cellular concentration from
particle masses 11
2 Area-Weighting interpolation for total velocity 13
3 Gaussian distribution with the approximate cell
concentrations 17
4 Macro flow chart of the PICFIC code 28
5 PICFIC advection simulation test with a = 2.5
Km, contours at 0.1, 10, 20, 40 and 60 .... 34
6 First-order finite difference advection test
with a = 2 .5 Km 36
7 Fromm fourth-order method advection test with
a = 2 . 5 Km 41
8 PICFIC advection test with a = 1 Km 43
9 Fromm fourth-order method advection test with
a = 1 Km 44
10 Crowley fourth-order method advection test with
a = 1 Km 45
11 Curve fitting with a polynomial. Pictured is a
Gaussian curve and an approximating cubic poly-
nomial. This picture corresponds to the poly-
nomial fitting used in calculating the fourth-
order flux scheme for the a = 1 initial distri-
bution. The shaded areas correspond to the
errors in the fit 46
12 Histogram of initial cell concentrations com-
pared to normalized Gaussian distribution ... 53
13 Renormalized 41st cycle of PICFIC calculation
(both advection and diffusion) compared to
normalized Gaussian distribution 54
14a Pollutant cloud advection in constant wind
field with no diffusion. Three time instants
are displayed together for comparison 56
v
-------
3SR-844
LIST OF ILLUSTRATIONS, contd.
Number Title Page
14b Same as Figure 14a. The contours are for
integer powers of 2 in pollutant concentra-
tion 56
ISa Pollutant cloud advection in constant wind
field with spherically symmetric diffusion . . 57
15b Same as Figure 15a. Contours at integer powers
of two 57
16a Pollutant cloud advection in constant wind
field with aspherical constant diffusion ... 58
16b Same as Figure I6a. Contours at integer powers
of two 58
17a Pollutant cloud advection in constant wind
field with aspherical, position-dependent dif-
fusion. Model of an inversion at the source
height 60
17b Same as Figure 17a. Contours at integer powers
of two 60
18a Pollutant cloud advection in shear wind field
with spherically symmetric diffusion, three
separate times are depicted as the cloud moves
to the right 61
18b Same as Figure 18a. Contours at integer powers
of two 61
19 Calculational grid for NEXUS CO simulation . . 69
20 Macro flow chart of NEXUS main program .... 71
20a Macro flow charts of NEXUS subroutines INPUT
and RESTRT 72
20b Macro flow charts of NEXUS subroutines SETUP
and OUTPUT 73
20c Macro flow charts of NEXUS subroutines DIFFUS
and PARCEL 74
vi
-------
3SR-844
LIST OF ILLUSTRATIONS, contd.
Number Title Page
20d Macro flow charts of NEXUS subroutines BORDER
and SOURCE 75
21 Supplied wind data (Ref. 16). The streamlines
are shown by solid lines. Isotacks in miles
per hour are dashed contours 79
22 Winds in NEXUS grid. The arrowhead points in
the direction of motion and the arrow lengths
is proportional to the wind speed 80
23 Initial (3:00 A.M.) CO surface concentration
(ppm) (Ref. 16) 82
24 Los Angeles County traffic distribution (un-
normalized relative car miles per cell) .... 83
25 Temporal distribution for Los Angeles Basin
traffic 84
26 Macro flow chart for SETUP 85
27 Isopleths for CO at 10:00 A.M 90
28 Plot of particle positions at 10:00 A.M. ... 91
29 Isopleths for CO (in ppm) at 10:00 A.M. as
simulated by Lamb's model (Ref. 16) 95
30 Surface winds at 10:00 A.M. showing convergences
near Downtown and Long Beach 96
31 Simulated hour-averaged CO concentrations com-
pared to measured concentrations at Downtown,
Los Angeles, Lennox, Long Beach and Hollywood . 98
32 Surface CO concentrations used for 3:00 A.M.
initial concentrations 99
33 At 3:00 A.M. the winds were calm and generally
toward the coast 100
34 High winds throughout most of the basin with
convergence near Burbank 102
Vll
-------
3SR-844
LIST OF ILLUSTRATIONS, contd.
Number Title Page
35 The morning high concentrations have been
blown into the San Gabriel Mountains and
Puente Hills but a new build-up is forming
in the Northwest 103
36 Evening traffic has produced a local build-up
of CO over Long Beach and a large region of
high concentrations corresponding to the sur-
face wind convergence in the Burbank-Pasadena-
Hollywood area , 105
37 Temporal distribution of traffic for freeways
and surface streets (Ref. 23) 109
38 Spatial distribution of freeway traffic (thou-
sand vehicle miles per day) (Ref. 23) 110
39 Spatial distribution of non-freeway traffic
(thousand vehicle miles per day) (Ref. 23) . . Ill
40 Hour-averaged CO (in ppm) calculated with an
inversion rising height is compared to the
standard results with a fixed inversion height
and to the measured concentrations 113
41 Macro flow diagram for computing photochemical
pollution 122
42a Initial distribution of primary pollutant A . . 131
42b Distribution of primary pollutant A after ad-
vection 132
42c Distribution of secondary pollutant, B, calcu-
lated from cell concentrations given in
Figure 42b 133
42d Distribution of secondary pollutant, B, calcu-
lated from Lagrangian concentration 133
43 Inversion height (feet above sea level) at 12
noon on September 30, 1969 (Ref, 17) 137
44 NO concentration in pphm at 6;00 A.M 144
Vlll
-------
3SR-844
LIST OF ILLUSTRATIONS, contd.
Number Title Page
45 N02 concentration in pphm at 6:00 A.M 145
46 N02 concentration in pphm at 10:00 A.M 146
47 0~ concentration in pphm at 10:00 A.M 147
48 HC concentration in ppm at 10:00 A.M 148
49 NO- concentration in pphm at 3:00 P.M 150
50 0, concentration in pphm at 3:00 P.M 151
51 HC concentration in ppm at 3:00 P.M 152
52 N02 concentration in pphm at 6:00 P.M 153
53 HC concentration in ppm at 6:00 P.M 154
54 0, concentration in pphm at 6:00 P.M 155
IX
-------
3SR-844
LIST OF TABLES
Table Title Page
I Calculated central values after 45 cycles for
different methods Ccorrect result: 69) .... 40
II PICFIC diffusion evaluation varying the number
of particles per cell using the area-weighting
procedure ... ......... . ...... 49
III PICFIC diffusion evaluation varying the cell
size using 1000 mass point particles ..... 51
IV Results of calculation compared to analytic
solution of Neuringer ............. 62
V 3:00 A.M. to 8; 00 P.M. average carbon monoxide
concentrations (£0 in ppm) , , ...... . . 92
VI Results of series of tests with NEXUS ..... 117
VII Rate coefficients for Eschenroeder rs model of
the hydrocarbon/nitric oxide mechanism .... 124
VIII Gaussian puff (a = 2Ax, a = 2Ay, a » 2Az)
central cell concentrations ..... ..... 135
IX Regional average of stations reporting .... 139
X Los Angeles APCD stations reporting most
plete data -- Station 1, Downtown Los Angeles . 140
Station 60, Azusa ........ 141
Station 79, Pasadena ...... 142
-------
3SR-844
INTRODUCTION
This paper reports the development and initial appli
cations of a new method for the solution of the turbulent at-
mospheric diffusion equation. The method, called PICK, is
based on the use of Lagrangian mass points and is one of a
family of Particle-In-Cell techniques for the solution of
partial differential equations. The turbulent atmospheric
diffusion equation is the standard description of the behav-
ior of pollutants in the atmosphere and solutions of the
equation have been used extensively for air pollution simula-
tion and prediction.
Other methods employed for the solution of the equa-
tion include analytic solutions and finite difference methods.
An analytic solution of the equation can only be obtained with
simplifying assumptions on the winds and diffusivities which
limit the generality of the solution. Finite difference
methods have difficulties in treating large gradients in the
solution and also may introduce a truncation dispersion of the
solution that masks the atmospheric turbulent diffusion. The
PICK method is the most accurate general technique for solu-
tion of the equations. With the PICK method a solution can
be obtained for any specified condition, and the method does
not suffer from the difficulties encountered with finite dif-
ference methods.
The purpose of this study was the development of the
PICK method (Section 1) and the demonstration of the method in
-------
3SR-844
the solution of evaluation test cases and actual air pollu-
tion problems. Test cases for the evaluation of feasibility
and accuracy and for comparison to finite difference solu-
tions were conducted with a two-dimensional computer code
PICFIC, as described in Section 2, For actual air pollution
studies, the PICK method was used in the three-dimensional
code NEXUS, The description of NEXUS and its application to
the simulation of CO in Los Angeles is given in Section 3,
The computer codes PICFIC and NEXUS were developed completely
by Systems, Science and Software CS3) as proprietary S3
computer software,
The description of photochemical smog requires terms
for the chemical reactions to be added to the turbulent at-
mospheric diffusion equations. These additional terms are
non-linear and couple the equation for each pollutant species,
The NEXUS/P code was developed to solve the equations with
the photochemical terms. In Section 4 the code NEXUS/P is
described and its application to photochemical smog in Los
Angeles is reported.
-------
3SR-84-4
1. DESCRIPTION OF THE PICK METHOD
1.1 THE DIFFUSION EQUATION
This report describes a method for the numerical solu-
tion of the equation of turbulent atmospheric diffusion using
Lagrangian mass points (called particles). The method is
called "PICK" to emphasize the Particle-In-Cell nature of the
technique as well as the K-theory approximations used in
developing the atmospheric diffusion equation. The atmos-
pheric diffusion equation is the standard mathematical de-
scription of the motion of pollutants in the atmosphere.
This PICK method has been applied to the simulation of air
pollution in two and three dimensions. The applications will
be described in the following sections.
Particle-In-Cell methods have been developed and used
in the PIC computer code for simulating compressible hydro-
dynamics and in the MAC (Marker And Cell) computer code for
2
incompressible fluid dynamics. These applications have been
reviewed by Harlow. In the PIC code, mass is associated
with each particle and as it crosses a cell boundary it trans-
ports cell averaged momentum and energy. Cell quantities for
density, pressure, velocity and temperature are calculated
from the particles and an equation of state for the fluid.
The particle positions are changed satisfying mass, momentum
and energy conservation during each time increment. In the
MAC code, the particles are Lagrangian markers that are used
to delineate fluid boundaries or free surfaces, or to tag the
-------
3SR-'844
fluid. The MAC particles are Lagrangian positions and do not
enter directly into any calculation, whereas the PIC particles
are used to represent the fluid.
Diffusion is not simulated in either the PIC or the
MAC code. These codes were developed using particles, in
part, to eliminate the dispersive errors inherent in other
methods. The PICK method introduced in this report also has
been developed to eliminate dispersive errors, those found in
finite difference solutions of the atmospheric diffusion
equation. The particles in PICK represent pollutant mass.
Changes in particle position are calculated to simulate pollu-
tant mass transport due to both advection and diffusion as
specified by the turbulent atmospheric diffusion equation.
The equation describing turbulent atmospheric diffu-
sion.
8c - TT9c vi£ - w9c + 9 *K 9c
9T ~ U97 V9y WSY 3x|Kx 3x
9C
9 c
equates the time rate of change of concentration c , •%•£ ,
_g„ gc _g _ ° t
to advective rate of change of c , ~ujT5T ~ v-g~ ~ wg~ > due to
the mean velocity field (u,v,w) plus the divergence of the
turbulent flux with the position dependent eddy diffusivity
(Kx
,K
' y
'V
'a_|
9x {
9c
x 9x
*fx(
Ky
i£
Since we will be considering general wind fields, there is no
assumption that the x-direction diffusion term is negligible.
The PICK method could actually be developed for a tensor eddy
diffusivity as discussed by Calder, but for simplicity and
-------
3SR-844
clarity the coordinate axis have been assumed to be the prin-
ciple axes for the diffusivity tensor. Discussions of the
derivation of this equation are given in several sources, e.g,
Sutton, Section 4 and of the applications to atmospheric dif-
fusion in Sutton, Section 8 and Pasquill, Section 3.
The diffusion equation can be re-written in a form
suitable for the PICK method by use of the following defini-
tions. The turbulent flux, for example, in the x-direction
is -K v°- , so an equivalent velocity, the "turbulent flux
X a X
velocity," can be defined such that
ufc = ~Kx H Cl,2aj
or
.. - Kx
In a similar manner, in the other directions, the turbulent
flux velocity will be
and
rH • C1'2d)
If a tensor eddy diffusivity were used, the turbulent flux
velocity could be written as
Kj ^ 9
*f
u,. = - - - , (1.2eJ
fi ' * c 9x. '
J
-------
3SR-844
where the subscripts are used for the three directions. This
is the novel feature at the heart of the PICK method - the
diffusion terms are incorporated into a velocity. With these
definitions, Eq . (1.1) becomes
It * =f
or by substituting the divergence of the mean wind fluxes,
(u~c, vc and w~c) :
l£ + Hue) 3 (vc) 3 (we) /3u . 9^ . 3w\ =
3t 3x 3y 3~z c\¥x 3y 3z/
Cu<0
3y 3z/ 3x
8(v£c) 3(wfc)
Tz
(1.3bJ
Now and in all further discussions we will restrict
our considerations to a divergence-free mean velocity (an in-
compressible mean flow) field which can be expressed as
That is, the mean velocity field (u,v,w) is solenoidal.
Thus, the term multiplying the concentration in Eq.. (1.3b)
is zero and the remaining terms can be regrouped to give
ac ^ 3(u+uf)c 3(v+v£)c 9(w+w£)c
Tt = " "33c "3y Jz '
It is apparent in this form of the diffusion equation that the
concentration is being transported by the sum of the mean and
-------
3SR-844
turbulent flux velocities. This sum, called the "total equiv
alent transport velocity," is
U = U + Mr
v = v + v£ (1.6)
W = W + Wr
The diffusion equation then reduces to the form
H * HP * !r * Ir1 - ° • (1-7
This equation is identical in form with the equation of con-
tinuity for a general compressible fluid, namely,
where p denotes the fluid density. It is this special fact
that is the basis for the PICK method as applied to the diffu-
sion problem. By the device of introducing the "turbulent
flux velocity" through Eq. (1.2) and the artifact of the "total
equivalent transport velocity," the original problem of turb-
ulent atmospheric diffusion is transformed into one describing
the advective changes of fluid density in a compressible fluid
moving in the fictitious velocity field (u,v,w) of total
equivalent transport velocities. Since the continuous field
variable p representing the fluid density can be visualized
in a discrete manner as the total mass of elementary fluid
particles in a unit volume of space, while the changes of
density may be considered in terms of changes in the number
-------
3SR'-844
of elementary mass particles occupying the unit volume, the
fundamental basis for the particle-in-cell treatment becomes
evident. The mass particles follow the fluid motion in the
fictitious velocity field (u,v,wj . That is to say, they are
Lagrangian particles in the non-solenoidal field of total
equivalent transport velocity. Their number in any cell
will determine the concentration of pollutant for the origi-
nal diffusion problem.
Eq . (1.7) will be the working equation for the PICK
method and will be referred to as the diffusion equation later
in this paper.
To complete the discussion, the boundary conditions
for the original problem must be transformed into conditions
within the field of fictitious total velocities. For the
original problem governed by Eq . (1.1), the boundary conditions
will involve specification of the domain boundaries of the
vector of total material flux, i.e., the vector f having com-
ponents :
Fx -
If
Evidently,
Fx = uc
Fy = vc (1.10)
= we
-------
3SR-844
so that the vector ? is directly proportional to the fic-
titious total velocity vector. For typical boundary condi-
tions the flux normal to a boundary is required to be zero
(simulating an impervious surface), or to equal a specified
value (simulating a transmittive boundary with flow into the
region) or to be proportional to the local concentration
(simulating reactions or absorption at a surface). The next
section describes how the diffusion equation is solved and
how these boundary conditions are handled in the PICK method.
1.2 THE PICK METHOD
The diffusion equation as written in Eq. (1.7) ex-
presses the conservation of matter for pollutant when consid-
ered as a compressible fluid being advected by the field of
fictitious total velocity. It expresses the fact that for any
small fixed volume of space the rate of increase of pollutant
concentration must equal the net pollutant flux that is pro-
duced by the fictitious total velocity through the boundary
surface.
In the PICK method the spatial distribution of pollu-
tant is represented by means of a large number of Lagrangian
mass particles, i.e., ones of constant mass of pollutant that
are simply advected in the fictitious field of total velocity.
Physical space is divided into cells of a fixed Eulerian grid
and the particles carry pollution from cell to cell as they
are moved by the fictitious velocity field. Since this field
is a non-solenoidal one, it will cause particles to move apart
or together, and thus uneven distributions may result. In
order to simulate satisfactorily the spatial distribution of
pollution, a sufficiently large number of particles must be
used in each grid cell. This requirement is examined and
quantified in a later section.
-------
3SK-844
Concentrations in the PIC method are computed as cell
averages. Two prescriptions for calculating cell concentra-
tions from the particle masses have been tried. In the
simplest prescription the cell concentration is the sum of the
masses of all particles in the cell divided by the cell volume.
Thus, when a particle crosses a cell boundary, the particle
mass is added to the cell entered and subtracted from the cell
it has left. The second prescription is equivalent to an area
(or volume) apportionment of a particle's mass between cells
based upon overlap considerations. For this purpose the
particle mass is regarded as uniformly spread over an area (or
volume) the size of a cell and centered on the particle posi-
tion. The overlap with adjacent cells then determines the
mass apportionment between them. These alternatives are illu-
strated for two dimensions in Figure 1. The first alternative,
being simpler, is faster to compute while the area-averaging
procedure has the effect of smoothing out the artificial gra-
dients caused by quantizing the pollutant mass into particles.
The results of calculations comparing these alternatives will
be discussed in Section 2. In both cases the values of the
cell concentration is associated with the cell center. Ini
tially , the particles are placed at random within each cell.
The number of particles in a cell is determined by the initial
concentration in the cell.
The fictitious total velocity field that is used to '
transport the particles in the PICK method is defined by its
values at the centers of the cells. For each such center,
the mean velocity vector (u,v,w) for the original real
velocity field must be specified and the turbulent flux
velocity calculated utilizing Eq. (1.2) with the derivatives
approximated by a finite difference algorithm. For example
in the x-direction the simplest finite differencing corresponds
to
10
-------
3SK-844
Cellular concentration in
Cell (2,2) = 2.0
—Particles having infinitesimal extent,
Cellular concentration in
Cell (1,1) = 0.04
(1,2) = 0.36
(2,1) = 0.29
(2,2) = 1.01
(1,3) = 0.10
(2,3) = 0.20
1
I .10
, .23
0.06-
0.04—*L
.20
.47
J"!
.541
(b) — Particles averaged over cell volume.
Figure 1 — Calculation of Cellular Concentration
from Particle Masses.
11
-------
3SR-844
/! l£
lc 8x
Ci + l
Ci
2Ax C-
(1.11)
where the concentration at the cell center x = x. is C^
and (x-+i x-_i) = 2Ax . Alternative differencing
algorithms are possible (such as directly approximating the
logarithmic derivative), but in view of the accuracy achieved
by this simple algorithm (as will be discussed in Section 2)
no others have been used.
The fictitious total velocity at each cell center is
then given by the sum of the original mean velocity and the
turbulent flux velocity. Once this has been calculated for
any elementary time-interval (t,t+At) the particles in the
cells are then moved for this interval with a velocity ob-
tained by linear interpolation according to the position of
the particle between the centers of adjacent cells, using the
velocity at each center. This is illustrated for two dimen-
sions in Figure 2. Using v to denote the total velocity
vector (u,v)
v(x,;
(y -
x.)
Vi,
Ax
v (x - x ~) + v fx
i+l,j i iji
+1 - X)l
Ax
Ay
(1.12)
12
-------
3SR-844
Figure 2 — Area Weighting Interpolation for Total Velocity
The shaded rectangular area is cell-sized and centered at the
particle position. Rearranging the terms in Eq. (1.12) gives
- x)(y -
x.)Cyj+1
x)(yj
(1.13)
is the
The product (x - x^)(y - y.) multiplying vi+1 -+1
shaded area overlapping cell (i+l,j+l) in Figure 2. Simi-
larly, the other products give the areas of the overlap with
13
-------
3SR-844
the other cells. In three dimensions the analogous volume
overlap is used.
In the PICK method the calculation to advance the
particle configuration in time proceeds in steps or cycles,
each of which calculates the desired quantities for time
t + At in terms of those at time t (an "explicit" time
advancement procedure)
x(t + At) = x(t) + uAt
(1.14)
y(t + At) = y(t) + vAt
The velocities are the fictitious total velocities determined
for the beginning of the time interval and interpolated to
initial particle positions. These are held constant through-
out any elementary time interval. It is evident that the
magnitude of the time interval At must be restricted, as
otherwise a particle could pass through many cells in a cycle
and well out of the region for which its velocity was inter-
polated. This could result in large inaccuracies and even
instabilities in the solution. Limiting the time step so that
no particle moves more than four-tenths of a cell within any
one cycle has been used as an empirical rule to avoid this
problem.
Every particle is advanced each cycle to a new posi-
tion using Eq. (1.14). Thus, the particle traces out in time
a trajectory for the pollutant mass. Movies of these tra-
jectories can be made showing the flow of pollutants.
Boundary conditions are introduced by modifications of
the fictitious total velocities. Thus, impervious barriers
are simulated by not allowing particles to be transported
14
-------
3SR-844
across these boundaries. Transmittive boundaries are permeable
to particles which continue to pass through them with the total
velocity.
This completes the cycle by which the particle method
solves the diffusion equation for a time step. In each cycle
the cellular concentrations are calculated from particle
masses, the fictitious total velocity for each cell is cal-
culated as the sum of mean wind velocity and the turbulent
flux velocity (which is calculated by finite differencing the
cellular concentrations), and the particle positions are up-
dated using an interpolated total velocity. Time is stepped
along cycle by cycle and the position of each particle traces
a trajectory of the pollutant mass.
1.3 SIMULATION OF EMISSIONS AND CHEMICAL REACTIONS
As pollutants are emitted into the atmosphere, addi-
tional pollutant mass must be added to the simulation. The
additional pollutant mass may be added to the mass associated
with a pre-existing particle at the source location, or a new
particle may be generated. The new particle would enter the
simulation at the source location with pollutant mass equal
to that physically generated during the time interval. Pol-
lutant generation is treated separately for each cell of the
grid. Point, line, area and volume sources can easily be
simulated.
The simulation of the chemical reactions begins with
a chemical mechanism — a set of chemical reactions and their
reaction rate constants — assumed to be an adequate descrip-
tion of the relevant pollutant reactions. The chemical mecha-
nisms can be translated into a set of chemical kinetics rate
equations — coupled first-order differential equations de-
scribing the time rate of change of each reactant and product
15
-------
3SR-844
pollutant species. This will be discussed in detail in Sec-
tion 4.2.
The solutions of the chemical kinetics equations
determine the cell concentrations but the changes in concen-
trations have to be related back to particle masses . Within
each cell the particle masses are re-weighted by the frac-
tional change in concentration due to chemical reactions,
m(t + At) - m(t)
For each chemical species and each cell the fraction is cal-
culated. If a pollutant species were first generated in a
cell by chemical reactions, there would be no pre-existing
particle to re-weight and a new particle would be introduced
into the cell.
1.4 ONE -DIMENSIONAL EXAMPLE OF THE PICK METHOD
In this section a solution of a diffusion equation by
the PICK method will be developed and its accuracy analyzed
by comparison with a known exact analytic solution. For this
purpose the diffusion equation to be solved has been simpli-
fied to describe one-dimensional constant advection, i.e., with
constant and uniform velocity u and Fickian diffusion (K =
JC
constant K) , without sources or boundaries, but including a
removal rate proportional to the local concentration, as given
by;
it * *
This removal rate can be thought of as simulating a chemical
reaction of the pollutant with a major atmosphere species
16
-------
3SR-844
whose concentration is essentially unchanged by the reaction.
The exact analytic solution, for constant u and K, and for an
initial concentration distribution that is Gaussian centered
at x = 0, of mean square deviation a2 and total mass Q, is
o
given by
c(x,t) =
/rr(2a2 + 4Kt)
exp
(1.17)
In the PICK method, the initial concentration distribution is
represented by a distribution of particles. For simplicity
we will use particles of equal mass. Each of the N particles
has an initial pollutant mass m(0) = Q/N where N must be
a large number for an accurate simulation.
The one-dimensional space is divided into cells of
size Ax , with the cells numbered as in Figure 3.
x/Ax
Figure 3 — Gaussian Distribution with the Approximate
Cell Concentrations.
17
-------
3SR-844
Thus, initially at t = 0, the number of particles in cell i
centered about x. = iAx is
x-+Ax/2
exp -
20:
dx
(1.18)
which is a histogram of the Gaussian distribution. Note that
n. is an integer approximation to the integral the accuracy
of which increases with N . If Ax is sufficiently small
this can obviously be approximated by
n (0)
20 2
0
(1.19)
which, of course, is also a Gaussian distribution. The steps
taken in the PIC method to simulate advection, diffusion and
chemical reaction will be described below using this initial
distribution of particles in cells.
First the mean turbulent flux velocity is calculated
for each cell. Following Eq. (1.12) this is given initially
by
u
fi
2Ax
- ni
n.
(1.20)
where the particle mass mQ has cancelled in the ratio of the
concentrations. The number of particles per cell n. can be
substituted from Eq. (1.19) to give
18
-------
3SR-844
U£i =
2Ax 1
<>i - *u>/2a: c*i
- e
x?
1-1
, (1.21)
Hence, by power series expansion of the two exponentials,
assuming (x2 + i ~ x|) and (x* ~ x^_i) are both small
compared to a2 ,
o
Ufi *
- JL!
(x? - x2 ) - (x2 - x2 )
v ; v -'
2a2
o
o K
- 2Ax
lxi+l " xi-lJLXi+l xi-lj
2a2
0
Kx.
~ ~
(1.22)
The cell-centered fictitious total velocity is
u- = u +
Kxi
(1.23)
This is the fictitious total velocity that appears in Eq. (1.6)
and is non-solenoidal. A particle at an initial position, say
x , between x. and x-+i will have a weighted velocity
given by the one-dimensional form of Eq. (1.13), namely
" x
X "
Ax
(1.24)
19
-------
3SR-844
or by simplifying
u(x) = u + §£ . (1.24b)
o
We note in the above that although the real fluid velocity u
is here postulated to be constant and uniform (and, hence,
solenoidal so that 9u/9x = 0), both the turbulent flux veloc-
ity and fictitious total velocity are time- and space-depen-
dent. For t = 0 , both increase linearly with distance from
the origin. In a time interval At the position of the
particle that is initially at x will be changed by u(xJAt
and the new position x1 will be given by
t (1-25)
so that
- UAt
x =
KAt
a^~~
0
It follows that the n.(At) particles lying in the i-th cell
at the end of the time interval At , i.e., whose positions lie
A -yr A v
between x' = x. - ~— and x1 = x. + •=— at time At , were
1 /r 1 L*
initially located between the positions
(x. - **) - uAt (x. + **
\ i 2 / \i 2
- uAt
and
_
. KAt KAt
~ 1 + ~
respectively, that is for a total distance, say A£ , centered
about the position (x. uAt)/(l + KAt/a2) and given by
1 o
20
-------
3SR-844
Ax
KAt
< Ax
(1.27)
But, initially, the number of particles centered in the line
segment of length Ax centered on position x. is given by
Eq. (1.19), from which it follows that initially the number
of particles contained in a segment of length A£ centered on
the position (x^ - uAt)/(l + KAt/a2) is given by
. . A j_ \ 2
i - uAt)
- uAt)
NAA
2a2 (1 + KAt/a2)2
27T02
0
KAt
(1.28)
It follows from conservation of particles that
K (At) =
NA&
2
- uAt)
(2cr2 + 4 KAt)
o
or by using Eq. (1.27) that
- uAt)
1
1 + KAt
a2
0
NAx
NAx c (2°l + 4 KAt^
/2Tta2
0
(x . - uAt) 2
2a2 + 4 KAt
0
A(2a2 + 4 KAt)
provided KAt/a2 is small.
21
(1.29)
-------
3SR-844
So far the effects of advection and diffusion in
Eq. (1.16) have been considered. To include also the effect
of the removal rate term simulating a chemical reaction, we
note that this by itself produces an exponential decay of con-
centration so that
c(t + At) = c(t) e"XAt . (1.30)
However, the changes in cell concentration have to be related
to particle masses. Within each cell the particle masses are
reweighted by the fractional change in concentration due to
the chemical reaction, as specified in general by Eq . (1.15)
and in particular by
m(At) = m(0) = m(0) e" , (1.31)
Finally, the cell concentration is given as the product of
the number of particles in the cell and particle mass divided
by the cell volume, so that
c(xi,At) = nt(At) m(At)/Ax
and hence for the above example
(xi - uAt)
c(x. ,At) = V exp
1 /Tr(2a2 + 4 KAt)
2a* + 4 KAt
- AAt
0, 32)
where Q = NmQ has been used. This is in agreement with the
analytic solution as given by Eq. (1,17).
22
-------
3SR-844
1.5 ACCURACY OF THE PICK METHOD
The PICK method is a hybrid of Eulerian and Lagrangian
techniques and has some of the advantages and disadvantages of
each. The Eulerian nature requires an assumption of homogeneity
within each fixed grid cell. The Lagrangian nature is intro-
duced by quantization of the pollutant mass (a continuous
variable) into particles (a discrete variable). The errors
introduced by the quantization into particles has been evalu-
ow-e.
ated operationally and >5" discussed in Section 2. The more
particles that are used, the less mass (on the average) is
represented by each particle, and the greater the variation
and resolution of concentration that can be simulated. A mini-
mum number of particles exists below which resolution and
accuracy are severely degraded while the maximum number of
particles is usually determined by computer capacity.
The Eulerian grid size determines the maximum spatial
resolution of the solution of the atmospheric diffusion equa-
tion. Though distributions smaller than grid size are visible
in the particle positions, only cell averaged concentrations
are solutions of the equation. Smaller than grid size phe-
nomena (for example, point emissions) can be represented by
particles, and the accuracy of such simulations has been
demonstrated by tests reported in Section 2.
The accuracy of advection simulation is limited by
the Eulerian grid size and the velocity interpolation prescrip-
tion. Mean wind patterns having a sub-grid scale, such as
smaller atmospheric eddies, can only be included in an averaged
sense as turbulent diffusion or by increasing the resolution.
If the wind field for advection at time t were known suffi-
ciently, the particle position could be expanded in a Taylor
series; for example in one dimension:
23
-------
3SR-844
x(t + At) = x(t) + u(x(t))At + hu^± At2 + ... , (1.33)
where u = 9x/9t and 9u/9t = u(9u/9x) have been used. The
first two terms on the right correspond to Eq. (1.14) (in
which u = u(x(t))) and the error to lowest order is given by
the last term, %u(9u/8x) At2 . Thus, in order to limit errors,
a limitation on the time interval is required. The time limit
employed has been chosen arbitrarily and empirically, and re-
stricts the maximum change in parcel position to four-tenths
of the cell dimension. Thus, the maximum relative error in
particle transport is
9u
u
_
9x 2 _ - 2u 9x _ Ax 9u
uAt 0.4Ax < u 9x '
As an upper limit 9u/9x <_ 2u/Ax , then this relative error is
bounded by 0 . 4 . In actual practice, this error is less than
10 percent.
Additional inaccuracies are introduced in the simula-
tion of diffusion — from the finite differencing used in re-
placing the diffusion term by the turbulent flux velocity,
Eq. (1.2) and Eq . (1.12). The approximation to the logrithmic
derivative in Eq. (1.12) can be expanded in a Taylor series to
show the inaccuracies incurred in the prescription in Eq . (1.12)
as follows:
24
-------
3SR-844
c * c
— - — - - 1- = . . \( p + r ' Av + ir"AY2
2Ax c 2Ax c (Lci i x 2i x
- '
- (ci - cjAx + c^'Ax - cV'Ax3 ...)
c.' c1." Ax2
— + -±z - + 0(Ax")
ci 6ci
That is, the lowest order error term in this approximation is
(T E cTx"5" ^x2 ' Thi-s has been acceptable in all applications
attempted to date since it is a small error in the diffusion
term, which itself is small.
Errors introduced by the simulation of reactions will
be discussed in Section 4.
1.6 SUMMARY OF THE PICK METHOD
In summary, the PICK method uses Lagrangian mass points
to simulate transport by advection and diffusion as described
by the differential equation for turbulent atmospheric diffu-
sion. An Eulerian grid is used for winds, emissions and values
of the concentration. The velocities used for the transport
are fictitious, non-solenoidal velocities which, when multiplied
by the local concentration, have a divergence equal to the
turbulent diffusion and advection in the original equation.
Boundary conditions of the differential equation become con-
ditions on the particle motion and are introduced through
modifications of the fictitious velocity field.
By use of the PICK method, three-dimensional solutions
of the atmospheric diffusion with temporally and spatially
25
-------
3SR-844
dependent winds and diffusivities are practical. The La-
grangian particles eliminate the problems of artificially in-
troduced dispersion that are associated with finite differ-
ence solutions of the differential equation. In applications,
air pollution can even be simulated for zero wind conditions.
The PICK method also lends itself to a new output form
(movies of pollution as represented by particle positions)
in addition to the conventional forms.
26
-------
3SR-844
2. TESTING THE PICK METHOD
2.1 DESCRIPTION OF THE PICFIC CODE
In order to test the concepts discussed in Section 1,
the PICFIC (Particle-I_n-Cell with FI_£kian diffusion) code was
developed (actually PICFIC treats non-constant diffusivities
and so is not restricted to Fickian diffusion). The equation
being solved is the two-dimensional analogue of Eq. (1.1) (for
convenience the z-direction was dropped). Boundary conditions
for the edges of the grid were also ignored and all simula-
tions were ended before any pollutants reached the grid edge.
These simplifications were used because PICFIC was developed
to test the advection and diffusion techniques proposed for
the PICK method. The results of these tests will be discussed
in Sections 2.2 and 2.3. Some additional simulations using
PICFIC will be presented in Section 2.4.
A .simplified flow chart of PICFIC is given in Figure 4.
A detailed computer-generated flow chart is given in Volume II,
Section 1, along with a listing of the coding.
PICFIC consists of three routines — the main program,
CALCP and CELL. The main program sets up the problem. CALCP
is the subroutine in which new particle positions are calcu-
lated using the total velocity field. CELL is an output sub-
routine for printing the cell concentrations.
In the main program, the 2-D velocity components U^.
and V.. and diffusivities K .. and K .. are specified
27
-------
3SR-844
MAIN PROGRAM
CALCP
SET WINDS U
DIFFUSIVITIES K
CALCULATE
PARTICLE POSITIONS
FOR INITIAL
DISTRIBUTION
(Transports
Particles)
CELL
(Prints Cell Con-
centrations)
CALCULATE
CELL CONCENTRATIONS
KIBAR
CALCULATE
TURBULENT FLUX
VELOCITIES V
TRANSPORT
PARTICLES BY
AREA WEIGHTED
TOTAL VELOCITY
CELL
INTEGERIZE
CELL CONCENTRATIONS
IKIB •= KIBAR
PRINT CELL
CONCENTRATIONS
Figure 4 -Macro Flow Chart of the PICFIC Code,
28
-------
3SR-844
(computer representations used are: u.. = U(I,J,1) , v.. =
U(I,J,2) , Kxij = D(I,J,1) and Kyij = D(I,J,2) . For the
problems used in testing PICFIC these values were initialized
and not changed during the simulation. The test problems used
a Gaussian distribution for the initial distribution of the
numbers of particles per cell. The UNIVAC 1108 operating sys-
tem available at S3 has a generator for normally distributed
random numbers. This was used to position particles so that
the cell concentration distribution was Gaussian.
The heart of PICFIC is CALCP which transports the
particles according to the PICK method. The PICK method re-
quires the calculation of the turbulent flux velocities, from
Eq. (1,2), which in turn use cell concentrations. As men-
tioned in Section 1, two prescriptions have been used, the
first of which treats particles as mass points, the second
uses an area-weighting procedure. In the first method, the
cell concentration is calculated by summing the mass associ-
ated with particles in the cell. With the area-weighting
procedure, each particle mass is divided among the four clos-
est cells, in the following steps:
(1) The particle at position (X,Y) is determined
to be in cell (N,M) by
N = TRUNCATED INTEGER(X/DX)
(2.1)
M = TRUNCATED INTEGER(Y/DY)
(2) The fractional distance of the particle from the
center of cell (N,M) is calculated from
29
-------
3SR-844
FX = X/DX - w+%)
(2.2)
FY = Y/DY - (M+%) .
(3) If these fractions are positive the particle is
in the forward portion of the cell and the other closest
cells have NN = N+l or MM = M+l indices; if negative the
indices are NN = N-l or MM = M-l .
(4) The particle mass is divided among the four cells
by the fractional overlap such that
KIBAR(N,M) = KIBAR(N,M) + (1-FX) • (1-FY) -KI
KIBAR(NN,M) = KIBAR(NN,M) + FX-- (1
(2,3)
KIBAR(N,MM) = KIBAR(N,MM) + (l-FX)'FY-KI
KIBAR(NN,MM) = KIBAR(NN,MM) + FX-FY'KI ,
where KIBAR is the storage location of the cell concentration
and KI contains the particle mass.
Once the cell concentrations have been calculated, the
difference algorithm, Eq . (1.12), is used to compute the cell
centered turbulent flux velocity for each cell. The last step
in CALCP is the moving of each particle. The particle's frac-
tional overlap with the four closest cells is determined as in
the above area-weighting, Steps 1-3. These are the fractions
needed in the area-weighting for interpolation of cell veloci-
ties to the particle position,
VX = (1-FX)«(1-FY).U(N,M,1) + (1-FX) • (FY*U(N,MM , 1)
FX «FY»U(NN ,MM,1)
(2.4a)
30
-------
3SR-844
VY = (1-FX)»(1-FY)«U(N,M,2) + (1-FX)•(FY-U(N,MM,1)
+ FX=(1-FX)«U(NN,M,1) + FX'FY-U(NN,MM,2) .
(2.4b)
The interpolated X and Y components of velocity are multi-
plied by the time step, DT , and added to the particle coordi-
nates :
X = X + VX«DT
C2.5)
Y = Y + VY-DT .
This completes the particle transport cycle as coded in CALCP.
The PICFIC code requires the following computer storage
(1) For each cell of the grid -- 7 variables;
(2) For each particle -- 3 variables .
In addition, 3,000 variables are used for storing the normally
distributed random numbers and a few hundred for miscellaneous
storage. For example, as compiled on the S3 UNIVAC 1108, the
PICFIC program size was 20,000 words, of which 3,000 were
particle variables, 9,500 cell variables and 3,000 random
numbers.
A fully commented listing of PICFIC, a computer-
generated flow chart and a sample problem output are given in
Volume II, Section 1. The tests made with PICFIC used the
code as described and compared the calculational results with
theoretical or analytic solutions.
31
-------
3SR-844
2.2 COMPARISON OF METHODS FOR PURE ADVECTION
The accuracy of the PICK method for advection was
analyzed using an off-axis wind field. The wind field was
chosen off-axis because neither PICFIC nor the finite differ-
ence methods used for comparison to PICFIC are rotationally
invariant. This corresponds to solving the 2-D advection
equation:
— + u— + v— =0 (2 61
9t 9x 9y l^toj
where c is the pollutant concentration, u is the x-compo-
nent of velocity (u = 10 meters/second), and v is the y-
component of velocity (v = 5 meters/second). The initial
tribution of pollutant concentrations is given by the Gaus-
sian:
c(x,y,0) = 711 exp -[(x-9)2 + (y-9) 2]/12 . 5 J , (2,7)
in which a = a = 2.5 Km. The cell averaged concentration
x y
of the center cell, cell (9,9) is 69. The exact solution at
time t is
c(x,y,t) = 711 exp -f(x-9-ut)2 + (y-9-vt) 2]/12. 5 , (2.8)
The calculation using the PICFIC code had a 30 x 50 grid with
each cell one kilometer on a side. The size of the grid was
selected to be large enough so that the deviations from the
exact solution would be seen and small enough that computer
costs would be kept to a minimum.
The initial distribution was constrained within the
grid boundaries by setting to zero cell concentrations over
32
-------
3SR-844
eight cells from the center of the distribution. Each of the
256 cells with a non-zero concentration had a single particle.
The results of this calculation are shown as a contour plot
of cell concentrations in Figure 5. The initial distribution
is shown in the first diagram for T = 0 . The contour
values correspond to 0 .1 , 10 , 20 , 40 , and 60. The squaring-
off of the outermost contour, the 0.1 value, is due to the
truncation of the initial distribution, i.e., cells over
eight cells away from the center of the distribution were set
to zero concentration.
The results in Figure 5 show that the PICK technique
can simulate advection to the accuracy of the interpolated
wind field or, in this example with a constant wind field, to
the accuracy of the computer. The accuracy of the PICK tech-
nique for advection is inherent in the Lagrangian nature of
the particles .
Similar advection tests were made using a first-order,
7 8
Fromm's fourth-order and Crowley's fourth-order schemes.
Each of these schemes uses a splitting technique to combine
g
integrations in each direction. The first-order scheme
uses the algorithm
= n uAt , n _ n v
ci Ax l i-l,j ijj
n+1 _ . n+J* vAt
-
for velocities u~ , v~ > 0 , and with superscripts referring to
time (t = nAt) and subscripts for position (x = iAx ,
y = jAy) . This "upstream" differencing originated with
Lelevier. Mahoney and Egan refer to this technique as a
forward time and backward space differencing scheme and report
33
-------
3SR-844
Tiiue = 0 seconds
Time = 1200 seconds
Time = 1800 seconds
Figure 5 - PICFIC Advection Simulation Test with
a = a = 2.5 Km, Contours at 0.1, 10,
y 20, 40 and 60.
34
-------
3SR-844
on an artificial diffusion from the differencing truncation
error having a diffusion coefficient
K =
uAt]
Ax J
(2,10)
for one dimension.
A criteria for time step limitation can be based on
Eq. (2.10). If the diffusivity were to be negative, spatial
gradients would be amplified each time cycle and the solution
would soon become meaningless. To prevent this instability,
uAt
At is chosen such that
< 1
In two dimensions the
same criteria must be met for each direction separately
uAt
Ax
< 1 and
vAt
Ay
< 1
If splitting were not used in an
explicit scheme (for example, as reported by Pandolfo
the restriction is
12
), then
uAt
Ax
vAt
Ay
< 1
With splitting the truncation diffusion is independent in
each direction and the same expression holds for two and three
dimensions, i.e.:
uAtl
Ax~J
vAtl
Ay J
C2.ll)
for two dimensions. The effects of the truncation diffusivity
in the first order scheme can be seen by the spreading of the
contours in Figure 6. The flattening of the outermost contour
35
-------
3SR-844
Time = 0 seconds
iiii
I!.!
I 11
111
Time « 1200 seconds
Time = 1800 seconds
Figure 6 — First-Order Finite Difference Advection Test
with a = a = 2.5 Km.
x y
36
-------
3SR-844
near the top of the grid is due to the distribution running
off the edge of the grid.
The diffusivity calculated from the results pictured
in Figure 6 was
KX = 2900 m2/sec K = 1900 m2/sec , (2,12)
whereas the expression for the truncation diffusivity,
Eq. (2 .11) , gives
Kx = 3000 m2/sec K = 2000 m2/sec (2.13)
with u = 10 m/sec , v~ = 5 m/sec and Ax = Ay = 1000 m ,
At = 40 seconds. These can be contrasted with atmospheric
diffusivities that have been measured in the range of 0.01
to 1000 m2/sec. Thus, this first-order donor cell scheme
would not be useful for simulating pollutant transport with
these winds and cell size because the truncation diffusion
would be greater than the atmospheric diffusion and would mask
the real atmospheric diffusion.
The truncation diffusion originates from errors made in
the truncation of the finite difference approximation of the
derivatives. With the first order method, the advection equa-
tion (in one dimension) is differenced
n+1 n n n
Ci " Ci = ~uCi ^i'1 . (2,14)
This can be expanded by a Taylor series about time n and
position i , for example
37
-------
3SR-844
3_\n a2_\n A 2
n n 3c\ *„ j. o c \ AX
c. , = c. -
to give the expression
,2^11 Av2
+ .../At =- |£) Ax.0r^f.../Ax
) ( 3x/i 9x /i Z )
C2.16)
The second derivatives can be related by differentiating the
advection equation
32c = _,,3 /3c
When this is substituted back into the expansion, the advec-
tion equation is recovered to first order and the coefficient
of the second order term is the truncation diffusivity; drop-
ping subscripts,
3c _ 3c . uAx/ uAt\ 32
3t - -
One way of reducing the truncation error is to use a
higker order difference schemet In any finite difference
scheme, the terms dropped can cause a dispersion of an initial
distribution during each calculational time step. Higher
order schemes in general increase the accuracy of the calcula-
tion and reduce the calculational dispersion. On the other
hand, higher order schemes give overshoot errors in the pres-
ence of large spatial gradients and require more computations
each time step. The same difficulty with large gradients is
encountered in polynomial curve fitting where a high order
38
-------
3SR-844
polynomial may introduce large fluctuations in approximating a
rapidly changing function.
Two fourth-order methods were compared to the PICFIC
calculation (the term "fourth-order" comes from the fact that
all terms up to fourth order are kept in the finite difference
expansion of the time and space derivatives). Fromm's fourth-
7
order method is a space-centered flux conserving formulation;
it can be expressed as follows:
:u *
cV.' - c". *
C2,19}
where
'24
C2.20)
and
a = vAt/Ay
Fn+is is obtained at time n+% using similar displacements
i-^J
in the i subscript, with a = uAt/Ax , This method was used
39
-------
3SR-844
to calculate the results shown in Figure 7. The slight dis-
tortion of the outer contour is of minor importance since this
is the 0.1 contour value. The major change between the first
order scheme shown in Figure 6 and the Fromm fourth-order
scheme is that the distribution remains peaked with the fourth
order scheme but is dispersed in the first order scheme as
shown by the maximum value at the center of the calculated dis
tributions in Table I. An examination of the calculated dis-
tribution revealed that
(1) the a2 did not spread as a linear func-
tion, as would be the case for a diffusion
error, and
(2) errors in a2 were less than 2 x 10 Km .
As mentioned previously, a higher order method would be ex-
pected to have increased accuracy over the first-order scheme
and to eliminate the problem of truncation diffusion in the
solution of the atmospheric diffusion equation.
TABLE I
CALCULATED CENTRAL VALUES AFTER 45 CYCLES
FOR DIFFERENT METHODS (Correct Result: 69)
A.X = Ay = 1
INITIAL VALUE AT CENTER
DISTRIBUTION METHOD OF DISTRIBUTION
a = 2.5 PICFIC 69
First-Order 29
Fromm Fourth-Order 68
a = 1.0 PICFIC 69
Fromm Fourth-Order 38
Crowley Fourth-Order 37
40
-------
3SR-844
Time = 0 seconds
Time = 1200 seconds
Time = 1800 seconds
Figure 7 - Fromm Fourth-Order Method Advection Test
= 2.5 Km.
with a =
J\.
y
41
-------
3SR-844
A more sharply peaked distribution was also used to
compare Fromm's method with the PICK technique. For an ini-
tial Gaussian distribution with a = 1 Km (normalized to the
center cell concentration = 69) , PICFIC was used to calculate
the diagrams of Figure 8. Again, using the PICK technique,
there is no distortion. The same initial condition using
Fromm's fourth-order method resulted in distortion of the dis
tribution as shown in Figure 9. Again, the distortion of the
outer contour (0.01) is not significant except to note that
between these contours are negative values of concentration
which are obviously unphysical. However, the major signifi-
cance of this figure is the smearing of the peak values at
center of the distribution. The center value is only 38 in-
stead of the correct 69.
Q
Crowley's fourth-order method was also tested with
the sharply peaked distribution. In Crowley's method, the
fluxes used in Eq . (2.19) are:
C2.22)
with notations and definitions as in, Fromm's method. The re-
sults obtained with Crowley's fourth-order method are shown in
Figure 10 for the a = 1 Km initial distribution. These re-
sults are essentially the same as those from Fromm's method
shown in Figure 9 .
42
-------
3SR-844
Time = 0 seconds
'^
/H
j
i
t
i
j !
i
1
1
1
I
!
|
•
t
i
i :
1 ;
1 ;
1 •
i |
1 !•
|i
j j
• t
.
r .;
1 j
• 1
•
•
r
i
1
i
• i
i
i
I I
I i
Time = 1200 seconds
Time = 1800 seconds
Figure 8 - PICFIC Advection Test with ax = a = 1 Km,
43
-------
3SR-844
Time = 0 seconds
ii
I
i i j ! ! i ' :
S ' ! ' i
i; : •
Time = 1200 seconds
Time = 1800 seconds
Figure 9 - Fromm Fourth-Order Method Advection Test
with a = a = 1 Km.
44
-------
I
3SR-844
Time = 0 seconds
mil
11
I !
1:1!
! i! !
! i
i
Time = 1200 seconds
Time = 1800 seconds
Figure 10 — Crowley Fourth-Order Method Advection Test
with a
Km,
45
-------
3SR-844
Both fourth-order methods show an inability to treat
a sharply peaked distribution. The cause of the error is de-
picted in Figure 11. In air pollution studies, such sharply
peaked distributions could result from large single sources.
Using a high order difference technique to eliminate the trun-
cation diffusion in an air pollution model would then require
special treatment of each large emission source. The PICK
technique, on the other hand, permits simulation of advection
irrespective of the nature of the distribution of concentration,
Computer time studies of each method were also made.
The computer used for these calculations was a UNIVAC 1108
operating in a multi-processing mode with the EXEC 8 system.
It was found for the illustrated calculations that the first-
order method took 15 seconds, both fourth-order methods re-
quired 19 seconds, while the PICFIC calculation was completed
in 12 seconds. The PICFIC calculation for which the results
are plotted used 200 particles initially placed at cell
centers within the initial distribution, Another calcula-
tion was performed with the PICFIC code in which 1,000
particles were placed at random within the initial distribu-
tion. The corresponding computational time increased to 18
seconds. The results of the PICFIC calculations for pure
advection were independent of the number of particles used.
2.3 EVALUATION OF DIFFUSION SIMULATION
The second test of the PICK technique addressed the
question of how many particles are needed for a given resolu-
tion. As mentioned above, the PICK technique accurately
simulates advection with as few as one particle per cell for
each cell having a non-zero concentration. On the other hand,
with the simulation of diffusion the velocity field is non-
solenoidal and requires that particles move apart or disperse,
46
-------
3SR-844
Figure
11 -
to the
Curve fitting withj
a Caussian curve and^
the errors in the fit,
47
-------
3SR-844
even within a single cell. Diffusion can be calculated only
in a collective sense, and a large number of particles are
needed for accuracy.
A test problem was generated for evaluating the accu-
racy of the PICFIC code in simulating diffusion. The grid
consisted of 50 by 27, 1 kilometer square cells. A wind
parallel to the (x) axis having magnitude u = 10 meters/
second and a symmetric diffusivity constant KX = K = 103
meters2/second were used. The wind field parallel to the
axis was chosen so that diffusion both, parallel and perpen-
dicular to advection could be studied, i.e., in the y-direc-
tion there is no wind so there is only pollutant transport by
diffusion, and in the x-direction advection is much greater
than diffusion so there is primarily advective transport. The
initial pollutant distribution was Gaussian with a = a =
x y
1 Km. The solution was calculated to t = 3200 seconds and
the mean position and standard deviation in both directions of
the distribution were determined. Analytically, these are
given by
x = x + ut y = y + vt
(2,23)
a2 = a2 + 2K t a2 = a2 + 2K t .
x o x y o y
The accuracy of the calculation was evaluated by Comparing
computed PICFIC results with the analytic formulas, Eq. (2.23).
"Effective" velocities and "effective" diffusivities (corre-
sponding to a linear least squares fit of the positions and
standard deviations at each time cycle) were calculated and
compared to the correct values (u = 10 m/sec, v = 0,
Kx = Ky = 10 m2/sec). Table II shows these values for four
calculations corresponding to 1,000, 500, 100, and 10 particles,
all using the area weighting procedure in the calculation of
cell concentrations.
48
-------
TABLE II
PICFIC DIFFUSION EVALUATION VARYING THE
NUMBER OF PARTICLES PER CELL USING
THE AREA-WEIGHTING PROCEDURE
CORRECT
1000
u 10.0 m/sec 10.0
v 0.0 m/sec -2,9 x 10"7
K 1.0 x 1Q3 m2/sec 0.93 x 103
JC
K 1.0 x lo3 m2/sec 0.93 x 103
EFFECTIVE
Number of
500
10.0
3.0 x 10"
0.92 x 10
0.91 x 10
Particles Used
10.0
10.0
7 -1.4 x lo-7
3 0.81 x 103
3 0.80 x 103
10
10.0
-3.6 x 10'5
0.42 x 103
0.30 x 103
Computational time
(seconds)
37
23
11
oo
-------
3SR-844
These results once again confirm the ability of the
PICK technique to simulate advection accurately. Even with
as few as 10 particles the effective velocities are the same
as the imposed wind velocity. The effective diffusivity does
depend on the number of particles used in the calculation.
The results show 93 percent agreement with 1,000 particles.
With 100 particles, the agreement has dropped to 80 percent.
There is only qualitative agreement in diffusivities with 10
particles. The correct final pollutant distribution has
a2 = 7.4 Km2 . Thus, 90 percent of the pollutant mass
(corresponding to 2a in each direction) is spread over a
90 Km2 circle or the corresponding 100 cells. For 10 or even
100 particles, this distribution requires the particles to be
spread thin, and, as a result, averaged processes may be
poorly represented. Significantly, the effective diffusivity
was less accurate towards the end of the calculation, when the
distribution was the most diffuse. In the case of 100
particles, the effective diffusivity dropped to 60 percent of
the correct value after the distribution was spread over 30
cells. This suggests that a lower limit of approximately
three particles per cell is needed for accuracy. This agrees
well with experience gained with PIC (the hydrodynamics code
using the particle-in-cell method with no diffusion) where 10
particles per cell are recommended.
Additional information was gained from these test cal-
culations ; an estimate of the speed of PICFIC and of the dis-
tribution of time between particle and cell related computa-
tions, Every time step requires 3-5 x 1Q~5 seconds to process
each cell in each dimension, and 1-3 x 10"1* seconds for each.
particle computation, On the S3 UNIVAC 1108, this corresponds
to 10-15 operations per cell per cycle per dimension and 40-100
operations per particle per cycle per dimension. Computer
storage requirements were described in Section 2.1,
50
-------
3SR-844
The alternate method of calculating cell concentrations,
using the mass point procedure, was also tested. It has ad-
vantages for the simulation of chemical reactions (as will be
discussed in Section 4). The area-weighting procedure pro-
vided a smoothing in concentrations , but by not using area-
weighting for concentrations, a computational time savings of
14 percent in the particle related computations was achieved
for the PICFIC code. Area-weighting is still used to inter-
polate the velocity field in the calculation of particle
velocities. This modification resulted in slightly poorer
statistics for the diffusion simulation as shown in Table III
for 1,000 particles and a 50 x 27 grid.
TABLE III
PICFIC DIFFUSION EVALUATION VARYING THE
CELL SIZE USING 1000 MASS POINT PARTICLES
u
V
K
X
K
y
CORRECT
25 x 14
10.0 m/sec 10.0
0.0 m/sec -6.8 x 10"6
1.0 x io3 m2/sec 1.4 x 103
1.0 x io3 m2/sec 0.85 x io3
EFFECTIVE
Grid Size
50 x 27 100 x 54
10.0 10.0
-6.5 x 1Q-6 -5.4 x 1Q-6
0.90 x IO3 0.76 x io3
0 .87 x io3 0.73 x io3
Computational time
(seconds)
14
32
107
51
-------
3SR-844
A set of computations was made to determine the effects
on diffusion simulation accuracy of doubling and halving the
cell size. The PICFIC code without area-weighting (using the
mass point procedure) was used for these test calculations.
These results are shown in Table III. The accuracy of the
advection simulation was unchanged by the grid size variation.
In the case of double-sized cells (grid size 25 x 14), the
accuracy of the simulated diffusion was degraded by lack of
resolution. Physically, this case would correspond to a
pollutant distribution localized to a variance of one-fourth
of a cell size. When the resolution was increased by halving
the cell size (grid size 100 x 54)t the simulated diffusion
was degraded. This degradation was due to an insufficient num-
ber of particles per cell to define the distribution accurately
as it spreads. Towards the end of this calculation, the dis-
tribution is spread over 300 cells. This corresponds to about
three particles per cell; with so few particles per cell, the
diffusion calculation becomes inaccurate as was discussed
previously. During the first 1,000 seconds of the calculation,
the effective diffusivities averaged 0.84 x lo3 meters/second
in the x and y directions, respectively.
These tests have shown the accuracy of the PICK method
by comparing the first (mean position) and second (mean square
deviation) moments of the concentration distribution calculated
by PICFIC with their correct values. Different distributions
may have the same first and second moments. To demonstrate
that the PICK method accurately simulates diffusion preserving
the initial Gaussian shape of the distribution, normalized x
and y distribution plots are shown in Figures 12 and 13. In
Figure 12, the initial x and y cell^, concentration distri-
butions used in the test are shown along with a normal curve.
After 41 cycles of PICFIC (with both diffusion and advection)
52
-------
3SR-844
-3
y/Ay
Figure 12 — Histogram of Initial Cell Concentrations Compared
to Normalized Gaussian Distribution.
53
-------
/n
A
A
\
\
\
2 -1
x/Ax
3SR-844
I
•3
-2 -1
I
1
y/Ay
Figure 13 - Renormalized 4lit Cycle of PICFIC Calculation
Advection and Diffusion] Compared to Normalized
Gaussian Distribution,
54
-------
3SR-844
the calculated distributions are renormalized (to c(0) = 1
and 0=1) and plotted with a normal curve in Figure 13.
The shape of the distribution is properly maintained by the
PICK method calculation. From these plots, it can be seen
that the accuracy after 41 cycles is as good as the initial
distribution.
2.4 PICFIC SIMULATIONS
In this section, a number of additional simulations using
the PICFIC code will be presented. First, for completeness, the
original PICFIC simulations are included even though they have
been reported elsewhere. These original simulations were made
to show the capabilities of the PICK method to reproduce the ef-
fects found in the atmosphere — anisotropic diffusion, inversion
layer and wind shear.
The PICFIC code permits two forms of graphically pre-
senting calculated data: plots of particle positions and con-
tour plots. Particle plots represent pollutant concentra-
tions by density of particles and thus are useful for qualita-
tive assessment of large amount of data. The contour plots
are computer-generated and show the usual angularity due to
computer contouring techniques. The calculational grid is not
shown but consists of 1,250 cells (25 * 50),
Figure 14 illustrates a single puff of pollutant
emitted on the left and carried downwind (to the right), and
displayed at three instants of time. Both types of data plot-
ting are shown for comparison. In this case no diffusion was
permitted and Figure 14 shows that the PICK method intro-
duces no artificial diffusion. When a constant diffusivity
is assumed, the usual Gaussian puff results as shown in
Figure 15. Figure 16 illustrates the effects of anisotropic
diffusion. The effect of an inversion at the emission height
55
-------
3SR-844
"
Figure 14a — Pollutant cloud advection in constant wind field
with no diffusion. Three time instants are dis-
played together for comparison.
Figure 14b —
Same as Figure 14a. The contours are for integer
powers of 2 in pollutant concentration.
56
-------
3SR-844
r«..:.
Figure 15a — Pollutant cloud advection in constant wind field
with, spherically symmetric diffusion,
Figure 15b —
Same as Figure 15a.
of two,
Contours at integer powers
57
-------
3SR-844
Figure 16a — Pollutant cloud advection in constant wind field
with aspLerical constant diffusion.
Figure 16b —
Same as Figure 16a,
of two,
Contours at integer powers
58
-------
3SR-844
is depicted in Figure 17. Lastly, the dispersion caused by
wind shear and by constant diffusivity on a single puff is
shown in Figure 18.
The most complex test case run to date was suggested
1 4
by Taylor: a particulate puff in a shear wind. The speci-
fications of the test were: an initial two-dimensional point
release at 10 meters above ground of particulate matter in a
gravitational field in the negative y-direction, having a
terminal velocity (v = 1 m/sec) in a linear shear wind
(u = y m/sec) and isotropic turbulent diffusivity (K
X
K = 1 m2/sec) . The equation being solved is more general
than the atmospheric diffusion equation, Eq . (1.1), due to
a term for the terminal velocity of the particulate matter,
9c
-v. T^T- . The equation can be written as:
t a t
3t + U8x
9c
3c
(2,24)
This has the same form as the atmospheric diffusion equation
with v" replaced by -v. and so the PICK method can be used
for the solution. The solution is just a special case of a
general analytic solution in two dimensions for particulates
in a linear shear wind reported by Neuringer. PICFIC was
used with a 75 * 16 grid, in which Ax = Ay = 1 meter, and
1,000 particles. The calculation was begun after the puff
had spread to a a > 1 meter, at t = 0.5 seconds. Table IV
shows the results of the calculation compared to the analytic
solution of Neuringer. The 0.5 sec initial conditions for
PICFIC have a 10 percent error in standard deviations due to
the random positions of the particles and the coarseness of
the cell resolution. By 3.0 seconds the puff had started to
approach the edges of the grid and the calculation was stopped,
This was after 75 cycles and so an adequate evaluation was
59
-------
3SR-844
:nfif
V. \'
¥'•
Figure 17a — Pollutant cloud advection in constant wind field
witL aspherical, position^dependent diffusion.
Model of an inversion at th.e source height.
Figure 17b —
Same as Figure 17a.
o f tvro.
Contours at integer powers
60
-------
3SR-844
r •..
Figure 18a — Pollution cloud advection in shear wind field
with, spherically symmetric diffusion, three
separate times are depicted as the cloud moves
to the right.
Figure 18b —
Same as Figure 18a.
o f two.
Contours at integer powers
61
-------
TABLE IV
3SR-844
Time
(sec)
0.5
1.0
2.0
3.0
PICFIC
x y
15.4 10.
19.7 9.
27.7 8.
34.7 7.
(meters) ANALYTIC (meters)
a a x y a a
x y 7 x y
00 1.14 1.11 15.4 10.00 1.04 1.00
50 1.67 1.46 20,0 9.50 1.63 1.41
50 3.06 1.98 28.5 8.50 3.05 2.00
53 4.87 2.34 36.0 7.53 4.86 2.41
possible. In comparison with the analytic solution at 3.0
seconds, the PICFIC puff had not spread as far in the verti-
cal, and perhaps had not responded to the high velocities at
the upper levels and so was retarded in mean x position.
Still these errors are less than the initial 10 percent error.
One final test was made in which the accuracy of repre-
senting a point source by a distribution with a spread of less
than Ax was evaluated. An initial Gaussian distribution with
standard deviation of one-eighth of a cell dimension was used
in this test. The wind and diffusivity used was the same as in
the diffusion tests described in Section 2.3 and the evaluation
was the same — a least-square linear fit to the mean positions
and standard deviations to calculate effective velocities and
effective diffusivities, The effective velocities were
correct to three significant places due to the Lagrangian
nature of the PICK particles. The effective diffusivities were
92 percent and 60 percent of the correct values in the parallel
and cross wind directions, respectively. Thus, while in the
62
-------
3SR-844
PICK method, strictly speaking, only the cell concentrations
are a solution of the atmospheric diffusion equation, the simu-
lation of sub-grid distributions is only slightly degraded.
The implication of this empirical result is that small scale
distributions, as for example due to point source emissions,
can be included in PICK simulations without any special con-
siderations .
2.5 SUMMARY
The two-dimensional PICFIC code has been used to evalu-
ate the accuracy of the PICK method for solution of the equa-
tion of turbulent atmospheric diffusion by simulation of advec-
tion and diffusion. The accuracy for advection was independent
of the number of particles per cell (as long as at leas^t one
per cell is used for each cell with non-zero concentration) and
was as accurate as the wind field specified at grid cell
centers. Diffusion is simulated using a non-solenoidal veloc-
ity field so particles must be sufficiently dense to simulate
the non-uniform transport. Test problems conducted with the
PICFIC code established that the simulated diffusion was less
than the correct diffusion when there were too few particles
per cell. In a simulation with PICFIC using ten particles per
cell, the effective diffusivity in the simulation was over 90
percent of the correct value. With three particles per cell,
the effective diffusivity in the simulation had dropped to
60 percent of the correct diffusivity.
PICFIC was used to simulate advection and diffusion of
a sub-grid scale distribution (a = Ax/8). Though PICK was
developed to solve the diffusion equation to resolution of a
grid cell, the sub-grid scale distribution was advected with
an effective wind velocity equal to the correct value to three
63
-------
3SR-844
significant places and was diffused with an effective diffusiv-
ities of 92 and 60 percent of the correct values in the parallel
and cross wind directions, respectively.
Comparisons of PICK and finite difference methods of
solution for pure advection showed the PICK method introduced
no artificial dispersion and could obtain an accurate solution
for sharply peaked distributions. The results of the first-
order finite difference technique were a rapid dispersion of the
test distribution, corresponding to a truncation diffusivity,
K = l . The fourth-order techniques tested ad-
L |_ AX J
vected a broad distribution (a = 2,5 Ax), but could not solve
the equation with a sharply peaked distribution (with
a = 1.0 Ax, negative concentrations were calculated).
The computer storage requirements for the PICFIC code
were seven variables for each cell, three variables for each
particle, and 3,000 variables used for random positions. This
resulted in a 20,000 word computer program for test calcula-
tions reported. Analyzing the computational time required for
the tests has shown the following time requirements for the S3
UNI VAC 1108:
• Cell-related calculations: 10-15 basic
computer operations or 3-5 x 10" 5 seconds
per cycle per dimension and
• Particle-related calculations: 40-100
basic computer operations or 1-3 x 10" ^
seconds per particle per cycle per dimen-
sion.
Using the PICFIC code, the PICK method was used to
simulate :
64
-------
3SR-844
Pure advection in constant wind field;
Advection in a shearing wind field;
Turbulent diffusion (with at least three
particles per cell) ;
Anisotropic and inhomogeneous diffusion
(two component diffusivity and inversion
layer); and
Settling of particulate matter in a gravi-
tational field.
65
-------
3SR-844
CTPus Page Intentionally Left Blank,)
66
-------
3SR-844
3. NEXUS SIMULATION OF CARBON MONOXIDE IN LOS ANGELES
After testing the PICK method in two dimensions with
the PICFIC code, the method was applied to solve the
three-dimensional atmospheric diffusion equation, Eq. (1.1).
The resulting code was named NEXUS, an acronym for Numerical
Examination of Urban Smog. NEXUS was developed to solve the
general atmospheric diffusion equation and thus to simulate
air pollution in three dimensions including temporally and
spatially variable sources, winds and turbulent diffusivities.
The NEXUS code was applied to the simulation of carbon
monoxide (CO) in Los Angeles on September 23, 1966. The par-
ticular day and pollutant were chosen to correspond to an air
pollution simulation reported by Lamb. The meteorological
•.x/c
data compiled by Lamb v^as used in this NEXUS application.
This minimized the effort required to make the calculation and
also permitted a comparison of both models. The meteorological
data obtained from Lamb "consisted of maps of hourly (3:00 A.M.
to 8:00 P.M.) isotacjcs and streamlines for surface winds, and
the diffusivities (K = K = 1.2 m2/sec, K = 0.6 m2/sec)
x y z
below the inversion height (1,000 feet) used in his work.
NEXUS requires specification of the winds in three dimensions
throughout the volume of the simulation. The winds aloft were
chosen in a divergence-free manner.
Carbon monoxide was considered inert for the seventeen
hour simulation. That is, the simulation contained no CO
sinks; the only CO depletion occurred by flow out of the basin.
67
-------
3SR-844
Since Lamb's CO source inventory was not available, other pub-
lished data were required to construct the emissions inventory.
The main source of CO in Los Angeles is traffic, contributing
97.7 percent of the 11,000 tons per day as reported by the
17
Los Angeles Air Pollution Control District source inventory.
1 8
Eschenroeder and Martinez published a Los Angeles County
19
traffic distribution for 1951 that they attributed to league
20
but has since been traced to Larson, et.al. Both spatial
and temporal distributions were reported. This 1951 traffic
distribution was used, but was normalized to the 1966 daily
CO emissions.
The initial CO concentrations throughout the volume of
the basin were required to initialize the simulation (at 3:00
A.M.). The only available data were ground level observations
at 12 LOS Angeles County Air Pollution Control District (APCD)
monitoring stations. These data were interpolated on the sur-
face by Lamb, aloft concentrations were arbitrarily assumed to
fall off to a background level of 5 ppm. The methods used to
enter the initial concentrations, source and wind data into
NEXUS will be discussed in Section 3.2.
The Los Angeles application of NEXUS used a grid of
22 x 16 cells in the horizontal, each 4.47 x 4.47 kilometers;
a map of the grid is given in Figure 19. The cells were 100
meters thick and four cells were used in the vertical, re-
sulting in a total height of 400 meters.
The boundary conditions used in NEXUS were an imper-
vious barrier at the ground with transmittive sides and top.
The bottom of the lowest layer of cells was assumed to follow
the terrain and be at ground level. The vertical velocity
used to transport the particles was always assumed to vanish
at ground level and vary linearly within the first layer of
cells. No particles were permitted to enter the ground. If,
68
-------
LOS ANGELES BASIN
(JD
Figure 19 - Calculational Grid for NEXUS CO Simulation.
-------
3SR-844
by a calculational overshoot a particle were below ground
level, it would be re-positioned above. At the transmittive
boundaries, when a particle moved into the outer half of the
boundary cell it was removed from the calculation, simulating
pollutant flow out of the region. Since the boundaries were
chosen arbitrarily, there is also the possibility of pollu-
tant flow through the side boundaries into the simulation
region. To treat this, NEXUS used large Eulerian cells out-
side the border. These border cells were outside the volume
of the simulation (the central grid) and only a gross treat-
ment was used. The border cells were initially (3:00 A.M.J
assumed to have above-background pollutant concentrations,
chosen so that pollutant concentrations were continuous
across the central grid boundary.
This work has been reported previously in preliminary
form.21'22
3.1 DESCRIPTION OF THE NEXUS MODEL
In this section, the structure of the NEXUS code will
be described. A macro flow chart is shown in Figure 20, with
code listings and detailed flow charts included in Volume II,
Section 2.
NEXUS and PICFIC both are based on the PICK method.
In developing NEXUS from PICFIC, the main changes .were the ex-
tension to three dimensions, the input of mean winds, diffu-
sivities, and source emissions frpm tape each cycle, and the
treatment of the borders and the upper and lower boundaries.
The main program of NEXUS directs the flow of the cal-
culations and calls the major subroutines. The INPUT sub-
routine is called first and reads the cards that specify the
problem to be run. The NEXUS code has the option of stopping
70
-------
MAIN PROGRAM
INPUT
3SR-844
NO
>( SETUP
\
OUTPUT
DIFFUS
PARCEL
BORDER
SOURCE
OUTPUT
YES
C STOP J
Figure 20 - Macro Flow Chart of NEXUS Main Program,
71
-------
INPUT
3SR-844
SUBROUTINE
INPUT
READ
SPECIFICATION
FOR RUN
-K RESTRT
RESTRT
READ
RESTARTING
VARIABLES
Figure 20a - Macro Flow Chart of NEXUS Subroutines INPUT and
RESTRT.
72
-------
3SR-844
SETUP
/ READ
INITIAL
I CONDITIONS;
DISTRIBUTE
INITIAL PARTICLES
AT RANDOM WITHIN
CELLS FOR
INITIAL CONDITIONS
OUTPUT
» PARTICLE
POSITIONS
Figure 20b - Macro Flow Chart of NEXUS Subroutines SETUP and
OUTPUT,
73
-------
3SR-844
DIFFUS
PARCEL
CALCULATE
DIFFUSION FLUX
VELOCITIES
r
CONVERT TO
CELL FRACTION/
TIME STEP AND
ADD MEAN
VELOCITY
ZERO
CONCENTRATION
ARRAY
CALCULATE
AREA-WEIGHTED
VELOCITY
PARTICLES
OFF GRID?
ADD PARTICLE
WEIGHT TO
BOUNDARY
FLUX
ADD PARTICLE
WEIGHT TO
CONCENTRATION
ARRAY
FINISH
ALL
PARTICLES
Figure 20c - Macro Flow Chart of NEXUS Subroutines DIFFUS and
PARCEL.
74
-------
3SR-844
BORDER
SOURCE
UPDATE
BORDER CELLS
BX1, BX2,
BY1, BY2
AMD FLUXES
CREATE NEW
PARTICLES
FOR SOURCES
ADD PARTICLE
WEIGHT TO
CONCENTRATION
ARRAY
Figure 20d - Macro plow Chart of NEXUS Subroutines BORDER
and SOURCE.
75
-------
3SR-844
and restarting a given problem from magnetic tape. The choice
of starting mode is specified to INPUT. Restarting is accom-
plished with a call to the RESTRT subroutine from INPUT. On
the other hand, if a new problem is to be generated, a call to
SETUP is made from the main program. SETUP reads the initial
concentrations from tapes. The cell concentrations are con-
verted to pollutant mass and assigned to particles. The
initial particle placement is random within each grid cell to
eliminate any systematic bias. Initial concentrations at
ground level within the central grid are known; the upper
level within cells and the borders are initialized with an
assumed exponential decrease from the known concentrations.
That is, the initial vertical profile of pollutant distri-
bution over background concentration is taken to be reduced
by one-half every 100 meters, c(i,j,k)=(c(i,j,1)-b)/2(kAz/100)+6
for k > 1 . Away from the central grid the concentration de-
creases a factor two for each 16 km of distance. These are
regions with no observational data and thus assumptions are
required for full specification of the initial conditions.
The function of the subroutine OUTPUT is to generate computer
edits to print concentrations.
The main loop of the main program calls subroutines to
calculate in order: the turbulent flux velocity (diffusion),
the transport of particles, the border fluxes, the source
emissions, and, if desired, edits of the results. The DIFFUS
subroutine reads from tape the cell velocities and diffusivi-
ties and the time step. The turbulent flux velocity is com-
puted and added to the mean velocity for each cell. This
total velocity is then converted to units of cell fraction per
time step. Finally, DIFFUS zeroes the cell concentration
array which will be recalculated in PARCEL. Historically,
the particles in the PICK method have been called "parcels"
76
-------
3SR-844
and the transport subroutine name is PARCEL. Since the total
velocity has already been calculated in DIFFUS, PARCEL begins
by calculating the area-weighted velocity associated with
each particle position and then moving the particles. If the
new particle position is outside of the central grid, the
particle contributes to a flux across the border that will be
treated in BORDER. Particles within the central grid have
their mass apportioned by three-dimensional area-weighting
among the eight closest cells and added to the cell concen-
tration array; this procedure is applied to all particles.
The function of the BORDER subroutine is the treat-
ment of pollutant mass outside the central grid. In these
regions the actual concentrations and sources are unknown
and a sophisticated treatment is unwarranted. Treatment of
these regions at all is only warranted if there is pollutant
flow into the central grid. The border treatment used in
BORDER assumes cells twice as high and three times as long
as central grid cells. In the border cells, advection is cal-
culated using first-order finite differences and diffusion is
ignored. Pollutant flow into the central grid contributes to
border fluxes that are treated as additional sources.
The upper boundary is placed such that the inversion
height is between the third and fourth layer of cells. The
diffusion of particles upward through the fourth layer is very
slow a'nd, as will be discussed in the next section, the advec-
tion is also minimal. Particles that do penetrate past the
middle of the fourth layer are removed from the calculation.
At the lower boundary (for particles below the middle of the
first level) , the three-dimensional area-weighting is reduced
to horizontal area-weighting because there are no cells below
the particle. The vertical component of the total velocity
is set to zero at the bottom of the grid and linearly inter-
polated to the particle position as described in the next sec-
tion.
77
-------
3SR-844
The SOURCE subroutine adds new particles to the calcu-
lation in simulation of source emissions and pollutant flow
across the borders. The horizontal coordinates of the new
particles are chosen at random within the source grid cell.
In the Los Angeles basin application, carbon monoxide is
emitted by motor vehicles which are assumed to be located at
the bottom of the first level of cells. In order to simulate
this low level source, the vertical random placement of newly
generated particles is restricted to the lowest one-third of
the bottom layer of cells — approximately the lowest 100 feet,
The new particle mass contributes to cell concentrations
through area-weighting.
This concludes the NEXUS calculational cycle. The
cycle is repeated until the time period for simulation is
completed. The position of each particle as stored in the
computer traces out the trajectory of the pollutant mass.
Each cycle may be plotted on a Stromberg-Carlson 4020 to
form a movie of CO patterns for the day.
3.2 INPUT DATA
The data required for a NEXUS simulation are the source
emissions, mean winds and diffusivities specified for each cell
throughout the grid during the time period being simulated.
The meteorological data were obtained from Lamb. The winds
were in the form of hourly averaged horizontal surface stream-
lines and isotack maps as shown in Figure 21. These were
interpolated and extrapolated to the cell centers of the grid,
as shown by wind velocity vector plots in Figure 22. Eight
directions [north, northeast, east, southeast, etc.) were used
for entering the winds into the computer which resulted in
discretized patterns in the vector plots. Wind maps for each
hour from 3:00 A.M. to 8:00 P.M. were used in the NEXUS simu-
lation. Initial surface pollutant concentrations for the
78
-------
3SR-844
rYs ^-i'' *C"' • A'^- / '^ ViCA^j
X^/ -^ /> -yA'^) >>&?H/]
4->^X\ ^>^^--^^s>"~>v< |\i
—\ >•%_"•• >-^1- ^lA J'\
Figure 21 — Supplied Wind Data (Reference 16). The stream-
lines are shown h)r solid lines. Isotacks in
miles per hour are dashed contours.
79
-------
LOS ANGELES BASIN
9 00
oo
o
r r A A A
i i t i t SAW* AD* HISi
A
i
oo
-fa.
Figure 22 —Winds at 9:00 A.M. in NEXUS Grid. The arrowhead points in the
direction of motion and the arrow length is proportional to
the wind speed.
-------
3SR-844
simulation were also those of Lamb. The 3:00 A.M. observa-
tions of carbon monoxide made by the Los Angeles APCD were
interpolated to generate the map of initial surface concentra-
tions shown in Figure 23. The diffusivities used by Lamb were
1.2 m2/sec and 0.6 m2/sec in the horizontal and vertical direc-
tions, respectively, and were constant in time and space.
These diffusivities were used in the S3 simulation below a
1000 foot inversion height.
The CO source emissions inventory used in the simula-
20
tion was compiled by Larson, et.al. The spatial traffic
distribution is shown in Figure 24 (after interpolation) and
the temporal distribution is shown in Figure 25. The total CO
emissions within Los Angeles County in 1966 were approximately
107 kilograms per day, as reported by the California Air Re-
17
sources Board. These total dail)
normalize the distribution so that
17
sources Board. These total daily emissions were used to
24 22 16
TDE = y^ y^ y^ si.(h) , 0.3,1)
h=l j=l i=l
where TDE is the total daily emissions, S^.(h) the hour
"h" emissions within cell (i,j) for the 22 x 16 cells in the
lowest level of the NEXUS grid,
A computer code, SETUP, was developed to translate the
data mentioned above into the computer tape format for the
NEXUS code. A general flow chart of SETUP is given in Figure
26 and a detailed listing, flow chart and sample problem input
and output is presented in Volume II, Section 2.
The main program of SETUP first calls the INITAL sub-
routine, which reads cards specifying the average CO mass per
81
-------
GO
-
«W)
n> \ xii ^.Vv=/-
Figure 23
K \ I
- Initial C3:00 A.M.) CO Surface Concentrations (ppm) (Re£. 16).
-------
oo
12
IS
8 13
10
13
16
22
22
11
10
10
26
23
20
16
12
40
36
24
10
21
19
19
15
13
12
17
12
13
14
12
10
12
21
10
10
1
2
U-l
i
OO
Figure 24 — Los Angeles County Traffic Distribution,
(Unnormalized relative car miles per cell.)
-------
EVENING RUSH HOURS
OO
8
7
6-
5-
PERCENTAGE OF
LOS ANGELES BASIN 4
TRAFFIC/HOUR
3-
-
-
MORN;
i i
LNG RUSH HOURS
i
"
i.
1 •
t
j
1 i l
• •'.' •'•'.' "!•••.'• '.';'•'• •.';?•.'
V -- ? " 'i, *" .'. .'.'- ; .'
.'."/--'• ; ' ' - ,'i ",''•'•••'
f^-.'V'.' •- .' ^iV-V ••
iVO^v;/^v;:'
••'-•",•'.•' •':•'.'• ;'v
:i'..%'.."./!;''.:v ;•'•:.:'•
•V '*"•• '.'''' '~ \
''•' ~r ' V1"' ''••"«• • '
.tTv'--'''. ';•-''-'•''• ''- "
V'V":'v;-:;'i'i':::.^'
.';•>'. •/-."'--'v" -)-'-'--'i'-
r-'ji-V-.''^-*, •;':{=;
j^v-Vv.J/1-;1-''::':.!-;.;
1 1 1
MIDNIGHT 2
8 10 NOON 14
TIME OF DAY
16 18 20 22 24
50
i
OO
Figure 25 — Temporal Distribution for Los Angeles Basin Traffic.
-------
3SR-844
VELOC
SETV
TIMSTP
VELOC
DIFF
SOURCE
YES
figure 26 -Macro Plow Chart for SETUP,
85
-------
3SR-844
particle, background concentration and the initial CO con-
centrations. The background CO concentration must also be
specified. The use of a fixed background concentration was
suggested by the APCD observations in which measured CO levels
never fell below 5 ppm during the simulation period. The ad-
vantage of a background concentration is a reduction of the
total pollutant mass that is tracked by the particles. With-
out a background concentration, either more particles or more
mass per particle would have been necessary. With a higher
background the observed low concentrations could not have been
simulated. The average CO mass per particle was chosen so
approximately 10 particles per cell would be used in the simu-
lation. Average CO concentrations were 10 ppm so each particle
represented one-half ppm (added to the 5 ppm background).
The main program of SETUP next reads cards specifying
temporal and spatial source emissions distributions. The
VELOC subroutine is called to read cards specifying hourly
winds. Two hourly wind card sets are read and linear inter-
polation is performed to compute the ground level, horizontal
winds at the required time.
The main loop of SETUP is a sequence of calls to sub-
routines SETV, TIMSTP, VELOC, DIFF and SOURCE. SETV uses the
ground level, horizontal winds to calculate the wind field by
an arbitrary construction that results in a divergence-free
field. The horizontal winds on the second level are assumed
equal to those in the first level. The continuity equation,
Eq. (1.4), is solved on the first level for the vertical com-
ponent of the wind at the second level as follows:
86
-------
3SR-844
wi - ! = ° »
- -3Az|Ui + l, j ,:
w.
2Ax 2Ay
i,j,2 = ij,l , vi i 2 = vi i 1 (3.2)
-Ljjj'- J-jJj-1-
The vertical component on the third level is assumed equal to
that on the second. The horizontal winds on the third and
fourth levels are set to be opposite in direction and of the
same speed as those of the first and second levels. The fourth
level vertical wind is set to zero, as shown below:
w. . _ = w. • 0
i,j ,3 i,j ,2
C3.3)
w, , , = 0 ,
By this construction, the interpolated wind field is made a
divergence-free solution of Eq. (1-4) for level 1 in the form;
87
-------
3SR-844
v^ ^ •1 -1 vi,j-1,1
2Ax
2Ay
, Wi, 1,3/2 ~ °
Az
2
where w. . _ ,~ is interpolated as w. . ,/7 = -- w. . ?
i ,3 ,3/2 ^ i ,J ,3/^ 3 i ,j ,z
For other levels the wind field is only approximately diver-
gence-free .
The time step is set in subroutine TIMSTP by requiring
0.4 Ax >_ uAt , 0.4 Ay >^ vAt , 0.4 Az >^ wAt , for all cells,
i.e., during one time step the greatest wind-driven motion is
restricted to be less than four-tenths of a cell dimension.
Next the time is updated, T = T + DT . At the new time VELOC
is called to calculate the new surface velocities.
The subroutine DIFF sets the diffusivities:
K = K = DH = 1,2 m2/sec
x y '
K DV ^'^ mVsec below the inversion height
z 0.0 m2/sec above the inversion height
Finally, SOURCE calculates the source emissions during the time
step for each cell from the normalized emissions distribution.
The calculated source emissions in each cell for the time step
is checked against the average particle mass (particles are re-
quired to be within 20 percent of the average mass). Since
only integer particles can be used in NEXUS, the cell emissions
are added to a storage array until the summed mass is sufficient
to create a particle in the cell.
88
-------
3SR-844
The calculated wind field, diffusivities and source
emissions for each cell are written on tape for use in NEXUS,
The main loop of SETUP is repeated until the time, T ,
reaches the final time.
3.3 RESULTS OF NEXUS SIMULATION
The NEXUS calculation, simulating the Los Angeles air
basin for 17 hours, used 12 minutes of central processor unit
time on the S3 UNIVAC 1108 and used 60K of computer memory.
The numerical results were in two forms: cell concentrations
and particle positions for each of the 163 cycles required.
The cell concentrations were hour-averaged and contour plots
(isopleths) were made for the ground level concentrations.
For example, the 10; 00 A.M. isopleths are shown in Figure 27.
The particle positions can also be plotted to illustrate the
NEXUS resultsj for the ground level at 10:00 A.M. these are
shown in Figure 28. The particle positions are not hour-
averaged, but are actual positions at 10:00 A.M.
The NEXUS results were compared to observed CO concen-
trations made by the APCD and compiled by Lamb, and to the re-
sults of Lamb's model. Lamb averaged the APCD observations at
each of the 12 stations for the 17 hours simulated. These
data, the computed results of Lamb.'s model, and the NEXUS re-
sults are shown in Table V. In these comparisons, the obser-
vations are taken to be "correct," even though the measure-
ments are known to be in error (in 1966 the APCD did not cor-
rect CO observations for water vapor and thus the measurements
can be expected to be 3-5 ppm too high). This error also is
present in the initial concentrations used to start the NEXUS
simulation.
Three statistical tests were used for comparing the
NEXUS simulation with observations:
89
-------
LOS ANGELES BASIN
10 00
to
o
c/4
1
oo
Figure 27 - Isopleths for CO at 10:00 A.M.
-------
LOS ANGELES BASIN
10 00
S^H^^
t: •'•".' -V^BiM*-'-^-.'•;. :.... •..' A'::•.••.'-.• "•' V'..:.. •> .
•.'.-'.' -'V:-:!:-:--V'v;'....'.•••.••:.'..••. •--"-:.AAZOSX" /V-->; 'V.
'£;-^ "• : •' "••*• •••>.••" j^ytv
,^:;;;:i-:./.v:-^;XX>.r:.v::: .. -V: -.:V •., V-. ;.^V-
--- • •- • • '•••/
LA HAM*
.£*
ENTEJ
PUENTE HILLS
/ft:
W^'AiC
A/^S/T^
ANA'MT? •
A"
i
oo
Figure 28 - Plot of Particle Positions at 10:00 A.M.
-------
3SR-844
TABLE V
3:00 A.M. to 8:00 P.M. AVERAGE CARBON MONOXIDE CONCENTRATIONS
(CO in ppm)
STATION
Downtown
Los Angeles
Azusa
Pasadena
Burbank
East Los
Angeles
West Los
Angeles
Long Beach
Hollywood
Pomona
Lennox
Anaheim
La Habra
OBSERVED
16
13
17
16
12
16
14
17
13
13
9
6
SIMULATED BY
THE NEXUS MODEL
18
9
17
17
16
12
12
17
9
10
9
8
LAMB'S MODEL
SIMULATION
(Ref. 16)
22
3
15
7
5
13
8
7
3
11
7
3
i
92
-------
3SR-844
(1) Mean difference:
N
6 = C°i - Si)/N
(2) Root mean square difference (RMSD)
N
CO. - S.)2/N , (3.6)
(3) Correlation coefficient:
N
>i - S)
p = •L"'L = , (3.7)
N
, - S)2/N
where 0. is the observation at the ith station, S. the
]_ j_ ^^
simulated concentration at the ith station, and 0 and S
the mean of 0- and S- , respectively.
The mean difference is a measure of the accuracy of the
model results averaged over the basin. Such basin- and day-
averaged concentrations depend on total emissions, net pollu-
tant flow through the boundaries, and average basin ventilation
both vertical diffusion and mean wind. The mean difference for
the NEXUS results was 0.73 ppm, well within the error intro-
duced in the measurements by variations of water vapor. The
calculations from Lamb's model had a mean difference of 4.8 ppm,
93
-------
3SR-844
This could be partially due to his assumption of CO-free air
off the ocean and perhaps his source emissions inventory.
The root mean square difference (RMSD) is a measure of
the deviations of calculated and observed concentrations at
each of the stations. The NEXUS results had a RMSD of 0.67
ppm, Lamb's model 4.8 ppm. Thus, NEXUS simulated the spatial
distribution of CO (day-averaged) more accurately. Another
statistical test of accuracy of the spatial distribution of CO
is the correlation coefficient. The NEXUS results had a cor-
relation coefficient p = 0.73 and Lamb's model p = 0.35 .
Many air pollution models have had to use additive and multi-
plicative factors to reduce the RMSD of model simulations and
observations. These adjustments of the results are a linear
transformation of the form:
S-' ~ A + B S- ,
where constants A and B are chosen to increase the accuracy
of the model to simulate the observations. However, linear
transformations do not change the correlation coefficient. In
the way it is used here, the correlation coefficient represents
the similarity of the spatial distributions of observations and
simulated results. Thus, NEXUS was better at the simulation of
high and low concentrations in the proper spatial distribution
than Lamb's model.
This may be explained on the basis of the 10:00 A.M.
isopleth plots. The results of Lamb's model for 10:00 A.M.
are shown in Figure 29, and can be contrasted with the NEXUS
10:00 A.M. predictions in Figure 27. Lamb's model predicts
highly peaked concentrations — up to 150 ppm. These are not
observed, but correspond to convergences in the horizontal
winds as shown in Figure 30. Lamb's model used two-dimensional
94
-------
to
en
to
>=>
i
oo
Figure 29 - Isopleths for CO (in ppm) at 10:00 A.M. as
Simulated by Lamb's Model. (Reference 16.)
-------
WIND DATA -
10 00
A «fto. r r r A^ A A, A SAN*
V A A
A r r 4 r A
$.„-»„*
S^NTA MONICA H»LITW^» A * A , A „
JAD^ilELT MT-S -t ^ 1
V
A- v'M
AIU»A .^
A £
A*
- /
PACiriC OCEAN
OO
Figure 30 — Surface Winds at 10;00 A.M. Showing Convergences Near Downtown and
Long Beach,
-------
3SR-844
winds that are not divergence-free. NEXUS uses a three-
dimensional wind field and so can avoid these unphysical ef-
fects .
The day- (or 17 hour) averaged concentrations have
smoothed over many of the significant fluctuations of the CO
concentrations. Hourly-averaged CO concentrations from both
the APCD and NEXUS are plotted in Figure 31. The shapes of
these curves are similar but in Downtown Los Angeles and Holly-
wood, the timing and magnitude of the morning maxima are differ
ent. This could be attributed to two assumptions in the
input data — the source inventory and inversion height. The
source inventory used a basin-wide average time distribution,
whereas Downtown Los Angeles and Hollywood could be expected to
be strongly influenced by early morning freeway traffic. The
constant inversion height was used mainly for comparison with
the work of Lamb, this also would have the effect of dispersing
the early emissions too rapidly. Each of these possible
sources of error is discussed further in Section 3.4.
The NEXUS simulation shows the realistic variation of
pollutant concentration with source emissions and meteorologi-
cal conditions through mean winds and turbulent dispersion.
The following sequence of figures (Figures 32 through 36)
illustrates the patterns and concentrations that develop dur-
ing the 17 hour simulation.
At 3:00 A.M. carbon monoxide is distributed throughout
the basin with higher concentration in the central basin area,
as shown in Figure 32. The 3:00 A.M. winds are represented by
arrows in Figure 33. The winds show a slow flow toward the
coast line. Carbon monoxide concentrations build up in the
slow-moving air throughout the morning rush hours.
97
-------
3SR-844
Hourly CO Concentrations
30-
20-
10-
Downtown Los Angeles
3 4 it, 7 8 9 10 II 12 13 H 15 16 17 IS 19 20
40-,
30-
Lennox
3 4 5 6 7 8 9 10 11 12 13 U 15 16 17 18 19 20
CO
(ppm)
Long Beach
3 4 5 6 7 8 9 10 II 12 13 14 15 16 17 18 19 20
Present Simulation
Measured Concentration:
I I
I
I
I
Hollywood
f=^
3 4 5 6 7
9 10 U 12 13 U 15 16 17 18 19 20
Time of Day
Figure 31 - Simulated Hour-Averaged CO Concentrations Compared
to Measured Concentrations at Downtown Los Angeles,
Lennox, Long Beach and Hollywood.
98
-------
LOS ANGELES BASIN
3 00
I
oo
Figure 32 - Surface CO Concentrations Used for 3:00 A.M. Initial Concentrations.
-------
WIND DATA
3 00
o
o
> . ^ > > -
PACIFIC OCEAN
CO
?3
i
oo
Figure 33 —At 3:00 A.M. the Winds Were Calm and Generally Toward the Coast.
-------
3SR-844
The winds from off the ocean increase during the morn-
ing hours; the 10:00 AM. winds are shown in Figure 30. There
is a flow towards Burbank, but the winds are light there. The
winds are also blowing pollutants eastward off the map towards
Riverside-San Bernardino. The winds form a vortex flow around
the Palos Verdes Hills, accompanied by low winds over Long
Beach. Heavily polluted air flows into the San Fernando Valley
in the Northwest corner of the map either past Burbank through
the Cahuenga Pass or over the Santa Monica Mountains. The air
then flows slowly across the Valley to the North. Air reach-
ing Long Beach via the vortex around the Palos Verdes Hills
has had sufficient time to accumulate high pollutant concentra-
tions . The low winds at Long Beach permit further increase of
concentration. Air coming off the ocean still retains an
above-background level of carbon monoxide due to pollution from
the previous day which was blown out to sea. The air passes
over the coastal traffic on its way towards Downtown Los
Angeles and additional CO is picked up. Thus, by the time the
air reaches Downtown Los Angeles (or even later, Pasadena),
pollutant concentrations are very high, Figure 27.
By 3:00 P.M. the sea breeze is fully established and
is sweeping across the basin, Figure 34. The low winds in the
Northwest corner are the first to be reversed. The Palos
Verdes-Long Beach area still has a low-speed vortex in the wind
pattern. The strong winds are now carrying air originating
further off shore. This air mass is probably less polluted as
it reaches the coast. Moving rapidly across the basin, the
air has little time to accumulate pollutants. Carbon monoxide
concentrations are reduced starting at the southwest corner of
the basin, Figure 35. In the San Fernando Valley, the wind
reversal and low wind speeds have resulted in high CO concen-
trations moving towards Burbank.
101
-------
WIND DATA -
15 00
.
A A
4 4 J
SANTA MONICA
* MTS
t
r * t f /
TBUHIANK
t / T"^ t /• /"
///* / / / /
Aj;>s/-T^ //•-»/• ^ | t / / / /
POMONA7
' /
/ Z-1
\ \ \ \ \ \
\ \ \
I
oo
Figure 34 — High Winds Throughout Most of the Basin with Convergence
Near Burbank at 3:00 P.M.
-------
LOS ANGELES BASIN
15 00
Figure 35 — The Morning High Concentrations Have Been Blown Into the San
Gabriel Mountains and Puente Hills but a New Build-Up is Form-
ing in the Northwest at 3:00 P.M.
oo
-------
3SR-844
The end of the calculational day at 7:50 P.M. is
shown in Figure 36. The high CO concentrations in the north-
ern portion of the basin are caused by both polluted air from
the Oxnard plane (northwest off the grid) and low wind speeds
along the wind front during evening rush hours.
3.4 TESTING OF NEXUS
The NEXUS code was extensively tested to determine
(1) sensitivity of the results to wind field and boundary as-
sumptions, to input data and to the initial random locations
of particles; (2) possible improvements through a new source
inventory and a non-constant inversion height; and (3) the ac-
curacy of variants of the PICK method and comparison of NEXUS
computing time with that required by finite difference methods.
The NEXUS calculation as described above will be referred to
as the standard calculation.
The first series of tests was directed toward evaluat-
ing the effects of the wind field construction, Eq. (.3.2) and
(3.3). The necessity of a divergence-free wind field was seen
in the NEXUS comparison to Lamb's model and will be shown
again in some of the variants of the PICK method to be de-
scribed below. The lowest two layers' winds were constructed
^, o 4>-
from horizontal ground level observations and the continuity
equation (requiring the winds to be divergence-free). The
winds in the upper two levels were chosen arbitrarily but ap-
proximately divergence-free. To test the effects of this choice,
we conducted two other test calculations. The first test had
all velocities on levels 3 and 4 set to zero. This resulted in
day averaged concentrations above 20 ppm in Pasadena, Burbank,
and East Los Angeles and higher concentrations than the
standard results everywhere. This can be explained in terms
of the vertical motion of the particles — they could be
104
-------
LOS ANGELES BASIN
19 50
o
l/i
50
i
CO
Figure 36 — Evening Traffic has Produced a Local Build-Up of CO Over Long Beach
and a Large Region of High Concentrations Corresponding to the Sur-
face Wind Convergence in the Burbank-Pasadena-Hollywood Area at
8:00 P.M.
-------
3SR-844
advected by convergent winds and build-up without the possi-
bility of escaping into the vertical as required by the con-
tinuity equation. Above the middle of the second layer, the
continuity equation would no longer describe the interpolated
wind field. The second test retained the vertical component
of upper level winds, but used zero for the horizontal com-
ponents on levels 3 and 4. In this test the results differed
by no more than 1 ppm from the standard results for day
averaging. This indicates the insensitivity of the results
with NEXUS to the form of the wind field above the layers
that affect the ground level (lowest two layers due to the
area weighting) as long as vertical transport is retained.
The testing of the boundary conditions consisted of
a calculation with the concentrations in the border cells set
to zero (only air at the background CO concentration, 5 ppm,
entered the central grid). The effects of this were minor
except during the late afternoon convergence. If the assumed
background concentration were lowered the borders could have
a greater effect. A series of tests was also conducted to
determine the effects of the choice of the background concen-
tration. In successive calculations the background concentra-
tion was lowered until, with a 5 ppm background, there was no
appreciable change in day averaged results. The number of
particles (and, thus, the average mass per particle) were
varied in another series of tests. An optimum for computer
time and accuracy was determined for 4,000 - 8,000 particles
(the number of particles varies as they are created by source
emissions and flow off the central grid) corresponding to ap-
proximately a metric ton of CO per particle.
Next, the sensitivity of NEXUS results to stochastic
variations in the required input data (winds and source inven-
tory) was evaluated. If no systematic errors are assumed in
106
-------
3SR-844
the data, the mean error in the data is zero and the errors
in the data could be expected to be normally distributed
around zero. That is, if the data used for the standard re-
sults were correct, then the data with stochastic errors
would be scrambled such that
Scrambled = Standard d + R) , (3.8)
where R is a normally distributed random variable with mean
zero. For these tests the standard deviation for R was
taken to be one-fourth. Thus, 90 percent of the time the
scrambled data falls within 50 percent of the standard value.
The first test used scrambled source emissions. Equa-
tion (3.8) was used with the source emissions spatial distri
bution for each cell. A second test used a scrambled source
emissions temporal distribution. The last test of this type
used scrambled winds — each component of the horizontal
ground level winds was scrambled separately by use of Eq.
(3.8). The results of these tests show that stochastic errors
tend to average out. In day-averaged concentrations, the
scrambled data resulted in deviations of less than 10 percent
from the standard results (the station-by-station comparison
is given in Section 3.5); the mean differences, root mean
square differences and correlation coefficients are also in
close"agreement. Hourly-averaged concentrations averaged
within 20 percent of the standard results.
The last of the sensitivity analyses was a comparison
of three calculations each using different random positions
for the placement of particles. Point particle procedure was
used in these tests. Eight hour averages were computed — these
were within 10 percent of each other. For hourly averages, the
maximum variation was under 35 percent and average deviation
107
-------
3SR-844
was approximately 10 percent. The use of the area-weighting
procedure would be expected to reduce this deviation because
the procedure smoothes out effects of particle position.
Two possible improvements in the standard input data
were tested: new source emissions inventory and a time-
dependent rising inversion. A source emission inventory was
23
recently completed by Roberts, et.al. and forwarded to us
by Roth. In Roberts' inventory, traffic is divided between
freeway and non-freeway, each having its own spatial and
temporal distribution as shown in Figures 37, 38 and 39. A
major difference between this new inventory and the one used
previously is the total carbon monoxide emitted per day. In
the previous inventory, the emissions were normalized to the
107 Kg/day reported by the California Air Resources Board.
Roberts' inventory has a total of only 0.6 x 107 Kg/day for
approximately the same region. Other differences are in the
spatial distributions of the two inventories. These differ-
ences may be attributed to traffic changes in the years be-
tween their compilation and the greater resolution of Roberts'
work. The temporal distributions also differ by their claimed
resolution. It was anticipated that the explicit treatment of
the freeway traffic in Roberts' inventory would lead to better
correlation with the observed early morning rise in CO at the
Downtown and Hollywood stations. This was partially borne out
by the NEXUS calculation; the morning rise in CO concentration
at Downtown between 6:00 and 8:00 A.M. was simulated (approxi-
mately an hour later than observed). However, the day averaged
results were poor — having a mean difference of almost 2 ppm,
a RMSD of 3.2 ppm and a correlation coefficient of only 0.55.
Since the correlation coefficient is invariant under a linear
transformation of the results, these poor results with Roberts'
inventory cannot be attributed entirely to the low value of
total CO daily emissions.
108
-------
,09 -
o
ID
freeways
surface streets
•0
0
o - AC
trf i— 1 * "~
4J 3
O
-H O
•J-l i-l
O Cn
-H -iH
4J 03
U w n"5
J^ JJj . U o _
-------
3SR-844
2S
23
21
2.1
20
11
IS
17
/.*
J
ft
13
IZ
i;
fi
S
i
$
3
2
1
I 13 <• £" 6 7 8 1 iti ll
59 !
1
1
1G9 ! 271
._)_
•
. . . j_- ..
i
— -i--- •
L_
t
.-..\ —
r !
i
94
155
242
J222
1
311 303 090
1
|J 352
J356
J242
1 i
""" k !
N. 101
jX
; 1
i 1
| I
I 1
1
' i
|
-p-j —
1
j
i
; 1
i I
1
i
119 !
.J!.
144
156
126
456
137
.
473
426
\. _,
141
295 511
295
' .,
N r
—
- ,
1 2. 3
Figure 38 - Spatial Distribution of Freeway Traffic
Ctkousand vehicle miles per day). (Ref, 23)
110
-------
3SR-844
25
«
23
n
11
2*
11
i?
17
it
if
*
13
12
1)
1C
•i
S
7
&
"^
3
2
1
' ?_.. ?__«:.
35 ! 12 34|149
_. L. ..!._ j. ... .
60 ' 155 1 202! 223
i |
137 ' 255 233; 238
1
239 . 250 i 249' 250
139 1 118 ! 107 194
25 . 30 17 12
9 ' 25 ! 0: 2
1 '
0 1 18 j 0 j 7
47 64 1 99:117
r "X
- ' 4-
i i ;
' r :
j
; i :
i ! •
i
!
i j '
i i
: i
i_ |_
.j ; " ~"
i :
i i ; '
T-*-, TT"1
£
146
252
266
264
182
59
30
98
231
476
\
"T
67 ? S
172 • 135 1 43 1 2S
. - . ., . t
22) \ 152 | 8S! 66
295 1 251 • 233 ( 54
414 328 ! 276 212
350 268 280 247
i
63 121 1 155 , 116
70 ' 32 | 260 ; 478
283 , 403 ; 549 512
512 ' 346 I 459 478
254 , 338 ' 261 412
305 255 255 311
\06 . 266 . 413 392
\
2i ! 261 : 296 397
\189 255 328
•V" " , "
•124 215 412
1 9\ 360 ' 324
; 1 J239 : 304
\^i( 121 ' 188
j\30 j 49 j 28
[V~y^»
i , i
-i-j-4-
it.
17
84
J
20
1 —
215
,180
25
1
375
572
468
413
478
331
285
358
322
231
249
266
142
176 "
-—
10
//
1
73
27
4
204
245
259
473
642
512
477
393
314
176
210
12.1
175
254
1^
J_
_._
_._.
-»-- — V
II
12 >3 /£ if
~l ------ -T-
01 0 : 4 ' 2
7 , 2 0.0
157 50 4 0
131 58 91 113
204 21 173 241
242 191 162 318
164 ' 116 190 236
225 217 ' 229 285
579 290 199 186
323 324 293 264
364 292 273 299
273 ' 322 305 265
254 292 255 386
^ ' i
265 ; 293 257 !298
121 270 . 312 J301
79 204 ' 283 430
J.
94 283 192 !282
i
157 413 | 386 317
y^r^^
! i
i i
, 1 !
--f— -j— f—
-r-H~
! -L. J
ii±J
12. '3 '% >f
It, /7 i$
L
200
0 1 0
21 1 2
212 156 91
308 241 218
227 ' 236 185
243 . 250 242
75 ' 108 ' 173
1
195 127 43
190 185 G9
253 ' 190 295
263 229 234
241 ' 230 135
93 96 104
204 133 - 186
165 .155 173
"1 .
162 . 75 . 116
175 ' 60 125
"^Jal 48 ; 114
_.;Xj-
— j_.p.
i i
" rt~~
TIT
it" i? "»8~
'
Figure 39 — Spatial Distribution of Non-Freeway Traffic
(thousand vehicle miles per day). CRef . 23)
111
-------
3SR-844
The observed early morning rise in CO could also be
explained by reduced dispersion during the morning rush hours.
The standard results were calculated with a fixed inversion
height for comparison with Lamb's work. Next a rising inver-
sion height was included with the standard input for a NEXUS
simulation. The same diffusivities as previously used were
retained (K = K = 1.2 m2/sec and K = 0.6 m2/sec or 0.0),
x y 2
but until 6:00 A.M. K = 0.0 was used for all levels. From
LJ
6;00 A.M. to 12 noon, K = 0.6 m2/sec was used below the in-
Z
version height and K =0.0 above. The inversion height
Zj
was obtained by a linear interpolation between ground level
at 6:00 A.M. and the standard inversion height (300 meters)
at 12 noon. The hourly-averaged carbon monoxide concentra-
tions at two stations are shown in Figure 40.
The rising inversion did have the anticipated effect
of limiting early morning dispersion and the calculated re-
sults for the early morning hours are in better agreement with
the measurements. The day averaged concentrations were higher
in this simulation than in the standard results of NEXUS. In
comparison with the APCD observations, the rising inversion
NEXUS results had a mean difference of -1.17 ppm, RMSD of
3.7 ppm but correlation coefficient of 0.7. Perhaps the
below-inversion diffusivities were too low leading to the after
noon concentrations higher than observed. This is suggested
10
by the work of Eschenroeder and Martinez in which a velocity
dependent vertical diffusivity is used. Their vertical dif-
fusivity would average about 6 m2/sec in the afternoon, com-
pared to the 0.6 m2/sec used here. These two examples illu-
strate some of the improvements possible within the framework
of NEXUS. Next a few modifications to the framework and to
the PICK method used will be described.
112
-------
CO
(ppm)
Downtown Los Angeles
30-
20-
10-
3 A 56 7 8 9 10 11 12 13 U 15 16 \7 18 19 20
30-
20-
10-
.Lennox
f—I—J.
Rising Inversic"
NEXUS Standard-
Measured Concentration
3 4 5 6 7 8 9 10 11 12 13 U 15 16 17 18 19 20
Time of Day
Figure 40 — Hour-Averaged CO (^in. ppm) calculated with an inversion rising height
is compared to the standard results with a fixed inversion height
and to the measured concentrations.
en
oo
-pi
-------
3SR-844
The NEXUS model as reported here has a very general
formulation, describing three-dimensional advection and diffu-
sion. However, by setting input variables (winds and diffusivi-
ties) it can simulate models having only two spatial dimensions.
It also has both Eulerian and Lagrangian features. Modifications
of NEXUS can be made to emphasize either, and thus evaluate the
advantages and disadvantages of Eulerian or Lagrangian methods.
The first modification was the assumption of a single
layer below a constant inversion height and use of only the
observed, horizontal, ground-level winds. This modification
assumes complete mixing below the inversion height. In this
25
respect it is similar to the work of Weisburd, et.al., with
a single cell moving along a Lagrangian trajectory. The
statistical comparisons of the day-averaged CO concentrations
calculated with the modified NEXUS and the observed concen-
trations showed a mean difference of -0.33 ppm, a RMSD of
3.39 ppm and a correlation coefficient of 0.49. This indicates
time- and space-averaged the concentrations were correctly
simulated, but the spatial dependence of the day-averaged con-
centrations were poorly correlated. The hourly averaged re-
sults had sharp peaks of very high concentrations. For
example, at Lennox 42 ppm was calculated for 8;00 A.M. when
18 ppm was observed. These unphysical maxima correspond to
convergences in the observed horizontal wind field.
Two approaches were suggested to alleviate this problem:
(1) "smooth" the observed winds to make a divergence-free,
horizontal wind field and (2) modify the calculated concen-
trations to take into account any non-zero divergence,
v-u $ 0 . The first approach was used by Browne.26 The second
approach was adopted for the next test. At each cycle the
divergence of the wind field is calculated and the concentra-
tions (and thus the particle masses) are modified by an arti-
ficial factor so that;
114
-------
c = c
Ax Ay(l + Atv"«u)
where the term in the brackets is a volume change due to the .
divergence V«u . Since this must be computed for each cell,
the divergence was computed by first-order finite differencing
With this modification, the correlation coefficient for the
resulting day-averaged concentrations was 0.74, but the mean
difference was 3.25 ppm and the RMSD 4.13 ppm. That is, the
results were lower than the observations but the spatial dis-
tribution concentrations were well correlated. This suggests
that the vertical mixing was too high in this test, but the
modification of the concentration to take into account the
non-divergence-free wind field was a step in the right direc-
tion.
The last series of tests evaluated the computer time re-
quired by finite difference techniques with the 3-D Eulerian
grid of NEXUS. The techniques used were similar to the first-
order and Fromm's fourth-order methods described in Section 2.2,
except that splitting was not used. The first-order method
introduced errors that are exhibited as truncation dif fusivities
which, for the grid used in NEXUS and a typical five minute time
step (assuming uAt/Ax -0.4 , vAt/Ay -0.4 and wAt/Az -0.4)
were: K = K = 800 m2/sec and KZ = 4 m2/sec. The fourth-
order method computed negative concentrations that had to be
ignored. The actual simulation results would be expected to
improve if splitting were used, but the time required for
computations would be almost identical. The computer time
115
-------
3SR-844
required for ,17 hour simulation with the standard input data
using the first-order method was four minutes, while the
fourth-order method required five minutes. These can be com-
pared with the twelve minutes required by the standard NEXUS.
3.5 SUMMARY AND CONCLUSIONS
The NEXUS code was developed to apply the PICK method
to problems of regional air pollution. NEXUS simulates ad-
vection and diffusion of a non-reactive pollutant in three-
dimensions and was applied to the simulation of CO concentra-
tions in Los Angeles. In the development of NEXUS, the PICFIC
code was extended to three dimensions, border cells and bound-
ary conditions at the ground were added. All input data was
written on tape by a SETUP routine and read by NEXUS. Using
meteorological and source data available in the literature,
NEXUS results for 17 hour-averaged CO concentrations at 12
stations averaged within 20 percent of Los Angeles Air Pollu-
tion Control District observations. The calculation of CO
throughout the region for 17 hours took 12 minutes on a
UNIVAC 1108. The NEXUS results were in better agreement with
observations than Lamb's model using the same meteorological
data.
The results of series of tests with NEXUS are sum-
marized in Table VI. The sensitivity of NEXUS to the input
data was tested by introducing random variations Cscrambling) in
the data. A new vehicle source inventory and a time-dependent
(risingj inversion height were tested by NEXUS calculation.
NEXUS was used to simulate a two-dimensional well-mixed-layer
model and to test a divergence-free modification. Two finite
difference methods using the NEXUS Eulerian grid were compared
to NEXUS for computing time. The finite difference methods
took one-third to one-half as much computer time as NEXUS for
the same 17 hour simulation, but encountered difficulties men-
tioned in Section 2.2.
116
-------
TABLE VI
DOWNTOWN LA
AZUSA
PASADENA
BURBANK
EAST LA
WEST LA
LONG BEACH
HOLLYWOOD
POMONA
ANAHEIM
LENNOX
LA HABRA
APCD
OBSERV.
16
13
17
16
12
16
14
17
13
9
13
6
Concentration
Coefficient
6
s*
NEXUS
STD.
18
9
17
17
16
12
' 12
17
9
9
10
8
0.73
0.67
2.68
LAMB'S
MODEL
22
3
15
7
5
13
8
7
3
7
11
3
0.55
4.8
6.66
NEXUS
NEW
SOURCES
13
10
'14
14
14
11
11
13
10
12
9
8
0..5S
1.92
3.20
NEXUS
RISING
INVER,
19
10
23
20
20
13
14
18
10
10
11
8
0.70
-1.17
3.67
SCRAMBLED SOURCES
TIME SPATIAL
DISTRIB. DISTRIB.
19
10
19
18
18
12
13
18
10
10
11
8
0.72
-0.33
2.86
18
10
19
18
18
13
13
17
10
10
11
8
0.73
-0.25
2.66
SCRAMBLED NEXUS
WINDS 2D
WINDS
18
10
18
18
18
13
12
16
10
10
11
8
0.71
0
2.68
14
10
22
17
17
12
14
13
10
14
15
8
0.49
-0.33
3.39
NEXUS 2D WINDS
5 V-u=0
MODIFICATION
12
10
11
12
12
9
9
12
10
9
9
8
0.74
3.25
4.13 ^
i
oo
-------
3SR-844
(This Page Intentionally Left Blank.)
118
-------
3SR-844
4. NEXUS/P SIMULATION OF PHOTOCHEMICAL SMOG
In Section 1.3 a technique to introduce the effects
of chemical reactions into the PICK method was briefly de-
scribed. With this technique the NEXUS code was extended to
simulate chemically reacting pollutants, forming the NEXUS/P
model (Numerical Examination of U_rban Smog with Photochemistry)
The major addition was the encoding of the technique by which
particles can change mass during the calculation. In NEXUS,
each particle has a fixed mass, set when the particle is intro-
duced into the calculation. Chemical reactions occurring in
the atmosphere have the effect of removing primary pollutants
and forming secondary ones. The net effects are changes in
the concentrations of the reacting and product pollutants.
The PICK method is based on using particles to represent pol-
lutant mass and so for reacting and product species the mass
represented by the particles must be changed and new particles
must be created.
The turbulent atmospheric diffusion equation being
solved with the NEXUS/P code can be written as
at
_3cs _.
3c
3c
3x
K,
3c
K
9y y 9y
3z
K
z 3 z
* s
r,t
Bsrt cr ct
r,t
(4.1)
119
-------
3SR-844
where the concentration subscripts have been used to emphasize
different chemical species and the two new terms represent bi-
molecular and photochemical reactions (B and P, respectively).
There is a separate differential equation for each species be-
ing simulated. The differential equations are no longer
linear and are coupled together through the chemical terms.
In the PICK method as used in NEXUS/P, the equations
are rewritten as
9c 3uc 3vc 9wc
1 = _ s^ s_ s s f4 21
9t 3x 9y 3z *s > ^ J
where the fictitious total velocity is used for pollutant
transport and the chemical terms have been grouped into a
source or sink for the "sth_" species, S . Eq. (4.2) is solved
assuming the effects of pollutant transport and of chemical
reactions are separable. Pollutant transport is simulated for
a time interval, then chemical reactions are simulated for the
same interval. That is, the transport and reactions are cal-
culated in sequence. The NEXUS code was used for the pollu-
tant transport and new subroutines were added for reactions.
The NEXUS/P code is described in Section 4.1 and photochemical
simulation is detailed in Section 4.2. The results of tests
made with the NEXUS/P code are given in Section 4.3. A photo-
chemical simulation of Los Angeles smog is described in Sec-
tions 4.4 and 4.5.
4.1 DESCRIPTION OF NEXUS/P
The NEXUS/P code was developed from the NEXUS code.
Simulating chemical reactions requires re-weighting the
particles based on cell concentration changes. A simple re-
weighting algorithm is possible with, the point particle
120
-------
3SR-844
procedure, Eq. (1.15). Thus, NEXUS was re-written for point
particles in developing NEXUS/P. In addition, new subroutines
were included to simulate the chemical reactions to calculate
the changes in cell concentrations and to introduce new
particles for product species.
The new subroutines (KEM, CHEM and NPART) are shown
in the general NEXUS/P flow chart, Figure 41. The KEM sub-
routine checks each cell to determine if chemical simulation
is necessary. CHEM is used to simulate the effects of chemi-
cal reactions as will be discussed in Section 4.2. Subroutine
NPART adds new particles. The re-weighting of particle mass
is completed in an addition to PARCEL.
The S3 UNIVAC 1108 has available 64K of central pro-
cessor computer core. This cannot accommodate the require-
ments of NEXUS with multiple species. Basically, NEXUS/P
requires four variables (x, y, z and mass) for each particle;
separate particles for each specie; and winds, diffusivities,
species concentrations and the chemical changes in concentra-
tion for each cell. For 2,000 cells, 10,000 particles and 5
species, the core required would be approximately 230K.
NEXUS/P instead uses bulk storage (magnetic drums) and shifts
data into the central processing unit for computations as
required. In developing NEXUS/P, NEXUS was re-written to
minimize the number of variables needed during each calcula-
tion. This has resulted in a requirement for the cell con-
centrations and chemical changes in concentrations for all
species, or, at another time, the winds, diffusivities, con-
centrations, chemical changes in concentrations, and particle
position and mass for a single species (a maximum of 54K for
the above example) in the central processing unit at a single
time. Two new subroutines, ROLLIN and ROLOUT, were developed
to handle the peripheral storage and manipulation for concen-
tration, chemical changes in concentration, and particle posi-
tion and mass.
121
-------
3SR-844
SETUP (separate computer
initial
conditions
1
time step
program)
winds and
diffusion
_ _ _l
sources
1
J
NEXUS
SETUP
initializes
program and
sets initial
conditions
L_
KEM
calculates
chemical
reactions
1. DIFFUS
read in
winds and
diffusion
2. NPART
calculates
new
parcels
3. PARCEL
moves
parcels
4. BORDER
calculates
mass flux
in borders
5. SOURCE
reads sources
calculate new
parcels
ROLLOUT
writes parcel
locations and
concentrations
Figure 41 -
Macro plow Diagram for Computing Photochemical
Pollution.
122
-------
3SR-844
An operation manual for NEXUS/P was written and is
included as Volume III. This manual includes descriptions
and listings of the subroutines of NEXUS/P and the input
data processing routine SETUP, all translated into IBM 360
FORTRAN IV.
4,2 PHOTOCHEMICAL SIMULATION
The simulation of reactions occurring in photochemical
smog begins with the choice of a set of chemical reactions or
chemical mechanisms. In the atmosphere the actual number of
elemental reactions can be very large, beyond the capacity of
present computers. For example, of the many hydrocarbons in
an urban atmosphere, at most three types of hydrocarbons are
represented in chemical mechanisms reported in the litera-
27 28 29
ture. ' ' Present chemical mechanisms are based on lump-
ing many similar species into a single "effective" species and
many similar reactions in a single "effective" reaction. The
effective hydrocarbon species used in present chemical mecha-
nisms may not be a real hydrocarbon, but instead react as an
average of real hydrocarbons. The reactions used may also
not be elementary reactions and may even be non-stochiometric.
This approach is referred to as the "Lumped Parameter" ap-
proximation .
The NEXUS/P model can use any method that predicts
the ch-anges in cell concentrations due to chemical reactions.
For the Los Angeles photochemical smog simulation reported in
Section 4.5, a set of reactions based on the lumped parameter
29
approximation and reported by Eschenroeder and Martinez was
used. Eschenroeder's chemical mechanism is given in Table VII
The set of chemical reactions can be translated into
a set of coupled, non-linear, first-order differential equa-
tions — the chemical kinetics rate equations. For example,
123
-------
3SR-844
TABLE VII
RATE COEFFICIENTS FOR ESCHENROEDER'S MODEL OF
THE HYDROCARBON/NITRIC OXIDE MECHANISM
(Reference 29)
(Stoichiometry imbalances may occur because of
lumped parameter assumptions.)
1.
2.
3.
4.
5.
6.
7.
8.
9-
10.
11.
12.
REACTION
(hv) + N02 -»• NO + 0
0 (+ 07.) (+ M) -> 0,, + M
LJ O
03 + NO -> N02 (+ 02)
0 + HC -* 2R02
OH + HC ^ 2R02
R02 + NO + N02 + 0.5 OH
R02 + N02 + PAN
OH + NO ^ HN02*
OH + N0? + HNO *
L, O
03 + HC -> R02
CH-0 +) NO + NO- -*• 2HNO **
z z z
(hv) + HN02 ^ NO + OH
RATE COEFFICIENTS
0.4 min"1
2.64 x 106 min"1
40 ppm" min
6100 ppm" min
80 ppm min
1500 ppm" min"
6 ppm min
10 ppm min
30 ppm" min
0.0125 ppm"1 min"1
0.01 ppm min
0.001 min"1
rc
Rate constant incorporates third body concentration.
ft
Water vapor incorporated into rate coefficient.
is used as the chemical symbol for a free radical
RO,
lumped species.
HC stands for the hydrocarbon lumped species.
Parentheses are used to indicate species concentrations
that are not chemically changing and are incorporated into
the rate coefficient.
The solar flux is represented by "hv," the maximum solar
flux is incorporated into the first rate coefficient.
124
-------
3SR-844
the inorganic loop of Eschenroeder's mechanism is:
kl
1. (hv) + N02 »> NO + 0
k2
2. 0(+02)(+M) - » 03 + M
3. 03 + NO
where hv represents ultraviolet photons of the solar flux
(which vary from sunrise to sunset), the third body, M , is
the total concentration of all species present (necessary for
conservation of energy and momentum in the reaction) , and the
k's are the rate coefficients. Species concentrations that
are not varying due to the chemical reaction are usually in-
corporated into the rate coefficient. In the second reaction,
the concentration of molecular oxygen and of M can be as-
sumed constant and incorporated into the rate coefficient
ki = k2[02][M] . The rate for each reaction is given by
TI = k1[N02]
= k[0][NO] , (4,3)
where the square brackets indicate the numerical value of the
concentration .
Thus, the chemical kinetics rate equations (differen-
tial equations for the total time rate of change of a species)
can be written as
125
-------
3SR-844
d[N09J
- *J[0]
kl[N02]
- k2[03J[NO]
= k2[0] - k3[03J[NO] .
(4.4)
These differential equations are non-linear and are
difficult to solve. The rate equations are all of the form
= F - Re
(4.5)
where F is a summation of the rate of all reactions forming
species c , and Re is a summation of all those reactions remov-
ing c . F and R are called the formation and removal rates,
respectively. For example with the full chemical mechanism,
Table VII, the time rate of change of
this form as :
RO
can be written in
d[R02]
dt
j
2k4[0][HC] + 2k5[OH][HC] + k1()[03][HC]j
kfi[NO] + k [N07]j[RO?] .
U i Li } L,
The non-linearities have been suppressed in this form, but
another difficulty, that of "stiffness," can be discussed.
For simplicity in this discussion of stiffness, R and p
are assumed constant. Then the solution of
126
-------
3SR-844
Eq. (4.5) is the sum of a constant term and an exponential
-Rt
term, e . For solutions decaying to zero (no constant
term in the solution, or physically, a small equilibrium
value), the solution is reduced by a factor of 1/e in a time
1/R . In typical applications the values of R for differ-
ent species concentrations may differ by many orders of
magnitude. In this case the times associated with signifi-
cant variations in different species concentrations differ
by orders of magnitude and the differential equations are
called "stiff." The standard methods for solution of first-
order coupled differential equations (Euler's method, Runga-
Kutta, Adams predictor-corrector) have a time step limitation
of the order of 1/R for the largest removal rate. In essence,
the calculation is limited by the fastest reaction even though
it may not lead to significant concentration variations. For
example, in Eschenroeder's mechanism the atomic oxygen is
rapidly removed by reactions with hydrocarbons (reaction 4)
and molecular oxygen (reaction 2), but usually is not varying
rapidly because it is in equilibrium with photochemical pro-
duction from NO- .
An approach designed to treat general stiff equations
was reported recently by Gear. S3 has developed a general
chemical reaction model, KEM-GEAR, using the Gear method for
integrating stiff equations. A test case was computed using
29
KEM-GEAR and compared to results published by Eschenroeder.
Eschenroeder used a Fade" approximant for the exponential solu-
tion of Eq. (4.5) and the two results agreed to the accuracy
of the reported concentrations. An advantage of the KEM-GEAR
approach was the computer time required. On the S3 UNIVAC 1101
the test case took 7.9 seconds; Esch_enroeder reported 50
seconds on a CDC 6400. The CDC 6400 is approximately 20 per-
cent slower than the UNIVAC 1108.
127
-------
3SR-844
Thus, KEM-GEAR was about a factor of five faster than the Fade*
approximant method. The disadvantages of KEM-GEAR are the
program size and core requirements since it is a multi-step
integration method. This motivated the development of a sim-
pler technique.
Some species are formed as rapidly as they are removed
so that their concentration does not vary appreciably. These
species are in a pseudo-equilibrium. This is not a true equi-
librium since the reactants in the forming and removing reac-
tions may still be changing. If a species were in this
pseudo-equilibrium,, the term dc/dt in Eq, (4.3) would be
small compared with F or cr giving
c « F/R , (4.6)
This is much simpler to compute than to integrate Eq. (4.3),
and so the results of KEM-GEAR calculations were analyzed to
determine if, when and for what species
0.999 <
CF/RJ
1.001 , (4.7)
The analysis showed that species 0 , R02 and OH satisfy
Eq. (4.7) and could be accurately solved using Eq. (4.6).
These are free radials and could have been expected to be in
pseudo-equilibrium due to their fast removal reactions, though
HN02, another free radical was found not to satisfy Eq. (4.7).
The terminating products (PAN, HNO,) are only products
of the reaction and so would not control the time step for the
standard methods. Of the other species, the three radicals
mentioned above can be solved from Eq. (4.6) and only the re-
maining five species (NO, N02, 03 HC and HN02) must be
128
-------
3SR-844
calculated by solving the differential equations. A small
code, CHEM, was developed in which these five remaining differ-
ential equations are solved. The solution is based on the
analytic solution of Eq. (4.2) for a time step, At , assuming
constant F and R :
c(t+At) = F/R + [c(t) - F/R]e~RAt . C4.8)
For a sufficiently small At , R and F are almost constant
and the solution is a good approximation. This method has
been used by Keneshea. In CHEM the species concentrations
are predicted as indicated by Eq. (4.8). Then F and R are
recalculated using the predicted concentrations, c , and new
"corrected" concentrations, c , are calculated again using
Eq. (4.8). A desired fractional error limit, ER , is used to
compute the appropriate time step, At , from a trial time
step At1 :
c
At = ER A c
CP ~ Cc
At' . (4.9)
An initial time step of }0 9 seconds is used to insure starting
accuracy. With an ER = 0.005, the CHEM results agreed with
those of KEM-GEAR to better than 10 percent at the N02 maxi-
mum concentration and, at worst, less than a 20 percent dif-
ference at the end of the calculation (150 minutes). CHEM
required only 1.8 seconds computer time for the same calcula-
tion KEM-GEAR used 7.9 seconds . . . better than a factor of
four savings in computer time or about twenty times faster
than Eschenroeder's solution.
The CHEM code was incorporated as a subroutine into
the NEXUS/P code for chemical simulation to form the quantity
129
-------
3SR-844
S in Eq. (4.2). An intermediate subroutine KEM was developed
to call CHEM for each cell in which chemical reactions are
significant. The CHEM subroutine calculates new cell concen-
trations due to chemical reactions for a time interval equal
to the NEXUS/P time step. The CHEM subroutine returns the
fractional change in concentration for each species:
_ cCt+At) C4 1(n
CF -- - ' L J
If the initial cell concentration of a species produced by
chemical reactions were zero, CHEM would return the negative
of the concentration produced:
CF = -c(t + At) .
The negative sign is used as a flag to indicate that CF is
the new concentration instead of the fractional change.
4.3 NEXUS/P TESTING
The accuracy of the NEXUS/P code was tested for (1) a
stationary pollutant puff with chemical reactions, (2) an ad-
vecting puff with chemical reactions, and (3) an advecting
and diffusing puff with chemical reactions .
The first test evaluated chemical reaction simulations
in the NEXUS/P formalism. The pollutants (NO, N02 , and HC)
were initially set in a cube of five cells on a side. The
winds and dif fusivities were set to zero, In NEXUS/P initial
concentrations were converted to particle masses. The CHEM
subroutine was used to calculate the fractional changes in
concentrations for each species and cell. Then the particle
masses were re-weighted to account for changes in mass due to
130
-------
3SR-844
chemical reactions. Finally, the new concentrations were
calculated from particle positions (which had not changed
because the winds and diffusivities were zero) and masses.
Since no advection or diffusion was permitted, the changes
in the concentrations depended only on the reactions and so
could be compared directly to the results of the CHEM code
above .
In the second test, the simulation of advection was
added to that of chemical reactions. In this test, the non-
linearities of the chemical kinetics rate equations are im-
portant. As the initial 125-cell pollutant puff is advected
downwind, the cell concentrations along the leading and trail
ing edges do not remain constant, and the non-linear chemical
equations introduce errors.
To illustrate this important point, consider a one-
dimensional problem as illustrated in Figure 42a.
Cell 1
0.0
1.0
1.0
1.0
0.0
0.0
Figure 42a — Initial Distribution of Primary Pollutant A.
Initially, cells 2-4 each have only primary pollutant, A, with
cell concentration equal to 1.0. The effects of advection for
a wind blowing four-tenths of a cell in time step At are simu-
lated first, resulting in the distribution shown in Figure 42b
131
-------
3SR-844
Cell I
0.0
0.6
1.0
1.0
0.4
0.0
Figure 42b — Distribution of Primary Pollutant A After
Advection.
If a quadratic chemical mechanism (A + A -»- B) is hypothesized
to illustrate the non-linearities in chemical kinetics rate
equations, the rate equations can be written as
Furthermore, if k is small, then k[AJ2 At « [AJ and the
effects of the chemical reaction on the concentration of A
can be ignored in order to focus on the secondary pollutant
production. With [A] approximately constant during the re-
actions the equation for production of [B] can be integrated
[B] = k[A]2 At .
That is, the production of the secondary pollutant is propor-
tional to the square of the concentration of the primary pollu-
tant. If, as an Eulerian technique would require, the cell
concentrations of A are used in this equation, the distribution
of B would be:
132
-------
3SR-844
Cell 1
0.0
0.36 kAt
1.0 kAt
1.0 kAt
0.16 kAt
0.0
Figure 42c - Distribution of Secondary Pollutant, B, Calcu-
lated from Cell Concentrations given in
Figure 42b.
On the other hand, the concentration of A inside the puff is
1.0 and outside 0.0. That is, in a Lagrangian frame the pro-
duction of B would be:
[B] = k[1.0]2 At = kAt ,
and would be distributed:
Cell 12 3 4 56
0.0
0.6 kAt
1.0 kAt
1.0 kAt
0.4 kAt
0.0
Figure 42d — Distribution of Secondary Pollutant, B, Calcu-
lated from Lagrangian Concentration.
This is the physically correct picture. The methodology used
in chemical simulation with Eulerian cell concentrations intro
duced errors of 40 percent in cell 2 and 60 percent in cell 5.
The random initial positions of the particles introduces even
further irregularities in the concentrations. The NEXUS/P ad-
vection test resulted in average concentrations for the center
27 cells of the 125-cell cube that were within 10 percent of
the correct values for the primary pollutants CNO, NO,., , and
HC) and within 20 percent for ozone after 40 cycles.
133
-------
3SR-844
A third test of NEXUS/P was conducted which also in-
cluded the effects of diffusion. In this test, the concen-
trations do change, even in Lagrangian coordinates, and so no
comparison to the CHEM results were made. Qualitatively, the
results were correct, but no adequate method for evaluation
was known.
These tests evaluated the changes made in NEXUS in
the development of NEXUS/P and demonstrated the difficulties
inherent in chemistry with Eulerian concentrations when sharp
gradients are present. A series of tests was also made with
an initial Gaussian distribution of pollutants with a = 2Ax,
.X
a = 2Ay and a = 2Az . These results are shown in Table VIII
y ' z
for the central cell concentrations and should agree with the
CHEM results. In these calculations, very few particles per
cell were used (5 particles per cell for primary pollutants
and an average of 1.5 particles per cell for Cu) and thus, the
comparisons represent the upper range of. the error introduced
by Eulerian chemistry for this smooth distribution. Even so
the errors in the central cell were within an acceptable
range O15 percent) and NEXUS/P was next applied to a real
photochemical smog problem.
4.4 INPUT DATA
The NEXUS/P code was used to simulate photochemical
smog in Los Angeles on September 30, 1969. This day was chosen
to correspond to data supplied by others. The grid used is
shown in Figure 43, a total of 22 cells in the x-direction
(cells 3 to 24) x 21 cells in the y-direction (cells 4 to 24)
x 4 in the z-direction with each cell 2 miles x 2 miles x
100 meters. Dickson and Start32 supplied surface winds. The
inversion height was obtained from Roth24 and the source
2 3
emissions inventory from Roberts, et.al.
134
-------
3SR-844
TABLE VIII
GAUSSIAN PUFF (QX = 2Ax, a = 2Ay, az = 2Az)
CENTRAL CELL CONCENTRATIONS
CYCLE POLLUTANT
0 NO
N02
HC
°3
5 NO
N02
HC
°3
10 NO
NO 2
HC
o,
CHEM
(Alone)
5
1
10
0
4
1
9
8.8
2
3
7
2.8
.3
.0
.0
.0
.4
.9
.3
x 10 3
.4
.6
.7
x 10 2
NEXUS/P
(No Winds
or Diffusion)
5
1
10
0
4
1
9
8.8 x
2
3
7
2.8 x
.3
.0
.0
.0
.4
.9
.3
io-3
.4
.6
.7
10-2
NEXUS/P
(Winds ,
No Diffusion)
5
1
10
0
4
1
9
9.3 x
2
3
7
2.6 x
.3
.0
.0
.0
.4
.9
.3
ID'3
.6
.2
.7
ID'2
135
-------
3SR-844
• The wind data from Dickson and Start consists of Los
Angeles APCD wind network observations interpolated, extrapo-
lated and smoothed to form hourly surface wind fields. The
inversion height data obtained from Roth are snown in Figure 43.
The inversion heights shown were used for the maximum heights
with inversion height rising linearly from ground level at
6:00 A.M. to the maximum height at noon. This is the same
procedure used in a NEXUS simulation of CO in Los Angeles. In
addition to the mobile source inventory shown in Figures 37-39,
Roberts compiled emissions from power plants, refineries and air
23
craft as well as the low level distributed sources. Hydro-
carbon emissions were separated into high and low reactivity.
However, Eschenroeder's chemistry mechanism only uses a single
hydrocarbon, so all hydrocarbons were lumped together. The
oxides of nitrogen data were not separated into NO and N02 so
an arbitrary 19 to 1 ratio was used for NO to N02 in emissions.
The initial concentrations of NO, N02, and HC were
interpolated from LA APCD observations at 6:00 A.M. The APCD
data were very sparse (only three stations reporting total
hydrocarbons), so a "1/r2" interpolation was used throughout
the grid. That is, the value at surface point x,y is given
by:
c.
f^f Cx - x.-)2 + Cy - y±
c(Xy) - 1 1 x x
m
£-• (x - x^2 + Cy - y±:
where the ith of m stations has concentration c..
The SETUP code for NEXUS/P is described in the NEXUS/P
operations manual, Volume III.
136
-------
3SR-844
? /«• /,-, if.. /7 l.V i
Figure 43 - Inversion Height (feet above sea level) at 12 noon
on September 30, 1969. (Reference 24)
137
-------
3SR-844
4.5 RESULTS OF THE NEXUS/P SIMULATION
The Los Angeles photochemical smog from 6:00 A.M. to
10:00 P.M. was simulated by NEXUS/P. The total computation
time on the S3 UNIVAC 1108 was less than one and one-half
hours for the 185 cycles. This was a five species (NO, N02,
HC, 0, and HN07J simulation on the 22 x 21 x 4 grid covering
O L
44 miles x 42 miles x 400 meters.
The prime concern of this study was to demonstrate the
feasibility of the PICK method for the solution of the atmos-
pheric diffusion equation and, in this section, for the simu-
lation of photochemical smog. The results of photochemical
smog simulation are very sensitive to the chemical mechanism
and source emissions inventory. Eschenroeder's photochemical
mechanism is still being developed and an expanded reaction
set was reported recently. The mobile source emissions in-
ventory was not satisfactory even for the CO simulation. For
these reasons, the results were given only cursory analysis
and comparison with LA APCD observations. In Table IX, the
calculated pollutant hour averages are averaged for the loca-
tions of all APCD stations reporting observations. This is
the closest possible comparison to a regional hour-averaged
concentration. The NEXUS/P results show a faster oxidation
of NO to N09 than is observed. The NO (the sum of NO and N0~)
^ X L
is removed too rapidly in the simulation by the product forming
reactions CPAN?HN02), Oxidant (CO concentrations usually re-
main low due to reaction 3, converting NO to N02> With NO too
low, 03 is too high. Both NEXUS/P results and APCD observa-
tions show a double maximum in N02; in the calculated results
at 7:00 A.M. and 2:00 P.M., and in the APCD observations at
9:00 A.M. and 1:00 P.M. The morning peak may be due to the
morning emissions being oxidized as the sun rises (the UV
radiation is required for reactions 1 and 12). The afternoon
138
-------
3SR-844
TABLE IX
»
REGIONAL AVERAGE OF STATIONS REPORTING
(observations are in parentheses)
TIME
6
7
8
9
10
11
12
1
2
3
4
5
6
:00 A.M.
:00
:00
:00
:00
:00
:00 Noon
:00 P.M.
:00
:00
:00
:00
:00
NO (pphm)
32
14
5
1
1
1
1
1
1
1
2
3
3
(33)
(30)
(22)
(11)
( 4)
( 3)
( 4)
( 5)
( 1)
( 2)
( 3)
( 3)
( 3)
NO 2
7
22
17
10
5
6
9
11
12
10
8
8
7
POLLUTANT
(pphm) HC (ppm)
( 6)
(10)
(17)
(18)
(14)
(ID
(15)
(18)
(16)
(10)
(11)
(10)
Cio)
8
8
7
6
6
5
5
5
3
2
2
2
1
( 7)
( 7)
( 4)
( 4)
( 4)
( 4)
( 5)
( 4)
( 3)
( 3)
( 3)
( 3)
C 4)
03 (pphm)
0
6
22
26
27
32
33
40
38
31
26
24
18
( 1)
( 2)
( 4)
( 6)
( 9)
(12)
(14)
(13)
(10)
( 8)
( 6)
( 4)
C 2)
13-Hr. Average 5 (10)
10 (13)
5 ( 4)
24 ( 7)
139
-------
3SR-844
maximum occurs because the only stations reporting are located
in high N09 regions. As expected, the NO results have the
Lt
proper time dependence — almost a monotonicly decreasing func-
tion from morning rush hours to evening rush hours, when new
emissions are large. The hydrocarbons (HC) also decrease until
evening rush hours. The poor comparison of average HC concen-
tration simulated by NEXUS/P at 6:00 P.M. may be due to the use
of only a three-station average as well as reactions with the
higher-than-observed 03 concentrations.
Only three LA APCD stations reported hourly averages
for HC. In Table X, the observations at these three stations
are compared to the NEXUS/P results. Even on this station-by-
station comparison the same conclusions can be reached: the
NO conversion to N09, removal of NO and oxidant formation are
L X
too rapid. For example, the peak in 0, simulated for the
Downtown LA station at 8:00 A.M. is related to the over-calcu-
lation of N02 at 7:00 A.M.
TABLE X
LOS ANGELES APCD STATIONS REPORTING MOST COMPLETE DATA
(observations are in parentheses)
Station: 1 Location: Downtown Los Angeles
POLLUTANT
TIME
NO (pphm) N09 Cpphm) HC (ppm) 0., (pphm)
6
7
8
9
:00 A.M.
;QO
:00
:00
48
21
2
2
(45)
(51)
(30)
UO)
5
41
24
3
C 2)
C 1)
( 7)
(11)
10
8
8
6
(10)
C 9)
0,1
8
50
23
C
C
(
(
1)
3)
7)
8)
140
-------
3SR-844
Table X, contd.
POLLUTANT
TIME
10:00 A.M.
11:00
12:00 Noon
1:00 P.M.
2:00
3:00
4:00
5:00
6:00
Average of
Hour
Reported :
Station: 60
NO
2
1
1
1
1
2
2
2
2
8
(pphm)
(10)
( 7)
( 3)
( 4)
( 2)
( 2)
(16)
Location
N02
5
12
17
7
7
9
10
14
15
14
(pphm)
(12)
( 8)
( 8)
(10)
( 8)
( 6)
( 7)
HC
5
4
4
4
3
1
2
2
1
4
(ppm)
( 4)
( 3)
( 2)
( 2)
( 2)
( 2)
( 2)
( 2)
( 4)
°3
29
45
48
67
51
28
36
37
23
33
(pphm)
(15)
(22)
(15)
(17)
(13)
( 6)
( 5)
( 4)
( 2)
( 9)
: Azusa
POLLUTANT
TIME
6:00 A.M.
7:00
8:00
9:00
10:00
11:00
12:00 Noon
NO
6
1
I
0
0
0
0
(pphm)
( 4)
( 8)
( 2)
( 1)
( 1)
( 1)
N02
8
12
7
5
4
4
4
(pphm)
C 6)
( 8)
(13)
(16)
(11)
( 9)
HC
8
7
6
6
5
5
5
(ppm)
( 8)
( 7)
( 5)
( 5)
( 5)
( 5)
( 6)
°3
1
7
14
21
20
23
21
(pphm)
( 1)
( 3)
( 7)
(13)
(17)
(18)
141
-------
3SR-844
Table X, contd.
TIME
1
2
3
4
5
6
:00 P.M.
:00
:00
:00
:00
:00
NO (pphm)
0
0 ( 1)
0 ( 1)
1 ( 1)
0
2 ( 1)
POLLUTANT
NO 2 (pphm) HC
4
5
14
13
10
11
( 4)
( 4)
C 4)
( 5)
C 7)
5
3
2
2
1
1
(ppm)
[
(
(
C
6)
4)
5)
4)
5)
6)
°3
31
24
45
37
30
9
(pphm)
(27)
(19)
(15)
(12)
(
C
7)
3)
Average of
Hours
Reported :
Station:
1 ( 2)
79 Location
8
( 7)
4
c
5)
22
(11)
: Pasadena
POLLUTANT
TIME
6
7
8
9
10
11
12
1
2
3
:00 A.M.
:00
:00
:00
:00
:00
:00 Noon
:00 P.M.
:00
:00
NO (pphm)
26 (23)
5 (19)
1 ( 3)
1 ( 1)
1 ( 1)
1 ( 1)
0
0
0
o c i)
N02
8
20
10
3
3
3
3
3
5
3
(pphm)
C 6)
C 9)
(11)
( 9)
U4)
(17)
C 5)
HC
8
8
8
6
6
6
5
5
4
3
(ppm)
C
C
(
C
c
c
c
c
c
c
5)
5)
3)
3)
3)
4)
5)
3)
2)
2)
°3
0
9
4
25
45
44
53
67
46
26
(pphm)
C
C
(
(
C
C
C
c
c
c
1)
1)
1)
1)
1)
2)
7)
7)
9)
4)
142
-------
Table X, contd.
3SR-844
POLLUTANT
TIME NO (pphm) N02 (pphm) HC (ppm) 03 (pphm)
4
5
6
:00 P.M.
:00
:00
1 (
1 I
1 1
: i)
: i)
: 2)
8
12
9
( 6)
(ID
(14)
2
2
2
( 2)
( 2)
( 3)
24
21
2
(
(
(
3)
2)
1)
Average of
Hours
Reported:
5)
8 (11)
5 ( 3)
33 ( 3)
The meteorological effects are similar to those dis-
cussed in Section 3.4 for the CO simulation. The photochemis-
try, a non-linear process, complicates the analysis, but
nevertheless the same conclusions about the effects of the
winds, inversion and horizontal wind convergences were evi-
dent in this simulation also.
The computer drawn isopleths calculated with NO and
N02 averaged over the initial hour simulated (6:00 A.M.) are
shown in Figures 44 and 45. At 6:00 A.M., the hydrocarbon
(HC) concentrations averaged above 5 ppm and oxidant (0,)
O
concentrations below 1 pphm and were both not plotted. By
10:00 A.M. (illustrated in Figure 46 - 48) the combination
of morning emissions and slowly increasing on shore winds
would have led to build-up of NO in the central basin, but
the chemical reactions have proceeded too rapidly converting
NO to N0?. NO is below 5 pphm and not plotted. N02 has also
been removed by the reactions and oxidants have been formed.
In Figures 46 and 47, N02 is below measured levels, but oxi-
dants are over predicted. In the afternoon the observations
143
-------
(I')
2.O
7
\
\
\
(37.)
N
1
oo
Figure 44 — NO concentration in pphm at 6:00 A.M.
-------
Cn
oo
Figure 45 — N02 concentration in pphm at 6;00 A.M.
-------
oo
Figure 46 — N02 concentration in pphm at 10:00 A.M.
-------
C/3
!*
i
oo
Figure 47 — 0, concentration in pphm at 10:00 A.M.
-------
oo
X
/ N_
\
fl-
A
\
7
•s.
\1
I
en
?o
i
oo
Figure 48 - HC concentration in ppm at 10;00 A.M.
-------
3SR-844
show N02 being depleted and so the NEXUS/P results are in
better agreement with observations (Figure 49). The oxidant
isopleths (Figure 50) are still too high, but show a spatial
distribution similar to the observations. The winds in the
afternoon have transported the major high levels of oxidants
toward the northeast. At 6:00 P.M. after the evening rush
hours, the NCU and HC show considerable similarities, but the
oxidants remain too high. The N02 distribution (Figure 52)
shows areas of higher concentrations near convergences in the
surface winds even though this is a secondary pollutant pro-
duced after conversion of emitted NO.
4.6 SUMMARY AND CONCLUSIONS
The NEXUS/P code was developed from NEXUS with the
addition of subroutines for chemical simulation and for
particle mass re-weighting to treat changes in concentration
calculated in the chemical simulation. The CHEM code was
developed to solve the chemical mechanism reported by Eschen-
9 Q
roeder and was 40 times faster than Eschenroeder's reported
technique.
The NEXUS/P was first applied to simulating a puff of
reactive pollutants. This simple application was used to
evaluate the compatibility of the PICK formalism as used in
NEXUS and the simulation of chemical reactions by the CHEM
code. In the first test, the simulated reactive puff was
stationary and the NEXUS/P results agreed with the CHEM simu-
lation alone. For the next test, the reactive puff was moved
with a constant velocity. At the boundaries of the puff,
errors were observed due to the use of Eulerian concentrations
for non-linear chemical simulation, as described in Section 4.3,
When concentration gradients were reduced by the use of a
Gaussian distribution for the advecting puff, the errors were
149
-------
tn
O
JO
i
oo
Figure 49 — NO- concentration in pphm at 3:00 P.M.
-------
I
oo
Figure 50 — 0, concentration in pphm at 3;00 P.M.
-------
on
t-o
\
\
\
\
7
\7
A
L
*-
O)
?3
I
oo
Figure 51 — RC concentration in ppm at 3:00 P.M.
-------
tn
0*)
\
\
\
\1
(«r
\
\
\
(
I* II \
vj
(w
V
fit
J—5
6
03)
en
oo
Figure 52 —NO- concentration in pphm at 6:00 P.M
-------
\
V
(3)
~ -+-
I
L
V
i
oo
Figure 53 — HC concentration in ppm at 6:00 P.M.
-------
Cn
Ui
cn
00
Figure 54 — CU concentration in pphm at 6:00 P.M.
-------
3SR-844
reduced. The magnitude of these errors for a realistic ap-
plication is not known at this time and should be investi-
gated. An advecting and diffusing reactive puff was also
simulated but a correct solution for comparison was not
known.
Photochemical smog in Los Angeles was simulated using
NEXUS/P. The 16-hour real time simulation required 1.5 hours
on a UNIVAC 1108. Approximately one-half hour was spent in
chemical simulation and one-half hour in data manipulation.
The photochemical smog was simulated using a mechanism de-
29
veloped by Eschenroeder with five calculated species (NO,
N02 , HC , 0., and HNCUJ and three species in pseudo-equilibrium
CO, R09 and OH). The source inventory used was reported by
23
Roberts, et.al, The required meteorological data was ob-
tained from Dickson and Start (winds) and Roth (inversion
height).
The results of the NEXUS/P simulation were compared
to Los Angeles APCD observations. In the simulation, NO was
converted to N02 too rapidly, total NO + N02 was depleted too
fast and too much 0, was produced. These deviations from the
observations occur basin-wide and are attributable to inac-
curacies in the chemical mechanism used in the simulation.
Some generally observed trends were reproduced — double maxi-
mum in N02 , NO + N02 and HC depletion during the day and 0.,
build-up and depletion. The calculated spatial distributions
of primary and secondary pollutants have the same general
features as observed but inaccuracies in the chemical mecha-
nism complicate any analysis.
156
-------
3SR-844
5, SUMMARY AND CONCLUSIONS
Tfie PICK particle^in^cell method for solution of the
atmospheric dispersion equation has been developed, tested
and applied to air pollution problems. The PICK method was
used in a two-dimensional computer code, PICFIC, for evalua-
tion of the method and comparison to finite difference
methods as applied to advection test problems.
In the simulation of advection, PICFIC is faster and
more accurate than the finite difference methods tested, and
PICFIC produced a solution for cases which the finite differ-
ence techniques failed to solve due to truncation errors.
PICFIC was also evaluated for accuracy of diffusion simulations
and for speed with different grid sizes and different numbers
of particles per grid cell. The accuracy of the PICFIC diffu-
sion simulation decreased with decreasing numbers of particles
per cell. Using ten particles per cell, the effective diffu-
sivity calculated with PICFIC was 0.90 of the actual value and
with three particles per cell was 0.60 of the actual value.
The computing time required was 30-50 micro-seconds per cell
per dimension and 100-300 micro-seconds per particle per di-
mension for each cycle on the S3 UNIVAC 1108. Thus, greater
resolution and accuracy may be gained at the least expense by
increasing the number of cells for a constant number of par-
ticles Cas long as .more than three particles per cell are
used). PICFIC was also used to simulate particulate matter
falling in a gravitational field from a line source in a shear
157
-------
3SR-844
wind; less than 10 percent deviation from the analytic solu-
tion resulted,
PICFIC was extended to three dimensions forming NEXUS.
NEXUS was developed with boundary conditions and input/output
routines appropriate for regional air pollution problems.
NEXUS was used to simulate carbon monoxide in Los Angeles.
The meteorological and source emissions data were obtained
from the literature. This first application had very good
correlation with Los Angeles Air Pollution Control District
observations. A correlation coefficient of 0.73 was obtained
by comparing NEXUS results with observations at 12 stations.
The 17-hour simulation took 12 minutes on the S3 UNIVAC 1108.
Tests were conducted with NEXUS to determine the effects of
different assumptions made for the boundary conditions, upper
levels winds and source emissions. The sensitivity of the
results to random errors in the meteorological and source
emissions data was also evaluated. A movie was made that de-
picts the flow of pollutants as simulated by NEXUS. In the
movie, dots are used to represent particle positions,
A second major application of the PICK method to an
air pollution problem was the simulation of photochemical
smog in Los Angeles, A technique for incorporating the ef-
fects of chemical reactions was added to NEXUS to form
NEXUS/P. A subroutine, CHEM, was developed to solve the chem-
ical mechanism reported by Eschenroeder and Martinez with
five calculated species CNO, N02, 03, HC and HNO-J . When
NEXUS/P was tested on an advecting puff of reacting pollu-
tants, errors were introduced by the use of Eulerian cell
concentrations in the non-linear chemical simulation. The
errors are small for small spatial gradients in concentration
as may be the case in some regional air pollution problems of
interest. For Los Angeles, the 16-hour photochemical smog
158
-------
3SR-844
simulation took 1,5 hours on the S3 UNIVAC 1108, The NEXUS/P
calculation shovred that NO is converted to N07 too fast and
Li
too much 0^ is produced, These features in the results are
attributed to the chemical reaction mechanism. The develop-
ment of a better chemical mechanism is central to an accurate
photochemical smog model. The errors introduced by the use
of Eulerian cell concentrations in chemical simulation should
be investigated to determine the significance of these errors
for photochemical modeling of actual air pollution problems
and to eliminate or minimize the errors.
This study has developed and demonstrated a new method
for simulating atmospheric advection and diffusion. The PICK
method has been shown to be more accurate and reliable than
the finite difference techniques for advection. PICK is the
preferred method for many air pollution problems where its
greater accuracy will result in a less expensive simulation
than finite difference methods using a smaller grid cell size,
The PICK method, being a hybrid Lagrangian-Eulerian method,
offers a unique capability to assess the importance of the
errors introduced by the use of Eulerian cell concentrations
in chemical simulation. Due to the accuracy of advection
with the PICK method, it may also make a significant contri-
bution to hydrodynamic modeling for meteorological applica-
tions .
159
-------
3SR-844
(This Page Intentionally Left Blank.)
160
-------
3SR-844
REFERENCES
1. A.A. Amsden, "The Particle-In-Cell Method for the Calcu-
lation of the Dynamics of Compressible Fluids," Report
No. LA-3466 (1966), Los Alamos Scientific Laboratory,
New Mexico.
2. J.E. Welch, F.H. Harlow, J.P. Shannon, and B.J. Daly,
"The MAC Method," Report No. LA-3425 (1965), Los Alamos
Scientific Laboratory, New Mexico.
3. F.H. Harlow and A.A. Amsden, "Fluid Dynamics: An Intro-
ductory Text," Report No. LA-4100 (1970), Los Alamos Sci-
entific Laboratory, New Mexico.
4. K.L. Calder, "On the Equation of Atmospheric Diffusion,"
Qtrly. J. Royal Meteorological Society (1965), 91, 390,
pp. 514.
5. O.G. Sutton, Micrometeorology, McGraw-Hill, New York (1953).
6. F. Pasquill, Atmospheric Diffusion, D. Van Nostrand, New
York (1962) .
7. J.E. Fromm, "A Method for Reducing Dispersion in Convec-
tive Difference Schemes," J. Computational Physics (1968),
3_, 176.
8. W.P. Crowley, "Numerical Advection Experiments," Monthly
Weather Review (1968), 96, 1, pp. 1.
9. W.P. Crowley, "Second-Order Numerical Advection," J. Com-
putational Physics (1967), 1_, 471.
ho® ro <-\f P
10. R.D. Richtmyer and K.W. M^ta^i^sz^ Difference Methods for
Initial-Value Problems, second ed.^Interscience Publishers,
New York (1967) . The first-order advection differencing
scheme described on pp. 292 is attributed to R. Lelevier.
11. B.A. Egan and J.R. Mahoney, "A Numerical Model of Urban
Air Pollution Transport," Conference on Air Pollution
Meteorology, AMS, Raleigh, North Carolina, April 5-9, 1971.
12. J.P. Pandolfo, M.A. Atwater and G.E. Anderson, "Prediction
by Numerical Models of Transport and Diffusion in An
Urban Boundary Layer," (1971), Center for The Environment
and Man, Hartford, Connecticut.
161
-------
3SR-844
13. R.C. Sklarew, "The Grid Model of Urban Air Pollution,"
APCA Paper No. 70-79 (June 1970), Systems, Science and
Software, La Jolla, California.
14. A.D. Taylor, letter to K.L. Calder, dated 26 May 1970.
15. J.L. Neuringer, SIAM J. Appl. Math. (1968), 16., 4, pp. 834.
16. R. Lamb, "An Air Pollution Model of Los Angeles," Masters
Thesis (1969), University of California, Los Angeles.
17. State of California, Air Resources Board, "Air Pollution
Control in California," Annual Report (1969), pp. 13.
18. A.Q. Eschenroeder and J.R. Martinez, "Mathematical Model-
ing of Photochemical Smog," AIAA Eighth Aerospace Sci-
ences Meeting, New York (January 1970), Paper No.70-116.
19. D.M. Teague, "Los Angeles Traffic Pattern Survey," Vehicle
Emissions I, SAE Progress in Technology (1964) , 6_, 102.
20. G.P. Larson, J.C. Chipman and E.K. Kauper, "Distribution
and Effects of Automotive Exhaust Gases in Los Angeles,"
Vehicle Emissions I, SAE Progress in Technology (1964) ,
^__
21. R.C. Sklarew, "Preliminary Report of the S3 Urban Air Pol-
lution Simulation of Carbon Monoxide in Los Angeles,"
(1970), Systems, Science and Software, La Jolla, California,
22. R.C. Sklarew, Systems, Science and Software, "Simulation
Models for Air Quality Control," IEEE Conference on Deci-
sion and Control, Miami Beach, Florida, December 1971.
23. P.J. Roberts, P.M. Roth and C.L. Nelson, "Contaminant
Emissions in the Los Angeles Basin — Their Sources Rates
and Distribution," Report No. 71 SAI-6 (1971) Systems
Applications, Inc.
24. P.J. Roth, private communication (1971).
25. M. Weisburd, L.G. Wayne, R. Danchick and A. Kokin, "De-
velopment of A Simulation Model for Estimating Ground
Level Concentration of Photochemical Pollutants," Report
No. TM-(L)-4673/000/00 (1971), System Development Corp.,
Los Angeles, California.
26. N.E. Browne, "Simulation Model for Air Pollution Over
Connecticut," J. Air Pollution Control Association (1969),
19, pp. 570.
162
-------
3SR-844
27. L.G. Wayne and T.E. Ernest, "Photochemical Smog, Simu-
lated by Computer," APCA Paper No. 69-15 (1969), New York
28. J.V. Behar, "Simulation Model of Air Pollution Photochem-
istry," Project Clean Air Research Reports 4_, (1970),
University of California. ~
29. A.Q. Eschenroeder and J.R. Martinez, "Concepts and Appli-
cations of Photochemical Smog Models,"
30. C.W. Gear, "A407 - DIFSUB for Solution of Ordinary Dif-
ferential Equations," (1971), Communications of the ACM
!U, 3, pp. 185.
31. T.J. Keneshea, "A Technique for Solving the General Reac-
tion-Rate Equations in the Atmosphere," Report No.
AFCRL-67-0221 (1967), Air Force Cambridge Research Labo-
ratories, Cambridge, Massachusetts.
32. C.R. Dickson and G. Start, private communication (1971).
33. J.R. Martinez, R.A. Nordsieck and A.Q. Eschenroeder,
"Morning Vehicle-Start Effects on Photochemical Smog,"
Report No. CR-2-191 (1971), General Research Corporation,
Santa Barbara, California.
34. R.C. Sklarew, "Computer Simulation of Carbon Monoxide in
Los iigeles," (1970), Systems, Science and Software, La
Jolla, California.
163
------- |