;PA-600/4-76-016c
  1976
Environmental Monitoring Series
     CONTINUED RESEARCH IN MESOSCALE All
             POLLUTION SIMULATION  MODELING:
                            Volume  III  -  Modeling
                       of Microscale Phenomena
                               Environmental Sciences Research Laboratory
                                   Office of Research and Development
                                  U.S. Environmental Protection Agency
                             Research Triangle Park, North Carolina 27711

-------
                 RESEARCH REPORTING SERIES

Research reports of the Office of Research and Development, U.S. Environmental
Protection Agency,  have been grouped  into five series. These  five broad
categories were established to facilitate further development and application of
environmental technology. Elimination of traditional grouping was consciously
planned to foster technology transfer and a maximum interface in related fields.
The five series are:

     1.    Environmental Health Effects Research
     2.    Environmental Protection Technology
     3.    Ecological Research
     4.    Environmental Monitoring
     5.    Socioeconomic Environmental Studies

This report has been assigned  to the ENVIRONMENTAL MONITORING series.
This series describes research conducted to develop new or improved methods
and  instrumentation  for the  identification and quantification of environmental
pollutants at the lowest conceivably significant concentrations. It also includes
studies to determine  the ambient concentrations of pollutants in the environment
and/or the variance of pollutants as a function of time or meteorological factors.
This document is available to the public through the National Technical Informa-
tion Service, Springfield, Virginia 22161.

-------
                                           EPA 600/4-76-016 C
                                           May 1976
     CONTINUED RESEARCH IN MESOSCALE AIR
       POLLUTION SIMULATION MODELING:
VOLUME III - MODELING OF MICROSCALE PHENOMENA
                 R. G. Lamb

     Systems Applications, Incorporated
             950 North gate Drive
        San Rafael, California  94903
                 68-02-1237
               Project Officer

            Kenneth L. Demerjian
     Meteorology and Assessment Division
 Enviromental Sciences Research Laboratory
Research Triangle Park, North Carolina  27711
    U.S. ENVIRONMENTAL PROTECTION AGENCY
     OFFICE OF RESEARCH AND DEVELOPMENT
 ENVIRONMENTAL SCIENCES RESEARCH LABORATORY
RESEARCH TRIANGLE PARK, NORTH CAROLINA  27711

-------
                                                                  11
                           DISCLAIMER
     This  report  has  been  reviewed by the Office of Research and
Development,  U.S.  Environmental Protection Agency, and approved for
publication.   Mention of trade names or commercial products does not
constitute endorsement or  recommendation for use.

-------
                                                                          m
                                CONTENTS
DISCLAIMER
LIST OF ILLUSTRATIONS	vi

LIST OF TABLES	ix

ACKNOWLEDGMENTS 	    x

   I  OVERVIEW  	    1

  II  ESTIMATION OF EDDY DIFFUSIVITY, WIND SHEAR PROFILES,  AND
      DISPERSION PARAMETERS IN THE PLANETARY BOUNDARY LAYER 	   12

      A.  Introduction	12
      B.  Description of the Numerical Experiments   	   16

      C.  Assessment of the Atmospheric Diffusion Equation   	   19
          1.  Wind Profiles from the Planetary Boundary Layer Model  ...   21
          2.  Theoretical  Forms for Eddy Diffusivities  	   24
          3.  Estimation of the Optimal Ky(z)  in the
              Atmospheric Diffusion Equation  	   25
      D.  Analysis of the Dispersion Parameters az  and ax
          in the Plume and Puff Models	34

      E.  Implementation of the Optimal Diffusivity and
          Wind Shear Profiles in an Airshed Model 	   50

 III -DEVELOPMENT OF A CLOSURE APPROXIMATION FOR THE CONCENTRATION
      FLUCTUATION TERMS IN THE GOVERNING EQUATIONS   	   58

      A.  Introduction	58
      B.  Derivation of an Approximation for Concentration
          Fluctuation Terms Such as A'B~r	62

      C.  Formulation of the Parameters c. and Cg	75
          1.  Release of A into a Uniform Field of-B	75
          2.  Release of Dilute Mixtures of A and B
              into a Homogeneous Turbulent Reactor   	   79
          3.  Release of Rectangular Slugs of Reactants
              into a Hypothetical, One-Dimensional  Reactor   	   81

      D.  Formulation of the Parameter A	-82

      E.  Testing the Closure Scheme Represented by Eq. (95)
          Using Observational Data	89

-------
                                                                           IV
IV  DEVELOPMENT OF A SCHEME  FOR PARAMETERIZING  THE  EFFECT  OF
    SUBGRID-SCALE  CONCENTRATION VARIATIONS  ON REACTION  RATE   	   98

    A.   Introduction	98

    B.   Statement  of the Problem   	98

    C.   Conditions Under Which  SSV  Effects  on Reaction  Speed
        Are Negligible   	100

    D. '  Parameterization Scheme for the  Subgrid-Scale
        Concentration Variations 	  106

    E.   Derivation of rAB in Terms  of 'Measurable  Parameters and
        Application of the Parameterization Scheme  to Specific Problems.  Ill

        1.   Random Distribution of  Point Sources  in a
            Three-Dimensional Space  	  114
        2.   Nonrandom Distribution  of Sources in  Three  Dimensions   .  .  .125
        3.   Isolated Point and  Line Sources	134

 V  DEVELOPMENT OF A SUBMODEL FOR RESTORING POINT SPATIAL  RESOLUTION
    TO GRID MODELS OF URBAN  POLLUTION	138

    A.   Introduction	138

    B.   Quantitative Illustration of the Subgrid-Scale
        Variations c"(r,t)	141
    C.   Derivation of a Submodel of c"  for
        Linearly Reactive Pollutants 	  146

        1.   Discussion of General Operational Problems  	  146
        2s   Derivation of the Working Equations for a Microscale
            Model  of Linearly Reactive  Pollutants  	  148
        3.   Design of the c" Model  Computer Program—
            Interaction with the Source  Inventory Data  Bank	153

    D.   Development of Mathematical Methods for Modeling Subgrid-Scale
        Concentration Variations c" of  Nonlinear  Pollutants   	  162

        1.   Derivation of a  Discrete-Time,  Continuous-Space
            Diffusion Equation  	  163
        2.   Development of a Microscale  Model of  Nonlinear Pollutants
            Based  on the Discrete-Time,  Continuous-Space Equation   .  .  .  168

VI  SIMULATION OF  BUOYANT PLUMES IN THE  PLANETARY BOUNDARY LAYER  .  .  .  .175

    A.   Introduction	175
    B.   Calculation of the Probability  Density  p(r,t|r',t')	176

    C.   Mathematical Basis of the Trajectory Simulation Model-
        Preliminary Design	178
        1.   The Particle Momentum Equation  	  178
      .  2.   The Particle Temperature Equation for Dry Plumes  	  183
        3.   The (T) Equation	186

-------
 VII   SUMMARY	188
      A.   Chapter I	188
      B.   Chapter II	189
      C.   Chapter III	191
      D.   Chapter IV	191
      E.   Chapter V	193
      F.   Chapter VI	194
VIII   FUTURE EFFORTS	195
      A.   Specification of the  Degree  of Spatial  Resolution
          Required in  Urban Pollution  Modeling  Studies	195
      B.   Evaluation  of the Actual  Spatial  Variability  of
          Concentrations at Receptor Sites  of Interest	196
      C.   Assessment  of the Need for Refined Mircroscale
          Transport  and Diffusion Formulas	198
APPENDICES
      A   FORTRAN PROGRAMS FOR  USTAR,  DK7,  AND  UBAR	199
      B   DERIVATION  OF EQS.  (193)  AND (195)	203
      C   FORTRAN LISTINGS OF THE CALC,  SIGMAX, AND  SIGMAZ FUNCTIONS.  .  .208
      D   THE MONTE  CARLO TECHNIQUE USED BY THE CALC SUBPROGRAM	212
      E   DERIVATION  OF rflR,  ?., AND ?R  FROM MULTIJET REACTOR
          DATA AND TOOR'SrtTHEORY.	218
      F   ESTIMATION  OFJHE ORDER OF MAGNITUDE  OF THE TIME SCALE
          OF DECAY OF  A"2	225
REFERENCES	230
FORM 2220-1	233

-------
                               ILLUSTRATIONS
 i


 2


 3

 4


 5


 6
Illustration of the Relationships Among the Solutions of
Various Governing Equations and Measurable Pollutant Concentrations

Illustration of the Relationships Among the Solutions of Various
Governing Equations and Measurable Pollutant Concentrations .  .  .  .
Schematic Summary and Index of Microscale Research
The Probability Density Function p(£,z,T/Z<---in Relative Units —
Derived From the Model  of Deardorff 	
Wind Profiles Given by Similarity Theory
for Neutral (Dashed Curves) and Unstable
(S)
and Deardorff's Model  (D)
 = 4.5 Conditions. .  .  .  .
 9

11


18


22
Comparison of the Optimal  K/ Profile (Solid Curve) with the Profile
Proposed by Shir and Shieh (Dash-Dot Line), the Profile given by
Similarity Theory (Dotted  Line), and the^Eddy Viscosity Calculated by
Deardorff (Dashed Line) in the Neutral  (h/L = 0) Boundary Layer .  .  .
                                                                               29
 7   Comparison of the Optimal KZ Profile (Solid Curve) with the Profile
     Based on Similarity Theory Proposed by Ragland (Dashed) in the Unstable
     h/L = -4.5 Boundary Layer	30

 8.  Comparison of Optimal Profile and Numerica-Empirical Solutions
     for the Neutral Case	32

 9   Comparison of Optimal Profile and Numerico-Empirical Solutions for
     the Unstable Case	33

10   Comparison of the Computed Nondimensionalized Root-Mean-Square Vertical
     Particle Displacement oz (Circles and Triangles) with the Pasquill-
     Gifford az Dispersion Parameter (Curves) for Neutral and Unstable
     Conditions as a Function of Downwind Distance from the Source 	 36

11   Comparison of the Computed Nondimensionalized Root-Mean-Square Vertical
     Particle Displacement az (Circles and Triangles) with the McElroy-
     Pooler az Dispersion Parameter Under Neutral and Unstable Conditions
     as a  Function  of Travel  Time t	38

12   Comparison of the Computed Nondimensionalized Streamwise Root-Mean-
     Square Particle Displacement ax (Circles and Trinagles) with the McElroy-
     Pooler Lateral  Dispersion Parameter ay for Neutral and Unstable Conditions
     as a Function of Travel  Time t	39

13   Comparison of the Predictions Made by the Puff Model and the Plume Model
     for the Case of Neutral  Stability	4]

14   Comparison of the Predictions Made by the Puff Model and the Plume
     Model  for the Unstable Case   	42

-------
                                                                             VII
15   Optimal aL (t) Profiles for use in the Puff Model Compared with the
     Empirical Data of McElroy and Pooler ..................  45

16   Optimal Streamwise Dispersion ox (t) Compared with the Empirical
     Lateral Dispersion Data of McElroy and Pooler ..............  46

17   Optimal az(x) Profiles for Use in the Plume Formula Compared with the
     Empirical Data of Pasquill and Gifford .................  47

18   The Fractional Error (lOOe) Resulting in the Puff Model and the Plume
     Model Using the Optimal o Profiles Shown in Figures 15 through 17 for the
     the Neutral Case ............................  48

19   The Fractional Error (lOOe) Resulting in the Puff Model and the Plume
     Model Using the Optimal a Profiles Shown in Figures 15 through 17 for
     the Unstable Case ............................  49

20   The Initial Slugs of Species A and B in a One-Dimensional Vessel ....  63

21   Fractional Error e in the Donaldson-Hilst Approximation £EQ.(59f)]  of
     Two Premixed Slugs of Material in a One-Dimensional Reactor .......  71

22   Mixing Turbulent Clounds of Reacting .Particle Species A and B ......  71
23   Comparison of Concentration Profiles in the Mixed Zone v/^g for Inert
     and Reactive Species ..........................  83

24   Functional Form of A(t) for Various Values of the Ratio a of the
25
26

27

28
29

30
31
32
Dependence of xmi-n on a for the Case of the Feed Ratio 8=1 	
Impact of Variations in the Feed Ratio 3 on the Functional Form of
X(t) for a = 700 	
Influence of Changes in the Initial Number and Separation of Reactant
Slugs on the Functional Form of A 	
First Comparison of Closure Model Predictions with Observational Data
Attempts to Simulate the Observational Data Using a Rate Equation that
Omits Concentration Fluctuation Effects 	
Second Comparison of Closure Model Predictions with Observational Data
System of Rectangular Sources 	
Description of the Volumes VA, VB, and VAB 	 : .
88

, 90

91
. 94

. 96
. 97
, 104
. 106
33   Ranges of the Vectors r" and r" +Ar Given Ar and Given r" Lies in v-j .  . 117

-------
34   Influence of Subgrid-Scale Concentration Variations  on
     the Chemical Rate of Change of A	123

35   Variation of the Dimensionless Parameter y as  a Function  of Rate
     Constant k for Various Atmospheric Stabilities 	  126

36   Uniform Network of Streets Illustrating the Calculation of the
     Stock Correlation Function G	129

37   Illustration of the Subgrid-Scale Concentration Variations Arising
     From an Instantaneous Point Source 	  144

38   Parameters Required to Specify a Finite Line Source   	  151

39   Flow Diagram of the c"-Model Routine Micro	155

40   The Contours of c" Arising From a Single Point and  Line Source,  as
     Computed by CALC	157

41   Horizontal and Vertical Cross Sections of the  Domain of Influence
     on c"(XR,YR)	159

-------
                          TABLES
Wind and Vertical  Eddy Diffusivity Profiles  in
the Surface Layer from Similarity Theory  	  23

Optimal  Parameters for Use  in  the Diffusion  Equation and in the
Gaussian Puff and Plume Models (Nondimensional Variables)  ....  52

-------
                        ACKNOWLEDGMENTS
     The calculations  of the  optimal diffusivity profiles reported in
Chapter II were performed in  collaboration with Drs. W. Chen and J. H.
Seinfeld of the Department of Chemical Engineering, California Institute
of Technology.   The author also wishes to thank Mr. Winston Shu of the
California Institute of Technology  for his assistance with the theory
reported in Chapter III and for making available excerpts from his
doctoral thesis, which is presently in preparation, for use in this
report.

-------
                             I    OVERVIEW


     In this report, the term "microscale" refers to all  phenomena that
have characteristic temporal  or spatial  scales that are too small  to be
resolvable in an explicit, deterministic.manner in urban-scale air pol-
lution simulation models.   Turbulence, concentration fluctuations, and
subgrid-scale concentration variations are examples of microscale  phe-
nomena.  In this chapter,  we  outline theoretically the way in  which each
of these important microscale features arises; in subsequent chapters  of
this volume, we propose techniques  for dealing with each  in air pollution
modeling studies.

     Consider the case of two pollutant species,  A and B,  that undergo the
reaction

                                 A  + B -> C

The fundamental equation governing  the concentrations of  each  species  is
the so-called mass continuity equation,  which in  the case  of A has the
form
where
          A    =  the instantaneous concentration of A at the point
                  iX-i a  Xn >  Xo»  ^/  >
          u.   =  the i-th  component of the instantaneous fluid
           1      velocity  at the same point,
          S.   =  the source strength (mass/volume/time),
           M
          K    =  the molecular diffusivity (of both A and B),
          k    =  the rate  constant of the chemical  reaction.

-------
     In an application of Eq.  (1) to a simulation of A in the atmosphere,
information concerning the instantaneous fluid velocity, u, is available
only from a few widely scattered meteorological  stations.  It can be shown
using the sampling theorem that from such a finite set of discrete obser-
vations, only the features of the velocity field with spatial scales
larger than about twice the average distance between sampling points can
be described explicitly [see Lamb and Seinfeld (1973)].  All  smaller fea-
tures are unresolvable and must be treated as stochastic variables.   Thus,
in an atmospheric simulation,  u. must be expressed in the form

                             u1 = Di + u:    ,

where u. is the mean velocity component, which embodies all of the large,
resolvable features of the wind field, and u! represents all  of the  unre-
solvable features.  The latter, which is usually called the turbulent
velocity component, represents a microscale phenomenon according to  our
definition.  Because of the stochastic nature of u!, both A and B are sto-
chastic; consequently, only their averaged properties have meaning.   In air
pollution studies, one of the most important averaged properties of  the con-
centration is the mean, or first moment, which we denote by angular  paren-
theses.  Taking the mean of Eq. (1), we find that <^A^> is governed by the
equation (assuming an incompressible atmosphere)
                                                                         (2)
In this equation, we have taken <( }to be an ensemble average,  which can be
regarded loosely as a time average.   The quantities A1  and B1  represent in-
stantaneous departures of the concentrations of species A and  B,  respectively,
from each one's. mean values.  Induced by the velocity fluctuations u!,  both
A1 and B1 are unresolvable and have  definitions similar to that of u!:

                               A1 =  A - 
                                                                         (3)
                               B1 =  B -

-------
     It is readily apparent from Eq. (2) that the microscale velocity vari-
ations u! affect the mean concentration <"A> both directly, through the
transport term , and indirectly, through the concentration fluctuation
term.  Both of these effects must be expressed in terms of (A), (B),
or other known variables before Eq. (2) can be solved,  because otherwise the
equation contains three unknown quantities:  , /u'A'\ and /A'B'V

     In the case where u! is caused solely by turbulence in the planetary
boundary layer, the most well-known method for approximating the fluctuating
transport term <(u.jA'N is the so-called gradient transport approximation, pro-
posed originally by Boussinesq in the nineteenth century.  Under this hypothesis,

                              " i   / ~ " ij 8x.i
where K.. is an empirical quantity called the turbulent diffusivity tensor.
        • J
Theoretical expressions (based on similarity theory)  are available  for some
of the components of K.., but these are restricted to the surface layer,
                      ' v)
which usually represents only a small  portion of the  layer into  which  pollu-
tants mix under typical  urban atmospheric conditions.  In Chapter II,  we
discuss this diffusivity in more detail and present new expressions for the
diffusivity—applicable to the entire  mixed layer—which we derived using a
new approach.

     In cases where u! is caused both  by ambient boundary layer  turbulence
and by perturbations of other origins, such as buoyancy forces  (when the
pollutant is hot) or turbulence generated in the wake of a large building
or hill, methods other than the simple transport theory [Eq.  (4)] must be
used.  The well-known plume rise formulas are examples of attempts  to  compen-
sate for buoyancy effects on u!.  In Chapter VI, we outline a technique for
obtaining refined estimates of the transport term<^u!A'^> in problems where
buoyancy forces play an important role.

-------
     Mathematical approximations for the fluctuation term/A'B/> have not
been used heretofore in urban-scale air pollution models; such terms have
simply been omitted.  It is known, however, that this simplification is
unjustified, especially for small-scale modeling and for fast chemical
reactions.  In Chapter III, we discuss some of the approximations that
previous investigators have proposed for these terms, and we develop and
test our own formulation.   For the time being, let us write

                            ~-  FAB    .                             (5)

Using Eqs. (4) and  (5), we obtain the model equation:
It is important to realize that this equation is actually a model of the
"true" equation [Eq. (2)J because it contains approximations of the terms
  to distinguish it from the true mean value
These two quantities are related as follows:
                                                                         (7)
where e is an error term whose value is a function of the model approxima-
tions used for /u!A'^> and<(A'B^> .   Note also that we have dropped the
molecular diffusivity term from Eq. (6) because it is known that relative
to turbulent diffusion, molecular diffusion is negligible insofar as the
mean <^A) is concerned.
     Figure 1, which summarizes the analyses presented up to this point,
     trates how A, <^A)>, and , and   would relate to actual data gathered from

-------
                                  • SOURCE

                                  MONITORING STATIONS'
                                                                  A,(t)
                   u
     ENSEMBLE OR TIME AVERAGING
     ASSUMPTIONS:   - - K..
                  i        ij
                     « K
                       rAB
                      See Figure 2
(a) Observed Instantaneous Concentra-
   tion at Stations  1 and 2
                                                                 ,Vtj
                                                      (b) Observed Mean Concentration
 (c) Analytic Solution of the
    Model Equation
FIGURE 1.   ILLUSTRATION  OF THE RELATIONSHIPS  AMONG THE SOLUTIONS OF  VARIOUS
         GOVERNING EQUATIONS AND MEASURABLE  POLLUTANT  CONCENTRATIONS

-------
    Since Eq.  (6) and its counterpart governing   are closed,  and hence
solvable, they can be used as the basis of an air pollution model.   How-
ever, because  these equations are nonlinear and because they contain vari-
able coefficients, numerical  integration methods must be used to  obtain
approximate forms of their solutions.  To implement numerical techniques
using a digital  computer, we  must first discretize the independent  variables.
But as we pointed out earlier, a finite set of samples of a continuous
function is sufficient to describe only certain large-scale features of that
function.  Consequently, before Eq.  (6) can be integrated numerically, it
must first be  "filtered" to remove all  small-scale variations, both in 
and in the independent parameters, such as K.. and u., that the grid system
cannot resolve.   The necessary filtering can be accomplished by averaging
the equation at each point over a volume equivalent to that of the  grid cells
and over a period, At, equal  to that of the time step used.

    Let the required space averaging operation be denoted by the  tilde (~)
and be defined by
                            x+Ax       y+Ay     /-z+Az
                          fff     ffj'.t)"*'  I*'  dz'
                          VAX       y-Ay     "Z-AZ
where x = (x,y,z) and where AX, Ay, and AZ represent the half-widths  of the
grid cell faces.  Averaging Eq. (6) in this manner and assuming,  as  is  usually
true in practice, that neither u.  nor K. .  possesses spatial  variations  of a
                                •        * \j
scale smaller than the grid, we obtain
                                                                         (9)
where
                               A" = m - m    ,
                                                                         (10)
                               B" =

-------
     The variables A" and B" represent variations in <(A>  and <(B)> ,  respec-
tively, that are of a scale smaller than the grid size.   Accordingly, A"
and B" constitute a new microscale phenomenon with which we must deal in
the development of an air pollution model.   Note that A" and B" affect the
chemical behavior of <^A)>  in much the same  way in which  the turbulent con-
                                                                        m
centration fluctuations A'  and B1  influence  [compare Eqs.  (9)  and (2)].
Hereafter we refer to A" and B" as subgrid-scale concentration variations
(SSVs).

     To complete our filtering of Eq. (6), we must next perform a  time
averaging of Eq. (9) over the period At of one time step.   In  practice,
the size of the time step At required to maintain computational  stability
and to minimize truncation errors  in the numerical integration of  the model
equations is much smaller than the characteristic time scale of any of the
variables entering in Eq. (9).  For this reason, the time  average  of <(A
is equivalent to<^A)>  itself.  The same holds for Sfl, K.., FflR,  and u..
                    m                              M   1J    AD       1
Consequently, we can regard Eq. (9) as having been averaged in both space
and time and as therefore representable in the required discrete form.

     Like its analogue [Eq. (2)],  Eq. (9) is not solvable  because  it con-
tains, in addition to the dependent variable <(A/> , the unknown quantity
A"B".  To our knowledge, this term has never before been considered in air
pollution modeling studies.  Apparently, this omission is  due  not  to the
belief that this term is small, but rather to a lack of awareness  of the
existence of the term.  In Chapter IV, we consider this SSV term in more
detail, both qualitatively and quantitatively, and we develop  a  mathematical
method for representing it in pollution models.   For the time  being, let us
write
                             G - A"?'    ,                               (11)

where G denotes a known function.   Equation (9)  now takes  the  form

-------
     Through a series of operations and approximations, each motivated by
necessity, we have arrived at an equation [Eq. (12)] that can feasibly be
used as the basis for a simulation model of urban air pollution.  A natural
question at this point is whether the solution <7V>  of this equation pro-
vides an adequate representation of the desired mean concentration .  To
answer this questions, we plotted in Figure 2(b)  the values of <^A>  that
Eq. (12) would yield at the two sample sites under the conditions of the
problem depicted in Figure l(a) [and 2(a)J.  After comparing these values
with the corresponding values of <(A^> shown in Figure l(b), one might con-
clude that   bears no resemblance to the desired quantity  and that,
consequently, an air pollution model based on Eq. (12) would be of no value
whatsoever.  However, a totally different picture emerges when the global
features of <(A>  are compared with those of .   Figure 2(c), which presents
this comparison, demonstrates that on scales that are large compared with the
grid scale, a positive correlation exists between <(A>  and ;   repro-
duces the overall features of the regional distribution of <(A^ quite well.
Only on scales comparable to and smaller than the grid size does Eq. (12)
break down.  This fact raises a second question:  What degree of spatial
resolution is required of an airshed model?

     The answer to this question must be inferred from the nature of the
task for which the model is intended.  For example, one such task is the
delineation of long period trends and regional  patterns in urban air qual-
ity.  For this purpose, a model based on Eq. (12) should suffice.  However,
models are also required that assess whether particular control  strategies
can achieve air quality standards.  To serve in this capacity, an airshed
model must possess point spatial resolution, since the standards for some
pollutants are currently expressed in terms of specific, short-period con-
centrations.   Thus, techniques for simulating urban pollution are required
that can reproduce the entire spectrum of spatial variations in the concen-
tration field.

     The basic model  equation [Eq. (6)] from which Eq. (12) was derived
possesses the desired range of resolution; but, as the reader will recall,
the integration of this equation with point resolution over an entire urban

-------
                                 Horizontal  Boundaries of One Grid Cell
                      I
             Continued from Figure 1
1
3 ' B + IT

»<*>«,
3X1

j 3Xj * SA - km - kFAB


                Averaging Over Time Step
                    and Grid Cell
Tt~TUi
                    S. - k 
                     A     mm
                                                                     t)\,
                                              (a) Same as Figure l(c)
                                                            (b) Model Equation Adequate for Com-
                                                               puter Simulation and Predicted
                                                               Values at Stations 1 and 2




	 i



/t
TT*
I ^
SOURCE

f
//
'
•^ISOPLETHS OF m

'STATION 1
XSTATI
i


ON 2


GRID CELLS

^
^
                                                             (c) Spatla) Distributions of
                                                                  and 
TIG'JRE 2.   ILLUSTRATION  OF THE RELATIONSHIPS AMONG THE SOLUTIONS OF VARIOUS
         GOVERNING EQUATIONS AND MEASURABLE  POLLUTANT  CONCENTRATIONS

-------
                                                                             10
region would require a prohibitive amount of computer time and storage.
This constraint is what motivated the spatial averaging, or small-scale
filtering, process that led to Eq. (12).  In Chapter V, we develop what
we call a subgrid-scale model that can be used in conjunction with an
airshed model based on Eq. (12) to obtain point resolution—in particular,
the quantity  --at any arbitrary point.  An alternative to the approach
presented in Chapter V (which we do not djjcuss in this report) is to cal-
                       —                s~~y
culate in addition to <(A>  the quantity A  .   By virtue of the definition
                         rn    i- *^~t)~\ i
of A" [Eq. (10)], the quantity [A" _r is a measure of the amplitude of de-
viations of <(A)>  f rom <(A)>  that one can expect within the grid cell where
/f~> o            ni         in
A"  is evaluated.  From the modeling standpoint, this approach has certain
operational advantages, but it does not provide the preciseness that is
often required.

     In summary, our research into microscale phenomena has focused on'.all
aspects of air pollution modeling that are attributable to mechanisms whose
spatial or temporal scales are too small to be resolvable in a simulation
model of an urban-scale region.   Figure 3 summarizes these research areas
schematically.  Because of the broad scope of this project, our efforts have
not yet produced working solutions to all of the problems that we have con-
sidered.  Our primary concern in this portion of our contract effort has
been to explore each of the microscale problem areas in some detail, to
assess the relative importance of each to the overall problem of air pollu-
tion modeling, and to design techniques for dealing with problem areas as
needed.  Only the implementation of our proposed techniques remains.  We
present a technical summary of the work described in this volume in Chapter
VII and a discussion of the work and problems remaining in Chapter VIII.

-------
                                  2
                                 ~~— + S. - KAXB>
                                 3X.3X.    A
     Chapter VI
     u]  Influenced
     by  pollutant
     buoyancy forces
    New  equation
3       3
                      Chapter  II

                      u{ due only to
                      boundary layer
                      turbulence
                  ,
                • "37
                              3
                                     SA - km - kFAB
                        Space Averaging
                                                               Fundamental
                                                               Equation
Fundamental
Model  Equation:
Point  Spatial
Resolution
     3_      3

                                   - kF
                                       AB



3 3m . . . \
If . , m. + «; j. Kfl> 
-------
                                                                            12
   II    ESTIMATION OF  EDDY  DIFFUSIVITY,  WIND SHEAR PROFILES,
         DISPERSION PARAMETERS  IN  THE  PLANETARY BOUNDARY LAYER
A.   INTRODUCTION

     The term u'.A'  entering  into  Eq.  (2)  represents the mean flux of the
material A arising  from turbulent fluctuations in the fluid velocity field.
The physical  significance  of this flux can be illustrated by the following
example.

     Suppose  that a fluid  is placed  in a  closed container and that turbulent
motions are subsequently induced  in  the fluid by some mechanical means.   Sup-
pose further  that numerous infinitesimal  velocity and concentration probes
are placed uniformly throughout the  fluid.  Let u (x,t) denote the functional
form of the expression that  best  describes the instantaneous signals gathered
from all the  velocity probes,  and similarly let c(x,t) denote the expression
that best represents the concentration measurements.  Let us assume that at
some instant  t~ a steady point source of  unit strength of the material  c is
introduced into the fluid  at the  point xn, and that no matter how many  times the
experiment is repeated, the  measured field c(x,t) is always found to be identical
to the solution of the continuity equation
where u. is the measured  velocity  field and where molecular diffusion does
not affect c.

     Suppose that many of the  velocity probes are now withdrawn and that the
velocity fields sensed by the  remaining probes is described by the function
U(x,t).   If the above  experiment were repeated, one would find that the
measured concentration field   c(x,t) spreads out more than is predicted by
the equation

-------
                                                                             13
A greater degree of spreading would be observed if still  fewer velocity
probes were used to obtain U.  This apparent spreading is the phenomenon
known as turbulent diffusion.

     Our example illustrates that turbulent diffusion is  actually an arti-
fact of our lack of knowledge of the true velocity field.  The example also
suggests that the magnitude of the diffusion phenomenon is related to the
density of velocity probes with which the fluid motions are monitored.
Since this aspect of diffusion and its potential significance to pollution
modeling have been discussed by Lamb (1971), and we do not elaborate upon
them here.  The primary concern in this section is to develop a mathematical
description of the diffusion phenomenon caused by turbulence in the planetary
boundary layer.

     There are basically two approaches to deriving such  a description—
Eulerian and Lagrangian.  The starting point of the Eulerian method is the
mass continuity equation:
where S is a known function describing the distribution of sources of the
scalar quantity.  This same equation was introduced above but in a slightly
less general form.  In a turbulent fluid, such as the atmosphere, the lack
of a total description of the velocity field u.  makes it necessary to treat
u.--and hence c--as stochastic variables.  These are customarily decomposed
into mean (and) and fluctuating (u! and c1) components.  Then, from
Eq. (13), one can see that  is governed by the equation
which is unsolvable because of the presence of the additional  unknown quantity

-------
                                                                             14
     This problem, known as the closure problem, is a fundamental  obstacle
that impedes the progress of all Eulerian-type approaches to the analytical
study of turbulent diffusion.   One of the oldest and best known attempts to
circumvent the closure problem is the so-called gradient transport hypothe-
sis, which dates back to Boussinesq in the nineteenth century.   In the con-
text of Eq. (14), this hypothesis states that
where K. - is a turbulent diffusivity tensor that must be evaluated using
        ' J
empirical data.  In Volume II, we discussed some of the attempts  made by
previous investigators to relate certain components of K. .  to measurable
                                                        ' J
properties of the boundary layer.  We review some of these  later  in this
chapter.

     Since Eq. (15) implies that turbulent diffusion behaves  like molecular
diffusion, the validity of this equation cannot be universal.  In recent.
years, several more advanced closure schemes have been developed, primarily
for the turbulent momentum and energy equations (see, for example, Hanjalic
and Launder, 1972; Lumley, 1970; Donaldson, 1969);  'but. none has yet
been demonstrated to be widely applicable to problems of atmospheric diffu-
sion.

     In contrast, the Lagrangian approach, starting from basic principles,
leads directly to the closed equation

                          rrt
              
-------
                                                                             15
fusivity of the scalar particles, but in practice the exact form of p is
virtually unattainable because of the extreme difficulty of tracking in-
dividual particles in turbulent fluid, especially the atmosphere.   Conse-
quently, just as the closure problem impairs the utility of Eq.  (14), the
lack of information regarding the exact form of p hampers the use  of its
Lagrangian counterpart, Eq.  (16).  And just as various schemes have been
advanced for rendering Eq.  (14) solvable, several hypotheses have  been pro-
posed for implementing Eq.  (16).   Among the latter,  the most well  known is  the
assumption that p is Gaussian.  Hith this approximation and the  assumptions
that the turbulence is isotropic and stationary, it  can be shown that the
widely used Gaussian plume and puff models are derivable from Eq.  (16).
Similarly, under the assumption that p describes a Markov process,  as would
be true in a study of pure molecular diffusion, it can be shown  that Eq.  (16)
reduces to the diffusion equation, represented by Eqs. (14) and  (15)  (see
Chapter V of this volume).   Unfortunately, none of these models  in  their
current forms provide a wholly adequate description  of atmospheric  diffusion.

     Recognizing the difficulties in obtaining actual  field data to test  and
implement turbulence theories, some investigators pursued the problem of sim-
ulating turbulent flows computationally (Deardorff,  1970, 1972;  Orszag and
Israeli, 1974).  In particular, Deardorff developed  a model that simulates
the turbulent planetary boundary layer, under various stability  conditions,
below an inversion base of constant height (Deardorff, 1970). This model
opens up new avenues, through Eq. (16), to applied studies of atmospheric
diffusion, because in the atmosphere, where molecular diffusion  has a negli-
gible effect on the distribution of the mean concentration ,  the proba-
bility density p that enters in Eq. (16) can be regarded as a function solely
of the turbulence velocity field and can accordingly be calculated  from the
numerically simulated turbulence fields.  In this chapter, we describe some
of our preliminary work along these lines using the  turbulence model  of
Deardorff described above.

     The main objective of our study was to use the  mean concentration pro-
files computed from the numerically derived values of p [referred to as the
numerico-empirical profiles] to assess the adequacy  of conventional atmos-
pheric diffusion theories.   First, we assessed the atmospheric diffusion

-------
                                                                             16
equation,  Eq. (14) with Eq. (15), by determining the profile of vertical
eddy diffusivity that produces the closest fit of the predicted mean con-
centrations with the numerico-empirical profiles.  Second, we compared the
numerico-empirical profiles with those predicted by (1) the Gaussian plume
formula with Pasquill-Gifford dispersion parameters and (2) the Gaussian
puff equation with McElroy-Pooler travel-time-dependent dispersion parame-
ters.  Additional studies of the plume and puff models that we carried out
are relevant to the work described in some .of the later chapters of this
volume.

B.   DESCRIPTION OF THE NUMERICAL EXPERIMENTS

     We considered the three-dimensional dispersion of marked particles in
an atmospheric flow between the ground (z = 0) and an elevated boundary of
                i\
constant height h as simulated by the planetary boundary layer model  of
Deardorff  (1970).  (In this chapter, we use the caret to distinguish dimen-
sional variables from those that have been nondimensionalized.)  In each of
                                      ^,
two experiments representing neutral (h/L = 0) and slightly unstable
(h/L = -4.5) atmospheric conditions, where L is the Monin-Obukhov length,
800 particles were released from points on a horizontal plane and followed
for a given period of time.  Using the computed trajectories of each particles
                                                     s\ /\   /*.  s\
we calculated the two-dimensional density function p(x,z,t|x',z',t')  in the
following way.  (The components of r and r1 were taken to be x, z.and x1, z1,
respectively.)

     First, we transformed the coordinate system of each particle to a moving
frame:
                             •(T) = XI(T) - [x! + U(ZS)T]    ,
                             .(T) = Z.(T) - Zs    ,                      (17)
                                i = 1, 2, ..., 800

where u(z ) is the mean value of the x-component of the wind at the level of
       x s
release of the particles and T = t - f.  Then we introduced a new function,
p    such that

-------
                                                                             17
                                  800h[2X(T)}|"1   ,     |
                                                                         (18)
                         0,    otherwise    ,

where N(C,C,T) is the number of particles at time T that lie in the range
C - 1/2 AC < CI(T)
-------
                                                                            18
  1.0
  0.8
  0.6
  0.4
  0.2
 2s—
-1.0  -0.5    0   _0.5


                (a) Unstable Conditions
                                 -1.0   -0.5     0    0.5    1.0    1.5
zf
  .0.3
   0.2
   O.I
        T=0.l
                                    = 0.4
 •1.0  -0.5    0    0.5   1.0     -1.0  -0.5
                    (b) Neutral  Conditions
                                                   0.5    1.0
  FIGURE 4.  THE  PROBABILITY  DENSITY FUNCTION p(£,z,T/Zs)-
   IN RELATIVE  UNITS—DERIVED FROM THE MODEL OF DEARDORFF

-------
                                                                             19
     We used the results of the above analyses for both the neutral and
slightly unstable cases in Eq. (16) to compute the mean cross-wind inte-
grated concentration, , of a chemically inert substance arising
from a continuously emitting point source.  For S(x,z,t) = M5(x - x )6(z - z ),
where M is the emission rate (a constant) and = 0, Eq. (16)
becomes
                        = lim M    p(x,z,t|x ,z ,t') dt1    .     (20)
Section C presents the results of the calculations.

C.   ASSESSMENT OF THE ATMOSPHERIC DIFFUSION EQUATION

     In the introduction, we mentioned that the implications of the gradient
transport hypothesis [Eq. (15)] regarding the nature of turbulent diffusion
are generally incorrect.  However, the so-called diffusion equation to which
this hypothesis gives rise through Eq. (14) is an attractive equation from
the standpoint of describing atmospheric diffusion.  Specifically, this equa
tion is simple and easily solvable using conventional finite differencing
techniques, and it lumps the effects of the turbulence into the single func-
tion K. . rather than into one or several  additional differential  equations.
       ' J
This last point is particularly significant in the context of photochemical
pollution modeling, where the number of differential equations involved is
already so large that the addition of more equations to handle the diffusion
terms would make the computational burden prohibitive.  Thus, a tacit con-
straint on our analysis was that the closure approximation we derived must
require no additional equations.

     These considerations therefore prompted the following question:   Does a
function K. .  exist that renders the differences between the solution  of the
diffusion equation and the solution of Eq. (16) acceptably small, and if so,
how does this "optimal" diffusivity compare with the expressions  in current
use?  As an added constraint, we insisted that the function K. . depend only
on spatial  coordinates and not on travel  time.  Without this restriction,

-------
                                                                             20
which contradicts observations, K. .  would become a function of the source
                                 ' J
distribution S and would therefore  acquire extremely unwieldly forms in
problems such as urban pollution modeling, where one must deal with a mul-
titude of widely scattered sources  of various shapes and sizes and tempor-
ally variable strengths.  In short,  we attempted to determine whether
turbulent diffusion can be described with acceptable accuracy in terms of
some local, hypothetical property of the flow field.

     To explore this question, we considered the two-dimensional steady-
state form of the atmospheric diffusion equation, in which turbulent
dispersion in the direction of the  mean wind is neglected, i. e.,
            u(z) ^ = 4    (z)       + M5(x - x )  6(2 - z )     .       (21)
                  8X    3z
Upon defining the dimension less variables,
                 -               K7               u(z )h
                                                               ,          (22)
      /
where h is a vertical length scale and u  is the friction velocity,  we obtain
the dimensionless form of Eq. (22):
  30 _ 8  L  3C
U TTT - — |K.7 T7
            Z. o Z

                                        6(x - xJ  6(z - zj     ,          (23)
                                 u    I         s         s
where 6(«) = 6(-)h and where x, z, and t are given by Eq.  (19b).   The bound-
ary conditions for Eq. (23) are

                          C(z,0) = 0                                     (24)

                          Kz|f = 0    ,    z = 0, 1     .                 (25)

     In contrast to the neutral case, in which the depth of the planetary
boundary layer is some fraction of u^/f, where f is the Coriolis  parameter,
Deardorff found that the unstable boundary layer extends up to the height z..

-------
                                                                             21
of the inversion base below which convective mixing is confined.   Conse-
quently, the proper choice of the lengtl
cases and h = u./f in the neutral case.
quently, the proper choice of the length scale h is h = z.  in unstable
     Before proceeding with the calculation of the optimal  diffusivity KZ,
we first consider the mean velocity profiles u that enter into Eq.  (23) and
also the diffusivity profiles K7 in current use.

1.    Hind Profiles from the Planetary Boundary Layer Model

     Figure 5 shows u(z) computed by the planetary boundary layer model  in
the two cases of h/L = 0 and h/L = -4.5.  It is of interest to compare these
profiles with those predicted by conventional  theories.

     According to the Monin-Obukhov similarity theory, the  mean velocity
gradient in the surface layer is given by (Monin  and Yaglom, 1971)
where
                  'z
               Hm \L
                                                                         (26)
                                    ni\ i /
                          9Z   KZ

                          1     ,     neutral  z < KU  /f    ;

                          /     -,-0.25
                          (l  -3 ft         ,     unstable -2  <  z/L  <  0     .(27)
The commonly accepted value of 3 is 15.   The mean surface layer velocity
profiles obtained by integrating Eq. (26) from z  to z + ZQ are given in
Table 1, where the profiles are expressed in both dimensional  and dimen-
sionless form.
     Above the surface layer, there is a change of wind direction with height
(the Ekman layer).  In the planetary boundary layer, v - 1/5 u at the geo-
strophic level, so that |v|  = u, even with turning.  In the unstable case,
i.e., z./L '< -1.5, both observational data and the numerical calculation
show that the change of wind direction is strongly suppressed, resulting in

-------
                                                                   22
28       30
34       36       38       40
                            IL
 FIGURE 5.   WIND  PROFILES GIVEN BY SIMILARITY  THEORY  (S)
  AND DEARDORFF'S MODEL (D) FOR NEUTRAL (DASHED CURVES)
         AND  UNSTABLE fh/L = -4.5) CONDITIONS

-------
                                            Table  1

     WIND AND VERTICAL EDDY DIFFUSIVITY PROFILES IN THE SURFACE  LAYER  FROM SIMILARITY THEORY
  Neutral:   z < 
-------
                                                                             24
a nearly unidirectional  flow in the planetary boundary layer.  Thus,  in the
diffusion calculations presented later,  we neglected the turning of the wind
with height in all situations.   Although several  approximate expressions
have been proposed to describe the wind  speed profile above the surface
layer, none of these agrees as well with the numerically calculated wind as
the profile obtained by simply extrapolating the  surface layer formulation
[Eqs. (26) and (27)] to the top of the boundary layer (see Figure 5).

2.   Theoretical  Forms for Eddy Diffusivities

     The surface layer vertical eddy diffusivities corresponding to the mean
wind profiles derived from similarity theory are  given in Table 1.   Consid-
erably less is known about the behavior  of K7 above the surface layer than
within it.  In two recent studies, expressions for K7 in the planetary bound-
ary layer above the surface layer have been proposed.
                                   ^,
     Ragland (1973) suggested that K-, be taken as constant above the  surface
layer at its value at the top of the surface layer.  In contrast, Shir and
Shieh (1973) assumed that in the planetary boundary layer under neutral con-
ditions K7 obeys the form

                          KZ - U#A(Z),    neutral    ,                   (28)

where £(i) = KZ exp (-4z/z-)> and that under nonneutral  conditions, K7 obeys
the form
                                                                         (29)
where z, is taken to be 10 m.,  In these computations, K7(z,)  is to be calcu-
                      /      \
lated assuming that z= 10 mjlies within the surface layer.   In the dimen-
sionless form adopted here, Eqs.  (28) and (29) become

                             Kz = Kze"4z                                 (30)

and

-------
                                                                             25
              Kz(z) =
                                 4(z-,-z)
   e
                              l
unstable    .           (31)
We postpone further discussion of these proposed diffusivity profiles until
after we have derived our "optimal" estimate of K7.

3.   Estimation of the Optimal K?j_z) in the Atmospheric Diffusion Equation

     In addressing the question of the existence of an optimal  diffusivity
«7, one faces the general problem of estimating the form of a functional
parameter appearing in a partial differential equation [namely, Eq. (23)],
such that the solution of that equation matches certain given data as
closely as possible.  If the data are available as continuous functions of
x and z, denoted here by C (x,z), then the customary criterion  to be mini-
                          JC
mi zed by the choice of K7(z) is

                         &.
                      = f °f
                        J  J
[C£(x,z) - C(x,z)]2 dz dx    ,               (32)
where £fi is the extent of x over which data are available.

     This problem is known as an inverse problem.  As Chen and Seinfeld (1972)
and Chen et al.  (1974) have shown, problems of this type can be solved effi-
ciently by techniques of optimal control theory.   The optimal  control  problem
to be solved is the following:  Determine the function K7(z) that minimizes J,
subject to Eqs.  (23) through (25) and to K7 > 0.   For this problem, the neces-
sary conditions for optimality assume the form of a two-point boundary value
problem:

-------
        C(0,z) = 0    ,                                                  (34)

         KZ|l= °    '    z = 0, 1    ,                                 (35)
                   IF
                      [Kzfl(u)l + 2^x>z) - C(X'Z)1    .             06)

       *U0,z) - 0    ,                                                  (37)

     K7 T7  (?) = 0    »    z = 0, 1     ,                                 (38)
      Z 9z
                     Jo
                       1C JL /jk
                       9Z 9Z \u
           ^L. = _ i   £ii.£_ iy>.\ HV = n                                  CM)
           5«      I   27 27 (,, I °X   U    '                             \^'
                   u
where 6J/6K7 is the functional  derivative of J with respect to K7 and ^(x,
           /-                                                    L
is the adjoint variable to C (x,z)  as shown,  for example, by Lions (1971).

     The two-point boundary value problem given by Eqs. (33) through (39)
cannot be solved analytically and must therefore be solved iteratively.
One straightforward means of determining the optimal value of K7 numeri-
cally is to use the method of steepest descent.  From the definition of
the functional derivative, we can rewrite Eq.  (39) as
                             1   «-
which is a direct expression for the effect of a perturbation in K7 on the
value of J.  The basis of the method of steepest descent is to choose 6K7
such that 6J is negative.  This can be accomplished by setting
                        6K7(z) =W(z)
                          L
I
                                                                         (41)

-------
                                                                             27
where W(z) is an arbitrary positive function of z.

     The algorithm proceeds as follows:

     (1)  Make an initial  guess of K,(z), calling it K?(z).   Select W(z)
                                    L-                 L                \
     (2)  Integrate Eqs.  (33) through (35)  from x = 0 to x = £fi using K7
     (3)  Integrate Eqs.  (36) through (38)  from x = H~ to x = 0 using 1C
          and C(x,z) from Step 2.

     (4)  Compute <5KZ from Eq. (41), and call  it <5K°.
     (5)  Revise the initial  guess of 1C by

                           Kz = K° + 6KZ

          If KZ < 0, set KZ - 0.

     (6)  Using 1C, return to Step 2 and repeat.

     (7)  Continue until (J  - J    )/J  
-------
                                                                             28
computing the new estimate.   The convergent value of Y\7 generally varies
depending on the initial  guess.  Thus, it is desirable to try several
different initial guesses to determine the variability of the convergent
1C profiles.  We carried  out this procedure in the present study, and we
substantially attained the profiles to be discussed subsequently, regard-
less of the initial  guess.  Nevertheless, results on the uniqueness of
profiles determined by using optimal  control theory are still unavailable.

     The resulting optimal diffusivity profile K7(z) for the neutral case
is shown in Figure 6, which also shows the profiles given by similarity
theory (listed in Table 1); the profile suggested by Shir and Shieh, Eq.
(28), and the eddy viscosity calculated by Deardorff.   The optimal  profile
agrees relatively well with that proposed by Shir and  Shieh and shows also
fairly close correspondence to the eddy viscosity profile computed  by
Deardorff.  However, the  diffusivity estimates given by similarity  theory
are consistently too large and differ from the optimal  value by nearly a
factor of 1C at the top of the boundary layer.

     The opposite situation is found in the unstable case (Figure 7), where
the optimal diffusivities are much larger than those given by similarity
theory.  When the latter  values were used in Eq.  (23),  the resulting errors
were a factor of 10 larger than those produced by the  optimal K-, profile.
For all common values of  h, the Shir-Shieh profile for K7 under unstable
conditions  (not shown in  Figure 7) gives smaller diffusivities than does
similarity theory.

     The fact that the optimal K7 profile does not drop to zero along the
broken line shown in Figure 7 cannot be attributed to  the algorithm that
was used to compute K7, because we found that the smaller diffusivity values
near the surface resulted in larger errors in the ground-level concentra-
tions predicted by Eq. (23).  A more likely cause is the larger truncation
errors in the finite difference model of Eq. (23).

-------
                                                                          29
        0,5
        0,4
        0,3

        0,2
        0,1
           0
0,01
0,03
0,04    0,05
FIGURE  6.   COMPARISON OF THE OPTIMAL Kz PROFILE  (SOLID CURVE) WITH THE
             PROFILE PROPOSED BY SHIR AND SHIEH  (DASH-DOT LINE), THE
            PROFILE GIVEN BY SIMILARITY THEORY  (DOTTED LINE), AND THE
              EDDY VISCOSITY CALCULATED BY DEARDORFF  (DASHED LINE)
                   IN THE NEUTRAL (h/L = 0)  BOUNDARY  LAYER

-------
                                                                           30
f.
u
    0,2-
                                                                  0.6
 FIGURE 7    COMPARISON  OF THE  OPTIMAL  KZ  PROFILE  (SOLID CURVE) WITH THE
                   PROFILE BASED  ON  SIMILARITY THEORY PROPOSED BY
             RAGLAND (DASHED)  IN  THE UNSTABLE  (fi/L  = -4.5) BOUNDARY LAYER
              (See the text for an  explanation of  the short  dashed  line.)

-------
                                                                             31
     Having described the differences between the optimal  K, profile and
some of the diffusivity formulas in current use,  we turn now to the im-
portant question of how well  the solutions of Eq. (23),  using the optimal
profiles, compare with the corresponding numerico-empirical  solutions of
Eq. (16).  The latter are displayed in Figures 8(a) and  9(a) for the
neutral and unstable cases, respectively.   Note that in  the  neutral case
the source height z  = 0.09, and in  the unstable case zg =  0.025.  To  fa-
cilitate comparison of these  solutions with those of Eq. (23), we plotted
the latter in Figures 8(b) and  9(b) in the form  of  a fractional error

                                     C(x,z) - C(x,z)
                            s(x.z)  =
where C denotes the solution of Eq.  (23)  and C~ represents  the  correspond-
ing value given by Eq.  (16).  These  figures give the values of  e  multiplied
by 100 so that they can be interpreted as plots of the percentage error.

     The figures reveal that the errors in Eq.  (23)  are nearly  randomly  dis-
tributed in space, except for the neutral case, in which rather large  errors
occur near the source.   The larger errors in this region are consistent  with
the conclusion drawn from theoretical  considerations that the gradient transport
hypothesis  should not hold when the  length scale of  the mean concentration dis-
tribution is comparable to or smaller  than the  Lagrangian length  scale of the
turbulence  (Lamb, 1973).   For the most part, however,  this  hypothesis  appears
to be reasonable.  Errors in the calculated concentrations  are  no larger than
about 20 percent at points farther than about 6u /f from the source in the
neutral case and at nearly all points  at  ground level  in the unstable  case.
Thus, if we find in future studies that the optimal  KZ profile  is insensitive
to the source height z   and to the number, shapes, sizes, and distribution of
                      o
particle sources, and that in such applications the errors  in the diffusion
equation do not increase significantly beyond those  shown in Figures 8(b)  and
9(b), then  we will have succeeded, through the concept of optimal diffusivity,
in transforming the diffusion equation into a model  of atmospheric diffusion
that is accurate enough to be useful  in a wide range of applied problems.

-------
                                                                              32
   0.4-
                                      X
       (a)  Nondimensionalized Cross-Wind Integrated Concentrations
           Computed from Eq.  (16)  Using the Probability Density p
           Calculated from Deardorff's Data for Neutral  (h/L  =  0)
           Stability
       (b) Concentrations Given by the Diffusion Equation [Eq.  (23)]
           Using the Optimum K^ Profile Shown in Figure 6 and the Wind
           Profile (D) in Figure 5.  Concentrations are expressed as
           the percentage deviation (i.e., lOOe) from those given in
           Part (a) above.
FIGURE 8.    COMPARISON  OF OPTIMAL  PROFILE  AND NUMERICO-EMPIRICAL  SOLUTIONS
                                FOR THE  NEUTRAL  CASE

-------
                                                                           33
  0.8-
A-
      0
16
                                    X/h
      (a)  Nondimensionalized Cross-Wind Integrated  Concentrations
          Computed from Eq.  (16)  Using the Probability  Density  p
          Calculated from Deardorff's  Data for the  Unstable
          (h/L - -4.5)  Case
       (b) Concentrations  Given by the Diffusion  Equation  [Eq.  (23)]
           Using the Optimum KZ Profile Shown  in  Figure  6  and  the  Wind
           Wind Profile (D)  in Figure  5.   Concentrations are expressed
           as the percentage deviation (i.e.,  lOOe)  from those given
           in Part (a)  above.
FIGURE 9.    COMPARISON OF OPTIMAL  PROFILE  AND  NUMERICO-EMPIRICAL  SOLUTIONS
                                 FOR THE  UNSTABLE  CASE

-------
                                                                             34
D.   ANALYSIS OF THE DISPERSION PARAMETERS a  AND a
     IN THE PLUME AND PUFF MODELS           Z      }
     Because we use the plume and puff models later, we attempt here to
optimize the performance of these models.  As previously mentioned, in both
the plume and puff models, we assume that the probability density function
p that enters in Eq. (16) has a Gaussian form that is completely determined
by the mean (x,z) and mean square (x2",!2") displacements of fluid particles
from their point of release.  Recall that these quantities are properties of
the turbulence.  Because of the great difficulty of tracking particles in
the atmosphere, past empirical studies of atmospheric diffusion have not at-
temped to measure the particle displacement statistics directly.  Rather, x^
and z2 have been treated merely as free parameters to be used in fitting a
Gaussian profile to concentration measurements made downwind from point
sources of known strength under various atmospheric conditions.  In this
manner, the so-called dispersion parameters of Pasquill and Gifford and
McElroy and Pooler were obtained.  If the probability density p were actually
Gaussian and if the assumptions involved in the empirical determination of
the dispersion parameters were correct, then these parameters should be
                                                      9"      9"
equivalent to the mean square particle displacements x^ and ^L .  Let us
                        —•
-------
                                                                             35
     The computed values of (•? - 2.^}^ as a function of x for both the
neutral and unstable cases are represented in Figure 10 by the open
triangles and open circles, respectively.  In this figure, the nondi-
mensionalized Pasquill-Gifford profiles of a  as a function of x for
three different atmospheric stabilities are indicated by the curves.
Since the positions, but not the slopes, of these curves depend on the
value of the length scale h that is used to nondimensionalize a  and x,
we have indicated the extent to which each curve would be displaced by
a different choice of h.

     Figure 10 clearly shows poor agreement between the calculated and
measured a's.  The chief discrepancy is a systematic shift of the calcu-
lated values to the right of the measured ones.   Thus, if we plot
(2? - z"2)  as a function of 0.24x, we obtain the closed circles and tri-
angles shown in Figure 10, and these agree much  better with the measure-
ments.  Note that in both the neutral and unstable cases, the slope of the
computed a  profile conforms for small  distances x to that of the Pasquill-
Gifford profile for Class C stability (i.e., slightly unstable); but that
at greater distances, the slopes become dissimilar, with the unstable case
becoming parallel to the Class B profile and the neutral case becoming par-
allel to the Class D curve.

     The closer agreement between the computed and measured data that results
when a smaller value of x, namely 0.24x, is used suggests that the Pasquill-
Gifford data pertain to particles released from a lower level than that used
in the numerical experiments.   This conclusion foll.ows from the fact that for
a release height z1, dx/dt - u(z'); and, thus, with a strong wind shear near
the ground (see Figure 5), x will increase with  z1 for a given value of the
travel time t.   This effect has apparently not been taken into account in rou-
tine diffusion estimates made with the plume formula (see Turner 1969).  Con-
sequently, one would expect those calculations to overestimate systematically
the ground-level concentrations resulting from elevated sources, if (as
appears likely) a  is less sensitive to z1 than  x is. .

-------
                                                                                36
 0,1
0,01
                                           fi=5000.'
                  fi=IOO
                                                          u*=0.l
                                                              o°=5000
                                             o   XA"
                                               x
                                            xx A-—Class B  )
                                         xox A    —Class C  \  Pasquill-Gifford
                                      x'         ---Class D  J
                                        A   o  Deardorff  h/L = -4.5
                                            A  Deardorff  h/l_=
                                            A  Deardorff
                                                        /-\
    o Deardorff  h/L =-4.5

J	:	I	I	i   i  i i  i I
    O.I
                    10
  FIGURE 10.    COMPARISON OF THE COMPUTED NONDIMENSIONALIZED ROOT-MEAN-SQUARE
                 VERTICAL PARTICLE DISPLACEMENT a,  (CIRCLES  AND TRIANGLES)
                 WITH THE PASQUILL-GIFFORD a?  DISPERSION  PARAMETER (CURVES)
                   FOR NEUTRAL  AND UNSTABLE CONDITIONS  AS A  FUNCTION  OF
                            DOWNWIND DISTANCE  FROM  THE  SOURCE

-------
                                                                             37
     The Pasquill-Gifford dispersion parameters are useful only in the
plume formula.  The puff model, however, requires o's that are functions
of travel time.  Figure 11 presents data of this type, which have been
collected by McElroy and Pooler in St. Louis.  In this figure, the data
ajre_ shown in nondimensional form, along with the computed values of
(z2 - z2)  for the same stability cases.  Again, the positions, but not
the slopes, of the empirical curves depend on the values of the length
scale h and the time scale h/u^ that are used to nondimensionalize the
data.  Figure 11 shows that for a value of u  =0.98 m/sec, which is not
unreasonable for the urban area in which the McElroy-Pooler data were col-
lected, the computed and measured a (t) profiles are in excellent agree-
ment in the unstable case and in good agreement for small t, under neutral
conditions.  (The downward turning of the computed a  profile at large t
is a result of the effect of the simulated inversion base at z = 1 in the
numerical model.)

     Figure 12 compares the computed streetwise dispersion a  = (x2 - x2)^
                                  	      .                  J\
                                  	
with the lateral dispersion a  = (y2 - y2)^ observed by McElroy and Pooler
                             J
as a function of travel time.  (We did not compute a  originally because it
is not needed in our problem.)

     The McElroy-Pooler data for a  show little variation between neutral
and unstable conditions, and the computed ov values exhibit the same behav-
                                           /s.
ior.  In fact, Figure 12 shows that when the empirical  data are nondimen-
sional ized using a friction velocity of about 0.45 m/sec, about one-half the
value used earlier with a , the computed and observed dispersion curves are
nearly superposed.

     Recall that if the turbulence were such that the probability density p
that enters in the generalized diffusion equation [Eq.  (16)]--and in the
puff equation—were actually Gaussian with joint moments xy = xz" - yT = 0,
as is assumed in the puff model, then the empirical a's would be equivalent
                                                   —7T J^  	— J^
to the root-mean-square second moments of p, i.e., x^ 2> z2  % and so forth,
as described earlier, and the concentrations predicted by the puff model
would agree exactly with those given by Eq. (16).  It follows, therefore,

-------
                                                                           38
    O.I
   0.01
                    i    i
              T—.—T  i  i i . |
                  /\
                                               u^0.8//
                                                £  '
                                                     .1^=0,98, h=IOOO

                                                    > = 500
                                                        •McElroy-PooIer
i  7 i    i   I  i  t i  i I
           I8
-------
                                                                              39
y,
   O.I
  0.01
                                 i  i .1  i
          O
          A
                   o
                   A
                                                        =0.55
                                                               ugO.5,' ft= 1000:
                                                               2000
                                                           = 0.35
                                                           ?0.44,h=u*/f
- McElroy- Pooler cry  h = 1000 meter

         neutral and unstable

 o  Deardorff ' o-x  fi/L=-4.5 (longitudinal o-x)
    Deardorff  crx
i  i  i i
                                                      = 0 (longitudinal crx )
                             i  i i i
         0.01
     O.I
      t
     FIGURE  12.    COMPARISON  OF THE COMPUTED NONDIMENSIONALIZED STREAMWISE
                   ROOT-MEAN-SQUARE PARTICLE DISPLACEMENT ax (CIRCLES AND
                   TRIANGLES)  WITH  THE MCELROY-POOLER LATERAL DISPERSION
                   PARAMETER  ay FOR NEUTRAL AND UNSTABLE CONDITIONS AS A
                                FUNCTION OF TRAVEL TIME t

-------
                                                                             40
that since the empirical a's are in close agreement with the numerically
calculated values of x2'5, z2'1, and so forth, major errors in the puff
model predictions reflect departures of the turbulence characteristics
from those consistent with the assumed Gaussian form of p.  Similar ar-
guments apply to the plume formula.  Let us examine these errors.

     For this purpose, we set up the standard plume model using Pasquill-
Gifford dispersion data, and the puff model of Lamb and Neiburger (1971)
using the McEl roy-Pooler data, to simulate the same two-dimensional  problem
that was treated earlier in Section C.  In each case, we nondimensionalized
the empirical data in the manner portrayed in Figures 10, 11, and 12 so
that each set was in closest agreement with the computed particle statis-
         ~Z
tics x , z , and so forth.  We should point out that the plume model  does
not account rigorously for the effects of the inversion layer present in
this problem.  However, in the trial calculations presented here, this de-
ficiency is not serious, because within the range of downwind distances
treated, the inversion layer has little overall  influence on the concen-
tration distributions.

     Figure 13 presents the results of the calculations for the neutral
case, and Figure 14 gives those for the unstable case.   As in the previous
figures, we have plotted the model predictions in terms of their fractional
departure e, Eq. (42), from the corresponding numeri co-empirical solutions
of Eq. (16).  The latter are given in Figures 8(a) and 9(a).  We consider
only the puff model  results here, because these results provide a better
assessment of the accuracy of the assumed Gaussian form of p due to the
closer agreement between the McEl roy-Pooler data and the computed particle
statistics.
     As shown in Figure 11, for the neutral  case the empirical  a  and the
           —^— i                                                   L-
calculated z2"5 profiles are nearly superposed for travel  times  t < 0.1,
which corresponds to distances x from the source in the range 0 < x < 4.

-------
II
 u*
       (a)  Nondimensionalized  Cross-Wind  Integrated  Concentrations
           Predicted by the  Puff Model, Expressed  as  the  Percentage
           Deviation (i.e.,  lOOe)  from Those  Shown in  Figure  8(a)
           Which Were Obtained from Eq. (16)  for the  Case of  Neutral
           Stability
      0.4-  -
        (b)  Results  of the  Plume  Model  Calculations  Expressed  as  in
            Part (a)
FIGURE 13.    COMPARISON  OF THE  PREDICTIONS  MADE  BY  THE  PUFF  MODEL  AND
             THE  PLUME MODEL  FOR THE  CASE OF  NEUTRAL  STABILITY.  The"
             puff model  used  neutral,  travel-time-dependent  McElroy-
             Pooler data;  the plume model used Pasquill-Gifford  Class  C
             data and an expanded x-axis  (see the  text).

-------
      0
16
      (a) Nondimensionalized Cross-Hind Integrated Concentrations
          Predicted by the Puff Model,  Expressed as the  Percentage
          Deviation (i.e., 100 )  from Those Shown in Figure  9(a),
          Which Were Obtained from Eq.  (16) for the Unstable Case
       (b)  Results of the Plume Model  Calculations Expressed as  in
           Part (a)
FIGURE 14.    COMPARISON OF THE PREDICTIONS MADE BY THE PUFF MODEL AND
             THE PLUME MODEL FOR THE UNSTABLE CASE. ,  The puff model
             used "slightly unstable" McElroy-Pooler  data;  the plume
             model  used Pasquill-Gifford Class B data and an expanded
             x-axis (see the text).

-------
                                                                             43
However, Figure 13(a) reveals that in this range the predictions of the
puff model  are considerably in error.  Concentrations are underestimated
on the centerline of the plume and are greatly overestimated at the
edges, especially along the plume edge that touches the ground.  This
behavior suggests that the diffusing particles behave as if a weak re-
storing force were preventing them from wandering too far away from their
level of release, i.e., the plume axis.  Thus, particles that have wan-
dered a distance z^2, say, from the plume axis have a larger probability
of moving back toward the plume centerline than of continuing on to still
more remote points.   In contrast, the Gaussian density implies that at
each instant the particles have equal probabilities of moving toward or
away from the plume axis.  In short, the Gaussian density has a much
larger kurtosis or flatness factor than the true probability density p
does under neutral conditions.

     Turning next to the unstable case, we see in Figure 11 that the measured
a  values and those computed for z" are in excellent agreement at all travel
times (and distances from the source), except those at the extreme upper limit
of the range treated.  Figure 14(a) reveals that the errors in the puff model
predictions are smaller in this instance than those in the neutral case but
that there is still a tendency for the predicted concentrations along the
plume centerline to be too small and for those at the edges to be too large.
One factor contributing to these errors is the assumption in the puff model
that the joint moment xz~ is zero.   Actually, this moment is strongly positive;
because of the rapid increase in wind speed above the level of release (namely,
            s*
Zpp.  = 0.025h) of the particles in the unstable case (see Figure 5), particles
that are displaced upward are also systematically displaced downstream.  Con-
sequently, at any distance x from the source, a significant fraction of the
particles that normally would be found at the plume edge are found farther
downstream, and in their places are particles from points upstream of x, which
thus have smaller vertical  displacements.  As a result, the effective width of
                                                                         ~2)^
the plume i.s consistently smaller than one would predict on the basis of z  2
alone.   This phenomenon is  not as  pronounced in the neutral case, because the
wind shear at the level of release (znri  = O.lu*/f) of the particles is much
smaller (see Figure 5).

-------
                                                                             44
     The analyses just completed illustrate the inadequacy of the basic
assumption underlying the plume and puff models that the probability
density p entering in Eq. (16)  is of a Gaussian form.   However, despite
these limitations, the plume and puff models remain attractive formula-
tions in applied studies  because of their mathematical  simplicity.   For
this reason, it is useful to inquire whether   profiles exist that  will
bring the errors in the predictions of these two models within the  range
of acceptability.

     To explore this question,  we used the variances a (T) and a (T) as
                                                      A.         L-
adjustable parameters for obtaining a least-squares best fit of the
Gaussian density function to the numerico-empirical expression p.  Figures
15 and 16 present the resulting "optimal" profiles of a  and a , respec-
                                                       X      £-
tively.  In each case, the corresponding empirical data of McElroy  and
Pooler, presented earlier in Figures 11  and 12, are shown for comparison.
Since the method used to obtain the empirical  data is similar to that used
here to.obtain the optimal a  and a  values, it is surprising that  the lat-
                            A      L.
ter  do not compare as well  with the empirical  data as do the computed root-
mean-square particle displacements shown in Figures 11  and 12.

     In a similar manner, we adjusted the dispersion parameter a (x) to ob-
tain the least-squares best fit of the plume formula to the numerico-empirical
concentration distributions given in Figures 8(a) and 9(a).  Figure 17 pre-
sents the results, along with the empirical Pasquill-Gifford data for compar-
ison.  We subsequently inserted these "optimal" dispersion coefficients into
the puff and plume models and performed new sets of calculations to compare
with the numerico-empirical  concentration profiles given in Figures 8(a) and
9(a).  The results, given in terms of the fractional error e used earlier,
are presented in Figure 18 for the neutral case and in Figure 19 for the un-
stable case.  Comparing Figures 13(a) and 14(a) with Figures 18(a)  and 19(a),
we find that although the optimal a  and a  profiles improved, the  overall
                                   A      i-
performance of  the puff model, errors in the predicted values remain intol-
erably large in certain regions.  This is especially true for the ground-level
concentrations predicted in the neutral  case [see Figure 18(a)].

     Figures 13(b) and 18(b) reveal that using the "optimal" az(x)  profiles,
the plume model fails by a considerable margin to achieve the level of accur-
acy that was attained in the earlier calculations.  The reader will recall
that in the previous computations, we expanded the x-axis to account for

-------
                                 o   Optimal  h/L = -4.5
                                     Optimal  h/L = 0
                                     I8
                                      i= ±O.OlJ
  0.05-
FIGURE  15.  OPTIMAL az(t) PROFILES FOR USE  IN THE PUFF MODEL
   COMPARED WITH THE EMPIRICAL DATA OF MCELROY AND POOLER

-------
                                                                46
0.5
O.I
                                                u*=0.5,h=IOOO

                                               u*=0.44th=u*/f
                                    o   Optimal  h/L=-4.5
                                    A   Optimal  h/L=0
                                        McElroy-Pooler
                                       (Neutral  S Unstable)
         0.5
0.
0.5
 FIGURE  16.  OPTIMAL STREAMWISE  DISPERSION ax(t) COMPARED WITH
  THE EMPIRICAL LATERAL DISPERSION DATA OF MCELROY AND POOLER

-------
                                                                         47
    0.5
cr2  O.I
   0.05
                                                        h=IOOO
                                                       = 0.5
                                        o  Optima!  h/L - -4.5
                                        A  Optimal  h/L=0
                                       ---- GLASS  B
                                       - CLASS
                                       -- CLASS
oj
Gifford
                                  _J	I
                                                     10
   FIGURE 17.  OPTIMAL az(x)  PROFILES  FOR  USE  IN  THE PLUME
          FORMULA COMPARED WITH THE  EMPIRICAL  DATA
                   OF PASQUILL AND GIFFORD

-------
                                                                        4'B
    0.4
    0.3

   0.2
    O.I
                 f- -0.10
                 -51/?
                 e2  = 0.77
                               8     10     12      14
                                                        16
                (a) Results of the  Puff  Model
         i=0.02
il
 u*
FIGURE 18.
                (b) Results of the Plume Model


      THE  FRACTIONAL ERROR (lOOe) RESULTING  IN THE  PUFF MODEL AND
      THE  PLUME MODEL USING THE OPTIMAL a PROFILES  SHOWN IN  FIGURES
      15 THROUGH 17 FOR THE NEUTRAL CASE.  Errors are  relative to
      numeri co-empirical  concentration distribution given  in
      Figure  8(a).

-------
                                                                             49
     A

     ~
     h
        0.8-
        0.2-
                                   8     10     12
                                    X
                      (a) Results of the Puff Model
                                                     14
                     (b) Results of the Plume Model
                                              16
FIGURE 19,
THE FRACTIONAL ERROR (lOOe) RESULTING IN THE PUFF MODEL AND
THE PLUME MODEL USING THE OPTIMAL a PROFILES SHOWN  IN  FIGURES
15 THROUGH 17 FOR THE UNSTABLE CASE.  Errors are relative to
the numerico-empirical  concentration distribution given in
Figure 9(a).

-------
 apparent differences in the mean particle velocities implicit in the empirical
data and in the numerical  turbulence model.   We also made arbitrary adjust-
ments of the atmospheric stability.   In contrast, we did not employ any of
these artifices in computing the optimal  a (x) values.   These facts, together
with the error values given in Figure 18(b)  suggest that the basic  plume
formula does not provide a wholly adequate description  of atmospheric
diffusion,  at least under neutral  stability conditions.   However,  under
unstable conditions, the accuracy is better, as indicated by Figure 19(b).

E.   IMPLEMENTATION OF  THE OPTIMAL  DIFFUSIVITY AND
     WIND SHEAR PROFILES IN AN AIRSHED MODEL

     Since  the optimal  diffusivity  profiles  derived in  Section  C are based
on a rather small ensemble of particle trajectories, the statistical sig-
nificance of these profiles is limited.  In addition, the question  arises
of whether  the profiles vary significantly with changes  in the  source
height.  (In our studies,  data from only  one release height were available
for each stability case.)   The data needed to resolve these problems are
currently being collected  using Deardorff's  newest model, and we plan to
continue studies similar to those described above but with a much wider
scope.

     For the present, however, we are developing the methods necessary to
implement these diffusivity and wind shear profiles in  operational  airshed
models.  For this purpose, we have  fitted fourth-order  polynomials  to the
wind shear  and optimal  diffusivity  profiles  derived earlier. Table 2 pre-
sents these in nondimensional  form.   The  only task necessary to use these
results in  a diffusion  model is to  convert the profiles  into dimensional
forms.  Thus, denoting  dimensional  variables by the caret ("),  we  have [cf
Eq. (22)]
                   K (--ju z,     ,     unstable case (z./L = -4.5)    ;  _  (43)
                    Z \ i- */  f I                          1
         Kz(z) = <
                                     neutral  case (z./L = 0)

-------
                                                                             51
Here KZ(Z)  is the nondimensional  diffusi vity, given in Table 2, at height
z (nondimensional) and u ,  z., f, and L are defined as before.   Similarly,
                        "ft   \
we have
                                                                         (44)
where u  is the friction velocity, u is the wind speed profile given
       "K        ^
in Table 2, and h = z. in the unstable case and u /f under neutral  condi-
tions.  Note that the wind profile u given in the table is for a particu-
                 A ^         A.
lar value of R = h/z~, where z~ is the surface roughness.   For different
values of R, the profiles must be adjusted as described below.
     In air pollution modeling studies, one normally has available only
                                                                      *
z., f, z~, and u,n, where the last term denotes the wind speed at 10 m .
Atmospheric stability data are seldom available from measurements.  Rather,
Pasquill's technique of categorizing the stability into Classes A through F
using ground level wind speed and insolation observations is frequently em —
ployed.  Recently, Golder (1972)  obtained approximate relationships between
Pasquill 's stability classes and the Monin-Obukhov length L.  Thus, we
assume that estimates of z./L can be obtained using Golder's formulas and
Pasquill's stability classes.  The friction velocity u^ must be determined
before the wind and diffusivity profiles can be converted into dimensional
forms .
     An often used formula for u  is
                                     = CDug
where Cn is a drag coefficient for the boundary layer and u  is the geostrophic
wind speed.   Lettau (1959)  proposed the formula
                                       Ii u \     ]-l
                                  Io9in f    ~1-8       '    neutral  case (46)
                                     IUVZ0/     J
 A technique for estimating the surface roughness z~ has been proposed by
 Lettau (1970).

-------
                                                     Table 2
                         OPTIMAL PARAMETERS FOR USE IN THE DIFFUSION  EQUATION AND  IN THE
                           GAUSSIAN PUFF AND PLUME MODELS  (NONDIMENSIONAL  VARIABLESt)
    Model
  Gaussian
  puff
Parameter
                                                                  Stability
                 Neutral   -  = 0
                                                        1.5  x  10
                0.0045 + 0.4r,  0.022 < t <_ 0.079
                0.026 + 0.129-t,  0.079 < T < 0.32
                                           Unstable  j1 =  -4.5   —- =  6.8  x  106
                                           .	I;	zo   •	
                                           =  -0.0023  + 0.65t -  1.89r2 + 5.32T3,
                                                               0.04 < T < 0.52
                  «X(T)
                1.85T,  0 < T <_ 0.045
                0.045 + 0.92-r,  0.045 < T < 0.32
                                           -  0.045 +  1.307-r,   0.05  < T <  0.7
  Gaussian
  plume
  oz(x)
0.032 + 0.0044X, 0.23 < x < 4.0
= 0.075 + 0.0167x,  0.5 < x < 5.5
- -0.065 + 0.0435X,  5.5 < x < 9.5
  Diffusion
  equation
              «  7.396 x 10"4 + 6.082 x 10"Zz
                 + 2.532z2 - 12.72z3 + 15.T7z4,
                                 0 < z < 0.45;
                                             -6.934-x 10"3 + 0.6113Z + 3.297z2
                                             -6.442z3 + 3.153z4,
                                              0.02 < z < 1.0
                                0,   z > 0.45.
                   u(z)-
                 29.82  +  213.2z  -  1989z'
                                             26.22 + 153.2z - 1428z':
                                    8743z3 - 14670z4, 0.022 < z < 0.21;        + 5541z3 - 7523_z4,_0.025 < z < O.J3
                                 u(0.21),  z  >  0.21.
                                                             u(0.3), z > 0.3.
t All lengths were made nondimensional  by  zi  »  Inversion  height  in  the  unstable  case and  by u^f in the
  neutral case.  Times were made nondimensional  by  Zj/u*  (unstable)  and 1/f (neutral).
5 The polynomial  given  for the neutral case is only qualitatively  accurate.  See Appendix A for a more
  precise formulation.
                                                                                                                     en

-------
                                                                            53
for neutral  conditions.   Our work with Deardorffs boundary layer model  has
led to a confirmation of this expression.   We have also found that in the
unstable case, the drag  coefficient is given approximately by
CD = 0.156
                                  / i \       I"1
                             loglo(z )~ 1>18       ,  unstable case     .   (47)
Because the geostrophic wind speed u  is  often  not known  in  air pollution
studies, Eq. (45) is difficult to use in  diffusion modeling  applications.
Consequently, we attempt below to obtain  u  from the 10-m wind  u,n  and  the
wind speed profiles that we have available.
     Consider first the unstable case z./L = -4.5,  which is  representative
of the mixed layer over urban regions during much  of the daylight period.
For the unstable case, we have from Table 2,

                   ~/L\
                   u zi                                  ?
                   ~^-LL  = u(z) = 26.22 + 153.2z  - 1428z^
                     *
                                   + 5541z3 - 7523z4
                                                                          (48)
                                     0.025 < z/z.  <  0.3
                          = 0(0.3)     ,     z =  z/z.  >  0.3
Note that this profile is valid only for z^ 0.025, R = Z./ZQ = 6.8 x 106,  and
z./L = -4.5.   For different values  of R, the wind speed profile U(Z,ZQ)  is (for
a given ZQ)
          u(z,zn) = u(z)  - p£n(6.8 x 106 zn)  ,     0.025  <  z < 0.3 ,      (49)
               U           K   V           U/       ZQ < .004    ;
where u(z) is given by Eq. (48) and k is the von  Karman constant  (-,35).   (Note
the restriction on the size of Z.)
     Let
                               Z10 =                                       (50)
                                      i

-------
represent the nondimensionalized 10-m anemometer height,
from Eq.  (49)
                         Then  we  have
                               u(z,n,zn)  =
                                           J10
 10'  0'    u
Therefore,
                                              u.
                                      u  =  .,10
                                       *    D(zio'zo)
                                        (51)
Note that the expression in Eq. (49) for u(z,zn)  is valid only for heights
z > 0.025.  If z,~ is smaller than this value, we are in the surface layer,
where similarity theory gives (see Table 1)
         u(z,zn) = |- 2/tan"1 x - tan^xj  + £n f;
              u    K|_ \               u/       yi
ZQ < z < 0.025
where
                               y —
                                                    ,1/4
     - 15(z
                                                                         (52)
                                               1/4
                                                                         (53)
In summary, we obtain the friction velocity u^ in unstable cases from the
formula
                                       "10
                                "* = "^^O^O7
where u,n is the wind speed at the anemometer height and u is given by
Eqs. (48) and (49) for z]Q in the range 0.025 < ZIQ < 0.3 and by Eqs.  (52)
and (53) when zn < z,Q < 0.025.   Note that the surface roughness z^, the
mixing height z., and the Monin-Obukhov length L must be known to obtain u

-------
                                                                           55
     In the neutral  case (z./L  = 0),  the evaluation of u^ is somewhat
simpler because the anemometer  height is always within the surface layer.
Here similarity theory gives (see Table 1)

                       1     /z + zn\        /z.
             u(z,zn) = r  £n - -     ,    h-L = OJ  neutral  stability
                             \   zr>   /        \L     /
                             \    0  /        \      /                     (54)

Thus, if Z-|Q = 10 meters and u,Q is the 10-m wind speed, we can assume
that
                                         >     neutral  stability    .      (55)
It is important to note that Eq. (54) is very accurate up to a height
z * 0.03 u^/f.   Since f « 10~4 sec"1  for latitudes  within the United
States and u  > 0.1 m/sec over urban  regions during most of the daylight
hours, it is easy to see that 10 m is a value that  is well  within the
range of validity of Eq. (54).

     Having established formulas for  the friction velocity u^, we can con-
vert the nondimensional wind speed and diffusivity  profiles of Table 2 into
physical forms  quite easily using the computer.   Our approach is to define
FORTRAN functions USTAR, DKZ, and UBAR, which calculate the friction velocity,
diffusivity, and wind speed, respectively, given the following data:

          uin    =  anemometer wind speed at height z,,.,
          zn     =  local surface roughness,
          z.     =  local inversion height,
          z./L    =  stability parameter (recall  that the Monin-Obukhov
           1        length L is assumed to be known),
          z      =  level at which the diffusivity  or wind speed is required,
          f      =  2fi sin   = Coriolis parameter.

-------
                                                                           56
     In the execution of the pollution model, the function USTAR is called
for each point on the horizontal  plane to compute the local friction veloc-
ity.  Thus, we have

                       u^ = USTAR(U10,ZO,ZIOVL,ZI,Z10)

where the arguments of USTAR are  given the proper local  values.   Once u  is
determined, the functions DKZ and UBAR can be referenced as follows:

                       K (z) = DKZ(Z,ZI,USTAR,F,ZIOVL)
                       u(z)  = UBAR (Z,ZO,ZI,ZIOVL,F,USTAR)

These functions merely use their  arguments to translate  the polynomials
given in Table 2 into dimensional values.  Note that the functions DKZ and
UBAR simply replace the diffusivity and wind speed variables in  the program.

     Listings of the FORTRAN programs of each of these three functions are
provided in Appendix A.  Note in  the listing of UBAR that in the portion of
the surface layer not covered by  the profiles given in Table 2,  we use well-
known analytical expressions for  u(z).  Equation (52) is an example of such
an expression.  Further information in this regard can be obtained from
Ragland (1973).

     At present, the three functions USTAR, DKZ, and UBAR are applicable only
to the stability cases z./L = 0 (neutral) and z./L = -4.5 (slightly unstable).
The data on which the profiles of Table 2 are based are  currently restricted
to these conditions.  We hope to  be able to extend the scope and accuracy of
this work in the future.  Deardorff (1972) has found, for example, that in all
situations where -z./L is greater than about 5, the proper velocity scale is

-------
                                                                           57
rather than u .   If this  is  true,  then  one K  profile should suffice  for
             "rT                              i-
all  unstable conditions  in the range z./L  < -5,  but this  remains  to  be
checked.   Also,  our K  profiles were obtained  for fixed source heights.
In principle, the diffusivity should be independent of source height;
but because of the heuristic nature of  the concept of a turbulent diffu-
sivity, this may not be  the  case.   Further work  is therefore needed  to
ascertain the influence  of source  height on the  optimal  K  profiles.

-------
                                                                         58
           III   DEVELOPMENT OF A CLOSURE  APPROXIMATION
             FOR THE CONCENTRATION  FLUCTUATION  TERMS
                    IN THE GOVERNING EQUATIONS
A.   INTRODUCTION

     Snapshots  of  smoke  plumes in the atmosphere always reveal  an  erratic
distribution of smoke  concentration within the local  vicinity  of the smoke
source.   In  contrast,  time exposures usually reveal  a smooth distribution
that closely resembles the classical Gaussian profile.   If two  such photo-
graphs of the same plume were superposed, the difference in apparent smoke
concentration at any point would represent the instantaneous magnitude at
that point of the  so-called turbulent concentration  fluctuation.   That is,
if A denotes the instantaneous concentration, portrayed by the  snapshot, and
A denotes the mean, revealed by the time exposure, then

                             A1 = A - A                          (56)

is the concentration fluctuation.

     In  situations where nonlinear chemical  reactions occur among  consti-
tuents of the plume or between the plume's constituents and chemical species
entrained from  the ambient atmosphere, the presence  of concentration fluctua-
tions can have  a profound impact on the behavior of  the mean concentrations
of all the reacting species.  This can be demonstrated mathematically by the
following simple example.  Suppose that two  species,  A and B,  undergo the
bimolecular  reaction
                             A + B

-------
                                                                          59
with rate constant k.  The equation governing the instantaneous concentration
of A is, therefore,

                              {£= -kAB    .                       (57)

If the two species are not homogeneously distributed, local fluctuations
A1 and B1 exist.  The mean concentration A" is then governed by the
equation
                                __    _
                            = -kAB - kA'B1    .                    (58)

This equation follows from Eqs.  (57) and (56) and the assumption that
A1 = B1 = 0.  Equation (58) indicates that, depending on their amplitude,
the concentration fluctuations can actually dominate the evolution of the
mean concentration distribution.
     Since the term A'B' is an additional unknown variable, Eq.  (58) cannot
be solved until another expression relating A'B' to A~ and B" or to some other
known quantities is found.   This is known as the closure problem.   The
search for such relationships, or closure approximations as they are generally
called, has been undertaken by several investigators.  Of these, the most
sophisticated work has been done by O'Brien (1968) and Hilst  et al.
(1973).  Both of these studies, which are similar in nature, follow along
the lines of the so-called invariant modeling approach.   We outline this
method below, using the work of Hilst et al. as an example.

     In the context of the problem represented by Eq. (58), the first step
is to formulate the equation governing A'B' and to add it to the system of
equations governing A and B.  We then have [repeating Eq. (58)]

                ^ = -kAB - kAW                                 (59a).


                ^:= -kAB - kFB^                                 (59b)

-------
                  '—- =  -k|_AB'2  +  BA'B1  +  BA'2  +  AA'B'
                                  	   	                   (59c)
                                +  A'B'2  +  B'A'2J
                                              2        2
The last equation  involves  the  mean  squares  A'   and  B1  ,  and  so  their
governing equations  must  also be  considered:

                   '2
                 dA'
                  dt
  = -2k[BA'2 + AA'B'  + A'2B'J    .             (59d)
                           21-  	          . 	1
                 „„ .            o"	       9"
                 ^-=  -2k[AB^ + BA'B1  + A'B'^J    .             (59e)
                                                 2
The presence  of third-order moments,  such  as  A'B1  ,  in  Eqs.  (59c)  through
(59e)  precludes a  solution  of the  system of equations represented  by
Eqs.  (59a)  through (59e).   Note  that  the problem here is  identical  to  the
one we faced  with  Eq.  (58), i.e.,  too many unknowns  for the  available  set
of equations.   It  is  apparent, therefore,  that  merely writing  down  the
governing equation for each new  unknown  that  appears does  not  in  itself
resolve the difficulty,  for this process produces an infinite  hierarchy
of equations.   At  some point, a  hypothetical  expression must be  introduced
to circumvent the  hierarchy and  to produce a  closed, finite  set  of  equations.
It is  at the  point represented by  the system  [Eqs.  (59a)  through  (59e)]  that
Hilst et al.  introduced their closure hypothesis in the form
2          '2                  M
            l2R,  _  AB     ,  ,  A
                 -

-------
where
             M =
                  0,  if (A'2/A2) (B'2/B2) < 1
                  1,   otherwise
     Equations (59f) and (59g) are not based on any physical  concepts or
on experimental observations.  Rather, they are simply mathematical  approxi-
mations that are designed to satisfy the following realizability conditions,
                 = 0
               if A'B1  = -AB
             A2B = A2B    ,     if A'2 = 0    ;
                   0
2B = A2B + A|2B +
             AB = AB + A|B + 2A'B'A
if A'2B'  = 0    ;
                                                                        (60)
We should note that Eqs. (59f) and (59g) are not the only relationships
that satisfy Conditions (60), nor are these conditions the only ones that
the statistics of A and B satisfy.
     Although Hilst et al. have found good agreement between the pre-
dictions of their model [Eq.  (59)] and the corresponding analytic solu-
tions of certain simplified chemical reactors, from the standpoint of
urban pollution modeling, their approach is not altogether acceptable.   The
objection is that in a system involving m nonlinear reactions among n species,
their technique introduced n  + m additional differential equations into the
existing set of n to be solved.  Thus, in a typical case where there are
five pollutants and five reactions (in which concentration fluctuations
are important), the number of differential equations to be solved is 15
instead of 5.  Needless to say, an increase in computational effort of
this magnitude cannot be accommodated in urban diffusion models given pre-
sent computer speed.

-------
                                                                                    62
             Since none of the existing closure approximations  for the turbu-
        lent concentration fluctuation terms  possess  both the degree of accuracy
        and the ease of implementation that are required  in pollution simulation
        studies, we set for ourselves  the task  of developing a  new approximation
        that would more closely satisfy these needs.   Our primary constraint was
        that the new scheme should  involve no additional  differential  equations.
        The existence of such a scheme has been suggested by our work in the pre-
        vious chapter with the analogous term u]c'.   In  Chapter II, we found
        that by a proper choice of  the functional  form of the diffusivity KZ,  the
        approximation
        could be made sufficiently accurate to be useful  in diffusion  modeling
        applications.  Moreover,  this  formulation does  not require that  additional
        equations be solved.   Using this work as  a guide,  we set out to  find  an
        expression of the form
                           A'B1  = f(A,B,known parameters)

        We describe the expression that we have  developed  and our attempts  to
        demonstrate its accuracy in the remainder of this  chapter.

        B.   DERIVATION OF AN APPROXIMATION FOR  CONCENTRATION FLUCTUATION TERMS
             SUCH AS A'B1

             The following simple example illustrates both the reasoning behind
        our approach and the essence of the closure problem itself.

             Suppose that  rectangular pulses or  slugs of the species  A and  B are
        present in a one-dimensional  system as shown in  Figure 20.  The mean con-
        centration A of A  is defined by
/
                                      XO+L
                                     xo

-------
                                                                          £3
            CONCENTRATION
         FIGURE 20.  TWO INITIAL SLUGS OF SPECIES A AND B
                    IN A ONE-DIMENSIONAL VESSEL
and a similar definition  applies  for  B.   These  results  follow  immediately
from Figure 20.   (Although  in  this  example  L  represents  an  averaging
length interval,  it  could just as well  denote time  or ensemble averaging
domains because  the  results  that  we derive  below  are the same  in  all  cases.)
Since A and B are not  uniformly distributed within  L, concentration fluctua-
tions exist.   At  points outside the slug  of A,  the  fluctuation A'  has  a
magnitude
                         A1  =  A  -  A  =  -
                                        VA
                               (62)
and at points  within  the  slug,
                          A'  •
                               (63)
     Similar expressions  pertain  to  the  fluctuations  B1
     Suppose  that  both  A  and  B  undergo  first-order  decay  reactions whose
rates are  given  by
                              dA
                              dt

                              dB
                              dt
                               (64)
-k.BB

-------
                                                                          64
The concentrations within the slugs then become

                            A = A^At    ,
                                                                  (65)


From these expressions and Eq. (61), we find that the mean concentration
A is
                     A =    A e-kAt =   e-kAt    t                (66)
and a similar equation holds for B.  Here, A,, denotes the initial value
of A~.  Upon differentiating Eq. (66), we find that
                                -  kA
                             dt   ~KAA    •                       (67)
Since this expression is not explicitly dependent upon £A or
the temporal behavior of A is independent of the way in which A is dis-
tributed within L.  We can conclude, therefore, that linear reactions
are unaffected by the presence of concentration fluctuations.

     Suppose, however, that A and B react with each other in the manner

                              A + B -* D

with rate constant k.  In the example shown in Figure 20, it is clear
that this reaction can occur only in the domain (of width £Ag) where A
and B are mixed.   Outside this domain, the concentrations of A and B
remain unchanged for all time.  Let Am and Bm denote the concentrations
within the mixed zone £/\B-  Then
                        dA    dB
                        dT ' dT - -kAmBm    '                   (68)

-------
                                                                      65
and from the definition given by Eq. (61), we have
                   A = U  V£A - £AB> + VAB
                                         Bm£AB
                                                                  (69)
Keeping in mind that all terms on the right side of Eq.  (69)  are con-
stants except Am and Bm, we obtain upon differentiating  Eq.  (69)
                                  =  . .     £AB
                                    "
                     dt    L   dt   "mm  L     '                (70)

and
                        = _kA B
                     dt    KAm m
These equations show that in contrast to the linear reactions  considered
earlier, the behavior of A and B is now explicitly dependent upon  SL/\Q.
In other words, because of the nonlinearity of the chemistry,  temporal
changes in A and B are sensitive to the manner in which the  species  of
A and B are distributed in space and time.

      This  fact is  obvious from  Figure  20, which  shows  that  the quantities of
A and B  that  are  consumed by  the reaction A + B  ->-  D are determined by the ex-
tent  to  which  the  slugs of  A  and B  overlap. Moreover,  this  figure indicates
that  the extent of overlap,  i.e. £/\B>  is not  expressible in terms of A and
F alone  because these  variables are independent  of the positions of the
slugs within  the  interval L.

     Herein lies the essence of the closure problem—the  parameters  that
control A" and  F, in this case Am,  Bm, and np^, are not expressible  in
terms  of only  the  variables  A and  B themselves.   The  origin  of this  diffi-
culty  lies in  the  loss  of  information about the  details of  the distribu-
tion of A and  B within  L when A and B are averaged over the  interval  L.

-------
                                                                           66
Consequently, in situations where these details play on explicit role
(namely, where nonlinear processes are involved), the closure problem
arises in equations governing mean values.   This suggests that solving
the closure problem means retrieving the information that is lost in the
averaging process.   This simple notion is the basic principle underlying
our proposed closure approximation.   Specifically, restricting ourselves
to chemical reactions that do not perturb the fluid flow (all air polu-
tion reactions fall into -this group), we plan to show that most of the
information contained in the term A'B'  is also present in its counter-
part A|Bj, which pertains to the case where the species A and B are
inert, and since AjBj is determined solely by those processes that con-
trol the spatial distributions of A and B, then AjBj, and hence A'B',
are expressible in terms of measurable properties of the system.
     To demonstrate this approach, we consider the problem described in
Figure 20.  Under the passive reaction assumption, the extents of the
regions 2,/\, SL$, and &AB are the same, regardless of the chemical  proper-
ties of the species A and B.  Consequently, if inert materials Aj and
BI are released into this system in initial concentration AQ and  BQ,
and if the same mechanisms that controlled the positions and sizes of
the slugs of the reactive species A and B also act upon AI and BI, then
we would observe
                             —     O^A
                             AT - —r»                      (71a)

and the second moments would be

-------
                                                                           67
                             A2    VA                            .__  .
                             A  =  -—    ,                        (71c)
                                   n
                           AIBI  =     L        '                     (71e)

All  parameters  on the right-hand side of Eqs.  (71a)  through (71e)  are
identical  to the corresponding terms  entering  into Eq.  (69),  and  all
quantities on the left side of Eqs.  (71a)  through  (71e)  represent  mea-
surable parameters.   If the latter are  regarded as properties of  the
system and are  specified along with the other  parameters that describe
the  given  problem,  then the five equations comprising the set can  be
solved for the  five unknowns Ag, Bg,  £^, £R, and £^g.  These  results
can  then be substituted into Eq.  (69) to obtain the  variables Am,  Bm,
and  £/\g that are required to achieve  closure of the  reaction  rate  equa-
tion [Eq.  (70)].

     The results  of the solution of Eqs.  (71a)  through  (71e)  are
                             AB
                                 J__
                                        5
                                 rA
                                                                  (72)

                                 *7
                                  BI
                             R   =  —

-------
where
                             AB " A R
                                  AIBI
                                                    (73a)
                                  A
                                                                  (73b)
                    BT
               *"*     J-


                    Fi
                                                                  (73c)
After substituting these results into Eq. (69) and solving for Am and Bm,

we can write Eq. (70) in the closed form
       dA
       dt
krArB
               LAB
       A - A  (1  -
AB
LAB
                                     [A
                                (74)
Since the right side of this equation is equivalent to
                       -kAB = -k[AB + A'B1]
                                                    (75)
our approximation for A'B' is given by
      A'B1 =
             FArB
              AB
     A - AT 1  -
                 AB
                 B 'J
         1 -
             LAB
             LA /J
  -  AB
(76)
We emphasize that this relationship is exact in the case of the problem

of two rectangular slugs of A and B shown in Figure 20.

-------
     Before we extend this technique to more generalized situations, it
is of interest to test the accuracy of the Hilst et al. closure [Eqs.
(59f) through 59g)] in the context of this simple problem.
To simplify the mathematics, let us assume that
                  = a, =
                                                                 „
                                                                   = £,
                                                                           69
 i.e., that the species A and B are premixed.  If A = B, then we find from
 the  definitions of the bar average ( ) and the fluctuation A1, Eqs. (61) and
 (62), respectively, that
                    A
                     ,3 _
+
+
                                                             (77)
 For  this  quantity,  the Hilst et al. approximation  [Eqs.  (59f) through
 (59g)j  gives

A'3
ADH
A3,
L
A3,
L
1 C ^ J
1.5 - L H
o c *• 4.
O - C> i T
2~
L2 J
2"
o ^
p
                                                £ < 1/2
                                                                  (78)
                                                L/2 < a < L
Figure 21  shows the behavior of the fractional  error of this approxi-
mation, i.e.,
                          e -
                                •5     -5
                               i J   n i J
                                    ADH
                           (79)
                                 A
                                  ,3
as a function of the normalized mean square fluctuation
                                      ,
                            A2  "'
                           (80)
The large size of errors  in this  approximation,  even for this elementary
problem,  emphasizes  the extreme difficulty of finding a closure scheme
based on  hypothetical  relationships  among statistical moments that is
accurate  for a wide  range of situations.

-------
                                                                           70
     Consider now the extension of the closure scheme developed above
for the simple problem shown in Figure 20 to the generalized situation
exhibited in Figure 22.  The disjointed volume v^g represents the total
volume  where A and B particles are found together; and the volumes v^
and vn denote those regions where A and B are found in isolation.  If
     D                                               	
the averaging volume V implicit in the bar operator ( ) encompasses all
of Vn, Vg, and v^g, then we have by definition
            I
                               A
             A dv
                           LVA        ~VAB
                         (81)
                     B = ^
                B dv +
             B dv
                                       'AB
                         (82)
                           J
     AB = -J-  I    AB dv
             VAB
                                    (83)
Let the following variables  be defined:
                   i_  f
                   VA J
A - ~-  \   A dv
= 77-     B dv
                                    (84)
               a  =
                   'AB
        J   A
                        'AB
dv
b =
'AB /
B dv
(85)
                                   'AB
In terms  of  these  variables,  A  and  B  become

-------
                                                                         VI
  oc
-o
  LU
    m o:
   ii
   u
FIGURE 21.   FRACTIONAL ERROR e IN THE DONALDSON-HILST APPROXIMATION
[EQ.  (59f)]  OF A'3  FOR THE CASE OF TWO PREMIXED SLUGS OF MATERIAL IN
     A ONL-PIMENSIONAL REACTOR  (plotter! as 9 function of the
            normalized mean square fluctuation
               FIGURE  22.   MIXING  TURBULENT CLOUDS OF
                 REACTING  PARTICLE SPECIES A AND B

-------
                                                                           72
                               [AvA + avAB]    ,                   (86)
A = v|AvA+ avAB
                           = l|^BvB + BvAB]    .                   (87)
Recalling that the concentrations within and the extents of the zones
VA and Vg are the same regardless of the chemical  nature of the reacting
species, we find that the quantities corresponding to Eqs.  (86) and (87)
for inert particle species ftj" and Bj are
                       Ai =    Av  +           >                   (88)
                       BI
where
          aT = —  f   AT dv    ,     BT  = —  f   BT dv    .     (90)
           1   VAB J     l             l    VAB J     l
                    VAB                         VAB
Finally, we define the two parameters y and y  such that
                     f   AB dv = yabvAB    ,                       (91)
                     VAB
                       ATBT  dv = yTaT6TvflR    .                    (92)
                        I  i        111  MD
                   VAB

-------
                                                                          75
efforts to relate X to indices  of both of the chemical  natures of the
reacting species and of the turbulence and molecular diffusivity.

C.   FORMULATION OF THE PARAMETERS ?  AND r
                                    M      D

     We consider three distinct situations, each of which illustrates to
some degree the physical  significance of these parameters and how each
is related to properties  of the flow.

1.   Release of A into a  Uniform Field of B

     Suppose that a cloud of A  particles of initial  concentration AQ and
volume VQ is released into a turbulent fluid  in which particles of species
B are uniformly distributed. We desire to estimate g\.

     In this situation, v^, and hence CA, are initially zero.  Although
particles at the edge of the cloud are near particles of the other species,
there is no region in which one species is immersed in  the other.   However,
with the passage of time, particles of A, excited by molecular-scale agita-
tions, migrate into the ambient fluid and thereby create a region v^g.
The flux of migrating particles at the cloud  edge can be represented in
the form
                                                                  (101)
where < is the molecular diffusivity of A and dA/dn is  the gradient of
concentration of A particles normal  to the cloud surface.   If S denotes
the instantaneous surface area of the cloud,  then at any instant the number
of A particles moving into the ambient fluid  is  given by
                      djl -  f
                      dt ~ Jc
  E •  n ds                             (102)
S ~   ~

-------
                                                                          76
Since the portion v^ of the original  cloud  that  is  occupied  only  by  A
particles is  completely surrounded  by a  region VAB  created by  the migra-
tion of A particles  into the ambient  fluid,  all  particles  that leave v^
move immediately into v^g.   Therefore,
                                   -  N
where NQ = AnVg is  the total  number  of A  particles  released  and  N  is  the
number still  in v^  at any given  time.   From  Eqs.  (103)  and  (102),  we
obtain
                     =      F  • n ds                              (104>
                dt          F      S
     The particle flux F averaged  over  the  cloud  surface  can  be  expressed
in the form
                             F  =  < -    ,                          (105)
where A and £  are concentration  and  length  scales,  respectively.   We  pro-
pose that the  ratio A/?  entering in  Eq.  (105)  can  be  represented  to good
approximation  using the  following expressions  for  A and 5 :
                                   -    ,                          (106)
                                  V0
                                •$
where e is  the  energy dissipation  rate  of  the  turbulence.

-------
                                                                          77
     The quantity on the right side of Eq.  (107) represents the smallest
length scale associated with spatial  variations in a scalar field in
turbulent fluid (see Batchelor, 1967).   It represents the length scale
associated with the concentration gradient that is required to produce a
molecular diffusive flux of particles that just balances the compressional
fluxes arising from the strain rate field of the turbulence.  This length
scale is thus a logical choice for £, at least initially.  With the passage
of time, however, ? increases, but in a manner that is not possible to
describe exactly.  We are proposing that the overall effect of this in-
crease (and the simultaneous decrease in A) on the flux F can be accounted
for by Eqs. (105) through (107); i.e.,
                                       -    .                      (108)
                                       0
Here, VQ is constant because it is  the combined  volume occupied by all
particles that initially were of species  A.   Combining Eqs.  (108), (103),
and (104) , we obtain
                   dcfl
     It is not possible to obtain a precise description of the rate of
change of surface area S of a cloud of marked  particles in turbulent fluid.
A reasonable approximation for the average surface area might be

                              S ^ a2    ,                          (110)

where a is the distance between a pair of  marked particles released a
distance £Q apart where

-------
                                                                          78
and Sg is the initial  surface area  of the  cloud.   Batchelor (1950)  showed
that for particles whose separation I is within  the inertial  subrange of
a homogeneous turbulence,
                         d£2
Combining this result with Eqs.  (110)  and  (111),  we obtain  the estimate
                   S = S0(l  + £2/3 SQ-2/3  t2)     ;                 (112)
and upon substituting this  expression  into  Eq.  (109)  and  solving  the  re
sulting differential  equation,  we  obtain  finally
                                              ,                    (113)

where

          a = 3K1/4 e'/4 ^     ,                                  (114)

          b - .'/" c11/12 *-7/3    .                               (115)

and £Q is the initial  cloud diameter.

     Under conditions  of strong convection  in  the  atmosphere  with  e  on
the order of 10~2 m^/sec3,  Eq.  (113)  predicts  that for  a  cloud  of  material
of Schmidt number on the order  of  1 and  initial  diamter of roughly 1  meter,
CA will  reach its limiting  value of unity within about  20 seconds  after  release.

-------
                                                                          79
2.   Release of Dilute Mixtures of A and B
     into a Homogeneous Turbulent Reactor

     In view of its  definition [Eq.  (96)], c^ is equivalent, in a homo-
geneous system, to the probability that any inert particle of A will  be
found at time t within some distance 3r, say, of a particle of B.  In other
words, implicit in our definition of v/\g is the condition that the smallest
distance between particles  of opposite species is on the average 3r or
smaller.  The magnitude of  3r is in  turn an implicit function of the chemi-
cal reactivity of A and B because reactions are assumed to occur within
and only within
     Let #|AQ represent the event

             E |A« = given particle A~  is  within  a distance
                    8r of at least one B  at time t    .            (116)

If the concentrations  of A and B are sufficiently dilute that the proba-
bility of finding more than one B particle with  a distance 8r of the given
A is negligible,  then  clearly
prob  B|A0]  = ////  P(r»t,r',t|rAO,tAO,rBO,tBO)SB(rBO,tBO)  dr'  dr dtBQ  dr
              0  0    8V
                                                                  (117)
where (r/\o»t^Q)  is  the  release  point of the given inert A particle,  8v is  a
volume of radius  8r centered  at the  location (r,t)  of a particular inert B
particle, and SB(rBQ,tgg)  is  the number density -of  B particles released per
unit time and volume at
     Equation  (117)  gives the  probability that a particular A particle ex-
periences  event  E.   However, CA is  equivalent  to the event

-------
                                                                           80
      E = any randomly chosen particle of A is within a distance
          6r of at least one B at time t.

If N/\ particles of A have been released prior to time t, then there are
N/\ mutually exclusive ways in which E can be realized, and it follows that
                  il     i  I       f     \
                = H~  II  Pr0'")  ^l^n  ^A^rAO'^AO^ ^AD ^rAO    ' ^^^
                   A    0
where S^(r^Q,t^Q) is the number density of A particles released per time and
volume at (rAO^o)-  Upon substituting Eq. (117) into Eq. (118), we obtain
finally (see Lamb, 1974)
                            AT(r',t) BT(r,t) dr1 dr    .           (119)
                  XA "^ ...  L ~      i ~     ~   ~
If the volume 6v is sufficiently small that Aj(r',t) Bj(r,t) is nearly
constant for all r' in 6v and for all r, then Eq. (119) reduces to
                              Aj(r,t) Bj(r,-
                    dr    .             (120)
A
If, moreover, the chemical reactor has volume V and the statistics of Aj
and BI are uniform throughout V--the homogeneity assumption—then Eq. (120)
gives finally

                                           .                       (121)

-------
                                                                          81
In the limit as  A and  B  became  thoroughly  mixed,  r^g  ->• 1  and  ^  ^ BI  5v<
Since by definition 0  <  r,A  <_ 1,  we  see  that  Eq.  (121)  holds only if  the
concentrations are so  dilute that

                             FT  6v  <  1                             (122)
3.   Release of Rectangular Slugs  of Reactants
     into a Hypothetical,  One-Dimensional  Reactor
     Imagine a one-dimensional  system such  as  that  shown  in  Figure  20  in
which the reacting clouds  have  rectangular  shapes.   In  this  hypothetical
system, we see that
                     A0£AB _  £AB            _  £AB
where AQ is the (uniform)  concentration  o^ inert  particles  in  the  cloud.
The following mean values  are thus  defined:
          A-
          A
           I      L

            	   A£n           	   R^j)
            .2   VA          2 _  VB
Solving these six equations for £y\g and a^ and substituting the result into
Eq. (123), we obtain

-------
                                                                           82
and, similarly,
where
                                                                  (125)
                               = —
                             A   -2
                                 AI
                               - -
                             B"-2
                                 BI
                            AB   ATBT
It turns out (as we show later) that when Eqs.  (124) and (125) are used
in the closure approximation [Eq.  (95)], good agreement with observation
results, even in turbulent pipe flow reactors.   This suggests that Eqs. (124)
and (125) may be sufficiently accurate for all  general  purpose applications
of Eq. (125), but this remains to  be definitely established.

D.   FORMULATION OF THE PARAMETER  A

     Being dependent upon the properties of both the turbulence and the
chemical reaction, A is more difficult to formulate than the mixing para-
meters C  ar)d C-   By virtue of its definition, i.e.,

-------
                                                                           83
                    AB dv
           1
'r    /A
                      dv
                  AB
  B dv
AB
                                             Bj  dv
                                                                     (126)
IjBj  dv
                                   AB
X is a measure of the effect of the chemical reaction on the correlation
of the concentrations of A and B within the mixed zone V/\B-  This can be
seen more clearly in Figure 23, which compares profiles of the concentra-
tions of inert species in v/\g with those that might be observed if the
species were reactive.  In interpreting this diagram, one should keep in
mind that v/\g is defined for the inert rather than the reactive particles.
As a consequence, portions of the zone defined as v/\g may contain only one
of the species rather than a mixture of the two.   In fact, if the reaction
is extremely fast, no part of V/\B contains a mixture of the reactants ex-
cept a very small zone where the reactant clouds  merge.  In this case, it
is clear from Eq. (126) that A -> 0.  However, at  the other extreme where
the reaction is exceedingly slow, the concentration profiles within v/^
are nearly identical for both the reactive and inert materials, and it follows
that A -»• 1.
                                                                    - X
      FIGURE 23.  COMPARISON OF CONCENTRATION PROFILES  IN THE  MIXED
                 ZONE VAB FOR INERT AND REACTIVE SPECIES

-------
                                                                            84
     The behavior of the parameter x  for reactions  of intermediate speed
has been studied by Shu (1975).   In his  analyses,  the differential  equa-
tions governing the combined reaction and molecular diffusion of two reacting
species in one dimension were solved  numerically,  and the values of X were
calculated explicitly for a wide range of rate constants, molecular diffu-
sivities, initial conditions, and the like.   The model  equations used were
                        9A
                        at
/A
     -  kAB
(127)
                        9t
                              A
                                  -  knAB
                                     (128)
                                                                  (129)
                       8t
                                                                  (130)
where k is  the reaction  rate constant,  D is  the  molecular  diffusivity,  and
n is the stoichometric ratio.   These  equations were  solved subject  to  the
initial and boundary conditions
           A(x,0)  =  Aj(x,0)  =  <
                                AO     ,     -6  <  x  <_ 0     ;
                                0    ,     otherwise;
                                     (131)
           B(x,0)  =  Bj(x,0)  =  <
                                      ,     0  <  x  <_ 6     ;
                                0     ,     otherwise;
                                     (132)
                            lim  A, B = 0
                                     (133)

-------
                                                                             85
The computed functional  form of x(t) is shown graphically in Figure 24 for
several  values of the dimensionless group

                             tD   k62nA
                         a = ^= ——-    ,                     (134)
                             LR     u

where tD is the characteristic time scale of the diffusion; i.e.,
                                                                  (135)
and tg is the (initial) characteristic time scale of the reaction,
namely,
                          tR = (knAor     .                       (136)
All curves in this figure are for the case of the feed ratio
                                                                  (137)
equal  to unity.

     It can be seen that in all  cases both the initial  value and the final
values of A are  1.   Initially,  before the reaction has  had time to consume
appreciable quantities of the reactants, the concentration profiles are
virtually identical to those in  the inert case.   In the long time limit, the
eventual  total mixing of the species results in  uniform spatial concentra-
tion distributions  that, as is  easily verified,  are characterized by x = 1,
regardless of the value of k.

-------
     It is important to note that for a < 1, the parameter x has approxi-
mately unit value for all time (see Figure 24).  This greatly simplifies
the functional  form of the closure scheme [Eq.  (95)].  In contrast, signi-
ficant temporal variations occur in X as a becomes large compared with
unity.  A fortuitous and fortunate feature of this behavior is that for all
values of a, X reaches its minimum value Xm-jn at approximately t = 0.5tp and
maintains this value for a period of about 100 tp.  Inasmuch as reacting
species that are initially unmixed can be consumed only as fast as the mole-
cular diffusion and turbulence can bring them together, for fast reactions
(a » 1), the characteristic time scale must be tp rather than t^.  In view
of this, it appears from the behavior of x(t) described above (and shown in
Figure 24) that a constant value of X, namely,

                             X - »mjn    ,                        (138)

should yield good results in Eq. (95) for all values of a.  With this pros-
pect in mind, Shu examined the relationship between Xm-jn and a and found
that for 6=1,

                                 1  T T2"

as shown graphically in Figure 25.

     An interesting consequence of Eqs. (138) and (139) is that for large
values of a, we have
Upon substituting this result into Eq. (95) and using the resulting ex
pression for JB in the rate equation for A, we obtain

-------
10
  -1
10
  -2
10
  ,-3
                                                                        a = 1
                                                                           To"
                                                                          100
                                                                        1,000
                                                                       10,000
               10
                 -4
Source:  Shu (1975).
1.0"
10"^        10"'
      t/tD = t/(«tR)
10
  0
10'
             FIGURE 24.   FUNCTIONAL FORM OF X(t)  FOR  VARIOUS VALUES OF THE RATIO a
                      OF THE DIFFUSIVE AND REACTIVE TIME  SCALES tD AND tR
                                                                                                                oo

-------
 10
•icr1
 10
   -2
 10
   -3
 10
   -4
    10
      -1
10'
      Source:   Shu (1375).
10'

 a
 FIGURE  25.   DEPENDENCE  OF  X  .   ON ex  FOR THE CASE OF THE FEED RATIO (3 = 1
                                                                                                     CO
                                                                                                     en

-------
                                                                            89
                 0
     The absence of the rate constant k from the right side of this equa-
 tion emphasizes the fact that fast reactions are diffusion-limited processes.

     Shu (1975) also considered the effects of variations in the feed ratio
 3  and the number and initial separation of the slugs of reactants.  Figure 26
 shows how x(t) is affected by variations in 3 when a is held constant at 700.
 The most striking effect is the more rapid onset of the rise of x toward its
 final value of 1.  This behavior reflects the more rapid exhuastion of the
 species (in these cases A) that is present in a smaller quantity.  Despite
 this behavior, the assumption
                         x = constant = A  .
 should still prove accurate during the major part of the period  in which
 the  reaction is in progress.  Figure 26 also shows that except for the case
 of 3 « 1, the order of magnitude of Am-jn is unaffected by changes in 3.

      Figure 27 shows the influence on x of changes in the number and sep-
 aration of the initial slugs of reactants.  In each of the three separate
 cases shown, the minimum value achieved by X is of the same order of magni-
 tude.  These results, those presented above, and other experiments not des-
 cribed here suggest that Eqs. (138) and (139) are perhaps accurate even in
 general situations.  In the next section, we examine this possibility in
 light of observational data.

E.   TESTING  THE  CLOSURE  SCHEME  REPRESENTED  BY  EQ.  (95)
     USING OBSERVATIONAL  DATA

     Few empirical  studies  of turbulent  reactions  have  been  reported  in"
which all  of  the  information required  to apply  and  subsequently to test
Eq.  (95) has  been measured.   Usually,  the  data  reported  are  not sufficient

-------
                                                                            10'
Source:  Shu (1975)
                FIGURE 26.  IMPACT OF VARIATIONS IN THE FEED RATIO B
                    ON THE FUNCTIONAL FORM OF x(t) FOR a = 700

-------
10° H
iif'h
10
  -2
10
  -3
10
  -4
   10
     -5
10
  -4
10"
   Source:  Shu  (1975).
10
  -2
10
-1

T_
Tn
10
10'
10-
           FIGURE  27.   INFLUENCE OF CHANGES  IN  THE  INITIAL NUMBER AND
            SEPARATION OF REACTANT SLUGS ON  THE FUNCTIONAL FORM OF X

-------
                                                                            92
to permit estimations of the parameters FAB> CA» and <;g that are properties
of the inert species concentration statistics.  However, using the theory
of Toor (1962),  which relates the concentration statistics of inert species
to those of extremely fast reactants in the same turbulence, Shu (1975)
succeeded in collecting experimental data from the literature that are ade-
quate to apply and test Eq. (95).  (See Appendix  E  for  details.)  We  report
some of his findings here.   For more details, we refer the reader to Shu (1975)

     All of the data used were gathered from the multijet reactor designed
by Vassilatos and Toor (1965).  In this system, the reactants, usually an
acid and a base, are fed through alternate jets in a head consisting of
100 small jets at the end of a pipe several centimeters in diameter.  Down-
stream from the jet head, the mean concentrations are constant on all planes
normal to the pipe axis.  Thus, this reactor approximates a one-dimensional
system in which downstream distance corresponds to elapsed time.

     Using Eq. (95) in the governing rate equations, Shu modeled this
system by

      dFA     kArAR__
       — -         fl R  TF    C\ - r \~\  TF    n   r M          ("\ A.~\ }
       P	 TT	 ATBT iT-n ~ U   t>i\)l  irD - U - CrJ J    ,     U^U
dt      rnrn  nlul L'A   v'   -°A'J l'B
where
                             FA =       ,                          (142a)
                                  AI
                             FR =       ,                          (142b)
                              D   n"
                                  BI
and where Aj and Bj are constants.  Since the reactions are of the form
A + nB -»• products, and since the mean speed IT of material down the pipe
is constant, Eq. (141) can be written in the equivalent form

-------
                                                                            93
      dFfl     k*r   nAT
      3T =  --^§~1 [FA - (1  - ?A>]  [FA - ]  + ^     '         043)
where x represents downstream distance and
is the feed ratio of the inlet jets.

     Shu used Eqs. (124) and (125)  to approximate ^ and eg.   Subsequently,
he determined the values of f^, fg, and f^g entering into these expressions
and into Eq.  (143) by using Toor's  theory, mentioned above, together with
relevant measured data on extremely fast reactions in the multijet reactor
(see Appendix E).   We should add that although  Toor's theory  has  been corro-
borated by several experiments, a more thorough and meaningful  test could be
performed if data were available that did not require this  theory.

     To apply the approximation [Eqs. (136) and (139)] for A, Shu evaluated
a using reported values of k, D, and  other terms, and the measured scalar
microscale for 6.   In all cases, Shu  found that a £ 1, and so he used the
constant value x = 1  throughout the model validation calculations.

     Figure 28 presents the first set of comparisons between  the observed
values of F^ and those predicted by the model equation [Eq. (143)] for two
different values of the feed ratio  B.  In both cases, the agreement between
the model predictions (solid line)  and the observations (circles) is very
good.   The degree to which the concentration fluctuations dominate the
rate of change of the mean reactant concentrations is revealed by the dis-
crepancy between the data points and  the dashed lines shown in both panels
of Figure 28.  These lines show the predictions of the rate equations in
which the concentration fluctuation term TrE"1" has been omitted.  The com-
parisons show that the effect of the  fluctuations is a retardation of the

-------
  CO
    t
  o
  CO
  o
X
  O
  OJ
    rr
        o DflTn OF MflO.  RE=2U30
        - D.  Z. MODEL
            k = 0.83 x TO4 I mole"1 s"1
                     o
                         o
         FLUCTUATION
         OMITTED
                              o
                                  I     I
    0.0
1 .0
  2.0
X,CM
  Source:  Shu (1975).
3.0
.o-
0.0
                                                        	j	j	}	r



                                                        o ORTfl OF MRO. RE=2H30;'
                                                        - D.  Z. MODEL
                                                                                          6 = 1.5
U.O
               FIGURE 28.  FIRST COMPARISON OF CLOSURE MODEL PREDICTIONS WITH OBSERVATIONAL DATA

-------
                                                                             95
speed of the reaction.   This result is caused by the segregation of the
reactants initially and the finite time required for the turbulence and
molecular diffusion acting together to bring the two materials into contact.
As shown in Figure 29,  this phenomenon cannot be described in terms of A and
F alone.  In this figure, we have plotted the solutions of the equation

                              rtT      	.
                              31= -ekAB                          (144)

for several values of the constant e.   Note that no value of e exists
that will bring the solution of Eq. (144) into even crude agreement with
the observations.

     Figure 30 shows a comparison of the theory and observations for a situa-
tion involving both a faster reaction and larger values of the feed ratio
than those presented in Figure 28.  As this figure indicates, the agreement
between theory and observation is better than it was in the previous case.
Additional comparisons  can be found in Shu (1975).

     The results reported above are very encouraging and provide the hope
that our closure model  will perform well under a wide range of conditions,
both in the atmosphere and in the laboratory.  Our immediate goal is to
develop expressions for r/\g, c/\> and Cg that permit Eq. (95) to be applied
to reaction processes in the atmosphere.  Toward this end, we are currently
preparing a numerical experiment in which these functions will be estimated
using the turbulence model of Deardorff.  The procedure will be to release
pairs of particles in the numerical fluid in each of many realizations of
the turbulent flow and to calculate the mean square concentration field from
the ensemble of particle-pair trajectories generated thereby.  In this
way, we plan to derive the forms of f^g, c/\> and eg as functions of the
atmospheric stability,  wind speed, release height, surface roughness, and
the like.  Provided that the resulting expressions have universal features,
the universal functional forms of these parameters can be used in conjunction
with those of the diffusivities KZ presented in Chapter II to study reactive
plumes of material in the planetary boundary layer.  Progress on this project
will be reported in Shu (1975) and in future reports on our EPA contract
efforts.

-------
7.0
                                                                  A   OBSERVATIONS
                                                                  k *  8320
                                                                  0 -  1.5
    0
            FIGURE 29.  ATTEMPTS TO SIMULATE THE OBSERVATIONAL DATA USING A RATE EQUATION
                            THAT OMITS CONCENTRATION FLUCTUATION EFFECTS
                                                                                                             CTl

-------
                          VflSSILRTOSClSSU)
                          K  = 0.12U.E 05'
                          BETfl = 1.Z6
   0.0
4.0
                                                         o
          CO
          o
                                                         CO
                                                         o
                                                       X
                                                         CM
                                                         O


                                                         O
                                                         O
                                                           0.0
                        .0
   VRSSILRTOSC19GU)
   K  = 0.12UE 05
   BETfl =3.98
  2.0
X,CM
3.0
4.0
Source:   Shu  (1975),
                FIGURE 30.  SECOND COMPARISON OF CLOSURE MODEL PREDICTIONS WITH  OBSERVATIONAL DATA

-------
                                                                             98
  IV    DEVELOPMENT OF  A SCHEME FOR PARAMETERIZING  THE EFFECT
               OF SUBGRID-SCALE CONCENTRATION  VARIATIONS
                             ON  REACTION RATE


A.   INTRODUCTION

     The previous chapter  dealt with the instantaneous  concentration fluctu-
ations  produced by turbulence and the effect they have  on  the  reaction rate
of nonlinear  chemical reactions.  In this chapter, we treat  the  analogous
problem of variations in the mean concentration that are of  too  small a scale
to be resolvable by  a grid model of an urban atmosphere.   The  similarity of
this problem  to that of the turbulent concentration fluctuation  is so great
that the closure scheme we developed in the last chapter can be  used here as
the basis of  a parameterization scheme for the subgrid-scale effects.

     Before proceeding with the development and application  of the scheme,
we briefly review the nature of the subgrid-scale variation  (SSV) problem
and then develop necessary conditions under which SSV is of  negligible con-
sequence.   The last  two sections of this chapter are devoted to  the deriva-
tion of the SSV parameterization scheme and to the application of it to air
pollution problems of applied interest.

B.   STATEMENT OF THE PROBLEM

     Consider two pollutants, A and B, that undergo the reaction

                             A + B -> D

to form a third product D.  The equations governing the concentrations of
these species in an  urban atmosphere are generally taken  to  be (for modeling
purposes)

-------
                                                                              99
where

          A   = some time averaged concentration of A,
          u-   = the i-th component of the mean wind speed,
          K. .  = the turbulent diffusivity tensor,
           ' J
          S,,   = the distribution and strength of sources of A,
          k   = the rate constant for the above reaction.

An equation similar to Eq.  (145) governs the concentration of B:

          ^D       f\D     r\       ^D
                                                                         (145b)
Because the analytic solution of this system of equations is not available,
Eq. (145) must be solved numerically.  For reasons given in Chapter I, this
necessitates the averaging of Eq.  (145)  over volumes equal  to that of the
grid cells into which space is discretized in the numerical technique.  The
spatial average is denoted by the tilde  (~) and is defined by (using the con
centration of A as an example)
                                   X+AX   y+Ay   Z+AZ

                                     x  L  L
where r = (x,y,z)  and (ZAX,  2Ay,  2Az)  denotes  the dimension of the grid cells.
Averaging Eqs.  (145a) and (145b)  in the manner of Eq. (146) and assuming that
u and K are nearly constant  over  distances comparable to the grid cell dimen-
sions, we obtain
where SA and AB are defined as in Eq.  (146).   A similar equation can be derived
for B.  Equation (147)  and the corresponding equation for B are used as the
basis of urban pollution models; therefore, it is A and B, rather than the
"point"  values A and B  that these models predict.   The differences

-------
                                                                            TOO
                          A"  =  A - A
                                                                        (148)
                          B"  =  B - B

are the  subgrid-scale  variations.  Using the definitions given by Eq. (148)
and assuming  that

                          A  =  A

                          B  =  B

and so forth,  we can express Eq. (147), and the corresponding equation for
B, in the  more useful  forms

                  ||+ LA = SA - kAB - kA^B"     ,                       (149a)

                  ||  + LB = SB - kAB - k^B"     ,                       (149b)

where for  brevity we have introduced the linear operator

                        L  -=  =1 157 - it "u it     •                 "50>
                                    I      I       J

     The quantity kA"B" on the right side of Eqs. (149a) and (149b) represents
the influence  of the SSV on the rate of decay of A and B.  Until  now, all
models of  urban pollution have tacitly ignored SSV effects, but as we show in
the next section, this omission is generally unjustifiable.

C.   CONDITIONS UNDER  WHICH SSV EFFECTS
     ON  REACTION SPEED ARE NEGLIGIBLE

     If

                             A"B" «  AB    ,                         .. (151)

then it  is  clear from  Eq. (149) that SSV has a negligible  influence on A and B
and can  be  ignored.  To derive the condition under which Eq. (151) is satis-
fiedv we begin by deriving the equations governing A" and  B".

-------
                                                                             101
     Subtracting Eq.  (149a)  from Eq.  (145a)  and using  Eq.  (148)  to  define
A", we obtain

           j.nll                -      _             ,	J
           — ^ LA"  = SA' -  k(AB"  +  BA"  +  A"B"  - A"B")     ,                (152a)

and, similarly,

           P.DII                ~      -,             ,	<•
           H + LB"  = Sg -  k(AB"  +  BA"  +  A"B"  - A"B")     ,                (152b)

where SA' = SA - SA and S£ =  SB  -  Sg.   Multiplying  Eq.  (152a)  by  B"  and
Eq. (152b) by A" and  adding,  we obtain
           3 fA"R"^  4.  rfA"nM  -    v   I 9 A   3B   ,  3A   9B
           -(A B )  +  L(A B  )  -  -  K. .  — -^- +  ^- ^^
                                     L-  '   J     J    I _
                                                              SA'
                                -  k(AB"2 +  BA"B" +  A"B"2  -  B"A"B"

                                +  AB"A" + BA"2  + A"2B"  -  A"^")

Upon averaging this  equation  in  the  manner  of Eq.  (146) and invoking assump-
tions such as  A =  A  made  earlier,  we obtain
                                .  r^
                                'ijL8xi
                                         ill    ~AII  ~nii
                Annil  I  rn'iDli  -    V   I on   OL>   '  °"   °t>
                 B   ~t~  I/A 15   -  -  l\ . .
                                1J
                                          J      J    U

                              -  k[AB"2 +  BA"2 +  /TB"  (A +  B) + A"B"2+ A"2B"]
                                                                         (153)

     The first  term on  the  right side of  this equation causes a  decay of  A"B",
even in  the  absence of  chemical  decay, i.e., when  k  =  0.   Thus,  if A and  B
were inert and  if  A"B"  were nearly  uniform in space,  the steady-state value
of A"B"  would be such that
                          L3_A"  3B"  ,  9A"  8B"   -  RHC" j.  /\'ic»
                          fe--37: +  ^r:^7   -B SA +  ASB
                          L.  '   J     J    ' -J

-------
                                                                             102
                                                                  r ^^~ "^
Let the steady-state value prescribed by Eq. (154) be denoted by (A"B")T.
It follows,  therefore, that the conditions that result in the satisfaction
of

                                 (f-B-'Jj « AB                           (155)

are sufficient conditions for the realization of Eq. (151).  Since these
conditions are not difficult to derive, we use them as  the measure that we
seek of the  importance of the SSV.

     We note first that the left side of Eq. (155) can be expressed approxi-
mately in terms of known quantities using Eq. (154) and the equation corre-
                                   f^~~i*>
spending to  Eq. (153)  that governs  A" , namely,
              aA"        ?        aA" aA"
              on   + rA"  = -9Y.          + A"c;" - KAA"R" + RA"  4- A"
               3t     LA      ^ 1J3X. 9X.      A   K(AA B    BA     A
                                                                         (156)
     In Appendix F, we show that the first term on the right side of Eq. (156)
can be approximated by

                                             \  '  9
                                               A      •                   <157'
where X is on the order of the smallest dimension of a grid cell and K^ is
the turbulent diffusivity in the corresponding coordinate direction.  Also,
it is not unreasonable to assume that A"Sj! is on the order of

                                                                         (158)
From this approximation and from Eq.  (157), we conclude from Eq. (156) that
in the steady state
                                        £    ,                          (159)

-------
                                                                             103
ignoring  the  chemical  decay (k = 0)  as before.  Employing assumptions such
as Eqs.  (158)  and (157)  in Eq. (154), we conclude finally that
     The quantities on the right side of Eq. (155) pertain to the concentra-
tions of the reactive species A and B, rather than to their inert counterparts
considered above in deriving Eq. (160).  Since the decay rates of A and B are
on the order of

                              -TV- = -rr '^ -kAB    ,                        (161)

and since the production rates of these quantities are SA and SB, respectively.
we see that in the steady state

                             AB * T^ (a? + ?)'     ,                       (162)
                                  K    B + A
where

                                  a =  ^    .                           (163)

Assuming that
we can reduce Eq.  (162)  to the form
                                    2SflSR
                              AB      A B
                                   k(sA + SB)

Upon combining this result with Eq. (160), we conclude finally that a
sufficient condition for ignoring SSV effects on the chemical decay of A
and B is

-------
                                                                              104
                       2K ?ASB
                                  ~
                                k(S
                                           (164)
When this condition is  not satisfied,  SSV effects may or may not be important,
depending on the particular physical  situation.   In the next section,  we  con-
sider the derivation of a  necessary condition under which subgrid-scale
concentration variations significantly influence A and B.

     To demonstrate the quantitative  significance of Eq.  (164),  let us consider
a simple one-dimensional problem in which a series of sources  of strength  S  and
width W are distributed with separation £ along  the x-axis,  as shown in Figure
31.  We might consider  this problem to represent a cross  section of a  system of
infinitely long, parallel  streets.
                    T
                    S
                    I
                  FIGURE 31
   HH
SYSTEM OF RECTANGULAR SOURCES
Suppose that S represents  the emission rate of species  B and that  aS  is  the
emission rate of species A.   Suppose further that a grid model  simulation  of
this problem is to be constructed using a grid mesh size Ax.  We desire  to
know whether subgrid-scale concentration variations can be ignored in the  grid
model simulation of the concentrations A and B.

Through a straightforward  series of analyses, we find that in this problem
                         S"S"  =
                         VB
and that
                               I
                                       1    -
                                         ~
                                         S  =
                aSw
                 I
                                           (165)
(166)

-------
                                                                              105
     Substituting  these  results  into  Eq.  (164), we  find  that  SSV  can  be
ignored  if
                                   2K?
                         kS  < -, - * - -     .                      (167)
Normally £»w  and  y  might  be  on  the  order  of  one.   Under  these  conditions,
Eq.  (167)  can  be simplified to

                                  <
                             kS  $  ~-    .                                (168a)
                                  A

For the hydrocarbon-ozone  reaction,  k ^  0  .1  ppnf min"  .  Also, hydrocarbon
emissions from motor vehicles are  on  the order  of 10  gin/mile.   Using these
values we find that  for a  roadway  carrying  1000 vehicles  per hour,

                             kS  *  3  x 10"5sec"2    .                      (168b)
                                               2
A typical  value of K  in the  atmosphere  is  10 m /sec.   Consequently, we  find
                    A
from Eq. (168a) that SSV effects are  negligible in  the  present  problem if

                              A ~ 50 m     .                              (169)

     In the problem  shown  in  Figure  31,  A  is  on the order of £.  Thus, if the
source configuration shown here  represented a network of  actual streets,  £
would certainly exceed 50  m,  and SSV  effects  could  not  be ignored.   In cases
where £>Ax, A  is on  the order of Ax.   In this case, SSV effects could be  sup-
pressed by making  the mesh size  Ax suitably small;  but  this is  not  practical
in many cases, and one is  left with  a problem in which  SSV effects  might  have
to be parameterized  to achieve meaningful  estimates of  the spatial  averaged
concentrations A and B.  In Section  E, we  develop a quantitative measure  of
whether parameterization is actually  necessary.

-------
                                                                              106
D.    PARAMETERIZATION  SCHEME  FOR THE SUBGRID-SCALE
     CONCENTRATION  VARIATIONS
     Having  derived  the  conditions  under which  SSV effects are negligible
in the modeling  of nonlinear chemical  reactions and having demonstrated
that these conditions  are  frequently not realized in actual  air pollution
simulation studies,  we turn  now  to  the development of a mathematical  scheme
for representing SSV.  The scheme we developed  here is based on the closure
approximation developed  in Chapter  III.

     Let V denote the  volume of  any grid cell  in the grid model of A and B,
and let v,, denote the  total  volume  within V in  which particles of species A
are found alone.  Similarly, let VB denote the  volume where  B particles are
found alone.   The volumes  vfl and VR are distinguished from a third volume
                          M     D
VAR in which  particles of  A  and  B are mixed sufficiently for chemical  reac-
tions to occur.   Figure  32 illustrates these three volumes.   By definition,
we thus have
A dv +
                                                A dv
                         (170a)
                                            'AB
                                / r          c
                         B = 1  [J B dv +  J
B dv,
                                     (170b)
                                            'AB
         FIGURE  32.   DESCRIPTION OF THE VOLUMES  v.,  VB>  AND  VAB

-------
                                                                            107
It is  clear  from  the  definition  of  v*  and  vg  that  if  A  and  B were  chemically
inert,  the concentrations within  these volumes would  be the same as  in  the
reactive  case.  It  is  only  in  the domain VAB  that  the inert and chemically
reactive  systems  differ.  This fact is one of the  key points upon  which our
parameterization  scheme  is  based.

     We define  the  following four variables:
             A EE 1    /Adv     ,          B  = -    /
                 VA  J                        VB   J
                      VA                             VB
B dv
             a = 1     /A  dv     ,        b  EE 1     /B  dv     .        (171)
                 VAB   J                       VAB   J
                       VAB                            VAB

     Let the subscript I  denote  the situation where  the given  particle ?pecies
is chemically inert.   Then  analogous  to  Eq.  (171), we define
           ^T  = TT    f   AT  dv     »     h  ET7    /"  BT  dv     •         (172)
            1    VAB  J     l                  VAB  y
                      VAB                          VAB

 In terms of these variables,  A,  B,  Aj,  and  B,  become

                                                                         (173a)
                                                                         (173c)

                                       AB]        '     '                "073d)

-------
                                                                             108
     Treating  v^,  Vn,  and VAB as  invariants in the inert and reactive systems
as we have  done  above  is  permissible as long as the fluid motions, which con-
trol  these  volumes,  are not perturbed by the chemical  reactions.  In the
problems of interest to us, this  condition is always satisfied.
     Next,  we define
                         'AB
                                                                         (174a)
          A2  = 1
          AI    V
A2 dv + f  A2 dv
                               v
                                AB
+ V!VAB]
                                                (174b)
 r  =
 rB -
B2 dv +
J B2 dv
VAB
.if
M
yRB2v
D
                    (174c)
(We discuss  the  p's  later.)   Just as  Aj  and  Bj  are determined completely
by the fluid velocity  field  and  source distribution,  so too are r^g, r/\,
and TQ.   Consequently,  these five quantities can be regarded as known
parameters  that  characterize the given problem.
     Rather than parameterize the term A"B" explicitly, we shall instead treat
the term
                              AB = A"B"  + AB
                    (175)
Using the variables introduced above, we can write
AB = 77
                                     AB dv =
                     (176)
                                 'AB
where x is a parameter, that can vary with time.  The desired parameterization
is then obtained by examining a, b, and vnc in terms of the seven known quan-
                                         Ab.
tities given by Eqs.  (173), (174), and by substituting the results into Eq.
(176).  To accomplish this, we must use two additional approximations because

-------
                                                                             109
there are two more variables in Eqs. (173) and (174) than there are equa-
tions.   We therefore write

                             aj = ?AA                                    (177a)

                             bz = ?BB     .                               (177b)

In terms of the new parameters 5. and cR5 AT and BT become
                                 
-------
                                                                              no
where
                         A ~ ~~0       ^D = ^o~        r /,D - —r^r     .       (181)
                         A   A2   ,    B   B2     ,    AB    AB
Using all  of the above expressions, we obtain finally
                                                                         (182)
where

                                 h = -*-     .                             (183)

Except for the parameters h and g, which we  discuss in more detail  later, all
of the terms on the right side of Eq. (182)  are known quantities.   The next
step, therefore, in implementing this parameterization of AB is to  express the
right side of Eq. (182) in terms of the source distribution, fluid  velocity,
and other variables that characterize the given system.

     This task is simplified somewhat in situations where the species A and B
are released from the same sources, as is often the case in actual  air pollu-
tion simulation studies.  Under this condition,  we have


                             'A = VB

and gA and gD reduce to
     H      D
                              11	             11 „ _
                                                                          (184)
 Using these results, we find that
                                         s*      ^     ^  \
                                                                          (185)

-------
                                                                             m
A similar expression can  be found for i\g,,  -  r.g.   Now,  from Eqs.  (177)  and
(178),  we conclude that the term in parentheses  in Eq.  (185)  is  identically
zero.   Thus,
                                   AD

which reduces  with the aid of Eqs.  (180c),  (174a),  (178a),  and (178b)  to  the
final form
                             f*^f     ^  n, ^
                             AB  - hrABAB     .                             (186)

This simple form of the parameterization  is  valid whenever  the reacting species
are emitted from the same source or in other terms,  when the reactants are pre-
mixed.  That most of the air pollution modeling problems of interest fall  within
this class of problems is justification  for  our continuing  the analysis of the
parameterization scheme with Eq.  (186) rather than  the general  expression  Eq.
(182).  Thus,  in the next section,  we  consider h briefly and then derive  ex-
pression relating f«R to measurable parameters.

E.   DERIVATION OF ?AB IN TERMS  OF MEASURABLE PARAMETERS AND
     APPLICATION OF THE PARAMETERIZATION  SCHEME TO  SPECIFIC PROBLEMS

     Before considering i\B, we  should point out that h is  also a function of
the given system and ranges from a value  of  zero for infinitely fast reactions
to a value of unity for inert materials.   Chapter III further analyzes this
parameter.  In the present analysis, we  assume that h is a  known constant.

     By virtue of its definition, fAB  pertains to- the inert species Aj and Bj
and can accordingly be derived from linear theory.   That is,
                                AB   : -
                                     Ari

-------
                                                                             112
where AT  is  the solution of Eq.  (145a)  with k = 0;  i.  e. ,
                           9AT
                             -+ £A  = S     ,                            (187a)
and, similarly,
                           3BT
                               +iBS     .                            (187b)
We assume that A and B are released from the same sources.   Hence,

                                 SA = aSB    ,                            (188)

where a is the ratio of the emission rates of. A and B.   If the initial  con-
centration of AT and BT are both zero, then it  follows  from Eqs.  (187)  and
(188) that
 and hence that
                                  /V>      ~
                                   2       2
                            fAR = 4--=-4    .                          (189)
                                  R      A
                                  BI      AI
     The solution of Eq.  (187a) is of the form

                   Aj(r.t)  = jfj   p(r,t|r',t') S(r',t') df dr'    ,     (190)

 where p is the Green's function of Eq. (187a) and where, for brevity, we
 have dropped the subscript from the source function S«.  All aspects of
 the transport, diffusion, and interactions of the particles with bound-
 aries are embodied in the kernel  p.  If the turbulence is homogeneous,
 then p has the property

                           P(r,t|r',t') = p(r-r',t,t1)     ,              (191)

-------
                                                                              113
which leads  to  the following simplifications of the expressions for AT and
""?

                       .t
                         p(rrr',t,t')  S(r',t') dt' dr' dr,    ,         (192)
                v(r)   u

               ff p(r-r',t,t')  ?(r',t')  dt1 dr'     ,                     (193)

                  t.t
                     )(r^r',t,t') p(r-r",t,t")  S(r',t')  S(r",t")  dt1  dt"  dr1  dr" dr],
                                                                          (194)

 ~        rr /*t  t
 Aj(r,t) =11 I  C  p(r-r',t,t') p(r-r",t,t") 6(r',r",t',t") dt' dt"  dr1  dr"
         JJ JQ JQ
                                                                         (195)

where
                                 >.AX  ,.,Ay_  ~Az
              6(r'fr",t'.t")  =1|  2  12  JT S(r'+c,t') S(r"+s,t") d5
                ~   ~              Ifi        ~~       ~~
                                J-AX_ J-AV_J-AZ_
                                   222
                                                                         (196)

and V = AxAyAz  is  the  volume  of  the grid  cell v(r)  centered at  r.   The  derivations
of Eqs. (193) and  (195) are presented in  Appendix  B.  Since p is  the  known  Green's
function  of  Eq.  (187),  the parameter  r\n  can be evaluated  using Eqs.  (193)  and
(195) once the  source  correlation function  G has  been prescribed.   (Note  that S
is the source function  that enters into the grid  model  of  A and B.)  We consider
below the forms  that G  acquires  in three  problems  of practical  interest.

-------
                                                                            114
1.    Random  Distribution  of Point  Sources  in  a
     Three-Dimensional  Space

     Consider the  situation in  which  "Point"  sources  of (small)  volume v are
distributed  at random in  space, and suppose that each source has a constant
emission  rate m.   Let the number density of sources D(r)  be sufficiently
large that on scales  comparable to that of the grid cells used in a numeri-
cal  calculation,  the  density D  can be treated as a continuous function.   The
problem just posed might  serve  as  a model  of emissions arising from home
heating units, or, through proper choices  of D, from clusters of major point
sources of pollution.

     The  source strength  function  that describes this array of sources is

                                   •
                          S(r,t) = ™ U(r)     ,                           (197a)

where
                !1     ,    if r  lies within any source    ;
                                                                         (197b)
                0    ,    otherwise

Thus,

                          z+Az    y+Ay   x+Ax
               S(r) = 1  f      f     J      S(r') dx' dy' dz'
                         Z-AZ  ^y-Ay   x-Ax

                    = mD(r)    ,                                         098);

where we have taken the averaging Volume V to have the dimensions 2Az •  2Ax • 2Ay
Similarly, we have
                  AZ   AX   Ay  >2
   G(r'>r")=}f    I     I     \U(r' + 0 U(rn + ?) dsy d?x dsz     .   (199)
                •'-A 7  -AX  -AV

-------
                                                                            115
Let
                           W(r,Ar)  = U(r)  U(r + Ar)     .                  (200)

It can be seen from the definition  of U [Eq.  (197b)] that W is a random
variable whose value is either 1  or 0 depending on whether both r and
r + Ar lie within sources.   If we let

                                Ar  = r1 -  r"     ,                        (201)

Eq. (199) then can be written in the form
                  G(r',r")  =
m
v2
                                     X"+AX  y"+Ay  Z"+AZ
J      J      J
•'v/" Aw «'i<" Aw */-*
W(r,Ar) dr
                                    x"-Ax  y -Ay  z"-Az
                                                                         (202a)
or
                                   • 2
                        G(r',r")  = ^W(r",Ar)    ,                      (202b)
                                   v
where W(r",Ar) is the Volume average of W at the point r".

     Provided that the source density D in the vicinity of r" is sufficiently
large, we can exploit the random variable nature of W described earlier to
write

               W(r", Ar) = (1)  P(l;r")  + (0) P(0; r")

                         =P(l;r")    ,                                  (203)

where P(n;r") is the probability that W(r",Ar) has the value n, where n = 0, 1.
If Ar = (Ar  . Ar , Ar )  is larger than  the dimensions of the elemental sources,
        v  x    y    z        3
i.e.,
                               Arx > vx     •
                               Ary > vy
                                                                         (204)

-------
                                                                            116
where v = v v v ,  then P(l;r") is just the product of the probabilities that
           x y i.
the vectors r" and r" + Ar fall within sources.  This is a result of the sta-
tistical independence of the source positions.  Simple reasoning is sufficient
to show that the probability that a point chosen at random will lie within a
source is just the fraction of the total volume occupied by sources in the
vicinity of that point.  Thus, we conclude that

        W(r", Ar)  = v2D(r") D(r" + Ar)    ,    |Ar| > dimensions of v    , (205)

and hence that

   G(r',r") ~ m2D(r') D(r")    ,    if |r'  - r" |  > dimensions of v    .    (206)

     Hhen  |r' - r" |   has a value smaller than  the source size, the determination
of W is slightly more difficult because the probabilities that r' and r" each
fall within sources are no longer independent.  Instead, we have
                                                                           (207)
where P....-(n,m) is the joint probability that U(r") has value n and that
simultaneously U(r" + Ar) has value m.  By definition,
                                                                           (208)
where Py,,(n|m) is the conditional probability that U(r") = n given that
U(r" + Ar) = m, and where Py.O) is the probability that U(r" + Ar) =  1.
From this definition and from Figure 33, we conclude after some reasoning
that
                                       i ,  i ..i-   ,,,, i \ n I r" \
                      PuU-Ollsr") = -v  + (v  "vvv } D(r }                 (209a)

where
                      V = (v  -  |Ar |)(v  -  |Ar  |)(v  -  |Ar  |)     .       (209b)
                             A      A    jr      y    i-      t-

-------
                                                                            117
To see this, note first that it is  given that one point, say, r",  lies  in a
source, for example, the source marked v-,  in Figure 33.   Since the separa-
tion vector Ar is also given, the second point r" + Ar must lie somewhere
in the volume denoted v2 in Figure  33.  If this second point happens to lie
in the shaded volume, i.e., v', then both  U(r") and U(r" + Ar) will  have
unit value; but if the point falls  outside the shaded region, then U(r" + Ar)
will have zero value, unless the point happens to fall within a neighboring
source.  The two events, r" + Ar lies in v'  and r" + Ar lies outside v', are
mutually exclusive.    Thus, Eq.  (209) follows immediately.
             FIGURE 33.   RANGES  OF  THE  VECTORS  r"  AND r"  + Ar
                   GIVEN  Ar AND GIVEN  THAT r"  LIES  IN v
                                                       1
     We found earlier that
                         = Pu(l;r" + Ar) = vD(r" + Ar)
 but since |Ar| < v , v , v  in the present analyses, we have
                                    = vD(r")
(210)
 Combining this result with Eqs. (209), (207), and (208), we obtain the informa-
 tion required to evaluate W, given by Eq. (203), and subsequently G:

-------
                                                                            118
G(r',r")  =
                                           -i
                       '  - (v2 - vv1) D(r")J    ,      |Ar| < dimensions of v   .
                                                                         (211)
where v'  is given by Eq.  (209b).  In summary,
  G(r'r") =  -<
              . * j
             m  D(r') D(r")
                                    if any component of r1 - r" is larger than
                                    the corresponding dimension of v     ;
                       '  -(v2 - vv1) D(r")]     ,    otherwise;
                                                                         (212)
                       where v1 = (YX -
                               v = vxvyvz
                                                 v  -  |Ar  |)(vz -  |Ar
Upon substituting this expression into Eq. (195), we can solve for the mean
                                        ~2
square inert concentration distribution Aj arising from the random array of
point sources.
     For simplicity, let us assume that
                            D(r) = D = constant.
                                                                        (213)
 If the diffusivities K,, and K7 are also uniform in space, then all spatially
                      n      ^,    ~o      ~
 averaged quantities, such as A,, Aj, and S, will be independent of spatial varia-
 bles.  Consequently, regardless of the speed of the mean wind, there can be no net
 advection of spatially averaged quantities, and we can therefore assume that u = 0
 without any further loss of generality.  In this case, we have
    p(r-r',t,t') -
                              exp
                           H z
(x - x1)' + (y - yT   (z - z')'
           2                ?
         2or,"              te
           H                z
                                                                          (214a)
                 = 2KH(t - t1)
                 = 2K(t - t1)
                                                                        (214b)

-------
                                                                           119
The problem thus posed is one of a uniform spatial  density of continuous
point  sources acting in a homogeneous turbulence  of arbitrary mean speed.
     From  Eq.  (195), we now have
   Aj(r,t)  = m
•'/'(Iff
   0         Ar  >v
                              1/3
D p(r-r",t-t")  p(r-r"-Ar,t-t') dAr dr"   (215)
                                                  p(r_r",t-t")
                   Ar  <  v
                         1/3 v
                 p(r-r"-Ar,t-t') dAr dr"\ dt1  dt
 where we have assumed that
                                 =  Vy =  vz = v1/3
                                                              (216)
 The  Gaussian form [Eq.  (214a)]  of the  kernels p permits a reduction of the
 first  integral within the braces  in  Eq.  (215) to
1
1 -
^
v1/3
Ka^2 +

aS2)
1/2
L
1

A(c
v1/3 1
i1 + a" )
                                                                           (217)
 where
                    ^ = 2KH(t  -  t1)     ,    a[)2 = 2KH(t - t")
 and so forth.  The second integral within the braces in Eq. (215) reduces to
                                          -  7vD)	
                                     (a'2  +  a"2)  (a'2 + a"2)1/2
                                                             (21;

-------
                                                                            120
                                                                   1 /3
In deriving  Eqs.  (217)  and  (218), we had to assume that a^, o^, » v   .  Fur-
thermore,  we assumed  earlier  that vD «1.  Consequently, Eq. (215) can be
written  in the  form
fa - *{(
D2 +

(2.)

3/2

2KH[(t -

t
D
')

+ (t -

t")]

3/2

(2K.

,)1/2
                                                                    dt1 dt1
  •2,A2 .     Dm2t1/2
= m D t  +  	o i /o—
                                        -  D
                                                                        (219)
Similarly, we find that
                         72,   ,x    '2n2.2
                          Aj(r,t)  =  m D t
 Making use of this result and Eq.  (219),  we  obtain  from  Eq.  (189)
                                                                        (220)

                                                          ,2/3
                    Dt3/2     1
                                                  t»
                                                      Y 2/3 „ 1/3
                                                      KH    Kz
                                                                        (221)
 As indicated, Eq.  (221)  is not valid in  the  limit  as  t->0.   This  limiting
 value is easily obtained, however,  when  one  notes  that  at  small  times  the
 material is confined to  a puff of volumes  v  equal  to  that  of each  source.
 Straightforward analyses then lead  to the  result
                                                                         (222)
     The significance of these  results  to the  behavior of A and  B  in  the
 problem under study becomes  clear once  Eq.  (221)  is  incorporated into our
 parameterization [Eq. (186)] and the result is used  in Eq.  (147).   We obtain
                        SA        •
                        -JJ-JT+ LA = mD -  k'AB
                                                                        (223)

-------
                                                                             121
where

                       k-  = kh(l  + 	JW-T77)                        (224)

and

                                   /2 -
                             a =
                                      I/2
The parameter k1  is an effective rate "constant" whose departure from the
nominal  value k is a measure of the SSV effects on the chemistry of A and
B.  In the case where h - 1, we note that

                           lim k1  = k,

and that
                               k1  - 	Uj-^    .                         (225b)
The first of these indicates that SSV effects die out with time.   From
Eq. (222), we find that

                           lim k1 -        .                              (226)
 In view of the second property,  Eq.  (225b) ,  the importance of SSV effects
 diminishes as the source number  density D,  the turbulent diffusivity K, or
 both increase.

     Let us examine the SSV effects  quantitatively.   Because of the constancy
 of Sn and SR, A is homogeneous in space.  This fact results in

                                 LA  = 0

 Assume also that SA = SB so that A = B.  Equation (223) then reduces to
                              dt
                                                                         (227)

-------
                                                                             122
the solution of which is
                            A(t) =
            mD
            k1
                                             - 1
                        (228)
where
                                 Y =
It can be seen from Eqs. (228) and (225a) that the steady-state value of
A is independent of the SSV:
                               - limA(t) =
                                                 ( 229)
     Using A^, we can define a reaction time scale TR:
TR = (kAj-1 = (mkD)"'
                                                                         (230)
This time scale is a rough measure of the time required for the chemical
reaction A+B->D to achieve a rate equal to that of the source emission rate.
Substituting TD into Eq. (224), we find that by the time the chemical decay
              K
rate and the emission rate are in balance, the effective rate constant k1
has a value
                                          a(mk)3/4  1
                           k1 = k
           1 +
                                         /4 K K 1/2
                                            KHKz
                        (231)
Since k1 has a maximum at t = 0 and decreases thereafter, it follows from
Eq. (231) that if
              mk
       y =  V4
           D ' K,,
                                         3/4
                                           1/2
• 1
                                                                          (232)
then SSV effects will  play a significant role in the behavior of A and B as
they approach their steady-state values.  This fact is illustrated graphically
in Figure 34, where we have plotted the ratio
                                      A
                                                                          (233)

-------
FIGURE 34.   INFLUENCE OF SUBGRID-SCALE CONCENTRATION  VARIATIONS
               ON THE CHEMICAL RATE  OF CHANGE  OF  .A
                                                                                         CO

-------
                                                                             124
for several  values of the dimensionless group y.  Here, ASSV denotes the mean
value of A,  which accounts explicitly for the SSV using k', and A denotes the
corresponding value obtained using k.  According to Figure 34, y must have a
value comparable ti
are manifest in A.
                                           q
value  comparable  to  or  larger than  about  10  before  significant  SSV effects
     It  is  of interest  now  to  apply  our  criterion,  Eq.  (232),  to  an  assess-
ment of  the SSV  effects  in  actual  problems  encountered  in  pollution  simula-
tion studies.  One  such  problem  is the simulation of large urban  areas
containing  strong point  sources,  such as  power  plants and  refineries.   The
question here is whether such  sources produce a significant  SSV impact  on
the predictions  of  an airshed  model.

     To  approximate the  conditions of a  problem of  this  type,  we  let D
approach the value  of one source  per grid cell.  In  a typical  case where
                              3     3
the grid cell  dimensions are 10   x 10  x  50 m,  we have

                           D  =  2  x  10"8  m"3

We assume that the  source emission rate  in  this case is  100  gm/sec.   This  is
a rough  value for the NO emission  rate of power plants.  We  also  consider  here
three conditions of atmospheric  stability.   These and the  corresponding values
                                         o
assumed  for Ku and  K are as follows (in  m  /sec):
             n      Z
Condition
Stable
Neutral
Unstable
KH
5
50
50
K
ID'1
1
25
Substituting  the  above  values  in  Eq.' (232), we  obtain  the  curves  of y  as  a
function  of k shown  in  Figure  35.
                     3
     Note that y  = 10  is  the  smallest  value  of y  shown  in the  figure.   Conse-
quently,  all  combinations  of rate  constant and  stability classification  shown
in this figure represent cases where  significant SSV effects  occur, at least in

-------
                                                                             125
problems  of the  type  considered  here.   For reference,  we have indicated the
values  of the  rate  constants  for some  of the  important photochemical  pollution
reactions.   We should add,  however,  that for  some  of these  reactions, the value
of y cannot be read directly from the figure  because we assumed that the
concentration  of the  two reactants are approximately equal.  In cases
where the reactant concentrations are in the  ratio a (a
-------
«»•
^s.
CO
  : ^J-
   ^.

   "o

  III

  a
           Note:  .0 = 2 x 10"8 m"3

                      2
                 m = 10  gin/sec
                                                     k—5,/mole/sec
                                              2
                 K,, and K, for source height of 10 m.
                                                                                                            k •>. 10
                                                                                                                 10
                                                                                                        10-
                 FIGURE  35.  VARIATION OF  THE DIMENSIONLESS PARAMETER y AS  A FUNCTION  OF
                            RATE  CONSTANT k FOR VARIOUS  ATMOSPHERIC STABILITIES
                                                                                                                              CT)

-------
                                                                           127
where N  is  the  number of "point" sources in the averaging volume.   Let  A
                                                                      X
and A  denote the separation distances between streets,  and assume  that
                            AX » A  » w    ,
                            Ay » A  » W
                             y     y
Then we have
                       N = Mk\ + /^\H
                            /AXAY\/XX +
                            I  WA Vy
Therefore,
                                   A  + A
                               n =
                               U
                                   WAZA A
                                       x y
     The source  strength function S that describes the network  of  streets  is

                           S(r,t) = ^U(r)     ,                         (236)

where m is  the mass  emission rate of each "point"  source that constitutes  the
network, v  is the  point source volume given by Eq. (234), and

                              !1 ,   if r falls on  a street,
                                                                        (237)
                              0,   otherwise.

The source  correlation function G [see Eq. (196)]  is therefore  given  by
                                           U(r' + c)  U(r"  + c)  d-c   d^  .d^
                                                  "    ~     ~     x  y    z
                                                                        (238)

-------
                                                                               128
   where V = AxAyAz.  Making the change of variables
                              r  = r  + Ar


   and  using the fact that G(r',r'+Ar) is independent of r' for all r' lying in

   the  slab 0< z          X      X
   Therefore,
                               and  |Ar
                         whA_xAy if  |Ar  | ^ nx
                           -v          XX
                               and  |Ar  I = mx
                         whDV  if  |Ar  | = nX
                                      X      X
                               and  |Ar  |
        G(r',r'+Ar) =  G(Ar) =
                                 A ^— + Dd)(Ar )  d>(Ar )
                      WX AZ   WX AZ     vx   X;  VV   y
(240)
                               if  0 < z'< h    ,   and z' + Az < h     ;     (241)
                     6(Ar) = 0   for all z' > h

-------
                                                                            129
                      AX
                                        -w
                                                                         -Av(r)
                       ~]
FIGURE  36.   UNIFORM  NETWORK  OF  STREETS  ILLUSTRATING  THE  CALCULATION
                OF THE  SOURCE CORRELATION  FUNCTION G

-------
                                                                                130
    In  this expression,
                                   0    ,  if |Ar  | *  nx      ;

                                               A      A
                                   1     ,  if |Ar  |  = nx   and  (AT  |-Amx
                                               xx        y     y
                                                                           (242)
    Similar  values apply for 4>(Ar ).   Substituting  these  results into Eq. (195),


    we obtain






t) =  I   ill  p(r-r',t,t') p(r-r'-Ar,t,t")  G(Ar,t't")  dAr dr1  dt1 dt"
                       p(c+Ar,t,t")  G(Ar)  dAr d§  dt'  dt"
     0  0
                                                 (243)
                            P
                                  K>'
                                    + a'2)
                            exp
                                              |Ar I ,mx
                                            ,2
dAry
              + term in  (Ar ) + term in cj>(Arx)(Ar• )} dt1  dt"     ,
    where
**•;?
                                     m2h
                                                                           (244)
    and
                                    = aH(t")
    and so forth.   Evaluating  the  expanded term in Eq. (243), we obtain

-------
     term in
                                   exp
                oi 2 ,   i2\
                2(aH + aH  ,
                                                       f
2 ,   ,2x1/2

H + °H }       m=-
                                       <(Ar )  and  has  the
                y                                             x
 form
     term in
                /
                (a
                                Zexp
                                 m=-o
                    (mx )'
                             w
                          2 .  rt,2x1/2

                          H    H '
                  exp
                         2(a
                                                                \
                                                                       (246)
where
                           A  = A  ~
                            Y    x Ay
                                               (247)
On evaluating the product term in Eq. (243) and combining the result with


Eqs.  (245) and  (246), we obtain finally
   Aj(t) =
                             m2h
            0"0
                 w
       / 2     ,2x1/2
       (aH + aH )
r 2^roexp
                                    -(nx )
                                     v  x
                                                                    2 -,
                           exp
                     dt1 dt1
                       (248)

-------
                                                                              132
     Determining the mean value Aj  is not quite as complex.  We find through
straightforward reasoning based on  the homogeneity of A, that
                      Aj(t)  =
Total Mass in  V
        V
                                   mDVt V
                                 AxAya_
                      Az(t)
jut
wa
X  + A
_x	y_
  A A
   x y
(249)
Using this expression and Eq.  (248),  we can determine the variable i\B [defined
in the present instance by Eq.  (189)  that is required in the subgrid-scale
parameterization scheme [Eq.  (186)J.   Unfortunately, Eq. (248) involves non-
elementary integrals that can  be evaluated only by numerical integration tech-
niques.   Such methods will therefore  be necessary to implement the present scheme
in an airshed model.  For the  present, however, it is useful to obtain c-.n order-
of-magnitude estimate of the  effect of the subgrid-scale concentration variations
on the chemistry.

     According to  our scheme,  the measure of the importance of the subgrid-scale
chemistry effect is  the degree  to which the parameter
                                     1
                               1AB
                                     A2
                                     AI
differs from unity during the period when the generation of material by sources
and the chemical  decay are unbalanced.   Through simple, straightforward analysis,
we find that initially
  r   .1 I
  FAB ^ ln~
                                                                         (250)

-------
                                                                              133
assuming that A  ^ A  -v, x.   Note that rAR  is  independent of the emission
               X    jr                  MD
rate.   Thus,  assuming an effective emission depth h = 2 m, a street width
w = 10 m, and an average street separation A  = 150 m, we find that in a
diffusion model that employs a grid spacing AZ = 20 m in its lowest level,
     We know from the work presented in the previous sections that r,,  decays
because of turbulent diffusion and ultimately attains a constant value of
unity.   The time scale of this decay is given roughly by

                                            X2
where it is assumed that both a  and a,, are proportional  to t   .   Since rflR
                               2_      n                                   MD
is initially much larger than unity, significant subgrid-scale chemistry ef-
fects will  arise if T  turns out to be comparable to or larger than the time
scale Tp, for example, required for the chemical decay and source emissions
rates to reach an equilibrium.   The time scale Tp is equivalent to the time
required for

                                  kA 2 ^ S

where S denotes the source emission rate.   In the present problem,
                                         •
                              S = mD = -|r-
                                       WAAz

and
                                       • 9 9
                                A2 ,  4m V
                                f\     O O  O    •
                                     W JTAz
Thus,
                                    /L.-M/2
                                                                         (252)

-------
                                                                              134
     In  the  present  problem,  we  find  that
                      Tr ^ £--v, 400  sec                                   (253)
                       1   i\i I
               o
(using  K,, ^  60m /sec),  and
                        ,v                     _                           (254)
                       E
Subgrid-scale  effects will  be  most  pronounced for fast reactions  in which  T
                                                             7
is small.   For the  case  of  the NO-0-  reaction, we have k -^ 10 £/mole/sec.
Thus,  if vehicles emit about  1  gtn of  NO  per kilometer, on a street network
           4
in which 10  vehicles pass  each point per day,
                           TE  * 200  sec     .                              (255)

From this  estimate  and  from Eq.  (253), we  conclude that in areas  with  street
network densities  less  than or comparable  to  that treated here (i.e.,  streets
150 m apart),  and with  traffic densities greater than  or equal  to 10   vehicles
per day, subgrid-scale  concentration variations  will  have a nonnegligible ef-
fect on the  simulation  of ozone and  nitrogen  oxide concentrations.  To account
for these  effects,  we must incorporate a scheme  such  as that developed here
into the simulation model.
3.   Isolated  Point  and  Line  Sources
     The  analyses  presented  in  the  last  two subsections  deal  with  homogeneous
source distributions.   The results  of those analyses  are therefore applicable
only in portions of an  urban  region where  sources  are fairly  uniformly distrib-
uted upwind  of a given  site.   In  this section,  we  briefly consider the magni-
tude of subgrid-scale effects  arising from isolated point and line sources,  such
as power  plants and major highways.

     For  the purposes of this  study, we  assume   that, in a point source plume,
the mean  concentration  is described with adequate  accuracy by

-------
                                                                              135
             c(x,y,z) = ^— exp - -^exp  - -^    ,                (256)
                        i . f~r i-c      •  —   / I     l^.£_E                   *    *
where Q is the emission rate,

                                2   ,,  x
                               °z - kz -    '

                               a2 = K  *    ,
                                y    y u
and where the x-axis lies on the plume center! ine with the mean wind parallel
to the x-axis.  Assuming that the grid volume V = AxAyAz has dimensions much
larger than the y and z dimensions of the plume, we find after integrating
Eq. (256) that
where d is the effective source diameter.  This value of fflR pertains to the
grid cell  that contains the source.  At a point nAx downstream from the source,
f^g can be found by replacing £n(Ax/d) by £n[nAx/(n + I)AX].  At each point in
space, the value of fAR is constant, because of the constancy of the source
and meteorology, but the value that emitted material "sees" as it moves down-
stream decays with time.

     In the context of the present problem, the effect of subgrid-scale concen-
tration variations on reaction rate is accounted for in our parameterization
scheme by  replacing the true rate constant k by an effective value k' given by
                               k1 = rABk    .                              (258)

The resulting perturbation in the space averaged concentration due to the subgrid-
scale variations is therefore
                                dc _
                                           ~   rr
                                     1 + (kcjt)  '

-------
                                                                              136
where t = AX/U is the residence  time of material  in  each  grid cell.   With
Eq. (259)  reduces to (assuming that kCjt > 1)
                                       AX/TIT
                                                                         (260)
Typically;  AX = Ay, AZ ^ 20 m,  and K K  ^ 100  n$sec .   Recalling  that  our
analysis assumes that the plume dimensions a  and  a  are  much  smaller  than
AZ and Ay,  respectively, we conclude that for  sufficiently  large  u,
                                    * 20    .                             (261)

This result indicates that the perturbation due to the  subgrid  scale  is  very
large for reactions in which (kCjAx/u)  > 1.  The NO-0-  reaction in  power plant
plumes is one example of a situation that satisfies  this  criterion.   When
kCj(Ax/u) becomes much smaller than unity, the subgrid-scale  effect diminishes.
These results agree with the analyses of randomly distributed point sources
presented in Section IV-E-1.

     Turning finally to the subgrid-scale effects produced  by isolated  line
sources, we assumed a line source plume concentration  distribution  of the form
                           c(x,y,z)  =   ^exp -(-%-)                    (262)
                                                \2oz /

where Q.  is  the emission rate per unit length  of source  and
Integrating  Eq.  (262)  and assuming  that Az»a2,  we obtain

-------
                                                                              137
                                  -1/2,
                           f   -  u   Az
                           FAB -   l/2
As in the case of the point source, the subgrid-scale-induced perturbation
in the concentration will  be on  the order of
                                      -1/2
                                      u   AZ
                           c    1-AB    A 1/2  K 1/2
                                            z

when
                                kCj  i~ > 1     .                           (265)
                                    u
Since the line source produces  a mean concentration level  CT  of

                                       Q,
                                                                         (266)
                                      UAZ

Eq. (265)  is  satisfied if
                                     -2
                               kQL >          .                            (267)

For the NO-CL reaction with k = 10  £/mole/sec,  Eq.  (267)  reduces  to

                            QL > 10"3 gm/sec/m                           (268)
                                        3
when u = 3 m/sec, Az = 20 m, and Ax = 10  m.   If NO  emissions  are  on the order
of 1 gm/km, Eq.  (268)  is  satisfied only by roadways  carrying 10  or more vehi-
cles per day.  Few highways rank in this category.   Even for those that do, r.R
is only on the order of unity; consequently,  little  if any subgrid-scale chemistry
effect occurs.

     We conclude  that  in  contrast to  the point sources, line sources of commonly
encountered strengths  produce negligible subgrid-scale chemistry effects, at
least insofar as  the N0-03 reaction is concerned.  All .other reactions that are
explicitly treated in  air pollution simulations  are  slower than this one and
accordingly are  less affected by subgrid-scale concentration variations.

-------
                  V   DEVELOPMENT OF  A SUBMODEL
             FOR RESTORING POINT SPATIAL  RESOLUTION
                TO GRID  MODELS  OF URBAN POLLUTION
A.    INTRODUCTION

     In  Chapter  I, we have discussed and explained the inability of grid
models  (i.e., those in which space and time are discretized) to resolve
features  in  the  concentration field smaller than the discretization interval.
Generally speaking, urban diffusion models of the grid type employ grid
networks  with a  horizontal mesh size of several kilometers and a vertical
mesh several tens of meters in length.  Since most major sources of air
pollution, such  as power plants, refineries, and highways, have scales much
smaller  than these, there is a great deal of fine structure in the concentra-
tion distribution that grid models cannot resolve.  Indeed, the locations
and intensities  of concentration maxima, the extent and location of the
zone of  oxidant  depression near roadways, and the pollutant exposures in
street canyons are examples of important pollution phenomena about which
urban-scale  grid models can provide no information.  Moreover, validation
studies  of grid  models use monitoring station data that represent time
averages  at  a fixed point rather than averages over the large volumes re-
presented by each grid cell.  Consequently, without an estimate of the
amplitude of small-scale concentration variations in the vicinity of each
monitoring station, the validation process can be hindered, because it would
be difficult to  ascertain whether discrepancies between the computed and
observed  concentrations were due to small-scale spatial variations or to
errors in the model.

     These weaknesses have been cited by some in arguments favoring the
development  of alternative pollution modeling approaches.  Spectral methods
are among the alternatives most often proposed.  In this chapter, we develop

-------
                                                                             139
a "microscale  model"  (not  based  on  spectral  methods)  that can be used in con-
junction  with  an  urban-scale  grid model  to achieve point spatial resolution
of the  concentration  distribution at  any specified point.  In developing this
microscale  model,  one fundamental constraint was  imposed:  The microscale
model must  operate totally independently of the urban-scale model.   That is,
we sought the  capability of examining the near-source small-scale concentra-
tion distribution  at  any point without having to  make multiple runs  of the
large-scale model. We also wanted  the microscale model  not to be structured
around  or coupled  to  the urban-scale  model.   If all  of these conditions are
met, the  microscale model  can be used with any grid  model,  and the  required
calculations can  be performed with  a  minimum of time and effort.

     Before proceeding with the  mathematical  development of the microscale
model,  we review,  for orientation purposes,  some  of  the basic concepts intro-
duced in  Chapter  I.
     Most air  pollution models in current use are based on  the premise that
the processes  governing the concentrations of each pollutant can be  described
by an equation of  the form

                                             S +  R    ,                (269)
                       .             ..
                 at    i  ax.    ax.   ij  ax.
                           I      «        J
where
           u-  =  the  i-th  coordinate  component of the mean wind,

          K. .  =  the  turbulent  diffusivity tensor,
           ' J

           S  =  the  strength of  all  systematic sources  of the pollutant,

           R  =  the  time  rate  of change  of the concentration due to
                chemical  interactions  among all  of the  pollutants present.

In adopting Eq.  (269),  one  tacitly assumes that if all  of the parameters  in
this  equation  were  known  precisely,  then the solution c(r,t) of  the result-
ing equation  (and  initial and  boundary conditions) would correspond exactly
to an error-free measurement of  the  concentration  at the point (r,t).

-------
                                                                          140
{Actually,  c  in  Eq.  (269)  represents  a mean value, which we can think of as
a time average over  an interval  of several  minutes, say, depending on the
nature of the diffusivities  used [see Lamb  (1971)].}

     If Eq.  (269)  is to be solved by  numerical  methods executed on a digital
computer, this equation must first be "filtered" to remove all  small-scale
variations  that  the  grid network cannot resolve.  One way to achieve such
filtering is  to  space average the equation  at each point over a volume V
equal  to that of the grid  mesh.   Thus, if the grid network has  mesh dimensions
of AX by Ay by AZ,  then it can be shown that any space averaged variable de-
fined by

                                x+Ax  y+Ay z+Az
                                                                      (270)
                                X-AX y-Ay Z-AZ

is essentially free of spatial  variations unresolvable by the grid network.
[In Eq.  (270), r = (x,y,z).]   Similar definitions can be written for S and R.
Averaging Eq.  (269) in the manner of Eq.  (270) and assuming that u and K are
nearly constant on scales comparable to the grid dimensions, we obtain

                 9C .  77  oC  _  a  i/   8C  ic in
                       i             u                .
                                         J
     \
Numerical  methods  can  then be applied to this equation to obtain an approx-
imate solution  for c.

     Thus, a  grid  model yields the quantity c(r,t) rather than the point
value c(r,t).   The difference between these quantities,  which we denote by

                      c"(r,t) = c(r,t) - c(r,t)     ,                   (272)

is what  we refer to as the subgrid-scale variations (SSV).  Upon subtracting
Eq.  (271)  from  Eq.  (269),  we  obtain the equation governing the SSV:
                     77      =      L'       + c u i  nil
                     u            k         s  + R
                at      i  ax.    ax.

-------
                                                                               141
where

                      S"(r,t) = S(r,t) - S(r;t)
                                                                      (274)
                      R"(r»t) = R(r,t) - R(r,t)

Since R is not a function of any spatially variable parameter other than
the concentration,  it is clear from Eq.  (273) that subgrid-scale variations
c" can arise only from subgrid-scale variations S" in the source strength.
Given the usual  dimensions of grid meshes used in airshed models (see the
earlier discussion), there are no urban  areas in which S", and hence c", is
everywhere identically zero.

     It is the mean value c(x,t), rather than the space averaged mean c(x,t)
given by the grid model, that possesses  the point spatial resolution and
hence all of the information that we require in air pollution studies.
According to Eq. (272), we can acquire this information from the grid model
by calculating in addition the subgrid-scale field c".  Thus, our approach
to the restoration  of point resolution to the grid model is simply to develop
a microscale, or subgrid-scale, model  based on Eq. (273) whose output c"(x,t)
can be combined with that of the grid model c(x,t) to describe the point
field c(x,t) at any arbitrary point.

     To illustrate  an important property of the subgrid-scale variation c"
and to give the reader a better physical feel for the nature of c", we
examine this field  quantitatively in an  elementary example in the next
section.

B.   QUANTITATIVE ILLUSTRATION OF THE SUBGRID-SCALE VARIATIONS c"(r,t)

     For simplicity, we consider a one-dimensional problem governed by the
equation

                                 2
                         |f = K-^-|+ fi(z) 6(t)                       (2-75)
                         9t       ^

-------
                                                                            142
with initial  and boundary conditions

               c(z,0)  = 0
           lim c(z,t)  = 0
                                             (276a)

                                             (276b)
where K is the turbulent diffusivity, assumed to be a constant, and 6 is the
delta function.   The solution of Eqs. (275)and (276) is
                   c(z,t) = (4TrKt)~1/2 exp f-z
                                             4Kt
                                              (277)
Now let
                                 z+Az
                   c(z,t) =
                            2Az
              c(z',t) dz1
                                Z-AZ
Averaging Eq.  (277) in this manner, we obtain
              c(z,t) =
                       4AZ
erf
                   - erf
where
                                         (278)
(279)
erf(x) = -2=  f
        •v7r  J
                                           -X
                                               dx
                                     0
is the standard error function.   Note that Eq. (279) is the solution of
the equation obtained by averaging Eq. (275); i.e.,
                                                                      (280)

-------
                                                                            143
where

                          I 1    ,    -AZ £ z < AZ
                   u(z) =
                          I 0    ,    otherwise

Let us assume that Eq.  (280) is a grid model  representation of Eq.  (275)
and that c, as given by Eq.  (279), is the model's output.  The subgrid-
scale variations are therefore [from Eqs.. (277) and (279)]
c"(z,t) = -= exp  -                f  L       „ erf
                                4AZ         4K
                                                                      (281)

We can see from this expression that the amplitude of the SSV decreases as
the discretization interval AZ in the grid model is made smaller.  More-
over, for fixed AZ, the SSV in this example decreases with time.  The latter
effect is shown graphically in Figure 37.  The figure also shows c for com-
                                                                            p
parison.   It can be seen that for the earlier travel time (i.e., t = 0.031AZ /K) ,
c" is up to three times as large as c and is therefore not a negligible
variation.

     To obtain a more concise description of the relative magnitude of the
subgrid-scale variation c", we note first that c" is largest at the source
location (i.e., z = 0).  Thus, the maximum amplitude of c" relative to c
can be expressed by


                                •             •
From the  expressions  for c"  and  c,  we obtain

                   p(t)  - j=.  [a*  erf (^)j "'


where
                                                                       (283)
                                  \/4Kt
                              * =     "

-------
                                                                                  144
  1.8




  1.7




  1.6




  1.5




  1.4




  1.3




  1.2




  1.1




  1.0




  0.9




  0.8
c
o


S 0.7
J-
+J


§ 0.6
c
o
LJ

  0.5




  0.4




  0.3




  0.2




  0.1




    0




  -0.1




  -0.2




  -0.3




  -0.4




  -0.5




  -0.6
     A
            ="f
                                C (a* = 0.25}
                           \
\    0.4 \  oTe     oTs      uo

 \         N.



  \            ^\
   \              ^-	

     \                     /

      \
                                                    1.2
                                     Amplitude
FIGURE 37.   ILLUSTRATION OF  THE SUBGRID-SCALE CONCENTRATION  VARIATIONS

              ARISING FROM AN INSTANTANEOUS POINT  SOURCE

-------
                                                                           145
That is,  o* is the approximate half-width of the pollutant cloud at time t
normalized by the grid mesh size Az.  It is evident from Eq. (283) that a*
is the key parameter that determines the relative size of the SSV in this
instance.   Upon evaluating Eq. (283) for several values of a*, we obtain
p = .10 when a* = 0.1, p = 0.34 when a* = 1 , and p = 0.08 when a* = 2.  From
these values, we conclude that, in the case of a point source, the SSV can
be neglected relative to c after the pollutant cloud has grown to a width
of about four grid intervals.

     These results take on added meaning when one realizes that the one-
dimensional equation considered here is an approximate model of a two-
dimensional steady-state plume.  In this case, time is related to the down-
wind distance x from the plume source by the relationship

                               x = ITt    ,                             (285)

where u is the mean wind speed (directed along the x-axis).  Translated
into terms of the two-dimensional plume, the conclusions reached above re-
garding the relative magnitude of the SSV indicate that SSV effects will
be significant within a distance of approximately
                                x =  j                                 (286)

 from the plume source.  In the grid model, this distance x will represent

                                                                  e*   (287)

 grid intervals from the source.  Thus, for fixed u" and K, the  number of
 grid cells affected by SSV actually increases as the grid size Az is in-
 creased, albeit both c" and c become insignificant sufficiently far from
 the source.

     Although we considered only the maximum value of c"  (i.e., c"(0,t) .in
 the analysis above, we can expect  (see Figure 37) that the SSV will become
 negligible everywhere in space on  the same time scale as that  of  the
 maximum.

-------
                                                                              146
C.    DERIVATION OF A SUBMODEL OF c" FOR LINEARLY REACTIVE POLLUTANTS

1.    Discussion of General  Operational  Problems

     As stated in the introduction to this chapter, we seek to obtain c"
from its governing equation [Eq. (273)] and to use this result along with
that given by the grid model  to obtain  the point mean field c(x,t).  Since
our earlier discussion of the equation  governing c", i.e., Eq. (273), was
very general, it did not reveal several of the operational problems that
arise in solving this equation in given situations.  We discuss these below,
within the context of a simple problem, before proceeding with the develop-
ment of specific models.

     To keep the mathematics  as simple  as possible, we consider here only
a single pollutant that decays at a nonlinear rate.  The equation governing
the concentration of this pollutant in  the urban atmosphere is assumed to be

                    H+ LC = -kc2 - k + S    ,                  (288)

where [_ is the operator

                      L s «1  SJq-- ^ KIJ alj    •                   <289>

u.  is the mean wind speed,  K-• is the turbulent diffusivity tensor, and
/  2
(c1 ) is the turbulent concentration fluctuation term discussed in Chapter III,
Although this example is mathematically much simpler than that for photo-
chemical air pollution, which involves  a system of coupled nonlinear equa-
tions, the problem posed by Eq. (288) has many of the essential character-
istics of the more complicated real situation.  Thus, techniques developed
for a pollutant that undergoes a second-order reaction can be applied to
more general situations through straightforward extensions.

     Suppose that a grid model for the  pollutant cited above is to be
developed.  As outlined earlier, the procedure for deriving the working
equation of the model is to perform at  each point a spatial average of

-------
                                                                            147
Eq.  (288)  over a volume V equal to that of the grid cells.  In this way,
we obtain  the working equation

                                   *"»     '—»
                 |f + [_c = -kc2 - kc"2 - k + S    ,              (290)
                 d U

where c and S are spatial averages and c" is the subgrid-scale concentration
variation.  Note that it is the spatial average mean square turbulent fluc-
tuation term <(c  )> that enters into the grid model equation.  Upon subtract-
ing Eq. (290) from Eq. (288), we obtain the equation governing c":
          9t
                        :"2 - 2kcc" + kc"2 - k" + S"    ,        (291)
         2
where " represents the subgrid-scale portion of the mean square fluc-
                 o
tuation field (c1  >.  The second term on the right side of Eq. (291) represents
the chemical interaction of the subgrid-and supergrid-scale concentration
fields.  Since the latter is a known quantity at each point (it is given by
the grid modol),  Eq. (291) contains only one dependent variable, c" (after
t^^>
c"^ and (c'^) have been parameterized), and this equation can thus in
principle be solved.

     The nonlinearity of Eq. (291) dictates that its solution can be approxi-
mated only by using numerical integration techniques.  It is to be expected,
however, that the implementation of these techniques will require the solu-
tion of operational problems unlike those associated with the modeling of c.
This requirement arises because point and line type sources give rise in S"
to delta functions, which are not amenable to treatment using discrete grid
numerical  methods.  We could, of course, attempt to evaluate Eq. (291) on
a fine grid, but this would only create new subgrid-scale problems.  To
circumvent these problems, we develop a new technique later in this chapter.

     In contrast,  in cases where the pollutant is linearly reactive, the
implementation of the microscale model is quite simple because the governing
equation,  Eq. (291), reduces to the much simpler form

                       3r "
                         -+ Lc" = -kc" + S"    ,                     (292)

-------
                                                                           148
and if the pollutant is inert, the equation reduces further to

                           ||^+ Lc" = S"    .                         (293)

     In the remainder of this section, v/e develop the known analytic solu-
tion of Eq. (292) into a working model of c" that can be used with any grid
model of CO, S0?, or other essentially linear pollutants.

2.   Derivation of the Working Equations for a Microscale Model of
     Linearly Reactive Pollutants

     Consider a pollutant that decays linearly at a rate described by

                              ^| = -kc    ,                           (294)

where k is the rate constant.  (Note that inert pollutants are represented
by k = 0.)  The equation governing c" in this case is Eq. (292), the solu-
tion of which is of the form  [see Lamb and Neiburger (1971)]
               t
c"(r,t) =  / /   p(r,t|r't') exp |~-k(t - t1)] S"(r',t') dt1 dr1
            'O


                +  f p(r,t|r't0) exp [-k(t - tQ)] c"(r',t0) dr'

                                                                      (295)

where p denotes the Green's function of Eq. (292) and its boundary conditions,
For  simplicity, we assume that c"(r,tQ) = 0. in which case the last term in
Eq.  (295) vanishes.  Our aim now is to derive formulas for c" pertinent to
general source configurations S" in conditions of open terrain where the
form of p is rather simple.  Consider first the source function S".

-------
                                                                          149
     The source  distribution S  can be expressed in the general form
                  Spa(t)  6(^ -  rJ  + 2  SL3(t)  6(Z - Z3}  6[y - S3(X)
                                     P                                    (296)



where S    denotes  the  emission  rate (mass/time)  of the a-th point source (of
       pa

which N  are present),  r  is  the position of the  a-th source, and 6 is the


delta function.   The line sources represented by the last term in Eq. (296)


have strengths S\_   per unit length and lie along the curves y - Sg(x) at


elevation Zp.   From Eqs.  (296), (274), and (270), we obtain



                N
     S"(r,t)  =     S  (t)  6(r - ra)  +     SLg(t)  6(z ' V
                                       Q. ..
                                        T ~} V
                                        t J IS.  it   /  \                      / OQ "7 \

                                      AxAyAz  ijk ~    '                 ^    '





where Q. .,  is the total  mass of pollutant emitted per unit time in the grid
       1 J K

cell (i,j,k)  and



                            !1,   if r lies  in the grid cell (i,j,k);


                                                                         (298)


                            0,   otherwise.





Upon substituting Eq.  (297)  into Eq.  (295), we obtain an expression for c"  of


the form


                        N               M             L
             c"(r,t)  =  r.  P  (r,t)  + >   L (r,t)  - >   V (r,t)     ,     (299)
                        '     a ~       .«•! *   p ~      i  ••   Y ~

                       a=l             3=1            Y=1


where P  is  the  contribution to c"  from the a-th point source, L. is the con-
       a                                                        p

tribution to c"  from  the  3-th  line  source, and where V  is the collective con-
                                                      Y

tribution of all  sources  in  the y-th grid  cell.   Below we derive the functional


form of each of  these terms, assuming that in the  open terrain and quasi-steady


flow regimes of  interest  to  us here  the kernel p in Eq.  (295) has the Gaussian


form

-------
                                                                           150
p(r,t|r',t')  =
where

          exp -
(x

- X1 - UT)
9 2
2o
X
2 (y

- y' -
9 2
2ay
v,)2]

exp
(2 - Z')2
9 2
z
+ exp
(z + z1)2]
9 2
2a
                   T = t - f,
                  °  = F(T)»
                                                                         (300)
                                                  (301a)
                                                  (301b)
           u,  v      = constants,
           w         = 0.
                                                  (301c)

                                                  (301d)
This form of p contains the implicit assumptions that the earth acts as a
reflective barrier to the pollutant particles, that the turbulence is homo-
geneous, and that no low-level  inversions exist.  Assumptions (SOlc) and
(SOld) are not generally valid, but are acceptable approximations over the
relative small areas, comparable to a grid cell, in which c" is finite.  To
maximize the accuracy of Assumption (301b), we plan to use the "optimal"
a (t) and a (t)  profiles developed in Chapter II.
 X         Z

a.   The Point Source Contribution, P

     Using the delta function's property
/  f(x)  6(x -  XQ)  dx = f(xQ)
                                                            a
-------
                                                                           151
                /•*

«(r-'V = Spa(t)J  Pxy(x'y'WT)
                                                  'T) dT    '
                                                               (302)
where
Pxy(x,y,Xa,ya,T) = [27rax(T) ay(i)
                                             exp
                                                     o

                                        (x - Xa - UT)



                                            2a2(T)
                                               v '
                               (y - ya -


                                   2a2(i)
                                      v
                                            2
                                                              (303)
and
 pz(z,za,i) =  /2T
-1
exp
r ( }
az
+ exp
(Z + Za) ]
2a2(r) J
                                                              (304)
We should mention that the assumptions of a reflective earth and the absence


of an upper level inversion can be relaxed by substituting the appropriate


expression for p  in Eq.  (302).  Equation (302) is the first of the set that


we need to model c".
b.   Line Source Terms
     We first assume that all  line source segments are straight lines lying at


ground level  (i.e.,  z  = 0,  3  = 1,2, ..., M) and that the strength of each seg-
                     P

ment is constant along its length.   Suppose that the 3-th line source has


length I   that it is centered at  (xn, yj, and that it makes an angle 0_ with
        p                           P   P                               p
the y-axis as shown  in Figure  38.
    FIGURE 38.   PARAMETERS  REQUIRED TO  SPECIFY  A  FINITE  LINE SOURCE

-------
                                                                          152
The line source term L   which represents this source in the c" formula
[Eq.  (299)J, can be shown to be [see Lamb and Neiburger (1971)]
                  0.                                    _               C
                  r     (   [(x - UT - x ) cos 0R+  (y - VT - y ) sin 0 ]
Mr.*) = Slo(t) /  exp {-
                                             2a
                  -1
                      erf
                           (y  -  VT  -  y  )  co.s  0  -  (x -  UT - x«)  sin e  +
                           	P	P	P	p
         - erf
                    •-                                         K
               (y - VT - y )  cos 0  - (x - UT - xj sin 0Q - -*•
           PZ(Z,0ST)
                                                                         (305)
where a = ax =  ay and  where pz is  given by Eq.  (304), with
the general  line source term needed to compute c".
                                                              = 0.  This is
 c.   Volume Source Terms
     The volume source terms V  arise from S, each volume source having the
 dimensions Ax • Ay • AZ of the grid cells in the urban-scale model.  It is
 not difficult to show that the volume source term V  is given by
 V  (r t) = 	^
 Y^' ^   8AV
                       ^-  X   -  UT  +  ^-\       /X - X  - UT
                     	1	^-j _  erf '	^
                               - VT
                    erf I
                             /2
                                          - erf
                                                fy - y  - VT -
                                                        /2
            Ferf(A) - erf(B) + erf(D) - erf(c)l
                                                                         (306)

-------
                                                                           153
where
              z - z
                /2 o                      /2 a

                                                                         (307)

                  Z-j    _____               7-4-7  -4-
                  -1  -  O                i. ~ i.  ~ ~7T
          B =
and (x ,y , z )  denotes the center of the y-th grid cell.  Equation (306) is
the volume source term required in Eq. (299) to calculate c".

3.   Design of the c" Model Computer Program--!nteraction
     with the Source Inventory Data Bank

     We have derived the basic equations required to calculate c".  The next
step is to put these results into the form of a subroutine that the airshed
model can call to compute concentrations c at any given point.  As previously
stated, the c" model must not interfere with the operation of the full urban-
scale model, and it must not require the latter to perform special calculations
In short, the c" subprogram must be capable of calculating c" given only the
following information (which is always available in the grid model itself):

         (u, v)      = horizontal components of the mean wind in the
                      vicinity of the receptor,
         z./L       = local stability parameter,
         AX, Ay, -AZ.= dimensions of the grid cells,
         u          = local friction velocity,
         z.         = mixing depth in the vicinity of the receptor,
         *c          = value computed by the grid model for the
                      receptor point,
         (x,y,z)    = coordinates of the receptor.

Thus, if MICRO is the name of the c" model, the call for it by the grid model
might look something like the following:

-------
                                                                            154
     CALL MICRO (X,Y,Z,T,ZlOVL,USTAR,ZI,UBAR,VBAR,r
                    DELTAX,DELTAY,nELTAX,CTILDE,CPP)

where CPP = c"(x,y,z,T).

     The MICRO program must also have access to the  source  inventory,  where
each point and line source segment,  or link, is stored separately.   In partic-
ular, the following source data are  needed to compute c":

          XPS(I)j
          YPS(I)>   = (x,y,z)  coordinates of I-th point source.
          ZPS(I))
          XLS(J)  )   _ (x,y) coordinates of the center of the  J-th  line
          YLS(J)  |     line source (ZLS = 0 assumed).
          THETA(J)   = angle of the J-th line source  (see Figure  38).
          XLNGTH(J) = length of the  J-th line source.
          SP(I)     = strength (mass/time) of the I-th point  source.
          SL(J)     = strength (mass/length/time) of the J-th line  source.

Using all of the  above information,  MICRO computes c" in a  manner  described
roughly by the flow chart shown in Figure 39.  We have completed the steps
labeled SIGMA and CALC in the flow chart; Appendix C  presents listings of
their source programs.  Steps  SEARCH, SORT, and DOMAIN have not  been completed,
and without these modules, the MICRO subprogram lacks the capability of search-
ing the source inventory itself.   However, if the user provides  the relevant
source information, the existing steps SIGMA and CALC can compute  c" at any
given point.

     In its present form, SIGMA consists of two FORTRAN functions,  SIGMAX(T)
and SIGMAZ(T), which are used  by the subprogram CALC.  The  SIGMA functions
calculate a  and  a  for a given travel time T and for given stability, mixing
           X      £-.
depth, and friction velocity conditions (which are passed through  common blocks
to the function programs).  To evaluate the integrals entering in  the  equations
derived earlier for c", the CALC subprogram uses a Monte Carlo quadrature
technique, described theoretically in Appendix D.  This technique  allows infi-
nite flexibility  in the ability of the CALC program to evaluate  c"  both near

-------
                                                                       155
               CENTER J
Using given values of USTAR, ZI0VL, and II,
convert the stored nondimensional  universal
functions OX(T ), ay(-r*), and a^t*) into
the dimensional  forms ox, ay, and  az.
                   \/
Using ox, ay, az, UBAR,  VBAR,  X5  Y,  and Z,
compute the boundaries of the  "domain of
influence," i.e., the region containing all
sources that might affect c"(X,Y,Z,T).
Search the source inventory for all  point
and line sources that lie in the domain of
influence.
                   \/
Store location and strength data for all
sources that may have a significant effect
on c"(X,Y,Z,T);  reject all  other sources.
                   \/
Calculate  c"(X,Y,Z,T)  using  Eqs.  (299),
(302),  (305),  (306), and  the source data
listed  by  the  SORT  step.	
               f RETURN J
SIGMA
DOMAIN
SEARCH
SORT
CALC
  FIGURE 39.    FLOW DIAGRAM' OF  THE  c"-MODEL  ROUTINE  MICRO

-------
                                                                              156
and far from point and line sources.   It also requires minimal  complexity
in programming and permits  control  over the accuracy of the integral eval-
uation.

     Test calculations performed with the MICRO model  in its present form
for the problem portrayed in Figure (40a) are presented in Figures (40b)
and (40c).   Although these  figures  are self-explanatory, we wish to call
attention to the size of  c"  in the immediate vicinity of the line source
in this calculation:  It is over two orders of magnitude larger than c at
the same point [cf.  Figures (40b) and (40c)J.

     One of the main problems that  we anticipate in the implementation of
microscale models such as MICRO into full-scale urban diffusion models is
the incompatability of the source inventory data bases that the two kinds
of models require.  That is, in grid models, such as that developed by SAI,
the source inventory consists of total emissions from each grid cell rather
than emissions from each point and  line source within a cell.  The MICRO
model requires data in the latter form, which is not derivable from the
existing source inventory.

     A logical solution to this problem is to collect raw data in the form
of point and line source emissions  and to store it in this form.  It is then
a simple matter to compute total emissions from any grid cell network for
use in a grid model.  This  approach was employed in a pollution modeling
study of the Los Angeles basin conducted by Lamb and Neiburger (unpublished),
Data from this study consist of the coordinates and strengths of each of ap-
proximately 6000 roadway links in the Los Angeles area for the year 1966.
These data are stored on discs and  can be converted, using simple programs,
into emissions from any desired grid cell network.  We have procured this
data bank and will use it in the final developmental phases of the MICRO
program.  Our chief interest is in  using this data to test the DOMAIN,
SEARCH, and SORT modules, which will allow the MICRO program to interact
with the data bank for the purposes of computing c" at any desired point.

-------
       u =   (2,3)
       *   SLIGHTLY
           UNSTABLE
          CONDITIONS

     MEAN WIND VECTOR
   10-METER
      STACK
 .GROUND-LEVEL
 LINE SOURCE

•-GRID CELL
 (200 x 200 x 20)
(a)   Source Configuration  Showing the
     Conditions  of a Test  Problem for
     Solution Using the CALC Routine
     Listed in Appendix C.
                                                                                                  0.02
                                                                                              0.07
                                                                   -0.01
                            (b)
c" (x 10J)  Computed by CALC for
the Problem Given in Part (a).
Note the large magnitude of c"
near the line source.   The
elevation of the point source
results in  positive c" only at
large distances.
(c)
c (x 10 )  as a Perfect
Grid Model  Would Compute
It for the Problem Shown
in Part (a).
         FIGURE 40.     THE CONTOURS OF c"  ARISING  FROM  A  SINGLE POINT AND LINE SOURCE, AS COMPUTED BY CALC
                                                                                                                   en
                                                                                                                   —i

-------
                                                                                158
The final  version of MICRO will  enable the user of a grid model to restore
point-scale spatial  resolution to his diffusion model at any desired point.
The following sections describe  the DOMAIN, SEARCH, and SORT modules that
we are presently developing.

a.   DOMAIN

     We demonstrated analytically in Section V-B that c" decays'with travel
time, even in the absence of chemical decay.  Let the characteristic time
scale of this decay  be T.  Thus, sources located a distance

                                I > [02 + v2]1/2 T                       (308)

upwind of the receptor point (XR,YR,ZR) have a negligible effect on c".   Fur-
thermore, because of the limited horizontal and vertical dispersion rates, all
sources downwind of (XR,YR,ZR) and those outside a wedge-shaped sector upwind
of (XR,YR,ZR) have negligible effect on c".  By ignoring all such sources, the
MICRO program can optimize its calculation speed.

     We showed in Section V-B that the decay time scale T of c" is approximately
equal to the time required for a point source cloud to grow to the- size of a
grid cell.  Thus, if <*x = axtbx> °.  = avtby» and az = aztbz> we have
                                                                         (309)
                                                                 )
The upwind dimension £ of the domain of influence is therefore

                                  L = VT    ,                            (310)
and T is given by Eq. (310).

     The width of the domain of influence is proportional to a (t), assuming
that a  - a .   We conclude that the domain of influence is the wedge-shaped
      x    y
region shown in Figure 41.

-------
                                                                         159
                                           (XR,YR  --:	
DOMAIN OF
INFLUENCE
                                                              GRID CELLS
                                                                   > x
            (a)  Horizontal View of the Domain of Influence
                                                         •AZ
                                             RECEPTOR
                          =  UT
  (b)  Vertical  Cross Section View A-A of the Domain of Influence
      FIGURE 41.  HORIZO.NTAL AND VERTICAL CROSS SECTIONS OF THE
                  DOMAIN OF INFLUENCE ON c"(XR,YR)

-------
                                                                           160
b.    SEARCH

     Once  the  boundaries  of the  domain  of influence have been established,
the next step  is  to find  and list all  point,  line,  and volume sources  that
lie within it.   As  indicated earlier,  the coordinates  of the i-th  point
source are stored in the  data bank as  XPS(I),  UPS(I),  ZPS(I).  Four param-
eters are  stored for each line source.   By testing  the coordinates of  each
source, we can determine  whether that  source  lies wholly or partially  within
the domain of  influence.

     The speed of the point- and line-source  searching operation  can be opti-
mized if sources  in the data bank are  catalogued according to the  grid cell
in which they  lie.   If this is the case, only  those blocks of data filed in
locations  corresponding to grid  cells  within  the zone  of influence need be
examined.   Determining the volume sources within the domain of influence 12
a simple matter,  since these correspond to the grid cells.  Note  that  even
if only the edge of a volume source (grid cell)  falls  within the  domain of
influence, the entire source must be treated.   The  extent to which that source
affects c" at  the receptor is automatically accounted  for fin the  c" equation.

c.   SORT

     The efficiency of the c" calculation can  be increased by considering only
those sources  that together have the largest  influence on c" at the receptor.
For example, the SEARCH step may find  a total  of 10 point and line sources
within the domain of influence,  but only three of these together may be re-
sponsible for  90 percent  of the  observed magnitude  of  c".  In this case, neg-
lecting the other seven sources  would  cut computation  time by 70  percent at
the cost of only a 10 percent error in  c".

     Because the  SORT process itself consumes  computing time, the  usefulness of
this step  is limited to those cases where there  are many point and line sources
within the domain of influence.   To determine  which sources to neglect, we need
an algorithm that permits rapid  estimates of the approximate effect on c" of a
given source.   The  following is  such an algorithm.

-------
                                                                          161
     To  each  point and line  source,  we  assign  a  number RANK.   For the i-th
point source,  RANK is  given  by
                                                                         (311)
where

                  d  -  [(XPS(I)  -  XR)2  +  (YPS(I)  -  YR)211/2     .           (-'31 2)

For the j-^th  line source,
                RflNK =  _.              -                     (    ,
                                " +
where
THETA(J) + TAN'           .           (315)
                  d =  [(XLS(J)  -  XP)2  +  (YLS(J)  -  YP)2]     ,              (314)
                  w =  XLNGTH(J)*SIN
     With a value of RANK assigned  to  each  point  and  line  source,  we  now  form'
the sum

                             RANKT  - SRANK     ,                           (316)

where the summation  is  over all  point  and line  sources  and where we arrange
the sources in  order of decreasing  values of RANK,  with point  and  line  sources
mixed as  their  values of RANK dictate.

     In the CALC  step,  the sources  are treated  in order of their RANK values,
the largest being treated first.  As we proceed through the list of sources,
we keep a running sum of the RANK values  of the sources treated.   Let this  sum
be denoted by

-------
                                                                           162
                                  K
                    RUNSUM(K)  =  ]T   RANK(k)     ,                        (317)
                                k - 1

where K represents the number of sources out of the total  of KMAX point and
line sources in the domain of influence that have been treated.   Note that,
by definition, RUNSUM (KMAX)  = RANKT.   Thus, when

                          RUNSUM(K)                                      ,.1R>
                            RANKT   - B                                  (6lti)

where 0<3<1, we terminate the calculation of the point and line  sources con-
tributions to c" and neglect  all sources with  labels between K and KMAX.   The
value assigned to 3 depends largely on the estimated overall accuracy of the
model; it is generally unnecessary to set 3 =  1 because errors in the model
itself do not warrant the added computing time required.

     A print-out of the list  of the sources actually treated and their RANK
values will  also be of aid in assessing the polluting potential  of each
source in the vicinity of the receptor.

D.   DEVELOPMENT OF MATHEMATICAL METHODS FOR MODELING SUBGRID-SCALE
     CONCENTRATION VARIATIONS c" OF NONLINEAR  POLLUTANTS

     In the  last section, we  derived solutions of Eq. (292), which governs
the c" distribution of linear pollutants.  Needless to say,  obtaining solu-
tions of the more general equation [Eq. (291)], which governs nonlinear species,
is a much more difficult process.   And modeling the complete system of photo-
chemical  pollutants is even more complex.  We  have  only begun to consider the
nonlinear problem.  We outline here a method that we have  developed and are
considering  for constructing  a microscale model o-f nonlinear pollutant species.
We demonstrate the method in  the context of Eq. (291) because it is mathemati-
cally much simpler than the equations  describing the complete, kinetic mechanism,
but yet it retains the essential nonlinear property of those equations.

-------
                                                                         163
     The  obvious  approach  to  solving  Eq.  (291)  is  to  express  the equation
in finite difference  form  and to  integrate  the  resulting  expression  numer-
ically.   However,  this  method only  creates  a  new subgrid-scale  problem:
because due  to  the  presence in  S" of  delta  functions  arising  from point  and
line sources,  subgrid-scale concentration variations  will  exist in any  dis-
crete analogue  of Eq.  (291).   Even  if it were possible  to  achieve a  grid
network  fine enough to  resolve  the  major  variations  in  the concentration
field, operational  problems would still remain.  For  example, because of
the arbitrariness  of  the locations  and orientations  of  point  and line
sources,  a variable grid network would be required to optimize  the accuracy
of the numerical  integration  process.   But  this might lead to difficulties
with computational  stability, boundary conditions, truncation error, and the
like, and would certainly  complicate  the  computer  program.

     To  circumvent all  of  the problems, we  develop in this section a rela-
tively new type of pollution  modeling equation.  This model is  not based on
any radically new theory.   Rather,  it is  based  on  an  equation that is an
intermediate step in  a  sequence of  mathematical operations that leads from
the general  Lagrangian  equation of  turbulent  diffusion  Eq.  (16), through the
well-known Gaussian puff and  plume  formulas,  to the  classical diffusion  equa-
tion [Eq. (291)J.   One  of  the unique  features of the  new  equation that  enables
it to avoid  the problems listed above is  that it is  continuous  in space  (like
the diffusion equation), but  discrete in  time.  In the  remainder of  this sec-
tion, we  derive the new model,  briefly discuss  its place  in the scheme  of
existing  modeling equations,  and outline  a  set  of  ways  in  which this model
can be implemented in  the  calculation of  c" of  nonlinear  pollutants.  The
actual implementation  process will  be undertaken in  a later study.

1.   Derivation of a  Discrete-Time, Continuous-Space  Diffusion  Equation

     The  starting  point of the  derivation is  the Lagrangian form of  the  dif-
fusion equation introduced in Chapter II, i.e.  ,

-------
                                                                           164
 = ff  p(r,t|r',t')  S'(r',t')  dt1  dr'  + /p(r,t r' ,tQ) dr'

              0                                                          (319)

where S is the source strength  function.   In  cases  where linear chemical  decay
occurs at a rate

                                     a£=-k(t)c    ,                     (320)

where k(t) is the rate "constant,"  it can  be  shown  [see, for example,  Lamb  (1971)]
that the proper form of the kernel  p entering in Eq.  (319)  is


                                        r   t          i
       p(r,t|r'Jt1)  = p'(r,t r',f)  exp - I   k(t")  dt"
                                        LJt,           J
                                                                         (321)

where p1 is the probability density that the  fluid  will advect  a particle
from (r',t') to (r,t).  Note that p1  is a  property  of  the fluid motions  alone
and is independent of the chemistry.   Through a series  of assumptions  and math-
ematical operations, it can be  shown that  the Gaussian  plume formula,  the puff
model, and the classical  diffusion  equation are all  derivable from Eq.  (319).
It is in the derivation of the  diffusion equation that  the  discrete-time,
continuous-space model arises.
     Thus,  let the time axis be divided into equal intervals of length At, and
    t  denote the time t - nAt.  For br
the  derivation by writing, for example,
 let t  denote the time t - nAt.   For brevity,  we  omit the reference to  time  in
Consider now the conditional  probability density

                             p(rnlrn_r  rn_2,...,r0)    .                "(322)

-------
                                                                          165
that a particle  is  at  r  at  time  t   given  that  it  was  at  r  ,  at time t  , ,
                      ~n          n  -                ~n- 1           n- 1
at r  ? at  time  t _?,  and so forth.   If  an  interval  At exists  such  that
    ML.          lie.
                       p(r  |r   ,,r   0,...,rn)  =  p(r  r  ,)     ,           (323)
                       p  ~n'~n-l  ~n-2'    '~0     K  ~n  ~n-l      '           v    '

then the particle  motions are  said  to  consitute  a  Markov process  in  the  dis-
crete time frame  because  any stochastic  process  whose probability density
satisfies a relationship  of the  form of  Eq.  (323)  is  called  a  Markov process.
An important property  of  a Markov sequence  is  that its probability density
p(r |r ) satisfies the Chapman-Kolmogoroff  equation
  ~n ~m

                      p(rjrnl)  =/p(rnlr£)  p(r£  rj dr£    .              (324)

where n>£>m are any integers.

     The physical  significance  of the  Markov  process  in the  context  of a tur-
bulent diffusion  problem  warrants comment.   Suppose that the turbulent eddies
that are responsible for  the random motions  of the pollutant particles have  a
time scale T.   A  reasonable estimate for T  is  the  Lagrangian integral  time
scale of the turbulence defined  by
                                       (t) v  (t + s) d£   '                (325)
                                  0
where v'(t)  is  the velocity  component at  time t  of  a fluid  particle.   (If the
turbulence is  stationary,  then  T. is  independent  of  the time t.)   The  time
scale T.  can be regarded  as  a  rough  measure  of the  lifetime of the so-called
energy-containing  eddies.

     Suppose now that  a particle is  released in  the turbulence and that its
velocity is  recorded at intervals  At = T. /10 apart.  We can expect that if
the velocity is, say,  positive  at  times t ,  t +, , and t +?, then there is a
greater chance  of  observing  a  positive velocity  at  time t + ~ than a negative
velocity;  and  conversely,  if the velocity is negative at the first three times,

-------
                                                                             166
then it  is  more  likely  to be negative than positive at time to-  In con-
trast,  if the  particle  velocity is recorded at intervals At = 10 .  apart,
then we  can expect  the  velocity at time t + ~ to be independent of that at
t  ? and all previous times.  The idea here is simply that a particle has
a "memory"  of  some  length T in the sense that the behavior of the particle
at any  time t  is  seemingly influenced by its history during the interval T
just prior to  the time  t.  This apparent memory is imparted to the particles
by the  organized  or systematic behavior that each eddy exhibits during its
lifetime.   It  is  the randomness of the birth and death processes of the
eddies  and  of  the transferral of a particle from one eddy to another that
causes  a particle ultimately to "forget" its distant past history.   This
type of reasoning leads  us to expect that the probability density of the
particle positions  will  satisfy the Markov condition [Eq. (323)] if the
discretization time interval At is chosen such that
                             At » T.     ,                              (326)
where T,  is  given  by  Eq.  (325).
       ,

     For reference, we  rewrite here the general equation [Eq. (319)] using
the notation  adopted  above:

                r        ,               r/"-'
    =/p(rrl r0) dr0 +11    p(r,lr'> S(r',t')dt' dr'
                                        ^
                                                p(r   r') S(r',t')  dt'  dr
                                                                         (327)
The reason for the split integral  in  this  equation  will  become  clear  later.

     For any time t1 <  t  ,,  we have  from  the  Markov  assumption
                        n~ i

                           p(rnlrn_,-r') =p(rJv,)     •

-------
                                                                         167
Thus, from the general  relationship
f(X  X )  =
   n  m
                                X£'Xm>
                                                          (329)
where X ,  X,,, and X  are any three random variables  with probability density
       11    XL-       111

functions  f, it follows  from Eq.  (328)  that
p(rn r1)  =
                                                                         (330)
Thus, p(r  r')  satisfies  a Chapman-Kolmogoroff type of equation,  even  through

r' can be the particle position anywhere on  the continuous  time axis  prior to


Vr



     Using Eq.  (324),  we  can express  the first term in Eq.  (327)  in  the form
           r0>
                                                                         (331)
Similarly, Eq.  (330)  permits  the second term in Eq.  (327)  to be written as
 jj"   p(rnlr')  s(r',t')  dt-  dr-  =  f p(rnlrn_-,)

   L/-N
                                          /-/n-l

                                        '         P(rn_-,  r')  S(r',t')  dt'  dr1


                                        ^ tn
                                                                         (332)
Combining Eqs.  (331)  and (332)  and comparing the result with Eq.  (327)  for


/c(r,t  ,)") ,  we conclude finally that,  by virtue of the Markov property,
\  ~  n-1  '

Eq. (327)  is  equivalent to

-------
                                                                             168
          ^p(r,tn|r',tn_1)dr'
                               P(r,tn  r',t')  S(r',t')  dt1  dr'     .        (333)
                          'n-1
This is  the  discrete-time,  continuous-space  diffusion  equation.
     From the  modeling  standpoint,  Eq.  (333)  has  several  very  useful  assets.
The first is that  it  allows  concentrations at time  t  to  be  calculated from
those at time  t-At and  from the  source  emissions  that  have  occurred during
the interval At just  prior  to  t.  This  is  in  contrast  to  Eq.  (319)  or to the
puff model,  for example,  which  requires  an integration over the  entire period
t~ - t to obtain the  concentrations  at  time  t.  The second  advantage  of
Eq. (333) is that  it  possesses  unlimited spatial  resolution.   As  we pointed
out earlier, this  property  is  essential  in modeling the  microscale  concentra-
tion distribution. The third  advantage  of Eq.  (333)  is  that finite differen-
cing techniques are not required to  evaluate  it.  In  fact,  all of the integrals
in this  equation can  be evaluated analytically when the  concentration distribu-
tion (c(r,t   -,)> is expressed  in series  form.   This advantage  circumvents the
     >  ~ n-i  /
problems, often associated  with  differencing  techniques,  of computational sta-
bility and of the  imposition of  artificial boundary conditions.

2.   Development of A Microscale Model  of  Nonlinear Pollutants Based
     on  The  Discrete-Time,  Continuous-Space  Equation

     Our aim is to use  Eq.  (333) as  the basis of  a  microscale  model.   In this
section, we  demonstrate how this might  be  accomplished,  using  as  the example
a single nonlinear species  c whose  microscale variations  c" are  governed by
Eq. (291).  The formal  relationship  between  Eqs.  (333)  and (291) is quite
complicated, and since  the  details  are  not important  to  the present analyses,
we outline the relationship only briefly below.

     First,  We note that Eq. (219)  is of the  form
                             Lc"  = -kfa + be" +  c"2 )+ S"    ,           (334)

-------
                                                                               169
where
                                a = -c"2 + "  ,                     (335a)

                                b = 2c  »                                 (335b)

Consequently, if S" were zero and if c" were initially uniformly distributed
in space (so that LC" = 0), then the chemistry alone would cause c" to change
in time, and we would have
                           C"(t) = -c + q f-Hr     >                    (336)
                                         LI ~ Q J

where

                               q = (~c2 _ 4a)1/2                          (337a)
                                              »

                                            TTT^-J
d - e-^L|J3 A ,A J                    (337b)
           C0
 and where c!! = c"(0).  Recall the form [Eq. (321)] that the kernel p acquires
 in situations where the scalar quantity decays at the linear rate
We anticipate, therefore, that if a nonlinear decay of the form found in Eq. (334)
occurs, then the kernel p entering in Eq. (335) can be approximated by


        p(r,t|r',t')  =  p'(r,t|r',t')  [-  |j- + ^ (H~ff)]    '              (338)

where, as in the case of Eq. (321), p1 is the solution of the equation

                                Lp' = 6(r - r') 6(t - t1)    .           -(339)

-------
                                                                           170
     We  now  state  without  proof that  the  equation
                                                               dr'
                             pc(r,tn|r',t')
fd +1>
T^Ty
:"(r'
                                           1-1
                                               dt1  dr1
(340)
becomes  equivalent  to  Eq.  (291)  in  the  limit as  At =  t  -  t _,  approaches  zero.
In Eq.  (340),  q  and d  are  functions  of  the  time  and space  variables  that appear
to the  right of  the vertical  bar in  the kernel  p.   Thus, we can say  that Eq.
(340)  provides a consistent representation  of Eq.  (291)  inasmuch as  the solution
of Eq.  (340)  can be made arbitrarily close  to  that of  Eq.  (291)  by making  the  time
step At  suitably small.  We should  hasten  to add,  however, that the  size of At
is restricted by Eq.  (326); so  the  equivalence of  Eqs. (340)  and (291)  can be
                                                                           *
realized only by scaling the time axis  with some value T and then making At  = At/T
suitably small by making T sufficiently large.   The implication of this is that
temporal  resolution is lost as  Eq.  (340)  is made to approach Eq. (291).  Never-
theless,  for all practical purposes, Eq.  (340)  will serve  well  as the working
equation  for a model  of c".

     Consider now the  evaluation of Eq.  (340).   We demonstrate  Our proposed
technique for solving  this equation  by  considering only the first integral on
the right side.   To simplify the analysis  further, we  assume homogeneous,  sta-
tionary  turbulence  and restrict  our  attention to a one-dimensional  problem.
Albeit  highly simplified,  this  case  is  adequate  to demonstrate  the essential
features  of our  technique.

-------
                                                                             171
     Under the conditions  that  we  have  imposed upon the turbulence charac-
teristics, the kernel  p1 becomes
                               (2*)
                                   '"
r
where
                                    x  =  x4  +  uAt
                  (341)
                  (342)
and u and a are  functions  of  x'  and  t   ,.
                                    n-1
     Turning next to  the  definition  [Eq.  (337a) ]  of the parameter q that enters
in Eq.  (340), we see  that.two  of the  variables  upon which  it depends,  namely,
        2
c and c" ,      spatially  averaged quantities.   Since the spatial  region of con-
cern to us  in the microscale model  is only as  large as  several  grid cells, it
is not unreasonable to  treat these  two quantities as constants.   This  assumption
is also supported by  observational  studies of  mean concentration  variations in
both the horizontal  and vertical  directions.   The third quantity  entering into
                                r\
the definition of q--namely, <^c'  )"--is variable  over the  microscale region,
but for simplicity we also regard it  as essentially constant.   Combining all of
the above assumptions in  Eq.  (340), we obtain
             c"(x,tn)  =  -
dx1
where
                                  l
                                                                         (344a)
                              a  -
                                                                         (344b)

-------
                                                                             172
The first  term  on  the  right side of Eq. (343) is clearly equal to c(x,t )


[Eq.  (333)].





     Now let  c"(x,t  )  be  represented by the power series
                      c"<*. V • ^ an>nx»    .                      (345)







The last term in  parentheses in Eq. (343) can now be written in the form





                                      OD

                                1    "C~"^          m


                              = j-   2^  dn-l,mx      '                (346)


                                    m - 0
where
             =  Vi,o +  c + 

(347) (348) m d , = aa , - f- > d , , a , . . (349) n-1 ,m n-1 ,m 3n /^ n-1 ,m-k n-1 ,k v ' u k = 1 Substituting Eq. (346) into Eq. (343), we obtain c"(x,tj + c(x,tj = c(x,tj = —iL— , «n_1>m ^0 1^iro m'= 0 • fx-mexp(- (x- X'/ DAt)' }dx' / - UAI; 2o* (350)


-------
                                                                           173
Changing the integration variable  to

                            5 = x' +  uAt  -  x    ,

we obtain
                                     E  Vi.
                             0    u  m  =  0        •>
                                                                          (351)

where

                                x =  uAt - x    .                           (352)

     Equation (351) can be written  in  the  equivalent form


          c(x,t ) = 	"	V™*  d  -,   *C~*   ryr-—-—r-yj I £  6        d £

                             m  " °        k  ~  °                            (353)

Evaluating the integral in this  expression,  we obtain
                               co
                /  *  \    M   v^  ,1       ^   I/ -\m-2k 2k
               c(x,t  ) = -£-   >    d   i    V  m!(-x).   P   ,?.  _ ni|
                    1     pn   •*-«'  II-I,IM ^^  7 p, \ 17     ?M i *•     'v--     >
                          u      _              ^iLisy:^iii6.i\y;
                             m =  0         k  =  0
                                                                           (354)
where
                             (2k  -  1)!!  = 1-3-5.,.(2k-l)

and where E1'1^ denotes summation over all  k that do not exceed m/2.   Upon
expanding the right side of  Eq.  (354)  and collecting like powers of  x,  and
upon writing the left side of this equation,  i.e.,  c(x,t ), in its series
form
                                    oo
                                   ^—-^^
                        f^ [ V  "f~ 1  —   s    f^   V  -t*f^{Y"f~)
                                  m = 0

-------
                                                                             174
and grouping these terms  with  those of the right side of Eq.  (354), we obtain
the set of equations  for  the coefficients  d   .   The desired  mircoscale field
c"(x,t )  is now in the  form
where
                                       -   a
                                n    m = 0  n'm
                              ansm    m
and the f  are known functions.   The modeling of c" is thus reduced simply to
the algebraic manipulation of coefficients in a series expansion.   Mo differ-
ential  equations or integrals require manipulation.

     As we mentioned earlier, this modeling approach has not been  developed or
tested heretofore.   Consequently, considerable work remains to be  done in pro-
ducing a working microscale model for photochemical pollutants.  We hope,
however, to make considerable progress toward that goal  in our continuing model
development program.

-------
                                                                             175
             VI   SIMULATION  OF BUOYANT  PLUMES IN
                THE PLANETARY BOUNDARY  LAYER
A.   INTRODUCTION

     Up to this  point,  all  of  the  work  presented in this volume has per-
tained to passive  scalar  quantities,  that  is, to scalars that do not inter-
act dynamically  with  the  fluid.  Under  this condition, the mathematical
problems of diffusion modeling are greatly simplified because the effects
of the turbulence  on  the  scalar can be  described in terms of quantities,
(namely, the turbulent  diffusivity) that are properties of the flow.  It
is this assumption that made it possible for us in Chapter II to derive
profiles for the diffusivity Kz that  were  independent of the character
of the diffusing substance.

     In air pollution studies  involving emissions from power plants, oil
refineries, and  other sources  in which  large quantities of both heat and
pollutants are discharged simultaneously,  the "passive scalar" assumption
is not valid on  the microscale,  i.e., over distances up to about 1 km
from the source.  The reason for this is that the heat exhausted along
with the pollutants produces buoyancy forces, which enchance the vertical
transport of the pollutant material.  The  result is the so-called plume
rise phenomenon, which  has been  studied analytically by Briggs (1963) and
others.  In most of these investigations,  attempts were made to describe
the plume rise using  semi-emm'rical formulas.  Some of these formulas
were found to work reasonabl   /ell under certain special circumstances,
but not in others. Apparent   „  there is no formulation that is applicable
to a wide range  of condition^.   In this chapter, we outline a model that
we hope will  fill  this  need.

-------
                                                                             176
     The starting point of our simulation is the Lagrangian diffusion
equation introduced in Chapter II [Eq.  (16)]:
          
-------
                                                                            177
Thus, p is a function of the horizontal projection of the distance to the
release point—but not of the positions of the source or receptor separately-
and of the travel time T—but not of the clock time t or release time t'
separately.  Consequently, we can calculate p in the following way.

     For a given release height z1 , say z1 - z-,, we release a particle
from each grid point on the horizontal plane at height  z, at the initial
instant t  when the data begin.  If there are N grid points in each hori-
zontal plane, then N particles will be released.  We allow each particle
to move forward in time under the action of the forces exerted upon it by
the turbulent fluid.  By tracking each particle, we create an ensemble of
N particle trajectories, each of duration T, from which the density function
p(£,n,z, z-, ,T) ,  0 < T < 'T, can be computed by conventional techniques [see,
for example, Bendat and Piersol (1968)].  Here, T denotes the length of
the interval during which data from Deardorff's model are available.  In
our early work, T corresponds roughly to 1 hour in real  time.

     For the function p derived by this process to be applicable
to a given dispersion problem—for example, the diffusion and transport
of pollutants from a stack of height H, diameter D, exhaust velocity V ,
                                                                      A
and temperature 0 --the ensembel of particle trajectories on which the
                 o
function p is based must adequately reflect the behavior of particles
released under such conditions.   In the example just cited, it is clear
that wind shear, buoyancy, and other phenomena have strong effects on
the particle motions.  Furthermore, as a result of the proximity of the
particles to one another at their time of release, -particle interactions
affect drag forces, heat exchange between the particles, and the atmos-
 phere,  and similar  factors.   Thus, unless  all  of  these  mechanisms  are  pro-
 perly modeled  in  the  simulation of the  ensemble of  particle  trajectories,
 the  density  function  p  based  on this  ensemble  will  lead to erroneous  es-
 timates of the  mean  concentration when  used  in  the  <(c>  equation.

     Hinze (1959) has reveiwed the theory relevant to the prediction of
the motion of discrete, bouyant particles in a turbulent fluid.  In the
remaining sections, we use this theory to formulate a model that, using
Deardorff's data as input, can simulate the trajectory ensembeles required

-------
                                                                         178
to study the diffusion of any type of pollutant (e.g.,  heat, moisture,
gases)  released under any conditions in the planetary boundary layer.
C.    MATHEMATICAL BASIS OF THE TRAJECTORY SIMULATION MODEL--
     PRELIMINARY DESIGN
     The Particle Momentum Equation
     Let (usv,w) represent the velocity V of the scalar (pollutant)  particle
with respect to a fixed Eulerian frame, and let (u ,v ,w )  represent the
average (spatial) velocity V  of the fluid in the immediate vicinity of
the given particle, i.e., the local  environmental  velocity.  We distinguish
between the "environmental" velocity V  and the ambient fluid velocity Vf,
                                     "** C                                "" I
because, as a result of the presence of large numbers of scalar particles,
the local velocity of the fluid that a particular particel  "sees" will
generally differ from Vf.  Also, let T and T  represent the particle and
environment temperatures, respectively.  Provided that the  particles are
sufficiently small, the momentum equation governing their motion can be
written in the form [see Hinze (1959)]
    dt
where
                    [—1
/(a,t)  -  Ve(x(a,t),t)   -
Y Hf Ye(X(a,t)
                                                                       (356a)
                        3   =
                                  9v
                                 +
                                     •')
                                                                       (356b)
                     (356c)
                                                                      '(356d)

-------
                                                                          179
and
          r  =  the particle radius,
          p  =  the particle density,
          v  =  the kinematic viscosity of the environmental fluid,
         pg  =  the density of the environmental fluid,
          a  =  the release points of the particle in question,
     X(a,t)  =  the position of the particle at time t
          g  -  the gravity vector,

and

                        at  •  It  + v • '    •                    <357>

The second term on the right side of Eq. (356a) represents a Stokes drag
force on the particle.  For this approximation to be valid, the Reynolds
number of the particle with respect to the particle radius r and relative
velocity |V - V   must be on the order of unity or smaller.

     The last term on the right side of Eq. (356a) represents the accelera-
tion of the environmental fluid relative to the axis of the moving particle.
In cases such as those of interest to us, where the scalar and fluid par-
ticles are dynamically quite similar, accelerations of the one relative to
the other are sufficiently small that the last term in Eq. (356a) can be
neglected.  We assume further that the response time a   of the particle
is so  small, relative to the time scale of changes in V  that the quasi-
steady-state relation

                              ~nrV(a,t)  -  0                      (358)

holds at all times.  Under these two assumptions, we obtain from EQ. (356a)

                     u  =  ue    ,                                 (359a)

                     v  =  v     ,                                 (35Qb)

-------
                                                                        180
                                                                   (359c)
In view of the  absence of horizontal  forces  on the particles  other than
those exerted  by the  fluid,  we assume that the horizontal  components  of
the velocities  of all  particles are equal  at all  times to  the corresponding
velocities of  the local  ambient fluid.   That is,  we assume that
                             u  =  uf
                                                                   (360)
                             v  =  vf    ,
where Uf and Vf are the velocities given by Deardorff's data.

     It is more convenient to work with potential temperature than with
densities because the former are given for the fluid by Deardorff's data.
Using the gas law and the definition of potential temperature, i.e.,
                                     R/cp
where p is pressure in millibars, we can convert Eq. (359c) into the
form
                                                                   (361)
 Later, we derive the equation governing the particle temperature 0 that
 can be used in conjunction with Eq. (361) to obtain the particle's vertical
 velocity component w.  However, to render Eq. (361) solvable, we must first
 express we in terms of the known ambient fluid velocity component Wf.

-------
                                                                          181
     Toward this end, we consider two limiting cases.  First, in the limit
as the travel  time becomes large and the initial cloud of scalar particles
becomes so dispersed that each particle is surrounded mainly by ambient
fluid particles, we have
                           lim wo = wf    .                        (362)
                           t-*»  e    T
Recall that we represents the velocity of the environmental fluid, com-
posed of both scalar and ambient fluid particles, in the immediate vicinity
of the given scalar particle.  At the other extreme, where the travel time
T is near zero, the entire collection of scalar particles can be envisaged
as a buoyant jet of diameter D (equal to that of the source diameter)
issuing into the ambient fluid.  The velocity of each particle can be as-
sumed to be equal to that of the jet and can be approximated by empirical
and theoretical means.  For the present, however, we assume that initially
the particles form a spherical cloud of diameter D.  This approximation
might apply under very unstable conditions.  The initial velocity of euch
particle is now assumed to be equal  to that of the entire cloud and is to
be calculated by assuming a drag force on the cloud of the form
                              -
                                3/CD\2
where CD is the drag coefficient of the spherical cloud and uro is the
speed of the cloud relative to the ambient fluid.  Invoking the quasi-
steady-state assumption as before, we can equate the drag force D and
the buoyancy force B, where
                                  0-0
                          \B\ = 9           ,                     (364)

-------
                                                                            182
to obtain

                                  2
                          (w - wf)
                                T
or
         w
                     = wf + sign(s) P^1         .                (365)
                        T             JL
Upon substituting this expression into the left side of Eq.  (361), we
can find the effective environmental velocity we in the limit as T -> 0.
We obtain
                 4Dg|0 - 0J\1/2   o 2  /0 - 0
f          -
Tim w  = w  + Sign(e - 0                   -       --      .    (366)
This, then, is the environmental velocity component we in the initial
instance where the cloud of scalar particles behaves as a sphere moving
through the ambient fluid.  From the two limiting cases given in Eqs.  (362)
and (366), we surmise that we can be expressed in the general form

                      we = wf + <|>(T) G(e,0e)     ,                 (367)
where
                         4Dg|0 - 0J\1/2   .2   /9 -
  G(e,ee) = sign(0 - e)                 -       -l           (368)
and where (0) = 1, Hm (T) - 0.   (In the general case,  G(0,0e)  represents
the velocity of the buoyant plume as it emerges from  its  source.)   The  rate
at which 
-------
                                                                          183
cloud of particles becomes dispersed,' or, in other terms, how rapidly the
population of particles surrounding any given scalar particle becomes satur-
ated with ambient fluid particles.  We attempt in a later section to derive
an explicit expression for 4. based on turbulence properties.  Thus, from
Eqs. (361) and (367), we arrive at the final form of the expression for w:

                                      ? 2
              w = wf + <(,(T) G(0,0e) + |yL (0   Qe)    .           (369)

The next task is to derive the equation governing the temperature 0.

2.   The Particle Temperature Equation for Dry Plumes

     As long as there is a difference between the particle temperature
0 and the ambient fluid temperature 0f, not only is a buoyancy force
exerted on the particle, but also a heat exchange between the particle
and fluid acts to eliminate any temperature differences.   This heat flux
can be expressed quantitatively by

                         q = -FA(0 - 0e)    ,                     (370)

where h is the heat transfer coefficient of the particle and A is the
particle's surface area.  The temperature change resulting from the heat
flux q is

                        d0     FA(0 - 0e)
                           =~               '                    (371)
where Vp is the particle volume and Cp is the specific heat capacity (in
units of energy/mass/°K) .   For spheres in the Reynolds number range
1 < Re < 25, F has the form (Kreith, 1966)

-------
                                                                             184
                           h - acu"p     ,                       (372)
where
                           ^r   (2ur)l/
                                     CO


Combining Eq's.  (372) and (371), we obtain
                                                                  (374)
     Here it is necessary to express Oe in terms of 0f.  For this purpose,
we consider the two limiting cases that we considered earlier in relating
we to Wf.  First, in the limit, as the travel time becomes very large, 0e
approaches 0f because each scalar particle "sees" an environment composed
almost totally of ambient fluid particles.  We can therefore write
                           lim ee = 0f    .                        (375)
In the other limit, when T - 0, we envisage as before that the entire
collection of scalar particles comprise a spherical cloud of diameter D
and that this cloud exchanges heat with its environment at the same rate
as a rigid sphere would under identical conditions.  (In the general case,
we envisage the particles as comprising a buoyant jet, but we do not con-
sider this case here.)  Moreover, we assume that internal mixing in the
cloud is sufficient to result in an equi -partitioning among all of the
particles of the total heat exchanged.  In this case, we have
                               (¥} <« - 'f>
limf =-h'"  ' •     3y     ,   .             (376)
               pcf'V'"3

-------
                                                                            185
or
                        o
                          dt
                                                                  (377)
For spheres with Reynolds numbers in the range 25 < Re < 105, the heat
transfer coefficient is given to good approximation by (Kreith, 1966)
                       h =
0.37kf /u D
~D
                                       0.6
                                                (378)
where kf is the thermal conductivity of the fluid (units of energy/time/
meter/°K).  Combining Eqs. (377) and (373), we obtain
                  1 im
                  T+0
                                       (0 - 9f)
    dt
        1.4  0.6
            v
                                                (379)
and by equating this expression and Eq. (374), we obtain the effective
environmental temperature 0e at the initial instant of release:
       lim 0  = 0
       T+0  e
1  +
                            2.2rkf(e - 0f)
    3aPC0(w-wf)°-4D1-4v°-6
                                                  -1
                                       (380)
     Since both we and 0e are uniquely determined once the populations of
scalar and ambient fluid particles surrounding any given scalar particle
have been prescribed, we assume that the transition of 0e from its initial
value [Eq. (380)] to its final form [Eq. (375)] is describable in terms of
the function cf> used earlier with we.  In particular, we assume that
           1 +
                     2.2rkf(0 - 0f)
                                           -1
                                 - (T)]0
                                                           f
                                           (381)

-------
                                                                          186
This equation, together with Eqs. (374),  (367), and  (369), formsia closed
system of equations that can be solved for the temperature 0 of  the par-
ticle and vertical velocity component w as a function of travel  time T
and the ambient fluid conditions 0f and Wf described by Deardorff's data.
In short, we now have model equations that we can use to simulate the tra-
jectories of dry particles of any type in a turbulent fluid.   (The deriva-
tion of the corresponding equations for vapor plumes is currently in progress.)
All that remains to complete the model is to derive a suitable functional
form for C|>(T) .

3.   The 4>(i) Equation

     As defined, (r) has the limiting values

                           lim ({.(T) = 1    ,                      (382a)
                           T+0

                           lim (T) - 0    .                      (382b)
                           T+0

In  physical terms, <{>(T) represents the fractional population of  scalar
particles within a distance of several radii of a given scalar particle.
This suggests that  is equivalent to the mean particle concentration in
an  expanding cluster.  Under the action of turbulence in the inertia! sub-
range, a cluster of particles expands at  a rate given by
                           = 20  2/3
                             30
where e is the energy dissipation rate per mass of  fluid,  £Q  is  the  ini-
tial diameter of the cluster, T is the travel time,  and 
-------
                                                                            187
The latter condition may not always be fulfilled in some problems of
interest to us where significant turbulent energy is created by the par-
ticle cloud itself.  Indeed, observations reveal that under neutral and
stable atmospheric conditions, the turbulent energy within the plumes pro-
duced by large power plants is often significantly greater than that in
the atmosphere just outside the plume.  Under these_conditions, it will
be necessary to resort to the available empirical and theoretical know-
ledge of buoyant jets to derive a saitable approximation for <(>(T).  Never-
theless, in those cases in which Eq. (383) holds, the mean cloud diameter
d at time T is given by


                 d(T) . To2 + 10 e2/3 D2/3  rfK    _
where D denotes the initial  cloud diameter.   This expression holds  only
as long as d < 10D.

     If the initial concentration of particles in the cloud is  unity,  it
follows from the considerations presented above that
                                  4   D3
                                  I* 8
                                                                  (385)
We assume that by the time <£/  is outside the inertial  subrange, or
d > 10D,  and Eq.  (385)  ceases to hold, either the mignitude of cj> has fallen
to a value near zero or ($2)^ ^ AX, where AX is a measure of the grid size
in Deardorff's model.  In the latter case, the growth rate of the cloud can
be estimated explicitly from the computed particle trajectories.

-------
                           VII  SUMMARY

     For the reader's convenience,  we summarize here the major points
presented in each of the six  previous chapters.

A.   CHAPTER I

     In this chapter, we introduced the term "microscale" to refer to all
phenomena whose space or time scales are too small  to be resolvable ex-
plicitly and deterministically in grid models of urban pollution.   Such
small  features result because of turbulence and finite differencing
techniques.   Those produced by the  former include the turbulent velocity
fluctuations themselves, usually denoted by u1, and the concentration
flucutations c' generated by  u1.  Both of these microscale features af-
fect the spatial  and temporal behavior of the evolution of the mean con-
centrations  of photochemical  pollutants.  No exact mathematical expres-
sions  have yet been found that can  describe these effects in generalized
situations.   Consequently, an important problem in pollution modeling
is the development of an approximate description of these phenomena
that,  under  most conditions of interest, will be as accurate as the
data upon which the model calculations are to be based.  Chapters  II,
III, and VI  address these problem areas.

     The use of finite differencing techniques in diffusion modeling
gives  rise to the microscale  phenomenon that we have termed the SSV, or
subgrid-scale concentration variation.  This feature occurs because the
spatial variations that exist in the concentration distribution near
point  and line sources, such  as power plant stacks or highways, are of
much too small a scale to be  resolved by the grid meshes used in urban
pollution models.  We showed  that the SSV can affect the grid-averaged
concentrations of photochemical pollutants in much the same way that the
turbulence-induced concentration fluctuations can  influence the time mean
concentration.  We also pointed out that the SSV can complicate the use

-------
                                                                              189
of pollutant concentration values predicted by a simulation model:
The concentration levels observed at a fixed point will differ from the
spatially averaged concentration predicted at that point as long as the
SSV is not zero.  Model validation studies and concentration extrema
forecasts are two examples of applications in which this problem arises.
Chapters IV and V addressed all aspects of the SSV microscale effects.

B.   CHAPTER II

     In this chapter, we demonstrated how Lagrangian diffusion theory can
be implemented using numerical turbulence models,  Using this technique,
we were able to derive expressions for the distribution of the mean con-
centration of passive material issuing from a continuous point source
in the planetary boundary layer,  We referred to these distributions,
shown in Figures 8(a) and 9(a), as the numerico-empirical (NE) solutions
of the Lagrangian equation.  Inasmuch as the data sets on which these solu-
tions were based were rather small, the results presented in Chapter II
are only tentative.  Future studies will attempt to achieve greater ac-
curacy and will address the important questions of how the concentra-
tions are affected by changes in source height and atmospheric stabilities
(other than those treated here) and by changes to other source types,
such as area or line sources.

     Proceeding under the assumption that the accuracies of the present
NE solutions are at worst comparable to those of available empirical
data, we used these solutions as standards for assessing the accuracies
of the three major diffusion models:  the diffusion equation and the
Gaussian puff and plume formulas.  We devoted subsequent work to optimi-
zing the performances of these models by adjusting the functional forms
of the free parameters in each model such that the resulting predictions
were in closest overall agreement with the NE solutions.  These  "optimal"
parameters, i.e., diffusivities K, and dispersion coefficients a  and a  ,
                                 Z                              Z.      A
are summarized in functional form in Table 2.  With regard to the accura-
cies of the optimal models, we note the following:

-------
                                                                             190
     >   Of the three models,  the  diffusion  equation  is  by far the  superior
        one.   Errors are generally on  the  order of 20 percent,  except at
        points near the source  under  neutral  conditions in which  case much
        larger discrepancies  are  observed  [see  Figures  8(b) and 9(b)].
     >   Relative to-the plume formula,  the  Gaussian  puff equation  is  a
        slightly superior model1,  but,  in quantitative terms,  neither
        provides an acceptable  description  of atmospheric diffusion
        under neutral  stability conditions,  at  least in the problem
        considered in Chapter II.   Both models  tend  to  overpredict
        ground-level concentrations arising  from elevated sources.  Er-
        rors  of 100 percent are prevalent,  and  in isolated areas  they
        are much larger.  In  the  particular  problem  treated,  the  ac-
        curacies of the optimal models  were  acceptable  under  unstable
        conditions (see Figures 18 and  19).

     Considering the Gaussian puff and  plume models  in  general, we note
the following points:

     >   Spatial concentration distributions  in  the planetary  boundary
        layer are decidedly non-Gaussian.
     >   The Pasquill-Gifford  data commonly  used in the  plume  forumla
        are not applicable to emissions from elevated sources unless
        some  allowance is made  for wind shear effects.   When  such  modifi-
        cations are made, they  greatly  improve  the accuracy of  the plume
        formula for elevated  sources  [see  Figure 13(b)s as compared with
        Figure 18(b)].

     The optimal diffusivity, dispersion coefficient, and wind  shear  pro-
files presented in Table 2 were implemented  in  the form of FORTRAN
function routines.  These routines can  be  used  in place of the  corres-
ponding variable names in existing diffusion models  to  achieve  results
that are compatible with the  predictions of  Deardorff's boundary layer
model.

-------
                                                                           191
C.    CHAPTER III

     Using empirical  data, we showed (see Figure 28) that concentration
fluctuations generated by turbulence can dominate the temporal behavior
of the mean concentration of materials that undergo nonlinear chemical
reactions.  Although  several previous investigators have suggested ap-
proximate mathematical expressions for describing these effects, none of
these equations is well  suited to pollution modeling studies because
each introduces too many additional differential euqations into the
system to be solved.   For this reason, we set for ourselves the task
of developing a new approximation that does not entail multiple equations.
The results of our efforts for the case of a bimolecular reaction are
represented by Eq. (95)  for the generalized case, and by Eq. (99) for
the case where the reactants are premixed.

     In Section C of Chapter III, we develop for various situations ap-
proximate functional  forms of the parameters £. and ?R, which enter into
the generalized closure  scheme expressed by Eq. (95).  [See Eqs. (113),
(121), (124), and (125).]  Tests—based on empirical data—of the accuracy
of the overall scheme produced excellent results (see Figures 28 and 30).

D,   CHAPTER IV

     Just as turbulent concentration fluctuations can affect the mean
rates of nonlinear chemical reactions, so can subgrid-scale variations
in the concentration fields simulated by numerical models.  In Chapter
IV, we derived a test of the significance of these effects  [see Eq. (164)],
When this condition is satisfied, SSV effects are negligible and can be
ignored; but, when it is violated, SSV effects may or may not be impor-
tant, depending on the particular physical situation.  To handle these
cases, we used the concepts that we employed in Chapter III to treat
turbulent concentration  fluctuation effects to develop a scheme for
parameterizing the SSV influences.  In its most general form, this scheme
is described mathematically by Eq. (182); and in cases where the reac-
tants are emitted from the same sources,  it takes the form of Eg. (186).

-------
                                                                             92
    We  applied  the  latter  to  three  problems  of  practical  interest:

    >   A  random spatial  distribution  of  point sources,  such  as  building
        heating  emissions.
    >   A  network of streets of  arbitrary separation.
    >   Strong,  isolated  point and  line sources.

From these applications,  we derived  a  dimensionless  number y,  defined
by Eq.  (232),.that provides a  quantitative measure of  the  magnitude  of
SSV effects on nonlinear  reactions  simulated  in  urban  pollution  models.
As portrayed in  Figure  34,  the magnitude  of the  SSV  impact grows as  the
value  of y increases.   Using representative rate  constants, source
strengths, diffusivities, and  the  like, we evaluated the parameter y for
several  of the important  photochemical pollutants  (see Figure  35).
From this  study, we  drew  the following conclusions:

    >   The NO-Oo reaction  is  the  only reaction  explicitly treated in the
        current  SAI  model that is  affected by subgrid-scale concentra-
        tion variations.  In this  case, the SSV  suppresses the effective
        rate of  ozone depletion.
    >   Effects  from freeways  carrying 10 or more vehicles per  day  are
        significant  when  the wind  is parallel to  the freeway,  but not
        when it  is perpendicular.   Smaller traffic volumes produce pro-
        portionately smaller effects,  and larger  volumes produce pro-
        portionately larger ones.
    >   SSV effects  from  networks  of city streets  carrying 10  or more
        vehicles per day  are significant  when (1)  the  streets  are 150
        meters or more  apart and (2) the  meteorology or source emission
        rates are in a  state of  transition.   Under steady-state  conditions,
        the SSV  effects from uniform street networks become negligible,
        but those arising from strong  isolated sources do  not, as described
        in the second point above.

    These findings  are only tentative.   In future work, we will attempt
to corroborate them  through numerical  experiments.

-------
                                                                             193
E.    CHAPTER V

     An objection frequently raised against grid point diffusion models
is  that in spreading emissions from point and line sources uniformly
throughout the local grid cell, such models produce distorted descrip-
tion's of the time concentration field within the immediate vicinity
of point and line sources.   This effect complicates the tasks of model
validation and concentration extrema predictions.   In Chapter V, we
developed a so-called rnicroscale model  that can be used as a subprogram
in any grid model of CO,  S02, or other inert or first-order reactant
to resolve the detailed concentration field at a given point.  This sub-
model has been implemented in FORTRAN and tested in several preliminary
calculations, as shown in Figure 41.

     The chief weaknesses of this microscale model in its present form
are its inability to handle pollutants that react  nonlinearly and its
use of Gaussian kernels.   The latter precludes applications to problems
of the following types:

     >  Those in which aerodynamic effects are important, such as
        downwash in the lee of an elevated roadway or building.
     >  Calculations of pollutant levels in street canyons.
     >  Estimations of pollutant levels on a roadway where vehicle
        wake turbulence is primarily responsible for the initial pollu-
        tant dispersal.

For the purpose of developing a microscale model for pollutants that react
nonlinearly., we derived a special discrete-time continuous-space dif-
fusion equation [Eq. (333)].  We showed that, when the concentration
field is expressed in a series form and line and point sources are rep-
resented in the usual delta function forms, the nonlinear microscale
model reduces to a simple set of algebraic equations (see the last
equation in Chapter V).  Tests of this equation and attempts to relax
the Gaussian kernel assumption implicit in-the current microscale model
versions will be undertaken in the next phase of the EPA contract work.

-------
                                                                             194
F.    CHAPTER VI

     In this chapter, we outlined a method whereby the boundary layer
model of Deardorff and the equation governing fluid particle dynamics can
be combined to simulate buoyant plumes and cooling tower exhaust in the
planetary boundary later.   The idea was to create an ensemble of particle
trajectories from which the probability density function p that enters
into the general  Lagrangian diffusion equation [Eq, (16)] can be derived.
The analyses presented in  Chapter VI are intended primarily as an illus-
tration of the technique we are planning to use to simulate buoyant
plumes.  Further theoretical  analyses will be required to develop the
final form of the model equations,

-------
                                                                           195
                        VIII   FUTURE EFFORTS

     The work reported  in  this  volume  has  attempted  both  to  emphasize
the aspects  of pollution modeling  that are affected  by microscale  phenomena
and to develop mathematical  tools  for  describing  these effects.  However,
two basic aspects  of this  project  are  still  incomplete:

     >  A thorough and  systematic  evaluation of the  magnitude  of
        microscale effects in  specific problem areas of pollution
        modeling interest.
     >  A thorough testing of  the  accuracies of the  various  tools
        that we have developed  for treating  microscale phenomena.

This last deficiency can be remedied  by means of  a series of straight-
forward numerical  experiments  in which each  of the microscale  modeling
techniques is examined  and exercised  under controlled conditions.   We
reported the results of a  few  such studies—the tests of  the closure
scheme presented in Chapter III is one example—and  others are planned
for the next phase of our  work.  It is the evaluation of  the magnitude
of the actual microscale problem area  itself that will  require further
planning and careful study.   Indeed,  the outcome  of  these evaluations
may well alter the future  course of our microscale modeling  work.   We
foresee three basic steps  in the design of these  evaluative  studies, as
discussed in the sections  below.

A.   SPECIFICATION OF THE  DEGREE OF SPATIAL  RESOLUTION
     REQUIRED IN URBAN  POLLUTION MODELING  STUDIES

     Air quality standards are  currently expressed in terms  of certain
time-averaged concentrations,  but  no  mention is made of the  extent of
the spatial  areas  to which these criteria  apply.  Considering  the  fact
that air quality standards are  intended primarily as safeguards of public

-------
                                                                           196
health,  one can infer that these standards apply to all  points where
plants or animals vulnerable to pollution damage are found.   Thus, for
assessments of air quality relative to cropland, livestock,  hospital
patients, and so forth, or for validations of a diffusion model using
data gathered by a pollution monitoring station, an urban-scale pollu-
tion model  possessing point spatial resolution (as described in Chapter
V) should be adequate.   However, to assess the pollutant dosages received
by highway patrolmen ,  bus and truck drivers, road maintenance men.,
and other similar occupational groups, a model should have not only the
capability of providing line-integrated concentrations,  but  also the ability
to simulate pollutant concentrations on the sources themselves.  The
last feature is outside the scope of present urban diffusion models,
because on the urban scale highways can be treated as line sources of
zero width within which concentrations are infinite.  Moreover, on a
highway or within a street canyon, the initial dispersion and the rates
of fast nonlinear chemical reactions are controlled by vehicle wake tur-
bulence and other aerodynamic effects that are not described by the con-
ventional Gaussian puff and plume models or by the commonly  used dif-
fusion equation.

     These considerations and the likelihood of observing intolerable
highway-integrated pollutant dosages point to the need for special micro-
scale models that can simulate pollutant levels in the motorist's
frame of reference.  One of the goals of the next phase  of our micro-
scale work will be to develop such a model to supplement the capa-
bilities of our urban-scale model.  We also plan to examine  further the
resolution required for various types of modeling studies.

B.   EVALUATION OF THE ACTUAL SPATIAL VARIABILITY OF
     CONCENTRATIONS AT RECEPTOR SITES OF INTEREST

     Having established the locations at which pollutant concentration
calculations are needed and the degree of spatial resolution required at
each site,  one is faced next with the following question:  Are local

-------
                                                                            197
spatial  variations in the concentration field so large that an urban-
scale model  alone is incapable of providing a representative estimate
of the true  pollutant levels at those sites?  In the case of the fixed-
point dosage estimates,  this question is essentially one of whether local
sources  are  responsible  for a significant fraction of the pollution ob-
served;  and  in the case  of the roadway-integrated dosage calculations,
the question is whether  background concentrations are so small--compared
with those on and produced by the roadway—that the urban-scale model
itself is needed.

     Continuous pollutant measurement data can help resolve these ques-
tions.  For  example, suppose that such a record is available for a
ground-level traverse of a city.   If moving averages of this record were
made to remove all variations that a grid model of that pollutant could
resolve, then the residual concentrations would represent the micro-
scale variations in question.  If these microscale deviations in populous
areas turn out to be only a small fraction of the total concentration
observed at  those same locations, then urban-scale models alone, i.e.,
without a supplementary  microscale module, should be adequate for assessing
urban air quality.

      Analyses of the type just described have been performed on unpub-
lished oxidant measurements made by Lamb and Neiburger in the Los
Angeles  basin in 1967.  Preliminary results of this work indicate that
significant  microscale oxidant variations occur on freeways themselves.
However, the amplitude of these perturbations decreases very rapidly with
distance, both upwind and downwind, from the freeway.  Microscale varia-
tions in oxidant levels  also occur on heavily traveled city streets and
in street canyons, but the amplitudes are less than those observed on
freeways.  Finally, on all streets with little or no traffic, microscale
oxidant  variations are negligible.

-------
     These preliminary findings  emphasize the need for a microscale model
applicable to roadways themselves.   They also suggest that, except for
sites on streets  or freeways,  it should be possible to obtain accurate
fixed-point assessments of oxidant  levels using an urban-scale model
alone (one in which subgrid-scale concentration variations have been prop-
erly parameterized).   We plan  to continue these empirical studies so
that we can determine where the  true microscale modeling problems lie.

C.   ASSESSMENT OF THE NEED FOR  REFINED MICROSCALE
     TRANSPORT AND DIFFUSION FORMULAS

     The foregoing discussions point toward the need for pollution simu-
lation models applicable to both freeway and street canyon environments.
Since within such regions transport, diffusion, and concentration fluc-
tuation effects are controlled by vehicle wake turbulence and other aero-
dynamic effects not normally considered in conventional  pollution modeling
theories, the question arises  of whether refined diffusion formulations
will be required  to develop the  microscale models needed.  We have al-
ready begun to develop a model of wake turbulence that can be used in
conjunction with  the theory presented in Chapter III of this report to
simulate chemical reactions in the  wakes of motor vehicles.  We also
plan to develop expressions for  the kernel p, which enters in the
Lagrangian diffusion equation  [see, for example, Eq. (295)], to render
this equation applicable to situations where aerodynamic effects, such
as building and elevated roadway wakes, are important.  This work will
be described in subsequent phases of our EPA contract effort.

-------
                                                    199
               APPENDIX A
FORTRAN PROGRAMS FOR USTAR, DKZ, AND UBAR

-------
                        USTAR =
                                   ROUTINE
U10 = wind speed at anemometer height,
Z10 = anemometer height,  in  meters,
ZI = inversion  height,  in meters,
ZO = surface roughness,  in meters,
ZI0VL = ZI/L =  stability  parameter.
   FUNCTION USTAR (U10,ZO,ZI0VL,ZI,Z10)

             IF (ZI0VL)   10,20,30

10             Z10=Z10/ZI
               ZO=ZO/ZI
             IF (Z10.GT.0.025.AND.ZO,LT.0.004)  GO  TO  15
               X=(1.-15.*Z10+ZO)*ZI0VL)**0.25
               XO=(1.-15.*ZO*ZI0VL)**0.25
               A1=ATAN(X)-ATAN(XO)
               A2=ALOG((X-1.)/(XO-1.))-ALOG((X+1.)/(XO+1.))
               U10BAR=(2.*AHA2)/0.35
               USTAR=U10/U10BAR
   RETURN
15
   RETURN
20
             IF  (Z10.GT.0.3)Z10=0.3
              UU-26.22+153. 2*Z10-14.28. *Z10**2+5541. *Z10**3
                    -7523.*Z10**4-ALOG(ZO*6.8E6)/0.35
              USTAR-U10/UU
               USTAR=0.35*U10/ALOG((Z10+ZO)/ZO)
   RETURN
C  EXPRESSIONS  BELOW  FROM  RAGLAND  PAPER
30             ZZ=Z10*ZI0VL/ZI
             IF (ZZ.GT.1.0)  GO  TO  35
               USTAR=0.35*U10/(ALOG((Z10+ZO)/ZO)+5.2*ZZ)
   RETURN
35             USTAR=0.35*U10/(ALOG((Z10+ZO/ZO)+5.2)
   RETURN
   END

-------
                   VERTICAL DIFFUSIVITY ROUTINE
Z = height of level  where KZ is wanted,
ZI = inversion height,
USTAR = friction  velocity,
F = Coriolis parameter  = 2ft sin 0,
ZI0VL = ZI/L = stability parameter.
   FUNCTION DKZ(Z,ZI,USTAR,F.ZI0VL)

             IF (ZI0VL)  10,20,30
10
   RETURN
15
   RETURN
20
  ' RETURN
25
  RETURN
30
35
RETURN

RETURN
END
  Z=Z/ZI
IF (Z.GT.1.0) GO TO 15
  DK=-6.934E-3+0.6H3*Z+3.297*Z**2
       -6.442*Z**3+3.153*Z**4
IF (DK.LT.0.0) DK=0.0
  DKZ=DK*USTAR*ZI

  DKZ=0.6123*EXP(-(Z-1.)**2/.02)

  Z=Z*F/USTAR
IF (Z.GT.0.45) GO TO  25
  DK=7.396E-4+6.082E-2*Z+2.532*Z*Z
       -12.72*Z**3+15.17*Z**4
  DKZ=DK*USTAR*USTAR/F

  DK=3.793E-3*EXP(-(Z-0.45)**2/2./4.E-2)
  DKZ=DK*USTAR*USTAR/F

  XL=ZI/ZIOVL
IF (Z.GT.085*XL) GO TO 35
  DKZ=.35*USTAR*Z/(1 .+4.7*Z/XL)

  DKZ=0.0

-------
                                                                           202
                    WIND SPEED PROFILE ROUTINE
   FUNCTION UBAR (Z,ZO,ZI,ZI0VL,F,USTAR)

             IF (ZI0VL)  10,20,30
10             ZO=ZO/ZI
               2=1/21
             IF (Z,GT.0.025.AND.ZO.GT.0.004)  GO TO  14
               X=(1.-15.*(Z+ZO)*ZI0VL)**0.25
               XO=(1.-15.*ZO*ZI0VL)**0.25
               A1=ATAN(X)-ATAN(XO)
               A2=ALOG((X-1.)/(XO-1.))-ALOG((X+1.)/(XO+1.))
               UBAR=(2.*AHA2)*USTAR/0.4
   RETURN
14
   RETURN
15
   RETURN
20
   RETURN
25
   RETURN
30
   RETURN
35
              IF (Z.GT.0.3) GO TO 15
                UU=26.22+153.2*Z-1428.*Z**2+5541.*Z**3
                     -7523.*Z**4-ALOG(Z0*6.8E6)/0.35
                UBAR=UU*USTAR

                UBAR=USTAR*(32.33-ALOG(ZO*6..8E6)/0.35)

                ZO=ZO*F/USTAR
                Z=Z*F/USTAR
               IF  (Z.GT.0.055.AND.ZO.LT.0.006) TO TO 25
                UBAR=USTAR*ALOG((Z+ZO)/ZO)/0.37

              IF (Z.GT.0.21) Z=0.21
                UU-29.82+213.2*Z-1989.*Z**2+8743.*Z**3
                     -14670.*Z**4-ALOG(ZO*1.5E7)/0.37
                UBAR=UU*USTAR

                ZZ=Z*ZI0VL/ZI
              IF (ZZ.GT.1.0) GO TO 35
                UBAR=USTAR*(ALOG((Z+ZO)/ZO)+5.2*ZZ)/0.37
              IF (ZZ.GT.6.0) ZZ=6.0
                Z=ZZ*ZI/ZI0VL
    WRITE (6,1000)
                UBAR=USTAR*(ALOG((Z+ZO)/ZO)+5.2)/0.37
    RETURN
1000 FORMAT  (1HO, 'NO ACCURATE WIND DATA ABOVE SURFACE  LAYER  IN  STABLE  CASE')
    END

-------
                                                    203
            APPENDIX B
DERIVATION OF EQS, (193) AND  (195)

-------
                                                                          204
                             APPENDIX  B
               DERIVATION  OF  EQS,  (193) AND (195)
     Equation (193)  states that


          Aj(r,t)  =  C f  p(r - r',t,t') S(r',t')  dt1  dr'
                       0

 and Eq. (195) reads  as  follows:
 A?(r,t)  =fff  f  p(r-r',t,t')  p(r-r",t,t")  G(r' ,r" ,f ,t") dt1 dt" dr1  dr"
  1  ~     I / /   I     ~ ~          ~~           ~~
         JJ JQ JQ
     In the derivation  of  these equations, we first want to  show  that
l)  d'  d  =
                                    'J
  Av(r)
                                                                      (B-l)
 where

          "              x-Ax   Ay   AZ

                        -Ax  -Ay  -AZ
                                                                      (B-2)
 Since each of the  three  integrals in  /dr' in Eq.  (B-l)  are  similar, we con-
 sider the one-dimensional case only, i.e.,
 /-X+AX  r                                      X+AX   •
J      J  P(x" - x1) C(x") D(x') dx1  dx"  =  /  /     p(x" - x1) C(x") D(X') dx" dx1   .

-------
                                                                            205
Let



                   x" = x - x1  + x1  =>dx" = -dx1




Then
         X TAX


     f f      p(x - x-j) C(x - x] + x1)  D(x') dx] dx

       w '  A v
                              y.X-.+AX

                              /        C(x - x]  + x') D(x') dx1  dx
                                                                 ]
                             X,-AX
Now let
                      t. = x'  - x,  $>  dx1  =
                r            r Ax
            =  /  P(x - x^  /     C(x + 0 D(x]  + 0  d? dx]
                               Ax




                            -AX




Upon  repeating  for  the y  and z  integrals, we  obtain  Eq.  (B-l)





     Now  we  want  to show  that
1  f \ ff p(r - r1)  p(r - r")  S(r')  S(r")  dr1  dr"  dr
v J   \JJ    ~   ~     ~~     ~      -    ~   ^ j  -

  AV  L
                      =  f f  P(r -  ^')  P(r - r")  S(r')  S(r") dr1  dr"
                        ^yy     ~~




where the product mean value  in  the  right integral  is defined as in Eq. -(B-2)


To prove this,  we start with  the one-dimensional  problem as before:

-------
                                                                             206
  x. X+&X


J     ff
                  - x'} p(c - x"}  S(x') S(x") dx'  dx1'
  X-AX
   l     rx+*x r
= 2& J     y  P(C - x') S(x') F(c)
                                                 dx'
                                                                      (B-3)
                     X-AX
where F(c)  = p(c - x")S(x")dx".   We can use Eq.  (B-l)  to write Eq.  (B-3)

in the form
                      P(C - x1) S(x') F(s) dx1
                                                                      (B-4)
where
S(x')
Now let
                     rAx  /•
                    J    JpU - x"  + 0  S(xr  + 5)  S(x")  dx"
                    -AX
                      x"  = X"  + r  =>  dx"  - dx"
Then
S(x') F(c) =
               2A
                                -  x")  S(x'  + 0  S(x"  + 0  dx"
                   -AX
             = f  P(C -  x")  S(x')  S(x")
                                       dx"
where the tilde average here is defined as in Eq.  (B-2).   Using this result

in Eq.  (B-4),  we obtain

-------
                                                                       207
/ iff  P(~ " ~'}  p(r  - r"} S
-------
                                                    208
             APPENDIX C

FORTRAN LISTINGS OF THE CALC, SIGMAX,
        AND SIGMAZ FUNCTIONS

-------
                                                                          209
   -SUUROUTINE  CALC(XALPHA,YALPHA,ZALPHA.XR„YR,2R.ISTYPE.XLNGTH,THETA»
   V      DELTAX.D EL TAY,DELTAZ,UB/iR,VeAR,RESULT)
    COMMON/M ICF<0/ IT F MAX , NCHFCK. N S A MR , X S AMP , F PSLON , I X , R I TF . S I GF A C
    co MM ON/T A RLE:/ ERFTC isos) ,EXPT{2005 ) .NERF,NEXP,ERFLIM,DERF.DEXP,
   ^  	.	tz XPL I M — •	•	-	- 	  ..._T  	- _,T _	.  ._
    DIMENSION TEST(SO)
    LOGICAL  HALT,RITE.CALL1
                  IOUT  =  ITRMAX/NCHECK •«•  1
  	SQ2 =—1 ,-4 14 21 3-5-	:	—
                  SQ2PI  = 2.5066283
                  T2PI = 6.2821653
                  DELTA-0.0
.___	IF_{ I STYPE )- 20,10,15-
 10               COSTH  = COS(THETTA)
                  SINTH  = SIN(THETA)
                  XLO2 = XLNGTH/2.
                 -DELT-A=XLC2	
               GO TO  20
 15               DXO2  = OELTAX/2
                  DYO2  = DELTAY/2
                     2-= DELTAZ/2,
                  OELTVe - DELTAX*DELT AYfOELT AZ*8
                  DELTA=DXO2
 20               XMXA  = XP-XALPHA
_____ I _ yMYA_=_YR-Y ALPHA ___
                  ZMZA  = ZR-Z^LPHA
                 ' ZPZA  = ZR-t-ZALPHA
                  ZMZA22 = ZMZA«*2/2.
                  ZPZA22 = ZPZA**2/2.
                  R =  SORT (XMXA<--XMXA+YMYA*YMYA)
                  VEL  = SQRTC UBARwOBAR-f VBAR*VBAR)
                  TO =  R/VEL
                  SIGX=SJGf/AX(TO)
                  SGDELT-SIGFAC*SIGX+DELTA
                  TUP=(R+SGDELT)/VCL
                  TLO=
-------
                                                                                           210
      60                   EXPARG=(XUHART*XU0ART-»-YVBART*YVBART)/2.0/SIGX?
                           PXPY=EXPGCEXPARG)/SIGX2/T2PI
      65
                       -GO  TO  75
     70                   Al = < XUBART+DX02) /SQ2SGX
                        •   A2=(XUBART-DXC2 >/S02SGX
                             = ( YV8ART-S-DYC2) /S02SGX _
                            ^= (YVBAKT-DY02 )/SQ2SGX
        *                  PXPY= = T22
                        IF( ITROUT-NSAMP.LT.O)  GO  TO  200
X -- CONVERGENCE- ROUT I NE.
                           SUMJ=0.0
                           JMIN=1TROUT-NSAMP-M
                DO  150  J=JM! K, ITRCUT
-- 150 -- SUM-J = SUM J-«-TEST(-J ) _____
                           SBAR=SUV J/XSAMP
                           SUMJ=0. 0
                DO  160  J=JMIN, ITROUT
-.-4..60 -- SUMJ = SUMJ+ { TFST( J) -SBAR )**2.
                           HALT=SQRT C SUMJ/XSAMP).LE.EPSLON*SOAR
                        IF(HALT)   GC  TO  205
                CONTINUE
                                   = 7RNGE*T22	
                        IF{.NOT.RITE)  PETURN
                DO  280  JTEST=1,ITROUT
                           ITRN-JTEST^NCHECK
 	280	WR-I-TE (-6, 1-04-2-) —I-TRN,T-EST-fJTES-T-)	
         RETURN
  1028 FORMAT(1H0.7E13.5 )
  1030 FORMAT(1H00'RECEPTOR  LIES ON  A  POINT SOURCE')
  104-2- FORM AT UH--.J AFTER' » I6.2X-. MTERAT IONS-TEST= •-» El 2. 5)
         END
         SUBROUTINE  TABLET
      — COMMON/TABLE/-£RFTt 1 505 )-,EXP-T-( 2O05 ) .NERF-.NEXP-sERFUI-H^DERF^OE XP,
        *          EXPLIM
                          Y=0.0
                          N02=NERF/2
                          DERF=ERFLI M/FLOAT(N02)
                          ERFTCN02P1 )=0 .0
                DO  15  N=1,NO2
                          ERE=ERF( Y)
                        •  ERFT( NQ2P1+N) =ERE
     IS                  ERFT{N02P1-N) =-ERE
     - - - DEXP=EXPL-J
                          X=-DEXP
         RETURN
         END
                DO  20  N=1,NEXP
                          X = X + DE:XP
                          EXPT=EXP<-X)

-------
                                                                          211
      FUNCTION EXPO(APG)
    ^COMMON/TABLE/ ERFT(1505),£XPT<2005KNERF..NEXP.ERFLIM.OERF^OEXP.

                    IE=IFIX(/»RG/DEXP)+1
                    IF(I5.GT.NFXP)  IE=NEXP
    -  	  EXPO-=EXPT< IE)- --	       	
      RETURN                                   		 -
      END

    -^-FUNCTION ERR OR (ARC)	.	
    ^COMMON/TABLE/ ERFT ( 1 505 ) , EXPT ( 20 05 ) . NERF , NEXP
     RETURN
     END—
              EXPL1M
                    IERF=IF IX( (ARG+EPFL IM)/DERF
             	-IF( IERF-.LT»t).-IERF=1   	
                    IF( IERF.GT .NERF)  IERF = NERF
                    ERROR=ERFT(IERF)
              TN  5 T_GM A7. ( T~P
       COMMON WS.TAR,USTAR,H,F
   10                TS= T*WSTAR/H
                   IF(TS.GT.Oc525) GO  TO  15
                     SIGMAZ=H*( 0 .65*TS-1 . 89*TS*TS +5
       PFTIIPM	,
   15              IFCTS.GT.l .1 )  GO TO 16
                     SIGMAZ=H*(-1. 11 1*TS*TS + 2.344*TS-0.33.4
       RETURN
     .  RETURN
   20                TS=T*F/0.45
                     USTARF=LSTAR/F
   - LF(TS.GT^O.OS) -GO -TG-2 ______
                     SIGMAZ=0.45*USTARF*0.64*TS
       RETURN
   21              IF(TS.GT.0.175)  GO  TO  22
  	SIGMA.Z=X1 *45JJ-USTJXRF*< 0 .Ol
       RETURN
   22                SIGMAZ=0.45*USTAPF*( 0.058*0 .12"9*TS)
                     IF( SIGMAZ. GT.0.31*USTARF)  S I GM AZ=0 . 3 1 *UST APF
             LIMIT  nnr^  rn M A x ( *.\r,^/u ) = .Q  IN IJM^TAPI F r.A.sF AND H= .
             NEUTRAL CASE.
       RETURN
   30           WRITE (6, 1003)
       RETURN
 1003  FORMATC 1HO.. 'NO  SIGMA DATA FOR STABLE  CASE1 >
       END	   	

            TI 0 N  S I G M A X ( T f t
       COMMON"~wsT~A~Dt USTAK.H. F
                   IF(WSTAR)  30,20.10
      	__TS=-TJf WSTAR/H	
                   IF(TS.GT.O.OS)  GO TO 15
                     SIGMAX=2.2*TS*H
       RETURN
   15
       RETURN
   .20                TS=T*F/0.45
                   IF(TS.GT.O .1 )  GO TO  22
       RETURN
   22                SIGMAX=( .45*USTAR/F)*(0.1-»-0.92*TS)
       RETURN
-- 30 - : -- WR-LTE (6-. 1 003 X -------------
                     S1GMAX=0.5*(USTAR/F)
       RETURN
 1003  FORMATC 1HO» "NO  DATA  FOR  STABLE CASE*)
    -_  END _______________ .

-------
                                                212
        APPENDIX D
 THE MONTE CARLO TECHNIQUE
USED BY THE CALC SUBPROGRAM

-------
                                                                          213
                             APPENDIX  D
                     THE MONTE  CARLO  TECHNIQUE
                    USED BY THE CALC  SUBPROGRAM


     We  propose to use a Monte Carlo technique to evaluate the  integrals
entering  in  Eqs.  (302), (305), and (306).  The technique is best  des-
cribed  by the mathematical analyses that are required  to prove  its  validity.

     Consider a function f(t) defined in the interval  t« < t 5  t- + T
as shown  in  Figure D-l.  Suppose we pick a sequence of numbers  t.,
i = l,2,...,n, at random in  the interval tQ < t 1 tQ + T.   We choose the
numbers  by a random  process  such that each number in the interval is
equally  likely to occur.  For eacl
a corresponding unique number f.:
equally  likely to occur.  For each number t.  in the sequence,  there  is
                            f. = f(t.O
                             i    v i

Thus,  the  probability of observing a value f in the range f.  <  f  -  f. +  f
is by  definition
                                               m(e.)
               p(f.) = prob(fi < f < f. +  f) = -^     ,           (D-l)

where  m(e.)  is  the measure of the set e. in which f.  < f(t) - f.  +   f
(see Figure  D-2)  in the domain (tQ, tg+T).

-------
f(t)
                                                     t0+T
  FIGURE D-1.  A  FUNCTION  f(t)  DEFINED  IN  THE  INTERVAL tQ < t <_ tQ + T

-------
 f(t)
  A
fi
                                        TOTAL  LENGTH  = m(e.)  =  MEASURE OF SET e.
                 FIGURE  D-2.  ILLUSTRATION OF THE  PROBABILITY  THAT f(t)  OCCURS

                  IN THE INTERVAL fn- TO f-j + Af OVER  THE  TIME t0
                                                                                            VT
                                                                                                          rv>

                                                                                                          en

-------
                                                                           21G
     Now, by definition of the Lebesque integral, we have
              ,t0+T
f(t) dt -
                             lim
                             Af->-0
                                                                      (D-2)
where the summation is over all intervals f.- in the range  -°° <  f  < °°.
                                           J
Let the mean value of the random sequence f(t-j), i = 1,2,...  ,m,  generated
above be noted by
                                                                      (D-3)
Then by definition the ensemble mean value of the random variable f(t-j) is
               f = lim f  = lim  2,
                        "   Af+0  j
                                                                      (D-4)
where p(fj) = prob(fj < f <_ f j + Af) and the summation is over all
intervals fj in the infinite range of f values.   From Eqs. (D-l)  and (317),
we obtain
                    f = lim  2,  f-;
                        Af+0  j   J
                                    m(e.
                                      T
Comparing this with Eq. (D-2), we see finally that
             f
                ,+T
f(t) dt = T li
                                 m
                                       n
1 2, f(ti
                                                                      (D-5)
                                                                    on
This result means that we can approximate the integral of any functit
f(t) over any domain by simply multiplying the length of the domain T by
the mean value of the random variable f(t-j) formed by picking points t-j
at random in the domain of integration.

-------
                                                                           217
     The speed of execution of this Monte Carlo integration technique
can be increased greatly by tabulating the exponential  and error func-
tions that appear in the integrands of several  of the integrals of
interest.   That is, rather than use the EXP and ERF FORTRAN routines to
compute the function value each time it is required, vie create a table
of the values of each function initially and look the value up in the
table as needed.  This procedure has been found to reduce computation
time by about one-third.  The subroutines and functions TABLET, EXPO,
and ERROR listed in the CALC program replace the EXP and ERF functions.

-------
                                                      218
              APPENDIX E
              )F rAB, rA,
MULTIJET REACTOR DATA AND TOOR'S  THEORY
DERIVATION OF PAB, i?A, AND ffi  FROM

-------
                                                                            219
                             APPENDIX E
               DERIVATION  OF TAB,  JT,  AND TR FROM
            MULTIJET  REACTOR DATA AND TOOR'S THEORY
     The  purpose  of  this appendix is to give a brief description of the
multijet  reactor  data and how we used the data to calculate the inert
                       /\    <*.       A
correlation  functions (r*n, l\, and rR) used in Chapter III.

     As noted  in  Chapter III, Toor (1962) developed a theory  that relates
the  concentration  statistics of two species undergoing a very rapid irre-
versible  reaction  in a  turbulent mixer to the concentrations  of inert
species in  an  identical mixer.  To test the theory experimentally, Vassilatos
and  Toor  (1965) designed an ideal one-dimensional tubular reactor having a
head made of some  100 small nozzles (Figure E-l).  Reactants  are fed through
alternate jets  to  simulate a cross-sectionally uniform concentration profile.
A modified  version of this reactor was later made and used by Mao (1969).
Extensive measurements  covering a wide range of reaction speeds were made
in these  reactors.   However, the results provided only a partial test of
Toor's theory  because no inert mixing data were taken.  Later measurements
of inert  concentration  statistics in the same reactors by McKelvey (1968)
and  Zakanycz (1971)  have corroborated Toor's theory.
             polymethyl methacrylate
                      epox.y
                                              :1.5
                        4.5
           FIGURE E-l.  SIDE AND END VIEWS OF TUBULAR REACTOR

-------
                                                                            220
     The theory presented by Toor (1962) states that for very rapid reactions
with stoichiometric feed,


                                        *
                                          ,  ,      2
                                        

where x = t is the axial position along the reactor length and XQ is
some reference point where the reactor has reached the state of cross-
sectional homogeneity.  Toor and his coworkers designed the multijet reactor
so that homogeneity is achieved virtually at the inlet of the reactor.  Thus,
we set x~ = 0.  Toor (1969) further derived a relationship between the con-
centration fluctuations of two unpremixed inert species fed into this reactor:
                                                                      (E-2)
                  2

                                                                      (E-3)
Since the species are not premixed,
                        = -      .                   (E-4)
Combining Eqs. (E-l) through (E-4), we obtain

-------
                                                                             221
                     (A B (X)>                   (E-5)
              ii'  —  *    x          /  v x          »•«••* \ /     w         \  ^ /
                      
-------
                                                                222
   o
    •
-------
                                                           223
                     VERY  RRPID RERCTION
                     STGICHIOMETRIC  FEED
                       CD  MRO,  RE=1460
                      —  CURVE FITTING
0.0
   FIGURE E-3.  DECAY OF CONCENTRATION WITH DISTANCE
                 IN MAO'S REACTOR

-------
                                                            224
                     VERY  RRPIO RERCTION
                     STOICHIOMETRIC  FEED
                       O   MRO, RE=2i430
                       —  CURVE FITTING
0.0
   FIGURE E-4.  DECAY OF CONCENTRATION WITH DISTANCE
                 IN MAO'S REACTOR

-------
                                                       225
              APPENDIX F
ESTIMATION OF THE ORDER OF MAGNITUDE OF
    THE TIME SCALE OF DECAY OF A"2

-------
                                                                            226
                             APPENDIX  F
              ESTIMATION  OF THE ORDER OF  MAGNITUDE OF
                  THE TIME SCALE OF  DECAY OF A"2
      When  the  first term on the righthand side of Eq.  (156)  is expressed
in the form (157), we find that, in the absence of any  mean transport,
chemical  decay, or sources of subgrid-scale concentration  fluctuations,
                                                 iV*"*)
an initially present mean square fluctuation field A"   decays according to
                                          A"      ,                 (F-l)
where A is  some  length scale and KA is the eddy diffusivity.   In other
words, the  time  scale of decay of A"? due to turbulent  diffusion is of
the order of
                                                                    (F-2)
We wish to  show  from an estimate of T that A is of the order  of  the size
                                                       -^
of the grid mesh upon which the mean concentration field A is defined.

      For simplicity we consider a one-dimensional problem governed by
the equation


                                   = K^-  + «(z)6(t)               (F-3)
                                       o Z

with initial  and  boundary conditions

                            c(z,0) = 0                              (F-4a)

                        lim c(z,t) = 0    ,                         (F-4b)

-------
                                                                            227
where K is the turbulent diffusivity, assumed to be a constant, and 5 is

the delta function.  The solution of (F-3) and (F-4a,b) is found to be
          c(z,t) = (4TiKt)"1/2 expf- fiTT- )    .
                                 V   KV
Now let
                         Z+AZ
          c(z,t) =
          2Az /
      c(z',t)dz'
                        Z-AZ
Averaging (F-5) in this manner we obtain
                                                           (F-5)
(F-6)
          c(z,t) =
                   4Az
 erf  ^^  -
                                            z-Az
                                                           (F-7)
where
erf(x)  =
/
                                -^
                               e   dx
is the standard error function.  Note that Eq. (F-7) is the solution

of the equation obtained by averaging Eq. (F-3), that is
                                                                     (F-8)
where
                  =  1, -Az <_ z < Az;
                ;
                     0, otherwise.

-------
                                                                            228
Let us assume that Eq. (F-8) is a grid model  representation of Eq.  (F-3)
and that c, as given by Eq.  (F-7), is the model's output.   The subgrid-
scale variations are therefore [from Eqs. (F-5) and (F-7)]
      c"(z,t)   =   -==Lr exp  j-   — .	
              .   V4i:Kt     V  4Kt/   4Az
erf f£y£
           _ era
             er9   .——
                                               •   (F-9)
We can see from this expression that the amplitude of the SSV decreases  as
the discretization interval  Az in the grid model  is made smaller.   Moreover,
for fixed Az the SSV in this example decrease with time.

     To obtain a more concise description of the  relative magnitude of the
subgrid-scale variation c",  we note first that c" is  largest  at  the source
location (i.e., z = 0).  Thus, the maximum-amplitude  of c"  relative to c
can be expressed by
                 -  c" o.t
                 =  „.,
                    c 0,t
                                (F-1Q)
From the expressions for c" and c we obtain
where
a*erf f -TT
                                         -  1
                               (F-ll)
             T*  =
                                                                          (F-12)
                     Az
In words, a* is the approximate half-width of the pollutant cloud,  normalized
by the grid mesh size AZ, at time t;   It is evident from Eq.  (F-ll)  that a*

-------
                                                                           229
is the key parameter which determines the relative  size  of the  SSV  in  this
instance.   On evaluating Eq.  (F-ll)  for several  values  of a*, we  obtain
p = 10 when 0* = 0.1; p = 0.34 when  0* - 1;  and  P  = 0.08 when a*  =  2.  Thus,
the time scale of decay of p  is of the order
                         T
                          P   4K     K                                 ^r~'°

Now since c (0,t),  which  enters  in  the  definition  of p  [in  Eq.  (F-10)J,  also
decays with time,  the  decay rate of c"(0,t)  can  be no faster than  that of p.
Furthermore,  since  c"(0,t)  represents  the  maximum  amplitude of  c^(z,t),  it
                                                                2
is not unreasonable to assume  that  the  time  scale  of decay  of c"  is  com-
parable to that of  c"(0,t).   From all  of these  considerations we conclude
that
and hence that the  length  scale  A,  defined  in  the  beginning of this appendix,
is of the order of  AZ.

-------
                                                                             230
                                 REFERENCES
Batchelor,  G.  K.  (1967),  "Small Scale Variation of Convected Quantities Like
     Temperature  in  a Turbulent Fuel, Part I," J. of Fluid t-lech. , Vol. 5, p. 113.

	 (1950), "Application of the Similarity Theory of Turbulence to Atmos-
     pheric Diffusion," Quart. J. of the Royal Meteor. Soc., Vol. 76, pp. 133-146.

Bendat,  J.  S.,  and A. G.  Piersol (1968), "Measurement and Analysis of Random Data,"
     John  Wiley and  Sons, New York, New York.

Briggs,  G.  A.  (1968), in  Meteorology and Atomic Energy, David Slade, ed.

Chen, W. H.,  and  J.  H. Seinfeld (1974), "Estimation of Spatially Varying Parameters
     in  Partial Differential Equations," Int. J. Control, Vol. 15, pp. 487-495.

Deardorff,  J.  (1972), "Numerical Investigation of Neutral and Unstable Planetary
     Boundary  Layers," J. Atmos. Sci., Vol. 29, pp. 91-115.

           (1970), "A Three-Dimensional Numerical Investigation of the Idealized
     Planetary  Boundary Layer," Geophys. F1uid Dyn. , Vol. 1, pp. 377-410.

Donaldson,  C. duP., et al. (1973), "Atmospheric Turbulence and the Dispersal of
     Atmospheric  Pollutants," Report EPA-R4-73-016a-c, Environmental Protection
     Agency,  Research Triangle Park, North Carolina.

Donaldson,  C. duP.  (1969), J. of AIAA, Vol. 7, pp. 271-278.

Golder,  D.  (1972),  "Relations Among Stability Parameters in the Surface Layer,"
     Boundary Layer Meteorology, Vol. 3, pp. 47-58.

Hanjalic,  K., and B. E. Launder (1972), "A Reynolds Stress Model of Turbulence
     and Its  Application  to Thin Shear Flows," J.  Fluid Mech. , Vol. 52, pp. 609-638

Hilst, G.  R., et  al. (1973), "A Coupled Two-Dimensional Diffusion and Chemistry
     Model  for  Turbulent  and Inhomogeneously Mixed  Reaction Systems," Report
     RA-73-016C,  Environmental Protection Agency,  Research Triangle Park, North
     Carol ina.

Hinze, J.  0.  (1959), Turbulence, McGraw-Hill Book  Company, New York, New York.

Kreith,  F.  (1966),  Principles of Heat Transfer, Second Edition (International
     Textbook Company, Scranton, Pennsylvania).

Lamb, R. G. (1974),  "Lagrangian Analysis of a Passive Scalar in Turbulent Fluid,
     Part 1:  General Theory," J.  Fluid Mech., submitted for publication.

           (1971),  "Numerical Modeling of  Urban Air Pollution" PhD. dissertation.
     Department of  Meteorology, University of California, Los Angeles, California.

-------
Lamb, R. G. (1973), "Note on the Application  of K-Theory  to  Diffusion  Problems
     Involving Nonlinear Chemical  Reactions,"  Atmos.  Env. , Vol.  7,  pp.  257-263.

Lamb, R. G.,  and J. H.  Seinfeld  (1973),  "Mathematical  Modeling  of  Urban  Air
     Pollution — General  Theory,"  Environ.  Sci.  Techno]^,  Vol.  7, pp. 253-261.

Lamb, R. G.,  and M. Neiburger (1971),  "An  Interim  Version of a  Generalized
     Urban  Air Pollution  Model,"  Atmos.  Environ. ,  Vol.  5, pp.  239-264.

Lettau,  H.  H.  (1959),  "Wind  Profile, Surface Stress,  and  Geostrophic Drag Coef-
     ficients  in the Atmospheric  Surface Layer," Advances in  Geophysics, Vol. 6,
     pp. 241-256.

Lions,  J.  L.  (1971), Optimal  Contro1_of  Systems  Governed  by  Partial Differential
     Equations (Sprinqer-Verlaq,  New York, New  York)  396  pages.

Lumley,  J.  L.  (1970),  "Toward a Turbulent  Constitutive  Relation,"  J. Fluid Mech.,
     Vol.  41,  pp.  413-434.

Lumley,  J.  L., and B.  Khajeh-Nouri  (1973), "Computational Modeling  of  Turbulent
     Transport," Second  IUTAM-IUGG  Symposium on  Turbulent Diffusion and  Environ-
     mental Pollution, Charlottesville,  Virginia.

Mao, K.  (1969),'  "Chemical Reactions with  Turbulent Mixing,"  Ph.D.  thesis,
     Carnegie-Mellon University,  Pittsburg, Pennsylvania  (173 pages).

McElroy, J. L.  (1969), "A Comparative Study of  Urban  and  Rural  Dispersion,"
     J.  Appl.  Meteor., Vol.  8, pp.  19-31.

McKelvey,  K.  N.  (1968),  "Turbulent  Mixing  with  Chemical Reactions," Ph.D. thesis,
     Ohio  State  University,  Columbus, Ohio (123  pages).


Monin,  A.^S.,  and A. M.  Yaglom (1971), Statistical Fluid  Mechanics, MIT  Press,
     Cambridge,  Massachusetts.

O'Brien, E. E.  (1968),  "Closure  for Stochastically Distributed  Second-Order
     Reactants," Physics  of  Fluids, Vol. 11, p.  1883.

Orszag,  S., and M. Israeli  (1974),  "Numerical  Simulation  of  Viscous Incompressible
     Flows,"  Ann.  Rev,  of Fluid  Mech., Vol. 6,  p.  281.

Ragland, K. W.  (1973),  "Multiple  Box Model for Dispersion of Air Pollutants  from
     Area Sources," Atmos.  Environ. , Vol.  7, pp. 1017-1032.

Shir, C. C.,  and L. J.  Shieh (1973), A Generalized Urban  Air  Pollution Model and
     Its Application to  the  Study of S02 Distributions  in the St. Louis  Metropolitan
     Area," Report RJ  1227,  IBM,  San Jose, California.

Shu, W.  R.  (1975), Ph.D.  dissertation, in  progress.

-------
                                                                            232
 Toor,  H.  L.  (1969)  Ind. Eng. Chem. Fundam., Vol. 8, pp. 655-659.
 	       (1962),  "Mass Transfer in Dilute Turbulent and Nonturbulent
     Systems  with Rapid Irreversible Reactions and Equal Diffusivities,"
     A^LCh.E. Journal, Vol. 8', pp. 70-78.

Turner, D.  B.  (1969),  "Workbook  of  Atmospheric Dispersion Estimates,"  Publi-
     cation  AP-26, Environmental Protection Agency, Research Triangle  Park,
     North  Carolina.

Vassilatos,  G., and H. L. Toor (1965), "Second-Order Chemical  Reaction  in  a
     Nonhomogeneous Turbulent Fluid," A.I.Ch.E. Journal, Vol.  11,  pp.  666-673.

Zakanycz, S.  (1971), "Turbulence and the Mixing of Binary Gases,"  Ph.D. thesis,
     Ohio State University, Columbus, Ohio (172 pages).

-------
                                                                                   233
                                   TECHNICAL REPORT DATA
                            (f'lcasc read Inunii'liun^ on l/ic /i i im1 be tore coinpli InifJ
    .'.-'iOO/.'4-76rQ.16_c_
    . ,  ,-.'. ^ bun 1 II LE
                  CONTINUED  RESEARCH IN MESOSCALE  AIR
   il.I.KIlON  SIMULATION MODELING.   VOLUME III.   Modeling
  !' v'i i-roscnle  Phenomena
             5. REPORT DATE
             _ May T9 7 6    	
             6. PERFORMING ORGANIZATION CODE
 R. G. Lamb
             8. PERFORMING ORGANIZATION REPORT NO.

                EF75-25
  >•(. ill OH VI.NT, ORG \NIZATION NAME AND ADDRESS
 SYSTEMS APPLICATIONS, INC.
 950  NORTHGATE DRIVE
 SAN  RAFAEL,  CALIFORNIA   94903
 I/ ;,»0
-------