Submitted by
N.V. Gilliani & Associates, Inc.
Huntsville, AL
to
Washington University
St. Louis, MO
as part of
EPA Cooperative Agreement
"Three-dimensional Monte Carlo Model for Interactive Analysis"
April 1996
LAGRANGIAN PARTICLE MODELING AND INPUT REQUIREMENTS
This report covers Year 1 of the work performed by N. V. Gillani & Associates, Inc. (NVGA) under subcontract to Washington University (WU) in the EPA project entitled
"Three-dimensional Monte Carlo Module for Interactive Trajectory Analysis".
There were two NVGA Tasks in Year 1, as follows:
1. Evaluation of mesoscale meteorological models, including a
description of the method of derivation of u', v', w' for use
in the CAPITA Monte Carlo model ( as Met Processors
10 based on available met model outputs; and,
2. Review of plume rise formulas for major point sources, and
selection of "best" choice for implementation in C-MC.
Sections 2 of this report is focused on the work of Task 1, and
Section 3 covers the work of Task 2.
Section 2 describes :
- Lagrangian particle modeling, in general;
- the methodologies for the derivation of u', v', w' , including the likely choice for C-MC;
- the met input requirements of C-MC; and,
- an overview of two candidate mesoscale met models (RAMS and MM4/5) for their suitability as met preprocessors for C-MC, and of output databases currently available and expected to be available in the near future.
Simulation of the transport and dispersion of pollutants over
mesoscale and regional distances from their sources may be performed
by Eulerian or Lagrangian approaches. The Eulerian approach involves
simulations of pollutant concentration changes over a fixed gridded
domain based on explicit dynamic numerical solution of the basic
mass/momentum/energy equations along with appropriate initial
and boundary conditions. This approach permits straightforward
treatment of nonlinear processes (e.g., atmospheric chemistry)
arising from interactions of pollutants from multiple and diverse
sources. The principal disadvantages of this approach are related
to limitations in the spatial resolution of the processes and
inputs, and to computational complexity. The spatial resolution
limitation is imposed mainly by the horizontal grid size which
is particularly coarse in regional models (typically 10 to 100
km). Sub-grid-scale processes must be parameterized in grid models.
This limitation is particularly serious with respect to the treatment
of point sources, because (a) the emissions from such sources
actually attain grid dimensions only after transport and dispersion
over many hours (up to a full diurnal cycle) and long distances
(often hundreds of kilometers), and (b) the chemistry of point-source
plumes is, in general, diffusion-limited and quite non-linear.
Since major point sources have been estimated to be responsible
(in 1985) for as much as 91% and 46%, respectively, of the anthropogenic
emissions of SOx and NOx in the US (Placet et al., 1990), this
limitation of Eulerian models is believed to introduce quite substantial
uncertainty in the results of regional Eulerian models. The Lagrangian
approach is much more suitable for the treatment of emissions
from individual point sources because, in it, individual emissions
can be tracked in a frame of reference which moves with the emission,
and at a level of spatial resolution equal to the current volume
attained by the initial emission. The emissions may be represented
as plumes, puffs or point particles. For plumes and puffs, their
growth initially is predominantly by turbulent dispersion for
which parameterizations of varying degrees of accuracy have been
derived. The significance of turbulent diffusion, however, decreases
during long-range transport, but field observations have shown
that plume spread continues, evidently as a result of mesoscale
spatial circulations as well as wind shear effects induced by
temporal oscillations produced by diurnal and inertial effects
(McNider et al., 1988). The spatial resolution of plumes and puffs
diminishes as they grow in size, introducing progressively increasing
uncertainty, particularly when there is significant wind shear
within the plume dimension. The applicability of Lagrangian plume
and puff models based on traditional parameterizations of dispersion
is therefore very questionable for mesoscale and regional transport
distances. These limitations may be circumvented in Lagrangian
particle simulations in which individual particles are advected
based on the local flow field. Such simulations are commonly called
Monte Carlo (MC) or random walk models; the CAPITA Monte Carlo
(C-MC) model is such a model. Interest in MC models has grown
in recent years, particularly because of their conceptual simplicity,
and their flexibility in applications to simple as well as complex
flow situations. The discussion of the Lagrangian particle (LP)
approach in this section is referenced to dispersion of non-buoyant
and dynamically passive pollutants. It is applicable to plumes
following plume rise.
In the present context, our main interest is in LP modeling of
inertialess particles moving in an Eulerian frame according to
the following basic finite difference equations, and based on
wind data generated by a mesoscale Eulerian grid model:
x(t+Dt) = x(t) + [ + u'(x,y,z,t) ] . Dt (2.1a)
y(t+Dt) = y(t) + [ + v'(x,y,z,t) ] . Dt (2.1b)
z(t+Dt) = z(t) + [ + w'(x,y,z,t) ] . Dt (2.1c) ,
where, (x,y,z) represents particle position in a Cartesian coordinate
system; t represents the current simulation time; Dt
represents the particle advection time step; , and represent the
mean grid-scale (i.e., resolvable/deterministic) components of
the wind field; and, u',v',w' represent the (unresolved/stochastic/parameterized)
sub-grid-scale fluctuations of the wind field. The application
of the mean wind field terms in the above equations achieve the
advection of the particles as well as their relative diffusion
by wind directional shear with height; the application of the
fluctuating wind field terms achieve further dispersion of the
particles by finer scale flow features. The application of the
mean wind field is straightforward because all mesoscale met models
output these wind components (at typical horizontal grid spacing
of 5 to 50 km, and with temporal data at hourly intervals). It
must be realized, however, that in grid models, energy is resolved
on scales of two grid spacings. Furthermore, because of filtering
and numerical diffusion, the minimum resolution is, in practice,
closer to four grid spacings. For a model using a horizontal grid
spacing of 25 km, the minimum wavelength of the energy spectrum
which is resolved is on the order of 100 km. Below the resolvable
scale, the met models generally represent the subgrid part of
the energy spectrum by employing turbulent closure to parameterize
boundary-layer turbulence which has a scale length comparable
to the height of the PBL(~1 km). Thus, in grid models, atmospheric
energy of wavelength between about 1 km and about 100 km is not
represented (McNider, 1989). This universal shortcoming of grid
models may be expected to lead to a general tendency toward underestimation
of plume dispersion as simulated by LP models using input winds
generated by a mesoscale met model. In Gaussian models of point-source
plume dispersion, resolved scale motion is generally incorporated
through hourly wind inputs, and unresolved energy is parameterized
in terms of Pasquill-Gifford-Turner (PGT) dispersion coefficients.
The representation of the unresolved energy is based on ideal
dispersion concepts that may not represent the local micrometeorology
or account for shear effects. Furthermore, a gap may exist between
the interval of the energy spectrum represented by the hourly
data taken from towers or routine NWS measurements (surface or
upper-air) and the interval represented by the PGT curves. If
so, then the energy in the gap is not being taken into account
in the regulatory models, and its effect on dispersion is ignored.
The accuracy of each of the two terms (grid-resolved and subgrid-scale)
in Eq. 2.1 depends on the nature of the met model and its outputs,
as well as on the time step of the particle model. The accuracy
of the resolved (deterministic) terms depends mainly on the spatial-temporal
resolution of the met model, on model assimilation of measured
data, on model treatment of a variety of physical processes (including
boundary layer and surface energetics, cloud and precipitation
processes and energetics, mesoscale complex systems, etc.), and
on the numerical advection scheme used in the model. The accuracy,
or lack thereof, of the subgrid term depends mainly on the size
of the unresolved energy gap, on the turbulence closure scheme
used in the model, and on subgrid-scale inhomogeneities (e.g.,
surface characteristics) which determine the amount and intensity
of subgrid-scale fluctuations. In the LES approach ("Large
Eddy Simulation"), the computational mesh is fine enough
to resolve all motions down to the scale of the large energy-containing
eddies of microscale turbulence, and the relatively small effect
of finer scale eddies is parameterized. The LES model, however,
is computationally extremely demanding and its applications are
generally limited to much smaller domains than those of mesoscale
met models.
A substantial energy gap remains unresolved and unparameterized
in current mesoscale and regional models. Research is under way
to develop parameterization of the dispersion due to energy in
the gap. In view of this situation, it is useful to separate ui'
into two components, thus: ui' ~ ui'
, gap + u' , sub-gap.
In the remainder of the report, only the sub-gap component is
addressed, and is generally referred to as ui'.
The need for future inclusion of ui' , gap
must be kept in mind, however. In the meanwhile, some facility
to calibrate simulations without its inclusion by observations
must be considered.
In the current version of C-MC, the subgrid (sub-gap, stochastic)
terms in Eq. 2.1 are not treated explicitly. Rather, the deterministic
motion of the particles is perturbed in two ways to introduce
the effect of stochastic processes indirectly : at every advection
time step (Dt ~ 2 hours), the height
of each particle is changed by the application of a random vertical
displacement, with the constraint that no crossover is permitted
across the mixing layer boundary; also at every Dt,
a random horizontal displacement, of magnitude confined to a length
scale proportional to a prescribed constant mesoscale eddy diffusivity,
is applied to each particle. A revision of C-MC is presently under
consideration in which the stochastic dispersion (specification
and application of u', v',w') would be based more explicitly on
fundamental principles.
In most Lagrangian particle models, the stochastic component is
commonly modeled as a first-order autoregressive or Markov process
using linear Langevin stochastic differential equations (e.g.,
Gifford, 1982; Sawford, 1984) of the form:
= - + , (i = 1, 2, 3) , (2.2)
where, ui' denote the three unresolved components
of the flow field, TLi denote the corresponding
Lagrangian time scales, and ui" denote
purely random (uncorrelated) velocity fluctuations arising from
white noise modeling of particle accelerations.
The above model is based on independent Langevin equations,
i.e. the motion of every particle is assumed to be independent
of the motion of all other particles. This leads to ensemble average
distributions of particles rather than specific distribution realizations
in time. Conse-quently, the system does not possess the ability
to observe such coherent motions as plume meander and other characteristics
of instantaneous snapshot realizations of particle distributions
(Yamartino, 1993). Such instantaneous realizations can only be
simulated with explicit representation of the 3-D field of turbulent
eddies of various sizes. Besides explicit dynamical simulations
(e.g., LES) which are too demanding computationally, an alternate
and computationally much more economical approach based on kinematic
simulation (KS) and currently still under development has recently
been claimed to have the potential to provide realistic representation
of eddy spectra (Yamartino, 1993). This approach consists of a
Fourier expansion of the turbulent flow into a spectrum of eddies
each satisfying conservation of mass. In it, the amplitudes and
phases of the various Fourier components of turbulence evolve
according to Langevin equations. This evolutionary nature of the
turbulence fields "causes the final spectrum of velocity
fluctuations to smear out into surprisingly realistic shapes.
In addition, due to the fact that the particles now follow the
local turbulent flow fields, rather than move independently, one
observes coherent motions, such as plume meander, and other characteristics
of instantaneous snapshot realizations of particle distributions
rather than simple ensemble averages." (Yamartino, 1993)
Zanetti (1984) proposed the following form of the general solution:
u'(t) = f1 .
u'(t-Dt) + u" , (2.3a)
v'(t) = f2 .
v'(t-Dt) + f3
. u'(t) + v" , (2.3b)
w'(t) = f4 .
w'(t-Dt) + f5
. v'(t) + f6 .
u'(t) + w" , (2.3c)
where fi are unknown
coefficients to be determined. The first terms on the right hand
side of the above equations are the Markov (memory) terms which
represent Lagrangian autocorrela-tions of the wind velocity fluctuations
for lag time Dt; the last terms are
the purely random, uncorrelated terms; the middle terms represent
cross-correlations between the velocity fluctuations. The method
of determining the unknown coefficients is given in Uliasz, 1992.
Uliasz also discusses five special cases of the above model with
varying degrees of complexity depending on whether the Markov
terms, the cross-correlation terms, and the horizontal turbulent
diffusion terms (u", v") are retained. Other applications
have made varying approximations also: for example, McNider (1981)
neglected all three cross-correlation terms; Boybeyi (1993) neglected
two of them (f3
= f5 Å 0);
C-MC effectively neglects all six f-terms,
and also does not treat turbulent dispersion explicitly as in
Eq. 2.3.
If the cross-correlations between the velocity fluctuations are
neglected in Eq. 2.3, then f3
= f5 = f6
Å 0. This simplification, however, leads to local flows
which are not constrained by mass conservation (flow non-divergence),
and results in a tendency for particles to accumulate in low turbulence
zones. To correct for this, it is necessary to introduce "drift
velocity" corrections (ud, vd,
wd) in equations 2.3a, 2.3b, 2.3c, respectively.
Most commonly, only wd is retained (ud
~ vd ~ 0) since turbulence inhomogeneity is
most significant in the vertical, and it may be expressed as (after
Sawford, 1985) :
wd Å ( 1 - Rw )
TLw . . . [ 1+ ()2
] , (2.4a)
where, Rw = Rw(Dt)
_ Å exp (-
) is the Lagrangian autocorrelation
function of w', and sw2
is the variance of w'. Legg and Raupach (1982) and Ley
and Thomp-son (1983) have shown that this expression can be approximated
effectively by
wd Å ( 1 - Rw )
TLw . , (2.4b)
and, for coarse Dt, simply as
wd Å Dt
. . (2.4c)
When the cross-correlation terms are neglected, it can be shown
(e.g., Uliasz, 1992) that the coefficients of the Markov terms
reduce to
f1 = Ru
, f2 = Rv
, and f4 = Rw
, (2.5)
and the solution of Eq. 2 becomes (with the wd
correction):
u'(t) = Ru . u'(t-Dt)
+ u" , (2.6a)
v'(t) = Rv . v'(t-Dt)
+ v" , (2.6b)
w'(t) = Rw . w'(t-Dt)
+ wd + w" , (2.6c)
Bearing the definition of Rui in mind (see following Eq. 2.4a), the Markov (autocorrelation) terms become negligible when Dt >> TLi ( TLw < TLu ~ TLv ~O[1 min] ). If the present time step of ~2 h is maintained in the application of the revised C-MC, then the Markov terms will be negligible, reducing the model to a pure random walk model (Monte Carlo) in which the current particle motion is devoid of any memory of past history. For small Dt (TLi ), the Markov terms are significant and must be retained. The corresponding homogeneous Markovian turbulence is, as will be shown, completely characterized by 3-D fields of velocity fluctuations of variances ( su2 , sv2 , sw2 ) and corresponding integral time scales (TLu, TLv, and TLw ). Such a representation still ignores the fact that real turbulence consists of a wide spectrum of eddy sizes, and that the longer wavelength eddies have longer characteristic time scales than the shorter wavelength components. One consequence of this is that near boundaries (e.g., the ground), which act as spectral filters, the modeler must conjure up a time scale profile (e.g., TLw (z) ) that represents a best compromise time scale (rather than going to a more fundamental formulation based on eddy size). Another practical matter to bear in mind is that, because the scale of turbulence is smaller for the vertical component, it is often more efficient to use time-splitting, with the use of a shorter Dt for vertical turbulent dispersion. Typically, time steps of 10-30 s are used for the vertical terms and 1-5 min are used for the horizontal terms in a Markov simulation.
In a homogeneous symmetric turbulence model, the random terms
(u", v" and w") are characterized by = == 0, and
variances ( su"2
, sv"2 ,
sw"2 ). They
can be estimated at each time step as
u" = su"
. hu , v"
= sv" . hv
, w" = sw"
. hw , (2.7)
where hu , hv
, hw are
random numbers from a standard Gaussian distribution with mean
= 0 and variance =1. With the cross-correlation terms neglected,
the variances can be shown to be given by
s= ( 1 - R) . s.
, (2.8)
and the model becomes
u'(t) = Ru . u'(t-Dt)
+ ( 1 - R)1/2 . su
. hu , (2.9a)
v'(t) = Rv . v'(t-Dt)
+ ( 1 - R)1/2 . sv
. hv , (2.9b)
w'(t) = Rw . w'(t-Dt) + wd + ( 1 - R)1/2 . sw . hw , (2.9c)
The problem now reduces to one of quantification of the variances
su2 , sv2
, sw2 and
the corresponding time scales TLu, TLv,
and TLw (which are needed to evaluate the
R's and wd). Since the LP model is to be used
in conjunction with the flow fields predicted by the met model,
the turbulent statistics must be deduced so as to be consistent
with the dynamically evolving boundary-layer flow. The determination
of these turbulence statistics will depend on the turbulence closure
scheme of the mesoscale met model and the corresponding outputs.
Evaluation of these turbulence statistics is the most difficult
part of Lagrangian particle modeling. Ideally, such evaluation
must be based on suitable Lagrangian measurements of the flow
field. Unfortunately, Lagrangian properties are difficult to measure
and, therefore, must be inferred from Eulerian measurements. Different
methods have been proposed for relating Eulerian and Lagrangian
statistics, but none of them has been found to be totally satisfactory.
Uliasz (1992) uses probably the most advanced method based on
a 21/2 order closure scheme involving solution
of prognostic equations which include explicit simulation of the
turbulent kinetic energy (TKE). Boybeyi (1993) also used a mesoscale
model with explicit TKE simulation, but he follows the approach
of Hanna(1982) who proposed a set of semiempirical formulae based
on a limited number of meteorological similarity parameters (zi,
z0, L, u*, w*; zi = mixing
layer height; z0 = roughness height; L = Monin-Obukhov
length; u* = friction scaling velocity; w* = convective scaling
velocity) to estimate the turbulence statistics at each elevation
(z). In this approach, all the parameters necessary to calculate
the turbulence statistics are provided by the PBL model assuming
that surface layer similarity parameters are valid for the entire
mixing layer. For heights greater than the mixing height, Boybeyi
uses expressions for sui
proposed by Ettling et al. (1986) which include values of the
mean TKE, and he uses constant values of 30 s for TLw,
and 600 s for TLu and TLv
.
McNider (1981) proposed a scheme for the parameterization of sui
and TLi using the turbulence length scales
(lui
). This scheme is outlined below, because it is based entirely
on variables outputted by both mesoscale met models considered
here, and on first order closure schemes.
Unstable Conditions (Daytime CBL) --------------------------------------------
(Eqns. 2.10)
Based on Kaimal et al. (1976), lu
Å lv Å
1.5 zi and
lw Å
Then, the Lagrangian time scales (TLi) are
derived from the Eulerian time scales (TEi)
based on the following relationships:
TLi = b TEi
, where TEi = 0.2 ( lui
/ ) and b = 0.6 ( / sw
) ; and,
the standard deviations of the velocity components (sui
) are given by :
su = sv
= u* [ 12 + 0.5 (zi/|L|) ]1/3
and
sw = , where Kz
is the vertical exchange coefficient and
A = 
.
If the mesoscale model is based on first order closure of turbulence,
Kz could be available in its output. In case
it must be computed, several models of Kz
are available for the convective boundary layer. The profile used
in RADM for vertical mass transfer in the CBL combined the standard
similarity theory-based profile in the surface layer ( = ku*z/f(z/L)
) with Kz(z) = kw*z (1 - z/zi)
in the bulk mixing layer above. Such a choice is reasonable if
the K-approach is chosen. However, the K-approach assumes turbulence
to be vertically symmetric in the CBL. This is not the case in
fact, with non-similar rapidly rising buoyant plumes and relatively
gradual compensating subsidence. Unlike eddy diffusion which only
permits exchange between adjacent layers in both directions, the
real updrafts in the CBL permit rapid upward exchange between
the surface and all layers above. In this sense, upward exchange
is non-local. Stull (1988) has described various non-local closure
approaches. Pleim and Chang (1992) have developed an asymmetric
non-local closure scheme which is used in RADM, and may be incorporated
into MM5 in the future. Recently, Hurley and Physick (1993) have
incorporated the asymmetry concept in a Lagrangian particle model
through use of a skewed probability density function representing
turbulence. This pdf is made up of two Gaussian distributions,
one representing updrafts and the other representing downdrafts.
These distributions are weighted based on the probability (p)
of the particle being in an updraft, and are characterized by
different mean and variance values of the particle velocity for
updrafts and for downdrafts. The determination of these turbulence
statistics requires not only the mean (= 0) and the variance (2
= s) of the vertical velocity distribution,
but also its third moment, 3 = S, which characterizes
the skewness of the pdf.
Neutral and Stable Conditions --------------------------------------------
(Eqns. 2.11)
su = sv
= 2.3 u* and sw
= 1.2 
[ ]0.58
| | , where
dV/dz is the local vertical wind shear and 
is a mixing length, simply taken as a linear function of z near
the surface and a constant value of 70 m (after Blackadar, 1979)
elsewhere, as follows :
= 0.35 z , z < 200 m
= 70 m , z ³ 200 m .
Ri is the local gradient Richardson number given by, Ri _ Å ,
and the critical Richardson number, RIcr _
0.115 (100 DZ)0.175 ,
where DZ denotes the local model grid
spacing over which Ri is evaluated.
The Lagrangian time scales (TLi) are derived
from the Eulerian time scales (TEi) as previously
described. The length scales are given by
lu Å lv
Å 0.7 (z/zi)1/2
. zi = 0.7 (z . zi)1/2
Simulations performed for complex terrain with different versions
of the Lagrangian particle model of linked to various met
models (Uliasz and Pielke, 1993) have indicated that in mesoscale
applications, horizontal turbulent diffusion of pollutants is
relatively unimportant and may be neglected. There is also no
need for more refined particle models (e.g., models including
wind covariances). A more serious simplification consists of replacement
of the Markov chain scheme with the random-walk scheme which can
be run with a much longer time step. The random walk particle
models predict the same concentration patterns as the more advanced
particle models, although some differences in maximum concentrations
are evident. These simplified particle models are a very attractive
tool for mesoscale studies since they can be an order of magnitude
faster than the Markov chain models. It should be noted, however,
that the random walk particle models with long time steps may
not be suitable for use with some physical parameterizations,
e.g., plume rise, dry deposition (Uliasz, 1994).
Based on the above observations, it would appear that the present
version of C-MC is probably not bad for applications to long range
transport. The following recommendations, however, are made for
its possible improvement:
1. Consideration should be given to reducing the time step from
2 h to 1 h or less. Such a time step will still make it unnecessary
to include the autocorrelation terms. While the mesoscale models
output the mean wind components () at hourly intervals, these
variables are calculated typically at five minute intervals (the
outputs are NOT one-hour averages). Thus, the available met data
from the mesoscale met models will be much better resolved than
2 hours and there is no good reason (related to input information)
to not use them at a significantly finer resolution than 2 h.
2. A more defensible approach (compared to the current C-MC approach)
needs to be used to determine the very important random vertical
diffusion (e.g., based on a length scale, lw
, calculated as done in the McNider model presented here);
3. Consideration should be given to include a vertical drift velocity
correction term, wd , while the cross-correlation
terms may continue to be neglected. This correction may be expected
to improve the simulation significantly for conditions corresponding
to strong vertical turbulence inhomogeneity.
4. A mesoscale horizontal "eddy" diffusion term (or
some other mechanism to augment horizontal dispersion) should
be retained in the model, to simulate the effect not so much of
turbulent diffusion, but of the neglected "energy
gap". I suspect strongly that application of C-MC or any
other current particle model, with the wind fields generated
by a mesoscale met model , will underestimate
the horizontal dispersion of a point-source plume during mesoscale
or regional transport (even if the model includes simulation of
the microscale turbulent diffusion terms). Of course, in the absence
of any real knowledge of the effect of the missing energy gap,
there is currently no objective way to represent such effect accurately.
This is an area where further basic research is needed.
Based on the McNider particle model (Eqns. 2.1, 2.9, 2.4b, 2.5,
2.10, 2.11), the required met variables at each time step are
as follows:
3-D (x,y,z) : , , ;
2-D (x,y) : w*, u*, L , H0 , zi
, z0 (H0 = sensible heat
flux)
The sensible heat flux is included in the above list (although
it is not specifically needed if u* and w* are available), because
it is viewed as a useful diagnostic quantity to have access to.
If the TKE solution is available, it is preferable to characterize
the turbulent dispersion terms based on TKE (as in the use Uliasz,
1992, for example) rather than on surface layer similarity parameters,
particularly for dispersion outside the PBL.
An important objective of the present project is the ability to
implement the revised C-MC with available met outputs from state-of-the-art
mesoscale met models. Two mesoscale met models stand out in this
respect: the Regional Atmospheric Modeling System (RAMS; Pielke,
et al., 1992) developed by Colorado State University, and the
Mesoscale Models MM4 ( Anthes et al., 1987) and MM5 (Dudhia, 1993;
Grell et al., 1994) developed jointly by Pennsylvania State University
and the National Center for Atmospheric Research. In this project,
these two model sets were examined for their utility as suppliers
of inputs to the upgraded C-MC. The two model sets are very briefly
overviewed here specifically for their utility in the present
context.
Both modeling systems (RAMS and MMx) are dynamic 3-D Eulerian
models with prognostic equations for pressure, the mean wind components,
temperature and moisture. Both systems use comparable schemes
for model initialization (based on ECMWF data), and for four-dimensional
data assimilation (FDDA) based on data of the NWS network of upper-air
soundings and the surface network. RAMS and MM5, in particular,
are quite comparable as preprocessors for a long-range LP model.
The physics treated in the two models is at a comparable level
of depth and breadth for our purpose. Both have the nonhydrostatic
option which permits almost arbitrarily fine horizontal spatial
resolution, as well as 2-way multiscale grid nesting capabilities
to permit a combination of fine spatial resolution and computational
efficiency. MM4 is limited to hydrostatic assumptions which limit
the horizontal spatial resolution to about 10 km or coarser, in
general. The nesting capabilities of MM4 are also much less versatile.
Almost all available outputs from MM4 are at horizontal resolution
of 20 km or coarser (mostly at 80 km). With MM5 replacing it now,
it is unlikely that there will be many new runs of MM4. However,
existing available outputs from MM4 runs over the US significantly
exceed those from RAMS and/or MM5 runs. In particular, EPA has
run MM4 at 80 km horizontal resolution (15 vertical layers) for
the entire year 1990, and the outputs are readily available. Currently,
these outputs do NOT include all the needed variables as listed
in Section 2.4 (particularly most of the 2-D outputs listed).
EPA, however, has plans to either repeat the runs for the entire
year or to derive the missing outputs from the currently available
outputs for the full year. This is expected to be done within
the next few months. Thus, a full year's inputs to drive runs
of the revised C-MC may be counted on to be available at 80 km
and 1 hour resolution. In addition, some MM4 runs at 80 km (a
subset also at 27 km) are available for selected short case study
periods corresponding to the FIFE measurement program of 1987
and the ACIDMODES measurement program of 1988.
RAMS has been around for a long time and has also been very flexible
in terms of options. Thus, it has been implemented for a wide
variety of scenarios. Since the addition of the FDDA capability
in recent years, its outputs including this option are potentially
useful for C-MC implementation. Also, some RAMS runs have included
explicit simulation of TKE. However, the user-friendly version
of RAMS is a proprietary system and has not been used much by
EPA. RAMS outputs are not easily available. MM5 is a very new
tool, by comparison. It is, however, a public-domain tool and
is expected to be used extensively by EPA and others in the coming
years. It will be the standard met model for the Models3 modeling
system of EPA. Already, EPA has the outputs from the following
nonhydrostatic runs of MM5 available for 1987 and 1988:
Period Region Spatial Resolution Length
of Simulation
Jul 8-11, 1987 C. USA 54 / 18 km; 30 layers 96 h
Aug 3-5, 1988 C. USA 54 / 18 km; 30 layers 48 h
Aug 3-5, 1988 E. USA 54 / 18 km; 30 layers 48 h
Aug 3-5, 1988 W. USA 15 / 5 km; 30 layers 48 h
Aug 28-30, 1988 E. USA 54 / 18 km; 30 layers 48 h
Aug 28-30, 1988 E. USA 36/12/4 km; 30 layers 48 h
Aug 28-Sep 2, 1988 E. USA 54 / 18 km; 30 layers 120 h
In addition, EPA is also planning to run MM5 at 36/12/4 km horizontal
resolution for the entire period of the recent SOS Nashville Field
Study (~ 5-6 weeks in June-July 1995) after the input data become
available early in 1996. The outputs of all MM5 runs will have
all the variables likely to be needed for the revised C-MC runs.
The TKE option is also expected to become available in MM5 in the very near future.
Emissions accompanied by significant thermal and/or momentum effluxes
are subject to significant plume rise. For power plants and moderate-to-large
industrial "point" sources, the major contribution to
plume rise is from the heat flux and buoyancy effects rather than
from momentum effects. Whenever the wind speed, Uï
, is less than 10 m/s in the lower troposphere (which is most
of the time), plume rise, Dh,
is greater than the physical stack height, hs
. In lighter winds and neutral or unstable ambient
conditions, it can be much greater. Since the effective stack
height, heff (= hs+Dh),
has strong influence on pollutant atmospheric residence time and
hence on plume transport, dilution, chemistry, removal, and ground-level
impact, plume rise is very important, especially for large elevated
point sources.
Plume rise depends on physical and emission variables at stack
exit as well as on the condition of the ambient atmosphere (the
latter depending on dynamic as well as surface characteristics).
There are a large number of plume rise formulas out there, and
results of plume rise calculations based on these generally do
not agree well with each other or with observations. Such calculations
often give results differing by as much as a factor of ten. Much
of the confusion is due to effectively different definitions of
plume rise. The term "final" plume rise is much in use
in practice. It presumes that the plume reaches some specific
maximum height during short range transport (a few kilometers).
Such a presumption is valid only under stable ambient conditions,
in which a plume tends to level off and subsequently maintains
nearly constant elevation for many kilometers downwind. In neutral
conditions, plumes have never been observed to level off. In unstable
conditions, plumes dilute rapidly and exhibit looping behavior,
"swinging up and down enough to throw a maddening degree
of scatter in the comparison with any model, no matter how elaborate"
(Briggs, 1975). In neutral and unstable conditions, plume rise
effectively terminates if the plume breaks up due to the action
of ambient turbulence. In convective conditions, plume rise may
also terminate as a result of plume touchdown when the plume is
caught in a downdraft. Recently, there have also been attempts
to develop Lagrangian particle transport formulations for the
plume rise stage. In this section, both approaches are outlined.
The focus in this section will be on emissions from large point
sources for which buoyancy effects are generally the most important,
but momentum effects can also be significant, particularly in
the early stages of plume rise. Also, it is presumed that plume
rise calculations will be based on meteorological variables supplied
by outputs of mesoscale met models such as MM5, which typically
have rather coarse horizontal spatial resolution and quite fine
vertical resolution, and a temporal resolution of 1 h.
This section is essentially an overview (not at all an exhaustive
review) of the state of the art. It is based mainly on three well-known
review articles (Briggs, 1975; 1984; Weil, 1988) which are collectively
quite encompassing of the state of the art. Immediately following
release from the stack, plume rise is typically dominated by momentum
effects for a short time. Subsequently, plume rise and growth
near the source are largely due to entrainment of ambient air
as a result of buoyancy-generated turbulence within the
plume. The strength of this self-generated turbulence in the plume
decreases as the plume rises, and ultimately becomes fully dissipated.
In a stable atmosphere, plume rise ceases at that time. In a neutral
or unstable atmosphere, if the ambient turbulence is strong enough,
then its effect may, at some point during plume rise, become dominant
over the effect of the decaying buoyancy-generated turbulence.
This may cause the plume to break up, thereby terminating further
plume rise. Under stable conditions, theory based on the assumption
that the role of ambient turbulence is negligible generally yields
results which are consistent with laboratory and field measurements.
Such theory is outlined first. Its applicability to neutral conditions
is examined next, and the effect of ambient turbulence is then
explored for neutral and convective conditions. After that, the
effect of the elevated inversion layer capping the daytime CBL,
and the methodology to account for penetration of this layer,
are outlined.
The basic theory of plume rise (based on the so-called Brigg's
formulas) is limited to assumptions of uniform wind field and
constant stability in the atmospheric layers through which the
plume rises. Various choices of these constant values of Uï
and the stability parameter have been used in different applications,
e.g. : local values at stack height or at some multiple of stack
height (e.g., 1.5*hs);
or, average values for the layer between the physical stack height
and a guesstimated effective stack height. In cases where met
inputs at high vertical spatial resolution are available from
mesoscale met model runs or from field measurements, the effect
of wind shear and vertical variations of atmospheric stability
can, in principle, be incorporated in plume rise models which
possess such spatial resolution. Turner (1985) suggested a method
which does precisely that. It uses a layer-by-layer approach with
decreasing "residual buoyancy", and variable wind speed
and stability with height. Turner's scheme is attractive for its
simplicity and superior use of available vertical information.
It has been applied in the multi-plume model TUPOS (Turner et
al., 1986) and in the regional model RADM (Byun and Binkowski,
1991). Since Turner's scheme automatically incorporates the effect
of local stability on buoyancy and plume rise, it does not require
separate special treatment of plume penetration of the elevated
inversion layer.
The plume also spreads during plume rise by the entrainment of ambient air, at first mainly due to buoyancy-generated turbulence, but progressively more by the effect of ambient turbulence. In most applications, the plume is initiated as a point source at z = heff and x = 0 (see Fig. 1); thereafter, its spread is treated as non-buoyant and due to ambient turbulence only (in some cases, also by wind shear with height). In other applications, initial plume spread at final plume rise is assumed to be equal to plume rise (e.g., Turner et al, 1986; Byun and Binkowski, 1991) or some fraction of it. Such a plume may then be initialized at z = heff and downwind location x given by the well-known "2/3-law" (Dh ~ x2/3, see Eq.14 below). In what follows, an attempt is made, as much as possible, to provide formulas for plume centerline trajectory as well as spread during the course of plume rise.
As mentioned earlier, this is the only case in which the plume
attains a "final" maximum rise over a short range and
plume centerline then remains at approximately that height during
subsequent stable transport. Also, under stable conditions, the
approximation that the effect of ambient turbulence is negligible
has been shown to be good, and is made here.
Figure 1 illustrates relevant physical variables for plume rise in a stable atmosphere for two scenarios:
a) vertical plume rise under calm ambient wind conditions (Uï < 1 m/s), and
b) bent-over plume for higher wind speed conditions.
(a) Vertical Plume (b) Bent-over Plume
For each plume, w denotes the vertical plume rise speed; for the bent-over plume, Us (= ds/dt) denotes the speed of plume material along plume axis (Us = w for the vertical plume). D and r are the local cross-sectional diameter and radius of the plume. Let (r , p, T) and (rï, pï, Tï), respectively, denote the (density , pressure, temperature) in the plume (assuming a "top-hat" profile) and in the surrounding ambient air. The buoyancy flux is given by
FB = = , (3.1)
and at stack exit (subscript "0") by
(FB)0 = = = ,
where QH and 0 (= p r02. w0 ) are, respectively, the stack heat and volume flow rates.
Implicit in the above equations are two assumptions:
i) p Å pï , which, when combined with the ideal gas relation p = r R T, gives r a (1/T); and,
ii) (rï - r0)/r0
<< 1 in the inertial terms, but not in the buoyancy term
(Boussinesq approx.).
Governing Basic Equations
For the shaded plume element shown (between s and s+Ds),
the governing equations of mass, momentum and energy balance may
be written as:
Mass Balance : = E = 2pr . vE
, where (3.2)
E is the rate of mass entrainment from plume sides, and vE
(= dr/dt) is the entrainment velocity. Based on observed linear
growth of plume diameter with plume rise, it is assumed that
vE Å b
w . (Taylor Entrainment Assumption) (3.2')
The proportionality constant b is generally
referred to as the entrainment factor. Eq. 3.2' represents a first-order
closure of the effect of the plume's self-generated turbulence.
With the assumption that the effect of ambient turbulence on momentum
and energy balance is negligible, we have the following equations
for momentum and energy balance:
Horiz. Momentum : = - p
r2 w , (Du
= u - Uï) (3.3)
Vertical Momentum : = p r2
g , (3.4)
Energy Balance : = - p
r2 w , (DT
= T - Tï)
(3.5)
The energy equation is also the conservation of buoyancy equation,
with dq
/dt = 0 (i.e., the plume is assumed to rise adiabatically;
q = potential temperature). Combining
the definition of FB (Eq. 3.1) with Eq. 3.5
and the relations r a (1/T) and qï
Å Tï yields
= - r2 w N2
, (3.6)
where N2 = is a stability parameter, and N is the Brunt-Vasala frequency.
For plumes, N-1 is the time scale of buoyancy
depletion and of maximum plume rise in a stable atmosphere ( N-1
~ O[1 min] ).
In addition to the above, the kinematic condition is also satisfied
by the relation
= (3.7)
The solution of the above set of equations for a stable thermal
environment (N2 > 0), with (calm
wind conditions (Us Å w , Uï
Å 0 Å _Uï/_z) yields
a "final" plume rise given by
(stable, vertical plume) (3.8)
Briggs (1969) and Weil (1982) have found good correlation of this
formula with field observations covering a broad range of heat
sources under calm, stable conditions.
, the solution for plume trajectory during rise is obtained, in
terms of initial momentum and buoyancy fluxes [(FM)0,
(FB)0], as (see, for example,
Weil, 1988)
= ()1/3 [N'
.( (FM)0/(FB)0
) . sin(N'x/Uï)
+ 1 - cos(N'x/Uï)]1/3
, (3.9)
and the plume spread rate is given by
r = r0 + b' z' , (3.10)
where,
b' = , N' = , (FM)0 = w02 r02 , (FB)0 = w0 r02 ,
and kv is a factor which accounts for the displaced mass of ambient air as the plume rises;
kv Å 1 for a circular plume cross-section.
The entrainment factor (b) and the related plume growth factor (b' ) play a very important role in plume rise determination. The following empirical assumptions are generally accepted:
i) the entrainment velocity vE a w , or equivalently, that plume growth rate r a z; and that
ii) the entrainment factor is much larger for bent-over plumes than for vertical plumes.
Briggs (1975) recommends use of b = bv Å 0.125 for the vertical plume, and
b Å 0.6 for the bent-over plume.
Assuming kv = 1, gives b' = b'v Å 0.09 for the vertical plume, and
b' Å 0.42 for the bent-over plume.
Based on an extensive survey of field and laboratory observations, Briggs (1975; 1984) concluded that the best formula for the "final" plume rise in stable conditions with Uï > 1 m/s is
(3.11)
We consider now the plume rise in a neutral atmosphere with a
uniform flow field. For now, we continue to ignore the effect
of ambient turbulence in the entrainment process. Then, equations
(3.1) through (3.7) continue to apply.
For neutral stability, N2 = 0. From Eq. 3.6
then, FB = constant = (FB)0
, i.e., the buoyancy is conserved and the plume rise, in principle,
can continue indefinitely, unless the conditions change downwind
(as they do, at least because of the effect of ambient turbulence).
Near the source, the plume trajectory can follow three different
scaling laws, each valid in a different x-regime (Weil, 1988).
The actual trajectory depends on wind speed (Uï)
and on source momentum and buoyancy fluxes [(FM)0,
(FB)0]. The three x-regimes
are delineated by two "critical" x values, xc1
and xc2, which were derived
by Hoult and Weil (1972) as :
= ()2 ()3
and = 2()2 , with
M and 
B,
respectively, denoting the momentum and buoyancy length scales,
and R denoting a velocity ratio. These three quantities are defined
as follows:
M _
(r0/rï)1/2
R r 0 , 
B
_ (FB)0
/ Uand R _
w0 / Uï
.
R is an important parameter determining when plume downwash can occur in the lee of a stack. For most stacks (excluding cooling towers), downwash is negligible for R > 1.5 (Briggs, 1984).
The three segmented solutions are outlined below (Weil, 1988):
= ()1/2 ()1/2
x < xc1 ; (3.12)
= ()1/3 ()1/3
xc1 < x <
xc2 ; (3.13)
= ()1/3 ()2/3
x > xc2 . (3.14)
Eq. 3.12 was derived by Hoult et al. (1969) for a plume dominated
by mass and momentum fluxes very near the source where plume axis
departs very little from the vertical. Eq. 3.13 was derived by
Hoult and Weil (1972) who used b =
0.6. Briggs (1984) used b = 0.4 + 1.2/R.
The trajectory represented by Eq. 3.14 applies in the bent-over
stage of a buoyant plume where both the buoyant flux and the entrainment
rate are important. It describes the well-known "2/3 law"
(Slawson and Csanady, 1967; Briggs, 1975).
Of all the simple plume rise formulas, probably the best documented
is the "2/3 - law". It has been shown to be in agreement
with many laboratory and field observations with b
= 0.6. This fact (same b for both lab
measurements with no turbulence and field measurements, under
near-neutral conditions) indicates conditions for which the assumption
of neglecting the effects of ambient turbulence is valid. This,
however, is not so when ambient turbulence is strong, as in the
unstable CBL.
Gaussian models generally break up a plume into two stages : a
buoyant plume rise stage in which only self (buoyancy)-generated
turbulence is important, and a non-buoyant stage in which only
ambient turbulence is important. Observations of asymptotic plume
rise and looping plume behavior suggest the involvement of ambient
turbulence (small-scale and large-scale, respectively). There
are two main ways in which ambient turbulence can affect plumes.
First, if the small-scale turbulence (with eddy size r) is sufficiently
intense, it can increase the plume growth rate beyond that given
by buoyancy-induced turbulence and lead to an asymptotic rise.
Second, if the large-scale turbulence (eddy size > r) is strong,
it can carry entire sections of the plume up and down, thereby
dispersing the plume by meandering. Whether one or both of these
processes occur is determined by the turbulence structure of the
boundary layer. In the CBL, measurements (e.g., Willis and Deardorff,
1974) and numerical modeling (e.g., Deardorff, 1972) show that
the turbulence is dominated by large-scale elements consisting
of thermal plumes (updrafts and downdrafts). These elements extend
from the surface to the CBL top, are long-lived, and have length
and velocity scales proportional to zi and
w*, respectively.
Briggs has advanced two schemes which separately include the effects
of small-scale and large-scale turbulence on plume rise. They
are, respectively, the plume break-up model (Briggs, 1975) and
the plume touchdown model (Briggs, 1984).
The breakup model assumes that the plume maintains its integrity as long as its decreasing buoyancy-generated turbulence remains greater than ambient turbulence. Thereafter, ambient turbulence in the inertial subrange would dominate the entrainment process and would break up the self-generated structure of the plume. The most important eddies would then be those of size comparable to r. vE would now depend on the ambient turbulence dissipation rate (eï) and on plume radius, as vE = (r eï)1/3 > bw. This represents a different closure assumption than the Taylor entrainment closure for buoyancy-generated turbulence. Based on dimensional analysis, Briggs assumed the turbulence dissipation rate prior to breakup to be given by
eB ~ h
w3/z', where h
is a constant. Based on Eq. 3.14 for z', and recognizing that
w = dz'/dt, he derived the following expression for z' as a function
of eB during the
buoyancy dominated stage:
z' = ()1/3
( )3/5 ()2/5.
Briggs defined an "effective" plume rise as the point
of breakup of the plume, where eï
= eB(z' = Dh),
and
Dh = ()1/3
( )3/5 ()2/5
. (3.15)
Briggs used this equation as the basis for deriving plume rise
formulas based on plume breakup for the two cases of ambient turbulence
characterized by small-scale turbulence (neutral and weakly unstable)
and by large-scale turbulence (strongly unstable).
Plume Rise Based on Plume Breakup in a Neutral Atmosphere
In the neutral case, eï
is given by u*3/kz
= u*3/k(hs
+ z'), where k is the von Karman constant
= 0.4. Substituting this expression in Eq. 3.15 and Dh
for z', we get an implicit equation for Dh
which may be solved iteratively. Briggs has suggested the following
explicit equation as a good approximation (Weil, 1988) :
, (3.16)
where, F*N _ and
u* _ ()1/3
.
H0 = sensible surface heat flux and L is the
Monin-Obukhov length scale.
Plume Rise Based on Plume Breakup in an Unstable Atmosphere
In the unstable case, eï
is given by 0.25w*3 /
zi . Substituting this expression in
Eq. 3.15 and Dh for z', we get (Briggs,
1984) :
, (3.17)
where, F*U _ and
w* _ ()1/3
Plume Rise Based on Plume Touchdown in an Unstable Atmosphere
Plume rise can also be terminated if the plume gets caught in a downdraft in the CBL. Weil and Hoult (1973) derived a plume rise formula for such a case based on the assumption that a strong enough downdraft could occur for sw' (of ambient turbulence) exceeded plume rise speed w. Csanady (1973) proposed a similar model using w = m sw', where m was an empirically chosen constant. Briggs (1975) carried this idea further by arguing that plume segments caught in downdrafts would eventually be carried to the surface when the rise velocity was less than the downdraft speed. For a downdraft speed of wD, he derived the following approximate formula for plume rise corresponding to the point when plume bottom just reached the surface (Briggs, 1984) :
DhD Å 2.9
()3/5 hs2/5 (3.18')
Briggs (1984) also showed that wD Å
0.4 w*. Substituting this in Eq. 18' gives
(3.18)
Eq. 3.17 gives a slightly lower value for Dh
than does Eq. 3.18 for the unstable case, and is the one recommended
by Briggs (1984).
During the daytime, emissions in the CBL may interact with the
elevated inversion layer capping the CBL. The inversion layer
(of thickness Dhi)
can be thin, being characterized by an inversion jump (Dqi),
or thick and characterized by the potential temperature gradient
(_qi/_z).
Depending on the strength of the inversion layer and the residual
buoyancy when the plume reaches it, the plume may penetrate it
partially or fully, or be fully trapped below it. There have been
various attempts to model partial penetrations. For bent-over
plumes, the models of Briggs (1975), based on buoyancy depletion
of the plume as it crosses the inversion, and of Manins (1979),
based on the plume density or temperature distribution when its
centerline reaches the inversion, represent the two main approaches.
Following a review of these, Weil (1988) recommends the following
combination of algorithms. For a thin inversion, if the
value of Dqi is
available, he suggests use of the model proposed by Manins (1979)
for the fraction (f) of the plume trapped by the inversion, as
follows :
f = - (P - 0.08) , < 0.4 (3.19)
where zi' _ zi - hs ,
P _ (FB)0 / (U_bi zi' 2) = a dimensionless buoyancy flux ,
bi _ g (Dqi/q_)
= Inversion layer stability parameter.
Plume rise is given by the equilibrium height (z) where FB
= 0. In terms of f,
z= zi' /(f + 0.5) (3.19')
For a thick inversion (Dhi/h'
> 0.4), if the value of _qi/_z
is available, he suggests use of the model of Berkovicz et al.
(1986), as follows :
f = 0 if zi' < 0.5 z ,
f = 1 if zi' > 1.5 z ,
f = - 0.5 if 0.5 z< zi'
< 1.5 z , (3.20)
where = [2.63 Ps + (2/3)3]1/3 , (3.20')
Ps _ (FB)0 / (U_Nzi' 3) = another dimensionless buoyancy flux ,
and N_ g (_qi/_z)/q_.
STABILITY Wind Ambient Turb. Effect Ambient Turb. Effect
Condition Negligible Not Negligible
Stable Calm Dhmax
Eq. 3.8
Uï > 1 m/s z'(x) Eq. 3.9
r(x) Eq. 3.10
Dhmax Eq. 3.11
Neutral Uï > 1 m/s z'(x) Eqs. 3.12, Dhmax Eq. 3.16 (Breakup)
3.13, 3.14
Unstable Uï > 1 m/s Dhmax Eq. 3.17 (Breakup)
Dhmax Eq. 3.18(Touchdown)
_________________________________________________________________
Partial plume penetration of Thin
elevated inversion layer Inversion Eq. 3.19, 3.19'
Thick
Inversion Eq. 3.20, 3.20'
Plume spread at the time of final plume rise (rF)
may be approximated based on rF Å b'
Dh, with b'
Å 0.09 and 0.42, respectively for vertical and bent-over
plumes.
Stack/exit variables : hs, r0,
T0, w0, r0
Meteorological variables : (3-D) , , ;
(2-D) w*, u*, L , H0 , zi , (H0 )
Briggs (1984) considered the case in which the meteorological variables are non-uniform in the vertical. The decrease in buoyancy flux with plume rise is given by the following modified form of Eq. 3.6 :
dF/dz = - N2 Uï
r2 = - N2
Uï (b'z')2
. Integration of this equation gives
,
which can be used to determine the residual flux on a layer-by-layer
basis during plume rise, with variable values of ambient wind
speed and stability from layer to layer. The plume rise is then
calculated using the residual buoyancy flux in the appropriate
formula. If the calculated effective stack height exceeds the
top of the current layer, then the calculations are repeated for
the next layer until the plume rise falls below the layer top.
Such an event signals that the final plume rise is in the current
layer. The final value of zgives the final plume rise. Plume spread
at this point can be estimated based on rF
Å b' Dh,
and the downwind distance from the stack (xF)
can be estimated based on the "2/3 law".
Turner (1985) applied such a layered approach. However, he did
not consider the effect of ambient turbulence and made other simplifying
assumptions which did not match the state of the art (Petersen
and Dowd, 1988). Byun and Binkowski (1991) included a more state-of-the-art
formulation in the RADM plume rise module. In this application,
only two plume rise formulas were used: Eq. 3.8 for the stable
case and Eq. 3.16 for the neutral breakup case. In the unstable
case (CBL), with one exception, the plume is assumed to spread
instantaneously over the whole CBL. The exception occurs when
the stack height is within 200 m of zi. In
this case, a partial penetration is calculated based on the Briggs
(1984) approach (depleting buoyancy model). With the exception
of the case of plume spread over the whole CBL, the plume spread
at final plume rise is taken as equal to the final plume rise.
The applications of Turner (1985) and of Byun and Binkowski (1991)
were not subjected to evaluations based on observations.
Most evaluations of plume rise formulas have been applied in the
context of the assumption of constant wind speed and stability
with height. Glendening et al. (1984) tested three approaches
(constant Uï and stability
at stack level; constant average values of Uï
and stability for the whole plume rise layer; and the layered
approach) against aircraft field measurements of plumes from a
power plant on the S. California coastline. He found the three
approaches to perform progressively better in the order specified.
Recently, Erbrink (1994) presented results of field evaluations
(lidar data) of the application of a more sophisticated layered
approach than Turner's and found model performance to be generally
satisfactory.
A layered plume rise formulation is recommended based on the plume
rise formulas summarized in Section 3.2.5.
The application of Lagrangian particle modeling to plume rise
is still in an embyonic stage. Two approaches have been proposed
and applied to varying extents, both involving modifications of
the equations involving the vertical transport of particles.
In one approach, initiated by Zannetti and Al-Madani (1984), an additional vertical velocity component (wb) representing the buoyant plume rise is added to Eq. 2.1c, thus:
z(t+Dt) Å z(t) + [+ wb
+ w'] Dt. In this manner, the effects
of buoyancy-generated and ambient turbulence can, in principle,
be considered together in an explicit manner. An expression for
wb is derived by time differencing any appropriate
plume rise formula for z' (a version of Eq. 3.14, the "2/3
law", was used by the above authors), i.e., wb
= [z'(t+Dt) - z'(t)] /
Dt. In this manner, the approach becomes
semi-empirical, but also subject to the errors of the plume rise
formulas. In principle, one can use different formulas according
to local stability and wind speed. The approach also, of course,
uses local height-varying values of wind speed and stability.
This approach was also applied by Anfossi et al. (1993), also
with the "2/3 law" formula only.
In the second approach, the Langevin equation in the z-direction
is itself modified to incorporate the effect of buoyancy. This
has been done in two ways. van Dop (1992) modified the Langevin
equation (see, for example, Eq. 2.2) for w' to the following form:
= - + B + , (3.21)
where B [_(g/T) (q
- qï)]
represents the rate of change of w' due to buoyancy-generated
turbulence. van Dop (1992) derived the following dynamic Langevin
equation for B :
= - - N2 w' + . (3.22)
TLw and TLB
are the appropriate Lagrangian time scales. van Dop (1992)
derived analytical solutions of the above formulation for the
case of purely-buoyant plume (i.e., with no momentum effect) and
assumed simple forms for the length scales (constant and proportional
to t). He demonstrated that his solution for stable calm conditions
corresponded asymptotically (t --> ï)
to the appropriate classical plume-rise formula. His model, however,
predicted plume spread rate to be faster than the plume rise rate
(b' > 1). His derivation of the
Lagrangian equations is also based on the work of Netterville
(1990) on plume entrainment, which invokes the mechanism of "ex
trainment" (negative entrainment, due to ambient turbulence).
van Dop's paper is a very interesting and original contribution
to Lagrangian particle modeling, but it is not immediately applicable
to practical plume rise problems.
Most recently, Hurley and Physick (1993) have adopted a different
form of the second approach (modified Langevin equation) for the
simulation of plume rise in the CBL. The authors start with the
more general form of the Langevin equation:
= a + (C0e)1/2
, (3.23)
where C0 is a universal constant, e is the rate of dissipation of turbulent kinetic energy, and "a" includes all effects, other than purely random, which impact the velocity fluctuations. For the special case of homogeneous turbulence, for example, "a" takes the form:
a = , where Pw is the probability density function (pdf) representing the turbulence.
The model uses a skewed pdf made up of three Gaussian functions, one representing updrafts (+), one representing downdrafts (-), and the third representing buoyancy-induced turbulence, as follows:
Pw = 1/2 (p G+ + (1-p)
G- + Gb) (3.24)
with G = G(m, s)
= (1/
s) .
exp [ - (w' - m)2
/ 2s2 ].
Here, p is the probability of the particle being in an updraft, m+ is the mean velocity in the updraft, and s+ is the velocity standard deviation in the updraft, and similarly for the downdraft (with "-" instead of "+"). The first three moments of w' can then be equated to the values
= 0, = swï2
+ swb2 , and =Sw3
, to give the following solutions
for p, m+ and m-
: p = 1/2 [1 - ( s2 / (8 + s2
)] 1/2 ,
m+ = 1/2 (swï2
+ swb2) (1-p)/p
, and (3.25)
m- = -m+
p/(1-p),
where, s = 2 Sk / [2 + (swb
/ swï)2
] 3/2 .
The homogeneous turbulence parameterizations used are
swï
= 0.6w* , Sk = (Sw/swï)3
= 0.4, eï
= 0.6w*3 / zi
, and C0 = 2.0.
The model results agreed well with observations in the middle
80% of the CBL.
The turbulence properties associated with buoyancy-generated turbulence
are as follows:
sub2 = svb2
= (b' wb)2
, swb2 = (b'
wb/2)2 , and
eb = 1.5 wb2/z'.
Note that in Gaussian (homogeneous) turbulence, the total turbulence
affecting a particle can be accounted for by simply adding the
variances of velocity fluctuations and by adding the eddy dissipations
due to ambient and buoyancy-generated turbulence, i.e.
s2 = swï2
+ swb2 and e
= eï
+ eb .
The above model was used to simulate various plume rise and entrapment
laboratory experiments of Willis and Deardorff under capping conditions
with an inversion lid. It was shown that the model could reproduce
the observations when a value of b
= 0.7 (instead of 0.6) was used. The justification given for this
higher value was that b represented
the entrainment effect of not only buoyancy-generated turbulence,
but also of ambient turbulence (quite substantial in the convective
case, as here).
In conclusion, it appears that it is perhaps a little premature yet to try to include Lagrangian particle modeling for plume rise in a practical regional model such as the CAPITA Monte Carlo model. At the right time in the future, however, this would be the natural extension of C-MC.
Anthes R. A., E. Y. Hsie and Y. H. Kuo (1987). A description of
the fourth-generation Penn State-NCAR Mesoscale Model (MM4). NCAR
Technical Note, NCAR/TN-282 + STR.
Blackadar A. K. (1979). High-resolution models of the planetary
boundary layer. In Advances in Environmental and Scientific
Engineering, Vol. 1, Gordon and Breach, NY.
Boybeyi Z. (1993). A mesoscale atmospheric dispersion modeling
system for simulations of topographically induced atmospheric
flows and air pollution dispersion. Ph. D. Dissertation, Dept.
of Marine, Earth and Atmospheric Sciences, NC State U., Raleigh,
NC.
Dudhia J. (1993). A nonhydrostatic version of the Penn State-NCAR
Mesoscale Model: Validation tests and simulation of an Atlantic
cyclone and cold front. Monthly Weather Rev. 121, 1493-1513.
Erbrink H. J. (1994). Plume rise in different atmospheres: a practical
scheme and some comparisons with lidar measurements. Atmos.
Environ . 28, 3625-3636.
Ettling D. J., J. Preuss and M. Wamser (1986). Application of
a random walk model to turbulent diffusion in complex terrain.
Atmos. Environ. 20, 741-747.
Gifford F. A. (1982). Horizontal diffusion in the atmosphere:
a Lagrangian-dynamical theory. Atmos. Environ. 16, 505-512.
Glendening J. W., J. A. Businger and R. J. Farber (1984). Improving
plume rise prediction accuracy for stable atmospheres with complex
vertical structure. JAPCA 34, 1128-1133.
Grell G.A., J. Dudhia and D.R. Stauffer (1994). A description
of the fifth-generation Penn State-NCAR Mesoscale Model (MM5).
NCAR Technical Note, NCAR/TN-398 + STR.
Hanna S. R. (1982). Applications in air pollution modeling. In
Atmospheric Turbulence and Air Pollution Modeling, Nieustadt
F.T.M. and H. Van Dop (eds.), D. Reidel Publ., Dordrecht, pp.
275-310.
Hurley P. J. and W. L. Physick (1993). A skewed homogeneous Lagrangian
particle model for convective conditions. Atmos. Environ.
27A, 619-624.
Kaimal J. C., J. C. Wyngaard, D. A. Haugen, O.R. Cote, Y. Izumi,
S. J. Caughey and C. J. Readings (1976). Turbulence structure
in the convective boundary layer. J. Atmos. Sci. 33, 2152-2169.
Legg B. J. and M. R. Raupach (1982). Markov-chain simulation of
particle dispersion in inhomogeneous flows: the mean drift velocity
induced by a gradient in Eulerian velocity variance. Boundary
Layer Met. 24, 3-13.
Ley A. J. and D. J. Thompson (1983). A random walk model dispersion
in a diabatic surface layer. Quart. J. Roy. Met. Soc. 109, 847-880.
McNider R. T. (1981). Investigation of the impact of topographic
circulations in the transport and dispersion of air pollutants.
Ph. D. Dissertation, Dept. of Env. Sciences, U. of Virginia, Charlottesville,
VA.
McNider R. T., M.D. Moran and R.A. Pielke (1988). Influence of
diurnal and inertial boundary layer oscillations on long-range
dispersion. Atmos. Environ. 22, 2445-2462.
McNider R. T. (1989). Mesoscale and long-range pollutant transport
and diffusion. Invited Review Paper presented at the 5th Scientific
Assembly of Int. Assoc. Met and Atm. Phys., U. of Reading, U.K.
Pielke R. A. et al. (1992). A comprehensive meteorological modeling
system - RAMS. Meteorol. Atmos. Phys. 49, 69-91.
Placet M., R. E. Battye, F. C. Fehsenfeld and G. W. Bassett (1990).
Emissions involved in acidic deposition processes. SOS/T Report
1, National Acid Precipitation Assessment Program, US Govt. Printing
Office, Washington, DC.
Pleim J. E. and J. S. Chang (1992). A non-local closure model
for vertical mixing in the convective boundary layer. Atmos.
Environ. 26A, 965-981.
Sawford B. L. (1984). The basis for, and some limitations of,
the Langevin equation in atmospheric relative dispersion modeling.
Atmos. Environ. 11, 2405-2411.
Sawford B. L. (1985). Lagrangian simulations of concentration
mean and fluctuation fields. J. Clim. Appl. Met. 24, 1152-1166.
Stull B. (1988). An Introduction to Boundary Layer Meteorology.
Kluver Academic Publ.
Yamartino R.J. (1993) The use of kinematic simulation techniques
for advanced mesoscale dispersion models. Paper presented at the
AWMA Int. Conf. on the Role of Meteorology in Managing the Environment
in the '90s, Scottsdale, AZ, 26-29 Jan. 1993.
Uliasz M. (1992). The Mesoscale Dispersion Modeling System, v.
2.15. Technical Report, Dept. of Atmospheric Science, Colorado
State University, Ft. Collins, CO.
Uliasz M. and R. A. Pielke (1993). Implementation of Lagrangian
particle dispersion model for mesoscale and regional air quality
studies. In Air Pollution, P. Zanetti, C. A. Brebia, J.E.G.
and G. A. Milian (eds.), Computational Mechanics Publ., Southampton,
pp. 157-164.
Uliasz M. (1994). Source- and receptor-oriented regional dispersion
modeling using Lagrangian particle model. Paper presented at the
2nd RAMS Users Workshop, Ft. Collins, CO, 15-17 February, 1994.
Zanetti P. (1984). A new Monte Carlo scheme for simulating Lagrangian
particle diffusion with wind shear effects. Appl. Math. Modeling
8, 188-192.
Zanetti P. (1986). Monte Carlo simulation of auto- and cross-correlated
turbulent velocity fluctuations (MC-LAGPAR II Model), Environmental
Software 1, 26-30.
Anfossi D., E. Ferrero, G. Brusasca, A. Marzorati and G. Tinarelli
(1993). A simple way of computing buoyant plume rise in Lagrangian
stochastic dispersion models. Atmos. Environ . 27A, 1443-1451.
Briggs G. A. (1969). Plume Rise. USAEC Critical Review
Series, TID-25075, NTIS (81pp).
Briggs G. A. (1975). Plume rise predictions, pp. 59-111 in Lectures
on Air Pollution and Environmental Impact Analyses, D. A.
Haugen (ed.), AMS, Boston.
Briggs G. A. (1984). Plume rise and buoyancy effects, pp. 327-366
in Atmospheric Science and Power Production, D. Randerson
(ed.), US Dept. of Energy, DOE/TIC-27601 (NTIS DE84005177).
Byun D. W. and F. S. Binkowski (1991). Sensitivity of RADM to
point source emissions processing. Paper presented at the 7th
Joint AMS/AWMA Conf. on Application of Air Pollution Met., New
Orleans, Jan 14-18 1991.
Csanady G. T. (1973). Effect of plume rise on ground level pollution.
Atmos. Environ . 7, 1-16.
Deardorff J. W. (1972). Numerical investigation of neutral and
unstable planetary boundary layers. J. Atmos. Sci . 29,
91-115.
Hoult D. P., J. A. Fay and L. J. forney (1969). A theory of plume
rise compared with field observations. JAPCA 19, 585-590.
Hoult D. P. and J. C. Weil (1972). A turbulent plume in laminar
crossflow. Atmos Environ ., 6, 513-531.
Hurley P. and W. Physick (1993). Lagrangian particle modeling
of buoyant point sources : plume rise and entrapment under convective
conditions. Atmos. Environ . 27A, 1579-1584.
Netterville D. D. J. (1990). Plume rise, entrainment and dispersion
in turbulent winds. Atmos. Environ . 24A, 1061-1081.
Petersen R. L. and R. M. Dowd (1988). Discussion on "Proposed
pragmatic methods for estimating plume rise and plume penetration
through atmospheric layers" by Turner (1985). Atmos. Environ
. 22, 2061-2062.
Slawson P. R. and G. T. Csanady (1967). On the mean path of buoyant
bent-over chimney plumes. J. Fluid Mech . 28, 311-322.
Turner D. B. (1985). Proposed pragmatic methods for estimating
plume rise and plume penetration through atmospheric layers. Atmos.
Environ . 19, 1215-1218.
Turner D. B., T. Chico and J. A. Catalano (1986). TUPOS: A multi-source
Gaussian dispersion algorithm using on-site turbulence data. EPA/600/8-86/010,
US EPA, Research Triangle Park, NC (155 pp.).
van Dop H. (1992). Buoyant plume rise in a Lagrangian framework.
Atmos. Environ . 26A, 1335-1346.
Weil J. C. and D. P. Hoult (1973). A correlation of ground-level
concentrations of sulfur dioxide downwind of the Keystone stacks.
Atmos. Environ . 7, 707-721.
Weil J. C. (1982). Source buoyancy effects in boundary layer diffusion,
pp. 235-246 in Workshop on the Parameterization of Mixed Layer
Diffusion , Phys. Sci. Lab., New Mexico State U., Las Cruces.
Weil J. C. (1988). Plume Rise, pp. 119-166 in Lectures on Air
Pollution Modelling, A. Venkatram and J. C. Wyngaard (eds.),
AMS, Boston.
Willis G. E. and J. W. Deardorff (1974). A laboratory model of
the unstable planetary boundary layer. J. Atmos. Sci .
31, 1297-1307.
Zannetti P. and N. Al-Madani (1984). Simulation of transformation,
buoyancy and removal processes by Lagrangian particle methods.
pp 733-744 in Proc. 14th NATO/CCMS Int. Tech. Meeting on Air
Pollution Modeling (Ch. Wispelaere, ed.), Plenum Press, NY.
4.2 Plume Rise
| Submit your comments, feedback, questions, and ideas pertaining this page. Your input will be automatically added to the existing annotations. In order to add a new comment, you must be registered with the CAPITA People's Page. |
