PB80-227580
Development of Mesoscale Air Quality
Simulation Models. Volume 1. Comparative
Sensitivity Studies of PUFF, PLUME, and GRID
Models for Long Distance Dispersion
Environmental Research and Technology, Inc.
Concord, MA
Prepared for
Environmental Sciences Research Lab.
Research Triangle Park, IC
Sep 79
•£. Department of Commerce
il Technical information Service
-------
September 1979
DEVELOPMENT OF MESOSCALE AIR QUALITY
SIMULATION MODELS
VOLUME 1. COMPARATIVE SENSITIVITY STUDIES
OF PUFF, PLUME, AND GRID MODELS
FOR LONG-DISTANCE DISPERSION
Arthur Bass, Carl W. Benkley,
Joseph S. Scire, and Charles S. Morris
Environmental Research S Technology, Inc.
696 Virginia Road
Concord, MA 01742
NOAA Contract No. 03-6-02-35254
Project Officer
Herbert J. Viebrock
Meteorology Laboratory
National Oceanic and Atmospheric Administration
Research Triangle Park, NC 27711
OFFICE OF RESEARCH AND DEVELOPMENT
ENVIRONMENTAL PROTECTION AGENCY
WASHINGTON, DC 20460
-------
3037?-tm
REPORT DOCUMENTATION i. "ipowt EPA/DFV8Q/0Q7d 1 '¦
PAGE EPA-60^-80-056 j
P 880-227580
4, T.t»* ami SubMte
DEVELOPMENT OF MESOSCALE AIR QUALITY SIMULATION MODELS.
Volume 1. Comparative Sensitivity Studies of PUFF, PLUME,
and GRID "lode!s for Long Distance Dispersion
5. Report Date
September 1979
6.
7. Author*,) Arthur Bass, Carl W. Benkley, Joseph S. Scire,
and Charles S. Morris
8. Performmg Organixatton Rept. No,
3. P*rlo»min| Or£an«ration Name and Address
Environmental Research & Technology, Inc.
696 Virginia Road
Concord, Massachusetts 01742
10. Project/Tasfc/Worfc Unit No.
N/A
11. Contract(C) or Grant(G) No.
(C> 03-6-02-35254
(C)
12, Sponsoring Organisation Name and Address
Environmental Sciences Research Laboratory
Office of Research and Development
U.S. Environmental Protection Agency
Research Triangle Park, North Carolina 27711
13. Type of Report & Period Covered
Final Report
14.
15. Supplementary Notts
For magnetic tape see PB 80- -2-2 7 Y 7
16, Abstract (Lsmit; 200 words)
This report, Volume 1 in a series entitled 'Development of Mesoscale
Air Quality Simulation Models,' provides detailed comparisons and sensitivity
analyses of three candidate models, MESOPLUME, MESOPUFF, and MESOGRID. This was
not a validation study; there was no suitable regional air quality data base for
the Four Corners area. Rather, the models have been evaluated by analyzing: V
model response to varying mesoscale flow regimes, model sensitivity to systematic
variations in meteorological and model design parameters; and model accuracy
versus computational requirements and cost.
colcr imwco)
I.N BL.iK a;.d white
17. Document Analyst* a. Descriptors
Air pollution
Atmospheric models
A1gorithms
Dispersion
b. ldent»fters/Op«n-Ended Terms
NATIONAL TECHNICAL
INFORMATION SERVICE
C. COSATI F«eld/Grouo
in dipmtmiit of contact
SWMftflMA 2ZAI
IS. A«a«laM*ty Statement
19. Security Class (This Report)
21. No. of Pages
RELEASE UNLIMITED
UNCLASSIFIED
221
20. Security Class (This Pane)
22. Price
UNCLASSIFIED
*
(See ANSI~m 181
See ln»truct*on« on Reverse
OPTIONAL FORM r*2 (4-;?)
trofmerly NTtS-3*)
Depart mr*t of Commerce
-------
iWUmNMENfUlOtStdRCMHTtCMNtSlOO* f
DISCLAIMER
Publication of this report does not signify that the contents
necessarily reflect the views and policies of the U.S. Environmental
Protection Agency, nor does mention of trade names or commercial
products constitute endorsement or recommendation for use.
iii
-------
TABLE OF CONTENTS
Page
LIST OF ILLUSTRATIONS vii
LIST OF TABLES xiii
1. INTRODUCTION 1
1.1 Overview 1
1.2 Mesoscale Modeling Approaches 4
1.3 Integrated Mesoscale Modeling System 8
1.3.1 MESOPAC 8
1.3.2 MESOPLUME 13
1.3.3 MESOPUFF 15
1.3.4 MESOGRID 15
1.3.5 MESOFILE 18
1.4 Organization of Report 18
2. MODEL COMPARISON AND MODEL SENSITIVITY TESTS 22
2.1 Requirement for Objective Criteria 22
2.2 Test Protocol 22
2.3 Initialization and Inventor)' Choices 25
2.4 Description of Base Cases 25
2.4.1 Model Design and Run Parameters 25
2.4.2 Meteorological Scenarios 27
2.5 Description of the Model Output Analyses 27
3. COMPARISONS OF MODEL BASE CASE RESULTS 45
3.1 Qualitative Comparisons for Single and Multiple
Point Sources 45
3.2 Quantitative Comparisons of Model Results 73
3.3 Model Simulations of Plume Fumigation 78
4. RESULTS OF MODEL SENSITIVITY TESTS 92
4.1 Sensitivity to Decay and Dry Deposition Processes 92
4.2 Sensitivity to Mixing Depth Variations 99
4.3 Sensitivity to Systematic Variation of Wind
Direction and Crosswind Growth Rate 105
4.4 Model Sensitivity to Variations in Wind Speed 112
4.5 Sensitivity of MESOGRID to Variations in 114
v
-------
TABLE OF CONTENTS (Continued)
Page
5. MODEL COST/SENSITIVITY ISSUES 120
5.1 Overview 120
5.2 Model Execution Cost Formulas 120
5.2.1 MESOPAv. Run Time 120
5.2.2 MESOPLUME Run Time 121
5.2.3 MESOPUFF Run Time 121
5.2.4 MESOGRID Run Time 122
5.3 Model Sensitivity to Changes in Time Step 124
6. SUMMARY, CONCLUSIONS, AND RECOMMENDATIONS 131
6.1 Basic Model Development Objectives 131
6.2 Model Verification 132
6.3 Comparative Model Performance 132
6.4 Model Sensitivity Results 133
6.5 Adequacy of Meteorological Data 135
6.6 Comparative Costs 136
6.7 Recommended Model Choices 137
6.8 Other Applications 137
REFERENCES 140
APPENDIX A MESOPLUME TECHNICAL DISCUSSION
APPENDIX B MESOPUFF TECHNICAL DISCUSSION
APPENDIX C MESOGRID TECHNICAL DISCUSSION-
APPENDIX D SCALE ANALYSIS FOR TERMS IN MASS CONTINUITY EQUATION
vi
-------
EWtRONWfNTAi «€S£AW:m* TfCHN€*.0«v NC
LIST OF ILLUSTRATIONS
Figure Page
1-1 Existing and Proposed (through 1987) Major
Energy Sources in the Four Corners Region 2
1-2 lite Variable Trajectory Plume Model Approach 5
1-3 The Variable Trajectory Puff Model Approach 6
1-4 Integrated Modeling System 9
1-5 Mesoscale Meteorology Package (MESOPAC) 10
1-6 Rawinsonde Stations in the General Four Corners
Area 11
1-7 Mesoscale Plume Model (MESOPLUME) 14
1-8 Mesoscale Puff Model (MESOPUFF) 16
1-9 Mesoscale Grid Model (MESOGRID) 17
1-10 File Management and Analysis (MESOFILE) 19
1-11 Examples of MESOFILE Graphical Displays 20
2-1 A Schematic Representation of Episodes 1 and
2 Flow Regimes 28
2-2 Episode 1 Wind Vector and Isotach (m/s)
Fields 29
2-3 Episode 1 Mixing Depth (m) Fields 33
2-4 Episode 1 PGT Stability Class Fields 35
2-5 Episode 2 Wind Vector and Isotach (m/s)
Fields 37
3-1 Episode 1 Day 167 Single Source Dispersion
Model Comparisons (S02) 46
3-2 Episode 1 Day 167 Single Source Dispersion
Model Comparisons (S0|) 49
vii
-------
LIST OF ILLUSTRATIONS (Continued)
Figure Page
3-7 Episode 2 Day 124, 125, 126, and 127 Single
Source Dispersion Model Comparisons (SO^)
50
53
3-3 Episode 1 Day 169 Single Source Dispersion
Mo^el Comparisons (SO,)
3-4 Episode 1 Day 169 Single Source Dispersion
Model Comparisons (S0|)
3-5 Episode 1 Day 167 Multiple Source Dispersion
Model Comparisons (SO,) 55
3-6 Episode 1 Day 169 Multiple Source Dispersion
Model Comparisons (SO,) 58
67
3-8 Days 167, 168, and 169 Comparison of the
Different Models Results (Correlation
Coefficient r) 74
3-9 Days 123, 124, 125, 126, and 127 Comparison
of the Different Models Results (Correlation
Coefficient r) 74
3-10 Episode l--Comparison of the Number of Points
in the Plume NP Computed by Each Model 76
3-11 Episode 2--Comparison of the Number of Points
in the Plume NP Computed by Each Model 76
3-12 Episode 1—Comparison of the Different Model
Results [Fractional Deviation of the Means
(FDM)] 77
3-13 Episode 2—Comparison of the Different Model
Results [Fractional Deviation of the Means
(FDM)] 77
3-14 MESOPLUME Time Sequence of 1-hour Average
S0? Concentrations (7 pm MST Dav 166 to 1 pm
Day 167) " 79
3-15 Time Sequence of Mixing Depths (7 pm MST
Day 166 to 1 pm Day 167) 82
3-16 MESOPLUME Mean Plume Concentration Cycle
(for 1-liour Averages) 85
viii
-------
€\vino^Nmf»ES€APCH*^CMNaocv»c
LIST OF ILLUSTRATIONS (Continued)
Figure Page
3-17 MESOPUFF Time Sequence of 1-Hour Average SO^
Concentrations (7 pm MST Day 166 to 1 pm
Day 167) 86
3-18 MESOGRID Time Sequence of 1-Hour Average S02
Concentrations (7 pm MST Day 166 to 1 pm
Day 167) 89
4-1 MESOPLUME Sensitivity to Decay Rate and Dry
Deposition Variations (Days 167-169) 94
4-2 MESOPUFF Sensitivity to Decay Rate and Dry
Deposition Variations (Days 167-169) 96
4-3a MESOPLUME Fractional Difference Field
Comparison for v, (S0„) = 0.1 m/s Case vs.
Base Case (Day 169) 97
4-3b MESOPLUME Fractional Difference Field
Comparison for v, (SO.) = 0.001 m/s Case
vs. Base Case (Day 169) 97
4-3c MESOPLUME Fractional Difference Field
Comparison for v, (SO.) = 0.1 m/s Case vs.
Base Case (Day 167)
4-3d MESOPLUME Fractional Difference Field
Comparison for v^ (S0_) = 0.001 m/s Case
vs. Base Case (Day 167)
98
98
4-4a MESOPLUME Fractional Difference Field
Comparison for k1 = -S.56 x 10"^ Case vs.
Base Case (Day le>9) 100
4-4b MESOPLUME Fractional Difference Field
Comparison for k = -5.56 x 10"7 Case vs.
Base Case (Day 1&9) 100
4-4c MESOPLUME Fractional Difference Field
Comparison for k.. = -5.56 x 10*^ Case vs.
Base Case (Day 167) 101
4-4d MESOPLUME Fractional Difference Field
Comparison for k. = -5.56 x 10"7 Case vs.
Base Case (Day 167) 101
ix
-------
LIST OF ILLUSTRATIONS (Continued)
Figure Page
4-5 MESOPUFF Sensitivity to Mixing Depth
Variations (Days 167-169) 102
4-6a MESOPLUME Fractional Difference Field
Comparison for Mixing Depth Multiplier =
2.0 Case vs. Base Case (Day 167) 103
4-6b MESOPLUME Fractional Difference Field
Comparison for Mixing Depth Multiplier =
0.5 Case vs. Base Case (Day 167) 103
4-6c MESOPLUME Fractional Difference Field
Comparison for Mixing Depth Multiplier «
2.0 Case vs. Base Case (Day 169) 104
4-6d MESOPLUME Fractional Difference Field
Comparison for Mixing Depth Multiplier =
0.5 Case vs. Base Case (Day 169) 104
4-7 MESOPLUME Sensitivity of r ( ) and FDM (---)
to Systenatic Variations in Wind Direction A0
and Plume Growth Rate do^/dt (Day 167) 108
4-8 MESOPLUME Sensitivity of r ( ) and FDM ( )
to Systematic Variations in Wind Direction AS
and Plume Growth Rate do^/dt (Day 169) 109
4-9 MESOPUFF Sensitivity of r ( ) and FDM (---)
to Systematic Variations in Wind Direction A0
and Plume Growth Rate do^/dt (Day 167) 110
4-10 MESOPUFF Sensitivity of r ( ) and FDM (---)
to Systematic Variations in Wind Direction A0
and Plume Growth Rate do^/dt (Day 169) 111
4-11 Effect of Wind Speed on Idealized Plume Shape 113
4-12a MESOPLUME Sensitivity of r to Wind Speed
Multiplier (a) Variations (Days 167-169) 115
4-12b MESOPUFF Sensitivity of r to Wind Speed
Multiplier (a) Variations (Days 167-169) 115
4-13a MESOPLUME Sensitivity of FEW to Wind Speed
Multiplier (a) Variations (Days 167-169) 116
4-13b MESOPUFF Sensitivity of FDM to Wind Speed
Multiplier (a) Variations (Days 167-169) 116
x
-------
LIST OF ILLUSTRATIONS (Continued)
Figure Page
4-14 MESOGRID Sensitivity to Variations in K
(Day 167) 117
4-15a MESOGRID Difference Field Comparison for
K Multiplier =0.1 Case vs. Base Case
(Say 167) 118
4-15b MESOGRID Difference Field Comparison (Limited
Grid) for Requested Time Step = 5 Minutes
Case vs. Base Case (Day 167) 118
5-la MESOPUFF CPU Run Time (t) and Sensitivity
of r to Variations in Sampling Rate and
Puff Emission Rate (Days 167-169) 125
5-lb MESOPUFF CPU Run Time (t) and Sensitivity
of FDM to Variations in Sampling Rate and
Puff Emission Rate (Days 167-169) 126
5-2 Variation of Correlation Coefficient with
Puff Emission Rate 127
5-3 Comparative Model Run Execution Times
(5-day run) 129
A-1 Schematic Representation of the Segmented
Plume 146
A-2 One Possible Arrangement of Sampling and
Basic Computational Grids 148
A-3 Calculation of the Trajectory of a Plume
Segment Endpoint 151
A-4 Calculation of the Concentration C(i,j) by
the Sampling Function 157
A-5 A Schematic Representation of Adjacent Plume
Segments in Rapidly Decelerating Flow 160
A-6 A Schematic Representation of Adjacent Plume
Segments in Highly Curvilinear Flow 160
A-7 Response of Two Plume Elements to Changes in
Mixing Depth 164
A-8 MESOPLUME vs. the Conventional Gaussian Plume 168
XI
-------
-------
FN'. V*, Mf S£:>: ,• n
LIST OF TABLES
Table Page
2-1 Model Sensitivity Test Protocol 23
2-2 Nominal Design Parameters 26
2-3 Statistical Measures Used 43
4-1 Values of the Deposition Velocities and Decay
Rates Used in the Sensitivity Tests 93
4-2 Matrix of Kind Shift and Growth Rate
Sensitivity Cases 106
5-1 Maximum Values of u and vs. Time Step At 123
6-1 Recommended Model Choices 138
A-1 Coefficients for Dispersion Parameter Formulas 155
A-2 Comparison of Cu/Q Values for MESOPLUME and
Turner Workbook-Computed Values 167
B-1 Coefficients for Dispersion Parameter Formulas 181
B-2 Comparison of Cu/Q Values for MESOPUFF and
Turner Workbook-Computer Values 192
-1 -1
C-l Inverse Monin-Obukhov Length L (m ) as a
Function of PCT Stability Class, for Roughness
Length = 25 cm (After Golder 1972) 208
C-2 Stability-Dependent Expressions for in the
Surface Laver and 3K
210
C-3 Maximum Values of u and K , K., vs. Time
Step At, for Nominal Model'Parameters 212
C-4 Determining when MESOGRID will Modify K_
Profiles 2 214
xiii
-------
f»*V'»*ONM£NT^ *«•
I. INTRODUCTION
1.1 Overview
In response to a growing national commitment to the greatly
expanded use of indigenous coal reserves to meet energy generation
demands, the Four Corners region of the southwestern United States will
experience increased development of coal-fired steam electric power
plants and coal gasification and oil shale recovery facilities. But, as
mandated by federal statutes and regulations reflecting the strong
nationwide support for environmental protection, further coal-based
energy resource development (ERD) will only be permitted where con-
sistent with the maintenance of human health, welfare, and environmental
quality. At the present time, the permitting of a new major coal-
burning facility will be most affected by the regulations governing
emissions and ambient air concentrations of criterion air pollutants, as
well as by their effects on other air quality-related values, e.g.,
visibility and acid rain.
For major new coal-fired steam electric power plants and other
coal-based major ERD point sources, the pollutant of most concern at the
present time is sulfur dioxide (SO?), especially considering the growing
body of evidence linking SO? emissions to increased concentrations of
sulfate aerosols (SO.J . Under the present statutes, the permitting of a
new (or expanded) major source of S02 emissions is subject to both the
existing National Ambient Air Quality Standards (NAAQS), as well as the
new source review requirements for the Prevention of Significant Deterio-
ration (PSD) of existing ambient air quality.
Thus, the siting and generation capacity of coal-based energy
facilities in the Four Corners region will be increasingly constrained
by the expected impacts of these facilities on both local and regional-
scale ambient air quality. To illustrate the regional nature of the
problem, Figure 1-1 displays one proposed scenario for the locations of
major energy-related SO2 point sources in the Four Corners region
projected to 1987.
Suitable simulation tools are clearly needed to assess the impacts
of different energy development scenarios on air quality throughout the
Four Corners region. To meet this need in the public interest, the
National Oceanic and Atmospheric Administration (NOAA) has sponsored a
study by Environmental Research & Technology, Inc. (ERT) to develop,
compare, evaluate, and exercise a number of alternative approaches to
regional-scale ambient air quality modeling in the Four Corners area. A
major objective of this study has been to provide a suite of air quality
simulation models that are both technically sound and computationally
practical for assessing regional-scale impacts of energy development
scenarios.
1
-------
7901110
f
Figure 1-1 Existing and Proposed (through 1987) Major Energy Sources
in the Four Corners Region
-------
fNv-«ONMIS-*A. »£<*«*->-»TtC«NT« DC- »
-------
^sea^sTfCHNOtoc*
-------
Figure 1-2 The Variable Trajectory Plume Model Approach
-------
f WBONWeV'V RgS£«.«C^» nEC«NOi.OG* -NC
r i\*\
i N« I
o
6^~
vJ
*
Figure 1-3 The Variable Trajectory Puff Model Approach
I
-------
that take into account the often-important variability of plume
trajectories, especially over regional transport scales.
Variable-trajectory models have some major advantages for regional-
scale impact assessments, as compared to the conventional straight-line
Gaussian models:
• Unlike the standard Gaussian plume model, the variable-
trajectory model does not assume that the wind field persists
uniformly over periods of many hours to days.
• On regional scales, the spatial and temporal variations in
transport wind flows can produce large meanders in plume
trajectories. For the models described here, meander effects
can dominate plume dispersion on regional scales.
• Straight-line models are invalid for describing plume dis-
persion under strongly recirculating or stagnating meteoro-
logical flow conditions.
On the other hand, variable trajectory models do have disadvantages:
• They require substantially more computation time and computer
resources than do straight-line models;
• The technical problem of distinguishing (a) plume dispersion
effects - resulting from meander of the mean wind from (bj smaller-
scale turbulent eddy-induced dispersion represented by the
conventional Pasquill-Gifford-Turner dispersion coefficients -
remains an open issue, and the choice of dispersion parameters
for such models is still very uncertain (Nappo 1978). On
balance, however, many experts advocate the use of variable-
trajectory approaches when warranted by spatially and tem-
porally varying meteorological conditions. The uncertainty
in plume dispersion is, in many cases, of secondary
importance - as this report will show - when compared to the
uncertainties in characterizing the wind field.
• To date, little validation work has been done with variable
trajectory models, principally for lack of adequate obser-
vational data bases; much more research, development and
verification work should be done; some is now under way.
• These modeling approaches are new and unfamiliar. If they are
to be used for regulatory applications, both the regulatory
agencies and the air quality modeling community must be made
familiar with their advantages and proper use for air quality
impact assessments over mesoscale transport ranges.
7
-------
(Vv*aONu«s<*. at<*«*-~•* 'V,
1.3 Integrated Mesoscale Modeling System
To make as clear as possible a comparison among the respective
advectivc-diffusive approaches, the original models were completely
restructured, so that MESOPLUME, MESOPUFF and MESOGRID are identical in
every aspect other than their advect ion-diffusion algorithms. Thus, the
three new models share, where possible, the identical modules for plume
rise, plume decay and deposition, plume growth with distance or travel
time, and treatment of plume fumigation by entrainment within the growing
mixed layer.
The dispersion models have been incorporated in an efficient,
easy-to-use "integrated mesoscale modeling system". This system is
schematically depicted in Figure 1-4. It consists of a closelv-
integrated set of independent models that perform discrete functions:
meteorological preprocessing; plume transport and diffusion; postprocess-
ing and analysis. The mesoscale meteorology preprocessor called MESOPAC
(Benkley and Bass 1979c) drives any of the three mesoscale transport-
diffusion models (MESOPLUME, MESOPUFF, or MESOGRID). In turn each of
these models communicates its results in identical manner to the post-
processing system, MESOFILE, responsible for file manipulation, display,
and statistical analysis of all model output fields. Each component of
this modeling system is described briefly below; complete descriptions
of each of the five components are provided in the respective model
User's Guides (companion reports, Volumes 2-6).
1.3.1 MESOPAC
Figure 1-5 illustrates the principal features of the mesoscale
meteorological package (MESOPAC). MESOPAC is used to construct hourly
time sequences of spatially variable, gridded meteorological fields.
The fields produced are:
• horizontal wind components (u,v) at a user-specified level,
• mixing depth, and
• Pasquill-Gifford-Turner (PGT) stability class.
The input data required by MESOPAC are routinely available from
National Weather Service (NWS) rawinsonde stations.* Figure 1-6 illus-
trates all the NWS rawinsonde stations in the Four Corners region.
Twice-daily observational data from all of these stations were used when
available. [Although the computational domain used for the dispersion
model calculations is bounded by the rectangle shown in the figure, all
of the rawinsonde stations shown were used in constructing the meteoro-
logical fields.]
*Use of MESOPAC is not restricted to the Four Corners geographic region.
8
-------
KmmmmttM, plmmx>i & tt-omAom mc
I
I
Figure 1-4 Integrated Modeling System
9
-------
rtpflumtfcffe ki ,s waft *«
Figure 1-5 Mesoscale Meteorology Package (MESOPAC)
10
-------
79011J1
• RAP
BOl
WMC •
• LBF
• DDC
• AMA
#MYF
•TUS • MAF
• ELP
Figure 1-6 Rawinsonde Stations in the Four Corners Region
-------
The MESOPAC model generates a one-level horizontal wind field,
based entirely on the observational station wind data supplied. Neither
terrain channeling influences nor variations in wind with height above
terrain are assumed in this model. As precedents for using single-level
wind fields to make regional-scale transport-diffusion calculations,
other published simulation studies of plume transport on the regional
scale have used advecting wind fields invariant with height [e.g., Liu
and Durran (1977) or Hales et al. (1977)]. Indeed, for the approaches
typified by MESOPLUME and MESOPUFF, the use of layer-averaged winds is a
practical necessity to avoid splitting plume/puff elements at each time
step. Little incremental benefit is gained by using multiple-layer wind
fields in simulations over rough terrain such as the Four Corners area -
regions for which reliable wind field data are sparse, and for which
any current interpolation scheme is conjectural.
In view of the chief objective - to provide practical modeling
tools for simulations over the Four Corners region - it was decided to
construct horizontal wind field" that might best describe the regional-
scale advection of plumes as averaged over the mixed layer. It is not
trivial to determine the appropriate level or the height interpolation
method to be used: local wind circulations, induced by the complex
topography of the Four Corners region, or by local thermal effects, may
generate near-surface flows that are effectively uncoupled from meso-
scale boundary layer transports. For application in the general Four
Corners area, the rawinsonde data for the 700-millibar (mb) (mandatory)
level was used to construct the horizontal wind field, because the
700-mb level frequently lies within the upper portion of the boundary
layer in the western United States; and, over most of the Four Corners
region, the surface pressure is less than the next mandatory level,
850 rab.*
In brief, MESOPAC interpolates gridded wind field values (u,v) to
model grid points at each time step, using an inverse-square distance
weighting procedure with wind data from one or more user-specified
radiosonde stations. For hours other than 0000 Greenwich Mean Time
(GMT) and 1200 GMT, the initial-guess wind field is synthesized by
linearly interpolation with time from the bracketing 0000 GMT and
1200 GMT wind fields. It is not suggested here that such a procedure is
technically preferred; only that, in the absence of a good prognostic
technique (for example, a limited area fine mesh numerical weather
prediction model), linear interpolation is probably a better assumption
than is simple persistence. A prospective user of one or another of
these models might consider the use of alternative schemes or data
sources to construct the wind fields.
To help ensure pollutant mass conservation by the advection-
diffusion scheme, the divergence of a wind field computed by MESOPAC can
be minimized to within a user-specified absolute tolerance factor
*A wind data adjustment procedure based on surface wind velocities
(e.g., Pack et al. 1978) was not adopted for the present study,
because the surface wind velocity in regions of complex terrain is
likely to be significantly different from that most representative
of transport averaged over the boundary layer.
12
-------
1
t rta+txaG* ®c
[e.g., <_ 110-,+1 (s_1)]. MESOPAC uses an iterative wind field relaxation
scheme adapted from those of Clark and Eskridge (1977) and Liu and
Good in (1976).
The gridded field of mixing depths is generated at each time step
using a technique that weights the relative importance of mechanical and
convective production of turbulence, while accounting for diurnal vari-
ations. The mixing depth is assigned [following Benkley and Schulman
(1979)] as the larger of two terms: (1) a mechanical value (Plate 1971);
(2) a convective value computed for daytime hours as the height of the
intersection of the morning temperature sounding with the adiabat
extending upwards from the observed surface temperature.
The gridded field of PGT stability class is generated at each time
step, using an algorithm closely paralleling that of Turner's (1970);
the essential difference is that MESOPAC uses as a measure of insolation
intensity the convective mixing depth value obtained above. The fields
of PGT stability class are used by each of the three dispersion models
in computing plume rise; and are also used by MESOPUFF and MESOPLUME to
compute plume dispersion rates at downwind distances to 100 km.
Stability values are also used by MESOGRID, in conjunction with mixing
depth and wind speed information, to produce the spatially and tem-
porally variable Kz profiles used in the vertical diffusion algorithm -
see Appendix C.
MESOPAC is described fully in the User's Guide to MESOPAC, Volume 6
of this series (Benkley and Bass 1979c).
1.3.2 MESOPLUME
The MESOPLUME model (Figure 1-7) is a multiple point source,
variable-trajectory Gaussian plume segment model driven by a Lagrangian
wind field. The modeling method used in MESOPLUME has been adapted from
that proposed by Hales et al. (1977). Unlike the conventional Gaussian
plume approach in which the wind is invariant along the plume axis at
any given time, MESOPLUME allows for deformation of the continuous plume
by a temporally varying wind field.
MESOPLUME divides the Gaussian plume into contiguous segments, each
segment describing a portion of plume behavior between successive time
periods (Figure 1-2). The end points of each segment are advected in a
Lagrangian sense. The lateral width of each segment end point depends
on its distance from the source along the instantaneous plume axis, not
on the distance traveled (or travel time) of the segment end point.
For the Four Corners regional study, MESOPLUME has been used with
26 by 26 horizontal grid resolution and a grid spacing of 40 km, for a
13
-------
I»*,«»* i ffi MNCXiX.* ¦»,
Figure 1-7 Mesoscale Plume Model (MESOPLUME)
14
-------
£Wv«« ft ifcCMNCXOl,* 9ft-
total computational domain of 106 fan2.* (Both the grid point resolution
and grid size are user-variable.) It accommodates multiple point
sources (up to 10 per run) and includes modules for Briggs* plume rise,
linear chemistry, deposition, growth, and fumigation.
MKSOPLUME is described fully in the User's Guide, Volume 2 of this
series (Benkley and Bass 1979a).
1.3.3 MESOPUFF
The MESOPUFF model (Figure 1-8) is a variable-trajectory, multiple
point source Gaussian puff model driven by a Lagrangian wind field. The
basic modeling algorithm has been adapted from that described by Start
and Wendell (1974). MESOPUFF simulates a continuous plume by super-
posing discrete puffs (Figure 1-3). Each puff is advected in a
Lagrangian sense; its time history is independent of preceding or
succeeding puffs. The dimensions of an individual puff are proportional
to its travel distance (or travel time).
The grid size and resolution of the MESOPUFF model can be made
identical to MESOPLUME. MESOPUFF also accommodates multiple point
sources (up to 10 per run) and shares the identical supporting modules
for plume rise, linear chemistry, deposition, growth, and fumigation.
With appropriate choice of input parameters, MESOPUFF reproduces to
a very high degree the results of the conventional straight-line
Gaussian plume model, if final (equilibrium) plume rise is assumed. For
the particular choice of inputs in Section B.9, for example, MESOPUFF
results and straight-line plume model results agree within a few
percentage points from 5 to 100 km from a point source.
MESOPUFF is described fully in the User's Guide, Volume 3 of this
series (Benkley and Bass 1979b).
1.3.4 MESOGRID
MESOGRID is a hybrid Lagrangian-Eulerian advection-diffusion model
based on the Egan-Mahoney method of moments (Egan and Mahoney 1972), and
adapted from the SULFA3D model (Egan et al. 1976) - see Figure 1-9.
The model accounts for horizontal advective transport and vertical eddy
diffusion, but not for horizontal diffusion. At each time step, pollu-
tant mass is advected and diffused in a Lagrangian sense; immediately
*If MESOPLUME is run with a much finer grid size (e.g., 5 km) and is
driven by a uniform, stationary wind field until quasi-steady-state
conditions are reached (that is, plume elements are created and fall off
the grid at the same rate), it will give near-field (<_ 100 km) results
within a few percentage points of the Turner Workbook results in the
limiting case of straight-line flow at distances for which the vertical
distribution of concentration can be assumed to be uniform (see
Section A.3).
IS
-------
£WV»«OW*NTAl Mm. >. * . V."* **
?
r»
Figure 1-8 Mesoscale Puff Model (MESOFUFF)
16
-------
Figure 1-9 Mesoscale Grid Model (MESOGRID)
17
-------
tNV*c *•».**! A. 4 •' V - !V
afterwards, a mass decomposition to a stationary Eulerian grid is
performed. The numerical method conserves the zeroth, first, and second
moments of the pollutant mass distribution; this method minimizes
pseudodiffusive errors associated with conventional finite-difference
approximations.
In common with any Eulerian finite-grid model, the MESOGRID model
will not conserve pollutant mass if it is driven by a wind field that
(a) has significant local divergence or convergence or (b) has sig-
nificant confluence, difluence or stagnation.
MESOGRID is described fully in the User's Guide, Volume 4 of this
series (Morris et al. 1979).
1.3.5 MESOFILE
MESOFILE is a flexible postprocessing package consisting of: file
management; file merging and manipulation; statistical analysis; and
graphical display components (see Figure 1-10). The MESOFILE package
records, catalogues, and archives all the relevant output from any of
the mesoscale diffusion models. Its file merging features allow the
user to aggregate easily the results files from any arbitrary number of
separate model runs; thus, different combinations of emissions (e.g.,
the combined effects of existing sources and proposed sources) can be
synthesized without rerunning the transport-diffusion models.
MESOFILE allows the user to specify the time averaging of the model
concentration fields over any multiple of the basic time step. Kind,
mixing height, PGT stability, and ambient air concentration fields can
be contour plotted. Figure 1-11 illustrates the kinds of graphical
displays that can be produced by MESOFILE with access to a suitable
plotter. These displays, used extensively throughout this report, help
greatly in describing and understanding the spatial and temporal vari-
ations in meteorology and ambient concentrations on regional scales.
MESuilLc is described fully in the User's Guide, Volume 5 of this
series (Scire et al. 1979).
1.4 Organization of Report
Section 2 of this report describes the model comparison and
sensitivity test protocol, including the nominal choices of model param-
eters for each of the models, as well as the base case meteorological
scenario data. Section 3 describes the results of the model comparison
runs made for both a single point source and for a set of ten point
sources, for each of two multiday meteorological episodes. Section 4
describes the sensitivity of the models to systematic changes in input
parameters and meteorological fields. These results are referenced to
actual comparative model run costs in Section 5 for different choices of
the relevant model execution parameters. Section 6 presents some
18
-------
BMM^^NCICMKMATKMNaOav *C
Figure 1-10 File Management and Analysis (MESOFILE)
19
-------
TtCHwotoG*
Mixing Height Field
Concentration Isopleths (fig/m*)
I
Figure 1-11 Examples of MES0F1LE Graphical Displays
20
-------
tmrnmrnwrn, «$« awcm « TiCHNOtos*
-------
fNVOONMENTJk flCSfWCH 4 •«*
2. MODEL COMPARISON AND MODEL SENSITIVITY TESTS
2.1 Requirement for Objective Criteria
At the present time, no ambient air quality monitoring network
operating in the Four Corners area seems suitable for verifying short-
term regional-scale transport-diffusion models. [Such a network would
be a well-distributed set of air quality monitors, situated to respond
to regional-scale impacts of the point sources in the general region and
not be unduly influenced by near-field major source impacts or back-
ground sources. SheLh et al. (1978) present a useful discussion of
criteria for the spatial distribution of sampling networks and sampling
intervals for verifying regional-scale transport models.] The lack of
suitable data in the Four Corners area for regional-scale model verifi-
cation is hardly unique - a comparable model development effort for the
northern Great Plains area (Liu and Durran 1977) faced the same diffi-
culty; as a result, many parameters in their model were chosen entirely
on the basis of sensitivity studies.
In the absence of a suitable ambient air quality data base for
model verification in the Four Corners region, the following criteria
were used to provide an objective, systematic basis for comparing model
performance and cost effectiveness:
• theoretical bases of the modeling assumptions,
• model performance under different realistic meteorological
flow regimes,
• qualitative (graphical) and quantitative (statistical)
measures of model sensitivity to different parameters,
• model sensitivity vs. execution cost tradeoffs associated with
changes in the model input variables,
• computer resource requirements, and
• ease of model use and interpretation.
2.2 Test Protocol
Table 2-1 outlines the test protocol for the more than 100 multiday
model sensitivity tests that were conducted. Each of the meteorological
input field variables model design parameters and model run parameters
were varied over the dynamic ranges indicated. It was assumed that all
variables and/or parameters in Sets 1,2, and 3 are uncorrelated. Of
course, certain variables are implicitly interdependent (e.g., the rate
of SO2 deposition influences the amount of S02 remaining for conversion
to SOu); and other variables are explicitly interdependent (e.g., the
vertical eddy diffusivity used by MESOGRID depends on the nnxing depth);
22
-------
CNVtBONMCMAt «r5E*«CH* ?£€M*«Os f'Xi*
TABLE 2-1
MODEL SENSITIVITY TEST PROTOCOL
Range of
Category Variation
Mesoscale Flow Regime
A) Zonal flow
B) Recirculating flow
C) Stagnating flow
Systematic Errors in Meteorological Input
A) Wind speed ± 2x
B) Wind direction ± 22.5°
C) Mixing depth ± 2x
Adjustable Modeling Design Parameters
A) Deposition velocity + lOx
B) Decay rate ± lOx
C) Horizontal plume growth rate ±0.5 m/s
(MESOPLUME, MESOPUFF)
D) Vertical Eddy Diffusivity (MESOGRID) ± lOx
Model Run Parameters
A) MESOPLUME
i) Time step 1-12 hrs
B) MESOPUFF
i) Time step 1-12 hrs
ii) Puff emission rate 1/3 - 12 hrs
iii) Sampling rate 1/3 - 12 hrs
C) MESOGRID
i) "Requested" time step S min - 1 hour
23
-------
. tKSI aht Mik rj c 'V.'
but as a matter of practical necessity, in order to limit the number of
sensitivity tests to be performed in this study, the sensitivity of the
models to simultaneous variations of two or more parameters was not
examined.
In the tests of model response to different mesoscale flow regimes,
the three models, representing three different plume transport-diffusion
modeling approaches, should be expected to perform with varying degrees
of realism under different mesoscale flow regimes. The first set of
sensitivity tests shown in Table 2-1 were designed to compare model
performance for (A) flow regimes in which the models might be expected
to compare well (e.g., highly zonal flow); and by contrast, (B,C) flow
regimes in which the models should differ substantially.
The next set of sensitivity tests shown in Table 2-1, explores the
influences of systematic errors in the meteorological input data. No
single choice of wind level is adequate to represent plume dispersion,
especially over the broad range of effective plume heights, mixing
height variations, ambient turbulence levels, and terrain conditions
encountered in a regional-scale simulation. Different choices of data
level will produce systematically different meteorological data bases.
Furthermore, a variety of different algorithms could be used to con-
struct mixing depth and wind fields. Subgrid scale resolution problems
in the data collection network will cause additional deviations in the
meteorological input data used with the models. It is questionable, for
example, whether interpolated grid values of mixing depth and wind speed
can be estimated to much better than a factor of two, or wind direction
estimated to much better than ±22.5°.
The third set of sensitivity tests shown in Table 2-1 addresses the
effects of present uncertainties in model design parameters: speci-
fically, decay rate, deposition velocity, horizontal plume growth rate
(for MESOPLUME and MESOPUFF), and vertical eddy diffusivity (for
MESOGRID). The ranges of variation assumed for the values of the
deposition velocity and decay rate reflect the large uncertainties in
these quantities. For example, estimates of SO,, decay rates vary over
approximately an order of magnitude, from M3.5% to ^5.0% per hour or
more (Wil son 1978). Quoted deposition velocities also vary over an
order of magnitude, from M3.2 to ^2.0 centimeters per second (cm/s) for
SO, (Husar et al. 1978). For the MESOPLUME and MESOPUFF models, the
long-range, time-dependent dispersion coefficients are highly uncertain;
beyond 100 km transport distance, the range of horizontal plume growth
rates (dov/dt) may be 0.5 ±0.5 m/s [the nominal choise is from Heffter
(1965)]. 'For the MESOGRID model, the algorithm used to create vertical
profiles of Kz probably yields estimates only to within an order of
magnitude at any specific location or height.
Finally, the fourth set of sensitivity tests examines the relative
sensitivity benefits vs. costs for different choices of model run
parameters: advective time step; plume segment or puff emission rates;
24
-------
f •mv-oONMSN'"*-. »ESE *»<> i tf cmncxOG * »•*'
sampling time steps. For MESOPLUME and MESOPUFF, the basic time step
used for the transport-dispersion calculation is varied from one hour to
twelve hours. For MES0GR1D, the "requested" time step (see Section C.4.4)
is varied from five minutes to one hour (by reducing Kz values to ensure
that the vertical numerical stability criterion is not violated). For
MESOPUFF, both the puff emission rate - which also controls the fre-
quency of puff trajectory update - and the puff sampling rate were
simultaneously varied. In the cost-benefit analyses made to explore
model performance with variat i ons in time step, a very stringent (short)
time step was used for the base case; then progressively longer time
step(s) were chosen; in each instance, the results were compared, in
terms of cost and reproducibility of important result features, with the
base case results of the same model.
2.3 Initialization and Inventory Choices
In keeping with the regional-;cale nature of the problems, each
model run experiment consisted of a five (or seven) consecutive day
simulation using hourly meteorology. The first two (or three) days of
each model run were considered to be an initialization period: the
analyses of model results were limited to the remaining three (or four)
days of the simulation. In these experiments, the two-day initial-
ization period is sufficient to allow the plume from any source to reach
the boundary of the computational grid; thus the unrealistic effects of
source start-up are largely avoided.
The model comparison experiments were performed using alterna-
tively: (1) a single point source only; and (2) a set of ten sources
chosen from the actual source distribution shown in Figure 1-1. Two
different meteorological episodes were run. To keep both computational
cost* and analysis efforts within reasonable bounds, the various model
sensitivity tests used only a single major point source of SO2, the
Navajo power plant in northern Arizona. This source was chosen for two
reasons: (a) For the two meteorological episodes considered, the winds
over most of the grid area were generally westerly the plume from a
source in the western portion of the grid would tend to spread eastward
through a significant portion of the grid; and (b) long-range transport
is the primary interest, and a source chosen for the sensitivity
analysis should have significant impact at large distances - the Navajo
plant has the highest emissions of any point source in the Four Corners
region.
2.4 Description of Base Cases
2.4.1 Model Design and Run Parameters
The nominal base-case values chosen for the adjustable model
design parameters are displayed in Table 2-2. These nominal choices are
the same as those used by Hales et al. (1977). Although no assertion is
*The costs of operating the MESOPLUME and MESOPUFF models depend on the
number of sources and a very large number of multiday sensitivity tests
were required.
25
-------
emvsfomiFNTAi. «£§£*ac« s rgowtcxiv »«c
TABLE 2-2
NOMINAL DESIGN PARAMETERS
S07 Decay Rate = 2%/Hour
SO^ Deposition Velocity = 1.0 cm s *
SO" Deposition Velocity = 0.1 cm s 1
4
Plume Growth Rate (MESOPLUME, MESOPUFF):
£ 100 km, PGT coefficients
do .
> 100 km, = 0.5 m s"
BASE CASE VALUES OF THE MODEL RUN PARAMETERS
MESOPLUME:
Time step = 1 hour
Meteorological fields update rate = 1 hour
Grid interval = 40 km
6 2
Computational mesh size = 26 x 26 (10 km )
MESOPUFF:
Time step = 1 hour
Meteorological fields update rate = 1 hour
Puff emission rate = 1 puff/hour
Puff sampling rate dynamically determined (see Section B.6)
Same computational mesh and grid interval
MESOGRID:
"Requested" time step = 1 hour
Meteorological fields update rate = 1 hour
Same computational mesh and grid interval
26
-------
KtJMnCHl TfCMMLOOV WC
being made that these are optimal values, they are near the means of the
range of values documented by the authors previously cited. The choices
of model run parameter base-case values are also shown in Table 2-2.
2.4.2 Meteorological Scenarios
Two very different multiday meteorological episodes (shown sche-
matically in Figure 2-1) were chosen for the model sensitivity tests.
Episode 1, June 14 to 18, 1978 (days 165 to 169), was chosen for
its wel1-defined regimes of both zonal and stagnating flows. Figure 2-2
displays the wind vectors and isotach contours every 12 hours starting
at 0000 GOT on day 167 (after a 2-day model initialization period). The
flow on days 165 to 167 was quite steady from the west and southwest ;
whereas the flow on days 168 and 169 was light and disorganized, especi-
ally in the southeast portion of the grid. By the end of the fifth day,
flow tended to return to zonal conditions near the source; thus, the
concentration patterns on the last day should be affected by both stag-
nating and zonal conditions. Figure 2-3 and Figure 2-4 display the
fields of mixing depth and stability class, respectively, computed for
this episode.
Episode 2, May 1 to 7, 1978 (days 121 to 127), was chosen tor its
distinct zonal and recirculating regimes. Figure 2-5 displays the wind
vectors and isotach contours every 12 hours starting 0000 GMT on day 124
(after a 3-day model initialization period). Two strong cutoff low
pressure systems moved through the grid, the first from late on day 121
to day 123, and the second from day 125 to early on day 127; at other
times the flow was basically zonal. Wind speeds were generally high
throughout the period, except in the centers of the cutoff lows.
2.5 Description of the Model Output Analyses
Because of the large volume of model output data generated, the
sensitivity analyses were designed to yield maximum information from a
limited number of measures or descriptors. Qualitative (graphical) as
well as quantitative (statistical) measures of the point-by-point and
bulk differences between two concentration fields - a "base" (or
"reference") field, and a "test" field - were used to determine the
effects of "test" parameter variations on the predicted concentration
field. The lowest order statistical moments - the mean and standard
deviation - are normally meaningful bulk indicators of the distribution.
Higher order moments car. also be quantified by bulk formulas but are
less meaningful - the higher order details in a concentration field are
best seen in qualitative examination of the concentration isopleth
patterns.
To compare two sets of concentration fields, it is useful, first,
to compare the lower moments of the two distributions, and of their
difference field. But since the model comparisons are also directed at
27
-------
i
» J
00
Figure 2-1 A Schematic Representation of Episodes 1 and 2 Flow Regimes
-------
CMMMONUCMTAt WSMMCH* rCCHMUOo* «C
00 167 78
, . , Tn . ,
¦ • - 4 - ¦
- ' « V » '
.......
::; i::
. , , - .
10
i
s • • L? .1.
" * - -
: r / " "
\ H " * ' ' '
_l..\ -T
i:^-r
•N; v ' -
' 'J I 1
. 10.
- - . „8
12 167 78
Figure 2-2 Episode 1 Wind Vector and Isotach (m/s) Fields
29
-------
mm & n ohmti oa» m:
00 168 78
12 168 78
Figure 2-2 Continued
30
-------
wvi&cmmMM mmuncH * >k>mi jg* *c
00 169 78
12 169 78
Figure 2-2 Continued
31
-------
ewwmwuww mm.m^^nctmoLmsr m.
32
-------
00 167 78 12 167 78
Figure 2-3 Episode 1 Mixing Depth (m) Fields
33
-------
00 169 78
Figure 2-3 Continued
34
-------
trA-mor*e«TAi aes&vCHk recwMuarv* *c
00 167 78 12 167 78
00 168 78 12 168 78
1
/
1 *
! t
1D
p-
j
—V 1
/ 5
r 1
i
•!
1
i (
;! '
J \
C '
B
s
A
\
| 1
| L"'"T" •
j j
L - -
J!
- r
\
h—f—
c>
h i
¦••i !
•
\ , •
»
7902217
\ |
° i
Figure 2-4 Episode 1 PGT Stability Class Fields
-------
oo mm
12 Iffl 78
00 170 78
I
Figure 2-4 Continued
36
-------
EMvMOMUFMTAt nCSt-AIICM * tECMMk-"*.* <*t
00 124 TO
12 124 78
Figure 2-5 Episode 2 Wind Vector and Isotach (m/s) Fields
37
-------
WlMWKTNMf~*lAt WFSfc/WCMft TF» HNC*or.V »«
00 12S78
12 125 TO
Figure 2-5 Continued
38
-------
nMKMKHTK MMMOOTtOWOLOBV
12 126 78
00 126 78
39
-------
mmmo*amMwtmo&v mc
mm 78
12 127 78
Figure 2-5 Continued
40
-------
t V. tftCMMENTAL RESEARCH ft TtCHNC* OtiV (NC
00 128 78
Figure 2-5 Continued
41
-------
««*.>» '{.WW it,'
the more detailed structures of the respective concentration fields,
point-by-point comparison measures of the two fields are also required.
The relative spatial locations of the two fields should also be com-
pared . Consider, for example, a single dispersing plume: if a base-
case concentration field is defined and a test-case field is created by
varying a model input parameter, the two plumes may or may not overlap.
This information is especially important to possible regulatory appli-
cations, for difficult interpretation problems would arise in situations
in which the test-case concentration is high and the base case concen-
tration is low (or conversely). In summary, the statistical and quali-
tative (graphical) measures employed for model comparison analyses
should yield the following information.
• How closely do the lower order moments of the base and test
plumes compare--are the plumes of the same size and do they
have the same mean concentration?
• How closely do the higher order moments of the base and test
plumes compare--is a particular feature of interest similarly
reproduced in both the base plume and test plume?
• How closely do the spatial location of base plume and test
plume or distinct plume features coincide--do the two plumes
overlap closely or do they tend to disperse in different
directions?
Table 2-3 summarizes the statistics found useful in this study. Both
the mean plume concentration C, and the number of points NP over which
the mean is computed, are useful indicators of gross plume size and
magnitude of impact. The plume sine NP is defined for MESOPLUME and
MESOPUFF as the number of points within ±3 Oy of the plume centerline;
and for MESOGRID as the number of points experiencing at least a thres-
hold concentration of 10-*4 ug/m~3. The standard deviation has not been
included as a bulk description of individual plumes; the standard
deviation was found to be more descriptive of the along-the-trajectory
plume distribution than of the lateral distribution of the concentration.
Several statistics have been fr nd useful for the bulk comparison
of two plumes: the average fractional deviation (AFD), the fractional
difference of the plume means (FDM), the fractional difference of number
of points in the plume (FDN), and the correlation coefficient (r).
The average fractional difference is most useful when every point in the
test plume is also included in the base plume. However, this statistic
can take on unrealistically large values when the base case concentra-
tion values at a single point become very small. Thus, this statistic
is really only useful for examining the effect of an input parameter,
such as decay or dry deposition, that does not affect the spatial
coverage of the plumes. For other input parameter variations, the FDM
or the FDN are methods by which the lowest order plume moments are
42
-------
CWVRONMEH'Av Sf *«CH 4 'I-Ch***. OG* *C
TABLE 2-3
STATISTICAL MEASURES USED
Descriptor
Definition
Individual pluae bu'k descriptors
U Mean pluae concentration, C
2) Nuaber of points in the pluae, NT
r.j ci.j /!fp
For MFSGPLUW- and MESOPUFF, NP is
the number of points within
•3 3V of the pluae center!ine.
For fcESQGRIO, VP ;• the nuaber of
points experiencing at least a
threshold concentration of
10"-* , g b~ '
i,} range over only those
points in the flu®e
II. TMO-plune bulk descriptors
1} Average fractional deviation
CAFD)
2) Fractional deviation of the plume
scans (FEW)
pB-fr;
S j
cB-cT
Calculated only over the
intersection of the base (B)
and test (T) pluaes
CB CT " Eg 't
3) Correlation coefficient (rl r * Calculated over the set of
fl2 ,2tp ,jr ,2j grid points in the union of
| B " B j[ T " T j the two pluses
4} Fraction deviation of number of
points in the plu&e (FDN)
III. Two-pluse point-by-point cosparisons
1) Difference field
IDF)
2} Fractional difference field
(FDF)
(OF)
M
(C-)
IO - (CT»M
(C 1 - (C )
(F0F). . - B ''J
' i.j
i,j range over all grid
points
i.j range over only those
grid points in the inter-
section of the two pluses
43
-------
ffcwnowewnk -%r
compared. The correlation coefficient defined in Table 2-3 is another
useful statistic for comparing two-plume bulk properties. Two identical
but disjoint plumes have a correlation coefficient of -1; whereas, if
the two plumes are rotated into coincidence, the correlation is +1. For
two nonidentical plumes spanning the identical set of grid points, r
approaches unity as the respective values of the higher order moments of
the two-plume concentration distributions approach each other.
No comparison of absolute maximum concentration values is included
here; absolute concentration maxima are nearly always a near-field
phenomenon, and the way that the long-range transport models considered
here were applied does not necessarily yield accurate results at short
ranges (see Section 3.1).
44
-------
3. COMPARISONS OF MODEL BASE CASE RESULTS
3.1 Qualitative Comparisons for Single and Multiple Point Sources
As the mesoscale flow patterns are varied, the three models will
sometimes intercompare very well, and sometimes not as well. The two
meteorological episodes chosen represent some typical variations of
mesoscale flow patterns that can be expected to occur over multiday,
regional-scale simulations.
The first set of base-case model comparisons, Figure 3-1, depicts
the isopleths of regional-scale 24-hour average ground-level concentra-
tions of SO2 from the Navajo power plant, as calculated by each mode 1
for day 16? of Episode 1 (the highly zonal day). Figure 3-2 depicts
concentration isopleths of SOt, from the MESOPLUME and MESOPUFF runs. In
each of the concentration contour figures shown, the contour values
increase inwards from the outermost (lowest) value respectively, 0.1,
0.5, 1, 2, and 5 ug/m3.
The overall concentration patterns produced by each model are
consistent with the zonal flow field on this day. Near the source, the
results from MESOGR1D are appreciably different from the results of the
other two models. Among the factors contributing to these differences
are (1) the initial dilution of the point source through a grid volume
in the MESOGRID model, (2) the very different treatment of plume mixing
in the vertical, (3) inadequate sampling of the lateral plume distribu-
tion by MESOPLUME and MESOPUFF when the horizontal plume size is much
smaller than the grid spacing, and (4) the inability of the particular
contouring package used to realistically handle point singularities.
On the regional scale, however, the three models generally inter-
compare closely. It should be emphasized here that unlike the other two
models, MESOGRID has no explicit horizontal diffusion - evidence that
plume dispersion on regional scales (at least as represented by each of
these models) is not sensitive to small-scale plume diffusion.
In the next set of base case model comparisons, the 24-hour average
concentrations for S02 (Figure 3-3) and SO^ (Figure 3-4) from the Navajo
plant are depicted for day 169 of Episode 1 (the day with highly dis-
organized stagnating flow conditions, particularly in the southeast
corner of the computational area). By comparison with the other two
models, MESOPLUME appears to produce an anomalous high SO? concentration
region far downwind--possibly a result of the inverse wind speed depen-
dence of the plume model (see Section C.4.3), magnified in a region
where the flow rapidly stagnates. The MESOGRID algorithm is also
subject to this nonphysical behavior, but because pollutant material
distributed in some arbitrary manner within an individual grid cell is
effectively spread uniformly over the entire grid cell for the purposes
of plume sampling, MESOGRID does not show anomalous concentrations as
high as those of MESOPLUME. Of the three models, the MESOPUFF results
appear most realistic in the stagnating region.
45
-------
08V18? MESQPUFF
Day 18? MESOPLUME
Figure 3-1 Episode 1 Day 167 Single Source Dispersion Model Comparisons (S02)
46
-------
eMnwaHWUcawKHtnoanaGv m
D«V 187 MESOGRID
Figure 3-1 Continued
47
-------
Day 187 MESOf* IFF
D;n« 16 ?
Figure 3-1 Continued
48
-------
iKWCNMtfrMfc
Day 107 MESOPLUME
Day 167 MESOPUFF
Figure 3-2 Episode 1 Day 167 Single Source Dispersion Model
Comparisons (SO*)
49
-------
Figure 3-3 Episode 1 Day 169 Single Source Dispersion Model Comparisons (S02)
50
-------
Day 169 MESOGRID
Day 183 MESOPIUWE
Figure 3-3 Continued
51
-------
Figure 3-3 Continued
I
52
-------
Day 168 MESOPLUME
Day 169 MESOPUFF
-------
i»rtoHNrt as* #c
The Multiple source base case model comparison runs using the
Episode 1 meteorology are shown in Figure 3-5 which depicts, for
day 167, the 24-hour average S02 concentrations for a set of ten major
point sources in the Four Corners area. It should be emphasized that
this set of ten sources was selected from the larger inventory of major
SO2 sources in the area to illustrate various aspects of model multi-
source performance for a distribution of source sizes and source
locations.
In general, the regional-scale impacts are, as expected, dominated
by the largest point sources whose plumes remain individually discem-
able at distances of several hundred kilometers. Note that the
0.1 jjg/m3 contour for each model approximately envelops the middle third
of the region. The plumes are aligned in about the same direction, as
dictated by the uniformly zonal flow that day. The intercomparison of
model results shows most strikingly that, notwithstanding differences in
basic model algorithms, the overall regional-scale behavior of the three
models is very similar in this case. The 0.1 yg/m3 contour envelopes
closely coincide everywhere.
The next set of model multisource comparisons, shown in Figure 3-6,
correspond to day 169 of Episode 1, the day of disorganized stagnating
flow. As is the case for the multisource results for day 167, the
0.1 ug/m3 contour envelopes compare very well on the regional scale.
But, as also seen in the single-source results for day 169, there are
important small-scale discrepancies evident in the concentration fields
produced by the three models under such adverse flow conditions. As an
alternative to the two-dimensional concentration field representation,
the three-dimensional perspective plots in Figure 3-6 provide a powerful
comparison of the magnitudes of impact computed by the three models.
These plots clearly illustrate, for example, that MESOPUFF and MESOPLUME
produce large-amplitude concentration 'spikes' in the near-field, but
HESOGRID does not.
Episode 2, it will be recalled, is characterized by strongly
recirculating flow situations, corresponding to the passage of two
migratory synoptic-scale lows throughout the region of interest.
Figure 3-7 displays, for each model as previously, the 24-hour average
model results for SO2 for days 124 to 127 for the Navajo plant. Again,
the gross consistency of regional-scale model results is apparent,
despite the great differences in modeling approaches. Small-scale
features also intercompare quite well on days 124 and 125 when the flow
is quite smooth and steady. However, for the more complex flow patterns
of days 126 and 127, there are important discrepancies in the higher
order moments of the plume distributions. It is believed that these
discrepancies are chiefly attributable to differences in the way each
model behaves in the recirculating flow region around the center of the
cutoff cyclone and in the stagnating flow region near the center of the
cyclone.
54
-------
C»*VI»0**'NTAi. «f S£A»C>« » TfCHNCXOO* SC
Day 167 MESOPUFF Day 167 MESOPIUME
Figure 3-5
Episode
(S02)
1
Day 167 Multiple Source Dispersion Model Comparisons
55
-------
*feoNNtatcXs-. #c
56
-------
tNMVMItNMLf«MM>*ftHaMNteaG« «n£
-------
Day MESOPUFF
Day 169 MESOPtUME
-------
mtmatmmHmMMOHtnomikOBi «c
D«v 109 MESOGRID Day 169 Ml .11:
Figure 3-6 Continued
59
-------
Day W MESOPUFF D«»v W MESOGRiD
Figure 3-6 Continued
60
-------
MESOPLUME 10 SOURCES S02
DRY 167
Figure 3-6 Continued
61
-------
t»•**>** HE«C
MESOPLUME 10 SOURCES S02 DRY 169
Figure 3-6 Continued
62
-------
»v,Mfc. «t • • r«M». 06 > «
MESOPUFF 10 SOURCES S02 DRY 167
Figure 3-6 Continued
63
-------
MESOPUFF 10 SOURCES S02 DRY 169
t <*>
*
%
I
£
fl£ a
V
0 «
1
S*
I"
Figure 3-6 Continued
64
-------
MESOGRID 10 SOURCES S02
r.'.vceK-VjMf RTSf 4 TfOM^cxoGv -v
DRY 167
Figure 3-6 Continued
65
-------
MESOGRID 10 SOURCES S02 DRY 169
a
a"
Figure 3-6 Continued
66
-------
Day 124 MESOPLUME
Dav 124 MESOPLUME
Day 124 MESOGRID
Day 124 MESOPLUME
Figure 3-7 Episode 2 Day 124, 125, 126, and 127 Single Source Dispersion
Model Comparisons (SO.,)
67
-------
Day 124 MESOPUFF
Day 124 MESOGRID
Figure 3-7 Continued
-------
EKvt»ONM£*?*t HSWO- * TtOMMUtOO* *C
Day 12S MESOGRID Day 125 ME50PLUME
{
\
Day 125 MESOPUFF Day 125 MESOGRID
Figure 3-7 Continued
69
-------
Day 126 MESOPUFF Oav 128 MESOPLUME
I 8
-
\
\
V ,
. . .
; i
' /\J^r J
Day 126 MESOGRID Day 126 MESOPLUME
I 1
R I
Figure 3-7 Continued
70
-------
Day 126 MESOPUFF
Day 126 MISOURID
Figure 3-7 Continued
71
-------
Day 12? MESOGRID VESOPi UME
Figure 3-7 Continued
72
-------
{*.. r»ONM?S**. *(<%£*» '•*, •» •«*».¦». CO V
3.2 Quantitative Comparisons of Model Results
As described in the previous section, the qualitative similarities
in gross model performance are encouraging evidence for the basic
validity (or at least consistency) of each of the mesoscale modeling
approaches. Qualitatively, however, there are clear differences in the
near vicinity of sources; and marked discrepancies in small-scale
features on the regional scale. When quantitative analyses are used to
compare model results, a picture consistent with the qualitative con-
clusions emerges. Figures 3-8 and 3-9 illustrate the correlation coef-
ficient (r) intercomparing the results of all three sets of single-
source model runs for each day during each episode.
In general, the MESOPUFF and MESOPLUME models tend to correlate
best; for zonal and recirculating mesoscale flows, r ? 0.9, but for the
highly stagnating conditions encountered during days 168 and 169 of the
first episode, the correlation coefficient is as low as 0.56. The low
values of correlation coefficient result, in part, because MESOPLUME
predicts anomalously high concentrations in an area, about 500 km down-
wind from the source, in which the flow decelerates to stagnation con-
ditions. The rapid decrease in plume segment length (associated with an
inverse wind speed dependence) dramatically increases the concentrations
computed by MESOPLUME. MESOPUFF does not exhibit this dramatic behavior
because puff size is independent of change in wind speed.
The MESQGRID results are sensitive to the onset of stagnation
conditions (probably because the two-dimensional flow field used to
drive MESOGRID is locally divergent), as is indicated by the decrease in
correlation coefficient with onsetting stagnation. Qualitatively,
though, the hot spots that result in the MESOGRID model are not as great
in magnitude as those in MESOPLUME, perhaps because in MESOGRID a pollu-
tant mass is considered to have a particular subgrid cell distribution
during the dispersion calculation, but, for the purpose of concentration
field sampling the pollutant mass is considered to be spread uniformly
over that grid cell.
The MESOGRID model requires a divergence-free wind field if it is
to conserve strictly pollutant concentrations (in the absence of sink
terms). The base-case meteorological data were processed with MESOPAC
with a maximum allowable absolute wind field divergence of 10"1* s~1.
This tolerance level was chosen because a substantially more stringent
level (10 5 s 1) was found (1) to result in a loss of important sub-
synoptic detail in the flow pattern, and (2) to increase the costs of
running the MESOPAC program by about a factor of 4. Under most cir-
cumstances, the 10 ** s * maximum divergence level is sufficiently low so
that no significant effect is noticeable when the MESOGRID results are
compared to either of the other models. It is possible under some
conditions for the divergence term to be important (see Appendix D).
73
-------
"Test" Case vs, "Bast" Case
a MESOPLUME vs. MESOPUFF
~ MESOGRID vs, MESOPUFF
I I MESOGRID vs. MESOPLUME
0.5
167 168 1®
Day
Figure 3-8 Days 16", 168, and 169 Comparison of the Different Models
Results (Correlation Coefficient r")
1.0
123
124
125
Day
128
127
I Figure 3-9 Days 125, 124, 125, 126, and 127 Comparison of the Different
I Mode Is Results (Correlation Coefficient r)
74
-------
f nv <«CWMf N-? At »f <4? « * * f <' >¦**"* ':>-'* •**
Another useful quantitative criterion for model comparison is the
effective size of the plume, as measured by the number of grid points
with significantly nonzero concentrations. This measure is shown in
Figures 3-10 and 3-11 illustrating 24-hour average plume size for each
day of the two episodes. Note that (1) the MESOPLUME and MESOPUFF
models calculate plume concentrations out to three horizontal standard
deviations from the centerline and that (2) MESOGR1D will set to zero
any calculated concentration below a minimum value of 10 14 micrograms
per cubic meter (ug/m3). The fractional difference of number of points
in the plume (FDN) is similar in all cases, with the notable exception
of the MESOPUFF results for days 126 and 127 of Episode 2. This may be
due to the longitudinal diffusion in MESOPUFF, tending to spread pollu-
tant mass tangentially outward from a curved plume trajectory - the
resultant MESOPUFF plume then encompasses more grid points than those of
the other models. The results suggest qualitatively that horizontal
dispersion is dominated by plume meander, except under special flow
situations - note that the MESOGRID model, with no horizontal diffusion,
yields results for plume size very similar to those of the other two
models.
Figures 3-12 and 3-13 illustrate the time variation of the frac-
tional difference of the plume means (FDM) for each 24-hour average time
period during the two episodes studied. During the first episode, flow
stagnation conditions in the eastern part of the grid during days 168
and 169 result in an unrealistically high concentration region in this
area for the MESOPLUME model, as evidenced by the large negative values
of the FDM when MESOPLUME is compared to MESOPUFF. Concentrations as
large as 5 ug/m3 occur far downwind in the region of lowest wind speeds;
this occurrence is most likely the result of the inverse wind speed
dependence of the MESOPLUME model. The 24-hour average plume concen-
tration (in ug/m3) for days 167-169 are 0.39, 0.41, and 0.40 for
MESOPLUME and 0.27, 0.25, and 0.22 for MESOPUFF. During the second
episode, MESOPLUME predicts high concentrations in the low wind speed area
surrounding the center of the cutoff cyclone (days 126-127). The
MESOPLUME-MESOPUFF 24-hour average plume concentrations are very close,
however, in periods of zonal flows and moderate wind speeds during both
Episodes 1 and 2.
The comparisons for Episode 1 indicate that MESOGRID produces nearly
identical average plume concentrations as MESOPLUME, and somewhat higher
concentrations than MESOPUFF. The absence of horizontal diffusion in
the MESOGRID model is evidently not a significant factor on the regional
scales considered here. The FDM results for Episode 2 are again similar
for the MESOGRID and MESOPUFF models, with the exception of days 126 and
127. It is believed that the significantly higher 24-hour average plume
concentrations for the MESOGRID model result from an area of convergence
with low wind speeds.
75
-------
n MESOPUFF
I I MESOPLUME
t I MESOGRID
Figure 5-10 Episode 1--Comparison of the Number of Points in the Plume
MP Computed by Each Model
CL
z
500
400
Figure 3-11 Episode 2--Comparison of the Number of Points in the Plume
NP Computed by Each Model
76
-------
CVviBi-friMFNUt BtSCARCMA "f CMhlOviVW r^C
+1,0
2
o
-1.0
"Test" Case vs. "Base" Case
I I MESOPLUMEvs. MESOPUFF
QH MESOGRID vs. MESOPUFF
n MESOGRID vs. MESOPLUME
167
1%
Day
1®
Figure 3-12 Episode 1--Comparison of the Di f fei >nt Model Results
[Fractional Deviation of the Means (FDM)]
Figure 3-13 Episode 2--Compar.ison of the Different Model Results
„ [Fract ional Devi at ion of the Means (FDM)]
1
77
-------
f -SMf K
The model comparisons shown in these last two sections are believed
to be the only such set of regional-scale comparisons of plume, puff,
and grid models that have been achieved in an objective, systematically
neutral manner, with uniform control over the models' external environ-
ment (that is, their meteorological, input/output, and analysis routines).
These results give initial indications that for some regional-seale
simulation purposes, the choice among a plume segment model, a puff
element model, or a grid model may not primarily depend on the nature of
the model results--the results for each model might be very similar in
gross properties. Rather, the decision may be based on other issues
such as relative ;ccut ion cost or the regulatory acceptability of one
or another modeling approach. For applications where accurate knowledge
of the location and magnitude of concentration maxima is necessary, the
fine-scale model behavior may be the dominant issue, especially in the
near field of sources. In complicated flow situations the variable-
trajectory approach used by the three models is more realistic, and
therefore likely to provide more credible result s, than the conventional
straight-line Gaussian plume approach used for most regulatory appli-
cations at present.
The comparison of relative model performance will be explored
further in Sections 4 and 5. Issues such as model sensitivity to vari-
ations in input parameters or time resolution will further illuminate
the relative advantages and disadvantages of each model for various
meteorological and source inventory scenarios.
3.3 Model Simulations of Plume Fumigation
As the mixing depth oscillates during a diurnal cycle, pollutant
mass emitted from an elevated source may either be emitted within the
turbulent boundary layer (leading directly to near-source surface
impact) or may be emitted above the boundary layer (and not intersect
the surface until considerably downwind of the source). On the regional
scale, plumes may be transported aloft over long distances (many
hundreds of kilometers) during stable nighttime hours, only to be
rapidly mixed down to the surface as morning convection deepens the
mixing layer. Thus, the regional-scale predictions of each of the three
models will be affected significantly by the treatment provided for
plume fumigation. This section provides qualitative illustrations of
the treatment of fumigation by each of the models.
Figure 3-14 depicts the progression of 1-hour average ground-level
SO2 concentrations as predicted by the MESOPLUME model for the Navajo
power plant during Episode 1 from 7:00 P.M. [Mountain Standard Time
(MST) on day 166] through 1:00 P.M. on day 167. A set of mixing depth
fields at the corresponding times is depicted in Figure 3-15. These
sequences of figures can be interpreted as follows. After sunset, the
mixing heights decreased to a relatively small value (a mechanically
driven equilibrium value dependent on the free stream wind velocity).
When the mixing height became ress than the effective stack height,
78
-------
7 pm MESOPLUME
f«¥*R£JNMIsNt*l HfSEAWO* I IfCHNOlOGv Nt
1
I
1 am MESOPLUME
7 pm MESOPLUME
' "
i
i
¦ f ¦
I
; . W ]
r
—4 -- -
rJ
(
|
Figure 3-14 MESOPLUME Time Sequence of 1-hour Average SO, Concentrations
(7 pm MST Day 166 to 1 pm Day 167)
79
1
-------
7 am MESOPLUME
1 am MESOPLUME
Figure 3-14 Continued
80
-------
1 pm MESOPLUME
9am MESOPLUME
Figure 3-14 Continued
81
-------
Day 166 7 pm
f
\
Tt
w
\
| j
^ — "i ¦?
/ 4
' {
/
i
i
i j
\
r
• - -
— —
' %
! '¦
m
Day 167 1 am
Day 167 7 am
Figure 3-15 Time Sequence of Mixing Depths (7 pm MST Day 166
to 1 pm Day 167)
82
-------
»mWM#TM.I«§£AiC»4*'f04Nat®GV MI
Day 167 Sam
Day 187 1pm
Figure 3-15 Continued
83
-------
newly emitted SO; was no longer mixed to the ground; the SO,, previously
mixed to the surface during the day continued to be advected, diffusing
slowly. Figure 3-16 illustrates the extent that dispersion and deple-
tion of SO2 within the plume reduces the average mass in the portion of
the plume contacting the ground surface during nocturnal regional-scale
transport.
After sunrise, the convective mixed layer grew rapidly and
entrained the entire plume by 9:00 A.M. local time (Figure 3-14,
9:00 A.M. MESOPLUME). The plume was rapidly mixed downward to the
surface. The highest concentrations of the day occurred at this time
while the plume was trapped within a growing but shallow mixing layer
(see Figure 3-15, 9:00 A.M. MESOPLUME). As the mixed layer height
continued to grow, the emissions < sre diluted through an even deeper
layer. After sunset, the mixed layer height collapsed, and the entire
sequence of nocturnal plume transport began anew.
As shown in Figure 3-17 the MESOPUFF model displays virtually
identical plume fumigation behavior as MESOPLUME.
For the MESOGRII) model, however, fumigation is much less dramatic
(see Figure 3-18) because vertical dispersion is determined by the K-
profile, the vertical resolution (currently three vertical layers), and
the depth of each vertical layer (in the present version of MESOGRID,
the cell value assigned to K, is set to a very small number whenever the
height of the grid cell exceeds the mixing height). In the MESOGRID
fumigation sequence, the "old" surface plume is advected and depleted, as
previously. Some continuous mixing to the surface does occur, however,
and the onset of fumigation is much less abrupt, because the continuous
profile of vertical diffusivity is specified by only four discrete
values and usually only the lowest two values are relevant to the fumi-
gation process. With the onset of convection in the early morning, Kz
begins to increase anew as the mixing depth increases and the stability
decreases; this increase leads to a more rapid but not immediate
transfer of mass to the surface.
84
-------
Local Time of Day
Figure 3-16 MESOPLUME Mean Plume Concentration Cycle (for 1-Hour Averages)
85
-------
7 pm MESOPUFF
Figure 3-17 MESOPUFF Time Sequence of 1-Hour Average SO- Concen-
trations (7 pin MST Day 166 to 1 ps Day 167)
86
-------
Figure 3-17 Continued
87
-------
•mmmcnui * am («c
1 pm MESOPUFF
9am MESOPUFF
Figure 3-17 Continued
88
-------
fciwmcMMf**'*, HtTu'zmx -*v» #c
7 pm MESOGRID
1 am MESOGRID
7 pm MESOGRID
Figure 3-18 MESOGRID Time Sequence of 1-Hour Average SO, Concen-
trations (7 pa MST Day 166 to 1 pm Day 167)
89
-------
fMWO»«»jrAi. & rtu^ejtoct i«c
7 am MESOGRID 1 am MESOGRID
9am MESOGRID
7 am MESOGRID
Figure 3-18 Continued
90
-------
1 pm MESOGRIO
9am MESOGRID
Figure 3-18 Continued
91
-------
FNWUNU|-K'*i *KM ft n < wf** ,». -u
4. RESULTS OF MODEL SENSITIVITY TESTS
This section presents some of the highlights of the more than
100 multiday model sensitivity simulation experiments performed in this
study. Recalling the sensitivity test protocol described earlier in
Section 2, Table 2-1 identifies the kinds of sensitivity tests con-
ducted. Model sensitivities to differences in mesoscale flow pattern
have already been discussed.
4.1 Sensitivity to Decay and Dry Deposition Processes
Plume depletion processes are important factors on regional scales.
Each model permits two different depletion processes: (a) dry deposi-
tion of S02 and SOlj, modeled using the "deposition velocity" concept;
and (b) conversion of S02 to SO^, assumed to be approximated by a linear
decay mechanism. The MESOPLUME and MESOPUFF model sensitivities to
variations in deposition and decay rates were evaluated in detail over
intervals representing the 1ikely uncertainty in these parameters. The
much higher cost of running the MESOGRID model meant that fewer MESOGRID
sensitivity tests could be performed. Table 4-1 displays the pairs of
values of the plume depletion parameters used in the sensitivity tests.
One effective way to display the results of an air quality model's
sensitivity to variations in model parameters is the format introduced
by Hilst (1970). Figure 4-1 (4-2), depicts the sensitivity of the
MESOPLUME (MESOPUFF) model to variations in plume decay rate and to
variations in SOj/SOt, deposition velocities in the Hilst format. The
abscissa variable is the relative fractional difference in value chosen
for the perturbed model run parameter compared to the parameter value
chosen for the base case model run. For example, in one case, the
perturbed case value of decay rate is chosen to be ten times larger than
the base value, and the abscissa is -9; if instead the perturbed case
value is chosen as one-tenth that of the base value, then the corres-
ponding abscissal value is +0.9. The ordinate values in this figure
describe the average fractional difference (AFD) in 24-hour average
ground-level SO2 concentration values computed in the base case and
perturbed case runs, averaged over all grid points and over the last
3 days of Episode 1. Each of the points plotted in the figure repre-
sents , therefore, a separate 5-day run made to display this variation.
As seen from the figure, an approximate five-fold increase in the
decay rate, or in the deposition velocities, leads to approximately a
two-thirds decrease in the average relative fractional difference
concentration field; a five-fold decrease in the decay rate, or in the
deposition velocities, corresponds to only about a one-third increase in
average relative difference field. This asymmetrical relative model
response to variations in input parameters is very typical. MESOPLUME
(MESOPUFF) also shows a somewhat greater sensitivity to changes in decay
92
-------
t*z¥tm>MWWt«t*we"*t*CMNOtOG* **.
TABLE 4-1
VALUES OF THE DEPOSITION VELOCITIES AND DECAY
RATES USED IN THE SENSITIVITY TESTS
Decay Rate Coefficient, Deposition Velocity, fro/s)
MESOPLUME and MESOPUFF
kl
JS
-5.56
X
10-5
Vd
fso2, sop
ss
(0.01, 0.001)*
kl
ss
-2.22
X
10~5
fl
Vd
(S02, sop
sr
(0.01, 0.001)*
kl
=
-5.56
X
10 *
Vd
(S02, S0°)
=
(0.1, 0.01)
kl
=
-5.56
X
io"b*
*•
Vd
(S02, sop
sr
(0.0316, 0.00316)
kl
s
-5.56
X
10 *
Vd
(S02. sop
as
(0.01, 0.001)*
kl
=
-5.56
X
10~6*
£L
Vd
(S02, sop
=
(0.00316, 0.000316)
kl
St
-5.56
X
10 *
fL
Vd
(S02, sop
s
(0.001, 0.0001)
kl
5=
-3.89
X
10"6
Vd
(S02, sop
35
(0.01, 0.001)*
kl
ss
-5.56
X
io"7
Vd
(S02, sop
35
(0.01, 0.001)*
MES0GRID
kl
=
-5.56
X
10"5
vd
CS02, sop
SS
(0.01, 0.001)*
kl
SS
-5.56
X
io"6*
vd
(SO,, sop
ss
(0.1, 0.01)
kl
SS
-5.56
X
io"6*
vd
(S02, sop
=
(0.01, 0.001)*
-6~
kl
ss
-5.56
X
10 *
vd
(S02, S04)
ss
(0.001, 0.0001)
k.
ss
-5.56
X
io-7
Vd
(S02, sop
ss
(0.01, 0.001)*
*Base case values.
93
-------
, ». , :»*, aWf-^Tfts, *H I. f t;run*.* OU- -V
MESOPLUME SENSITIVITY TESTS
(DAYS 167-169)
PBASE = Base Value
Pjest = (Perturbed) Test Value
«
i
I
Figure 4-1 MESOPLUME Sensitivity to Decay Rate and Dry Deposition
Variations (Days 167-169)
94
-------
g H£S£A»tMS TtTXNCXOC.v
rate than to changes in deposition velocity. This sensitivity results
fro* the diurnal fumigation cycle and the models' treatment of dry
deposition and decay. For a given deposition velocity, the amount of
pollot ant removed by dry deposition depends on the ground-level con-
centration of the pollutant and the mixing depth. For the base case
choices of SO decay and deposition rates, the decay term dominates the
deposition term whenever the plume is mixed through more than about
1,800 m in the vertical (because of the inverse-M dependence of the
deposition term); when the plume is mixed through less than 1,800 m, dry
deposition is the dominant process in the model. Dry deposition,
however, occurs only when the plume is in contact with the ground,
whereas decay processes occur continuously, even when the plume is above
the mixed layer.
As seen in the previous section, the mixed layer is generally
shallowest at night; new emissions may remain above the mixed layer
until the plume fumigates the next morning, prohibiting dry deposition
of new emissions for up to 12 or more hours each day. After plume
fumigation, when dry deposit ion does occur, the mixing height values
usually increase to greater than 1,800 m in the spring and summer in the
Four Corners area, and decay processes dominate the removal.
Because the MES0PUFF model uses the identical decay and deposition
algorithms, it shows virtually the same sensitivities to systematic
variations in kj, vj(S02), and vjfSO^) as MES0PLUME (Figure 4-2). The
more 1imited set of experiments with MES0GRID show si milar results.
Figures 4-3 and 4-4 display another way in which to compare model
performance (in this case, MES0PLUME) with different choices of input
parameters. The contours in each figure describe the fractional dif-
ference, (Cg-C.j.)/C , in concentration field (FDF) values between the
24-hour average ground-level SO2 concentration values in the base case
(C ) and the test case (C_). In each case, grid point concentration
values less than 0.1 ug/m* were ignored. The first pair of figures
(Figures 4-3a and b) show the effects of a ten-fold increase (vj = 0.1),
and a ten-fold decrease (v^ = 0.001) in the values of S02 deposition
velocities, relative to the base case (vj = 0.01) for the Navajo single-
source case on day 169 (stagnation flow). In both instances, the
fractional differences increase with distance from the source and
indicate greater deviation from the base case for the choice of = 0.1
than for the choice of vj = 0.001. In contrast to Figures 4-3a and b,
Figures 4-5c and 4-3d illustrate the fractional difference field for the
cases vj = 0.1 and v^ = 0.001 on day 167 (strongly zonal flow). It is
noted that the fractlonal di fferences in this zonal flow situation, in
which plume segments have shorter grid residence times, are considerably
smaller than those of day 169. For example, for the choice of vj = 0.001,
there are test case values in the far-field that are more than twice
those of the base case on day 169, whereas on day 16? there are no test
case values even 50% greater than those of the base case.
95
-------
MESOPUFF SENSITIVITY TESTS
(DAYS 167-169)
Figure 4-2 MESOPtSH; Sensitivity to Decay Rate and Dry Deposition
Variations (Pays 167-169)
96
-------
tWMUMUC***! «*«*#*»•* 'HXNUl.».* *fc
¥,j - 0.1 m/s
Day 169
Figure 4-3a MESOPLUME Fractional Difference Field Comparison for
v, (SO.,) » 0.1 m/s Case vs. Base Case (Day 169)
a i
vj - 0.001 m/s
Day 169
Figure 4-3b MESOPLUME Fractional Difference Field Comparison for
v, (SO,) = 0.001 m/s Case vs. Base Case (Day 169)
d 2
97
-------
vd - 0.001 m/s
Day 167
Figure 4-3c MESOPLUME Fractional Difference Field Comparison for
vd (S02) = 0.1 m/s Case vs. Base Case (Day 167)
Figure 4-3d MESOPLUME Fractional Difference Field Comparison for
vd (S02) * 0.001 m/s Case vs. Base Case (Day 167)
98
-------
This sequence of figures emphasizes that the meteorological regime
can exert a great influence on the sensitivity of the model to different
choices of deposition velocity. Although only the MESOPIUME model is
illustrated here, the other models also show comparably wide variations
in deposition velocity sensitivity as the meteorology changes.
The next set of figures (Figures 4-4a through 4-4d) show the
corresponding fractional difference results of the MESOPLUME model for
ten-fold increases and decreases in decay rate kj, relative to the base
case, both for the stagnating condi t ions (day 169) and the zonal con-
ditions (day 167). The strong dependence on meteorological flow regime,
evident in the deposition velocity tests, is also seen here.
4.2 Sensitivity to Mixing Depth Variations
The sensitivity of the models to systematic variations in the
mixing depth was evaluated by varying the mixing height by simple multi-
plicative factors--the values used were 0.5, 0.67, 1.5, and 2.0. If the
plume is assumed to be entrained by the mixed layer under all base and
test conditions and if dry deposition is ignored, a doubling of the
mixing depth corresponds to a halving of the average concentration
(I-'DM = 0.5) with perfect correlation of base and test cases (r = 1.0).
Likewise, a halving of the mixing depth means a doubling of the average
concentration (EDM = -1.0) and again r = 1.0. In actuality, the
behavior of I-DM and r are very different, as illustrated in Figure 4-5
for the MESOPLUME 24-hour average SO - concentrations (Navajo single
source). For Episode 1, the FDM values seem only weakly correlated with
the mixing depth multiplier, and r becomes much less than unity for the
largest mixing depth multiplier. These results may be explained by
differences in the nocturnal behavior of the plume in the base and test
cases. If in the base case the plume is above the mixed layer, then
increasing the test case mixing depth may cause plume entrainment;
therefore, rather than the average plume concentration halving when the
mixing depth doubles, the ground-level concentration could actually
increase, which did occur on day 168.
The fractional difference fields (Figure 4-6 a, b, c, and d) on
days 167 and 169, show that the sensitivity of the models (here,
specifically, MESOPLUME) to mixing depth variations is highly episode-
specific. For the examples shown, raising (lowering) the mixing depth
lowers (raises) the plume concentrations in the very far field, because
plume elements far from the source will have experienced at least one
diurnal mixing depth cycle and will have therefore been entrained.
[This process occurs even in the case of the 0.5 nixing depth multi-
plier, because in the Four Corners region the summer mixed layer is
generally very deep (recall Figure 2-3)]. In contrast, no such definite
relationship can be assumed in the near field because of the episode-
specific nature of plume entrainment.
99
-------
smmom&mtt «Ese*RCw «, !io«a og* m
k1 - 5,56 E-5
Day 169
Figure 4-4a MESOPLUME Fractional Difference Field Comparison for
k1 = -5.56 x 10-5 Case vs. Base Case (Day 169)
k, = -5.56 E-7
Day 1&)
~[
--
1
!
i
i
I
s
!
, "t
/' !
Vl\ i
/ L /"
_.A
r
rJ
1
X
\ f
\
(
j
/
Figure 4-4b MESOPLUME Fractional Difference Field Comparison for
kj = -5.56 x 10"7 Case vs. Base Case (Day 169)
100
-------
k1 - -5.56 E-5
Day 16?
Figure 4-4c MESOPLUME Fractional Difference Field Comparison for
kj = -5.56 x 10"5 case vs. Base Case (Day 167)
—
i
]
I
i
* X"'
1
! /
OS |
/ J
L /
/
i
r
<-:i''...
ry
i
\
i
i
i
Fieure 4-4d MESOPLUME Fractional Difference Field Comparison for
kj = -5.56 x 10'7 Case vs. Base Case (Day 167)
101
-------
(MVWDNMIWTAt . *,
Mixing Depth Multiplier
Figure 4-5 MESOPUFF Sensitivity to Mixing Depth Variations (Days 167-169)
102
-------
PBLOEL - 2.0
Day 167
1
4
Figure 4-6a MESOPLUME Fractional Difference Field Comparison for Mixing
Depth Multiplier =2.0 Case vs. Base Case (Day 167)
Figure 4-6b MESOPLUME Fractional Difference Field Comparison for Mixing
Depth Multiplier ™0.5 Case vs. Base Case (Day 167}
103
-------
PBLDEL - 2.0
Day 169
Figure 4-6c MESOPLUME Fractional Difference Field Comparison for Mixing
Depth Multiplier = 2.0 Case vs. Base Case (Day 169)
PBLDEL = 0.5
Day 169
Figure 4-6d MESOPLUME Fractional Difference Field Comparison for Mixing
Depth Multiplier = 0.5 Case vs. Base Case (Day 169)
104
-------
SKviACSMFMAi *<
4.5 Sensitivity to Systematic Variation of Wind Direction and
Crosswind Growth Rate
A series of sensitivity tests were performed with the MESOPLUME
and MESOPUFF models to compare the relative importance of (1) horizontal
crosswind plume growth and (2) systematic wind direction variations.
This set of sensitivity tests was performed by adding a constant angular
shift factor to the base case wind di rection at each grid point and,
simultaneously, using a multiplicative scaling factor to adjust the
horizontal, time-dependent plume growth rates [used for the plume or
puff elements at travel distances beyond 100 km only]. (Both MESOPI.UME
and MESOPUFF use the PGT sigma curves for travel distances _< 100 km.)
These tests have three specific objectives:
(a) they investigate the sensit ivity of the model to changes in
plume growth rate. The results suggest that plume meander
(resulting from large-scale fluctuations in the mesoscale wind
field) dominates small-scale turbulent dispersion (as repre-
sented by plume diffusion coefficients). If this is generally
true, more exact parameteri zations of small-scale dispersion
may not be necessary.
(b) they try to determine how much systematic wind direction
variations influence the orientation and shape of a given
plume, because large systematic errors can occur when the
driving wind field is initialized. For example, the direction
of the 700-mb station wind used in this study as input to
MESOPAC could differ significantly from the wind direction
most representative of plume altitude for the particular
dispersion problem at hand. This wind direction initiali-
zation error could produce concentration patterns very
different from the real situation.
(c) they examine the relative importance of plume growth rate vs.
wind shift. This could be useful in determining whether, for
example, the effect of lateral plume growth is great enough to
override the effect of a systematic wind shift. The matrix of
sensitivity cases run with each model for meteorological
Episode 1 is shown in Table 4-2. For each model, the base
case run is defined by n£ wind shift and a plume growth rate
of 0.S m/s.
Note that the sensitivity of the models to random variations in wind
direction (such as those caused by improper initialization of the
subsynoptic scale wind field) was not studied, at least not directly.
Because the primary result of random wind variation is the enhancement
of plume diffusion, it was felt that the effect of imposing random wind
variations can be adequately represented by increasing the crosswind
plume growth rate.
105
-------
TABLE 4-2
MATRIX OF WIND SHIFT AND GROWTH RATE SENSITIVITY CASES
Angular Wind Shift
Plume Growth
Rate (m/s) 5.625° 11,25° 22.5°
0.0 x x x
0.5 xxx
1.0 x x
106
-------
f NVtfKlNMIST* flCSC AftCMft rtCMNCfcOC.* **
The two statistical measures found to be most useful in evaluating
the results were the correlation coefficient (r) and the fractional
deviation of the plume means (FDM). The correlation coefficient, r, is
useful because it is a measure of the degree of overlap among the con-
centration fields of two plumes. The FDM statistic is also useful
because it is a measure of the gross first-order moment difference
between plume structures, without regard to differences in their spatial
orientations.
The next set of figures shows the results for the zonal day (167)
and stagnating day (169) , both for the MESOPLUME mode 1 (Figures 4-7 and
4-8) and for the MESOPUFF model (Figures 4-9 and 4-10). Each figure
contains subject ivelv sketched curves of approximately constant corre-
lation coefficient (r) (shown as sol id lines) as we 11 as the FDM (shown
as dashed lines). These curves are seen to be more or less orthogonal,
because two plumes with very comparable mean values (FDM ^0) and high
correlation coefficient when the directions coincide (r >0.75), rapidly
become uncorrelated with only very modest angular shifts in wind
direction.
As illustrated next, the sensitivity to variations in wind direc-
tion is quite dependent on the meteorological episode. A comparison of
Figure 4-7 with Figure 4-8, and Figure 4-9 with Figure 4-10, shows that
during the moderate wind speed zonal flow conditions (day 167), a small,
systematic wind direction shift results in a test plume uncorrelated
with the base case plume; whereas, the light variable winds during the
stagnating flow period (Day 169) lead to results that are appreciably
less sensitive to the imposed variations in wind direction.
Wind direction shifts had very little effect on the FDM, especially
during the stagnating flow conditions. Systematic variations in the
wind direction do cause plumes to pass through different sections of the
grid and, therefore, experience somewhat different mixing height and
stability fields, but, in this case, the slight FDM dependence on wind
direction variations probably reflects the fact that these fields are
homogeneous.
The sensitivity of the MESOPLUME and MESOPUFF models to variations
in the plume growth rate, doy/dt, also shows a dependence on the meso-
scale flow pattern. For growth rates less than the base case (0.S m/s),
the FDM is much more sensitive to changes in the growth rate, when flow
conditions are zonal (relatively steady), than when stagnating (vari-
able). The results also indicate that once a nominal growth rate is
reached, further increases of the growth rate do not substantially
change the mean plume structure. The plume growth rate appears to be
important only until the plume reaches a size large enough to be influ-
enced by the mesoscale flow structure. Thereafter, the plume dispersion
is dominated by the larger (meso-) scales of the flow; horizontal plume
dispersion from small-scale turbulent diffusion is no longer of primary
107
-------
'I'.***.'-.* »*C
0.5
1.0
dov
^ (m s"1)
Figure 4-7 MESOPLUME Sensitivity of r (—) and FDM (---) to Systematic
Variations in Wind Direction A0 and Plume Growth Rate do /dt
(Day 16?) y
108
-------
HMMOt MM.MKM ftCWttOIJ*
0.5
dOy
_F
(m s-')
1.0
Figure 4-8 MESOPLUME Sensitivity of r (—) and FDM (—) to Systematic
Variations in Wind Direction A6 and Plume Growth Rate da /dt
(Day 169) y
109
-------
22.5"
<1
11,25°
5.625°
0.5
1.0
dov
dT1"""1
Fieure 4-9 MESOPUFF Sensitivity of r ( ) and FDM ( ) to Systematic
Variations in Wind Direction A8 and Plume Growth Rate dew dt
(Day 167)
110
-------
22.5*
as
* 11.25"
5 625°
0.5
d
aou
-g Cm s 'J
1.0
Figure 4-10 MESOPUFF Sensitivity of r ( ) and FDM (---) to Systematic
Variations in Wind Direction £0 and Plume Growth Rate do /dt.
(Day 169) y
111
-------
importance. The correlation coefficient is remarkably independent of
variations in dcv/dt - even in the extreme case of no plume growth
beyond 100 km (dby/dt = 0), the correlation coefficient is at least
0.75.
Thus, for regional-scale transport and dispersion as simulated by
these models, the specific nature of the wind field is critical -
chiefly because the mesoscale flow features appear to play the dominant
role in plume dispersion on these scales.
4.4 Model Sensitivity to Variations in Wind Speed
Just as in the case of wind direction, the initialization scheme
can introduce appreciable systematic errors in transport wind speeds.
To understand the effect of a systematic variation of wind speed on the
plume or puff model results, consider the idealized situation repre-
sented in Figure 4-11, Case A. Assume that the wind field is hori-
zontally uniform and steady, and that it takes t hours for the emissions
from a point source to reach the edge of the grid. Doubling the wind
speed at each grid point will halve the travel time to the end of the
grid (Case B). The mass of pollutant integrated over the entire plume
is twice as large for Case A as for Case B (assuming no removal),
because the plume in Case A contains t hours of emissions (assumed
constant) and the Case B plume contains only t/2 hours of emissions.
The centerline concentrations, however, should be approximately the
same in both cases; this may be understood as follows. The MESOPLUME
model assumes the centerline concentration C is inversely proportional
to uo :
y
C
-------
Case A
OaseB
Figure 4-11 Effect of Wind Speed on Idealized Plume Shape
113
-------
The wind speed sensitivity tests were conducted by adjusting the
wind speed at each grid point with a multiplicative factor in each model
rim. The factors used were 0.5, 1.5, and 2.0, relative to the base
case. Mote that, for these sensitivity tests, the MESOPAC-generated
¦echanical mixing depths and stability classes were not recomputed
whenever the wind speeds were changed.
Figures 4-12 a and b display the correlation coefficient r as a
function of the wind speed multiplier for wind speed sensitivity runs
MESOPLUME and MESOPUFF. The results are similar for both models;
meteorology-dependent behavior is seen in the day-to-day variations of
the results. Increasing or decreasing the wind speed by a factor of two
causes about a two-fold reduction in the correlation coefficient. As
was seen in the analysis with the correlation coefficient, day-to-day
variations of the fract ional difference of the means (Figures 4-13 a
and b) are highly episode-dependent. Large FDM values are sometimes
produced when a change in the wind speed mu11 ip1ier leads to a change in
the plume trajectory that exposes the plume to radica 11y di fferent wind
vectors, mixing depths, and stabilities. But even with the occurrence
of these large deviations, no well-defined FDM/wind speed multiplier
trend can be extracted from the analysis. Implicit dependencies on the
spatial and temporal variations in meteorology override any explicit
trends in the data.
4.5 Sensitivity of MESOGRID to Variations in
A series of tests has been performed to determine the sensitivity
of MESOGRID to variations in k. values. The concern here is to illus-
trate how well test case concentration fields, produced by altering some
or all of the K, values, reproduce a base case concentration field
generated using a predetermined K, criterion. The base case K,
criterion was chosen such that the MESOGRID "requested time step" is
one hour. Thus, all overly-large values of K2 are truncated to the
threshold value, which ensures numerical model stability in the vertical
for a one-hour time step (see Section C.4.4 for a complete discussion).
Two methods of varyir? K2 field values were adopted in the present
sensitivity study. The first method involved the direct reduction of K,
values by a K, multiplier (less than unity). In Figure 4-14, this
method was used to produce the portion of the sensitivity curve that
lies to the ri.ht of the ordinate axis. It is seen that reducing all K,
values by as much as 90% causes a significant reduction in the average
concentration of the plume. Figure 4-15a supports this contention--
when K2 values are reduced by an order of magnitude, ground-level
concentrations are decreased [Fractional Difference (FD) field values
are positive) everywhere except in the very far field. (Note that, in
Figures 4-15a and b only, MESOGRID plume concentrations as small as
114
-------
. \. aif <4 ««*¦>* 'f 'v *G • •«#
1.0
*- 0.5
0.5 0.67 1.5 2.0
Day 167 a
fvl Day 168
Day 169
Figure 4-12a MESOPLUME Sensitivity of r to Wind Speed Multiplier (a)
Variations (Days 167-lt>9)
IIS
-------
+0,5
S
a
-0.5
Wflh Day 167
1 1 Day 168
M Day 169
Figure 4-13a MESOPLUME Sensitivity of FDM to Wind Speed Multiplier (a)
Variations (Days 167-169)
Figure 4-13b MESOPUFF Sensitivity of FDM to Wind Speed Multiplier (a)
| Variations (Days 167-169)
i
116
-------
Figure 4-14 MESOGRID Sensitivity to Variations in (Day 167)
117
-------
a ,:ri> #*_
-4
('
Figure 4-15b MESOGRID Difference Field Comparison (Limited Grid) for
Requested Time Step » S Minutes Case vs. Base Case (Day 167)
118
-------
eWRONMENTAi R£S6*RC«*tFChMDiOG» «NC
10"" lig m"3 are included.) Apparently, reducing K, values directly
decreases the rapidity of vertical mixing; pollutant material in the
upper two grid levels will take longer to diffuse to the lowest
layer, thus resulting in lower surface concentrations. But, reducing
the surface concentration reduces the rate of removal by dry deposition,
leaving more pollutant mass in the atmosphere; thus, after long travel
times, ground-level concentrations in the "reduced-K," cases may
actually become greater than those in the base case.
The second method of varying K- involved the indirect modification
of overly large values of K: by reducing the "requested time step."
This method was used to produce the portion of the sensitivity curve to
the left of the ordinate axis in Figure 4-14. In pract ice, in part s of
the Four Corners study, K, values up to 1,000 m2/s were found to be
produced in strongly convective conditions by the algorithm discussed
in Section C.4.3, These very large values, unless truncated, would
force the model to use a time step even smaller than 5 minutes. Values
of K-. that exceeded the limits of the requested time steps (respectively,
30 minutes, 15 minutes, and 5 minutes) were automatically truncated.
Because of the truncation, there is some uncertainty in the values of
K-. at each grid point. The FDM fields also reflect the effects of
truncation, which are expressed as estimated error bases on the K2 values.
It is curious that allowing K, values to increase apparently
results in reduced surface concentrations—one would rather expect
stronger vertical mixing to increase the rate at which material from
above is mixed to the surface. The fractional difference field in
Figure 4-15b compares the 5-minute "requested" time-step test case to
the 1-hour "requested" time-step base case. This comparison is a useful
adjunct in determining the cause of the Kz versus plume concentration
relationship. (Note that a somewhat smaller computational grid has been
used here to reduce the cost of a MESOGRID simulation for a 5-minute
requested time step.) It can be seen that near the plume centerline,
base case values exceed test case values (FD is positive), whereas away
from the centerline, the opposite is true (FD is negative). Apparently,
when a much shorter time step is used, pseudodiffusion effects resulting
from multiple time steps are appreciable and tend to spread the plume
laterally. For a 5-minute "requested" time step, this effect is enough
to reduce mean plume concentrations by about one-third.
In the example shown, the resultant tradeoff of the two competing
mechanisms, horizontal pseudodiffusion and vertical explicit diffusion,
yields the most conservative plume concentrations when a "requested"
time step of approximately one hour is used.
119
-------
5. MODEL COST/SENSITIVITY ISSUES
5.1 Overview
The model sensitivity information presented in Section 4 is
necessary but not sufficient to base a choice of the respective modeling
techniques most appropriate to a given situation - particularly in the
context of regulatory applications, serious issues of relative time and
cost are also encountered. The use of an expensive, complicated method
with extensive time and cost requirements should be questioned if a
simpler, less costly method can produce similar results for the given
application.
To gain more insight into these issues, a set of sensitivity tests
was performed to study individual model performance vs. execution time
requirements and to compare the cost/sensitivity results for the three
models. The following questions, among others, were addressed.
• Which parameters control the cost of a given model simulation?
What are the cost factors? How do these vary with meteorology,
emissions inventory, and mode 1 run parameters?
• How do the model simulation costs compare for identical input
conditions?
• Which factors and parameters can be modified or relaxed to
produce substantial savings in model simulation costs, yet
with only minimal deviations in output concentration fields?
• How do the costs of running each of the three regional-scale
dispersion models compare to the costs of collecting and
preprocessing the meteorological input data and the costs of
postprocessing and analyzing the model output concentration
fields?
5.2 Model Execution Cost Formulas
5.2.1 MESOPAC Run Time
Execution cost of the MESOPAC model (Benklev and Bass 1979c) is
determined principally by:
• the number of time steps,
• the number of grid points, and
• the wind field divergence "tolerance factor".
120
-------
»€S€*«Cm* TptMNCXOGv
-------
• the number of sources S, and
• various meteorological factors that control the average
number of puffs resident on the grid.
The MESOPUFF model allows for the specification of puff release and puff
sampling substeps (see Section B.9). These substeps can be chosen to
ensure that discrete puffs emitted consecutively from a source will
overlap sufficiently to simulate a continuous plume. For a single
source, the following relationship for MESOPUFF execution time was found
useful for estimating run time for typical meteorological conditions in
the Four Corners area:
t(sec) = 20 * O.SN + 0.003 (n -1) N2 + 0.002 (n -11 N2 (5-2)
P s
Note that while both ru, and ns can be varied to simulate a continuous
plume, frequent sampling is often less costly than frequent puff
releases. For one puff release per time step (np=l), and with sampling
frequency dynamically determined (see Section B.6), MES0PU1F execution
time (on an IBM 370/158) over the Four Corners computational grid and
typical meteorological conditions used is given by:
t (sec) = 0.38 (S) (N) + 0.43 (S) (5-3)
For runs with more than a few sources (so that the second term on the
right hand side of (5-3) is negligible), the MESOPUFF model is about
twice as expensive to run as MESOPLUME.
5.2.4 MESOGRID Run Time
The run time required for a MESOGRID simulation is determined by
the number of time steps, the number of grid points, and various mete-
orological factors that control the time step At satisfying the hori-
zontal and vertical computational stability criteria (see Section C.4.4).
Note especially that MESOGRID execution time vari'.s directly with the
number of grid points but is independent of the number of sources.
MESOGRID becomes cost-effective compared to the plume and puff models
only when many sources are to be included. The execution time (on an
IBM 370/158) of MESOGRID for the Four Corners computational grid (and
typical meteorological conditions) is described by:
w , 300 N* „
t(sec) = —— (5-4)
where N* is the total number of time steps (of duration At minutes) in
a MESOGRID simulation. Table 5-1 lists the maximum values of u and Kz
consistent with various choices of At, assuming a horizontal grid
interval of 40 km and vertical layer thicknesses of 500, 1,000 and
2,500 m. The values of Kz shown are the averages of the values of Kz
at the adjacent integer levels, e.g., Kz (1.5) « 0.5 (Kz(l) + Kz(2)).
122
-------
TABLE 5-1
MAXIMUM VALUES OF u AND VS. TM STEP fit
Maximum Kz(m2/s) Allowed by
the Vertical Diffusive
Advective Stability Stability Criteria
Criterion u
At (Bin) (m/s) Level l.S Level 2.5
5 - 625 2,916
15 - 208 972
30 22.2 104 486
§0 11.1 52 293
123
-------
5.3 Model Sensitivity to Changes in Time Step
As seen in the preceding section, for a specified number of
sources, grid si:e, and meteorological conditions, model execution time
varies approximately inversely with length of the time step used. It
may be expected, however, that as time steps are lengthened to save
execution time, a model will no longer reproduce the results obtained
when a very small time step was used. This section examines the repro-
ducibility of model results as the time step is changed. For model runs
made with different time steps, the 24-hour average S02 concentration
fields are compared.
Figures 5-1 a and b illustrate in highly condensed form the results
of numerous cost/sensit iv: ty tradeoff experiments made with MESOPUFF.
The figures display the execution cost (in IBM 570/158 CPU minutes) and
either the correlation coefficient (rI or the absolute value of the
fractional difference of the mean (|FDM|1, as puff sampling rate and
puff emission rate are relaxed from "special" base case values of
20 minutes each. Here, the correlation coefficient r (given as the
average of the three correlation coefficients computed individually for
each of the last 3 days of the simulation) is used as a measure of the
extent to which a "test" case would reproduce the results of the
"special" base case. The execution time of the base case was about
6 minutes for the 5-dav Episode 1 run with the Navajo power plant. In
contrast, a test run made with sampling and puff emission rates of
2 hours took only 43 seconds [an 8-fold savings in run time); yet,
qualitatively at least, the test run reproduced the base case results
very closely (r = 0.88, !FDM, = 0.06). This suggests that computa-
tionally "optimal" values of puff sampling rate and puff emission rates
can often be chosen (at least for this episode).
In Figure 5-2, the MESOPLUME and MESOPUFF models are compared for
departure from their respective base case results when only the plume
segment (or puff) emission rate is varied; the sampling rate is dynami-
cally determined and therefore nearly constant in each case. The
correlation coefficient r is again used to measure reproducibility of
model results (expressed as 24-hour average SO; ground level concen-
trations). The figure indicates that as the plume segment or puff
emission rate is increased (from a base case rate of 1 hour), the
MESOPUFF model, for the most part, reproduces its base case results with
more precision than MES0PLUME. Indeed, the MESOPUFF results for a
4-hour time step show a correlation of r = 0.9 with the results for the
base case; by contrast, for a 4-hour time step, MESOPLUME shows a corre-
lation of only about 0.6. Model execution times are, of course, sig-
nificantly reduced when using a longer time step--for the 4-hour time
step simulations, MESOPUFF required 40 seconds (compared to 97 seconds
for its base case), and MESOPLUME required 36 seconds (compared to
76 seconds for its base case).
124
-------
Puff Emission Interval (Hours)
figure 5-la MESOPIIFF CPU Run Tine (t) and Sensitivity of r to Variations
in Sampling Rate and Puff Emission Rate (Days 167-169)
125
-------
12
4
2
1
12
Puff Emission Interval (hours)
Figure 5-lb MESOPUFF CPU Run Time (t) and Sensitivity of FDM to Variations
in Sampling Rate and Puff Emission Rate (Days 167-169)
126
-------
t* 8f •« A rkV»>* V
Plums Segment or Puff Emission Rate (Hours)
Figure 5-2 Variation of Correlation Coefficient with Puff Emission Rate
127
-------
*t «***>•*•»
For the MESOGRID model, the only way to increase the time step
(thereby decreasing simulation cost) is to reduce the maximum values of
K. at any time step, so that no K_ values exceed the maximum allowable
value for a certain "requested" time step. Section 4,5 demonstrated
that maximum K- values could be reduced to be consistent with a 1-hour
requested time step without significantly affecting reproducibility of
results on regional scales. HESOGRID execution costs can be minimized
by requiring that the model time step be solely determined by the advec-
tive stability criterion (the Courant condition). In the first episode,
the maximum wind field velocities at each time step forced the model to
use a 30-minute time step during 80°. of the simulation and a 60-minute
time step during the remaining 20° of the simulation. Total execution
time was 18 minutes (IBM 370/158) for the 5-dav meteorological episode.
In all cases, the costs associated with archiving, combining,
analyzing, and producing hard copy plots of output results (using the
iMESOFILE postprocessing system) are negligible when compared to the
costs of running one of the dispersion models.
Figure 5-3 compares the execution time required by each model to
simulate Episode 1, days 165 to 169, 1978, as a function of the number
of sources included in a simulation. The MI.SOPLUME, MESOPUFF, and
MESOGRID curves are based on Equations 5-1, 5-3, and 5-4. Thus, for
this episode, MESOGRID becomes competitive with MESOPUFF when about 40
sources or more are included in a simulation, but about 80 sources must
be included before MESOGRID becomes competitive with MESOPLUME. [Note
that the MESOPLUME and MESOPUFF models as currently configured will
actuallv accommodate a maximum of 10 sources per run (to save on com-
puter storage requirements). For either model, the costs of performing,
for example, 4 simulation runs with 10 sources each and then merging the
results files using the MESOFIEE System are only negligibly greater than
the costs of actually running the models with a 40-point source inven-
tory; and restriction to 10 sources produces substantial savings in
needed core storage. MESOGRID is current 1y restricted to a maximum of
50 sources per run, but this restriction could be easily changed.]
The comparative run costs shown are based on time steps and other
model run parameters chosen to optimize the reproducibility of results
of base case runs made with very small time increments. The costs shown
are therefore somewhat higher than necessary; in particular, the time
step(s) used by the plume and puff models can be increased by a factor
of two or even four with minimal_loss of reproducibility. [Note that if
modeling of only SO2 (and not SOi,) is desired, it is possible to run
MESOGRID for SO7 alone and thereby reduce MESOGRID execution time by
almost a factor of two.]
Finally, the possibly substantial cost of generating the necessary
meteorological inputs to the dispersion models should be emphasized. In
practice, the cost of running MESOPAC per simulated day of hourly time
128
-------
f w«c'mmn-'K % *« c ««**¦>..>•>» v
Figure 5-3 Comparative Model Run Execution Times (5-day run)
129
-------
nc»Aint>k •vc+txoc* «*.
steps is about one-third that of MESOGRID, or comparable to a 20-source,
1-day simulation by one of the other models. If a lengthy set of mete-
orological data is to be generated by hand and then used for only a few
simulations by one or another of the dispersion models, the labor costs
associated with manually reducing and concatenating the raw rawinsonde
data will probably dominate the total cost (labor plus computer time) of
the entire simulation. In such an instance, using MESOGRID rather than
MESOPLUME or MESOPUFF would yield only a negligible increase in total
simulation cost.
130
-------
6. SUW1ARY, CONCLUSIONS, AND RECOMMENDATION'S
6.1 Basic Model Development Objectives
As shown by the many model comparison, model sensitivity, and model
cost analyses presented or summarized in the preceding sections, the
basic model development objectives have largely been realized. When
used and interpreted appropriately, each of the new models appears to be
a practical tool for simulating regional-scale air quality impacts of
energy-related major point sources of S02 in the Four Corners area.
These different models provide the user with a range of performance vs.
cost tradeoffs from which to choose; the choice of model best suited to
a specific application depends on (1) size of inventory, (2) charac-
teristic or dominant meteorological regimes, and (3) importance of
small-scale plume dispersion and plume fumigation to regional scale
impacts.
Each of these variable trajectory models takes into account the
important mesoscale variations in the meteorological fields governing
long-range plume dispersion in the Four Corners region. The meteoro-
logical input data are readily constructed from available rawinsonde
data using the MESOPAC preprocessor; use of alternative meteorological
preprocessing schemes (in place of MESOPAC) is easily achieved because
of the uniform meteorological interface structure built into each of the
three transport-diffusion models.
The contiguous plume segment approach adopted in MESOPLUME, and the
puff superposition approach adopted in MESOPUFF, have each been pre-
sented as generalizations of the conventional straight line Gaussian
plume model, useful for situations dominated by spatial and temporal
variations in wind fields, mixing heights, and atmospheric turbulence
levels. In the limit of uniform, steady flows, both MESOPLUME and
MESOPUFF can reproduce well the results of the straight-line model.
This feature of these models should aid in their general acceptance by
the air quality conmmity.
Other important features that contribute to the usefulness of these
models for regulatory and other applications include:
• the ease with which model design and run parameters and
algorithms can be changed to facilitate use for a range of
model validation research and development projects;
• the ease with which the models can be used to simulate short-
term average ambient ground-level concentrations, and the ease
with which the postprocessing system can analyze and display
model results; and
131
-------
• the complete freedom to change the geographic extent and
resolution of each mode 1 and to modify (replace) the meteoro-
logical preprocessing system.
6.2 Mode 1 Verification
lor lack of a suitable regional-scale data base in the Four Corners
region, the performance of these new models has not been verified to
date against actual ambient air quality in the region. However, the
basic computational method used in the MllSOC.RIP model has been pre-
viously verified to a limited extent on regional scales by MESOGRlD's
parent model SULFAS!) in the first phase of the llectric Power Research
Institute's Sulfate Regional Experiment f FPRI SURF) .
6.3 Comparative Model Performance
All three models compare extremely well on regional scales when
meteorological flow conditions are such that each of these modeling
approaches is appropriate. Fven when the nature of the meteorological
flow regime puts MESOPLUMF (and, to lesser degree, MF.SOGRID) at a
serious disadvantage relative to MLSOPUFF, the models remain in good
agreement at the lowest (most disperse) limits of the concentration
field calculations (0.1 ug/ffl") .
However, when stagnation conditions are encountered in a modeling
simulation, the inverse wind speed dependence of the MESOPLUMF. modeling
approach can lead to spurious results - as evidenced by several
instances when concentrations increased substantially by large downwind
distances. MLSOPUFF is not as st?r i ou; ly affected in this manner.
Similarly, under conditions of disorganized flow or where signifi-
cant local convergence occurs, the MESOGRID model also shows spurious
increases in concentration at large downwind distances, but apparently
to a lesser degree than MESOPLUML. This spurious behavior can arise
even when the two-dimensional wind field is nondivergent, because
parcels of material are advected at each time step by successive, inde-
pendent, one-dimensional operations in the east-west and north-south
directions, and one-dimensional wind fields may have appreciable local
divergence.
The significant differences that are apparent among results of the
different models in the near field of sources can be partially explained
in terms of the differences in the way the respective modeling algorithms
deal with (or ignore) the small-eddy diffusion effects. By contrast, at
large distances from the source, plume dispersion (as represented by
these models) is dominated by the advective effects of large-scale
edd ies in the mean wind field, and the differences in the way each model
treats small-scale diffusion are much less important.
132
-------
is-, rf ,-m^oick.v **:
On the larger regional scales, at least, the effects of numerical
pseudodiffusion in the MESOGRID model can be ignored--as can be judged
by the overall similarity of the grid model results (for small concen-
tration isopleth values) with the results of the plume and puff models,
neither of which is subject to pseudodi ffusive effects. Time-
discretization errors in all three models can be minimized by conser-
vative choices of advective (or diffusive) time steps. Space-
discretization errors in MESOPLUME can be serious in cases of highly
variable wind fields.
6.4 Model Sensitivity Results
In interpreting the detailed results of the model sensitivity
experiments shown in Sections 4 and 5, their principal limitations
should be borne in mind.
First, it will be recalled that the model sensitivity analyses were
performed as if the effects of varying one mode 1 parameter could be
considered as entirely independent of all the other model parameter
choices; in reality, many of the model parameters are strongly coupled,
either explicity or implicitly. For example, the advective transport,
diffusion removal, and chemical transformation processes are treated as
formally separable, although they are often strongly (if implicitly)
coupled--obviously the mass of SO; already converted to SO" is not
available to be deposited nor to be transported as SO-, for subsequent
transformation further downwind. Another but subtler case of parameter
coupling is the relationship between the diurnal variat ions in mixing
height (hence, in the onset of plume mixing to the ground), and the
resulting downwind pattern of surface deposition.
Second, the interpretation of the sensitivity results is compli-
cated by the implicit, but often strong, interaction among meteoro-
logical variables that govern large scale dispersion (e.g., mixing depth
and wind velocity) and model parameters that control the actual finite-
difference calculations (e.g., grid size or time step). When the values
of meteorological variables change rapidly in time and space, spatial
and temporal discretization errors in the calculation scheme may be
expected to cause significant errors, as, for example, in the calcu-
lation of individual plume element trajectories. Therefore, values used
for the time step or spatial interval should be sufficiently small to
permit good approximations to the "exact" plume trajectories in the case
of analytically-specified wind fields. In this study, computer core
size restrictions precluded the study of space discretization errors.
However, it has been shown for MESOPUFF and MESOPLUMP, (see Figures 5-1
and 5-2) that time discret izat ion errors do not appear to unduly
influence the calculated regional-scale 24-hour average S02 concentra-
tions, if a time step not larger than 2-4 hours is used.
133
-------
Next, it must be strongly emphasized that, for a number of reasons
concerned with model practicality, the physical parameterinations used
in these models are very simple, and do not therefore take into account
some of the more detailed (but often only fragmentary) information
available. Examples include:
• the use of constant, spatially uniform rates for in-situ plume
conversion of SO; to SO*--no provision has been made for the
sensitivity of the conversion rate to humidity, or to the rate
with which a plume entrains surrounding cleaner air;
• no provision for wet removal processes has been incorporated;
further, a complete and explicit model of wet removal would
require consideration of spatially and temporally variable
precipitation patterns and associated rainfall rates;
• the use of constant rates of dry deposition, with no provision
for factors influencing atmospheric and surface resistances to
SO2 deposition. Such factors may be spat ially variable
(e.g., soil type, vegetation type and percent coverage) or
diurnally variable (e.g., atmospheric stability, plant
stomatal resistance).
With the simple parameterizations used in this study the models are
shown to be very sensitive to rates of decay and dry deposition (see
Figures 4-1 and 4-2). A five-fold increase in either of these removal
mechanisms causes about a two-third decrease in Episode 1 24-hour
average SO; concentrations; in contrast, a five-fold decrease causes
about a one-third increase in concentrations. This behavior is highly
episode specific, however (see Figures 4-3 and 4-4). Obviously, even
when using very simple (e.g., linear, first order) removal rate pre-
scriptions, the model user must choose, with care, values for decay rate
and deposition velocity that are appropriate for the regional-scale
simulation at hand.
Lastly, it is emphasized that these test results, although pro-
viding much valuable information, are not definitive, in view of the
strong dependence of the results on the particular choices of meteoro-
logical and emission source scenarios.
With these limitations in mind, the detailed set of results pre-
sented in Sections 4 and 5 support the following useful general con-
clusions :
• Uncertainties in the meteorological data will usually exert
the dominant influence upon model sensitivity. The tests
provide clear evidence that the regional scale results are far
more sensitive to typical variations in the larger (meso-)
134
-------
which the surface winds are interpolated upwards o represent the layer
averaged wind field. This approach was not adopted in the MESOPAC
model, in part for lack of any interpolation method deemed adequate to
handle the strong topographic influences and convective conditions that
often dominate itesoscale flow fields in the Four Corners area. If these
new dispersion models are to be used in other regions, however, it might
be appropriate to modify MESOPAC to use supplementary surface wind data.
As another possibility, these models could be driven by a 1;pited-area,
Numerical Weather Prediction (NWP) model that has demons' ¦' prog-
nostic skill in the geographic region of interest over a .. ^d of, for
example, 12 to 18 hours.
6.6 Comparative Costs
By far the largest factor in the cost;- of doing regional-scale
simulation with one or another of these models is the size of the emis-
sions inventory. As shown in Section 5.3, only for large numbers of
sources, for example, 20 to 80 or more, does the grid model (tie
simulation costs of which are independent of source inventory size)
become cost-effective relative to the puff and plume models (the simu-
lation costs of which are proportional to the number of sources). Thus,
to model the incremental air quality impact of one or a few point
sources on existing ambient air quality in a region, either MESOPUFF or
MESOPLUME is appropriate (depending on the meteorology), and MESOGRID
would be extremely impractical. If, however, the emissions inventory
consists of, for example, 100 sources (and near-field plume resolution
is not critical) MESOGRID is likely to be the better choice (again,
depending to some extent on the meteorology).
The detailed nature of the meteorological scenario also plays an
important role in determining model costs. For MESOPUFF and MESOPLUME,
episodes dominated by stagnating or recirculating flows will cost more
to simulate than episodes of largely uniform flow conditions:
• For sufficiently stagnating or disorganized flows, MESOPLUME,
the model of choice for the more zonal flows may break down
and use of the somewhat more expensive (MESOPUFF) model is
mandated.
• For the plume and puff models, the execution costs per time
step are functions of average number of plume or puff volume
elements resident on the grid; under strongly recirculating
conditions, volume elements may remain on the grid for several
days. For MESOGRID, however, stagnation episodes may be less
expensive computationally whenever the Courant condition
permits the model to use a longer time step because wind
speeds are low.
Because the costs of running these models are so dependent on the
specific episode and experience with the models is as yet limited to the
136
-------
ew»o»««KTAL Tfo«
-------
! •• . ^ ft *F . f> • »s
TABLE 6-1
RECOMMENDED MODEL CHOICES
Criterion Preferred Model(s)
Number of Point Sources
• Small (X < 20)
• 'tedium (20 - S ' 80)
• Large (X > 80)
Dominant Meteorology
• Zonal flows, moderate to
high wind speeds
• Stagnating or recircu-
lating flows
• Variable flow regimes
¦ MESOPUFF
• Choice is scenario-dependent
• MESOGRID
• MESOPUFF or MFSOPLUME
(MESOGRID*)
• MESOPUFF (MESOGRID*)
• MESOPUFF (MESOGRID*)
Model Sensitivity
• Near-source concentration • MESOPUFF
estimates, including sen-
sitivity to plume diurnal
fumigation cycle, small
scale diffusion, and plume
rise.
*For large inventories.
138
-------
pkocwmiw* «wk>i i cx«iooi *
• regulatory applications, e.g., new source review under the PSD
provisions of the 197? Clean Air Act where (1) the chief
concern is the incremental impact on regional air quality of
one or several increment-consuming new sources (use of
MESOPUFF or MESOPLUME) or (2) where the regional background
concentration fields aust also be simulated, starting from a
sizable emissions inventory containing both point sources and
area sources (use of MESOGRID).
139
-------
s»lk -I.^' V
REFERENCES
Benkley, C. K. and A. Bass, 1979a. Development of Mesoscale Air
Quality Simulation Models. Volume 2. User's Guide to MESOPLtlME
(Mesoscale riumc Segment) Model. EPA 600/7-79-XXX, Environmental
Protection Agency, Research Triangle Park, N'C, 141 pp.
Benkley, C. h. and A, Bass. 1979b. Development of Mesoscale Air
Quality Simu1 at ion Models. Volume 5. User's Guide to MESOPUEF
(Mesocale Puff) Model. EPA W)0/""-"9-XXX. Environmental Protection
Agency, Research Triangle Park, NC, 118 pp.
Benkley, C. I\. and A. Bass. 19~9c. Development of Mesoscale Air
Quality Simulation Models. Volume b. User's Guide to the MESOPAC
(Mesoscale Meteorology) Package. EPA 600/7-79-XXX. Environmental
Protection Agency, Research Triangle Park, NC» 85 pp.
Benkley, C. W. and L. L. Schulman. 19">.». Estimating Hourly Mixing
Depths from Historical Meteorological Data. . Appl. Meteor.
18:~~2-~S0.
Blackadar, A. K. 1962. The Vertical Distribution of Kind and Turbulent
Exchange in the Neutral Atmosphere. J. Geophys. Res. 67: 3095-
3102.
Boris, J. P., and D. L. Book. 1973. Flux-corrected Transport--I.
SHASTA, A Fluid Transport Algorithm That Works. Jour. Comp. Phys.
11:38-69.
Briggs, G. S. 1975. Plume Rise Predictions. Lectures on Air Pollution
and Environmental Impact Analyses. American Meteorological Society,
Boston, MA, p. 67-114.
Clark, T. L. and R. E. Eskridge. 197". Non-Divergent Wind Analysis
Algorithm for the St. Louis RAPS Network. EPA 600/4-77-049.
Environmental Protection Agency, Research Triangle Park, NC, 63 pp.
Draxler, R. R. 1977. A Mesoscale Transport and Diffusion Model. N0AA
Tech. Memo ERL ARL-74, 31 pp.
Egan, B. A. and J. R. Mahoney. 1972. Numerical Modeling of Advection
and Diffusion of Urban Area Source Pollutants. J. Appl. Meteor.
11:312-322.
Egan, B. A., k. S. Rao, and A. Bass. 1976. A Three-Dimensional
Advective-Diffusive Model for Long Range Sulfate Transport and
Transformation. Seventh International Technical Meeting on Air
Pollution Modeling and its Application, Airlie, VA, September 7-10,
1976, p. 697-714.
140
-------
I Nv **ONMT»«tAt ART* 4 'f r «*>. 00» V
REFERENCES (Continued)
Colder, P, 1972, Relations Among Stability Parameters in the Surface
Layer. Bound. Layer Meteor. 3:47-58.
Hales, J. M.» D. C. Powell and T. D. Fox. 1977. STRAM-An Air Pollution
Model Incorporating Non-linear Chemistry, Variable Trajectories,
and Plume Segment Diffusion. EPA 450/3-77-012. Environmental
Protection Agency, Research Triangle Park, NC, 147 pp.
Heffter, J. L. 1965. The Variations of Horizontal Diffusion Parameters
with Tine for Travel Periods of One Hour or Longer. J. Appl.
Meteor. 4:153-156.
Heffter, J. L. and A. D. Taylor. 1975. A Trajectory Model, Part I. A
Regional-Continental Scale Transport, Diffusion and Deposition
Model. NOAA Tech Memo. ERL ARL-50, 28 pp.
Hilst, G. R. 1970. Sensitivities of Air Quality Prediction to Input
Errors and Uncertainties. Proceedings. Symposium on Multiple-
Source Urban Diffusion Models. APC No. AP-86. Environmental
Protection Agency, Research Triangle Park, NC, pp. 8-1 to 8-40.
1-tusar, R. B. , J. P. Lodge, and D. J. Moore, eds. 1978. Workshop:
Sulfur in the Atmosphere. Atmos. Environ. 12:7-23.
Liu, C. Y. and W. R. Goodin. 1976. An Iterative Algorithm for Objec-
tive Wind Field Analysis. Mon. Wea. Rev. 104:784-792.
Liu, M-K. and D. D. Durran. 1977. The Development of a Regional
Air Pollution Model and Its Application to the Northern Great
Plains. EPA-908/1 - 77-001. Environmental Protection Agency,
Denver, CO, 281 pp.
Morris, C. S., C. W. Benkley and A. Bass. 1979. Development of
Mesoscale Air Qulity Simulation Models. Volume 4. User's Guide
to MESOGRID (Mesoscale Grid) Model. EPA 600/7-79-XXX Environmental
Protection Agency, Research Triangle Park, NC, 105 pp.
Nappo, C. S. 1978. Workshop on Long-Range Trajectory-Puff and Plume
Modeling of Cont inuous Source Emissions. NOAA Tech. Memo. ERL-
ARL-72, 18 pp.
Nuber, J. A., A. Bass, M. T. Mills and C. S. Morris 1977. A Review of
Regional-Scale Air Quallty Models for Long Distance Dispersion
Modeling in the Four Corners Region. EPA-600/7-78-066, Environ-
mental Protection Agency, Research Triangle Park, NC, 82 pp.
0' Brien, J. J., 1970. A Note on the Vertical Structure of the Eddy
Exchange Coefficient in the Planetary Boundary Layer. J. Atmos.
Sci. 27:1213-1215.
141
-------
I'MVMBBTMMFV1'**. i t -'X- *4.
REFERENCES (Continued)
Pack, IK M., C. J. Eerber, J. L. Heffter, k. Telegadas, J. K. Angel 1 ,
K. H. Hoecker and L. Machta. 1978. Meteorology of Long-Range
Transport. Atmos. Environ. 12:425-444.
Plate, E. J. 1971. Aerodynamic Characteristics of Atmospheric Boundary
Layers. U.S. Atomic Energy Commission [MIS Ref. No. TID-25465],
190 pp.
Scire, J. S., J. E. Beebe, C. K. Bcnkley and A. Bass. 1979. Develop-
ment of Mesoscale Air Quality Simulation Models. Volume 5.
User's Guide to the MESOFI EE Postprocessing Package. EPA 600/7-
79-XXX. Environmental Protection Agency, Research Triangle Park,
NC, 67 pp.
Sheih, C. M., G. P. Hess and B. B. Hicks. 1978. Design of Network
Experiments for Regional-Scale Atmospheric Pollutant Transport
and Transformation. Atmos. Environ. 12:1745-1753.
Start, G. E. and L. L. Wendell. 1974. Regional Effluent Dispersion
Calculations Considering Spatial and Temporal Meteorological Cal-
culations. NOAA Tech. Memo. ERL-ARL-44, 63 pp.
Turner, 1). B. 1970. Workbook of Atraosoheric Dispersion Estimates. U.S.
Dept. of H.E.W, Public Health Service, Publ. 999-AP-26, 88 pp.
Wilson, K. E. 1978. Sulfates in the Atmosphere: a Progress Report on
Project MISTT. Atmos. Environ. 12:537-548.
142
-------
t -,\«r ~.'»> &* sf *>« M* -K'.v .v
APPENDIX A
MESOPLUME TECHNICAL DISCUSS IOC
A.1 Basic Mass Conservation Equations
MESOPLUME, a variable-trajectory version of the conventional
straight-1ine Gaussian plume Model, is designed to take into account the
spatial and temporal variations in the horizontal advection, diffusion,
transformation, and removal mechanisms governing the dispersion of a
pluoae on regional transport scales. In MESOPLUME, a continuous plume is
modeled by subdividing the plume into a number of contiguous 'plume
segment' elements, as illustrated in Figure A-1. The conservation of
pollutant mass over a plume segment of length is is expressed by the
¦ass balance equation:
m
j / G(s,r,z) dr dz
is
(A-l)
/ / u C dr dz
S* As
- / / u C dr dz
where s, r, and z are the longitudinal, lateral, and vertical plume
directions, G(s,r,z) (g m~3s~1) is the rate of change (gain-loss) of
pollutant concentration C(s , r, z) (g m 3) bv conversion and removal
processes, AQ (g s J) is the resultant rate of change of pollutant mass
and u (m s 1) is the wind speed. In the MESOPLUME model, G(s,r,z) and u
are considered to be constant from s to s + As, where s is the current
distance of a plume segment endpoint from the emitting source, measured
along the plume axis.
MESOPLUME permits the user to specify two possible vertical
distribution functions; (1) a vertical Gaussian profile, ignoring any
effects of the mixing lid H; or (2) a uniform vertical distribution
below the mixing lid.
For Case 1, the ground-level axial plume concentration C(s,r,0) is
defined at the upwind edge of a plume segment by the expression:
C(s,r,0)
21*2
it u cjs) 0v(s)
exp
( 1
-r~
exp
f
j
*>
2a 1
2o/
y
(A- 2)
Preceding pageu blank
-------
146
-------
where Q(s) is the pollutant mass flux and o,ls) , o.(s) are the lateral
and vertical Gaussian plume dispersion coefficient distributions at
downwind distance s. Full reflection from the ground is assumed.
For Case 2, if the plume altitude (z) lies below the mixed lid H,
the ground-level axial concentration is expressed at the upwind edge of
the plume segwent by the expression for uniform vertical mixing:
("(s,r,0)
Of s>
/Jn u II 0 (s)
m v
exp
-r
(A-3)
where 11^ is the maximum mixing depth encountered by the plume segment
(see Section A.8). If, rather, the plume centerline lies above the
nixing lid, no ground-level concentrations are calculated. At the
downwind edge (s+As) of the plume segment, the ground-level axial
concentration (s+6s,r,0) is expressed as:
C(s+ds,r,0)
(Q(s)+(dQ/dt)}At
/T? u H a (s*&s)
my
exp
-r
. 2
2 o
v
C A - 4 l
The MESOPLUME model solves the mass conservation Equation A-1
independently for both sulfur dioxide (S02) and sulfate (SO,,J. The gain
functions include terms for the loss (gain) of SOz (S0U) by linear decay
of SO • to SOl] and terms for dry deposition of either species (see
Sections A.5 and A.6).
A.2 The Grid System
The original coordinate system used in the STRAM model has been
replaced in MESOPLUME by a simple Cartesian coordinate system. All
spatial model input data (emission source inventories and meteorological
fields) are referenced to the same grid, called the "basic computa-
tional" grid. To improve the resolution of the plume sampling function
(see Section A.4.3), MESOPLUME uses a sampling grid that is a subset of
the basic computational grid. The origin of the sampling grid may be
placed anywhere on the basic computational grid (but not on the northern
or eastern grid boundaries). The resolution of the sampling grid is a
¦ultiple of the resolution of the basic computational grid. Currently,
the maximum allowable size of both the basic computational and sampling
grids is 40 by 40 horizontal grid points. Figure A-2 illustrates one
possible arrangement of basic computat ional and sampling grids. Note
that the grid index of the origin (1,1) is assigned Cartesian coor-
dinates x = 0, y = 0. The sampling grid in this example originates at
(i = 2, j = 2) or (x = 1, y = 1) on the basic computational grid, and
has a spatial resolution four times finer than the basic computational
grid.
147
-------
y-5.0
j-«
y-4.0
,-5
v « 3,0
j-4
y-2.0
j = 3
¥=1.0
' i = 2 "
Ay = &d
y = 0.0
" i-1
* = 0.0
i = 1
I I
- + + + -
- + + + -
-+ + + -
I i i
l I F
I- H 1
_ _| | |
i i i
1 I I
1 1 I _
"h++-
" + ++-
i 1 1
i t t
- + + + -
- + + + -
I 1 i
I i i
¦w ann^HiB M^iaB mm
-+++"
" + + +"
1 1 1
1 I !
- + + +1
-+++-
—1—h + "
1 1 1
l i i
" |- ""I" ||11"^""1 •—
1 1 I
1 I 1
-hH- + "
1 1 1
1 1 1
__J—1__|—
-+ + + -
i i i
x -1.0
i-2
* = 2.0
t = 3
* = 3.0
i = 4
x = 4.0
x = 5.0
i-6
[~i* = 4d-*j
Figure A-2 One Possible Arrangement of Sampling and Basic Computational
Grids
I
f.
148
-------
twwHWMOmnaiHMKXt TKWNOl00» »C
.4.3 Specification of Model Inputs
MESOPLUME is, at present, driven by the meteorological fields
produced by MESOPAC, a MESOscale meteorology PACkage (Benkley and Bass
1979c), but an alternative meteorological preprocessing system with
appropriate grid resolution could be substituted by the user. MESOPLUME
requires the following fields (usually at hourly intervals) interpolated
to each model grid point:
• horizontal (u, v) wind components,
• nixing depth, and
• Pasqui11-Gifford-Turner (PGT) stability class.
Use of MESOPAC model output by the MESOPLUME model is straightforward -
the Cartesian coordinate system is identical in both models. Note that
the horizontal (u,v) wind components are considered invariant with
height. Vertical wind components are not used. The fields of PGT
stability class are used by MESOPLUME in doing plume rise and plume
growth calculations.
A.4 Plume Trajectory, Dispersion and Sampling Algorithms
The computational scheme of the MESOPLUME model has three distinct
functional elements: (1) a Lagrangian plume trajectory function, (2) a
plume dispersion function, and (3) a plume sampling function. The
Lagrangian trajectory function is used to advect the endpoints of each
plume segment during a basic time step; the resultant distance between
consecutive endpoints defines the length of each plume segment. The
widths at the upwind and downwind ends of each plume segment are deter-
mined by the plume dispersion function. Given the size and location of
each plume segment, the plume sampling function computes the concen-
tration exposure received during the time interval at each grid point
that lies within the plume segment.
2.4.1 Lagrangian Trajectory Function
This section describes how the endpoints of a plume segment are
advected during a time step; it has been adapted (but largely verbatim)
from Hales et al. (1977, pages 15-18).
Let the components of the wind in the x and y grid directions along
the trajectory be designated by u[t;x(t),y(t)] and v[t;x(t),y(t)] ,
respectively,
where
, t
x(t) = I u[t';x(t'),y(t')]dt'
149
-------
) 1
y(t) = -jj / v[t' ;x(t' ) ,y (t' ) ldt'
O
x(0) = x0 ; y(0) = yp
u(0;x(0),y(0)) = u(x0,y0)
v(0;x(0), y (0)) = v(x0,y0)
and x0,y0 are the initial coordinates (origin) of a trajectory of total
length s*. (The grid space unit Ad is used in the program to convert to
nondinensional spatial units on the computational grid.) The vectors
u(t;...) and v(t; ) provide a Lagrangian description of the velocity
field. Increments of advection over a further time interval At are
therefore given by
. t+fit
Ax = ig- / u [ t' ;x(t' ) ,y (t') ] dt' ; (A-5a)
. t+At
&y = _ J v [t' ;x(t') ,y (t') ] dt* . (A-5b)
After the new coordinates x(t+At) and y(t*At) are calculated for each
incremented endpoint of a given plume, a check is made to see if any of
these values are off the grid. If so, both the affected plume increment
and all plume increments emitted from the same source previously are
deleted from further consideration.
The integration of the wind components over the trajectory is
approximated by a two-step iteration method involving bilinear inter-
polation in space and in time (see Figure A-3).
The various positions in Figure A-3 are best defined by equations:
Xj = x(t) + u[t;x(t),y(t)] At ; (A-6a)
yx = yft) ~ v[t;x(t),y(t1] At ; (A-bb)
x, = Xj + uft+At;Xj,y ] At ; (A-6c)
y, = yj ~ vffAtjXj, >j ] At ; (A-6d)
*ln practice, however, the exact definition of s in MESOPLUME is the
current along-the-variable-plume-axis distance of a segment endpoint
from a source. In time-varying flows, it may be somewhat different
from the actual segment endpoint trajectory length.
ISO
-------
F*fv«M3NMifwTfti. rtewNOiOGv wt
Figure A-3 Calculation of the Trajectory of a Plume Segment Endpoint
151
-------
xffSt) = 0.S|x(t) ~ xj ; (A-6e)
yft+fit) = 0.5(y(t} * v,| . (A-6f)
Thus, the new position (xft+At),yis found by adding the
vector average of two advect ion increments to the former position
(x(t),y(t)). The first advection increment is calculated by advecting
the plume increment for the entire interval At using the wind effective
at (x (11, y (t) ] at time t. The addition of this increment to the posi-
tion Ix(t),y(t)] yields the position (xj,yj). However, (xj,yi) is not
assumed to be the actual end of the advection increment because the wind
may change along the trajectory. A second advection increment is cal-
culated using (xj,y}) as the starting point and the wind at that point
effective at the end of the ti«e increment. The addition of this advec-
tion incresent to position (xj.yj) yields position (xa.ya). Then the
new pluae increment position (x(t*At),y(t*At)) is taken to be the point
halfway along the line from (x(t),y(t)) to (x2,y2)•
The bilinear interpolation by which the effective wind components
u(t), and v(t) are calculated works as follows. Let tn and tn+j be the
effective times of the two gridded wind fields closest to time t. Time
interpolation weights tj and t:- are defined by:
t - t
n+1
- t
t < t < t ,
n — n+1
(A-7a)
tj = 1 - t7. (A-7b)
Next, the location of (x(t),y(t)) is noted on Figure A-3. The coor-
dinate values are between 3 and 4 for x(t) and between 2 and 3 for y(tJ.
Accordingly:
X = x (t) - 3 ;
V
y(t) - 2 ;
1 - Y .
(A-8a)
(A-8b)
(A-8c)
(A-8d)
152
-------
mviticmMWM. »TfeMNatco* #c
This yields
"(t) = S *p V uttn;3,2) ~ t, Xp Yp U(tn+1;3,2)
* lI *q Vp ul;4'2) (A"9)
* ll Xp * X2 Xp Yqu(tn+1:3»3)
+ *1 \ \ * h \ \ UfW4'3) •
Similar equations hold for v(t)» u(t+At), and v(t*At). Therefore, the
advection Equation A.5 is calculated from
Ax = tuft) ~ u(t*At)}At; (A-lOa)
Ay = (v(t) ~ v (t ~ At) }At. (A-lOb)
The subsequent positions are given by
x(t +At) = x(t) + Ax ; (A-lla)
v(t+At) = y(t) + Ay . fA-llb)
A.4.2 The Plume Dispersion Function
This section, describing the evaluation of the plume dispersion
parameters, has been taken largely verbatim, with modifications, from
>tales et al. (1977, pp. 23-24). The plume dispersion parameters Oy and
az are calculated for distances out to 100 km using plume growth
formulas fitted to the curves of Turner (1970). For distances greater
than 100 km the plume growth rates given by Heffter (1965) are used.
The growth of a plume with travel time t or along the plume trajectory
distance s from the source is represented by
do
oy(s*As) = o (s) + As I (A-12}
I s-»As/2 ;
153
-------
. (t ~ At)
+ '.t
dt
t+M/2
(A-13}
Similar equations for oz are used. These terms allow for spatial and
temporal changes in stability class to be included, without violating
the entropy principle (centerline concentrations cannot increase with
downwind distance).
The integral formulas for ov and o, for travel distances less than
100 km are of the following form's
0 q
o ,(s,a) = Y s (meters); (A-14)
V* Cl
b
c^(s,a) = s a (meters). (A-15)
Here
is a stability index designated as A, B, C, I), E, or F for the
PGT stability categories. The coefficients Y#
and bri are given in
¦a. .
Table A-1 as a function of stability index a and yield cv and a~ in
meters when s is specified in meters. Because the integral formulas are
not valid if the stability class changes over the travel distance s, the
derivative forms are actually used to carry out the computations
ds
= 0.9 Y
-0.1
(A-16)
do_
ds"
b A s
a a
(b -1)
a
(A-17)
MESOPLUME assumes that ov, o,
0 at the origin (s = 0).
The derivative formulas used to grow o.. and c- for travel distances
greater than 100 km are those given by Heffter (1965):
do
dt
y , -1,
Cm s )
0.5
(A-18)
do_
dt
-i i /"> -i n
(m s ) = 0.5 (2KJ t '
(A-19)
154
-------
TABLE A-l
COEFFICIENTS FOR DISPERSION PARAMETER FORMULAS
Pliaae Growth Coefficients
Stability
Index a Y Z b
a a a
A 0.36 0.00023 2.10
B 0.25 0.058 1.09
C 0.19 0.11 0.91
D 0.13 0.57 0.58
E 0.096 0.85 0.47
F 0.063 0.77 0.42
155
-------
,4,4,3 The Plume Sampling Function
Each plume segment resident on the grid (at the end of every time
step At ) is sampled, using the plume sampling function, to evaluate the
average concentration experienced at each sampling and intersection
during the previous ti«e step. For example, consider the hypothetical
plume segment depicted in Figure A-4. The plume segment centerline at
time t ~ At extends from {x,y) to (x+Ax,y+Ay) , the positions of the two
consecutive plume sofwent endpoints at t »'t. The plume segment length
is is = (Ax2 + Ay2 j1/'. The lateral extent of a plume segment is con-
sidered to be truncated at t.yiv. This is a reasonable simplification--
much less than 1» of the area under the Gaussian distribution function
lies beyond ±5cy from its center. In this example, plume segment radii
of size Sov(s) at the upwind edge (x,y) of the segment and of size
3ov(s+As) at the downwind edge (x+Ax,y+Av) are indicated. At time t+fit,
the grid point intersections (21,12), (22,12), (22,13), (23,13), and
(22,14) are each impacted by the hypothetical plume segment; each grid
point is assigned a certain average concentration C(i , j) resulting from
the presence of the plume segment over At. This evaluation is illus-
trated as follows: Suppose the grid point concentration C(i,j) is to be
calculated at (22,12). First, the ground-level concentrations C(x,y)
and C(x+Ax, y*ty) are computed using F,quations 2-3 and 2-4 for the case
of uniform vertical distribution. Next, a point (x' ,y' 1 is found such
that the 1ine segment of length r const rueted from (x' ,v') to the grid
point (i, j) is perpendicular to the plume segment centerline. The
effective source strength C(x',y') is then computed by linear inter-
polation :
C(x',y')
C(x,y) As^ + C(x+Ax, y+Av) As j
CA-20)
ASj + As,
The lateral plume dispersion coefficient o (x.y) is similarly inter-
polated from: '
oy(x',>'•)
av_(x,y) As, + a^(x+Ax, y+Av) As^
CA-21)
As. + As,
Finally, the grid point concentration C(22,12) is computed as:
C(22,12) = C(x',y') exp
2ov2(x',y«)
(A-221
156
-------
EV.«ONMf*'V *fCHNCXOr.* *iC
Figure A-4 Calculation of the Concentration C(i,j) by the Sampling Function
I
157
-------
If a grid point lies within nore than one plume segment at the end of a
tine step (for example, because plumes of multiple sources overlap), the
total concentration at the grid point C^.(i,j) is computed as the sum of
the individual contributions:
C-T(i,j)
mT
T
¦=1
(A-23)
where m-j- is the total number of segments impacting the grid intersection
{i,jl during the tine step of interest. Note that a plume segment is
only sampled beginning with the second time step of its history and at
every tine step thereafter.
The MESOPLUME model has two peculiarities intrinsic to the basic
advect ive-d i f fus ive scheme of the model, which can present problems when
the sampling function is applied. These problems arise from (1) the
inverse-wind speed dependence of the concentration algorithm (Equation A-2)
and (2) the way in which adjacent plume segments are actually juxtaposed
in a curvilinear flow. As an example of the inverse-wind speed depen-
dence problem, consider the case illustrated in Figure A-5, for which
the wind flow along the plume axis rapidly decelerates, and suppose ua,
the average wind speed over segment A, is twice uB, the average wind
speed over segment B.* Neglecting removal processes, MESOPLUME evaluates
the concentration C(xj,yi) at the downwind edge of segment A according
to Equation A-4:
C(Xj.Vj)
,/Fua Cv(V>V H,
(A-24 1
but can also evaluate the concentration C(x! ,y i)L at the same point
us ing Equation A-3 as appropriate to the upwind edge
Ug instead. Since u^ = 2Ug,
C(xl .yp L = 2C(x],y1)
ge of segment B using
(A-25)
If the grid resolution is such that severe discontinuities in the wind
field occur, the MESOPLUME model will indeed predict plume centerline
concentrations increasing with downwind distance. This prediction, a
violation of the second law of thermodynamics, represents a major
shortcoming of the plume segment model--it cannot treat stagnation flows
correctly.
*Th is does not imply that the two-dimensional wind field is significantly
divergent--the local flow about segments A and B may be exactly non-
divergent if the mass fluxes through the lateral boundaries are "correct."
158
-------
To address the second problem, that of juxtaposing adjacent plume
segments, consider the situation illustrated in Figure A-6. In this
figure two adjacent plume segments are advected in a highly curvilinear
flow. In this case, for each plume segment, the lateral lines that
define the plume segment widths (at the upwind and downwind end points)
arc constructed to be perpendicular to the local plume segment axes.
Therefore, if the average wind direction used to advect segment B is
very different from the average wind direction used to advect segment A,
the plume segments will not be perfectly contiguous, as shown. In the
example, there is a shaded region in which grid points will receive two
concentration doses during one time step and also a region where the
plume may actually "pass over" a grid point without impacting it at all.
"litis illustrates a second potentially serious shortcoming of the model -
it cannot treat strongly sheared flows (e.g., recirculating flows)
correctly.
A.5 Conversion of Sulfur Dioxide to Sulfate
MESOPLUME represents the conversion of S02 to SO^ as a simple
linear rate function--the changes in mass of each pollutant in one time
step St are given by the time-discretized equations:
AQn(SO,) = kj Qn(SO,) At (A-26a)
AQn(SO*) = -1.5 kj Qn(SO,) At (A-26b)
where Qn(S02) is the mass of S0? at the beginning of the nth time step
and ki(s 1) is the conversion rate constant. MESOPLUME uses the nominal
value kj = -5.56 x 10~6, suggested by Hales et al. (1977), unless
otherwise specified by the user.
Note that the conversion of S02 to SOi^ represented by Equation A-26
is (1) independent of the vertical or horizontal distribution of S02
within a plume and (2) independent of whether a plume lies above or
below the mixing height.
For the solutions to be adequate representations of the continuous
exponential equations (e.g., Q(S02) = kiQ(S02)), Q should not change
appreciably over the time interval At. Since this is frequently not the
case, the time interval At is automatically subdivided into as many
subintervals At* as are required to ensure that no more than a specified
fraction AQ^ of the mass Q is removed by decay during a subinterval At'.
This process is described further in the next section.
159
-------
;v,-«c•¦wv**, • • • • •-
(K»Y.I
Figure A-S
(*.,V.)
A Schematic Representation of Adjacent Plume Segments In
Rapidly Deceleratinc Flow
Figure A-6
A schematic Keprcsentation of Adjacent Plume Segments in
Highly t.urv 11 incur Flow
160
-------
A.6 Dry Deposition of Sulfur Dioxide and Sulfate
The changes of mass AQn{gJ of each pollutant resulting from dry
deposition are given by:
AQn(SO,) = - (v j (SO,) Qn(SO,J At! / Hb
(A-27a)
AQn{SO*) = -(vd(SO*) Qn(SO^) At) /
(A-27b)
where Qn(SO?) and Qn(SO^J are the masses of S02 and SO^, respectively,
in the plume segment at the beginning of the time step, v^fSOj) andk
vd(SOr() are the deposition velocities of each pollutant, and "n, is the
vertical depth of the plume segment (see Section A.8). MESOPLUME uses
the nominal values vjfSOj) = 0.01 (m/s) and vd(SO,J = 0.001 (n/s), as
suggested by Hales et al. (197"*'), unless otherwise specified by the
user.
As described before, MESOPLUME automatically ensures that the mass
Q does not change within one sub interval At* by more than a user-
specified fraction AQf. Taking into account both SO? removal mechanisms
(decay and deposition), the expression used by MESOPLUME to specify the
subinterval at' is
In other words, the subinterval size At' is chosen automatically by
MESOPLUME at the beginning of each basic interval At such that no more
than a specified fraction AQf of SO- mass is removed by either mechanism
(or, no more than 2AQf by both mechanisms combined) within At'. It
should be emphasized that the conversion and removal rates are not
affected by this restriction. Also, individual puff trajectories are
not updated at each subinterval At'; they remain unchanged over the
basic interval At.
A.7 Plume Rise
The effective emission height he of each plume segment is computed
as hg = hg + hp, where hs is the stack height (m) and h„ is the plume
rise (ra). The plume rise equations used by the MESOPLUME model are
those described by Briggs (1975) for equilibrium (final) plume rise.
For unstable and neutral conditions with h' <_ H (i.e. , the plume does
not rise into an elevated stable layer)
/
'•Q f (SO ,)
At' = MIN ¦ r —
1
AQ (SO,) H
t m
CA-28)
hp = h' = 1.6 F1/3 (3.5 x*)2/3 u"1
(A-29a)
161
-------
For unstable and neutral conditions with h' > H (i.e., plume penetrates
into on elevated stable layer)
1 -1 -1 1/3
h « MIX |h' » c + 18,75 F u S ) } ; (A-29b)
p OB
for stable conditions when u _> 1.37 a/s
h • 2.6 F1/3 S"1/J u"1/3 ; (A-29c)
P
for stable conditions when u < 1.37 m/s
h =5.0 Fl/4 S-3/8 . (A-29d)
p
where:
F = buoyancy flux (m4/s3)
x* = 34.49 F0-4 for F > 55
= 14.0 F0-625 for F < 55
S = (g/T) (30/Bz)
g = 9.8 m s "
T = 2908K
36/3z = 0.0137°K a"1
H = mixing depth (is)
zfe = H - hs (.)
u = wind speed (« $ S
u = MAX {u , 1.37}
m
A.8 Treatment of Plume Fumigation
When the "uniform vertical distribution" option is elected,
MESOPLUME takes the spatial and temporal variations in the mixing height
into consideration to determine the extent at ground level pollutant
impact:
(a) For a plume element centerline less than the mixing height,
MESOPLUME arbitrarily assigns the element a height of zero
and immediately mixes the mass of the plume element uniformly
through the mixing depth, with the mixing lid acting as a
perfect reflector. Equations A-3 and A-4 are then applied
to compute around level concentrations.
162
-------
mMMm «w««e» «tto«ot,«w «c
(b) For a plustc el event centerlinc greater than the mixing height,
no impact of the element is felt at the ground. However, if
subsequently the nixing height becomes greater than the height
of the plume element centerline, the entire pluae element is
immediately nixed uniformly through the nixing depth, and
Equations A-3 and A-4 again apply.
The nixing depth encountered by a plume element is likely to change
over time and space. MESOPLUME assumes that a plume element residing in
the nixed layer is always uniformly mixed throughout the maximum mixing
depth H,,, encountered by the element in its progression through the
computational space-time grid.
The restriction of uniform pluae element properties through a
character!zable depth such as Hg is imperative for a one-layer model
such as MESOPLUME to ensure that parcels of material, once entrained,
are not bifurcated in the vertical when the height of the mixing lid
changes. Therefore, in the MESOPLUME model a plume segment, once
entrained, is uniformly mixed through a realistic height, consistent
with the following assumptions:
(a) When convection ceases in the late afternoon, turbulent energy
in the afternoon mixed layer dissipates. Once the energy of
the large convective cells is completely damped, the remaining
small scale energy is not sufficient to yield a vertical
entrainment rate comparable to the daytime entrainment rate
produced by large convective eddies.
(b) The stable lid capping the convective layer most likely
maintains some integrity after convection ceases and acts to
prohibit any additional vertical entrainment of plume material.
(c) A modeling approach that redistributes mass into a shallower
layer when the mixing depth decreases will violate the second
law of thermodynamics.
Since no mechanism remains for turbulent transfer of pollutant
material either upwards or downwards once the mixing depth decreases, it
is asserted that the height most appropriate to define the extent of
vertical mixing of a segment is the maximum mixing height that the
segment encounters during any portion of its travel.
A hypothetical example is given in Figure A-7 to schematically
illustrate the interaction of the MESOPLUME vertical distribution
algorithm and the mixing depth progression algorithm used by MESOPAC
[Benkley and Bass (1979c), Benkley and Schuiaan (1979a)J. Here consider
two instantaneous parcels of material released at different times from a
source with an effective stack height of 250 m; the first release is
made at 0300 GMT on the first day; the second release is made at
163
-------
rntms
&•
18 21 OZ 3
Time of Day
21 OZ
Figure A-7 Response of Two Plume Element s to Changes in Mixing Depth
-------
0300 GMT on the second day. Assume that the hypothetical source is
located in the western United States - so that 0000 GMT corresponds to
late afternoon -the usual tine of saximum mixing depth.
The first parcel is entrained at 0600 GMT on Ray 1 as the nocturnal
hour ary layer oscillates slightly upward in height. MESOPIUME immedi-
ately nixes the parcel uniformly throughout this mixing depth. This
first parcel is nixed uniformly through 900 m by 0000 GMT on Day 2, and
remains nixed through 900 m subsequent to that tine, even though the
height of the mixed layer collapses. The second parcel, emitted at
0300 OH on Day 2 is fumigated at 1500 GMT and is mixed uniformly
through TOO m by the end of the afternoon. Therefore, the first parcel,
en it ted at 0300 GMT on Day 1, remains uniformly mixed through a deeper
vertical layer than the second parcel; this is consistent with the
greater extent of nixing on the first day as compared to the second day.
This mixing depth/fumigation scheme exhibits at least two short-
comings.
• When a plume segment becomes bifurcated in the vertical by
collapse of the mixing lid, how should the plume segment
lateral dispersion rate be characterized? Because MESOPLUME
is a one-layer model, the same set of horizontal dispersion
coefficients must be applied to both the above-mixing-lid and
the below-mixing-lid portions of the bifurcated segment.
• The dry deposition algorithm (see Section A.6) removes S0? and
SO^ at a rate proportional to the surface concentration--this
loss is then redistributed throughout the entire vertical
extent of the plume. It would be more realistic to distribute
the loss only through the current mixing depth, but this is
not possible in a one-layer plume segment model.
It is felt that these limitations of the scheme are necessary for a
method that is both computationally practical and yet realistic for use
with a one-layer model. The inclusion of a fumigation cycle is very
important; the marked effects of plume fumigation on the distribution of
regional-scale ground-level concentration patterns are illustrated in
Section 3.3.
A.9 Comparison to the Conventional Gaussian Plume Model
From the point of view of possible regulatory applications, one of
the most attractive features of MESOPLUME is its ability to recreate the
results of the conventional Turner Workbook plume model (using PGT
coefficients). When MESOPI.UME is run with appropriate uniform steady-
state meteorology and fine spatial resolution, it can reproduce the
165
-------
Turner workbook results to a high degree. This can be demonstrated, for
exanple, under the following test conditions:
• The horizontal grid spacing of the model is reduced, e.g., to
5 km.
• The meteorological fields (wind direction, wind speed, mixing
height and POT stability class) used to drive MESOPLUME are
taken to be spatially uniform and constant in time.
• MESOPLUME is run with one source until quasi-steady-state
conditions are well established; that is, plume segments are
created at the same uniform rate at which they disappear off
the edge of the grid.
Under these special test conditions, MESOPLUME simulates very well
the Turner Workbook plume concentrations out to distances of 100 km from
an emission source. (The PGT cv and az curves, represented in the model
by piecewise linear power law fits, are only defined to 100 km.)
Table A-2 and Figure A-8 provide a representative illustration of how
closely MESOPLUME results compare to those obtained for the conventional
Gaussian plume model using the PGT "0" curve for Oy, and a uniform
vertical distribution. Overall, the comparison is exce1 lent --the very
small fractional differences (typically less than 2°) are attributable
mostly to inexact fitting of the PGT ov curve by the power law function
used in MESOPLUME. The somewhat larger deviations near the source are
expected - the results shown here were obtained with a wind velocity and
time step corresponding to plume segments about 5 km in length - so,
near the source, the linear averaging done by the sampling function does
not approximate accurately the power law behavior of the plume center-
line concentrations. It is expected that with use of smaller time steps
and correspondingly shorter plume segments, the comparison of MESOPLUME
to the turner Workbook plume in the very near field would be even
closer.
In its present version, MESOPLUME does not include an option for a
reflected Gaussian vertical distribution, so that at distances where the
plume does not yet approach a_ uniform vertical distribution, MESOPLUME
will not give correct results. Under stable flow conditions, for
example, plumes may not approach a uniform vertical distribution for
many kilometers downwind. This version of MESOPLUME should therefore
only be used at and beyond distances for which the assumption of uniform
vertical mixing is appropriate - often only nt distances of 100 km or
more. [It would be very simple to modify the present version, however,
to incorporate an optional Gaussian reflected vertical distribut ion so
a* to make MESOPLUME suitable also for near-field computations.]
166
-------
10
15
20
25
50
35
40
45
50
55
60
65
70
75
80
85
90
95
TABLE A-2
COMPARISON OF Cu/Q VALUES FOR MESOPLUME AND
TORMER WORKBOOK-COMPUTF.D VALUES
(u - 2.78 ¦ s~l. PGT Class = D, H = 1000 m,
unifora vertical distribution)
PGT
a (¦)
>'
550
780
1,000
1,220
1,420
1,620
1,820
2,020
2,200
2,400
2,600
2,775
2,975
3,200
3,375
3,550
3,700
3,850
4,000
Turner Workbook
(C u Q"1)
MESOPLUME
CC u Q"1)
7.25x10
5.11
3.99
3.27
2.80
2.46
2.19
1.97
1.81
1.66
1.53
1.44
1.34
1.2S
1.18
1.12
1.08
1.04
1.00
-7
7.77x10
5.94
4.11
3.47
2.83
2.50
2.19
1.97
1.78
1.64
1.50
1.42
1.31
1.25
1.17
1.11
1.06
1.00
0.95
-7
Fractional
Deviation
0.07
0.16
0.03
0.06
0.01
0.02
0.00
0.00
-0.02
-0.01
-0.02
-0.01
-0.02
0.00
-0.01
-0.01
-0.02
-0.04
-0.05
167
-------
Distance (km)
Figure A-8 MESOPLUME vs. the Conventional Gaussian Plume
168
-------
Iv^OMVk TtoWW* ¦»
APPENDIX B
MESOPUFF TECHNICAL DISCUSSION
B.1 Basic Mass Conservation Equations
MESOPUFF, a discretized variable-trajectory version of the
conventional straight-line Gaussian plume model, is designed to take
into account the spatial and temporal variations in the advection,
diffusion, transformation, and renoval mechanisms governing plume dis-
persion on regional transport scales. In MESOPUFF, a continuous plume
is nodeled by subdivision into a sufficient number of discrete "puffs"
of circular horizontal cross section, as illustrated in Figure B-l.
The conservation of pollutant mass in a puff transported a distance As
is expressed by the mass balance equation:
fiQ » / / / G(r, i, z) dr d0 dz
...» (B_1}
h I / /
C dr d8 dz
s+As
* h I i I C dr d6 dz
where r, 9, define points relative to the puff center in cylindrical
coordinates, G(r,9,z) (g m~3 s"1) is the rate of change (gain-loss) of
pollutant concentration C(r,6,z) (g af3), &Q (g s"1) is the resultant
rate of change of pollutant mass, Mid u(m s"1) is the wind speed. In
the MESOPUFF model, G(r,0,z) and u are constant for s to s ~ As,
where s is defined as the total distance a puff has traveled since it
was emitted.*
For a discrete puff lying below the mixing height H, the circularly
symmetric ground-level puff concentration C(r,0;s) is defined as
C(r,0;s)
2* 0y2(s) gj(z)
exp
2 ay2(s)J
g9(z)
(B-2)
where Q(s) is the pollutant mass flux and ay(s) the "radial" Gaussian
plume dispersion coefficient at distance s. The use of a "radial"
Gaussian dispersion coefficient is a convenient computational device,
*By contrast, in the MESOPLUME (plume segment) model, s is the current
distance of a plume segment endpoint from the emitting source, measured
along the plume axis; for temporally varying flows, the two definitions
can yield significantly different values of s.
Preceding pageuNank 171
-------
r (% *\
i N • j
u
*
Figure B-l The Variable Trajectory Puff Model Approach
172
-------
TCmnOlCXV* *tC
nothing Bore, The functions gj (:) and g? (:} are dependent upon the
vertical distribution of concentration in the puff. Replacing H by Hm,
the Mxinuai Mixing depth encountered by a puff (see Section B.8),
MESOPUFF perwits the user to specify one of two possible algorithms for
the distribution function gf:), namely
1) a uniform vertical distribution algorithm within Hm, such that
gl (i) = Hg, and gp(z) - 1.0, and
2) a Gaussian, multiple reflection algorithm where:
• if < 2 H , gj (z) = and g,, (;) is a function that
accounts fo? multiple reflection effects and
if o >
z —
H
gj(-) = H, and g2(i) = 1.0.
For regional-scale transport, e.g., at distances from 100 to
1,000 km fro« a source, either algorithm will produce substantially
similar results, as a rule, because at travel distances M00 km a,
is likely to be greater than 2H„,.
Using the uniform vertical distribution function (1), the ground
level puff concentration C(r,0;s) at distance s is
C(r,0;s)
Q(s)
2tt o ~(s) H
y i
exp
2 o "(s)
y
(B—3)
At distance s+As, the ground level concentration becomes
C (r, 0; s+As) =
Q(s-As)
2r a (s + As) H
y m
exp
-r
2 oy*"(s+As)
(B-4)
To ensure that a series of discrete puffs overlap along the variable
plume axis with sufficient density to approximate a continuous plume, an
individual puff (as described by Equation B-2), should not travel a
distance As greater than 2 oy in one time step At. Where wind speeds
are large enough that this condition would be violated, it is necessary
to subdivide further the time step used to interpolate puff concentra-
tions to grid points.
The MESOPUFF model solves the mass balance equation (B-l) indepen-
dently for both sulfur dioxide (SO2) and sulfate (SOiJ . The "gain"
functions for each species includeterms for the loss (gain) of S02
(SO,) by linear decay of S02 to SOk, and also include terms for dry
deposition of either species--see Sections B.5 and B.6.
173
-------
B.2 The Grid System
To facilitate the interaction of the MESOPUFF model with a number
of input and output routines and pre- and postprocessors, a simple
Cartesian coordinate system has been adopted for MESOPUFF. All spatial
model input data (.emission source inventories and meteorological fields)
are referenced to the same grid, called the "basic computational" grid.
To improve the resolution of the plume sampling function (see
Section B.4.3), MESOPUFF uses a sampling grid that is a subset of the
basic computational grid. The origin of the sampling grid may be placed
anywhere on the basic computational grid (other than on the northern or
eastern grid boundaries). The resolution of the sampling grid is a
mult iple of the resolution of the basic computational grid. Currently,
the maximum allowable size of both the basic computational and sampling
grids is 40 by 40 horizontal grid indices. Figure B-2 illustrates one
possible arrangement of basic computational and sampling grids. Note
that the grid index of the origin (1,1) is assigned Cartesian coor-
dinates x = 0, y = 0. The sampling grid in this example originates at
li = 2, j = 2) or (x = 1, y = 1) on the basic computational grid, and a
spatial resolution four tines finer than the basic computational grid.
B.3 Specification of Model Inputs
MESOPUFF is, at present, driven by the meteorological fields
produced by the MESOPAC MESOscale meteorology PACkage (Benkley and Bass
1979c), but an alternative meteorological preprocessing system with
appropriate grid resolution could be substituted by the user. MESOPUFF
requires the following fields (usually at hourly intervals) interpolated
to each model grid point:
• horizontal (u,v) wind components,
• mixing depth, and
• Pasquil1-Gifford-Turner (POT) stability class.
Use of MESOPAC model output by the MESOPUFF model is straightforward-the
Cartesian coordinate system is identical in both models. Note that the
horizontal (u,v) wind components are considered invariant with height.
Vertical wind components are not used. The fields of PGT stability
class are used by MESOPUFF in doing plume rise and plume growth calculations.
B.4 Puff Trajectory, Dispersion and Sampling Algorithms
The computational scheme of the MESOPUFF model has three distinct
functional elements: (1) a Lagrangian puff trajectory function, (2) a
puff dispersion function, and (3) a puff sampling function. The
174
-------
V-6.0
i-«
*-4.0
j-S
f "3.0
i-4
V - 2.0
j-3
T
Ay-Ad
_L
V" 1.0
1-2 '
y-0.0
i-1
*-0.0
i-1
I
r 1 i"
|__
- + + + -
- + + + -
i i i
1 I t
1
j j |
,|, „j ,
i i i
i i i
-+++-
j L_j__
J
i r i
- + + + "
_ _j_ -f* *4~ ~
~ ~h "4" "f* ~
i i i
1 1 i
—I i 1—
n
1 1 I
1 1 " I
- + ++"
_ i i J
TTT
1 1 i
i i i
mm mahout- mutual mm
mm hhi^-ib- am^iiiiiiM. mm
-+ + + "
i 1 1
I I I
1__| J--
-4-++-
-+++-
i i i
1 1 i
" + + + ¦
—|__J—
- + + + -
I I l
*-1.0
i-2
*-2.0 *-3.0 *-4.0 *-5.0
i -3 i-4 i-5 i-6
Ax-Ad-*-
Figure B-2 One Possible Arrangement of Sampling and Basic Computational
Grids
I
175
-------
Lagrangian trajectory function is used to advect the centerpoint of each
puff during a basic tine step. The radius of each puff is determined by
the puff dispersion function. Given the size and location of each puff,
the pufF sapling function computes the concentration exposure received
during the tine interval at each grid point by summing up the individual
puff contributions at each grid point.
B.4.1 Lagrangian Trajectory Function
This section describes how the centerpoint of a puff is advected
during a time step; it has been adopted (but largely verbatim) from
Hales et al. (1977, pages 15-18).
Let the components of the wind in the x and v grid directions along
the trajectory be designated by u(t;x(t),y(t)] and v[t;x(t),y(t)] ,
respectively.
where
t
x(t) • / u11' ;x(t' 1 ,y (t' 1 ]dt'
l 1
y(t) = ^ I v [ t' ;x (t' ) ,y (t ' ) ]dt'
x(°) = x0 ; y(0) = y0
u(0;x(0),y(0)) = u(xQ,y0)
v(0;x(0),y(0)) = v(x0,y0)
and XQ.yo are the initial coordinates (origin) of a trajectory of total
length s. (The grid space unit Ad is used in the program to convert to
nondisensional spatial units on the computational grid.) The vectors
u(t;—) and v(t;...) provide a Lagrangian description of the velocity
field. Increments of advection over a further time interval At are
therefore given by
Ax
bd
t*'.t
j u[t' ;x(t') ,y (t' )•] df ;
t
(B-5a)
t+At
/
t
&y = I v[t';x(t'),y(t*)] df . (B-5b)
176
-------
mt n.c***#Ck.oG* *k
After the new coordinates and y(t*At) are calculated for each
incremented puff centerpoint, a check is made to see if any of these
values are off the grid. If so, the affected puff is deleted from
further consideration.
The integration of the wind components over the trajectory is
approximated by a two-step iteration method involving bilinear inter-
polation in space and in tine (see Figure B-3).
The various
positions in Figure H-
-3 are best
defined by equations:
*1
•st
*{t) ~ u(t;x(t),y(t)]
At ;
(B-6a)
>'l
s
y (t) ~ v|t;x(t),y(t)|
At ;
(B-6b)
*2
s
*j * u(t*At;Xj,yjJ At
f
(B-6c)
y2
as
* v [ t ~ At;xj,yj] At
'
(B-6d)
x(t*At)
=
0.S(x(t) + x?] ;
(B-6e)
y(t*At)
sr
0.5[yft J ~ y,] .
(B-6f)
Thus, the new position (x(t*At),yft~At)) is found by adding the
vector average of two advection increments to the former position
(x(t),y(t)). The first advection increment is calculated by advecting
the plume increment for the entire interval At using the wind effective
at [x(t),y(t)] at time t. The addition of this increment to the posi-
tion (x(t1,y(t)] yields the posit ion (xj,yj). However, Ui ,y i) is not
asswaed to be the actual end of the advection increment because the wind
nay change along the trajectory. A second advection increment is cal-
culated using (xi.yj) as the starting point and the wind at that point
effective at the end of the time increment. The addition of this advec-
tion increment to position (xj.yj) yields position (x2,>'2). Then the
new puff position (x(t+At),y(t+At)) is taken to be the point halfway
along the line from (x(t),y(t)) to (x2,y2).
The bilinear interpolation by which the effective wind components
u(t), and \(;) are calculated works as follows. Let tn and tn+1 be the
effective times of the two gridded wind fields closest to time t. Time
interpolation weights tj and t2 are defined by:
t - t
" t < t < t , ; CB-7a)
Vl * *n n - n+1
tj = 1 - t2. (B-7b)
177
-------
I
Figure B-3 Calculation of the Trajectory of a Plume Segment Endpoint
178
-------
Next, the location of (x(t),y(t}) is noted on Figure B-3. The coor-
dinate values are between 3 and 4 for x(t) and between 2 and 3 for y(t).
Accordingly:
Xq - xCt) - 3 ; (B-8a)
X « 1 - X ; (B-8b)
p q
Yq - yft) - 2 ; (B-8c)
Y - 1 - Y . (B-8d)
p q
This yields
u|t) - tj Xp Yp u(tn;3,2) ~ t2 Xp Yp u(tn+1;3.2)
~ h \ Yp ~ t2 Xq Yp uCtn+1;4.2) (B-9)
~ *1 Xp YqU(V3,3) ~ t2 Xp Yq»(tn+1;3fS)
~ *1 Xq Yq »fV4'3) - t2 Xq Yq u(Vl;4,3) .
Similar equations hold for v(t), u(t»At), and v(t+At). Therefore, the
advection Equation B-5 is calculated from
Ax « (u(t) ~ u(t+At)}At; (B-lOa)
Ay = ^v(t) + v(t+At)}At. (B-lOb)
Hie subsequent positions are given by
x(t*At) = x(t) + Ax ; (B-lla)
>'(<:~ At) - y(t) ~ Ay . (B-llb)
179
-------
B.4.2 The Puff Dispersion Function
This section, describing the evaluation of the puff dispersion
parameters, has been taken largely verbatim, with modifications, from
Hales et al. (19??, pp. 23-24), The puff dispersion parameters ay and
0- are calculated for distances out to 100 km using plume growth
formilas fitted to the curves of Turner 1,1970). For distances greater
than 100 k> the plune growth rates given by Heffter (1965) are used.
'The growth of a puff with travel time t or along the plume trajectory
distance s froa the source is represented by
o (s»As)
y
da
a (s) * hs ~|—*•
y ds
(B-12)
s+fis/2 ;
o (t«-At)
o (t)
do
t+At/2 .
(B-13)
Similar equations for az are used. These terms allow for spatial and
temporal changes in stability class to be included, without violating
the entropy principle (puff centerpoint concentrations cannot increase
with downwind distance).
The integral formulas for ov and o, for travel distances less than
100 km are of the following forms
o (s,a) = Y s°'9 (meters);
y ®
(B-14)
c,(s,a) = s (meters),
(B-15)
Here a is a stability index designated as A, B, C, D, E, or F correspond-
ing to the PGT stability categories. The coefficients Za, and ba
are given in Table B-l as a function of stability index a and yield
values of ov and o2 in meters when s is specified in meters. Because
the integral formulas are not valid if the stability class changes over
the travel distance s, differentiated forms of equations B-14 and B-15
are actually used to carry out the computations.
do
_£
ds
0.9 Y
-0.1
(B-16)
180
-------
smfmmmmm.mmmcftrfeMmm-' «t
TABLE I-I
COEFFICIENTS FOR DISPERSION PARAMETER FORMULAS
Puff Growth Coefficients
Stability
Index a Y Z b
a a a
h 0.36 0.00023 2.10
i 0,25 0.058 1.09
C 0.19 0.11 0.91
0 0.13 0.5? 0.58
E 0.096 0.85 0.47
F 0.063 0.77 0.42
181
-------
dc,
2
_
b Z s
a a
-------
T*Ch»*XOC.v **'
-------
where SLj is the total number of puffs impacting the grid intersection
(i,j) during the tine step of interest.
B.5 Conversion of Sulfur Dioxide to Sulfate
MESOPUFF represents the conversion of SO; to SO,, as a simple
linear rate function--the changes in mass of each pollutant in one time
step at are given by the time-discretized equations:
AQn(SO,)
k Q (SO,) At
£ B - 2 2 a)
AQn(SO°) = -l.S k Qn(SO,l it (B-22b)
where Qn{SOy) is the mass of 50; at the beginning of the nth time step
and k 5(s~!) is the conversion rate constant. MESOPUFF uses the nominal
value kj = -S.56 x 10"* (that is, 21. per hour) suggested by Hales
et al. (1977), unless otherwise specified by the user.
Note that the conversion of SO; to SO,, represented by Equation B-22
is (1) independent of the vertical or horizontal distribution of SO;
within a plume and (2) independent of whether a plume lies above or
below the mixing height.
For the solut i ons to be adequate representations of the continuous
exponential equations (e.g., Q(S0;) = kjQ(S0;)), Q should not change
appreciably over the time interval it. Since this is frequently not the
case, the time interval St is optionally subdivided into as many sub-
intervals At' as are required to ensure that no more than a specified
fraction AQf of the mass Q is removed by decay duriig a subinterval it'.
Thi s process is described further in the next sect ion.
R.6 Dry Deposition of Sulfur Dioxide and Sulfate
fl
The changes of mass AQ (g) of each pollutant resulting from dry
deposition are given by:
AQn(S0,) = -(vd(S0,) Qn(S0:) At) / (B-23a)
AQn (SOj) = -(vd(S0=) Qn (S0°) it] / llm (B-23b)
184
-------
6^v«CWM6Nt/k«e»SA»Kw4^Ci-«CtCsG* «c
where Qn (SO. ) and Qn(SOiJ are the masses of S02 and SO,,, respectively, in
the puff at the beginning of the nth time step, vj(S02) and vjlSOjJ are
the deposition velocities of each pollutant, and H is the vertical depth
of the puff element (see Section B.8K MFSOPLIIMF uses the nominal
values vj(SO;) = 0.01 (m/s) and vjfSOiJ = 0.001 (m/s), as suggested by
llales et al. (197?), unless otherwise specified by the user.
As described before, MESOPUFF can optionally ensure that the mass Q
does not change within one subintervai At' by more than a user-specified
fraction &Qf. Taking into account both SO; removal mechanisms (decay
and deposition), the expression used by MLSOPUFF to specify the sub-
interval at* is
JAQ-(SO,) AQf(SO,) M
• -i^rispr
(B-24)
in other words, the subintervai size At' can optionally be chosen by
MESOPUFF at the beginning of each basic interval At such that no more
than a specified fraction AQf of S02 mass is removed by either mechanism
(or, no more than 2&Qf by both mechanisms combined) within At'. It
should be emphasized that the conversion and removal rates are not
affected by this restriction. Further discussion of the opt iona1
sampling mechanism is contained in Section B.9.
B.7 Plume Rise
The effective emission height he of each puff is computed as
he = h + hp, where hs is the stack height (m) and hp is the plume rise
(m). The plume rise equations used by the MF.S0PUFF model are those
described by Briggs (1975) for equilibrium (final) plume rise. For
unstable and neutral conditions with h' H (i.e., the plume does not
rise into an elevated stable layer)
hp = h' = 1.6 FI/J (3.5 u~' ; (B-2Sa)
For unstable and neutral conditions with h' > II (i.e., the plume
penetrates into an elevated stable layer)
-1 1 1/3
h = MIN'h' , (=£ ~ 18.75 F S ) 1 ; (B-25b)
185
-------
for stable conditions when u >_ 1.37 m/s
h . 2.6 f"3 S'"3 U-"3 ; (R-25c)
P
for stable conditions when u < 1.37 m/s
h =5.0 F1/4 S-3/S . CB-2Sd)
where:
F = buoyancy flux (m^/s^)
0 4
x* = 34.49 F for F > 55
f)
= 14.0 F " for F < 55
S = (t/T)(3e/3z)
g = 9.8ms
T = 290"K
36/3z = 0.137°K uT 1
H = mixing depth (hi)
z. = H - h (m)
b s
u = wind speed (m s *)
u = MAX (u , 1.37}
B.8 Treatment of Plume Fumigation
When the "uniform vertical distribution" option is elected,
MESOPUFF takes the spatial and temporal variations in the mixing height
into consideration to determine the height extent of ground level
impact:
(a) For a puff element which is emitted at an effective stack
height which is less than the mixing height, MESOPUFF assigns
the element a height of zero and immediately mixes the entire
mass of the puff uniformly through the mixing depth (the
mixing lid acts as a perfect reflector). Equations B-3 and
B-4 are then applied to compute ground level concentrations.
(b) For a puff element which is emitted at an effective stack
height which is greater than the mixing height, no effect of
the puff is felt at ground. Howe 'er, if subsequently the
mixing height becomes greater than the height of the puff
centerpoint, the entire puff element is immediately mixed
uniformly through the mixing depth, the puff height is set
equal to zero, and Equation B-3 and B-4 again apply.
186
-------
»t Sf!M*'I f.MNfM.rv.,'. <•*
Tlie mixing depth encountered by a puff element is likely to change
over time and space. Ml.SOPUFF assumes that a puff element residing in
the mixed layer is nixed through the maximum mixing depth Hm encountered
by the puff in its progression through the computational space-time
grid.
The restriction of a uniform puff element through a characterizable
depth such as H,,, is imperative for a one-layer model such as MESOPUFF -
so it is imperative for a one-layer model such as MESOPUFF - so that
parcels of material, once entrained, are not bifurcated in the vertical
when the height of the mixing lid changes (as sensed by a parcel in
transport]. Therefore, in the MESOPUFF model, a puff, once entrained,
is uniformly mixed through a realistic height, consistent with the
following assumptions:
(a) When convection ceases in the late afternoon, turbulent
energy in the afternoon mixed layer dissipates. Once the
energy of the large convective cells is completely damped, the
remaining small scale energy is not sufficient to yield a
vert ical entrainment rate comparable to the daytime entrap-
ment rate produced by large convective eddies.
lb) The stable 1 id capping the convective layer most 1ikelv
maintains some integrity after convection ceases, and acts to
prohibit r. y additional vertical entrainment of plume material.
(c) A modeling approach that redistributes mass into a shallower
layer when the mixing depth decreases will violate the second
law of thermodynamics.
Since no mechanism remains for turbulent transfer of pollutant
material either upwards or downwards once the mixing depth decreases, it
is asserted that the height most appropriate to define the extent of
vertical mixing of a puff is the maximum mixing height that the puff
encounters during any portion of its travel.
To schematically illustrate the interaction of the MESOPUFF
vertical distribution algorithm and the mixing depth progression
algorithm used by ML SOP AC (.Benkley and Bass (1979b), Benkley and
Schulman 1979)) consider the hypothetical example given in Figure A-7
of Appendix A.
Here, consider two puffs of material released at different times
from a source with an effective stack height of 250 m; the first puff
release is made at 0300 GMT on the first day, the second puff release at
0300 GOT on the second day. Assume that the hypothetical source is
located in the western United States so that 0000 GMT corresponds to
late afternoon--the usual time of maximum mixing depth.
187
-------
nMK'V **t . • »w.
The first puff is entrained at 0600 <1MT' on Day 1 as the nocturnal
boundary layer oscillates slightly upward in height. MESOPUFF immedi-
ately mixes the puff uniformly throughout this mixing depth.
This puff is mixed uniformly through 900 m by 0000 GMT on Day 2,
and remains nixed through 900 m subsequent to that time, even though the
height of the nixed layer collapses. The second puff, emitted at
0300 GMT on Day 2 is fumigated at 1500 GOT and is mixed uniformly
through 700 m by the end of the afternoon. Therefore, the first puff,
emitted at 0300 GMT on Pay 1, remains uniformly mixed through a deeper
vertical layer than the second puff; this is consistent with the greater
extent of mixing on the first day as compared to the second day.
This mixing depth/fumigation scheme exhibits at least two short-
comings.
• When a puff becomes bifurcated in the vertical by collapse of
the mixing lid, how should the puff radial dispersion rate be
characterized? Because MESOPUFF is a one-layer model, the
sane set of horizontal dispersion coefficients must he applied
to both the above-mixing-1 id and the below-mixing-1 id portions
of the bifurcated puff.
• The dry deposition algorithm (see Section B.6) removes SO; and
SO^ at a rate proport ional to the surface concentrat ion--this
loss is then redistributed throughout the entire vertical
extent of the plume. 11 would be more realistic to dist ribute
the loss only through the current mixing depth, but this is
not possible in a one-1aver puff element model.
It is felt that these limitations of the scheme are necessary for a
method that is both computationally practical and yet realistic for use
with a one-layer model. The inclusion of the fumigation cycle is very
important; the marked effects of plume fumigation on the distribution of
regional-scale ground level concentration patterns are illustrated in
Section 3.3.
B.9 Accurate Simulation of the Continuous Plume
Because MESOPUFF simulates a continuous plume by means of a series
of discrete puffs, care must be exercised so that puffs overlap suf-
ficiently to allow the plume to appear smooth and continuous to the
sampling function. There are two methods by which this can be accom-
plished: (1) by subdividing the distance As that a puff travels in one
time step, and (2) by subdividing effective distance As' between
consecutive puffs.
The first method allows for sampling each puff on the grid at a
rate ns times more frequently than that with which the puff is advected.
For ns > 1 MESOPUFF will choose a sampling time step uts = At ns~!
188
-------
where 4ts is the basic model time step. Therefore, the distance iss
traveled by a puff between consecutive plume samples is Ass = is n~ 1.
The second method allows for releasing puffs at a rate n times
more often than the basic tirae step; each puff released will nave an
initial mass of SO^, assigned as Qp{502) = Q(S02) np~ 5 where Q(S02) is
the mass of SO; emitted by a source in one 'basic' time step; the
corresponding distance iSp between two consecutive puffs is ftSp =
As* Hp'1 . MESOFUFF updates the Lagrangian trajectory of each puff
resident on the grid at each 'puff release' time step.
In practice, both options appear about equally effective in simula-
ting a smooth and continuous plume (See Section 5.3J. The puff release
method is slightly more expensive but may yield more accurate puff
trajectories in strongly curvilinear or temporally varying flows.
Both options can be exercised jointly; however, the larger of
(iip,ns) should be an even multiple of the smaller. As an example,
suppose the basic time step is 1 hour, with Hp = 2 and ns = 4; puffs
will be released every half-hour from each source, puff trajectories
will also be updated every half-hour, and puffs will be sampled every
quarter-hour, or twice per puff release step.
If the dynamic sampling step option is invoked in MESOPUFF (see
Section B.6) the actual samplings step n„ used by MESOPUFF during any-
basic model time step wi11 be taken as the larger of (a) the user-
specified np and (b) the dynamically-computed np.
Some guidance in selecting the appropriate puff release and
sampling rates for a MESOPUFF simulation is provided by Figure B-S,
which depicts the strong effect of increasing the distance between
adjacent puff centerpoints upon the ability of MESOPUFF to recreate a
continuous plume. It is seen that when adjacent puff centerpoints are
separated by no more than 2 a _, concentrations computed by MESOPUFF will
differ by no more than about ?2# from concentrations as computed by
MESOPUFF in the limit as puff separation distances go to zero; if the
centerpoints at adjacent puffs are separated by about 4a., deviations
can be a large as 50% or more. *
Because puffs are assigned zero width when released, it follows
that for a given combination of puff release rate, sampling rate, and
wind speed, a certain amount of time must elapse before adjacent puffs
have grown wide enough to provide the degree of overlap necessary to
simulate a continuous plume. But, because high rates of puff release
and puff campling quickly become cost-prohibitive (see Section 5.2.3J,
some near-field inaccuracies must be tolerated if long-distance trans-
port is to be modeled in a cost-effective manner. Therefore, if one
wishes to model plume impact at distances greater than 100 km, for
example, the puff relase rate and puff sampling rate should be chosen
so that, for the maximum expected wind speeds, adjacent puffs wi11
overlap sufficiently within, say a downwind range of 50-100 km.
189
-------
Puff Separation Distance (in units of oy)
Figure B-5 MESOPUFF Model: Percentage Deviation with Puff
Separation Distance
190
-------
£ «s w w *'t c *•**>. ¦ -
-------
Distance
(km)
5
10
15
20
25
30
35
40
45
50
55
60
65
70
75
80
85
90
95
100
TABLE B-2
COMPARISON OF Cu/Q VALUES FOR MESOPUFF AND
TURNER WORKBOOK-COMPUTED VALUES
(u = 2.78 m s"1, PCX Class - D» II = 1000 m,
uniform vertical distribution)
PGT Turner Workbook
oy Cm) (C uQ"1!
300 11.80 x 10"*
550 7.25
780 5.11
1,000 3.99
1,220 3.27
1,420 2.80
1,620 2.46
1,820 2.19
2,020 1.97
2,200 1.81
2,400 1.66
2,600 1.53
2,775 1.44
2,975 1.34
3,200 1.25
3,375 1.18
3,550 1.12
3,700 1.08
3,850 1.04
4,000 1.00
MESOPUFF Fractional
(C u Q-') Deviation
12.01 x 10"
0.02
7.47
0.03
5.18
0.01
4.02
0.01
3.28
0.00
2.80
0.00
2.43
-0.01
2.15
-0.02
1.94
-0.02
1.77
-0.02
1.62
-0.02
1.49
-0.03
1.39
-0.03
1.30
-0.03
1.22
-0.02
1.15
-0.03
1.09
-0.03
1.04
-0.04
1.00
-0.04
0.96
-0.04
192
-------
DtStonce ikm i
figure S-6 MF.SOPUFF vs. the Convent ional Gaussian Plume
193
-------
APPENDIX C
MESOGRID TECHNICAL DISCUSSION
C,1 Model Formulation
The horizontal advection, vertical diffusion,_1inear decay, and dry
deposition of sulfur dioxide (SO2) and sulfate (SO,,) species on regional
scales are represented in MESOGRID by a discrete-level numerical repre-
sentation of the continuous equations describing the mass conservation
of the respective pollutant species:
!5
at
9Cj
3x
Sv
3_
3z
3C
C —-
2 32
-klCl -f1C15(Z)H,Azk) ~ —
(C-l)
ac,
4
3t
aC2 f fS + 3_
3x % 3y 3z
3C,
(
z 3z
~ j kJC1 -f2C2S(z,H,Azk)
(C-2)
where
x is the east-west horizontal coordinate (m);
y is the north-south horizontal coordinate (m);
z is the vertical coordinate (m);
t is the time (s);
Cj, C, are the ambientconcentrations of sulfur dioxide (S02)
and sulfate (SOl.) respectively (g m 3);
is the source emission rate of SO2 (g m 2s :) within a
vertical cell of height Azj, (the SOij emission rate is
assumed to be zero);
u(x,y), v(x,y) are, respectively, the x and y components of hori zontal
wind velocity (m s *);
is the vertical eddy diffusivity (m2 s"1);
kj is the rate (s"1) of linear decay of S02 to SOu;
fj, f-, are the dry deposition rate functions (s"1) of S02 and
SOS, respectively. Dry deposition is considered only
when the height z of a parcel of pollutant is below the
mixing height H; and
6(z,H,Az^) = 1 for z <_ and k = 1; 5(z,H,Azk) = 0 for 2 > H or k ^ 1.
Preceding pag&Hank
197
-------
The set of equations used by MI.SOO.RII» invokes these important
assumptions:
• the aclvectr.t wind field is two-dimensional (hori zontal) only;
thus if the -t h!» i is to conserve mass strictly, it must he
driven by a horizontally nondivergent wind field.
• horizontal diffusion is negligible compared to advective
transport,
The MESOGRID model solves the set of disc retired advection-
diffusion equations on a three-dimensional grid with either two or three
vertical levels (as desired). The Egan-Mahoney (.1972) method of moments
is used to minimze pseudodiffusive errors by conservation of the lower
order moments of the pollutant distribution. The grid cell resolution
is user-specifiable (to a maximum of 26x26x3 cells). The horizontal
grid size is uniform; the vertical cell dimensions may vary between
levels. The choice of grid cell dimensions will reflect the following
considerations:
1) the overall size of the area to be modeled, and required
spatial resolution of concentration fields,
2} the advective and diffusive numerical stability criteria, and
3) the relative model execution costs as dictated by items 1 and
2 above.
C.2 The Grid System
MESOGRID is referenced to a simple Cartesian coordinate system for
easy interaction with a variety of input/output (I/O) routines, as well
as pre- and postprocessors. The MESOGRID grid cell system is illus-
trated in Figure C-l. Each grid cell is denoted by an ordered triple of
indices (i,j,k) where i is the east-west (x-direction) index (i=1 is the
westernmost cell), j is the north-south (v-direction) index (j=l is the
southernmost eel 1), and k is the vert ical (z-d i rect i on) index (k= 1 is
the lowest level). For example, the center of the southwesternmost and
lowest grid eel 1 (with indices 1,1,1) corresponds to the physical loca-
tion (x=0, y=0, z=Azj/2). All emission source data and meteorological
data input to the MESOGRID model are referenced to this grid system.
C.3 Specification of Model Inputs
The meteorological fields required to drive the MESOGRID model are
generated by the MESOPAC model (Benkley and Bass 1979c), a meteoro-
logical preprocessing package developed specifically for this purpose.
MESOPAC produces gridded hourly interpolated fields of:
• horizontal (u,v) wind components (m s~1) ,
-------
ir^ w « * "t -><•*% r>"» * **
| (i,j) «indicas kteotifyiof c«i number in horizontal
I
j 11,41
3 +
1
1
<2,41
+
J <1,31
2 +
!
I
<2,31
(3,3)
4-
+
J 11,21
1 •§¦
1
I
12,21
(3,2)
(4,2)
+
+
~
t
1 (1,1)
1
0 i—~
0
(2,1)
(3,1)
(4,1)
1
= Act
1
1
I
I
+ CNI
I
1
I
I
I
•
-~Ax = Act*-
I
Figure C-l A Schematic Representation of the ME50GRID Cartesian Coordinate
System
199
-------
at *9 **+¦*'*
• nixing depth (»), and
• Pasqui11-Gifford-Turner {PUT) stability class.
Interaction of the MESOGRIP model with the MESOPAC preprocessor is
straightforward--the coordinate system is identical in both models and
the grid point values of meteorological variables produced by >"'SOPAC
apply directly to the center of each MESOGRIP cell. At present, the
horizontal wind components generated by MESOPAC are treated as invariant
with height, although the MESOGRIP model itself can accommodate ver-
tically layered wind fields. The mixing depth and stability fields are
used by MESOGRIP to determine eddy diffusivity {) vertical profiles
(see a detailed discussion in Section C.4.j»).
Emission source locations are user-specified through a three -
dimensional Cartesian coordinate system fx,)*,; coordinates], and
MESOGRIP assigns the appropriate (i,j,k) grid cell, for example (see
Figure C-l), a source with x,y coordinates (O.h, 1.4) is assumed to fall
within a grid cell with horizontal indices (2,21. Acceptable ranges of
x and y are (-0.5 <_ x <_ iMX -0.51 and (-0.5 <_ y jmax -0.5) ; if x or y
values are specified beyond those limits, a MESOGRIP run will terminate
with the error message SOURCE NOT IX GRID. Source locations are
assigned vertical indices (k) as illustrated in i'igure C-2. Thus, a
source with an effective release height, ho = hs ~ hp = 600 m, where hs
is the stack height and hp the plume rise (see Section C.8), is assumed
to fall within the second (k=2) vertical layer.
Point source emissions Q(g s ') are assumed to be distributed
uniformly throughout a grid cell; an equivalent volume emission source
rate is defined by:
q = 2__ (C-3)
azk(M)^
where Q is the eel 1-averaged source emission rate (g m~ ?s~1 , Q is the
point-source emission rate (g s 1), A;j. is the vertical depth (m) of the
k-th cell, and Ad is the size of a horizontal grid cell. As an example,
if Q 1 103(g s"1),. &Z2 = lO3 m, and ad = 4 x 101* m (tvpical values) ,
then Q = 6.25 x 10"10 (g m~V>).
C.4 Solution of the Advection-Piffusion Equations
The tracer equations for each pollutant (Equations C-l and C-2)
are solved in finite difference representations. The solution proceeds
in two sequential steps. In the first step, pollutant mass is advected
by the wind field. Immediately thereafter, in the second step, pollu-
tant mass is diffused in the vert ical at a rate proport ional to the
layer values of the K_ field. These computational steps are described
next.
200
-------
p>v»:w.-Ts, «t se aoc*. »-fkv\ '.-o «
fa» = 2500f
Ait = 1000 m
Az, «500m
r» = 3
n = 2
n = 1
Height (ml
4000
2750
1500
1000
500
250
0
figure C~2 Vertical Structure of the ME50GR1D Computation Grid Assuming
Default Values
201
-------
?¦ mw***"-# *?>-•¦«-.- *»¦
C.4.1 Horizontal Advcction
Conventional finite-difference approximations to the two advection
term?, u:C/Sx and v?C/3y in Equations C-l and C-2, produce truncation
errors; these errors can introduce spurious, pseudodiffusive effects in
the predicted concentration terms. This artificial diffusion by the
numerical scheme can be, in some instances, orders of magnitude larger
than that resulting from the real atmospheric mixing process. Various
special finite-differencing methods have been proposed for control of
numerical diffusion (e.g., SHASTA--see Boris and Book 1973). In the
present model, the approach adopted is a unique mixture of both Eulerian
and Lagrangian approaches, called the Method of Moments (Egan and
Ma honey 1972). This method substantially reduces the pseudodiffusion
associated with numerical advection term by using one or more statis-
tical moments of the concentration distribution within a grid element
scheme. This section describes this procedure in detail. To keep the
explanation as simple as possible, the discussion of horizontal advec-
tion will be initially limited to the description of advection in the
east-west direction only.
Defining C as the average concentration within a grid element, the
forward-in-time, backward-in-space, finite-difference approximation to
the horizontal advection for a uniform wind field u > 0 (neglecting
diffusive, source, and sink terms) is
C™ = (1 - o) Cl ~ aC1 (C-4)
in,k m,k m-1,k
th
where T denotes the T tine step such that t = TAt, m and k denote the
rath (jth or jthj horizontal and the kt'1 vertical grid elements and
o = uAt/M, the ratio of the advection distance per unit time step to
the grid element dimension. The Courant advective criterion stability
requires that a <_ 1. After many time steps, the resulting horizontal
advection can be likened to advection with velocity u and upwind-
downwind mixing with a pseudodiffusive coefficient of magnitude
uAd(1 - o)/2. A substantial reduction in this diffusive effect can be
achieved if the first and second moments of the material within a cell
are calculated after each time step and used to adjust the amount of
material advected downwind during the next time step. Specifically, the
procedure calculates the first and second moments of the concentration
distribution after material has been advected into and out of a grid
cell. These moments are then used to define a simple rectangular
concentration distribution in the grid cell having the same amount of
total mass and identical first and second moments.
Suppose Cm denotes the relative displacement of material within the
Bth Ceil from the center of the eel 1, such that ranges from -0.5 at
the left hand to *0.5 at the right-hand extreme boundary of a cell. The
202
-------
zeroth, first, and centered second noments of the grid cell concentra-
tion distributions, corresponding to the mean concentration, the center
of mass, and the scaled distribution variance are given by
0.5
C(c ) d£
II HI
-0.5
0.5
CU) c dC /C
Rl Hi TO FR
CC-5)
-0.5
0.5
C(C ) (C - F )2 dCJC
n m m m m
-0.5
For simple rectangular concentration distributions, these integrals
are readily evaluated by summation for each grid element in terms of the
material distributions of the portions remaining and newly advected in
after each successive time step. It may be seen that the nature of the
downwind transfer depends on the value of the "portioning" parameter,
P = (F * a * 0.5 R - 0.5)/R (C-6)
s m m m 1 3
For Pm < 0 none of the material is advected into the m+1 cell. If
PB ^ 1 all of the material is advected into the downwind cell. For
1 > PB > 0, an amount of material PmCm is advected to the downwind cell,
and (1 - PH)Cm remains in the mttl cell. Considering now the general
case for the mth cell with inflow from the m-1 cell, outflow to the m+1
cell, and continuous source addition during the time step, the forward-
time computation procedure at grid element m can be represented as
~T-1
r ~a
C_ = C_ + C_ + Q^AT
CT+1 FT+1 = C F + C F ~ 0 At(a/2)
m m tv a a in v '
C.+1 (r2^+1 " CrlRr2 * 12(FI+1 " Fr)2] (C-7)
+ C [R2 ~ 12(FT+1 - F )2]
a1 a m a J
* Q„At[l ~ 12(F^+1 - 0/2)2],
203
-------
where subscripts r and a indicate quantities remaining and newly
advected in, respectively. This procedure completely eliminates pseudo-
diffusive effects for a simple rectangular concentration distribution.
Stall diffusive errors will remain when more complicated distributions
are advected.
In the MESOGRID model, the scheme described above has been general-
ized to two-dimension horizontal advection. figure C-3 illustrates in
top view the advection of a block of material with a uniform concen-
tration distribution. In analogy to the one-dimensional advection
procedure, one nay define (suppressing the horizontal grid element
subscripts (i,j)
P. = (F, + o. + 0.5 R. - O.S)/R,
ill i l
P. = (F. ~ n. ~ 0.5 R. - 0.5)/R.
111 J 1
re-si
Then (imagining a newly emitted parcel of material contained within one
grid cell), for the case illustrated in Figure C-3, where 0 < Pj,
P| < 1, the contribution to the new concentrations at each of the four
grid cells sharing the material after advection is
c?; . =
1+i.j
T
= c: .
i.j
>yl
cT+1 , =
l.J+l
- cT .
(1 - p
cT+; . , .
1*1,J+l
. cT .
I.J
p.p.,
1 3
cT+1 =
1,J
. CT .
1..1
(1 - P
(C-9)
The rules for advection for P < 0 or P > 1 follow those derived for
one-dimensional advection, as do the expressions for calculating the
first and second moments. In general the moments are evaluated from the
concentration distributions advected into a cell from more than one
adjacent cell. The computation procedure determines which neighboring
eel Is contribute to the moment calculations a-id computes IX , IC F ,
EC R 2, and then the new values for each element: m m m
m m
T+l
c = rc (c-io)
_ . IC F
,.T+1 mm %
= Tj— (C-11)
204
-------
»Vv-*;'
-------
rr
rc R " , T ,
" « _ 12 fT+1
,T+1
12
,i»l
(C -12)
€.4.2 Vertical !>i f fusion
The vertical diffusion terms S/3z(K,3C/3:) in Equations C-l and
C-2 are simulated by a conventional forward-tine, centered-difference
technique modified so that the vertical grid spacing can be variable.
(The MESOGRID model is easily modified to include more vertical levels
if computer resources and budgets permit.) In regions where parameters
or concentrations change rapidly with height, resolution and accuracy
can be improved with smaller vertical grid spacing. The vertical
diffusion calculation may be represented by
Cl = d - *k + 1 - Vl' Ck + \+l Ck+1 * Vl Ck-1* where the
horizontal subscript m is now suppressed, and
"V
(Kk * *WAt
k+1 - izk (Azk + A;k+1)
tK, + K )tt
vk-l = A-k (Mk ~ fiz^)
The diffusivities are calculated at the center of each grid
element, and denotes the depth of the kth grid element. The diffu-
sive computational stability criterion requires that the y's are less
than 0.5. The first and second moments of the horizontal distribution
are maintained in the diffusive transfer by two additional equations:
..T+l
T+l
Yk+1 " Yk-1
y C F
n-1 n-1 n-1
rT cT • rT trT
Ci Fi + Yi ^ C, , F. .
k k k+1 k+1 k+1
T-l
i+12 K - pr
1 " Yk+1 " Yk-1
" Yk+1 Ck+1
k+1
12
Ft . fT+1
k+1 k+1
(C-14)
rT
Vk-1 k-1
k-1
+ 12
FT - FT+1
k-1 k-1
^/cj+1 (C-15)
206
-------
t NV*IONMt ***. «f » **C>« * f C»*l£X (X, V *<
C.4.3 Calculation of Spatially Variable K, Profiles
MF.SOGRID computes spatially variable profiles of vertical eddy
di ffusivity for use in the vertical dispersion algorithms. The scheme
uses the nixing depth, wind speed, and stability class information
provided by MESOTAC to compute K, profiles in each of two altitude
regimes: (1) the surface layer of depth S = 0.1H and (2) the Ekman
layer extending from the top of the surface layer to the mixing height H,
To develop the profiles of K- in the boundary layer, it is assumed
that atmospheric pollutants diffuse in the same manner as heat. The
vertical eddy diffusivity in the surface boundary layer is then given
by:
0.35 u. 2
(C-16}
where
u# is the friction velocity,
z is height above ground,
L is the Monin-Obukhov length,
is the nondimensional potential temperature gradient.
By definition,
surface layer.
For stable conditions
j-l equals 1 for neutral conditions U- =
in the
1 ~ 4.7 ;
(C-17a)
and for unstable conditions (z/L > 0)
- 1/2
1 " is r
(C-17b)
Substitution of Equations (C-17a and C-17b) into Equation (C-16) yields
for stable conditions
0.35 ku„;
1 ~ 4.7 f
(C-18a)
207
-------
and for unstable conditions
-I1/2
O.SSu,; 1 - 15 fi
(C-18b)
The value of u. is taken to be 0.035 uf, where uf is the "free stream"
wind velocity (Blackadar 1962). If the MI-SOPAC model is used to compute
wind fields, the winds may not be significantly less than the free
stream wind speeds as long as the reference level chosen to describe the
wind fields falls typically within the upper two-thirds of the boundary
1aver,
If the Monin-Obukhov length L is assigned a value dependent on the
PGT stability class estimated by MESOPAC, the surface layer K- profiles
becotae stability-dependent as well. Table C-l depicts values of L 1
as a function of PGT stability class for a typical roughness length
zg = 25 cat, "as suggested by Colder (1972).
TABLE C-l
INVERSE MONIN-OBUKHOV LENGTH l."1 (m"1) AS A FUNCTION OF
PGT STABILITY CLASS, FOR ROUGHNESS
LENGTH z = 25 cm (AFTER GOLDER 1972)
PGT Stability Class
B CD
-0.11
-0.05
-0.01
0.01
0.05
To develop the Kz profiles in the Ekman layer, MESOGRID uses an
interpolation method suggested by O'Brien (1970) to determine diffu-
sivities between the top of the surface layer (S) and the top of the
nixed layer (H). Taking into account the physical requirement that the
first derivative of Kz be continuous with height in the Ekman layer,
O'Brien formulated the second-order equation
(z - H)'
(h - sr
S ~ c= - s)
(C-19)
3K
2
W
* 2
H - S
208
-------
(Ssmyw-'K AWrnfFCMMOv.^.*
where
z is the height at which K_ is to be determined
H is the height of the mixed layer (from MESOPAC)
5 is the height of the surface layer
K is the value of K at the top of the surface layer (from
Equation C-16) 2
3K,
i—^ is the derivative of K evaluated at S.
3 f *7
Here, it is assumed that Kz approaches zero at the top of the mixed
layer. Three stability-dependent equations for 3K,/3z|c are formulated
by differentiating Equation C-16 and evaluating at height S:
cor unstable conditions
3K
2
az
0.35u.
1 - 22.S
i-i/:
i - 15 r
(C-20a)
for neutral conditions
3K
2
3z
0.35u* ;
and for stable conditions
(C-20b)
8K
2
3l
0.35u#
1 +4.7r
(C-20c)
Table C-2 contains formulations for surface-layer K, (Equation C-16)
and 3K2/32|s (Equation C-20) appropriate to each PC.T stability class,
using the stability-dependent values of L"1 from Table C-l.
Equations C-16 through C-20 are used by MESOGRID to calculate
values of Kz at the mean height of each vertical grid element when-
ever the neteorological data is updated by MESOPAC. Because of the way
the MESOGRID model characterizes vertical diffusion the mixing depth, H
cannot be used directly as the depth through which pollutants diffuse in
the vertical. However, the effect of a mixing lid can be simulated
209
-------
TABLE C-2
STABIL1T¥-DEPENDENT L.TRESSIONS FOR Kz IN THE
SURFACE LAYER AND 3K
2
3z
PCT
Stability
K, (surface layer)
3K
2
32
0.35u#z{l + 1.65 z)I/2 0.35ut(1 ~ 2.48 S)(1 ~ 1.65 S)
-1/2
0.35u4z(1 + 0.75 z)1/2 0.35u.(l ~ 1.13 S)(1 ~ 0.75 S)
-1/2
0.35u#z(l ~ 0.15 z)1/2 0.3Su4(1 ~ 0.23 S)(1 ~ 0.15 S)
-1/2
0.35u#z
0.35u.
0.35u.z(l + 0.05 z) 0.35u.(l + 0.05 S)"
0.35u,z(l ~ 0.23 z)"1 0.35u„(l + 0.23 S)'
210
-------
*««•>* ffCHNtxcxv* *<-
approximately—MESOGRID sets K» in a grid cell essentially to zero if
the aean grid cell lies above the nixing height. Referring to the
exaaple in Figure €-2, if the nixing depth H lies between 1,000 and
2,750 meters, K- is set to zero in the uppermost layer but is left
unchanged in the two lower layers. But because diffusive transfer of
¦ass by two contiguous vertical grid cells is made proportional to the
average of the K- values for each cell (Equation C-13), without further
¦edification of the Kz profiles, mass transfer would take place between
Levels 2 and 3 even though Level 3 lies entirely above the mixing depth.
MESOGRID will correct partially for this behavior; the average of Kz for
Levels 2 and 3 will be set equal to zero if 1,000 <_ H < 1,500 m, but
will be unchanged if 1,500 <_ H < 2,750 m.
€.4.4 Numerical Stability Criteria
As described before, computational stability of the computational
scheme is ensured only if:
~ < 1 (C-2 1)
and
(Kk * Sil3 At <_0.5 (C-2 2)
AZk(iZk + Azk^
Equation (C-2 1), the advective stability criterion (Courant condition)
Beans that a parcel of mass must not be advected, in either horizontal
direction, through a distance greater than one grid space Ad within a
time step At. Equation (C-2 2) means that the effective vertical diffu-
sion velocity (K^ + Kk±l)/(Azk + Azk±l^ cannot be such that a parcel is
displaced vertically (by diffusion) more than a distance Azfc/2 within
At.
Given the maximum values of u and Kz at the beginning of each hour,
MESOGRID calculates the maximum time step allowed by the more stringent
of these stability criteria, except that the actual time steps used by
MESOGRID are constrained to be even fractions of one hour; thus the
possible choices for At are 5, 6, 7.5, 10, 15, 20, 30, and 60 minutes.
Note that, notwithstanding the above, MESOGRID does not currently allow
selection of a time step shorter than 5 minutes.
Table C-3 lists the maximum values wind velocity u and Kz values
(expressed as averages between adjacent levels, e.g., Kj 5 = (Kj + K2)/2,
for the nominal choices (Ad » 40,000 n, Azj = 500 m, A22 = 1,000 m, and
A13 = 2,500 m). Under stable conditions, the advective criterion will
normally dominate, whereas for unstable conditions (resulting in com-
puted Kz values up to 103 in test simulations), the vertical diffusion
211
-------
WILE C-3
MAXIMUM VALUES OF u AND Kj 5- *2 5 VS.
TIME STEP At» FOR NOMINAL MODEL PARAMETERS
Advective Diffusive
Model Tine Stability Stability
Step Criterion ___ Criterion
At («iin) (a s"1} Kt g|m2 s"!) K, 5(rn2 s"1)
5 135 625 2,916
15 45 208 972
30 22 104 486
60 11 52 293
Assumptions: M = 40,000 a
42.j = 500 m
= 1,000 m
iZj = 2,500 m
212
-------
stability criterion may dominate. Specific choices of horizontal and
vertical grid cell dimensions can stability criteria significantly
different from those of the hypothetical example. For instance, if the
desired near-ground vertical resolution dictates choice of a small
value, the maximum time step may be determined by the diffusion criteria,
whatever the stability class.
In the MESOGRII) model simulations, it was often found that for
unstable conditions, would be so large that the time steps necessary
to ensure computational stability were smaller than the minimum value
deemed reasonable for regional-scale transport (.',t = 5 minutes].
Because such small time steps result in very expensive computing costs,
a method, explained next, has been provided to give the user a degree of
control over the diffusion-specified time step by svstematically
decreasing very large values of K-. Such a procedure may be justifiable
if the MESOGRID model is used for regional-scale simulations (as dis-
tinct from near-source impacts). Because, in a grid model, the mag-
nitude of K. values determines the rapidity with which vertical mixing
will occur bv diffusion--as long as the reduced K- values lead to
uniform vertical mixing within travel distances that are small compared
to the regional transport scales of interest, the effect on regional-
scale concentration patterns should be minimal (unless ground-level
removal processes are of dominant importance).
The approach taken here forces the model to reduce large values of
K; to within a range consistent with both the advective stability time
step Ata (Equation C-2 l) and a user-specified "requested" time step Atr.
There are six distinct cases delineated in Table C-4 for modification or
noranodification of K-; these six cases are the six possible orderings of
Atr, .'.ta and the diffusive stability time step tj (Equation C-2 2) by
increasing magnitude. No modification of K- is required in Cases 1, 2,
and 3, because the advective stability criterion is at least as stringent
as the diffusive stability criterion, thus the model selects a time step
shorter than that required by the largest K. values. In Case 4, the
diffusive stability criteria is more stringent than the advective
stability criteria, but no Kz modification is necessary because
Atr
-------
TA1LE C-4
DETERMINING KHF.V P«SOGRID WILL MODIFY K PROFILES
Case Kuaber Possible Tine Step Orderings Modify K??
1 tt < at, < At No
a — d — r
2 at < At < At, No
a - r - d
3 At < fit < At . No
r — a — d
4 at < At, < at No
r — a — a
5 at, < At < at Yes, in accordance
d r ~ 3 with At
r
6 At, < &t < At Yes. in accordance
d a r with At
214
-------
C,5 Boundary Conditions
The MESOGRID model permits horizontal boundary concentrations of
S02 and S}~ to be specified independently for each vertical level and at
each east-west and north-south boundary face of the computational
region. The boundary value concentrations are assumed to occupy ficti-
tious grid cells adjacent to the boundary of the regular grid. These
concentrations are then advected into the grid using the procedure
described in Section C.4.1. These boundary values, which have default
values of zero, are assumed to be time independent.
The boundary condition describing mass transfer through the "top"
of the computational region is specified in terms of a reflection coef-
ficient. The reflection coefficient is user-specified and can range
from complete reflection of mass from the upper boundary (value of
unity) to total loss of material through the upper boundary (value of
zero). The boundary condition describing mass transfer out of the
bottom of the grid is user-specified through the 'deposition velocity'
concept. Section C.7 gives a more detailed description of the bottom
boundary condition.
C.6 Conversion of Sulfur Dioxide to Sulfate
MESOGRID represents the conversion of SO2 to SO* as a simple
linear rate function. The changes in mass of each pollutant in one time
step fit are given by:
&(f(S02) = kj Qn(S02) At (C-23a)
AQ(SO~) - -1.5 k2 Qn(S02) At (C-23b)
where Qn(S02) is the mass of S02 at the beginning of the nth time step,
and kj (s *) is the conversion rate constant. MESOGRID uses the nominal
value kj = -5.56 x 10~6, as suggested by Hales et al. (1977), unless
otherwise specified by the user.
It should be noted that the conversion of S02 to SO^ represented by
Equation C-23 is (1) made independent of the vertical or horizontal
distribution of S02 within a plume, and (2) independent of whether a
plume lies above or below the mixing height.
When Equation C-23 is applied to an individual grid cell after the
advection calculation and priorto the diffusion calculation, the changes
in concentrations of S02 and S0U resulting from linear conversion are:
ACn(S02) = kj Cn(S02) At (C-24a)
&Cn(S0p - -1.5 k2 Cn(S02) At (C-24b)
215
-------
•-. >•» ¦*** *,-K O* ~ \ * - V
where, for example, Cn(SO^) is the grid cell concentration of S0,> at
time step n, following the horizontal advection operation.
C.7 Pry Deposition of Sulfur Dioxide and Sulfate
The changes of mass M}nfg) of each pollutant at time step n result-
ing from dry deposition are given by:
AQn(SO^)
[vd(SOp qjfsojl Atj/Az
(C-2Sb)
where Qn(SOj) and Qn(SO~) are the total masses of SOj and S0„, respectively,
in a vertical column at the beginning of time step n CQi is the mass in
the lowestcel1); v. (S02) and v^(SOZ) are the deposition velocities of
S0Z and SO" respectively; and :j is the vertical depth of the the
lowest grid layer. MESOGRID uses the nominal values vjfSOj) =
0.01 (it s 1) and vj(SO^) = 0.001 (m s 1), as suggested by Hales et al.
(1977), unless otherwise specified by the user.
There is an important difference in the way mass is depleted by
MESOGRID as compared to the plume and puff models; although the amounts
of SO/ and SO^ removed by dry deposition are the same in each model for
vertically uniform distributions, MESOGRID depletes mass in only the
lowest of its vertical layers rather than throughout the entire mixing
depth (as do the other models). However, since MESOGRID exercises the
dry depostion algorithm before the vertical diffusion algorithm, mass
loss in a cell in the lowest level is redistributed through cells in the
upper layers (as long as they are below the mixing lid). Therefore, if
the effective "diffusion velocity" (see Section C.4.4) is large, the dry
deposition algorithm in MESOGRID yields the equivalent results as the
algorithm in the other models.
C.8 Plume Rise
The effective emission height hc of each point source is computed
as he = hs + hp, where hs is the stack height (m) and hp is the plume
rise (m). The plume rise equations used by the MESOGRID model are those
described by Briggs (1975) for equilibrium (final) plume rise.
For unstable and neutral conditions with h' < H (i.e., the plume does
not rise into an elevated stable layer).
h
P
, , , , 1/3 . *,2/3 -i
h' = 1.6 F (.>. 5 x*) u
(C-26i)
216
-------
EVv^Owevr*,. «t»**o* * ChNOuOO* *c
For unstable and neutral conditions with h* > H (i.e., the plume
penetrates into an elevated stable layer)
h = MIN
P
h«» (z^ + 18.75 F ii"1 S"1)
1/3
for stable conditions when u > 1.37 a s
-1
; (C-26b)
2.6 F "3 S -l/3U "1/3;
(C-26c)
for stable conditions when u < 1.37 a s"
S.O F "" S -3/S;
(C-26d)
where:
T
30/3z
H
2b
u
4 -3
buoyancy flux (a s )
34.49 F0*4
14.0 F0-625
Cg/T) (96/3z)
9.8 m s'2
290 °K
for F > 55
for F < 55
-1
0.0137 K
nixing depth (ib)
H - h (¦)
S -1
wind speed (m s )
MAX {u , 1.37}
217
-------
V, .>t -v» i,.. ,. *.
APPENDIX I)
SCALE ANALYSIS FOR TERMS IN MASS CONTINUITY EQUATION
Consider the mass continuity equation with decay:
— = -?• (vC) - kC + CV- v
3t
where
3C
3t
is the time rate of change of concentration,
-V-(vc) is the advective term in flu* form,
-kjC is the decay term, and
CV« v is the divergence term.
Using worst-case scaling values, an order of magnitude estimate
for the time scale of the advective term is approximately:
L 2 x 1O^m _ ,
T , 1- 77 r—i ^ J hours
ad U 2m/s
where L is a length scale typicallv containing most of the lateral
concentration gradient of the plume, and U is a typical stagnation wind
speed velocity. By comparison, the time scale for decay is much longer:
d ' In(0.98) % 50 hours
the e-folding time for a rate of decay of 2% per hour.
But, for the divergence term taken as 10 ** sec *, the time scale
for mass change resulting from the absolute divergence of the horizontal
wind field is given by:
T,. 3 hours
div
or of the same order as the advective time step. Therefore, the results
of the MESOGRID model (in common with any Eulerian grid model) can be
significantly influenced by divergence presented in the wind field.
Preceding page^blank
------- |