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