HYNHallis
In dip-kalita Model fop tap-Stagnant
Ralsttii Happison
Has 17, 1987
-------
Report No. EPA-910/9-87-179
WYNDvalley
An Air-Quality Model for Near-Stagnant
Flows in Constricted Terrain
Halstead Harrison
WYNDsoft, Inc.
6333 77th Avenue SE
Mercer Island, Washington 98040
May 17, 1987
Prepared for
U S Environmental Protection Agency
Region 10
1200 Sixth Avenue
Seattle, Washington 98101-3183
-------
On Air-Quality Modeling in Complex Terrain,
with Emphasis on WYNDvalley,
A Dispersion Model for Near-Stagnant Flows
Halatcad Harrison
WYNDeoft Inc.
6333 77th Avenue SE
Mercer Island
Washington 98040
(206)-232-1810
The mathematics of air-quality modeling is discussed with special attention to problems
associated with complex terrain. Practical solutions to some of these problems are suggested and
an implementation is described as WYNDvalley, a dispersion model adapted to the study of trace
species concentrations in near-stagnant conditions within complex valley drainage systems. The
model is illustrated by application to the Methow, Twisp, and Chewack valleys of Okanogan
County, Washington.
I. Introduction
To guide wise management of the air we live in, we wish to develop models and
simulations to permit accurate estimates of the concentrations of air pollutants,
appropriately averaged over time and space. This is a difficult task. The
motions of the air are characterized by erratic vascillations that display temporal
scales between seconds and seasons, and spatial scales between millimeters and
thousands of kilometers. Limitations of cost affect the density of observations
that we may make to specify the initial state of the systems we wish to model,
and to follow the evolution of the winds and emissions throughout the period
that interests us. Perhaps surprisingly in these times of exponentially improving
computer technology, we are also limited by the computer speeds and memory
that would be needed for near-rigorous solutions to the equations we pose. For
these reasons practical progress requires approximation.
In the following section I shall explore alternate choices for these approximations,
starting from common assumptions and with moderate rigor. Here I shall
emphasize practical implementation of computer models to be useful under condi-
tions when the flows are dominated by constraints imposed by local topography.
In section III I shall illustrate these choices in more detail through a specific
model, WYNDvalley, and discuss some results from a specific simulation, the sen-
sitivities of these results to the data and to the model's assumptions, accuracies
-------
- 2-
and limitations, ways the model might be improved, and the likely costs and
benefits of effort to effect these improvements. I conclude with a summary and
documentary appendices.
II. Some Mathematics of Tracer Dbpersion in the Atmosphere
Readers anxious to get on with specific details of WYNDvalley may turn immedi-
ately to section ID.
My purpose in including the present discussion is to make explicit a rather long
series of approximations that must be introduced along the way towards solutions
of equations over which ... as it seems to me ... unnecessary arguments sometimes
erupt about the merits of alternate and disparate seeming algorithms. The point
is that excessive partisanship over 'Gaussian Plume1 or 'Box' models is ... or
ought to be ... somewhat vitiated by an appreciation of the many severe approxi-
mations that are common to both these useful approaches.
In the sub-sections that follow I shall highlight those approximations that are
common to both Gaussian and box models, and those separate to each. I shall
conclude this section with a discussion of the virtues and limitations of the two
approaches to air-quality modeling, with emphasis on matching models to the
questions that are to be asked of them.
A. The Dispersion Equation: General Approximations
The central concept of tracer dispersion is, very simply, that mass must be con-
served. The governing equation is:
In [l] C(x,y ,z,t) is the concentration of a trace chemical species whose motions
and dispersion we wish to model as functions of the spatial coordinates * ,y ,z,
and time f; y is a generalized representation of the gradient operator which, in
Cartesian coordinates, becomes d/dx + d/dy + d/d«. The wind velocity is
u(x ,y ,z ,t) and Q (x ,y ,z ,t) is a source or sink. The bold fonts emphasize that v
and u are vectors [while C and Q are scalars] and the central dot between the
vectors denotes the scalar product.
-------
- 3 -
One can imagine solving [l] directly by numerical integration: set up a three
dimensional grid, specify C(x ,y ,z u(x ,y ,z ,t) and Q (* ,y ,z ,t), and then use
Euler's rule .. or other convenient algorithm .. to solve for C(x,y ,x ,tt +dt).
Problems arise, however, because very small increments must be specified in each
of the coordinates. We have all seen near-neighbor wind flags blowing in
disparate directions, and we have all experienced wind speeds and directions that
change abruptly within a few seconds. Suppose then, for illustration, that we
specify a grid with 10 meter resolution over a 10 km x 10 km x 1 km box. This
would require 108 grid elements just to specify the instantaneous concentration
field. I understand that there are parallel machines that can do this, barely, but
for our purposes the most apt comparison may be to the CRAY I, that in one
common configuration can address 'only' 7 x io6 random-access memory locations.
Ignoring this, just for the illustration, with 10 second time resolution for the
simulated fields, and 10 Mhz computational clock speed with 10 cycles per 4-byte
multiplication, this brute-force integration would require approximately 10 com-
puter minutes per simulated hour. That is, the computer would be running only
six times faster than the world, for a 10 km x 10 km x 1 km box that often would
be too small for our interest, anyway.
This illustration is silly for a further important reason: we simply cannot afford
the density of measuring instruments necessary to specify the initial concentra-
tions and winds in those 10® grid elements, nor to verify the precision of our cal-
culations, after the fact.
For these reasons, let us proceed to condition equation [1] for approximations
that are forced upon us anyway by the intrinsic variability of our noisy atmo-
sphere.
First let's try to stretch out the time unit for each integration step. Notice that
in the brute-force example above, a ten second time step is limited to wind
speeds of less than one meter per second, across the grid lengths of ten meters.
But most measurements will report winds and tracer concentrations that have
been averaged over tens of minutes, or hours, and often at speeds exceeding 1
m/s.
-------
- 4 -
To spread out the time steps, let's define:
+r/2
fa ,v ,*)s — / fa >y >* J W [2]
r -r/2
where { may be either the concentration, C ,the source, Q, or any component of
the wind velocity, u. The time, r is a somewhat imprecisely defined subscale over
which we want to smooth our equation variables to match more closely the real
averaging times of our measurements.
With [2] we next define primed quantities such that:
0{x,y ,',t) s 0{x,y,z) + e {x,y,z,t) [3c]
Q(x,1,i,t)ea Q(x,y,z) + J (x,y,z,t) [3q]
+ [3u]
V (x ,y ,z,t)ss v(x ,y,z) +r/ (z ,y ,z ,t) [3v]
u>(x,y,t,t)s3iv{x,ylz)+u! (x ,y ,z ,t) [3w]
In [3] the barred quantities are time-averaged over r and the primed quantities
express the more slowly varying fluctuations at times greater than r, « , v, and w
are the three scalar components of the vectorial wind, respectively parallel to x,
y, and z. Note that by definition:
Fr=qr=u7=tr=w7=0 [2']
Substituting [3u] back into [l], while temporarily neglecting the coordinates y
and z:
T? —T, [(5 + t' Xr + * >]+<« + * ' IU1
Expanding and with a second overbar to denote averaging the operators of [la]
over r,
+ ) + Q + 9" [lb]
The middle two terms within the [ ] average to zero by definition [2'], as does q7.
Then dropping double overbars and noting that HC = uC,
—^[i®+r?l + 5 M
-------
- 5 -
The effect of this manipulation is to lump into the ««• all the temporal fluctua-
tions that occur in times that are short with respect to r, the measurement time.
Since by definition this term is no*-measured, it has to be modeled. What shall
we do with it?
We could set it equal to zero: that after all is a model, too. But this vigorous
approximation excessively degrades the accuracy of what we're trying to do. One
useful approach is to write additional differential equations of the form
du -c „ [—5- -j — 1
— = F[n-J
but this expansion rapidly increases the complexity of the equation set, and clo-
sure problems intrude with higher-order moments.
By far the most common approach is to model the ire7 in terms of the quantities
that are measured at the time scale, r, namely a and C. This is the natural 'next
higher approximation'. Drawing heavily upon analogies to the diffusion of heat in
rigid bodies [an imperfect analogy that will causes trouble later, because the air is
not rigid], Batchelor, Prandl, Reynolds, Richardson, and many others have
explored the utility of postulating:
For obscure reasons this is called the 'exchange' approximation (or the 'austasch'
approximation by those who wish to make it even more obscure). The Kt
defined by [4] is ultimately an empirical coefficient, called a diffusivity. It has
dimensions of lengtlAime"1, and it may be thought of as the product of a root-
mean-square 'eddy velocity' times a scale length over which these eddies are
effective in transporting momentum. The subscript 2 denotes «,hat K, is a scalar
component of the vector K. Some authors generalize [4] still further by defining
K as a tensor, but this is really just mathematical exhibitionism: remember that
the components of K are ultimately empirical and can be measured in ways that
force a vectorial interpretation.
Now we can proceed by substituting [4] into [lc], hereafter dropping the overbar
notation as excessively busy; remember, however, that all the measurements and
modeled variables remain implicitly as averaged over t. We get:
-------
- 6-
*£s=[_l + _L + ±)
dt [ dx dy dzJ
-UC+K— + K— 4-K —
uC + K' 17+ K' 9F + ^~
dz
+ C) = Q
dx dx dy dy dx dz
[5a]
and
dx dy dz
„ dC „ iC „ dC]
' dx ' dy ' dz
if dC j, dC jr dC
jV« iV¦ TJ"" A*
1 dx ' dy ' dz
j. dC „ dC j, dC
' dx ' dy ' dz
is diagonal
(5b]
Condition [5a] is just a statement that the iPiF, Pc7, and Ore7 do not affect the
orthogonality of the * -y -z coordinate system. [This is not equivalent to stating
that the turbulence is isotropic.] Condition [5b] again states that the diffusivities
coordinate system, and [5b] again states that the diffusivities need not be treated
as tensors. In comparison to other approximations that follow, neither [5a] or
[5b] are very severe.
With these approximations [id] simplifies to:
dC a(nC) d(vC) d(ivC)
dt dx dy dz
[le]
d_
dx
dc
K*~di
dy
k¦ dC
dx
j, dC
K'-d7
Q
In [le] the K,,Kt, and K, have so far been left within the exterior spatial
differential operators to emphasize that they are in fact functions of scale. This
property of the atmosphere makes eddy diffusion fundamentally different from
that of heat in rigid solids. Measurements over widely differing scales show that
atmospheric diffusivities are larger for bigger puffs than smaller [See for example:
Gifford, 1982].' At scales larger than about 30 km eddy-diffusivity coefficients
scale roughly as X2, where X is a curvature scale of a diffusing puff, that may or
may not be isotropic.
-------
- 7 -
At smaller scales K,,Kt, and K, vary with diminishing powers of X,, Ar> and X,,
but they do not approach scale independence until puff sizes of millimeters, or
less.
To proceed with any elegance beyond [le] therefore requires differentiation at all
scales, a need perhaps best assuaged by spectral expansions of the variables in
Fourier series. Alternately, with less elegance but perhaps more common sense,
one may ignore the scale dependence of the K( at this step, and proceed with the
differentiations as if the they were scale independent; then at a later step empiri-
cal scale corrections can be imposed a postertori. It is further a common approxi-
mation to neglect the spatial gradients of «, v, and w at this step, though in fact
this is not strictly required ... nor is it a very good idea. However, for simplicity
and because Gaussian models do require this further approximation, I shall make
it here to achieve a 'final1 working version of this series of equations [1]:
dC
dt
—1»
dC
dC
- XD
dc
dx
,dy
bz
+ *'-0 + *-0 + *-0+« M
Let us pause for a moment to reflect on what... if anything ... useful that we
have done by recasting [1] into [if]. Recall that the central idea was to subsume
all the temporal fluctuations that occur over time scales smaller than the meas-
urements' averaging times, r, into terms of the form JFc7, and then through the
'exchange' approximation to parameterize these into K(. This may permit us
both larger time steps and a wider spatial net, because at constant wind speeds
each parcel travels further during the longer times. Unfortunately, however, this
latter advantage is not automatic, as a new and serious difficulty intrudes when
we attempt to integrate [If] numerically with a relaxed grid.
This trouble is perhaps best illustrated with a sketch, as in figure 1. Suppose
that an initial puff concentration is uniformly spread between two adjacent grid
points, as with the prism at the left of the figure. With a wind blowing to the
right, then after a time increment dt the puff might be translated to the position
of the center prism. But with a finite grid the mass distribution of this puff
would be interpreted as split between the two adjacent grid elements, as with the
pair of prisms at the right. Suddenly the puff has 'diffused' to greater width,
purely as an artifact of the finite grid. Unless the translation fortuitously places
the puff exactly over a neighboring grid, which is hard to achieve generally when
-------
uotsnj.^iQ xeaxaauinN
U
U 1
i i r
BB
^
T
3^inE> x j
-------
- 8 -
the winds vary in direction and speed, the numerical diffusion may easily exceed
the physical diffusion. One may minimize this artifact by requiring tight grids,
but this is a step backwards in the whole approach from [l] towards [if].
What is to be done about this new problem? Two approaches are common: brute
force, or a clever transform. Both have virtues and defects. The WYNDvalley
model that will be described in the next section is essentially a brute-force
approach: set up a moderately fine grid (100's of meters) with moderately small
integration steps (100's of seconds), accept moderate numerical diffusion, but
compensate for it by diminishing the external diffusion parameters, K(. The
advantages and limitations of this approach will be described in appropriate
detail later. Because the transform approach is the more commonly used, how-
ever, it is proper here to continue its development for a few more paragraphs, to
illuminate the comparison.
B. Gaussian Models: Special Approximations
In this sub-section I amplify briefly on an additional set of approximations
inherent in Gaussian solutions to [If].
If «, v, and w are approximately constant over short times and spaces, then [if]
is a linear second order differential equation. Several of the approximations lead-
ing to [If] where consciously adapted to achieve linearity, so as to take advantage
of a convenient property of linear equations, that two independent solutions of
[if] may be added to one another to make a a third. Thus if we can derive an
analytical solution to [if] for an infinitesimal puff, then we may superimpose puff
solutions to achieve a general solution. This approach leads to the justly famous
'Gaussian7 models.
While the derivation of a Gaussian solution to [if] is straight forward, it is also
rather lengthy, and it would overbalance the present text. Interested readers
should therefore consult other references [see for example: Pasquill and Smith.,
1983; Longhetto, 1980]. For present purposes permit me to sketch the derivation,
only, for the more common 'plume' version. A more complex 'puff version may
be derived similarly, with somewhat fewer restrictions.
-------
- 9 -
1. Assume a steady-state for times less than r. Temporal changes at times
greater than r will be accounted for by specifying slower variations in the
winds and diffusivities.
2. Neglect diffusion parallel to the wind: that is, set K, = 0.
3. Set dC /dt — 0 and seek a solution for the remaining partial derivatives in
2, p, *, of the form X(x)Y(y)Z{z).
4. Substituting XYZ into [if], and rearranging one may separate the partial
differential equation into three total differential equations, one for x, one for
y, and one for z. These have particular solutions of the form
X(s)=Cle-b
Y(y) = <7asin(nf y) + Cjcos(n, y)
Z (*) = C^sii^n, x) + Cecos(nt x)
5. Set C3 and CA equal to zero, to force transverse symmetry.
6. Set the lateral and upper boundary conditions at infinity.
7. Evaluate the general solution by integrating over n, and n, with Fourier
transforms. The expansion coefficients of these transforms are evaluated by
assuming a profile for Y(y) and Z(x) at a specified z.
8. Choose Gaussian functions for these profiles. This choice takes advantage of
a unique property of Fourier transforms, that the cosine transform of a
Gaussian is a Gaussian also. [Other choices are possible, but most require
solutions involving sums of infinite series.] Then if the real plume is a Gaus-
sian at any x, it will be a Gaussian at all other x.
9. Evaluate k by conservation of mass.
The result of these assumptions and manipulations can be cast in several closely
related ways. Perhaps most quoted of these [Turner, D3., 1967] is:
In [7] Q has a slightly different meaning from that in [if]: here it refers to a total
current [kg/s]; there it referred to a current per unit volume. H is a plume
center-line height above the surface (at r = 0); x as before is a downwind dis-
tance from the source of Q; y is the transverse distance across the plume, and av
-------
- 10-
and a, are dispersion lengths equal respectively to \/2Kf s /o and y/2K, x /« .
The coordinate z is measured upwards from the local surface, and the term in
x +H accounts for perfect reflection of the plume from that surface. Usually ar
and |