EPA-R4-73-030e



July 1973
ENVIRONMENTAL  MONITORING SERIES




-------
                                        EPA-R4-73-030e
URBAN AIR  SHED  PHOTOCHEMICAL

     SIMULATION  MODEL  STUDY
 VOLUME  I  -  DEVELOPMENT  AND  EVALUATION
     Appendix  D  -  Numerical  Integration
           off  Continuity Equations

                      by
                  S. D. Reynolds
              Systems Applications, Inc.
               9418 Wilshire Boulevard
             Beverly Hills, California 90212
               Contract No. 68-02-0339
             Program Element No. 1A1009
          EPA Project Officer: Herbert Viebrock

               Meteorology Laboratory
         National Environmental Research Center
       Research Triangle Park, North Carolina 27711
                  Prepared for

        OFFICE OF RESEARCH AND DEVELOPMENT
        U.S. ENVIRONMENTAL PROTECTION AGENCY
              WASHINGTON, D.C. 20460

                   July 1973

-------
This report has been reviewed by the Environmental Protection Agency and




approved for publication.  Approval does not signify that the contents




necessarily reflect the views and policies of the Agency, nor does




mention of trade names or  commercial products constitute endorsement




or recommendation for use.
                                 11

-------
                                        EPA-R4-73-030e
URBAN AIR  SHED  PHOTOCHEMICAL

     SIMULATION  MODEL  STUDY
 VOLUME  I  -  DEVELOPMENT  AND  EVALUATION
     Appendix  D  •  Numerical  Integration
           of  Continuity Equations

                      by
                 S. D. Reynolds
              Systems Applications, Inc.
               9418 Wilshire Boulevard
             Beverly Hills. California 90212
               Contract No. 68-02-0339
             Program Element No. 1A1009
          EPA Project Officer: Herbert Viebrock

               Meteorology Laboratory
         National Environmental Research Center
       Research Triangle Park, North Carolina 27711
                  Prepared for

        OFFICE OF RESEARCH AND DEVELOPMENT
        U.S. ENVIRONMENTAL PROTECTION AGENCY
              WASHINGTON, D.C. 20460

                   July 1973

-------
This report has been reviewed by the Environmental Protection Agency and




approved for publication.  Approval does not signify that the contents




necessarily reflect the views and policies of the Agency, nor does




mention of trade names or commercial products constitute endorsement




or recommendation for use.
                                 11

-------
                          TABLE OF CONTENTS


                                                                  Page

INTRODUCTION	D-l

I.    THE METHOD OF FRACTIONAL STEPS APPLIED TO THE INTEGRATION
      OF THE MASS CONSERVATION EQUATIONS	D-6

      A.  Transformation of the Governing Equations	D-6
      B.  The  Computational Grid	D-ll
      C.  Nomenclature	D-ll
      D.  Finite-Difference Equations	D-13

II.   METHODS FOR SOLVING THE LINEAR SYSTEMS OF EQUATIONS.  ...   D-33

      A.  The Algorithm for the Tridiagonal System	D-34
      B.  The Algorithm for the Block-Tridiagonal System ....   D-35

III.  DISCUSSION OF THE NUMERICAL METHOD	D-38

-------
.INTRODUCTION

     In this Appendix we present a  comprehensive  discussion  of  the
methodology used  to  obtain  the numerical  solution of  the  governing
airshed equations.   During  this second phase  of model development
and  evaluation, we have carried out further studies of the numerical
scheme described  in  Appendix  D of Roth et al.  (1971)   and implemented
several important changes.  We shall discuss  the  motivation  for
making these changes and include a complete  description  of  the finite-
difference equations that were finally adopted.   All  simulation results
included  in the Final Report  were obtained using  the  numerical  methods
described in Sections I and II.

     We begin by  writing the  governing airshed equations, which express
the  conservation  of  mass of each chemical species in  a turbulent fluid
in which  chemical reactions occur.   Employing the so-called K-theory, these
equations can be  written as:
                                 -
      3t     3x      3y      3z      3x
                                   I =  1,  2,  ...,  p
 for
                    h(x,y)  <_ z £ H(x,y,t)

-------
where
          x,y = horizontal coordinates
            z = vertical coordinate
        u,v,w = the three deterministic components of the wind vector
           c^ = mean concentration of species A
        K^,IC. = horizontal and vertical turbulent diffusivities
           S^ = rate of emission  of species  £  from elevated sources
           Rj, = rate of production of species  &  through chemical
                reaction
    fXpfyq »vvr = west, east, south, and north model boundaries
       h(x,y) = terrain elevation
     H(x,y,t) = inversion base elevation
     Equation  (1) must be solved subject to initial and boundary
conditions.  The initial condition is
                       c£(x,y,z,to) = fA(x,y,z)                     (2)

where fj;, is the initial concentration distribution of species  &  .   The
boundary conditions that apply at the ground and inversion base  are:
 (1)  z = h(x,y)
where
n,  is the unit vector normal to the terrain directed into the atmosphere,
and Qj is the mass flux of species  £  at the surface.*
*  We indicate vector and tensor quantities through the use of single
   and double underbars, respectively.
                               D-2

-------
(2)   z = H(x,y,t)





           ~^-l*c£) 'U  = Z fy (x,y,z,t) -n_H          if
                       = 0                          if V-nH> 0
                                                                      (4)
where V is the advective velocity of pollutants relative to the moving

inversion base,




                                           9H
                       V =  ui  +  vj + (w - -T- )  k
                       —    —     —        dt   —






n_   is the unit vector normal  to the surface defined by the inversion

base and directed out of the modeling region,  and g^ is the mean concen-

tration of species  i  aloft  (just above the inversion base).



      The conditions that apply  along the horizontal boundaries are:
 (3)   x = x  or
            (uCjl  -  KJJ-;^)  = ug£                      if U-n _< 0


                                                                      (5)


                 -  K ——  = 0                       if U-n > 0
 (4)   y = ys  or  yN
            (vc*  '  ^^  = vg                      if
                                                    if O.n > 0
                                 D-3

-------
where U_ = ui_ + vj_, n_  is the outwardly directed unit vector normal to
the horizontal boundary, and g, is the mean concentration of species
H  just outside the airshed boundary.

      These equations, which are a repetition of Equations  ( 5) -  (9)
in the Final Report, are coupled since they share a common argument  in
the chemical reaction terms R^ (c.. ,c_ ,... ,c ).  In the present version
of the model, we simulate the dynamic behavior of six chemical species
—reactive and unreactive hydrocarbon, NO, N0_, O,, and CO.  The extent
to which "unreactive" hydrocarbon and CO participate in chemical reactions
is generally small, and hence we have set R., equal to zero in Equation  (1)
for these two species.  Consequently, we must solve a system of four coupled,
nonlinear partial differential equations involving reactive hydrocarbon,
NO, NO  , andO_, and two sets of linear partial differential equations
involving unreactive hydrocarbon and CO.

      In Appendix D of Roth et al. (1971), we formulated a numerical
scheme based on the method of fractional steps for the solution of
equations similar to (1) - (6), the  only difference being that horizon-
tal diffusion was ignored.  The method of fractional steps, discussed
by Yanenko (1971), was especially developed to treat multidimensional
problems arising in mathematical physics.  The characteristic feature
of the method is that the solution of a multidimensional problem is
replaced by the successive solution  of several problems having fewer
dimensions.  For example, a four-dimensional equation in (x,y,z,t)
space may be split into three two-dimensional equations, in  (x,t),
(y»t), and (z,t)  space.

      One of the primary considerations regarding the use of finite-
difference methods for the solution  of Equation (1) is the manner in
which the advection terms are to be  approximated.  In Appendix D of
Roth et al. (1971), we discussed this problem and proposed to study
two high-order techniques—a second-order scheme given by Crowley
(1968), and a fourth-order scheme devised by Fromm (1969) .  After fur-
ther examination of these two schemes, we found that both had the unfor-
tunate property of producing negative concentrations in areas on the
grid where large concentration gradients occurred.  Thus, to avoid nega-
tive concentrations and yet maintain the accuracy of a high-order tech-
nique, we implemented another second-order scheme, described by Price
et al.  (1966), which we found to be  capable of handling gradients with-
out producing negative concentrations.  We discuss the application of
this method in more detail in Section III.

      In the next section we describe the algorithms used to obtain  the
numerical solution of Equations (1)  - (6).  First, we perform a change of
                                D-4

-------
variables that transforms the region between the ground and the inver-
sion base from an irregular to a rectangular region.  In so doing, we
write the transformed equations in a way that makes it possible to use
a conservative finite-difference representation.  Second, we formulate
all finite-difference schemes in a conservative manner.  And finally,
we employ a Newton iterative procedure to solve the nonlinear equations
which result from the implicit finite-difference treatment of the (z,t)
fragment.

      In Section II, we include the recursive formulas used to solve
the systems of linear equations which result from the application of
the implicit finite-difference methods described in Section I.

      The last section of this Appendix is devoted to a discussion of
the numerical method.  We explain the difficulties encountered in the
application of advective difference approximations, and the manner in
which they were overcome.  To evaluate the numerical treatment of the
nonlinear chemical reaction terms, we have conducted two tests of the
numerical method.  The results of these tests are presented and dis-
cussed in Section III.
                                 D-5

-------
I.    THE METHOD OF FRACTIONAL STEPS APPLIED TO THE INTEGRATION  OF THE
      MASS CONSERVATION EQUATIONS

      In this section we discuss in detail the numerical solution of
Equations (1) , six partial differential equations in three spatial di-
mensions and time.  The problem is complicated by the fact  that  four
of the equations are coupled and nonlinear due to the appearance of
the chemical reaction term R£.  To apply the method of  fractional
steps, we fragment the original four-dimensional equation  (x,y,z,t)
into three two-dimensional partial differential equations in  (x,t) ,
(y,t), and (z,t) respectively.  We include the chemical  reaction  and
elevated source terms, R^ and S^, in the (z,t) fragment.  The solution
of the (x,t) and  (y,t) fragments is explicit, while the solution of
the (z,t) fragment is implicit.  To compute an integration  time  step,
we successively solve the equations representative of the  (x,t),  (y,t),
and (z,t) fragments.  Before presenting the details of  the  finite-
difference equations, however, we first discuss certain aspects  of the
equations and their solution.

A.    Transformation of the Governing Equations

      The modeling region, as defined, is bounded in the vertical by
the ground and the base of the elevated inversion.  We  designate the
elevation difference between these two boundaries as the mixing  depth.
One of the characteristic features of the mixing depth  is that it
varies with location  (i.e., with x and y at a given time t).  This
variation is due in part to variations in terrain elevation,  and in
part to the general increase in the elevation of the inversion base
with distance from the coast.  We also note that the upper  boundary
changes position with time as the inversion base moves  up,  then  down,
during the course of a day.  To implement the finite-difference  scheme
on a "fixed" rectangular grid, we perform the following change of
variables:

                                T = t                              (7a)

                                5 = x                              (7b)

                                n = y                              (7c)

                                       z-h(x,y)                    ._  .
                                p ~ H(x,y,t)-h(x,y)                v a;
                                 D-6

-------
      Under this change  of variables,  Equation. (1)  becomes
         (AHC,) + £ (uAHC^) + J- (vAHc,)
                     KHAH
                       3p v  A
                           3p
                  3n
                          3c
                       3n
                  3p
                                3AH
                         fe+  m\H*
                         \dr\  p 9n /  3 n
                        S^AH
                              l,2,...,p       (8)
where
AH -




 W
                                                                (9)
    initial and boundary conditions become:
                                                               (10)
                              D-7

-------
(1)    P « 0
      Q4(5,n,T)
AH
                         AH  3n     AH
                                          9n
                                                                         (11)
(2)   p = 1
        AH
                                  AH
                                                    AH
                                                                           if W < 0
              0 =
                           AH / 3p
                                         3n
                                                                          (12)




                                                                          if W > 0
 (3)   5 -  Cg   or  5W
UCA -
                      AH
         r!!*   i
         l 35  " AH
                  35    AH  35     35 I 3p
                                                             if U-n < 0
                                                            if U-n > 0
                                                                          (13)
                                  D-8

-------
(4)   n = r)  or n
           S     N
          !!£   -i-fe,   l£H\!!i
           3n  " AH I 3n    P 3n  J  3p  l~
      In the application of Equations (8) - (14) to the computation of
pollutant concentrations in the Los Angeles Basin, several terms in the
above equations are generally small with respect to other terms  and
therefore can be neglected.  Changes in the terrain and inversion base
elevation with location are gradual, indicating that the derivatives
3h/3£, 3h/3n,  3AH/35 and 3AH/3n  are considerably smaller than one.
With the assumption that terms in Equation (8) containing these deri-
vatives can be neglected,* Equation  (8) becomes
                   3  /  A   8cA   3  /     3cA   3 /Kv
                   K ^VH IT/ * 3n \VH ~3Ty + "37 \te
                                     £=i,2,...,p
We note that, whereas the assumption is valid that the derivatives  3h/3C»
3h/3n, 3AH/85 and 3AH/3n are small in the-flat portions of the Los Angeles
airshed, it does not hold in the mountainous areas, such as the Palos
Verdes Hills and the Santa Monica Mountains.  However, since the terms
omitted in Equation  (15) involve horizontal diffusion, which is generally
less important than horizontal advection, we do not expect to incur ex-
ceptionally large errors in using Equation  (15).  In the application to a
   To properly account for convergence and divergence effects,  W is  still
   given by Equation (9),   See Equations (71)  - (76) for a description
   of the evaluation of  W ,
                                 D-9

-------
different airshed  in which topographical effects are of prime  importance,
it would be desirable  to use Equations (8) - (14) as the basis of the
model.

      Under the  assumptions invoked in obtaining Equation  (15) ,  the ini-
tial and boundary  conditions become:
                                                                     (16)
 (3)
 (1)   P = 0


                               K., 9c

(2)   p = 1


                ..      Kv
                                                 ifw  0
                  AH 8p
                   UC n - K_. -r=- = ug o            if U-n  <  0
                     x-iiac     ^ *•               	
                                                                     (19)

                   K.. -TT- =0                    if U-n  >  0
                                  D-10

-------
(4)   n =  n_ or  n.
            S      N
                    - K.. -r—  =  vgp               if U'n < 0
                       H  of)        *                  — — ~~

                                                                    (20)

                    ——  =  0                      if U-n > 0
                     3n                                	
      Thus we shall use Equations (15) - (20) as the basis for the
simulation model of photochemical air pollution in the Los Angeles
airshed.  We note that these equations can be solved on a fixed rec-
tangular grid with spacing A£, ATI, and Ap.

B.    The Computational Grid

      Originally, a 25 x 25 x 5 grid network was adopted for the Los
Angeles simulation study, where  A?  = An  = 2 miles, and Ap = 0.2
( £ and n  can be scaled by the length of the airshed, but this only
introduces additional arithmatic operations into the finite-difference
equations).  Upon examining the 25 x 25 horizontal grid layout, we
found that 198 of the 625 grid squares were over the Pacific Ocean and
the San Gabriel Mountains—two areas having no emissions sources.  Con-
sequently, we eliminated these grid squares from the original grid to
form the region over which the actual computations were performed.  This
amounted to a savings of approximately 30% in the required computing
time.  The grid is illustrated in Figure 1.  Thus the computational
grid consists of 2135 grid cells (427 horizontal grid squares x 5 verti-
cal levels).

C.    Nomenclature

      The exposition of finite-difference schemes generally requires the
use of subscripts and superscripts to denote spatial and temporal aspects
of a particular variable.  This application is no exception, and, in
fact, is somewhat more complicated than most.  We shall denote the con-
centration distribution on the grid by:

                              m n
                                 D-ll

-------
                                                                               25-
•  I  i
                                                                            25-
                        Figure  1.   The Modeling Region
                                   D-12

-------
where      £    is  an  index denoting chemical species, 1  < £  < 6

       i,j,k    are indices designating the grid cell location
                in  the £-i":-p coordinate system.

                    *W - 1 - XE
                    Js .5 3 _5 JN
                     1 _< k j< K

                 where Iw,  IE, J0, and  JN indicate  the column or row
                 number of  cells adjacent to  the  west,  east,  south,  or
                 north borders, and  K is the  number of vertical  levels.

           n    is  an  index in time, n  ^_ 0

           m    refers  to the iteration nuinber,  m >_ 0

We show subscripts  and superscripts  on  other  variables in an analogous
manner.

D.    Finite-Difference Equations

      The method of fractional steps is applied  to the solution  of
Equations (15) - (20)  as follows.  First, Equation (15)  is split
into three two-dimensional  equations in (£,T) ,  (n^)  and (P,T).   The
three equations are given by:
Step I:    -^  (AHCjl) + -^  ^x^,  -  gg i •>„«« -35-1                  (21)



Step II:   —  (AHCp) + —  (vAHcp)  =  T- (K AH -r—)                  (22)
           9x      *    ori       *    9ri V  H   3r)  I


Step III:  A  (AHcj) + -j~  (WCo)    =  -j- (7^ -r-^ )+ R.AH + S .AH      (23)
           ^—      -r    dp     *      op \ fin  dp  /    J6      X/
           O I                           \       r



The calculation of   0C.  .     given values of   C. .    requires three
 .                    x*i/"i/jc                    x/ i / n »jc
steps:

      Step I:    Integration  in  the   5   direction
      Step II:   Integration  in  the   n   direction
      Step III:  Integration  in  the   p   direction

We now present the actual  finite-difference expressions employed in the
three-step integration procedure.
                                 D-13

-------
STEP I:  Integration in the  £  Direction


      We advance the concentration distribution at the n   time level,
  n                                               *
 c.    , to the first intermediate time level, 0C.  . .  , by integrating
* ^ »3 »k                                        x.  i , j ,K

Equation (21) over the time interval AT.  The explicit, seconds-order


finite difference expressions adopted for the spatial derivatives are

based on the results of Price et al, (1966) .
       r*     =  r"      +    AT
                         *                                             (24)
where

                          n
                                        rn
                                        C
      \
-2,j,kl
                             2         £i-l,j,k   £ i-2,j,k           (25)
                            n         If*          f^    \
                            i-i,j,k   \fi-a, j,Jc   Jri,j,k/    if  aj^^^  >.
or
                          n

      £l   Pn          = ^i:

                                          -l,j,k-  £Ci,j,k)   ifaU,jfk<0
                                 D-14

-------
                       u.  i. .  AH. i   . AT
       an  ,         =   V'*'J'k	i-i'J	                               (27)
        i-i,j,k                A?
  Y? ,  ,  .     =     *•-*•!•*•	                             (28)
   i-i»3»k                A?2

     In Equation (24) we have made the following assumption with regard
to the evaluation of AH at the first intermediate  time  level:
             *       n
           AH.  . = AH.  .

Since the effects of temporal mixing depth  changes on pollutant concen-
trations are best treated in Step III, we defer movement of the inver-
sion to that point in the computational  scheme.
     When applying Equation  (24) in grid cells adjacent to the western
or eastern borders  (£„ or E._), we employ the  expressions for the
                     W     E
boundary conditions given in Equation  (19).   Thus, grid cells adjacent
to the western border  (£_7) are treated as follows:
                        W

 t I..i jfk        i. I .i.k + . _. n    UFT -4-.-i.v-  ~ oFr  xJ. -i v)     *29)
                                                                   \
                                                              *'j'V
where
                                    n                        n
                                                         if
                                                                        (30)
     AT  ji            V*'3^ /   n
                                                n
                                            "
                                                         if a"  ,  .    <  0
                                 D-15

-------
and
            n
                                      n
                                                  n
                                                                       (32)
      rn
,k " £CI
                                                            \
                                                     w+l,j,kJ

                         n
      AT
                                                                        (33)
            w
                                         1£
                                     n
                                                    n
 In Equation  (32) we have replaced  fcciw-l,j,k by  !fllv-h,i,K
 the value of the concentration outside the modeling region.
 The treatment at the eastern border (C )  is analogous to that given
                                       £
 for the western border above.  It should be noted that in moving
 from grid cell to grid cell in the E, direction, we need only compute
                  A? I i
                                 D-16

-------
for grid cell  (i,j,k),  since
is available  from  the  calculations performed in grid cell  (i-l,j,k).
STEP II:  Integration in the   n  Direction


      In this step, we  advance from the first level intermediate  con-
                *
centrations  .C.  . .  to the  second level intermediate concentrations
  **         *•  1t'3tK-
.C.  . ,   by the finite-difference analog of Equation  (22). These  equations
*• i f D rK
are the  n  direction counterparts of Equations  (24)  -  (33).



            **            *          AT   /  *           *       \
           f        —    f      4-         IF        —  T?        I       ttA\
          0 ,' -5 V   ~   ff-i-ilr       n   ln
             '3'           f11'     AnAH  .
                                      1 »D


where


                        Rn
       AT   *           pi,i-i,k /„ _*
          ,F
       An ri,j-i,k        2
                                                                        (35)
                                         if
                                 D-17

-------
£L  F*          =    I*j-*A
An «. i,j-i,k            2
                                        *         *      \
                                                 r       i
                                        1,:)^   i i,j+l,k)
                                                                       (36)
                                         if 6? .  ,  .   < 0
and
      _n                 v.  .  _,    AH.  .  i  AT

          -               1>J        1>]" 	
                        K" .        AH1}  .  ,  AT
      ,n                 H.  .  ,      i/D-5!

       i.j-i.k        -   lf3~*'K	:—                         (38)
Using the expressions for the boundary  conditions as given in Equation

(20), grid cells adjacent to the  southern  border (nc)  are treated as
                                                    s

follows:
      *Ci,Js,k        '   £Ci,Js,k
                                            i,J                          '

                                                                       (39)
where
           *          _ /R0        \   "
           i,j -J5,k   ~ I pi,j  -J5,k ) £gi,j -is,
              S         \   S    /       S
Al  F
An ri,j_->»,k     ri,j_-j5,k; £yi,j.-i5,k       if r; T.,  ,. > o (40
                                                       s
                                  D-18

-------
      An «,i,J ->5,
              s
 1"V15'k  /    *          *


	2	   3ACi J  k " *Ci J +1 1
    2      ^  X, 1,J ,k   * l,Jg+l,J
                                            if
      Aj_   *               i,Js+*i,k                 n       .


      An ri,Js+»5,k   =      2       ( \Ci,Jg/k " Hgi,Js->5,kJ
                                                   *        \

                                                   i,J +l,k )
                                                      s     /
                                            if  •.j
                                                    S
                          e"
      A    j.               if*J"^";J/K.   j.
      AT   *                  Q      /   *
      a L  p            _       s      /3 p         -  C

      An ^ i,J +h,k    "      2       I  £ ifJ +l,k    1 i
                             i,Js+2,k)
                             n
                                        *           *
  6^ T ^ JlC*i T  V * "C^  -  -- '-I           (42)
   lf s ^'H   ' s'k
                            •i^^i.,.* - »ci.Vi.xj
                                           if  B   , _,_, .  <  0
                                                i,J +h,k
                                                   S
The treatment of  grid cells adjacent to the northern border (nN)  is

analogous to that given for the southern border above.
                                 D-19

-------
    STEP III:   Integration in the  p Direction
                                         **
         Having established values of  .C. .     from Step II, we now complete
                                       *• i, j ,k
    the integration procedure by solving a finite-difference approximation

    of Equation (23).   These equations have been formulated in an implicit

    (Crank-Nicolson) manner to avert stability problems that might arise in

    the treatment of the diffusion term when the grid spacing is small due

    to a shallow mixing depth.  In addition,  the implicit formulation is also

    more suitable for treating the nonlinear chemical reaction terms.
n+1      _   **     ""i,J        AT
         '
    /  n+1

n+1 V i'j'k">
                                                              p
                                                               **
                                                          _n+l   A
                                                          Ri.n.kA"i.i        (44)

                        **      n  \     AT   /  n+1     n+1
                              Aun  ^ .L —£1	/ SnT^   AH. 7
                                               * i»j »k   i»j
where

                     An+1
      |L  ^+1     =  i^k"lg (£C^+^ k 1 + £C^+^ J                    (45)
                                                           =  2,3
      ^1 F**      =  i'P'k"ls [  c**      +  nC**    ]                    (46)
      Ap £*ifjfk-»i       2    U i,j,k-l    £  i,j,k^
                     +
                        n
                                    **           **
                                  D-20

-------
                                             AT
                                        Ap
                                                                       (47)
                         n
                                              AT
                                                                        (48)
The boundary conditions  at  the ground are:
(1)   P = 0
£L  Fn+1    -  AT
Ap £rifjf«s     Ap
                                                                        (49)
                           **
                                          n
                                  -
                      Ap  £  i,j,>5     Ap
                                                                        (50)
where k = h is equivalent  to p  ~ 0
 (2)   p = 1
AT  _n+l
                               i   n+1

                                  9
                                                                        (51)
AT

Ap
                               ' f
                               ir
                                                                  0     (52)
                                 D-21

-------
                                                                       (53)
where  k = K+'s  is equivalent to  p = 1.


      At this point we need to discuss the following three topics:
      1)   Step III applied to inert species
      2)   Step III applied to chemically reactive species
      3)   The calculation of W.


The Treatment of Inert Species

      For those species which are treated as being chemically inert,

we set  0R.  .   = 0.  The difference equations for Step III reduce to
          1/D'
a system of linear equations involving the values of  0C.  .   .  We can
                                                      *• 1 1 D rK
write the linear system for species  H  as follows :
                                 D-22

-------
                     Vl Vl Vl




" cn+
£ 1'
£ i,
*<
»<
1
1
1
1




~ ~
yl
V2
Vl
yK :
                                                                       (55)
where the elements of the matrix and the column vector on the right-hand
side of the equation can easily be determined using the finite-difference
expressions given in Step III.  We describe the algorithm used to solve
this tridiagonal  system of linear equations in Section II.

      The Treatment of Chemically Reactive Species
      The computation of vertical transport coupled with chemical reaction
is the most complicated step in the numerical procedure.  We will thus dis-
cuss this aspect of Step III in more detail.  We recall that reactive
hydrocarbon (HC), NO, N00, and 0, are all coupled through the chemical
                        4b       <3
reaction terms. Therefore, we are confronted with solving a system of nonlinear
algebraic equations in 4K unknowns, where K is the number of vertical levels.
In the present simulation effort, we employ five vertical levels, and so
we must solve twenty nonlinear equations in twenty unknowns, where
the unknowns are the values of ,C? .   , 1 ^ k _f 5, for HC,NO,NO2» and O,.
                               * i i j »*
                                  D-23

-------
     To solve the system of nonlinear equations, we use a Newton iterative


procedure.  We begin by introducing the iterates ^C.  .  .  Using Equation  (23),

                                    o n+1          i»3fk
the initial values of the iterates, 0C. .    , are given by the following
                                    * i/ 3 »k

slightly modified Taylor series expansion
      °cn+1
      £ i,
                        Au                /                     \

                  **      i/J      AT    /   **          **     \

                  i,j,k AHn+l    Ap AHn+l UFif j,k->j " £ i,j,k+»J
                   ?  . ,, AH?  .)
                   i,D,k   i,]^
                                                           AH
                                                             n+1
                                                                       (56)
The difference expressions for the terms in Equation  (56) are given in


Equations  (45) -  (54).   Higher iterates are obtained by setting
              m+1 n+1

                £ i,j,k
£ i,j,k
                                            m6Cn+1
                                            £  i,j,k
                                        (57)
Next, we substitute the right-hand side of Equation (57) into Equations


(44) - (54) in place of  iCjk  and neglect all quadratic  terms
involving
                    .   We can write the resulting system of linear
equations in the following form:
                                   K-l
                             \   \




I"6X1
m.
m.
*2
{K-1
m
A Y

—


m
pl
m
P2
K-l
K
                                                                       (58)
                                   D-24

-------
where A, ,  B  and E.  are 4x4 matrices, and  6X.  and  e  are column
vectors.   The coefficient matrix is block-tridiagonal, and systems of
equations in this form can be solved efficiently.  The chemical reac-
tion rate expressions are given in Appendix B of this report.  We de-
fine below the elements of each matrix and column  vector.
                                             k =  2,3,...,K           (59)
where
                                                                    (60)
and  I  is the 4x4 identity matrix

      \ = ek I                              k  =  1,2,...,K-1        (61)

where
      \=  —TTT  v  ~rjr"   -   ^'ik^y                      (62)
       *•    ~ i .JIT-J.  \     ^           1. I.JVT"! /
     \
                                D-25

-------
where
      1+\
               n+1     \
              Pi,j,k+'s 1
                   ,n+l
                                                        k = 1
      1 +
          2AH.  .
,n+l
                               ,n+l
        + 2AHX \
                     ,n+l
                                                        K  ^  
-------
where
                   k!2°H + k!3°3
      d!2 ' d21 = d!4 = d41
      d22 - k3°3 + k6N°2 + V°2 +  k!4R°2
                 - k
      d31 = k!3°3
      d32 * k3°3
      d34 * k4°3
      d42 = k6N°2 ' k3°3 *  k9H°2  '  k14R°2



      d43 = k4N02-k3NO




      d44 = kl + k4°3 + k5N°3 + k6N° + k!0H°2 + k!5R°2
and
       °3
      NO, =  C    .
        2   4 i,],k
                                 D-27

-------
      The expressions for the steady-state species O, OH,  HO  ,  R02/


and N0_ appearing in the matrices   D,  (see Equation  (65))  are  evaluated
                        -Li
using the values of   C.  .   .  See Appendix B  for the  expressions re-
                     *• i / 3 /k

lating the concentrations of the steady-state species to  the  concentrations


of HC, NO, NO2, and 0.,.  The remaining terms in Equation  (58) are defined by:
              m.
                                m
                                  . ^
                                  6c.  .
                                m  n+1

                                3  i,j.k
                                                                        (66)
                                                                         (67)
 where
    ID—
and  P.
                                                m—
                                                m—
                                                m—
                                                m—

                                                4%
                                                                         (68)
                                   D-28

-------
and
           **
£Pk   =  £Ci,j,k
               AT
                                               **
                                                           **
                                                      ~  £Fk,j
                                        » J
      \
,j,k+»5y
                                                   (69)
                 2AHn+1
+ 1
                                                AH
                                                  n+1
     m—
     £pk
                                                                        (70)
                   AT
                 2AHn+1
We note that each iteration involving the solution of Equation  (58)
requires only the additional evaluation  of  mD,  and  "p   before the
system can be solved,  since all other terms are invariant  from one
iteration to the next.  Thus we establish the matrices A. ,  b, I,   and
E , and the column vectors  P.  .  To perform an  iteration, we  evaluate
 k          	                k
 D    and   P   , and using Equations  (63) and  (67)  we complete the
definition of the linear system given by Equation  (58).  We describe
the method used to solve this block-tridiagonal system of linear
equations in Section II.
                                 D-29

-------
      Upon solving Equation  (58)  for values  of    6c.  .   , we  use
Equation  (57) to obtain an improved estimate of
. C.  ...  We terminate the iterative process when
JC 1 , J , K
6
                              n+l
                              i,J,k
                            m n+l
for all  H  and  k  in each column of grid  cells.   For  the  present
simulation study,  e = 0.01.

      'Fhe Evaluation of W
      One of the assumptions commonly made  in  airshed modeling is
that turbulent atmospheric flow  is incompressible.   Upon making this
assumption, we can write:
                            -
                           3x     3y    3z
                                                                        (71)
or, in  (£,ri,p) coordinates
                           3uAH    3vAH
where
We note that  W ,  as defined in Equation (9) , is given by
                                                                        (72)
                                                                        (73)
                                                                        (74)
                                  D-30

-------
Thus we vise Equation  (72)  to solve for W and Equation (74) to evaluate
W.
      Equation  (72) is  integrated subject to the following boundary
condition :

                           W  = 0        at p = Q                           (75)

This condition restricts  the wind to flow parallel to the terrain at
ground level.  To obtain  values  of W  from Equation (72) , we write the
following finite-difference  equation:

+  £fi-  ["(vAH)J  . ,  .  -  (vAH)?  ..,..1
   An  I     i,D-*ifk        1,3+^Jfk I
                                                                          (76)
                                            for k = 1,2, ... ,K
where
      w"     - o
      "l.J.* "

Note that Equation  (76)  requires knowledge of the vertical  structure of
the horizontal wind components.  In the  present simulation  study,  u
and  v  are  taken to  be invariant  with height.
                                D-31

-------
      This completes the presentation of Step III of the numerical
integration procedure.  We now turn our attention to discussing the
stability characteristics of the difference equations.  The basic
statement that we can make with regard to stability is that if each
step is stable, then the entire procedure is stable.  Thus we ex-
amine each step individually to determine the conditions for stabi-
lity, and we begin with Step I.  To facilitate the analysis, we
assume the u, AH and K  are constants, and using the stability cri-
teria given by Price et al. (1966) , we obtain
                                          f
                                  ~>    -*

We obtain a similar expression for Step II after making the analogous
assumptions

                            ViT   W       1
                             4"  * 4,»     -  2

Typically, the left-hand sides of these expressions are less than 0.3
for a grid spacing of two miles and a time step of four minutes.
      Step III is more difficult to analyze since the equations are
both coupled and nonlinear.  We can, however, state the conditions for
stability when the equations are linear, that is, when R. = 0.  Since
we have employed the Crank-Nicolson treatment in this step, the solu-
tion is stable for all values of  AT.  We must rely on numerical ex-
perimentation and experience to ascertain if the treatment of the full
nonlinear set of equations is stable.  We can report that, as of this
writing, no evidence of instability has ever been observed.
                                 D-32

-------
II.   METHODS FOR SOLVING THE LINEAR SYSTEMS OF EQUATIONS

      In Step III of the numerical integration procedure, we must
solve two systems of linear equations, each of the form

                                Ax = y

where the characteristics of the matrix  A  distinguish the first
system from the second.  The two forms assumed by matrix  A  are
tridiagonal and block-triagonal, respectively.  A tridiagonal matrix
of order  n  can be written in the following manner:
                 1,1
1,2
                                                         n-l,n
The block-tridiagonal matrix assumes the same form, but the distin-
guishing feature is that each element, A.  ., is itself an m x m matrix.
Thus the order of the matrix in this case'is n-ra .   We note that the tri-
diagonal matrix is actually a block-tridiagonal matrix with m=l.  The
algorithm employed to solve the block-tridiagonal system is equivalent
to that used to solve the tridiagonal  system, but we shall treat each
separately for two  reasons.  First, the recursion formulas for the tri-
diagonal system are particularly simple to formulate; whereas, the re-
cursion formulas for the block-tridiagonal system involve matrix inver-
sions which need to be discussed in further detail.  Second, we wish to
formulate the algorithms using notation similar to that found in the
computer codes described in Volume II.
                                 D-33

-------
A.    The Algorithm for the Tridiagonal System

      As the algorithm we have adopted for the solution  of a tridia-
gonal system of linear equations is described in  detail  by Conte
(1965), we only summarize the recursive formulas  here.   We wish to
solve the following system of linear equations:

o
a . b , c
n-1 n-1 n-1
• a b
n n





\
x.
x
r
;
l-l
X
n





yi
"*
yn-l
y
n
      The solution vector  x  is obtained from the  following four-step
procedure:
      (1)
                   = yl/tol
(2)   Compute  tu.   and  z.   from

                       C,  ,/0).
                wi = bi '
      (3)  Set: a)  =b  -ac  n/u  ,
                 n    n    n n-1  n-1
      (4)  Solve for the solution vector
                x  =  (y  - a z  ,) A)
                 n    •'n    n n-1 ' n
                                                    i =  2,3,...,n-l
                                D-34

-------
                Vi = Zn-i
                                                       i = 1,2,...,n-l
This algorithm is the basis for the code in subroutine SIMDIG, which
is part of the Atmospheric Pollution Simulation Program described in
Volume II.


B.    The Algorithm for the Block-Tridiagonal System

      The algorithm described in this section is fashioned after that
given by Isaacson and Keller (1966).  Consider the following block-
tridiagonal system of linear equations:
                                Ax = f
                                                                       (77)
which we write as
         B,
                            , B  , C  ,
                          n-1  n-1  n-1
                               n
                                   B
                                    n





xl
x.
X
I
i-l
X
n





fl
f2
fn-l
f
n
where  A , B., and C.  are m x m matrices, and  x.  and  f.
m-component column vectors.  We write the coefficient matri:
                                                            are

product of lower and upper triangular matrices in the following manner:
                                D-35

-------
A = LU =
                                             i  r.
                                                                   i   r
                                                                      h-l
      If we equate each element in  A  with each element in the matrix  formed
      by the product  LU, we can write the following identities
               rl • DI1 cl
                     Di  = Bi -  Ai Fi
                                                                                (78)
                                                i — 2,3,... ,n
                                                        " 2,3,...,n—1
                                                                          (79)
Defining a new vector  y ,  we note that Equation (77)  is equivalent to


                               Ly = f

                               Ux = y

We proceed by solving Equation  (80) for  y,  and then  solving Equation
(81) for x, which is the desired solution.   The first step  can  be  written
as
                                                                                (80)

                                                                                (81)
                                                                                (82)
                                  = Di  l'£i " A
                                                                          (83)
                                      D-36

-------
and the second step is given by
                 n
                     SS   V
                          n
                                                                        (84)
                Vi ~   Yn-i "  n-i n-i+1
                              - F  .x
                              i =  1,2,...,n-l
The algorithm given above is formally correct,  but we can improve  upon
it computationally in the following way.   We write the following equi-
valent expressions for Equations (78) and (82), and Equations  (79)  and
(83):
                                                                        (85)
M")
                          =  c.'.f.  - A
iyi-l)
                              i  =  2,3,...,n
where the parentheses and dotted line  indicate a partitioned matrix.*

      We now summarize the algorithm in its final form.
(1)   Set
      B
(86)
(2)   Solve  for T  and y
(3)   Set
D.  = B.  - A.F.
 ill i-l
                       fi " fi ' Aiyi
                                                     i  =  2 , 3,. .. ,n
*  The partitioned matrix (r,iy,)   is  formed by combining the m x m matrix
   F   with the  m x 1  matrix  y    to form a  m x  (m + 1)  matrix.
                                  D-37

-------
 (4)    Solve  for   F.   and
                                                    i —  2,3,... ,n
 (5)    Set               x  = y
                         n   •'
                             yn-i *   rn-ixn-i+l
      We solve  the m+1  linear systems in steps  (2)  and  (4) by
Gaussian elimination.   The  five-step procedure  given above serves
as the basis  for the  code in subroutine BLCKTD, which is part of the
Atmospheric Pollution Simulation Program described  in Volume II.
III.  DISCUSSION OF THE NUMERICAL METHOD

      In Sections I and II, we haye described the computational procedures
employed in the solution of the governing airshed equations.  In this sec-
tion we recount our experiences in the formulation, testing, and applica-
tion of the finite-difference expressions, as these experiences have had
a direct bearing on our final choice of methods.  We begin by discussing
the horizontal advection-diffusion approximations which form the basis of
Steps I and II.  Of primary concern is the treatment of the advection
terms,  3(uc)/3x  and  3(vc)/3y, since the effects of truncation error in
the approximation of these terms may be important.  Initially, we examined
second and fourth-order schemes given by Crowley  (1968) and Fromm (1969),
respectively.  After conducting tests involving the advection of hypotheti-
cal concentration distributions, we concluded that there was reasonable
justification for adopting the more accurate fourth-order scheme.

      Upon concluding these initial tests, we initiated an investigation
to determine what effect the choice of method would have on the magnitudes
of the predicted concentrations.  Meteorological and source inputs used
were typical of those employed in actual  simulations.   Unfortunately,
                                 D-38

-------
at the outset of this study, we observed that, during Steps I and II
of the integration procedure, negative concentrations were predicted
in the vicinity of coastal power plants by both the second and fourth-
order finite-difference methods.  When the vertical diffusion/chemical
reaction calculations were carried out in these areas (i.e., Step III),
the iterative procedure usually failed to converge.  We concluded that
the centered, high-order difference schemes may predict negative con-
centrations in areas on the grid where concentration gradients are
large.  One such area exists near coastal power plants since NO concen-
tration levels in grid cells located over the ocean may differ substan-
tially from those located just onshore.

      We considered two alternatives that would eliminate this short-
coming in the numerical method.  The first was to adopt another difference
formulation of the advection terms which would not produce negative con-
centrations.  A difference scheme that has this desirable characteristic
is the simple first-order approximation.   Although positive  concentrations
are predicted using this scheme, the method is not as accurate as the
second and fourth-order methods.  The second alternative was to specify
a lower limit on the value of the predicted concentration in any grid
cell.  Thus, whenever a negative concentration is predicted, that con-
centration is set equal to some small value.  Adoption of this alternative
would have the effect of adding additional "sources" of pollutant to the
model.  During the time we were evaluating these two alternatives, we
discovered the uncentered, second-order technique which we eventually
adopted.  The two important features that this scheme possesses are:
1) it is more accurate than the first-order scheme, and 2) it is capable
of predicting positive concentrations in areas where  relatively large  concentration
gradients exist.  We performed advection tests using this uncentered second-
order scheme and found that it was more accurate than the first-order
scheme, but less accurate than the centered second-order scheme given by
Crowley (1968).  We finally adopted the uncentered second-order technique
as a compromise between achieving accuracy on the one hand, and avoiding
the possible calculation of negative concentrations on the other.   Since
the vertical concentration profiles  are generally  rather  smooth and
turbulent diffusion is the primary vertical transport  mechanism, we
have encountered no difficulties in  approximating  the  advection and
diffusion terms in Equation (23)  with finite-difference expressions.
The difference representation of  3 (Wc£)/8o ,  as given in Equation  (45),
is a centered, second-order approximation.


      In addition to studying advection approximations for use in Steps I
and II, we have also conducted a test of the finite-difference formulation
used in Step III.  The primary objective of this test is to evaluate the
numerical treatment of the nonlinear chemical reaction terms.  First we set
all wind velocities and surface emissions fluxes equal to zero in a parti-
cular column of grid cells.  We then established a uniform distribution of
each pollutant in the column of cells and integrated the equations.  The
governing airshed equations reduce to the following form:
                                D-39

-------
                                                                        (87)
which we recognize as the ordinary differential equations employed in
simulating smog chamber experiments.  We then solved Equations  (87)
using the method developed by Gear  (1971) and compared the results.
Using the following initial conditions:


                             HC = 40 pphm
                             NO = 20 pphm
                             O  =  0.5 pphm
                            NO_ =  5 pphm
                             CO = 10 ppm

and rate constants typical of those used in simulation studies of the Los
Angeles airshed, we found that, after two hours, paired sets of results
agreed to within 2% for all species.  Our results were obtained using a
four-minute time step and a relative error tolerance of 1%, while the
results from the Gear routine were obtained using a one-minute time step.

      To determine the effect that the size of the integration time step
has on the predicted concentrations, we performed two simulations, iden-
tical in all respects except for the size of the time step.  We chose the
time period 1100-1300 PST since the photolysis rate constants are largest
at this time of day.  We carried out a normal simulation starting at 0500
PST and terminated it at 1100 PST.  We then used the results obtained at
1100 PST  as the initial conditions for the test.  The two integrations
were performed using one-minute and  four-minute time steps.  For purposes
of comparison we calculated the percentage difference between the two
results in each grid cell at 1300 PST,
                   % Difference =
                                   *  **
                                  C -C
                                     *
                                    C
                                          x 100%
where                        C  = result using a 1-minute time step
                            C   = result using a 4-minute time step.
                                D-40

-------
We then averaged the percentage differences in all grid cells and also
found the largest percentage difference.  The results are presented below:


         Species     Average % Difference    Largest % Difference

          RHC                0.9                     3.1
          NO                 2.2                     6.9
          03                 1.8                     6.3
          NO2                1.4                     2.1
          CO                 0.8                     4.9
          URHC               0.6                     2.9
      Based on these results, we adopted a four-minute maximum time step.
In some instances, however, it is necessary to reduce the size of the
time step.  For example, if a negative concentration is ever predicted,
then the size of the time step is decreased, and the integration sequence
is reinitiated.  We note that negative numbers may be predicted by the
horizontal advection scheme if the time step is too large.  The time
step size is also reduced whenever the iterative procedure in Step III
fails to converge in a reasonable number of iterations (say 20).  The
relative error tolerance was chosen as 0.01 to insure that the Newton
iterative process had converged, or at least that the concentrations were
not changing by more than one percent.

      The final observation to be made with regard to the numerical scheme
deals with the order in which Equations (21) - (23) are processed compu-
tationally.  We note that there is no particular reason to integrate in
the  £  direction first; in fact, it is just as reasonable to integrate in
the  n  direction first.  To reduce the possibility of introducing a bias
into the predictions that might result from employing a fixed integration
sequence, we apply the following rule regarding the order in which Equations
(21) - (23) are integrated:

                                              Computational Order

                                     Step I         Step II        Step III

      odd-numbered time steps:   £ ~ direction  n - direction  p - direction
     even-numbered time steps:   n - direction  £ - direction  p - direction
                                D-41

-------
Thus the integration sequence is  £-n-p, n-C-P i C~n-p» and etc.
In Section I, we give the algorithm for the  C-n-P  sequence.  The
algorithm for the  n-£-p  sequence differs from the   £-n-p sequence
only in that values of
are determined from the finite-difference representation of Equation
(22) , and values of
                               **
are determined from the finite-difference approximation of Equation (21)
                                 D-42

-------
                             REFERENCES
Conte, S. D., "Elementary Numerical Analysis," McGraw-Hill, New York
     (1965).

Crowley, W. P., Monthly Weather Review, 96, 1   (1968).

Fromm, J. C., Physics of Fluids Supplement II, II-3  (1969).

Gear, C. W., Communications of the ACM, 14, 176  (1971).

Isaacson E., Keller, H. B., "Analysis of Numerical Methods," John Wiley
     & Sons, New York  (1966).

Price, H. S., Varga, R. S., Warren, J.E., Journal of Math. Physics, 45,
     301 (1966).

Roth, P. M., Reynolds, S. D., Roberts, P. J. W., Seinfeld, J. H.,
     Development of a Simulation Model for Estimating Ground Level Con-
     centrations of Photochemical Pollutants, Final Report and 6
     Appendices, Systems Applications, Inc., Beverly Hills, California
      (1971).

Yanenko, N. N., "The Method of Fractional Steps," Springer-Verlag,
     New York (1971).
                                D-43

-------