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^ ------- 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 |