-------
performed. Because the puffs were not, in general, normally
distributed in the vertical, azj does not define the puff as
uniquely as aT or a, T do.
x j. y •*-
The values that Nickola obtained for the puff dimensions
are displayed in Table 2. Note that for puff P5 no az-r value
was given at the 800 m arc. Examination of the data shows that
the vertical concentration profiles bears no resemblance to a
normal distribution at the observation tower. However, a value
for ozj is needed in order to determine the rate at which the
puff is spreading in the vicinity of the 200 m arc. Since the
puff diffuses very slowly iri the vertical under all conditions,
an average azj was computed from the azj's determined for the
other puffs. This value was used for puff P5 on the 800 m arc.
Since it takes a finite time for a puff to pass an arc,
the instantaneous width of the puff will be less than the width
of arc intercepted by the puff, due to the effect of meander of
the wind. However, Nickola (1971) points out that, for a puff
released under the unstable or neutral conditions considered
here, there is essentially no difference in these two widths,
due to the relatively high wind speeds and resulting short peri-
od of puff passage at an arc. Therefore, it is assumed that the
instantaneous and integrated standard deviations are identical.
Since most of the data was collected at the 1.5 m level,
it is at this level that ax and av are valid. However, useful
standard deviations are in the plane passing through the puff
centroid, which may be at a different level. Since the puffs
were released at the surface, it is assumed that their centroids
remain at this level, but that the standard deviations are those
measured at the 1.5 m level.
The a's and b's used to define the standard deviations in
(3) may now be determined using the axi» CTyi and ozj computed by
Nickola at 200 and 800 m from the source. A system of two equa-
tions in two unknowns may be set up for each coordinate direc-
tion. In the downwind direction,
In o = In a, + b. In t2QQ
X200
ln a = ln al + bl ln ^00 ' <43>
Similar equations may be derived in the other directions.
The a^, etc., obtained yield the measured standard devia-
tions at 200 and 800 m from the source, but ax, Gy and az deter-
mined when the puff centroid is not on a measuring arc are valid
only if the time rate of change of puff size is correct.
The puff's diffusion rates are determined using (34). The
22
-------
Table 3. Puff diffusion rates.
time
37.6
40.0
66.4
71.2
experiment
P5
P7
experime
P5
P5
P7
P7
al
2.034
1.247
-»t
mt 2 dt
32.6
33.9
11.7
12.2
a2 33
.02452 .2704
.01544 1.573
1 dav
~2 dt
3.09
3.54
2.21
2.55
bl b2
.8131 1.612
.7718 1.537
!^2
2 dt
.547
.568
.232
.227
b3
.8070
.3446
23
-------
results are displayed in Table 3 for the travel times on either
side of the time the puff centroid crossed the 200 m arc. The
diffusion rates all increase with time except for vertical dif-
fusion in experiment P7. In this case the diffusion rate de-
creases, since ~b% is less than 1/2. This cannot occur under the
assumptions stated in Section 2. It is shown in Section 3 that
the diffusion rate for a Gaussian puff increases with travel
time until the puff is far from the source. However, the assump-
tions used in the previous sections are not applicable in the
atmospheric surface layer.
It was shown in Section 2 that the nature of the function-
al dependence of the puff growth rate upon time changed as the
travel time increased. It must be assumed that this is true for
the growth rate of the analyzed puffs also. Since the puff's
dimensions were measured on only two arcs, only one set of a's
and b's are obtained for each experiment. Therefore, only one
functional dependence is computed using (43), so the puff dif-
fusion rates obtained cannot be valid over the entire travel
time of the puff. The diffusion rates obtained are therefore
a reasonable estimate of the true diffusion rates for some time
interval within the time required for the puff centroid to reach
the 800 m arc. It is assumed that this time interval includes
those observation times when the puff centroid is near the 200 m
arc.
24
-------
SECTION 5
DETERMINATION OF THE ENSEMBLE AVERAGED
CONCENTRATION ESTIMATE
Since the array of samplers is inadequate to completely
define the concentration distribution in a tracer puff, conven-
tional analysis techniques cannot be used here. The data used
at the six times closest to the time the puff passed the 200 m
arc are given in Tables 4 and 5. The analyses are most likely
to be correct at times when the puff centroid is near the mea-
suring arc. However, analyses are produced at fourteen times
in test P5, and ten times in test P7, so that the temporal bound-
aries will be far from the analyses of interest.
A modified analysis, based upon Sasaki's (1970a, b and c)
variational method is performed. First, a concentration model
is fitted as closely as possible to the available data. The
model values are replaced by Nickola's observations at the appro-
priate grid points. The results are then filtered in space and
time using the variational technique. The concentration model
is modified and refitted to the resulting analysis in order to
obtain an estimate of the ensemble averaged concentration.
THE INITIAL CONCENTRATION MODEL
The basic concentration model utilized is defined by (4).
However, the Gaussian model does not account for two processes
which occur in the atmospheric surface layer. Wind speed shear
is present in any boundary layer. Also, since the SS^r tracer
is inert, scavenging occurs only at the surface. These proces-
ses tend to deform the puff from a normal distribution.
In order to account for these processes, the Gaussian model
is modified by making u a function of height such that
u - UQ + G(z - ZQ) , (44)
where uo is the effective transport wind speed, computed by
Nickola (1971) at the observation level zo = 1.5 m. A linear
wind profile is assumed because both surface scavenging and wind
shear are accounted for by this method. This produces a dis-
tinctly different concentration distribution than would be
25
-------
Table 4. Available data for test P5 when the puff centroid is near
the 200 m measuring arc. (From Nickola et al. (1970b))
a. Concentrations
Time 104 106
32.8 1.5 4.7
35.2 .9 5.6
37.6 .7 2.7
40.0 1. 12.5
42.4 .4 2.0
44.8 .3 1.8
b. Concentrations
(nCi/m ) at the 1.5 m level
Azimuth (degrees)
108 110 112 114 116
63.8 120.
62.3 162.
37.2 191.
55.0 166.
11.5 93.
9.7 77.
(nCi/m3)
Time
32.8
35.2
37.6
40.0
42.4
44.8
0 167.9
1 210.9
5 176.7
7 157.4
3 100.1
7 114.4
at the
.8
214.4
161.6
130.0
121.4
110.0
78.0
203.4 64.6
132.5 48.7
127.4 26.8
117.2 14.6
87.3 24.6
58.8 50.7
114 azimuth
Level (m)
3. 6
208.5 191
96.1 54
71.8 29
84.8 26
76.1 34
34.4 16
118 120
9.2 4.5
5.7 3.4
4.1 1.0
2.0 .3
1.1 .7
6.2 1.3
.1
.7
.5
.3
.5
.3
.9
122 124
2.5 .2
2.1 .3
1.1 .1
1.5 0.
.4 0.
.6 .1
26
-------
Table 5. Available data for test P7 when the puff centroid is near
the 200 m measuring arc. (From Nickola et al. (1970b))
a. Concentrations (^Ci/m ) at the 1.5 m level
Azimuth (degrees)
Time 104 106 108 110 112 114 116 118 120 122
56
61
66
71
76
80
b.
.8 .2 0. 2.2
.6 .2 .1 4.
.4 .3 12.3 22.2
.2 2.1 12.7 25.9
.0 1.2 5.8 25.2
.8 0. 1.1 6.7
34.1
15.
49.9
53.4
63.2
32.3
Concentrations (|jCi/m ) at
Time
56.8
61.6
66.4
71.2
76.0
80.8
.8
237.3
240.4
181.2
132.6
90.2
46.8
134.7 212
85.4 213
90.7 160
111.6 107
70.2 71
30.5 44
the 114*
Level
3
221.2
214.0
158.4
75.9
43.1
37.3
.5 126.9 24.8 15. 6.
.6 100.6 58.7 27.6 10.4
.9 68.7 50.9 20.2 5.2
.5 32.2 36.7 21.2 3.2
.5 32.5 22.1 4.8 .8
.9 19.3 10.6 5.1 1.5
azimuth
6.1
234.0
78.2
61.0
24.3
17.8
6.8
124
.3
.9
.2
.3
.1
.8
27
-------
obtained using a traditional power law vertical profile, such
as the one employed by Yeh and Huang (1975) . Since the centroids
of the horizontal concentration distributions are displaced down-
wind with height, the maximum concentration in the forward part
of the puff is aloft, while the maximum concentration in the
rearward portion is on the surface. This is in agreement with
the observations. The G is then chosen such that the model con-
centrations are as close to those observed as possible. it is
implicitly assumed here that while tracer is scavenged only at
the surface, small turbulent eddies with time scales smaller
than the incremental averaging time tend to smooth our discon-
tinuities produced by the surface scavenging and the effect of
wind shear. Tl>e concentration model therefore has a normal dis-
tribution in the horizontal, but the vertical concentration dis-
tribution is not Gaussian.
The variables which must be determined such that the model
fits the observed data as closely as possible are; A(t), aj_, a2/
33, bi, b2/ b3, u(z), x(t) and y(t) . The alf a2, aj, blf b2 and
b3 are obtained by the method described in the previous section
and A(t) is defined in (2).
The x(t) and y(t) may be computed when the azimuth, cp(t) ,
of the puff centroid is known. The azimuth must be determined
at each observation time to ensure optimal fit to the data, so
an interval-halving technique is utilized for this purpose. An
azimuth interval is chosen which contains the optimalAaximuth.
The interval is then halved until the quantity
.
mized, where ^ is the measured concentration and -^ is obtained
from (4) . When the quantity is a minimum the optimal azimuth is
obtained.
Since, data useful in determining the optimal azimuth is
contained only in the measuring arc at 1.5 m elevation, other
data may be ignored in this interval-halving technique. Nichola
computed the transport wind speed at the arc level, so u(1.5 m)
is known. Since the interval-halving technique minimizes the
difference between the model value and the data, it is not nec-
essary to optimize A(t) at this point. Therefore, A (t) is taken
to be 1 without loss of accuracy in determining the optimal azi-
muthal angle.
The optimal azimuth displays some variability, as seen in
Table 6, for two reasons. As the puffs move with the flow, wind
fluctuations cause the puff centroids to meander. Furthermore,
the actual puffs are somewhat non-Gaussian, as seen in Fig. 6
to 9. As the model adjusts to the tracer observations on the
200 m arc, shifts in the optimal azimuth will occur. Therefore,
the variability of the azimuth is not due solely to changes in
the centroid ' s azimuthal position, it is partly due to irregu-
larities in the puff's concentration distribution. This intro-
duces an error in the estimate of the puff's azimuth.
28
-------
Table 6. Parameters optimized to fit the model to observations
Travel time
23.2
25.6
28.0
30.4
32.8
35.2
37.6
40.0
42.4
44.8
47.2
49.6
52.0
54.4
52.0
56.8
61.6
66.4
71.2
76.0
80.8
85.6
90.4
95.2
B(t)
a.
.7708
.4261
2.260
6.854
4.271
1.639
-2.151
-1.922
- .4363
-1.142
.04
- .2584
- .4419
- .577
b.
5.586
2.669
-1.827
.5847
.0433
- .5869
- .4587
- .7582
- .2361
.01481
A(t)
Test P5
1.207
1.008
.8978
.7309
.6701
.5632
.5420
.5794
.4675
.5036
.4494
.5275
.4654
.4514
Test P7
.9881
.8061
.5983
.5047
.4790
.5051
.4292
.4267
.3247
.1775
CD(t)
115.8
112.6
112.6
113.3
112.8
111.8
111.6
111.4
112.1
112.3
111.5
110.9
112.2
113.6
113.7
113.9
114.3
113.8
112.8
112.4
113.0
113.5
113.7
113.4
tp(t) = optimal azimuthal angle
A(t) = coefficient obtained from least squares fit
B(t) = constant obtained from least squares fit
Optimal G for release P5 = .1875
Optimal G for release P7 = .1016
29
-------
*
1 OJ
o
CM
'/I = t
/<•/
*€/•
.
v^\
ID
_
•J
(VI
O
CO
o
CD
O
LU
IU
o:
o
UJ
o
i
i_
rs
s
M
<
u
V
in
O
»;)•
n
*-
A
\1 o
in
-P
CO
-P
O
r-t
0)
(U
iH
E
in
ooooo
oo ooooo ooooo oooo ooooo
— O 0>00^~
-------
D 5
8
?
o I o
O
X
4J fC
W
(U C
•4-1 <]J
W to
O -P
g.
•H C
N O
fO -H
-P
o rtj
rH Q)
CO
(U ,Q
5°
ct
O -H
(1)
rH 4J
•H (C
m
O 0)
M rH
a Q)
>
C 0)
O rH
Q)
-P 0) rl
(0 -P O
M m MH
•P U r)
C -H fl)
0) Ti ft
u c
C -H QJ
O M
U w (0
M to
• re a)
r- g to
tji U fO
•rl -H C
fa EH 10
31
-------
2 |
O*«CM — oio*iotNj —
ui
U) _|
, co m
0)
•P
rd
en
C,
0
-H
-P
(0
S-l
-p
c
(U
o
c
o
u
00
32
-------
M
• O
r^\
CM 13
C
-P (0
10
04J
m
(U
r£ 5^
-p m
•H C
N O
(T3 -H
-P
Q)
cn
O -H
/X^"" H
// :
/ /
/// :
/
o
f
8
8
o
CO
t "
1 CO
"E^:
V «D
5 <°
a "
*•*- *-
>< 0
<^
^
O S < UJ 1 l_
0)
rH -P
•H 10
4-1
O 0)
M rH
a a)
> •
C (U T!
O M 0)
•H g
-P 0) i-l
fO -P O
^ ro m
-POM
C -H Q)
CU T3 ft
U C
C -H (U
O M
U
-------
SOURCE
OBSERVATION SITE
— PUFF CENTROID
X-Dt I
Fig. 10. Relationship between observation and model
coordinates. The cp represents the aximuthal angle to
the puff centroid, while cp1 is the azimuthal angle to
the observation site and ro is the radial distance to
the observation site.
MODEL
ANALYSIS
3000
2000
1000
922
600
700
600
500
f 400
•P 300
^ 200
O
a.
* '
60
50
40
30
20
0 OBSERVATION
10
20.B
23.2 25.6 280 30.4 32.8 352 37.6 400 424 44.8 47.2 496 52.0 34.4
t (sec.)—
Fig. 11.
8, 3.
Concentration profile at grid point 13,
34
-------
Once the azimuth is evaluated, the position of the obser-
vation sites with respect to the model coordinates may be com-
puted. Therefore, the x and y distances needed in (4) may be
obtained. As shown in Fig. 10, the optimal azimuth is the direc-
tion from the source to the puff centroid. Since the azimuth and
range of the observation site is known, x and y may be determined
from the geometry depicted.
The data on the_vertical tower is utilized to obtain the
mean transport wind, u, as a function of height. G is determined
by the interval-halving technique using data from all the avail-
able times. Not only is it undesirable to obtain G at each time,
it is impossible. Since A(t) has not been optimized, when data
at only one time is used, the quantity I(x' ~ X)^ ^s minimized
by forcing G to be either very small or very large, depending
upon whether the puff centroid has reached the measuring arc or
not.
A(t) is the final parameter to be determined. It may be
found using a least squares fit at each observation time. The
regression equation used to obtain A(t) is
X = B(t) + A(t) x' , (45)
where x' is given by (4) with Q1 = Q and B(t) is the value of x
on the ordinate. A comparison of Table 6 with Tables 4 and 5
shows that B(t) is very small in comparison with the observed
concentrations in the vicinity of the puff centroid, and may
therefore be ignored. As seen in Table 6, A(t) generally de-
creases with travel time. Due to the nature of A(t) (see Eq.
(2)), the decrease is most likely a result of the puff losing
mass as it proceeds downwind. However, if the puff were Gaus-
sian, residing on the surface, and totally scavenged upon en-
countering the ground, A(t) could not be less than one. When
the centroid of a Gaussian puff is on the surface, the mass flux
must be away from the surface, towards lower concentrations.
Only the concentration in the lower half of the theoretical puff
is considered scavenged. If no absorption occurs, then the
tracer is reflected from the surface and A(t) = 2.
Since the optimal A(t) is generally less than one, the
Gaussian model with the puff centroid on the surface cannot fit
the data very effectively because there must be some mass flux
towards the surface. As seen in Figs. 6 and 8, the tracer puff
is reasonably Gaussian in the horizontal. However, in the verti-
cal, the puff is somewhat non-Gaussian. The data shows that in
the forward portion of the puff there is a concentration gradient
towards the surface. This is reflected in the model, and is one
of the reasons it was decided to depart from the Gaussian verti-
cal profile.
The optimal A(t) in Table 6 decreases in a somewhat
35
-------
irregular manner and at times actually increases. Due to irreg-
ularities in the tracer puff, A(t) is forced to behave errat-
ically in order to fit the data as closely as possible at each
observation time using the least squares procedure.
Figs. 6 to 9 show that the model does a good job of fitting
the data in the horizontal; an acceptable fit in the vertical is
also obtained. Since the puff is not normally distributed in the
vertical, it is more difficult to model the vertical concentra-
tion profile. For example, in experiment P5, the data shows that
until t = 32.8 sec, when the forward portion of the puff resided
over the 200 m measuring arc, concentrations aloft consistently
exceed those at the .8 m level. However at t = 35.2 sec, the
concentration at the .8 m level greatly exceeds those aloft (see
Table 4). The model is not capable of changing the vertical pro-
file in such a rapid manner. However, at t = 40 sec in Fig. 7,
the .8 m model concentration has become larger than those at
higher levels. At longer travel times the concentration at the
.8 m level is always the greatest. Therefore, the concentration
in the central and rearward portions of the puff is largest near
the surface while in the forward portions of the puff the concen-
tration is largest aloft.
COMBINING MODEL DATA AND OBSERVATIONS
In order to reduce some of the irregularities in the model,
an analysis combining both model and measured data is performed.
The resulting concentrations are to be utilized in the determi-
nation of ensemble averaged concentration estimates. These esti-
mates, along with estimates of the tracer fluxes will be forced
to conform to the diffusion equation, Eq. (5), in which the local
time rate of change of concentration is one of the most important
terms. It is therefore imperative that the grid point concentra-
tions change smoothly in time, if reasonable fluxes are to be
obtained. Since the puff centroid's position fluctuates, and
A(t) changes in an erratic manner, the model concentrations at
a grid point exhibit a high amplitude, short period fluctuation,
as shown in Figs. 11 and 12. Therefore, it is necessary to im-
pose a low pass filter in time on the data in order to damp out
the short period fluctuations.
As seen in Table 6, there are times when A(t) actually in-
creases rather than decreases. This suggests that some source
has added to the tracer mass in the puff at that time. This,
however, is unreasonable, since the only source of tracer is at
the initial time. Smoothing in time eliminates the bogus
sources.
Figs. 6-9 show that the model does not conform to the ob-
servations exactly. Therefore, it is desirable to design the
filter such that there is a tendency to force the analysis to-
wards the collected data. This may be accomplished by using a
36
-------
MODEL
ANALYSIS
300
200
t 100
ST 90
E 80
5 ro
^ 60
X 5°
40
30
20
10
© OBSERVATION
20.8 23.7 256 280 30.4 328 35.2 376 400 42.4 44.8 47.2 496 52.0 54.4
t(sec.)-*
Fig. 12. Concentration profile at grid point 16,
7, 5.
37
-------
higher observational weight at grid points where data is avail-
able. However, it is more important that the resulting analysis
be smooth, rather than conform to the data. Therefore, filter-
ing weights are chosen such that the analyses do not always con-
form .to the data exactly. However, Figs. 6-9 show that in gen-
eral the analyses are closer to the data than the model is.
The concentration model cannot account for many of the
mechanisms by which a puff disperses. For example, a backing or
veering wind will tilt the puff. This produces a disparity be-
tween the model and the tracer puff, which can exist for a sig-
nificant length of time. Since the amount of data available is
minuscule, it is desirable to equip the analysis with a "memory",
so that model-data differences in a certain region of the puff
are "remembered" for a short time before and after the observa-
tion time. This may be accomplished using a quasi-Lagrangian
coordinate system, moving with a mean wind speed. The time fil-
tering alters the analysis at a grid point, which maintains its
relative position in the puff, at times other than the observa-
tion time. This is particularly evident in Fig. 12, where the
dotted line represents the analysis with a larger weight attach-
ed to the observational constraint when data is available. The
effect is also present, although not so noticeable, in Fig. 11.
This effect cannot be produced in an Eulerian coordinate system,
because the puff moves through the grid points, so that a speci-
fic grid point does not retain its relative position within the
puff.
Since transport wind speed shear is accounted for in the
model, the puff remains stationary within the Lagrangian coor-
dinate system only at the 1.5 m level. However, because the
shear is small and the time span over which the "memory" acts is
generally short, the effect of the speed shear upon the analysis
is negligible. Furthermore, the puff is free to move laterally
within the grid system. This will also shift the puff centroid
with respect to the grid points. However, the large time filter-
ing does not allow the puff centroid to shift very rapidly.
Therefore, the time filtering has the additional effect of forc-
ing the Lagrangian coordinate system to maintain its relative
position with respect to the puff centroid, and the inclusion of
observations into the analysis tend to distort the analyzed puff
from its Gaussian shape as well as dictate the position of the
puff centroid. It is important that the puff centroid remain
reasonably stationary with respect to the coordinate system,
otherwise it will be impossible to separate the effects of dif-
fusion on the puff concentration from the effects of advection.
Since measuring devices are located on vertical towers on
arcs of 200 and 800 m radii from the source, a cylindrical coor-
dinate system seems appropriate, because observation sites may
then be placed on grid points in an Eulerian system. In a La-
grangian system, however, the observations can be placed closer
38
-------
to grid points using a rectangular coordinate system, except
near the time the puff centroid is at an observation site. As
the origin of the cylindrical coordinate system moves, the range
and azimuth to an observation site changes. Therefore, distances
between grid points in the vicinity of a stationary measuring arc
are continuously changing. At the end of the grid farthest from
the origin, the grid spacing is much larger than the distance be-
tween observation sites, while at the near end the distance is
much less than that between the sites.
The narrow tracer puffs may be confined to 11 observation
sites on the 200 m arc. The analysis is shown observations only
from these sites. Since the arc from which observations are
utilized is so short, it can be approximated by a straight line.
Using a properly oriented Lagrangian cartesian coordinate system
with a grid row always located at 199 m from the source, the dis-
tance between the nearest grid point and an observation site are
much less than a meter in the y direction, and between .3 and
1.7 m in the x direction. Therefore, observations may be assign-
ed to the nearest grid point with little loss of validity.
In order to spread the effects of the observations and mesh
them with the concentration model, smooth spatial gradients are
required. Therefore, spatial as well as temporal filtering must
be accomplished. The analysis system is incorporated in the
variational functional:
J = / Wx- X)2 + a
V 't
,- - 2
+ a4(ff) + a-5(W 1 dv dt . (46)
The functional is minimized and appropriate boundary conditions
are applied to obtain the second order linear partial differen-
tial equation
-2 7 + -3 + «« + «5 - «1 + "1 X - •
dx dy dz dt
Since this Euler equation is elliptic, it may be finite-differ-
enced and solved numerically using an over-relaxation technique.
Because the measuring devices are not equally spaced on the tow-
ers, the finite-difference scheme must take into account the non-
equal grid spacing. The finite-differencing schemes used in this
research are reviewed in Appendix B, and the derivation and solu-
tion of the finite-difference form of (47) is discussed in Appen-
dix C.
FILTERING CHARACTERISTICS
The filtering characteristics of (47) may be examined by
39
-------
defining a response function, R, where
R = A/A . (48)
A is the amplitude of the filtered analysis, and X is the ampli-
tude of the observations.
Using
- - i(kx + % + *
-------
e
m
(N
J
4J
U
0)
(0
(0
(0
I
O
CM
CM
C5
C
•H
(0
O «
-H O
-P in
u
£ II
tn •!*
c a
O
CO C
0) (d
P5
in
•
. CN
(C
ro II
(•oas)
•H ro
fa c?
41
-------
g
ID
E
to
CO
<0
•J-
N
O
i-q
-P
0)
u
(0
m
-H
fa
CO
I
CO
ro
CD
(M
O
1C
n
foes)
« IO •*
(-38S)
g
IT)
M
O
CD
o
m
in
•H
fa
CQ
(0
0)
-H
fa
42
-------
O
>v
-J
o
= 6
•f
to
Q
(0
lO
in
(N
N
H
O
i-q
-P
(U
U
0)
(0
•H
fa
W
(0
1
(0
en
(•oas)°i
tn
•H
fa
E
in
OJ
d
•H
U)
C ir
O •
-rH CM
4J
U II
c
0) ^
w d
c
o -d
OH C
W fO
-------
e
to
O
>1
E
10
CD
O
X
0)
en
•H
fa
CO
fO
C/3
TJ
« OOCJU>O'4ra>t>Jtt>
^
o
J
4J
0)
CJ
CD
fO
CP
•rH
u>
10
IO
(•oes)
44
(0
0)
u
Cn
•H
fa
-------
smaller than the large scale distribution. This allows the anal-
ysis to conform to the observations more closely. A comparison
of Figs. 14a-d show that as the transverse wavelength decreases
the filtered amplitude also diminishes. This is desirable, since
the observed data shows short wave fluctuations which should be
damped to a certain extent in order to produce a smoother anal-
ysis.
BOUNDARIES
The concentration boundaries utilized in this analysis
technique are composed of generated data. The top and side
boundaries cause no problem, because the puff is completely con-
tained within these boundaries. However, the temporal and bottom
boundaries have a noticeable effect upon the analysis. The ef-
fect of the initial time boundary, computed at a travel time of
20.8 seconds, may be seen in Figs. 11 and 12. The heavy temporal
filtering forces the analysis away from the model value, towards
the boundary value. This is one of the reasons that further
analyses are performed at times as far from the temporal bound-
aries as possible. The high amplitude, short period fluctuations
observed in Fig. 11 near the initial time indicates that the op-
timal azimuth obtained when the puff centroid is far from the
measuring arc may not be a good estimate of the centroids true
position. A considerable amount of puff meander would be re-
quired to produce this effect. This is another reason for not
using the analyses close to the boundaries.
ANALYZED CONCENTRATIONS
The analyzed concentrations are reasonably close to those
observed during the time period of interest. This may be seen
by a comparison of Tables 7 and 8 with 4 and 5, respectively.
However, the azimuth of maximum concentration changes more rap-
idly in the data than in the analysis in test P5, a result of
the temporal filtering. The analyzed concentrations are closer
to those observed in test P7. One reason is the aximuth of max-
imum concentration changes more slowly in this case.
The analyzed concentrations differ from the model values
in several respects. The shape of the distributions in both the
horizontal and the vertical are very similar, but differences due
to the inclusion of observations are evident. Fig. 15 shows the
analyzed concentration on the 1.5 m level, where most of the ob-
servations were taken, for the six central time periods in test
P5. In each of these analyses, observations were available on
one entire grid row. At t = 32.8 sec, the data is on the third
grid row downwind of the puff center. At t = 44.8 sec, the data
is on the third grid row upwind of the puff center. The grid
length is chosen such that data is assigned to consecutive grid
rows at each new time. Since the weight assigned to data obser-
vations is much larger than that assigned to model observations,
45
-------
Table 7. Analyzed concentrations for test P5 when the puff centroid
is near the 200 m measuring arc.
a. Concentrations (nCi/
'm3)at
m ; at
the 1.
, 5 m level
Azimuth (degrees)
Time
32.8
35.2
37.6
40.0
42.4
44.8
104
1
1
1
2
2
2
106
5
6
6
13
7
6
108
48
51
40
51
23
20
110
115
151
171
150
95
75
112
184
220
194
167
114
106
114
190
144
134
121
93
63
116
63
56
41
30
33
45
118
9
8
7
6
6
9
120
3
3
1
1
2
3
122
2
1
1
1
1
1
124
0
0
0
0
0
0
b. Concentrations
Time
) at the 114° azimuth
Level (m)
.8 3 6.1
32.8
35.2
37.6
40.0
42.4
44.8
193
159
137
123
105
76
191
111
89
91
77
40
164
61
37
31
34
18
46
-------
Table 8. Analyzed concentration for test P7 when the puff centroid
is near the 200 m measuring arc.
a. Concentrations (yCi/m ) at the 1.5 m level
Azimuth (degrees)
Time 104 106 108 110 112 114 116 118 120 122
56
61
66
71
76
80
b.
.8 0 0 3 35
.6 0 1 6 24
.4 1 11 22 52
.2 2 12 26 55
.0 2 7 25 59
.8 1 3 9 32
Concentrations (^Ci/m ) at
T ime . 8m
56.8 216
61.6 223
66.4 174
71.2 126
76.0 86
80.8 47
134 212 123 26 13 5
101 216 109 55 23 8
101 165 80 49 18 4
109 110 43 36 18 3
70 72 36 22 6 1
33 44 22 12 6 2
the 114° azimuth
Level
3m 6 . 1m
223 231
214 89
158 65
79 22
45 18
36 7
124
0
1
0
0
0
1
47
-------
o
o
in
CN
.
CM
m
MH II
o »
4-1 01
-d
rH C
(U O
> U
0) (U
rH Cfl
e co
IT) .
• CN
rH OH
-p -
O
•H
-P ^1-
C m
0)
u w
C -H
0
u d
o
T3 -H
(1) 4J
N (0
>i M
rH -P
(0 C
in
in
u cr>
c •
o <£>
u
u
(13 C
e (0
Q
(0
•r-t
o
•rt
-P
(0
M
-P
C
(U
U
C
O
U
e
i
•H
X
(0
e
Q) .
U U
X TJ
(1) C
0
m u
m a)
rH CO
•H in
in II
(0
•P
0)
M ^3
(0 C
en (0
O
48
-------
o
en
•H
C
O
•rH
C
0)
O
C
O
u
g
i
•rH
X
(0
g
4-1
ft
0) •
u w
X ^d
0) C
o
(0 O
IT) 0)
r-l tQ
•H r-
fa n
in II
(D
4-J
(U
e-o
(0 C
»^
m -rH
rH U
a.
•H rH
fa (N
AV
x
<
§
^
o
Q
I I
(0
•H
C
o
-H
4-)
ffl
^
-P
C
0)
u
C
o
u
g
I
•rH
X
ffl
g
0) •
u w
X T3
(u d
o
(0 U
in qj
rH 03
• o
tn •
•rH O
fa <*
Q)
00
in -rH
rH U
•H
fa
49
-------
AV
x
<
Q
2 J
C
O
•H
-P
(0
M
-P
C
-------
there is a tendency to induce short waves in the analysis. How-
ever, inspection of Figs. 15a-f shows that with the weighting
system used this effect is not very noticeable.
Fig. 16 contains the analyses for the six central times at
the 6.1 m level. A comparison with Fig. 15 shows the effect of
wind speed shear. In Fig. 16 the maximum concentration is dis-
placed considerably windward. A comparison of the two figures
further shows that the maximum concentration at the 6.1 m level
is displaced towards the right with respect to the maximum con-
centration at the 1.5 m level. This is due to the extremely high
concentration observed at the 6.1m height on the 114° tower at
t = 32.8 sec. Fig. 12 shows that this observation forces the
analyzed concentration to be much larger than the model concen-
tration at times later than 32.8 sec at this quasi-Lagrangian
grid point. However, Table 4 shows that at later times the con-
centration is much lower at this observation site. Therefore,
at positions nearer to and upwind from the puff centroid, the
concentration is much smaller, so the combined observations pro-
duce a distinct bulge in the right front quadrant of the puff.
These analyses are consistent with the concept under which
the analysis scheme was devised: namely, that the abrupt changes
in concentration distribution are due to a combination of meander
of the puff centroid and irregularities in the puff distribution.
Fig. 16 indicates that the analysis exhibits a distinct departure
from the Gaussian distribution, especially at the shorter travel
times. Furthermore, the puff centroid is still allowed to mean-
der, as shown in Table 9a. Comparison with Table 6a shows that
the meander is considerably damped by the analysis scheme, how-
ever.
Figs. 17 and 18 show that the analyses for the six observa-
tion times closest to the time that the puff P7 centroid passes
the measuring arc, at the 1.5 and 6.1 m levels, respectively.
Note that the tracer is more nearly normally distributed in this
case. As with puff P5, the concentration aloft is greatest for-
ward of the puff centroid. This is demonstrated again in Ta-
ble 5. The observed concentration at the 6.1 m level is much
larger at t = 56.8 sec than at the following times. The combined
effect of wind shear and surface scavenging has again distorted
the puff so that its distribution is no longer Gaussian in the
vertical.
Even though puff P7 traveled for a much longer time to
reach the 200 m measuring arc, its concentration is comparable
to that of puff P5, because the puff is diffusing more slowly in
the former case. The puff appears to be losing mass due to more
rapid scavenging in test P7, but this is partly due to the fact
that the time interval between observations is twice as large as
in the other case.
51
-------
AV
x
Q
Z
o
in
ft to
T)
M C
O O
m u
0)
rH 03
0)
> 03
01 •
f-l CM
CO
•-H II
(1)
4J-H
HI U
0-*
-Hȣ(
-P rH
(0
il CO -
•P -H ;.
C g-'
0) C 10
U O CTi
C-H .
O 4J U3
U (0
Vl II
•O -P
0) C >i
tsl 0) i U
•-H C T3
03 O C
C U (fl
(0 -H •
VD X CM
iH (0 i-H
e
II
(Ji (U
••-I fl X
AV
x
<
°-\
O
2 J
Q
(0
C
o
-H
-P
to
M
0)
O
O
U
6
•H
X
m
e
0)
0) •
U CO
XtJ
nj u
-------
AV
x
-------
AV
X
<
o
Q
cn
•H
O
-H
-P
(0
M
•P
d
0)
u
c
O
u
e
-H
X
m
0)
O •
X co
-------
2
•rH
o
•H
4J
fO
0)
O
§
o
g
•H
X
(0
g
Q)
0) •
U CO
X T3
0) d
o
(0 U
r» o)
•H W
in .
-H •-!
CO ||
(d
•P
fO C
CO (0
-H
U
Cn oo
•H n
fa CN
55
-------
o
z
CO
•H
O
•H
-M
tO
0)
U
c
o
u
X
(0
g
-p
a.
0) •
U CO
x -a
id u
r^ a)
iH CO
-H
fa
CO
rO
r- -H
•H U
•H
EM
I I
AV
O
Z
z
o
w
•H
o
-H
4J
fO
n
•p
d
0)
u
§
o
g
-H
X
(0
g
-P
0,
-------
*1
o
•H
-P
(0
0)
u
c
•o
u
g
•H
X
rt)
g
0)
u
X
-------
AV
O
0) U
r-( 0)
0}
CD
\D m
(1) II
-p »
mm
O -
-H
MrH
•P CO
C CN
(U
U IQ
°§
•O -H
(U -P
N m
>i ^ •
•H -P g
(0 C
C d) m
< U cri
G •
O vo
• U
fO II
CO g
H 3 >,
g <1
. -H
tn X T3
•H US C
fa g (0
oo ° 6
^•00 tVJ Jo
w
-H
C
O
•H
•P
(0
M
4J
C
0)
U
C
O
U
X
a
e
0) •
U
-------
r i^ iii
AV
en
•H
C
O
•rH
-P
(0
M
4J
C
0)
O
C
O
U
g
•H
X
fO
g
0)
-p
a,
0) •
O CO
x -a
a) c
o
re u
oo a)
t-H en
• sf
CO
(0
<0
CO 03
rn
O^
s-d
•
Cnr-
•H ^H
fo rH
x
<
M
-H
d
o
•H
-P
a
o
c
o
u
g
i
-H
g
0)
0)
o
X
(U
en
10 O
03 O
H 0)
CQ
§
CJICM
•H •
ca
(0 II
0) -P
CO C
•co
TS g
00 \
r-t -H
O
• n
en
•H iH
fe 00
59
-------
o
z
ll
to
•H
•c
o
-H
-P
m
n
-M
d
c
o
u
X
a
e
-p
ft
0)
0) •
CO
(0 T)
00 C
rH O
o
JL <"
tn en
•H
en
rd u
0) -P
(C T)
w C
to
•n
0) g
oo \
i-l -H
U
•H CO
T i i
o
(0
•H
C
O
•H
4J
C
(U
U
C
O
u
(0
0)
0)
o •
X
-------
The parameters optimized to fit the initial concentration
model to the analysis are shown in Table 9. The optimal cp(t)
and A(t) in Table 9 are obtained in the same manner as those in
Table 6, except the analyzed concentrations have been utilized.
The scavenging coefficient, p, is on the same order for both
tests, but is more variable in test P7. Since B(t) is very
small, it may be ignored.
The analysis presented here using the indicated weights
has many desirable features. Filtering in time produces smooth
temporal changes in concentration and tends to allow model-data
differences to be preserved for short periods of time. The use
of larger observational weights when measured concentrations are
available forces the analysis toward the observed puff and away
from the model. The spatial filtering reduces the amplitude of
short waves and therefore helps blend the model and data together.
The analysis scheme is designed to yield concentration pat-
terns with the high frequency waves filtered out. This is con-
sistent with the observations used, which were averaged over the
collection interval. However, the analysis cannot be used di-
rectly as an estimate of the ensemble averaged concentration.
CONCENTRATION ESTIMATES
The ensemble average concentration estimate is based on
the concentration model developed at the beginning of the sec-
tion. The assumption that the ensemble average concentration is
approximately normal in the horizontal is reasonable since puff
irregularities are eliminated when puffs diffusing under identi-
cal conditions are averaged, due to the random nature of turbu-
lence. In the vertical, the ensemble average concentration pro-
file will still be affected by the combined effects of wind shear
and scavenging occurring only at the surface. As seen in Figs.
16 and 18, the positions of the concentration maxima at the 6.1 m
level change very slowly with time. Therefore, for the six cen-
tral times in each test, the average position for the concentra-
tion maximum at each height is computed with respect to the posi-
tion of the maximum at the 1.5 m level, using (44) with uo = 0
for the effective transport speed. This allows for the genera-
tion of the ensemble average estimates in a true Lagrangian coor-
dinate system, since the azimuthal variation of the puff centroid
may also be eliminated using an ensemble average model.
The quantity A(t) has already been computed using the anal-
yses and is shown in Table 9. This represents an improved esti-
mate for use in an ensemble average over estimates given in
Table 6, since A(t) decreases with time in a reasonably smooth
manner. However, the scavenging coefficient, p(t), which is
computed using A(t), displays a somewhat erratic nature. There-
fore, for an ensemble average, further smoothing of A(t) may be
desirable. Because so little is known about the nature of the
61
-------
Table 9. Parameters optimized to fit an ensemble averaged esti-
mate to the analysis.
Travel time
cp(t)
A(t)
B(t)
P(t)
32.8
35.2
37.6
40.0
42.4
44.8
56.8
61.6
66.4
71.2
76.0
80.8
B(t)
a.
1.072
.5611
.2351
.0338
.0354
.0516
b.
.1309
.1854
.2930
.2965
.1952
.1095
A(t)
Test P5
.6168
.5855
.5608
.5438
.5163
.5017
Test P7
.7479
.6156
.5234
.4769
.4574
.4153
cp(t)
112.5
112.3
112.1
111.9
112.0
112.1
113.9
114.0
113.8
113.2
112.9
113.0
f3(t)
.0198
.0154
.0172
.0168
.0372
.0266
.0140
.0144
= optimal azimuthal angle
= coefficient obtained from least squares fit
= constant obtained from least squares fit
= scavenging coefficient
62
-------
scavenging in the diffusion puffs or the effect of wind shear,
further smoothing is forgone.
In spite of a lack of knowledge concerning the precise man-
ner in which wind shear and scavenging affect an ensemble aver-
aged puff, reasonable estimates of the ensemble averaged concen-
tration may be obtained using the model described above. The
effect of wind shear and scavenging was examined for the data
available; Figs. 7 and 9 show that the modeling used to describe
these effects, though crude, is reasonably effective. The assump-
tion that the ensemble averaged concentration is nearly normally
distributed in the horizontal is a good one, because of the ran-
dom nature of turbulent flows. Puff P5 shows some displacement
in the transverse direction with height, possibly due to wind
direction shear. However, this effect is not included in the
ensemble averaged estimate, because, on the average, wind direc-
tion shear is negligible in the atmospheric surface layer. Con-
centration estimates obtained above will be used in the varia-
tional formalism developed in the next section.
63
-------
SECTION 6
COMBINING FLUX AND CONCENTRATION ESTIMATES
A variational technique has been developed to combine the
ensemble average concentration estimates and the analytical flux-
es such that the diffusion equation is satisfied. This is not
the first variational technique to use the diffusion equation.
Wilkins (1971, 1972) developed a technique for the purpose of
optimizing objectively contoured patterns of air pollution con-
centrations. He used the diffusion equation as a dynamical con-
straint in a numerical objective analysis technique to remove
pattern inconsistencies between successive isopleth pattern pre-
sentations.
THE VARIATIONAL FORMALISM
The diffusion equation is used here as a strong constraint,
so it is satisfied as closely as the numerical methods allow, by
adjusting both the concentration and flux estimates. This ap-
proach is similar to that employed by McFarland (1975) when he
used the continuity equation to obtain dynamically consistent
wind fields. As shown by Sasaki et al. (1977), the proper fi-
nite-difference scheme must be used in order to satisfy the
strong constraint.
In a Lagrangian coordinate system the diffusion equation
may be written as
dF dF dF
+ ^ + ^ + ;r-^+Bx = 0 • (51)
at
3y
where Fx = x'u'» Fy - X'v'» Fz = x'w' an<^ P i-s defined in Section
2. The fluxes and concentrations all represent ensemble averaged
quantities .
The variational functional which incorporates the concen-
tration, fluxes and diffusion equation is
J =
/ > (
3F
V3(Fy- V2
dF
0X)IdVdt .
(52)
64
-------
The
0 . (57)
The derivation of the finite-difference form of these equations,
a discussion of boundary conditions and the solution of the equa-
tions are included in Appendix D.
The above equations may be combined to obtain a second
order elliptic partial differential equation in X, the Lagrange
multiplier. Once X is determined •£, Fx, Fy and Fz may be obtain-
ed by substitution into (53)-(56). These ^uan^ities satisfy the
diffusion equation and are as close to •%, Fx, Fy and Fz, respec-
tively, as the weight ratios Y2/Yi» Y3/Y1 an<^ Y4/Yi allow, except
at grid points adjacent to or on the boundaries. Therefore, the
outermost three grid points in each spatial dimension are not
available for the final concentration analysis, since the con-
sistent finite difference scheme for second order partial differ-
entials involves two rows of boundary values. in Appendix D,
second order differentials of the concentration estimate appear
in the equation for X. Because it is desirable to obtain good
results at the 1.5 m level, where most of the data is located,
the first three levels of the ensemble averaged concentration
estimate are generated below the surface. The levels at which
concentration estimates are computed are -8.1 m, -4.9 m, -1.7 m,
1.5 m, 4.7 m, 7.9 m, 11.1 m, 14.3 m, 17.5 and 20.7 m. The grid
spacing was chosen to permit convenient computer output. In
order to reduce the truncation error in the finite-difference
scheme, the grid intervals used here are smaller than those used
in the preceding section. In the downwind and crosswind direc-
tions the grid intervals used are 8 and 4.8 m, respectively.
65
-------
ANALYSIS OF A GAUSSIAN PUFF
To analyse the results obtained from the variational tech-
nique, two different cases are considered for test P5. in Case 1,
•)( is generated such that it is as close to the ensemble averaged
concentration estimate as possible, except that no wind shear or
non-uniform scavenging is allowed so that the concentration dis-
tribution is Gaussian. It was shown in Section 3 that the esti-
mated fluxes and concentrations satisfy the diffusion equation
exactly in this case. Therefore, errors in the final analysis
will be numerical, such as truncation errors in the finite-
difference schemes. These numerical errors are not small, al-
though they have been reduced considerably by using the smaller
grid spacing. A measure of the finite-difference errors is the
variance of the residual of the diffusion equation (51), VR. The
variance is .533 for test P5, Case 1. This variance is unaccept-
able, because it indicates that the residuals are of the same
order of magnitude as the terms in the diffusion equation. How-
ever, when the variational technique is applied the variance is
greatly reduced to .000008. This indicates that the residuals
were reduced more than two orders of magnitude. While this is
a satisfactory reduction in the residual, it is not reduced to
the limit of accuracy of the computer. If the last term on the
left hand side of (51) is omitted, then the residual is further
reduced.
The concentrations and fluxes for the two central times
have been examined. They are displayed for the 1.5 m level at
the 37.6 sec travel time in Figs. 19-22. At this level the ver-
tical flux is much less than the horizontal. However, at the
top of the puff the two are the same order of magnitude. The
vertical concentration profile is shown in Fig. 23. Since there
is no wind shear or non-uniform scavenging allowed in Case 1, the
concentration distribution is Gaussian in the vertical as well as
in the horizontal.
The analysis performed above could also have been applied
to test P7. However, as shown in Section 4, the vertical spread
of the puff in this test occurs so slowly that it violates one
of the boundary conditions which is shown in Section 3 to be
valid for a Gaussian puff. Therefore, it makes little sense to
generate a Gaussian puff which conforms to the diffusion rates
of puff P7.
Diffusivities may be computed from the final fluxes and
concentrations. The averaged results for grid points where rea-
sonably large concentration gradients exist in Case 1 are shown
in Table 10. Note that, at the 37.6 and 40.0 sec travel times,
the averaged diffusivities are very close to the true diffusivi-
ties, which are given by (25) for the non-sheared, uniform sca-
venging case. This is to be expected, since the flux and concen-
tration estimates satisfy the diffusion equation exactly in this
case.
66
-------
AV
in
(0 -r-l
OJ U
4-) H
i-H O
CM
Q)
CO (0
(0 -H
U
^ O
O -H
4-1 -P
O -P 00
•rf C •
-MO)'*
td u
M C II
-P O
C U >i
0) <3
U 6
§11
U -H (0
•o to g
0) g
Cn o
fO Q) •
n ^ oo
(U H
> II
(0
• X
(U rH <3
r-l 0)
rQ > M
g 0) T3
Q) ^H C
W O
c g u
W Q)
in ro
(Ti
tn
•H -P
in
0)
g
o
00
CO
4-1
• rH ||
O g
i
,C g
-------
ID g
04 -P •
CD 00
<1> .
-PUN
0)
- to X!
H 1
^ g 1
U 4H <3
•r4
m
0) g O
tjl •
ro cu oo
^ r^
0) EH II
(0 X
• CQ
B CU T)
CU rH C
CQ O
ego
H (U
in CQ
CM
CN
VD
• 4-1
tr> II
-H 4-1
CM (C 4J
68
-------
I I I r
zv
/\
(U .
w r-~
m n
U
II
n
O -P
M-l
«H -H
O U CM
M i •
ft m
en
C o II
O CN
•H N
-P W <3
(C -H
^1 TJ
-P C C
com
0) -H
u -P »
c AS 6
O M
O -P CO
c •
na d) <*
i
(U <
> g
(d ^ »
g g
H
X O
Q)
g oo
W Q) II
C ,£
W EH X
<3
n m 03
w o
•H (1) 0)
fa -P w
AV
CM
Q J
in «.
CMTO
CO -H
(U U
•P dL
CN O
(N
0)
Ul CO
(0 -H
U
C
M O
O -H .
«W -P g
(0
C H CO
O -P •
•H C,
c u .
m
(0
X
(U H
•H 0 »
rQ > M
e
-------
The weight ratios are chosen so that the variance of the
diffusivities, Vx, Vy and Vz, respectively, remain at least an
order of magnitude smaller than the diffusivities, while still
forcing the final concentration to be very close to the esti-
mate. For Case 1 we know that the diffusivities are constant in
space, so the variances should be zero. Using the weight ratios
Y2/Y1 = .2, YS/YI = 1. and Y4/Yi = 5.9, the variances displayed
in Table 10 are sufficiently small. The correlation between the
final and estimated concentrations is almost perfect. Spot
checks were made to ensure that the magnitudes of the two were
also nearly identical.
If the variance of the diffusivities is small enough, then
the fluxes obtained using the average diffusivities and final
concentrations must be highly correlated with the final fluxes.
As seen from Table 11, the correlation between the fluxes is very
high in Case 1. This result is anticipated, since we already
knew that in this case the concept of down-gradient flux is a
good one .
When the correlations in Table 11 are good and the vari-
ances are small in Table 10, then the fluxes are proportional to
the concentration gradients, and the constants of proportionality
are the mean diffusivities. If the mean diffusivities are such
that (25) represents a good approximation, then the diffusivities
are the diffusion rates for the puff. The results displayed for
Case 1, a situation where we have an a priori knowledge of the
diffusivities and fluxes, indicates that the fluxes are propor-
tional to the concentration gradients and the diffusivities are
the diffusion rates for the puff.
ANALYSIS OF PUFFS P5 AND P7
In Case 2 (y_) is generated as described in the last section,
so that it is an estimate of the ensemble average concentration
for puffs P5 and P7. The same type of analysis of the results of
the variational technique used in Case 1 is performed here. How-
ever, in this case we do not have an a priori knowledge of the
results.
The variance of the residual of the diffusion equation is
caused by finite-difference errors and possibly because the flux
and concentration estimates do not satisfy the diffusion equa-
tion exactly in this case. However, VR is not increased over
that obtained in Case 1, being .601 for test P5 and .106 for test
P7. The residuals are again the same order of magnitude as the
terms in the diffusion equation. When the variational technique
is applied, VR is again greatly reduced, to .000008 for test P5
and to .00001 for test P7. Therefore, the residuals are reduced
more than two orders of magnitude for both tests.
The concentrations and fluxes have been examined for the
70
-------
Table 10. Diffusivities and their variances.
Time Kx K^ Vy ky
37.6 32.6 32.8 .053 3.09
40.0 33.9 34.2 .093 3.54
37.6 32.6 32.0 1.33 3.09
40.0 33.9 33.2 2.01 3.54
66.4 11.7 10.7 1.25 2.21
71.2 12.2 11.0 1.38 2.55
A i dcr _ "";_-, * Vx*
Kx = 2 dt 'Kx - JL
A 1
-------
Table 11. Correlation between fluxes.
Time
37.6
40.0
37.6
40.0
66.4
71.2
Cx
1.000
1.000
.9995
.9993
.9960
.9939
Cy
.9994
.9988
.9986
.9979
.9966
.9959
Cz
.9992
.9984
.9976
.9965
.9892
.9834
Case
1
1
2
2
2
2
C = correlation between F and K v y
X X X X
C = correlation between F and K v y
C = correlation between F and K v y
72
-------
-p
(0
in
00
-H ||
*"•§ *
to g O
£ C
» +J fO
CN
-g
0) U
CO (U O
(0 CO •
U CO
CN
M g II
0\
MH -rl X
CXI
X i
rH m en
MH sf T3
•O w O
C! -H u
-H 0)
> X w
to d
CO r-4 V£>
O «H •
o g n
3
T3 g II
<1) -H
tP X -P
rO fd
n g »
0) U
> Zl
C 0)
W r-l n
i^> in w
CN • -H
rH
X
tn d) 3
•H ,C rH
fc -P MH
73
-------
two central times in each test. They are displayed for the 1.5m
level at the 37.6 and 66.4 sec travel times in Figs. 24-33.
Since the variational technique leaves the concentration virtu-
ally unchanged, Figs. 24 and 19 are nearly identical. The down-
wind flux for puff P5, in Fig. 25, is virtually the same as the
flux in Fig. 20. The cross-wind fluxes in Figs. 26 and 21 are
again nearly identical. The horizontal fluxes and concentrations
also compare very well for the other travel times and levels ex-
amined. Of course the flux and concentration patterns at higher
levels in Case 2 are displaced downwind from those in Case 1, due
to the effect of wind shear and surface scavenging. Some deteri-
oration in the comparison of downwind flux was observed at the
highest level examined, 11.1 m, near the top of the puff.
While the horizontal fluxes did not change much from Case 1
to Case 2, the same cannot be said for the vertical fluxes. A
comparison of Fig. 27 with 22 shows that the maximum flux more
than doubled for Case 2 at the 1.5 m level. Furthermore, the
pattern is completely different. Rather than a symetric upward
flux with the maximum at the center of the grid, the maximum up-
ward flux is displaced well upwind in Case 2. The downwind por-
tion of the puff is dominated by downward flux in this case.
The vertical flux pattern has changed completely from Case 1
to Case 2 because the vertical concentration profile is very dif-
ferent in the two cases. The vertical profiles for Case 2 at the
37.6 and 66.4 sec travel times are in Figs. 32 and 33. A compari-
son of Figs. 32 and 23 shows how wind shear and surface scaveng-
ing change the vertical profile. While the maximum concentration
remains relatively unchanged, the largest concentration in the
downwind portion of the puff in Case 2 is aloft, while in Case 1
it is at the surface. Therefore, the concentration gradient in
the downwind portion of the puff at the 1.5 m level is directed
downward in Case 2, but upward in Case 1. Since the flux is
down-gradient, Fz also changes direction in this portion of the
puff. In the upwind portion of the puff, the concentration gra-
dient is larger in Case 2 at the 1.5 m level. This accounts for
the smaller upward flux in Case 1 in this region.
An examination of Figs. 24 and 28 shows that while the tra-
vel time for puff P7 is much longer, the concentration in this
puff is only slightly less than in puff P5, because this puff is
diffusing at a much slower rate. This may be seen from a com-
parison of the fluxes. While the shape of F^ is similar in Figs.
25 and 29, the magnitude is much larger in Fig. 25. The magni-
tudes of Fy in Figs. 26 and 30 are closer, which accounts for
the difference in the shape of the puff. Puff P5 is elongated
more than puff P7 because the downwind flux is relatively larger
in test P5. The vertical fluxes in Figs. 27 and 31 have a simi-
lar pattern, because the vertical concentration profiles in Figs.
32 and 33 are very similar.
74
-------
g
00
•H ||
C
-P Q) T>
W rC C
Q) -P ffl
-P
- g
» U
CN (U O
0) 00
COCN
(« g II
*rH ?S
M U CO
X "O
3 CN C
rH -rH O
4-1 CJ
[Q 0)
rH -H CO
(0
O X ^O
-H 3 •
0)
m
II
CD -H
in x «
(0 rO U
H g Q)
0) co
> Q)
EH g
0) \
rH -H
r^
C Q) in
W rH •
m
g I
•
t^ m w
CN • -H
rH
X
in 0) 3
•H & rH
fa -P 4-1
AV
o o S
*r co £!
Q
i
z
I I I I
(U
A I
-P U
V
-P CO
re
-P II
CO
l<4
»"(
ft fft *&
g C
rH Cl)
0)
CO
(2
W
•
00
CN
•
>
rH
Q)
(U
i — |
g
0
•
00
II
X
^
»
CO
75
-------
The mean diffusivities and the diffusivity variances for
Case 2 are given in Table 10. The correlation between the fluxes
obtained using the average diffusivities and the final fluxes for
Case 2 are shown in Table 11. The variances in Table 10 are low
and the correlations in Table 11 are high. Therefore, the fluxes
are nearly proportional to the concentration gradients, and the
constants of proportionality are the mean diffusivities. This
result could be anticipated for the horizontal fluxes, because
the ensemble average concentration maintains a near Gaussian dis-
tribution in the horizontal. Furthermore, it has already been
shown in Section 2 that the diffusivity concept has a physical
validity for the vertical flux near the surface. However, this
conclusion is not valid for all fluxes; for example, it may be
shown that the heat flux can be counter-gradient in the planetary
boundary layer (Businger, 1973).
Table 10 shows that the mean horizontal diffusivities for
Case 2 are close to the diffusivities computed using (36). There-
fore, Kx and Ky may be reliably estimated using the horizontal
diffusion rates for the puffs. However, this is not true for the
vertical diffusivities. In Case 2 Kz is consistently smaller
than Kz, and for test P7 Kz decreases with increasing travel time.
This result is not surprising since, when the puff is not Gaus-
sian, there is no physical reason why Kz should be the vertical
diffusion rate.
The magnitude of Kz varies as the weight ratio Y4/Yi chang-
es. If Y4/Y1 is made infinitely large, then Kz will conform to
KZ exactly, provided the concentration distribution remains un-
changed. Since the weight ratios used are very reasonable, be-
cause the concentration is relatively unchanged, and the expected
diffusivities are obtained in Case 1, the Kz's computed for Case 2
must be reasonably accurate. This is born out in a comparison of
the vertical diffusivities in Fig. 1 and Table 10. An average of
the analytical diffusivities in the lowest levels compares favor-
ably with the Kz's, which may be regarded as average diffusivi-
ties over the first few meters of the surface layer.
ENSEMBLE AVERAGE CONCENTRATIONS
The ensemble average concentrations obtained are very close
to the analyzed concentrations. A comparison of Fig. 24 with
Fig. 15c and Fig. 28 with Fig. 17c shows that at the 1.5 m level,
where a majority of the data is found, the concentration patterns
displayed agree very well. At levels where the analyzed concen-
trations are not so close to a Gaussian distribution, larger dis-
crepancies exist.
The comparison of the ensemble averaged concentration and
the data may be made in Figs. 6-9. The ensemble average concen-
trations are at grid points which do not coincide with the obser-
vation sites. In general, grid points are close enough so that
76
-------
I I t I
x
<
Q .
|| I i 1
0)
-P rH .
(0 M-l g
•> g 00
-H
-P d II
03 -H
0) g >i
•P <
0)
ȣ 'O
CN -P d
rfl
0) -
03 U g
fO 0)
U w o
VlCN CO
X U X
t-H
m GO .
• CO
•O •<* T3
Cnc:
•H o
£ to cj
g -H 0)
ox™
a) ID
dig
fO 3 II
Jt C
(U -H -P
> X
ro ro «.
g O
0) 0)
rH (U 0)
§
U
Q) n
a) co
n
I
GILO
-H « tn
fa rH -rl
0)
03
C
cn
CN
fl
AV
Q -
III I
-P
(0
r- g «
ft |g
-P -H CO
W C •
<1) -.H 0)
™ ,C >i
-P <
0)
W »T3
03 u C
U Q) (0
Cfl
* 6
O(N
in go
X-.H co
3 U
-HUH
M-l
C rH
•H CN »
£ (0
TO U3 T3
m -H d
O O
M X U
u 3 0)
rH 01
Q) Tt
dig •
fO rj V.O
0) -H
> X II
(0 (0
g I)
(U
rH 0)
-H
• (U U
O H H
oo
g r-
« CN
•rH
fcl
77
-------
AV
X
<]
IW g
g co
04 g
-p'd
-H
g
0)
a
0) »
to u g
rd x
> d ^
0)
tn
(0 (0 »
g o
(1) .0)
rH 0) CD
(U g
to \
C • -H
W rH U
-l -H •
0) -P CO
> (0
m ^ u
-P
cu c
0) <1
u
o c
U rO
CM X
f> rO
< g
tr> a)
CO
<
78
-------
(X
o
-p •
w oo
<0
-P II
- X
CM <]
CQ in
m ti
U C
O
M U
O 0)
M-t CO
0) -sf
•H .
•H i£>
m ^o
o
00°
^ -H
-P U
C 3.
0)
o r^
c m
o -H
u
tn
ro O CN
M -H •
(0
(0 W ||
-p
0) C N
•H o) <
6 ft T3
DOC
tn u ro
C
W g -
on X
g
•
tn
-------
discrepancies caused by this procedure are small. The concentra-
tions at the 1.5 m level in Figs. 6 and 8 are very close to the
observations. Fortunately, grid points and observation towers
coincided very well in both tests, and the vertical concentration
profiles also agreed with the data in Figs. 7 and 9. By using
the concentration estimate described in the last section, the
final concentration is constrained to fit the data.
The results obtained above are not affected by the use of
an average centroid position at each level. if such an approach
were not taken, then an advection term would have to be included
in the diffusion equation, as in (5). However, for the wind
shear observed in tests P5 and P7, this term is an order of mag-
nitude smaller than the other terms at most levels, and may there-
fore be neglected.
80
-------
SECTION 7
CONCLUSIONS
The ultimate value of the K-theory approach discussed in
the previous sections rests on the physical validity of the con-
cept of diffusivity. Lacking this, the prediction of concentra-
tion distributions, which is the ultimate goal of K-theory, re-
duces to rather arbitrary formula-fitting. It is evident that a
puff will be acted upon by a whole spectrum of turbulent fluctua-
tions. Only those fluctuations which are small compared with the
existing distribution of material can be represented by the dif-
fusivity concept. The fluctuations which are on a scale similar
to or greater than that of the puff itself will exert effects
ranging from distortion to bodily movement of the tracer distri-
bution. Therefore, the diffusivity concept alone cannot account
for the dispersion of an individual puff released from an instan-
taneous point source.
In light of this, the research performed has answered sev-
eral important questions, the most important being under what
conditions the diffusivity concept does have physical validity.
It is shown in Section 3 that if the diffusivity is considered
to be proportional to the product of a characteristic length and
a characteristic velocity, the physical interpretation is that
the tracer is spread over the characteristic length by the tur-
bulent eddies. It is further shown that the diffusivities sat-
isfy the diffusion equation and a meaningful set of boundary
conditions when the puff distribution is defined by the Gaussian
model. Therefore, the concept of the tracer fluxes being equiva-
lent to the negative of the product of the time dependent diffu-
sivity derived in Section 3 and the concentration gradient is
valid under the conditions for which the concentration has a nor-
mal distribution; a stationary homogeneous flow in the absence
of external boundaries. The concentration under these conditions
is the ensemble averaged concentration one would expect to ob-
serve in a non-sheared turbulent flow with uniform scavenging
near the surface.
A second important question is whether or not the diffu-
sivity can be identified with some other physical quantity which
may be determined from easily measured atmospheric variables. It
was shown in Section 3 that, under the described conditions, the
diffusivity is the rate of spread of the puff. This is fortunate,
since much work has been accomplished in the determination of
81
-------
diffusion rates, by Lin (1960), and Smith and Hay (1961), for
example. There is, of course, too little data currently avail-
able to determine if these results are directly applicable; how-
ever, a valuable foundation for further research has been laid.
A stationary homogeneous flow, required for a normally dis-
tributed puff, is unrealistic in the atmospheric surface layer,
and there exists an external boundary in the form of the earth's
surface. Therefore, an ensemble averaged concentration analysis
is developed which is decidedly non-Gaussian in the vertical. A
comparison with the data shows that the vertical concentration
profile adequately describes both the effect of wind shear and
also that of surface scavenging. Furthermore, the near normal
ensemble averaged concentration distribution in the horizontal
is shown to conform to the available data. Therefore, the flow
conditions implied by the ensemble averaged concentration, i.e.,
stationary, horizontally homogeneous turbulence, appear to be a
reasonable enough approximation of the true ensemble averaged
flow in the surface layer. This is not surprising since the as-
sumptions of horizontal homogeneity and stationarity are rou-
tinely used in the study of the structure of the surface layer.
Due to the assumption of horizontal homogeneity, the hori-
zontal diffusivity may still be identified with the rate of
spread of the puff. Therefore, horizontal diffusivities which
are time dependent only are likely to have the same physical
validity as in the case of a vertically homogeneous flow for the
shallow puff considered here.
When the concentration is not normally distributed in the
vertical, equating diffusivity to the diffusion rate, as in (25),
has little physical validity, since the vertical concentration
profile can no longer be specified by a standard deviation.
Therefore, it is not surprising that the Kz's obtained in Sec-
tion 6 do not agree with their estimates, even though a large
weight ratio is employed for the vertical flux. The variational
technique employed in the preceding section guarantees that the
ensemble averaged concentration and fluxes satisfy the diffusion
equation (Eq. 51). Since the concentration and horizontal fluxes
are shown to be close to their estimates, as expected, it may be
assumed that the vertical fluxes are also correct. Therefore,
the vertical diffusivities obtained are reasonable approximations
of the average vertical diffusivities in the layer occupied by
the puff.
Because of the good comparison between the KZ's and the
diffusivities in Fig. 1, it may be that the vertical diffusivity
is well defined by (14). The physical validity of this form for
the vertical diffusivity has already been established for a puff
released on the surface. However, the variance of Kz is smaller
than it should be if (14) represents the actual diffusivity. As
seen in Fig. 1, the diffusivity varies rapidly with height close
82
-------
to the surface. Furthermore, in the neutral case the diffusivity
is less than that in the unstable case. However, Kz was found to
be less in the unstable than in the neutral case. Eq. (14) does
not allow a time dependence for Kz. It may not be necessary,
however, since the dependence upon time obtained from the avail-
able data is ambiguous. In test P5 Kz increased with time, while
in test P7 Kz decreased. It is, therefore, evident that this
research did not reveal the precise nature of the vertical dif-
fusivity in the atmospheric surface layer.
Since ambient data upon which the profiles in Fig. 1 are
based is very limited, these profiles are subject to change. It
may be possible to obtain improved profiles by using (14) to
estimate vertical fluxes, employ the variational technique des-
cribed in the last section, and determine KZ at each available
level.
The diffusivity concept appears to be valid for an ensemble
averaged tracer puff diffusing on the surface. Furthermore, the
diffusivities can be identified with other physical quantities
which may be determined from measurable atmospheric variables.
The diffusivities may be used to estimate fluxes, which,
together with the concentration estimates, can be employed in the
variational technique to obtain fluxes and concentrations which
satisfy the diffusion equation. These quantities are assumed to
be the true ensemble average concentrations and fluxes.
83
-------
SECTION 8
BIBLIOGRAPHY
Businger, Joost A., 1973: Turbulent transfer in the atmospheric
surface layer. Workshop on Micrometeorology, American
Meteorological Society, Boston, Mass.
Gifford, F. A., Jr., 1957: Relative atmospheric diffusion of
smoke puffs. J. Meteorol., 14(5), 245-251.
, 1968: An outline of theories of diffusion in the lower
layers of the atmosphere. Meteorology and Atomic Energy,
1968, TID-24190. Clearinghouse for Federal Scientific and
Technical Information, 65-116.
Hildebrand, Peter H., 1977: A radar study of turbulent diffusion
in the lower atmosphere. J. Appl. Meteor., 16, 493-510.
Lewellen, W. S., M. Teske, and C. duP. Donaldson, 1974: Turbu-
lence model of diurnal variations in the planetary bound-
ary layer. Proceedings of the 1974 Heat Transfer and Fluid
Mechanics Institute. Stanford University Press.
Linn, C. C., I960: On a theory of dispersion by continuous move-
ments. Proc. Natl. Acad. Sci., 46, 566-570.
and W. H. Reid, 1963: Turbulent flow, theoretical aspects,
Handbuch der Physik, 8, Part 2, 438-523, Springer-Verlag,
Berlin.
McFarland, M. J., 1975: Variational optimization analysis of
temperature and moisture in a severe storm environment.
Ph.D. Dissertation, University of Oklahoma, Dept. of Meteo-
rology, Norman.
Mood, A. M., and F. A. Graybill, 1963: Introduction to the Theory
of Statistics. McGraw-Hill, New York.
Nickola, P. W., 1971: Measurement of the movement, concentration
and dimension of clouds resulting from instantaneous point
sources. J. Appl. Meteor., 10, 962-973.
, j. v. Ramsdell, Jr., and J. D. Ludwick, 1970a: An inert
gas tracer system for monitoring the real time history of
a diffusing plume or puff. J. Appl. Meteor., £, 621-626.
84
-------
, 1970b: Detailed Time - Histories of Concentrations Re-
sulting From puff and Short Period Releases of an Inert
Radioactive Gas; A Volume of Atmospheric Diffusion Data,
BNWL - 1272. National Technical Information Service.
Pasquill, F., 1970: Prediction of diffusion over an urban area-
current practice and future prospects. Proceedings of
Symposium On Multiple-Source Urban Diffusion Models, Air
Pollution Controle Office Publication, AP-86, Sec. 3.
, 1974: Atmospheric Diffusion. John Wiley and Sons, New
York.
Richardson, L. F., 1926: Atmospheric diffusion shown on a dis-
tance-neighbor graph. Proc. Roy. Soc. (London), AllO,
709-737.
Roberts, J. J., E. S. Croke, and A. S. Kennedy, 1970: An urban
atmospheric dispersion model. Proceedings of Symposium on
Multiple-Source Urban Diffusion Models, Air Pollution Con-
trol Office Publication, AP-86, Sec 6.
Roberts, 0. F. T., 1923: The theoretical scattering of smoke in
a turbulent atmosphere. Proc. Roy. Soc. (London), A104,
640-654.
Sasaki, Y., 1970a: Some basic formalisms in numerical varia-
tional analysis. Mon. Wea. Rev., 98, 875-883.
, 1970b: Numerical variational analysis formulated under
the constraints as determined by longwave equations and
low-pass filter. Mon. Wea. Rev., 98, 884-898.
, 1970c: Numerical variational analysis with weak con-
straint and application to surface analysis of severe storm
gusts. Mon. Wea. Rev., 98, 899-910.
, P. S. Ray, J. S. Goerss, and p. Soliz (1977): Errors due
to inconsistency of finite difference schemes. Submitted
to Monthly Weather Review.
Smith, F. B., and J. S. Hay, 1961: The expansion of clusters of
particles in the atmosphere. Quart. J. R. Met. Soc., 87,
82-101.
Sutton, O. G., 1953: Micrometeorology. McGraw-Hill, New York.
Taylor, G. I., 1921: Diffusion by continuous movements. Proc.
London Math. Soc., (2) 20, 196-202.
Tennekes, H., and J. L. Lumley, 1972: A First Course In Turbu-
lence. The MIT Press, Cambridge, Mass.
85
-------
Van der Hoven, I., 1968: Deposition of particles and gases.
Meteorology and Atomic Energy, 1968, TID-24190. Clearing-
house for Federal Scientific and Technical Information,
65-116.
Wilkins, E. M., 1971: Variational principle applied to numerical
objective analysis of urban air pollution distributions.
J. Appl. Meteor., 10, 974-981.
, 1972: Variationally optimized numerical analysis equa-
tions for urban air pollution monitoring networks. J. Appl.
Meteor., 11, 1334-1341.
Wyngaard, John C., 1973: On surface layer turbulence. Workshop
on Micrometeorology. American Meteorological Society,
Boston.
Yeh,Gour-Tsyh and Chin-Hua Huang, 1975: Three dimensional air
pollutant modeling in the lower atmosphere. Boundary-Layer
Meteorology, 9, 381-390.
86
-------
APPENDIX A
DERIVATION OF THE EXPRESSION FOR DIFFUSIVITY
The derivation of the expression for diffusivity,
is shown in order to reveal the nature of the relevant character-
istic length and eddy velocity. Mixing length theory could be
used to define the diffusivity. However, as Pasquill (1974)
points out, there is a considerable element of vagueness in the
whole idea of a mixing length. Therefore, this approach is not
used since it does not contribute to a clear understanding of the
characteristic length and velocity.
It is well known (e.g., Pasquill, 1974) that, based on di-
mensional grounds, the diffusivity may be determined by the pro-
duct of a characteristic eddy velocity and a characteristic
length scale. Since, as shown in Section 3, the diffusivity is
proportional to the perturbation velocity variance at distances
far from the source, it is assumed that the characteristic ve-
locity is the standard deviation of the velocity perturbation.
Pasquill (1974) states that only those eddies of a size
similar to or less than that of the puff are effective in dif-
fusing it. It is therefore assumed that the characteristic
length in the downwind direction, lc, is related to ax. Since
Ox grows indefinitely, but lc's growth may be limited by the size
of the largest eddies diffusing the puff, the relationship be-
tween lc and ax is taken to be
lc = C(t) ox . (A.2)
At long travel times C(t) must decrease so that 1 remains con-
stant.
Since the diffusivity may be determined by the product of
the characteristic length and velocity,
? ? ^
/"iTi-i1**"!^ / •* i \
x - C]_LU axJ , (A.3)
87
-------
where GI is a time dependent proportionality parameter. The pro-
duct may be replaced by the covariance of £x and u', £xu', since
the correlation between these quantities may be defined as
(A.4)
When tracer puffs are released many times under identical condi-
tions, and ensemble average concentrations are obtained, the re-
sulting distribution of particles must be identical to the con-
centration distribution. Therefore, the variance of £x is equi-
valent to the variance of the concentration as
Substitution of (A. 4) and (A. 5) into (A. 3) yields
Kx = C3 Axu' , (A. 6)
where C3 = C1/C2-
The relationship between &x and u1 must now be examined,
in order to express the diffusivity in terms of measurable quan
tities. Using the fact that the usual laws of differentiation
may be applied to the mean values of fluctuating variables and
their products,
It may be seen from a comparison of Figs. 2 and 3 that
x'(tf) ^ £x. However, the position of the puff centroid with re-
spect to the abscissa, a, does not change as rapidly as the posi-
tion of an individual particle does. Therefore, the time rate of
change of a particle's position with respect to the puff centroid,
C, is nearly equivalent to its rate of change with respect to the
coordinate system, or
dx' d<£x + a> d*x
d~t~ dt at ' >
At short travel times the puff is small and the approximation is
not a good one since the location of the puff centroid changes
rapidly as large eddies transport the entire puff. However, as
the puff becomes larger at longer travel times the eddies tend to
diffuse rather than transport the puff and the approximation be-
comes better. Therefore, at longer travel times,
88
-------
u1 - d* /dt . (A.9)
In those cases where (A.9) is valid (A.7) becomes
, ^l
V1' = 2 dt- ' (A-10)
if (A.5) also is used.
The diffusivity may be related to the rate of spread of
the puff by substituting (A.10) into (A.6), to obtain
C3 dax
Kx = iTdt--
bl
Using the definition a = a,t the diffusivity becomes
3t -L
2 ^l'1
Kx = ax b-jCgt . (A. 12)
The constant 03 may be evaluated by equating (A.12) with (12)
when b~L = 1/2, the condition under which (12) is valid, to obtain
C3 = 1. Therefore, (A.11) is equivalent to (A.I), so that, ex-
cept close to the source, the diffusivity expressed by (A.I) is
in fact proportional to the product of the assumed characteris-
tic eddy velocity and characteristic length scale. When the as-
sumption (A.8) is valid, the characteristic velocity is the stan-
dard deviation of the velocity perturbation, and the character-
istic length scale is shown in (A.2). The product of the pro-
portionality parameter and C(t), Cl, is equivalent to the cor-
relation coefficient, C2. Therefore, the diffusivity may be
uniquely determined from Vu1 , 0X and C2-
89
-------
APPENDIX B
CENTERED FINITE DIFFERENCE ALGORITHMS
The form of the finite-difference algorithm for non-equally
spaced grid intervals may be obtained from a Taylor series expan-
sion. Let P represent any dependent variable, while i, j, k and
£ are the grid indicies along the x, y, z and t axes, respective-
ly. For derivative evaluation at i, j, k or &, only those sub-
scripts different from i, j, k, or H are identified.
Taylor series expansions to approximate P]^_i and P]^+^ are
(Z - -*2
P = P -
k-l - - k-l z 2! --•-
and 9
(zk+l - Z)
Pk+l = P +(zk+l - Z> *zP + - 21 - 7zzP + •••
Neglecting higher order terms (resulting in a truncation error) ,
(B.I) and (B.2) may be solved for vzP. After some manipulation,
" z z ~ zk-l
-
-------
+ P - 2P
"U--1
'
Similar results may be obtained in the x, y and t axes. in Ap-
pendix C, algorithm (B.6) is used on the z axis, since the grid
intervals here are non-equally spaced, and the equivalents of
algorithm (B.7) are used for the other axes.
Significant errors can arise in numerical computations
when a set of differential equations is written in an inconsist-
ent finite-difference form (Sasaki et al. , 1977) . The inconsist-
ent form, exemplified by using (B.5) and (B.7) together, results
in a lack of satisfaction of the governing equations which may
be equal to that obtained when the variational technique is not
employed. A finite difference analog of the second derivative
operator of the form
- 2P + P
VZ(V P) = 7 P = - - 5 - — '
-------
APPENDIX C
NUMERICAL SOLUTION OF THE ANALYSIS FUNCTIONAL
The variational functional in Section 5 may be written in
finite-difference form as
2 2 2
J = T Ca^x - x) + a2(v X) +
-------
Since the Euler equation is elliptic, it may be solved using a
relaxation algorithm. Substituting the finite-difference analogs
(B.6) and the equivalents of (B.7) into (C.4) and rearranging
yields
cc^ a-D oc-3 ao a^ ct-i
__, p «^ ^ j iii i ") ~"
7 (AY)2 Zl Z2 (At)2 2
/ ^ • "~ \ *^ / ~™ ™~ \
2
(Ayr
2
(At)
+ ttl\] = R » (C.5)
where zl == (zk+1 - z)(z - z^.j^) + (zk+1 - z) 2
and z2 = (z - z]c_1) 2 + (z]c+1 - z) (z - z]c_1) .
The over-relaxation factor F is predetermined to speed the re-
laxation process. The residual Rv at the v-th iteration indicates
that the ^'s do not satisfy (C.5) exactly on that iteration. On
the succeeding iterations the residual may be reduced by choosing
a more appropriate ^, such that
-v+1 -v _ Rv __ ._ ,..
v = v + - . (C-6)
^ *•
(Ax) (Ay) Z (At)
This recursion equation may be iterated until the residuals at
all grid points are less than some arbitrary value. The root-
mean-square (R.M.S.) of the residual is a measure of the relaxa-
tion procedure's success in obtaining a solution to the partial
differential equation. For test P7 the R.M.S. was reduced from
79.94 to .01477 in 18 iterations, while in test P5 it was reduced
from 1709 to .006688 in 2,4 iterations.
93
-------
APPENDIX D
NUMERICAL SOLUTION OF THE FUNCTIONAL COMBINING
FLUX AND CONCENTRATION ESTIMATES
The finite-difference form of the variational functional
in Section 6 is
2 A 2
V + Y4(FZ - Fz>
7 F +VF + tf F +
V V \T \7 <7 *7
x. x y y z z
The observational weights yi» Y2' YS an<3 Y4 are prespecified con-
stants in this case. The Lagrange multiplier, X, is one of the
dependent variables.
As in Appendix C, the stationary value of the functional
is found by equating the first variation to 0. Therefore,
6J = 0 = E {2Yl(x - §
§1!
«•%*'*• fl»
2*a*
*Ii?
3 5 58'
n ** ^
«. O tt
o
5-
3-
a
8-
<«
t)
r
0
<
rn
m
C
CO
O rn
S' <
0) 33
O § o
S §2?
0) W Q)
— « 2
Q) O
0)
T)
D
3
^
o
o m
3S
O §
s i
*
o O
>
O
m
z
o
c
en
m
Z
<
X %
o o
z to
2 ^
<" «
Z 2
H f1"
> >
U
u
01
0 ro
H m
m 0)
3 5
a
m
Z
O
-------