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

-------