NUMERICAL INTEGRATION.OF THE
CONTINUITY EQUATIONS
Appendix D
of
Development of a Simulation Model
for Estimating Ground Level Concentrations
of Photochemical Pollutants
Prepared by
Systems Applications, Inc.
•Beverly Hills, California 90212
for the
A1r Pollution Control Office
of the Environmental Protection Agency
Durham, North Carolina 27701
-------
NUHERICAL INTEGRATION OF THE CONTINUITY EQUATIONS
Appendix D
of
Developr.1ent of a Si!'\ulation Hodel
for Estimating Ground Level Concentrations
of Photochemical Pollutants
Steven D. Reynolds
Report 7l-SAI-l2
June 1971
Prepared by
Syst~s Applications, Inc.
Beverly Hills, California
90212
for the
Air Pollution Control Office
of the Environmental Protection Agency
Durha.'n, North Carolina 27701
under Contract CPA 70-148
-------
TABLE OF CONTENTS
INTRODUCl' I ON
. . . . .. . . .
.....
. . . . . .
. . . .
I.
SELECl'ION OF A FINITE DIFFERENCE TECHNIQUE
. . .
. . .
II.
THE METHOD OF FRACTIONAL STEPS AND ITS APPLICATION
TO THE INTEGRATION OF THE CONrINUITY EQUATIONS. . . .
III. EVALUATION AND DISCUS.5ION OF THE NUMERICAL ~tETHOD . . .
A.
B.
C.
Stabi1i ty . . . . . . . . . . . . . . . .
Accuracy. . . . . . . . . . . . . . . .
computing Storage and T~e Requirements.
. . . . .
. . . . .
.....
~ERENCES . . . .
. . . . . . . . .
......
......
~
0-1
0-4
D-6
0-16
0-16
0-17
0-22
0-24
-------
INTRODUCTION
The simulation model for estimating ground level concentrations
of photochemical pollutants is based on the equations of continuity
for a turbulent fluid in which chemical reactions occur. These
equations are given bYI
aCi
- +
at
ac
u --!. +
ax
"" a ( aCi)
a;c Kx ik
aCi
v- +
ay
aCi
w-
az
+ aay (Ky :;i) + ;z (Kz ::i) + Ri (el' C2'''.' cp)+ Si
1 .. 1,2,...,p
for
XW~X~XE
Ys ~ Y ~ YN
h(x,y) < Z < H(x,y,t)
- -
t ~ to
x,y "" horizontal coordinates
z = vertical coordinate
u,v,w .. three components of average
wind velocity vector
ci m time-averaged concentration of species
K , K , K "" turbulent eddy diffusivities
x y z
S1 "" rate of emission of speciesi from
elevated sources
R1 - rate of production of species i
through chemical reaction
'Xw' ~, Ys' YN "" west, east, south, and north
model boundaries respectively
h(x,y) "" terrain ele~ation
H(x,y,t) = elevation of the inversion base
where
i
Assuming that turbulent diffusion in the x
be neglected, Equation (D-l) simplifies tOI
and
y
directions can
aC1 aCi
-+u- +
at ax
aCi
v-
ay
aCi
+ w-
az
a ~ a~i )
1:1- K -
3z z az
+ Rt + S1
1 - 1,2,..., P
D-l
(D-l)
(D-2)
-------
The initial and boundary conditions are:
initial
ci(~fY'z,to}. fi(x,y,.z)
boundary (1)
, ,
, ' ac '
-~ azi - 2i (Xty,t)
,at
z - h (x,y)
if W < 0,
ac
then -Kz -! . 0
az
then W gi(x,y,z,t)
aCi
.. W ci - K -
Zaz
at
z - H(x,y,t)
(2) ,if W ~ 0,
'aH aH
. where W co w .. u - - v-
ax ay
3H
- -
at
(3)
Ci . gi(x,y,z,t)
at x - 'if (or XE)
y - ¥s (or YN)
where
x and y are at boundaries through which the
prevailing winds ~nter.
and where
Qi(x,y,t)
fi(x,y,z)
CJi(x,y,z,t)
.. surface flux of species i
.. initial concentration distribution of species
.. function expressing the concentration of
species i on the boundary at points of inflow
i
These equations, which are a repetition of Equations (2) in the main text,
are coupled since they share a common arqument in the texms ~ (cl' Ci'...' cp).
We are thus confronted with the problelllS of solving p coupled, nonl near
partial differential equations in four dimensional space (x,y,z,t).
J!'inite difference methods have enjoyed widespread use in the
numerical integration of partial differential equations. In particular,
they have been successfully applied in the solution of a variety of
three dimensional (x,y,t) problems in fluid mechanios and diffusion.
Unfortuna tely, however, there is a dearth of reported experience in the
integration of 'IIOre complex systems of equations, such as Equations (D-l).
. It was thus necessary to investigate the use of modified versions of
techniques developed for simpler systems of equations.
The approach adopted for thb study, a variant of the method
of fractional steps, is based on the repeated but separate solution of
three two-dimensional problems--in (x,t), (y,t), and (z,t) space--at
, each dine .step. The solut:Lon procedure is impllcit and iterative in
(z,t) ~ce, explicit and dbect in (x,t) and (y,t) spaCe. The 'chemical
re'adtlon 'a"tld sour(:8 terms, R1 "and 'Si, ~:te included in the (z,t) portion
of the integration. Truncation errors incurred due to the finite difference
II
D-2
-------
approximation have been examined in tests of second and fourth-order
methods, as applied to the horizontal advective portions of the
fractionated equations--the (x,t) and (y,t) problems. The fourth-
order procedure appears to be significantly more accurate in treating
horizontal advection. It is still unclear, however, whether inclusion
of the higher-order approximation in the overall model is warranted
since a second-or~er method is used in the vertical, or (z,t), portion
of the calculations. Tests of the two approximations are continuing.
Finally, our experience thus far indicates that the procedure (with
second or fourth-order approximation included) is stable and poses no
particular computational difficulties.
In Section I of this Appendix, we discuss the basis for the
selection of the numerical technique adopted. A detailed exposition
of the method is given in Section II. We concluC1e the discussion,
in Section III, with an evaluation of the procedure and a sununary of
our experienoe to date with its use. I
0-3
i
I
I.
-------
I.
SELECTION OF A FINITE DIFFERENCE TECHNIQUE
Equations (D-2) have characteristics of both hyperbolic and parabolic
equations. Many methods have been reported in the literature for the
solution of these classes of equations. In evaluating available methods,
one must consider the following characteristics of each.
a)
stability
b)
accuracy
c)
computer storage requirements
d)
computing time requirements--expressible as the
ratio of computing time to simulation time.
e)
ease of implementation-adaptability of the method
to the solutions of Equations (0-2)
In this section we briefly examine several classes of finite difference
methods with these criteria in mind.
Finite difference techniques are often classified as explicit or
implicit, depending on whether each succeeding integration step is
direct or is based on the simultaneous solution of difference equations.
While implicit techniques are computationally more burdensome than
explicit methods, they do offer the advantage of larger ranges over
which they are stable. Other finite difference methods exist which
are difficult to classify in this way. Typically, these techniques have
the characteristics of implicit methods yet, because of some unique
aspect of the particular technique, involve less burdensome calculations
than are norwally eXPected with an implicit method. Two such techniques
were considered for use in the solution of Equations (D-2), the method
of fractional steps and the method of alternating directions.
The numerical integration technique that was ultimately adopted was
selected by elimination. The reasoninqvas as follows.
a)
Classical implicit methods involve the inversion of extra-
ordinarily large matrices when applied to Equations (D-2),
matrices on the order of 25,000 in size for a grid of
25 x 25 x 10 with four coupled equations. Even with a
grid of substantially fewer nodes, the computational
times would be excessive.
b)
Explici t schemes were ruled out because of the highly
nonlinear nature of the reaction rate tems. When
coupled reaction rate equations are solved for conditions
under which chemical rate processes (and not transport
processes) are limiting, and these equations are thus
expressible as ordinary differential equations, special
D-4
-------
methods must be employed in their solution to ensure
stability (See Appendix B, Section IV). The need for
such methods is, if anything, more pronounced in the
treatment of the coupled, nonlinear partial differential
equations. .
c) The method of alternating directions was attractive, as
it possesses the stability characteristics of classical
implicit techni~es and yet involves the manipulation
and inversion of tridiagonal matrices. A major
s~plification and savings in computational t~e is
thereby realized, when compared to classical implicit
methods, (see Richtmyer and Morton (1967) for further
details). Unfortunately, knowledge of the concentration
distribution (25 x 25 x 10 x 4 values) at three levels
of time (cqrrent and two previous) are required. In
order to store this information, use of external memory
would be required, resulting in ext~nded computation
times.
d) The method of fractional steps involves the fragmentation
of an n dimensional partial differential equation (with
time as a variable) into n-l two-dimensional partial
differential equations. An implicit formulation is
commonly used for those fragments which are parabolic
(i.e., contain second partials or diffusion terms),
an explicit formulation for those which are hyperbolic,
(l.e., advection o~1y). Since diffusion is included only
in the vertical, and since we consider a maximum of ten
strata, the implicit portion does not constitute a
particular computational burden (inversion of a 10 x 10 '
tridiagonal matrix). .
The method of fractional steps thus appears to be preferable to other
classes of techniques for the particular problem under consideration.
The method of fractional steps was developed in the Soviet Union
in the late 1950's. Since difference equations can be written in several
ways, a variety of formulations of this method have been developed
(see Richtmyer and ~~rton (1967) for a list of references). We have
adopted the Crank-Nicolson, treatment for the (z,t) fragment and have
considered both second and fourth-order difference methods given by
Crowley (1968) and Fromm (1969) for the (x,t) and (y,t) fragments.
In the next section we present details of these schemes.
D-S
-------
II.
THE METHOD OF FRACTIONAL STEPS AND ITS APPLICATION TO THE
INTEGRATION OF THE CONTINUITY EQUATIONS
The problem under consideration is the numerical solution of
Equations (D-2), six partial differential equations in three spatial
dimensions and time, with advection in all three directions and diffusion
only 1n the vertical, or z, direction. Four of the six equations are
coupled through the reaction rate terms. In applying the method of
fractional steps, we fragment the partial differential equation into
three separate two-dimensional partial differential equations, in
(z,t), (x,t), and (y,t), with the inclusion of the reaction and elevated
source terms in the (z,t) fragment. The solution in (z,t) is implicit,
while the solution in (x,t) and (y,t) is explicit. Each step in time
is comprised of the solution of each of the three fragments of Equations
(~2). Before presenting the details of the method, however, we wish
. u
to discuss certa~n aspects of the equations and their solution.
O~ep of the Nethod
Both second and fourth-order finite difference approximations
have been considered for the two horizontal integrations, in (x,t)
and (y,t), in which only advection is included. A second-order
scheme was selected for the (z,t) fragment, as a fourth-order
approximation was not justifiable due to the limited resolution
in the z direction. The advantage of a second-order method in
(x,t) and (y,t) is that it involves less computation than a higher-
order method, but at a price of reduced accuracy. As our experience
with both second and fourth-order methods is not yet sufficient to have
made a choice between them, we present both. See Section III, however,
for a more detailed discussion relating our experience with these
alternative formulations.
~nsformation of .the Governing Equations
In Appendix C we discussed at length the unique nature of the
temperature inversion that is so often present over Los Angeles.
One of its most notable features is an increase in elevation of the
base of the inversion with distance from the coast. We wish to
include this feature in our finite difference formulation of
Equation. (~2), and at the same time specify the base of the
inversion as the upper boundary of the airshed. Unfortunately,
since ground elevation and the height of the inversion base vary
with x, y and z,the region being modeled is not rectangular (in
three dimensions). It is thus convenient to perform a change
of variables to a coordinate system in which this region is rectangular.
To do this, we .set .
t -
x
xE - Xw
(1)-3a)
1>-6
-------
nD Y
YN - YS
(D-3b)
. ~ z - h(x,y)
P H(x,y,t) - h(x,y)
.(D-3c)
Under this change of variables, Equations (D-2) become
aCt + u aCi + v aCi +..!.. aCi .. 1 a ~ aCi)
- a["" at) ap OR"Z ap Kp ap
at xE - ~ YN - Ys OR
.)
+ It! + si for i" 1,2,...,p (D-4)
where
~H = H(x,y,t) - h(x,y)
w=w-
u
6h + P !!!!.)
~a~ at
v
(~ + ~) - aH.
~a" P al1 Pat
~-xw
y .- y
N S
for
t > t
- 0
o~t~l
O~I'I~l
o < P < 1
- -
.The initial and boundary conditions arer
initial
Ci(~'I'I,p,to) ~ f1(~'O,p)
boundary (l)
-Kp a C1
!H ap- - Q1(t,l'I,t)
at p .. 0
-~ ac
(2) if W>O, then P i
- or::- - .. 0
~H. 3p
at p .. 1
if W
-------
where
t and n are at boundaries through which the prevailing
winds enter.
Henceforth, we will be concerned with the solution of Equations
(D-4) in a rectangular grid with spacing 6~, An, and Ap. For a
25 x 25 x 10 grid, these values are 0.04, 0.04, and 0.1 respectively,
independent of time.
TelflPOl'a1 Resolutior-
The temporal resolution is associated with the frequency with
which the primary variables are altered. The resolution of available
data is the determining factor. Wind speed and direction are
reported hourly, air quality data are similarly reported, and the
emissions data, based on the emissions inventory reported in Appendix A,
a180 have a temporal resolution of one hour. It thUB seeI!\ed
appropriate to carry out calculations of t~fl upper wind field,
aircraft emissions, and power plant emissions only once an hour,
rather than at every time step. This simplification necessitates
fixing the .height of the inversion base for a period of one hour,
then instantaneously moving it to a new level (over each of 625
columns of air) where it remains fixed for the following hour.*
A~line interpolation procedure is employed to redistribute
the concentration field from the old grid to the new grid.
Nomerzc latul'e
Under the best of circumstances, the nomenclature associated
with the exposition of a finite difference formulation is clumsy.
Our application is no exception, and, in fact, is rather more
burdensome than is common. With apologies, then, the concentration
distribution on the grid is denoted by m n
. tCi,j,k
Ipn m
lCi,j,k = tC[(i-;)6C, (j-;)~n, (k-;)AP,n4t]
where
and
..
is an index denoting chemical species.
i,j,k
are indices related to location in the ~-n-p . coordinate
system.
n
is an index in time.
m
refers to the iteration number.
Returning now to a presentation of the procedure, we remind
the reader that the calculation of n+l given n requires
C . C
three . steps : .. i,j,k I. iit,k
. Computati~nal considerations also affected the decision to adjust
th e inversion base at hourly intervals. The .repetition of the upper
wind field calculation at each time step would greatly increase computa-
tion time.
0-8
-------
Step I:
Step II:
Step III:
Integration in the
Integration in the
Integration in the
p
~
1"1
direction
direction
direction
We first present the calculations for the formulation in which a
second-order approximation is applied in each of the three steps. We
then give the equations for Steps II and III using the fourth-order
formu~ation. As a second-order implicit calculation is carried out
in the p (or z) direction in either case, Step I is unchanged.
Step I - Integ:mtion in the
p
dirBot1-on
en
R. i,j,k.
difference scheme:
We advance
to
e*
R. i,j,k.
by the fOllowing iterative
~. - en
L i,j,k R. i,j,k
At
W (me. . +
+ ih~'k. R. i,j,k+l
AHi,j
cf
R. i,j,k+l
m . n )
- c - C
R. i,j,k.-l R. i,j,k.-l
4 Ap
1
IC 2(AHn Ap)L
i,j
{ "i.j.k-~(~C~.j.k-l + ,C~.j.k-l)
-I. ~ K \ (me * + en )
\Kpi,j,k.-~ Pi,j,K+J,j t i,j,k t i,j,K
+ IL (me *. + cf ) } +.!. ( mR, * + R,n \ +
-~i,j,k.+~ R. i,J,k+l 1 i,j,k+l 2 R. i,j,k R. i,j,~ .
n+1:I
s
R. i,j,k
(D-5)
The boundary conditions are given by: .
(1)
-Kpi,j,;
AHi,j
(R.e~,j'l - 1e~,j,O),.
Ap
n
..Qi,j
(2)
if W ~ 0, then
cf - cf
L i,j,NOZ 1 i,j,NOZ+l
0-9
-------
if W < 0, then
r/: W ( .n n )
Wt i,j,NOZ+~ - 2 tCi,j,NOZ + tCi,j,NOZ+1
- Kfli,jANoz+~
. AHi,jAp
(tci,j,Noz+1 - 1C~,j,NOZ)
where NOZ refers to the mesh point just under the inversion base.
Equations (D-S) must be solved simultaneously for reactive
hydrocarbon, nitrogen dioxide, ozone, and nitric oxide. One possible
scheme is the solution of the fort, nonlinear equations (10 layers x 4
species) by Newton's n\ethod. To avoid this, however, we have chosen
to write the terms
m *
R
.t 1,j,k
in such a way that each interation only involves the inversion of a
tridiagonal matrix for each species--a significant calculational
8implification. (See von Rosenberg (1969) for a description of Thomas'
algorithm for the inversion of tridiagonal matrices). As an example,
consider the term
m *
R
1 1, j,k
where .t = 1 refers to reactive hydrocarbon.
The rate of formation of hydrocarbon through chemical reaction is
given by
d(RH) - -kg(RH) (0) - klO(RH) (03) - kll(RH) (OH.)
dt
which we write:
~* m * m-l * m * m-l *
. 1 i,j,k - -kg 1Ci,j,k RADll,j,k - klO 1Ci,j,k 2 Cl,j,k
m * m-l *
-kil lCi . k RAD2i j k
,J, , ,
where t . 2 refers to 03' RAel to atomic oxygen, and RAD2 to OH. .
E>-10
-------
The only unknown in
~*
1 i,j,k
is
m *
C
1 i,j,k
,
and thus each iterative calculation for reactive hydrocarbon, as with
all other species, reduces to the inversion of a tridiagonal matrix. The
iteration is started using the concentrations at the previous time step,
0c* .. cn
I. i,j,k I. i,j,k
It is terminated when
m+l * m *
C - C '
." i,j,k I. i,j,k
11\+1 *
C
1 i,j,k
< e:
Step II - Integration in the
~
dil'ection
*
Having established values of "ci . k from Step I, we now
employ a second-order method incons4t~ative form (i.e. we write
!£. + u!£ as!£. + 3 (uc) - c 3u) as given by Crowley (1968) ,. to advance
at ax at ax ax'
*
C
I. i,j,k
to
**
C
I. i,j,k
c ** .. c * [1 + fit
I. i,j,k I. i,j,k ("E - ,~)6E:
(UH~,j ,k
Ui-~,j,k) ]
+ fit
(xE ,- XW) fI~
( * *)
F - F
"i-;,j,k I. i+;,j,k
(D-6 )
where
At
(~E - Xw)fI~
F* = ai-;,j,k ( c* + * )
I. i-;,j,k 2 "i-l,j,k I.ci,j,k
a
i-;,j,k
(ai-;,j,k)2( * * )
- C - C
2 I. i,j,k I. i-l,j,k
= ~j/!'!t fit.,
(xE - Xw) fI~
D-ll
0'::'0;. :
." .
. .':"
-------
The boundary conditions are written as
At
(xE - ~)A~
*
F ...
R. ~,j,k
u
~, j,k At *
(~ - Xw)At . R.9'~,j,k
for ~,j,k > 0
At
(~ - ~)At
..
F s:
R. NOX+~,j,k
~OX+~, j,k At:
(~ - Xw)At
..
t.9'NOX+~,j,k
for ~OX+~,j,k < 0
where NOX+~
refers .to the XE boundary.
At points of outflow w~ use a one-sided difference approximation.
particular,
In
**
C 8::
.. l,j,k
'C~,j+l'k + lC~,j-l,k)~
+ Ul,j,k At
(~ - ~)At
[( * . * )0 * J
C + C 2 - C
R. l,j-l,k R. l,j+l,k I. 2,j,k
for Ul,j,k ~ 0
**
C 8::
R. NOX, j, k
'C:OX,j-l,k + R. C:OX,j+l,k);4
UNOX j k At
+ "
(~. - Xw) At
[tC:OX-l,j'k - (tC;OX,j-l,k + tC;OX,j+l,k)"h]
for '\loX,j,k ~ 0
It should be noted that in moving from node to node in the
only
t direction,
At
(~ - ~)A~
*
F
R. i+~, j,k
D-12
-,
-------
need be calculated, since
At
(~ - Xw) Af;
*
F
g, i-l..I,j,k
is available from the previous calculation.
Step III - Integration in the
11
direction
The final step in the integration is analogous to Step II.
tci:~.k a tC~:j.k [1 + IrS ~tYs)A" (~i.j+~.k - ~i.j-~.k)]
+ At
(YN - ys)A'l
( ** **)
g,Gi,j-l..I,k - g,Gi,j+l..I,k
where
At **.
G ..
(YN - ys)A'l g, i,j-l..I,k
ai j-I..I k ( ** **)
'2 ' tCi,j-1,k + tCi,j,k
(B1 j-~ ~'1~ ** **)
- '-' C - C
2 I, i,j,k g, i,j-l,k
Vi . L kAt
B - ,J-~,
i,j-l..I,k (YN - ys)A'l
The boundary conditions are also treated in a manner analogous to Step II.
We now present the fourth-order approximation as given by Fromm (1969).
Steps II' and III', in t and 11 respectively, are as fo1lowsl
step II' ~
Integration i.n the f; direction
We advance
*
tCi,j,k to
**
C .
t i,j,k
by:
** * 1
C -C 1+
R. i,j,k R. i,j,k
At
24 (xE - ~)
.. )
27 u .;. U
[ (i+~~, j,k i-I..I, j,k
-(Ui+3/2,j,k - Ui-3/2,j,k)] I
+ At ( (4)F* - (4)F* )
(XE - Xw)A~ R. i-l..I,j,k R. i+l..I,j,k
D-13
-------
where the left superscript (4) refers to the fourth-order difference operator
given by
At
("E - ~) M;
(4) * .
F
1 i-~,j,k
=C&i-~,j,k 17( C* + C. )
12 1 i-1,j,k 1 i,j,k
- 'C;-2,j,k'+ 1C:+1,j,k) I
(ai-~,j,k) 21 ( *
+ 24 -- 15 tCi-1,j,k
. )
- C
t i,j,k
( * * )1
- C - C
1 i-2,j,k 1 1+1,j..~
(a )3 1 (
+ i-~,j,k - c* +
12 1 1-1,j,k
( * *) 1
+ C + C
1 i-2,j,k 1 1+1,j,k
:(a1-~,j,k)1+ ~ (* *)
+ 24 l- 3 tCi-l,j,k - tCi,j,k
(* * )1
+ C - C
t i-2,j,k . 1 i+1,j,k
1 C ~, j ,k)
a ...
i-~, j ,k
Ui-~,j ,kAt
(xE - ~) A~
where
For mesh points adjacent to the boundary, we use the second-order scheme
given in Equation (D-6). The boundary conditions are the same as those
presented in the second-order treatment of Step II.
Step III'.
Integration in the n direction
We advance
of Step II'.
**
C
1 i,j,k
to
cf+l
t i,j,k
by the ~ direction counterpart
1>-14
-------
cn+1 Ie C** {1+ 6t [27(V' - v )
L i,j,k L i,j,k 24(YN - Ys)~~ i,J+~,k i,j-~,k
- (Vi,j+3/2,k - Vi,j-3/2'k)]}
+ 6t (4)G** - (4)G** )
(YN - YS)6~ 1 i,j-~,k L i,j+~,k
where
6t
(yN - ys)6r',
(4) ** .ai,j-~,k { ( ** **)
G = 7 C + C
1 i,j(-~:: 12 ** )1}i,j-l,k . 1 i,j,k
- c.+ C
& i,j-2,k ~ i,j+l,k
,I
(6 )2 { ( )
+ i,j-~,k 15 C C
24 & ~:j-1,k - 1 ~:j,k
~ **. ** )}
- c - C
1 i,j-2,k t i,j+1,k
( $ i j_1 k) 3 { ( ** **)
+ ,~, - C + C
12 1 i,j-l,k 1 i,j,k
( ** **)}
+ C + C
1 i,j-2,k 1 i,j+l,k
( B i , j -~ , k) It { - ( * *
+ 24 3 1Ci,j-l,k
** )
1Ci,j,k
( ** ** ).}
+ C - C
1 i,j-2,k 1 i,j+l,k
The boundary conditions are treated in a manner analogous to Step II I.
D-15
-------
III. EVALUATION AND DISCUSSION OF THE NUNERICAL ME'lHOD
As was pointed out in Section I, alternativef.Ld.te difference
methods may be evaluated on the basis of stability, accuracy, and
computing time and storage requirements. .These characteristics have
been examined. as they pertain to ~1e use of ~e method of fractional
steps in the solution of Equations (D-4). W~ present our findings
in the discussion that follows.
A.
Stability
A theoretical stability analysis is a difficult undertaking
for nonlinear equations or systems of equations. As Equations (D-4)
are nonlinear, stability considerations were examined in two
sequential steps:
1) Set Ri = Si = 0 in Equations (D-4). Examine the
stability of the resulting linear equations.
2) Examine the stability of the full Equations (D-4)
through numerical experimentation.
Thus, a range of acceptable step sizes is estimated in S~ep 1.
Step 2 is obviously necessary to assess the effect of. the nonlinear
terms on stability. While this approach does not en~ure that
stability probl~s will not be encountered over certain ranges
of the independent variables--as ~e range studied in numerical
experimentation is necessarily limited--it will provide a good
indication as to whe~er problems of instability are likely.
We turn now to the application of the method of fractional
steps to Equations (D-4). In integrating forward in time from
n
C
.. i,j,k
to
n+l
C
R. i,j,k
stability is assured if each fractional step is stable. Thus,
we must examine~e stability of the numerical procedures that
constitute each of the three steps.
Step I is basically the solution of the equation
aCi
at
w aCi
+--
68 ap
=
1 a ( aCi )
682 a; Kp ap
+ Rt + S1
Since this equation is nonlinear, we first set ~ c Si C O.
The procedure commonly applied to the solution of the resulting
linear equation, the Crank-Nicolson finite difference scheme, is
stable for all values of 6t. We must rely on numerical experiments
to ascertain if the reaction and source terms give rise to
C-16
-------
instability over the range of step sizes (in time) of interest--
on the order of one to five minutes. In preliminary studies, no
evidence of instability has been observed. We plan to carry out
additional tests of the accuracy of numeric~l techniques in the
future, these same tests will be useful in further assessing
stability.
Step II is represented by the numerical integration of
aCi
at
+ ( u )
xE - XW
aCi
-
at
=
o
Crowley (1968) has shown that a condition for stability is that 101<1,
where
L
o =
u6t
(xE - Xw) 6t
Our experience thus far has indicated that a is usually less than
0.1. Step III, in the ~ direction, is the direct counterpart of
Step II, in that lal
-------
The effect of grid spacing on computational accuracy can be assessed,
however, rather easily. One need only perform a validation run
for a small portion of the Basin, using first the 2 mile x 2 mile
grid, then a 1 mile x 1 mile grid. The differences observed in
predicted values of concentrations give an indication of the
magnitude of the truncation errors. While this test has not been
carried out, it is planned for the beginning of the next phase of
our studies. .
The effect of the order of the method on accuracy has been
assessed for a simplified form of the continuity equations,
3c
at
+ ( u ) 3c
xE - "w '5"'f +
( v ) 3c 0
YN - yS an."
(D-7)
under conditions of uniform air flow. Following Crowley (1968),
we first examined the effect of transporting a conical distribution
of concentration across a 2 mile x 2 mile grid (see Figure D-l).
The fOllowing conditions applied to the test:
Wind:
Uniform, 612 mph,
(moving northeast
thus U" 6 r.\ph,
from the southwest
at 45°)
v .. 6 mph
Grid Spacing:
XE - Xw = YN - Ys .. 50 miles
~x .. ~Y .. 2 miles (or ~~ .. ~~ .. 0.04)
Initial Concentration
Distribution:
Right circular cone with concentration of
2 ppn at apex, 1 ppm at base, and 1 ppm
at all other points (see Figure D-1)
Time Step:
10 minutes
The combination of wind speed and direction and size of the time
step conspire to move the distribution exact~y one-half a grid
spacing per tirrLe step in both the ~ and ~ directions. In Table D-l,
we-report the errors in the peak concentration that accumulate over
a five-hour period, the time required to traverse 60\ of the Basin.
Note that the fourth-order approximation, as applied to the numerical
solution of Equation (D-7) proved to be substantially lOOre accurate
in predicting peak concentrations than the second-order method.
Furthennore, using the fourth-order approximation, after one-hour
all other predicted concentrations were accurate to within 2\, and
after five hours, to within 4\.
In a second test, we investigated the effects of the same wind
f!~ld on the 0500 PST carbon monoxide concentration distribution
shown in Figure D-2. A fourth-order approximation was applied to
Equation (D-7) in this test, with a time step of ten minutes. After
D-18
-------
1
I ,
:1
2S
1.4
23
21
11
lD
Iii
I~
17
.,
- -- _.-
- - ---1---
!
t. 0
i--
I
I
-.. --1--+
, --T'
-- -1-1-- -- .-. .-r-i---
_dn .-,-~Tr-
.-- --'. -..- .-
S
4
-. -- -...-
3
l.
2. 3 - ~
s
Eo
7
s q
10 II 11 '3 Ie, 11" .17 ,
'. 3,0 L'
12 2.3 2.. 25
Figure D-l.
TEST DISCRIPTION AND ITS PATH OF TRANSPORT
"
D-19
-------
1 ,S .6 7 _t L.!.t} .!1. 12 '~,..!..(. 11£ ..1l...E,:.S 19,2D 1.1 22 ", Jf 2-
~~ M~~ ;- -;= . ~ -~- -:;-~- ::~i-:t~- ~~r-tH-~r~ -~-_:~t~~; :
23 6 ; 7 7 8 9 lO:!l ~ _:~~-~. .1~.~~.~~-~~~.~0 ~.-i-~-~ 7 7 -~W 6 ~ 5 L!. 23
21 7 7 7 8 9 10 I .1~ ~-~~-~ -~ I ~~~l.~~~ ~~: 10 i-~~. 8 8. -~ 7 ~: 'j-~ 21
11 7'"T"77 T 8 . 9 9: 11. 12 ! 14 . 16 16 116 ; 15 i 14 : 13 12 ~ : 10-r;-r 8 8 I 8 \ 7; 7 \ 6 2.1
le 6 I 6 7-i-~+.~- :~.J~~-~~'~!'~T~:' .~~+~~t~:~-~~~~~. ,~l~4~.L~~t.~- _.!.L~L~ 7 7 20
1'1 6 I 6 7+-,-.L.~_. ..9_..~~~, _1~L1~t' ~~- ~~~l~+~~.l~-L~~. 1: ~-~~J~:-i~l-~ .~*,-~1_71.7 1'/
IS ,~,.~--~. ...7+~~. -~~ -l~h"'~ :~-t.1:'L1~~5 ~,:4 ~~-1-~_.~-~.~0!!_._9.. -~l..~I-.?.~~i-' /8
" 6 i 6. 6 :_~~ ~-' 8 9: 11 I 13 .~. 17...c?.: 16 ~... ,13.l.~~.1.~-~!.. 7. ,j 6 i '~~. 17
6 ; 6 I 6' i 7 8 r; 11. ~ 15 17 ; 19 : 18 r~T~5 13 i 12 i 11 i 10 8 8 7 I 6: 5 I 4"
" : I !' I : I
',:~ ~-t._~;-~.-- 9: 10 111 12 14 :6._/~--~.17 116 : 1:- -~._~~~-: 9 '_~ 7~. ~~ 4 IS"
T 5 I 6: 7' . i' ~O--~Ii~l-~: ~~ 15 : 16 ! 15 I 15 : U 12 t1 ; 10 h 8 ~.t 6 -:..l-, 5' _.~ /~
13 ',; l-;'I-'~:-";-:~~-' 1 _12 112__I._~t;.. U ~~: 1~~;j';;r1'2 ~-11~'~~:-; 1""'; 7 7 6! 5 5 13
-.. +"-T-""~'-r- -t- I : I I
~ ~! : I :: : ":::j:' L::; :: I :: ::- :::r;t~; ':f ';-1;- -; r : : --:- -; t-: ~;;
16; 7 : 8 , 10 '10! 10 10 10_1 .9_J~-TI ~ -~ ,~-+'.~L~ -~- .?. ...'. ' 6: 6 6'0
-~~T,'~'-;'~''';- ~O.-t'~.,-~~' , f 8! 7! 7 7 7' 7"! 7 ----+---.
:~~tEIECl~tft-~ .t'J~t ..~4~rT .~. ~~: ~. -~~ ~t -~ ~
, 5 : 5 6: 6 j 6 7 I 8 I 8 I , ; ;, ~.. - '6-1-~ .. - ~ - .r-~' 6 6 6 6 6 6 6 £
I : I
5 ~- ~ ~-II.-~-i-~.-t-~. ~~+- 7 ~- '~I-;~- : .. :.-' -~'--~'~tl~I':-:' -; _..~ .:' _:
" ~-;._~..--S-'i_5-+~ ~.,+~-+-~ ~ .~- _:'~u_L5 ~
: :.:: :::~-> ~,.ti~ ~ -:-~.:~..: ~-I-~'.':+-< i-:.- ;- ; -: ~ : .: ::
4 14 -;! 4 ! 4" 4 5 -;-. -;-~-;. --~-t:.t; ;-t-~ ~ I 4-t.~- 5 --;. 5 - 5 - 5 5
_.~
2. 3 ,.. 5 , 1 8 q 10 II 12. '3 ,,, ,,. " '7 I I q 10 L'
6 5
6 4
Figure D-2.
co CONCENTRATIONS, PPM.
30 September 1969, 500 PST (600 PDT)
D-20
-------
TABLE D-l.
RESULTS OF TESTS OF THE NUl1ERICAL TECHNIQUE
Second Order Fourth Order
Time (hr.) Peak Cone. (ppn) Error (ppm) Peak Cone. Error
o 2.0 0.0 2.0 0.0
1 1.83 -0. 17 1.87 -0.13
2 1.78 -0.22 1.86 . -0. 14
3 1.74 -0.26 1.86 -0.14
4 1.68 -0.32 1.85 -0.15
5 1.69 -0.31 1.84 -0.16
D-21
-------
four hours, the maximum di~crepancy between predicted and
correct values of concentration was 1.5 ppm, with the vast
majority accurate to less than 1 ppm.
While these tests confirm the accuracy of the fourth-order
method of .fractional steps' as applied to horizontal convective
transport (Equation (D-7», the degree to which these resplts are
reflective of the accuracy that can be expected when the method
is applied to the full model is unkno\o.'1'1. The nonuniformity of
flow typical of the meteorology of Los Angeles, the effect of the
nonlinear terms, the coupling with diffusion and advection in the
vertical, the relatively smooth and even distribution of
concentration that prevails over most of the Basin are all factors
that affect accuracy. Also, it is unclear whether the greater
accuracy of the fourth-order approx~~tion in the horizontal
(Steps II and III) is warranted when a second-order approximation
is used for calculations in the vertical (Step I). In other words,
the gre~ter accuracy of the higher-order method may he forfeited
when approximations of differing orders are applied in thesarne
calculation. Thus, a number of questions remain unanswered with
respect to the accuracy of the procedure. .
In future studies, we plan to compare second and fourth-order
methods for Steps II and III when applied to the full airshed model.
If the differences in predicted values do not vary qreatly, this may
be an indication that either errors in the vertical calculations
are of the order Of those associated with the second-order
approximation in the horizontal, or that nonuniform flow character-
istics or the nonlinear terms serve to lessen the difference in
accuracy between secono and fourth-order methods. Conversely,
if the differences in prediction are substantial, the fourth-order
approximation will be adopted for use despite the added computing.
time that this method entails. In the validation studies reported
1n the main text, we have used the second-order method, feeling that
this would prove sufficient for initial evaluation of the model's
performance.
C.
Computing Storage and Time Requir~ents
The full airshed model currently requires a storage capacity
of 315,000 bytes. It may be possible to reduce the storage
requirements in several wayss
1. Exclude a portion of the nodes lying OVer the ocean,
for instance, those more than four miles from shore. This
would result in a 15-25\ reducUon in nodes and thus a
20,000-40,000 byte reduction in storage. Computation time
would also be reduced.
2. Reduce the number of grid points in the vertical direction
from the present 10 to 5-7. Present results indicate that
such a reduction would not have an appreciable effect on the
accuracy of the calculation. However, specific tests are
needed to substantiate this point. .
0-22
-------
3. Use peripheral storage and overlay portions of the
program.
We will investigate each of these possibilities in future work.
As the computer program is currently constituted, two seconds
of computation time (on an IBM 360-75) are required per time step
for the prediction of carbon monoxide concentrations over a
25 x 25 x 10 grid. Thus, if ~t ~ 2 minutes, one minute of
computation time is needed to simulate one hour of actual time.
It should be mentioned that, while we are now using a two-minute
time step in validation runs for carbon monoxide, a longer time
step may well be adequate--perhaps on the order of five minutes.
We plan to investigate the effect of variations in the length
of time steps in the near future.
Preliminary tests of the model with all six species included
indicate that for two-minute time steps the ratio of simulated
time to computing time is approximately eight. Thus, an eight-hour
simulation would require one hour of computation. However, more
computing experience is needed with reactive species included in
the model to better establish this ratio.
,
0-23
-------
REFERENCES
Crowley, w. P., Monthly Heather Review, ~, 1 (19,68).
Fromm, J.E., Physics of Fluids Supplement II, 11-3 (~969).
Richtmyer, R.D., K.' W. Morton, "Difference Methods for Initial-Value
Problems," Second Edition, Wiley-Interscience, New York (1967).
von Rosenberg, D. U., "f.lethods for the numerical 501 ution of Partial
'Differential Equations," American Elsevier Publishing Company,
New York (1969).
I
I
I
D-24
------- |