United States Atmospheric Research and Exposure EPA/600/8-89/041
Environmental Protection Assessment Laboratory March, 1989
Agency Research Triangle Park, NC 27711
Research and Development
v>EPA
User's Guide to the
Complex Terrain
Dispersion Model Plus
Algorithms for
Unstable Situations
(CTDMPLUS):
Volume 1. Model
Description and User
Instructions
-------
vxEPA
United States
Environmental Protection
Agency
Atmospheric Research and Exposure
Assessment Laboratory
Research Triangle Park NC 27711
EPA/600/8-89/041
March 1989
Research and Development
User's Guide to the
Complex Terrain
Dispersion Model Plus
Algorithms for
Unstable Situations
(CTDMPLUS):
^
Volume 1. Model
Description and User
Instructions
-------
EPA/600/8-89/041
March 1989
USER'S GUIDE TO THE COMPLEX TERRAIN DISPERSION MODEL
PLUS ALGORITHMS FOR UNSTABLE SITUATIONS (CTDMPLUS)
VOLUME 1: MODEL DESCRIPTION AND USER INSTRUCTIONS
by
Steven G. Peny*
Atmospheric Sciences Modeling Division
Atmospheric Research and Exposure Assessment Laboratory
Research Triangle Park, NC 27711
Donna J. Burns and Lucy H. Adams
Computer Sciences Corporation
Research Triangle Park, NC 27709
Robert J. Paine,** Michael G. Dennis,** Michael T. Mills**
Environmental Research and Technology, Inc.
Concord, MA 01742
David G. Strimaitis,** Robert J. Yamartino,** Elizabeth M. Insley**
Sigma Research Corporation
Westford, MA 01886
On assignment from the National Oceanic and Atmospheric Administration,
VS. Department of Commerce
Authors of the original CTDM user's manual (EPA/600/8-87/058a), much of
which remains in this manual for CTDMPLUS.
ATMOSPHERIC RESEARCH AND EXPOSURE ASSESSMENT LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
US. ENVIRONMENTAL PROTECTION AGENCY
RESEARCH TRIANGLE PARK, NC 27711
-------
NOTICE
All the information in this document has been funded wholly or in part by the U.S. Environmental
Protection Agency; the portions on the original CTDM were funded under Contract No. 68-02-3421 to Envi-
ronmental Research and Technology, Inc. This document has been subjected to the Agency's peer and
administrative review, and it has been approved for publication as an EPA document. Mention of trade
names or commercial products does not constitute endorsement or recommendation for use.
-------
ABSTRACT
The Complex Terrain Dispersion Model PLus Algorithms for Unstable Situations(CTDMPLUS) is a
refined air quality model for use in all stability conditions for complex terrain applications. It contains, in its
entirety, the technology of CTDM for stable and neutral conditions. However, CTDMPLUS can also model
daytime, unstable conditions, and has a number of additional capabilities for improved user friendliness. Its
use of meteorological data and terrain information is different than current EPA models; considerable detail
for both types of input data is required and is supplied by preprocessors specifically designed for
CTDMPLUS. CTDMPLUS requires the parameterization of individual hill shapes using the terrain pre-
processor and the association of each model receptor with a particular hill (except for receptors in flat ter-
rain, which CTDMPLUS can also model).
In modeling stable to neutral conditions, a central feature of CTDMPLUS is its use of a critical
dividing-streamline height (Hc ) to separate the flow in the vicinity of a hill into two separate layers. Flow in
the upper layer has sufficient kinetic energy to pass over the top of the hill, while streamlines in the lower
layer are constrained to flow in a horizontal plane around the hill. Two separate components of
CTDMPLUS compute ground-level concentrations resulting from plume material in each of these flows:
LIFT handles the flow above Hc , and WRAP handles the flow below Hc .
In modeling unstable (convective) conditions, the model relies on a probability density function
(PDF) description of the vertical velocities to estimate the vertical distribution of pollutant concentration.
Terrain distortions'of plume parcel trajectories are accounted for, as are deflections of the daytime mixed
layer height.
Hourly profiles of wind and temperature measurements are used by CTDMPLUS to compute plume
rise, plume penetration, convective scaling parameters, the value of Hc , and the Froude number above Hc .
In stable/neutral conditions the profiles of turbulence data ( fle or a,, and « values) are used to compute
plume a y and a z values at plume height.
The model calculates on an hourly (or appropriate steady averaging period) basis how the plume tra-
jectory (and, in stable/neutral conditions, the shape) is deformed by each hill. The computed concentration
at each receptor is then derived from the receptor position on the hill and the resultant plume position and
shape.
This user guide is divided into two volumes: Volume 1 describes the model and how to use it, while
Volume 2 contains code listings. Auxiliary user manuals describe the terrain and meteorological preproces-
sors. The CTDMPLUS technical descriptions of algorithms for stable/neutral conditions remain virtually
unchanged from those of the original CTDM manual (Paine et al., 1987). However, due to improvements in
user friendliness, sections of this CTDMPLUS manual pertaining to user instructions and model structure
have changed significantly from those for CTDM.
ui
-------
TABLE OF CONTENTS
Abstract in
Figures vii
Tables ix
Symbols and Abbreviations x
Acknowledgements xvi
1 INTRODUCTION 1-1
1.1 BACKGROUND 1-1
1.2 MODEL APPLICABILITY AND TECHNICAL LIMITATIONS 1-2
1.3 OVERVIEW OF CTDMPLUS SOFTWARE COMPONENTS 1-3
13.1 Terrain Preprocessor 1-3
1.3.2 Meteorological Preprocessor, METPRO 1-5
1.3.3 Receptor Generator, RECGEN 1-6
1.3.4 Complex Terrain Dispersion Model PLus Algorithms for Unstable Situations,
CTDMPLUS 1-6
1.3.5 Interactive Set-up Program, SETUP 1-6
13.6 Concentration Postprocessors 1-7
1.3.7 Menu Driver Program 1-7
1.4 CTDMPLUS COMPUTER CONSIDERATIONS 1-7
1.4.1 Software Considerations 1-7
1.4.2 Hardware Limitations 1-8
2 TECHNICAL DESCRIPTION 2-1
2.1 OVERVIEW OF THE CTDM PLUS ALGORITHMS FOR UNSTABLE
CONDITIONS (CTDMPLUS) 2-1
2.1.1 Pollutant Releases into Stable/Neutral Layers 2-1
2.1.2 Pollutant Releases into Unstable Layers 2-6
2.2 CALCULATION OF METEROLOGICAL VARIABLES 2-11
2.2.1 Mixed-Layer Height 2-12
2.2.2 Wind Speed 2-13
2.2.3 Wind Direction 2-14
2.2.4 Potential Temperature Gradient 2-15
2.2.5 Turbulence Intensities 2-16
2.3 PLUME RISE CALCULATIONS 2-17
23.1 Momentum Rise 2-17
2.3.2 Neutral/Unstable Buoyant Final Rise 2-17
2.3.3 Neutral/Stable Buoyant Final Rise 2-18
2.4 DISPERSION PARAMETERS FOR STABLE/NEUTRAL CONDITIONS 2-19
2.4.1 Plume Spread Due to Turbulence 2-19
IV
-------
TABLE OF CONTENTS (continued)
2.4.2 Allowance for Source-Induced Effects 2-22
2.5 THE LIFT COMPONENT 2-23
2.5.1 Overview of LIFT 2-23
2.6 THE WRAP COMPONENT 2-31
2.6.1 Overview of WRAP 2-31
2.6.2 WRAP Structure 2-31
2.7 RECEPTORS NOT TREATED BY LIFT OR WRAP 2-37
2.7.1 Receptors Not Influenced by Hills 2-38
2.7.2 Receptors Upwind of the Impingement Point 2-38
2.8 PLUME PENETRATION INTO ELEVATED STABLE LAYERS 2-40
2.9 THE VERTICAL DISTRIBUTION FUNCTION 2-41
2.9.1 The Probability Density Function (PDF) 2-41
2.9.2 Terrain Influences On the PDF 2-47
2.10 TERRAIN INFLUENCE ON MIXING HEIGHT 2-53
2.10.1 Estimate of the Elevated Temperature Gradient 2-55
2.11 THE HORIZONTAL DISTRIBUTION FUNCTION 2-56
3 USER INSTRUCTIONS FOR RUNNING CTDMPLUS 3-1
3.1 OVERVIEW 3-1
3.2 INPUT DATA REQUIREMENTS 3-6
3.2.1 General Program Specifications 3-7
3.2.2 Terrain Data 3-11
3.2.3 Receptor Data 3-11
3.2.4 Meteorological Profile Data 3-11
3.2.5 Meteorological Surface Data 3-15
3.2.6 Meteorological Rawinsonde Data 3-15
3.2.7 Hourly Emissions Data 3-15
3.3 CTDMPLUS OUTPUT FILES 3-15
3.3.1 CTDMPLUS Output Listing 3-15
33.2 Concentration File 3-19
3.4 ADDITIONAL COMPUTER NOTES 3-20
3.5 CTDMPLUS SUBROUTINE STRUCTURE 3-22
4 USER INSTRUCTIONS FOR AUXILIARY PROGRAMS 4-1
4.1 RECEPTOR GENERATOR 4-1
4.2 INTERACTIVE SETUP PROGRAM 4-2
43 GRAPHICAL CONCENTRATION DISPLAYS 4-2
43.1 CHIRET 4-11
-------
TABLE OF CONTENTS (continued)
4.3.2 CHIDIS 4-11
4.33 CONTOUR 4-15
5 INSTRUCTIONS FOR USING THE MENU DRIVER 5-1
5.1 INTRODUCTION 5-1
5.2 USER INSTRUCTIONS 5-1
5.2.1 File Naming 5-1
5.2.2 Directory Setup 5-1
5.23 Initiating the Program 5-5
5.2.4 Terrain Preprocessor Programs 5-5
5.2.5 Receptor Generator 5-9
5.2.6 Meteorological Preprocessors 5-9
5.2.7 Running CTDMPLUS 5-12
5.2.8 Concentration Postprocessors 5-15
References Ref-1
Appendix A Details of LIFT and WRAP Algorithms A-l
Appendix B Test Case for CTDMPLUS B-l
VI
-------
FIGURES
2-1 Idealized picture of how the dividing-streamline plane and the plane of stagnation stream-
lines "cut" into a plume 2-4
2-2 Depiction of plume behavior as it passes over a hill as seen from the side 2-5
2-3 Depiction of plume behavior as it is deflected around a hill as seen from above 2-7
2-4 Depiction of plume behavior as it passes over a hill as seen from above 2-8
2-5 Idealized stratified flow about hills indicating domain of individual CTDMPLUS compo-
nent algorithms 2-25
2-6 Illustration of the relationship between the crosswind-average concentration profiles at SQ
and s , and the plume from one of the many point-source elements representing the flux of
material across the plane at s0 2-27
2-7 Illustration of terrain effect on the vertical distribution of plume material above Hc as
modeled in LIFT 2-29
2-8 Typical streamline patterns in two-dimensional flow around an elliptical cylinder 2-32
2-9 Top view of a plume in two-dimensional flow around a hill 2-34
2-10 Sketch of the flow around an ideal cylinder of elliptical cross-section 2-36
2-11 Treatment of receptors upwind of the impingement point 2-39
2-12 Illustration of the use of image sources to model possible paths 2-43
2-13 Illustration of the bi-Gaussian distribution of Pv 2-45
2-14 Comparison of normalized crosswind integrated concentration contours 2-46
2-15 Illustration of the intercept method used for finding a second estimate of the vertical
velocity that carries a parcel in a deformed flow over a hill 2-49
2-16 Illustration of the weighted-average approach 2-50
2-17 Illustration of the linear interpolation method used iteratively 2-51
2-18 Fluid modeling results from Perry and Snyder (1989) showing deflections in the z,
streamline 2-54
2-19 Typical profile of mean potential temperature 2-57
3-1 Interaction among components of the CTDMPLUS system 3-2
3-2 Outline of the main program CTDMPLUS 3-23
3-3 Outline of the subroutine SEQMOD 3-28
3-4 Outline of the subroutine DAYCALC 3-30
4-1 Portion of a sample interactive session from the RECGEN program 4-3
vu
-------
FIGURES (continued)
4-2 Receptor display from the interactive RECGEN program 4-4
4-3 Example of a receptor file generated by the interactive RECGEN program 4-5
4-4 Display from a sample session of SETUP: "OPTIONS" input file for METPRO 4-7
4-5 Display from a sample session of SETUP: "PROFILE" input file for METPRO 4-8
4-6 Display from a sample interactive session of SETUP: "SURF1" input file for METPRO 4-9
4-7 Display from a sample interactive session of SETUP: listing of files required by
CTDMPLUS : , 4-10
4-8 Sample interactive session for program CHIRET, a portion of the concentration display
postprocessor 4-12
4-9 Sample screen display from the CHIDIS program showing user options for concentration
displays 4-13
4-10 Sample screen display from the CHIDIS program showing previously displayed receptors
plus the current one and the concentration value 4-14
5-1 Directory setup 5-4
5-2 Screen display from the menu driver system, showing the main selection menu 5-6
5-3 Screen display from the menu driver system, showing the terrain preprocessor menu 5-7
5-4 Sample display from the menu driver system, showing the FITCON options screen 5-8
5-5 Sample screen display from the menu driver system, showing the HCRIT options screen ... 5-10
5-6 Screen display from the menu driver system, showing the meteorological preprocessor
menu 5-11
5-7 Sample screen display from the menu driver system, showing the selection of a file for use
with METPRO 5-13
5-8 Sample screen display from the menu driver system, showing the list of files selected for
use with CTDMPLUS 5-14
5-9 Sample screen display from the menu driver system, showing the concentration postpro-
cessor menu 5-16
5-10 Sample screen display from the CONTOUR program, showing the concentration isopleths
and hill contours 5-17
vui
-------
TABLES
3-1 Format of the "FOPTIONS" FITCON options file 3-3
3-2 Format of the "HOPTIONS" HCRIT options file 3-4-
3-3 Contents of the "CTDM.IN" file 3-8
3-4 Surface roughness length (meters) for land-use types and seasons 3-10
3-5 Format of the "TERRAIN" input data file 3-12
3-6 Format of the "RECEPTOR" input data file 3-13
3-7 Format of user-created "PROFILE" file 3-14
3-8 Format of the "SURFACE" file 3-16
3-9 Format of the "RAWIN" file 3-17
3-10 Variable emissions input format 3-18
3-11 Format of text concentration output file 3-21
5-1 Components in the CTDMPLUS menu driver 5-2
5-2 Standard CTDMPLUS file extensions 5-3
IX
-------
LIST OF ABBREVIATIONS AND SYMBOLS
ABBREVIATIONS
CTDM Complex Terrain Dispersion Model
CTDMPLUS Complex Terrain Dispersion Model PLus Algorithms for Unstable Situations
CWIC Crosswind Integrated Concentration
ERF Error Function
EPA Environmental Protection Agency
GLC Ground-Level Concentration
HDF Horizontal Distribution Function
PC Personal Computer
PDF Probability Density Function
RTDM Rough Terrain Diffusion Model
VPTG Vertical Potential Temperature Gradient
SYMBOLS
a
B
B0
B i , Bz
C
d
dQ
de/dz
f
F.
Fr
F,Fb
Major axis length of an ellipse, m (also denoted as La )
Term used in Fz formulation
Minor axis length of an ellipse, m (also denoted as Lb )
Factors involved in Bl,B2 terms in WRAP calculation
Function oih/L used in computation of wind direction adjustment with height
Parameter used in analytical expression for convolution integral /
Factors involved in vertical term in WRAP concentration calculation
CTDMPLUS modeled concentration ( ng/m3)
Stack-top inner diameter, m; also crosswind distance from plume centerline (m)
Source strength of a point source element (g/s)
Vertical potential temperature gradient (K/m)
Coriolis parameter = 1.458 x 10"1 sin(latitude), sec"1
Buoyancy flux parameter
Hill Froude number
Buoyancy flux (mVs3)
Vertical distribution function in Gaussian plume equation
Acceleration due to gravity = 9.8 m / s2
-------
LIST OF ABBREVIATIONS AND SYMBOLS (continued)
G Soil heat flux, W/ m 2 ; also, a constant used in determining vector wind speed from scalar
wind speed and oe; also denotes Green's function
Gm Maximum probability of Pa
h Stable mixed layer height, m; also final plume height, m; also hill height function (m)
h R Height of receptor above reference base height (m)
h g Modified height of receptor (relative to plume height) due to terrain effects (m)
hs Height of the stack top above the common stack base (m)
H Height at the top of the hill above common stack base, m; also the sensible heat flux
(WAn2)
H c Critical dividing streamline height (m)
i y Horizontal component of the turbulence intensity
i z Vertical component of the turbulence intensity
/ Convolution integral of Green's function with the hill height function
k von Karman constant (0.4)
K y Eddy diffusivity, horizontal com ponent (m 2 / s )
K z Eddy dif fusivity, vertical com ponent ( m 2 / s )
I Mixing length, m; also crosswind distance from plume centerline to receptor (m)
/ Mixing length in neutral conditions (m)
/ s Mixing length in stable conditions (m)
L Monin-Obukhov length (m)
L a Major axis length of an ellipse (m-also denoted as a)
L Minor axis length of an ellipse (m-also denoted as b)
L , L x Parameters used in analytical expression for convolution integral /
L x Length scale of hill in along-wind direction (m)
L y Length scale of hill in cross-wind direction (m)
m Modified inverse length scale =n(l + l^/lj)1/2; also trajectory slope in PDF
n Ratio of N/u (or1)
N Brunt-Vaisala frequency (s"1 )
N0 Value of TV atS0 (flat terrain value)
P Plume penetration factor
XI
-------
LIST OF ABBREVIATIONS AND SYMBOLS (continued)
P Probability density function of vertical lift
p Perturbation pressure (Boussinesq flow)
Q Emission rate (g/s)
r Ratio of ellipse major axis length to minor axis length, r - a /b
R Displacement from a reference point, m; (R2 = x2 + y2 + z2)
R , Richardson number
R L Parameter used in analytical expression for convolution integral /
s, -S Square of the Brunt-Vaisala frequency, stability parameter ( s"2 )
s0 Downwind distance at which the plume impinges upon or is deflected by terrain (m)
ST Speed of the incident flow at the tower location ( n T, v T ), (m/s)
S. Speed of the incident flow approaching a hill (m/s)
t Travel time (s)
Ta Ambient air temperature (K)
T Factor for streamline distortion hi the vertical (height correction factor)
T i Distortion factor due to terrain effects in the horizontal direction
TL Lagrangian time scale (s)
T Lo Value of TL at s0 (flat terrain value), (s)
T s Stack gas temperature (K)
T Factor for wind speed distortion due to terrain influences
T y Terrain-effect factor equal to ratio of a'y/a'y.
T z Terrain-effect factor equal to ratio of c * / c \
Tau Terrain-effect factor equal to ratio of ou/ouo
Taa Terrain-effect factor equal to ratio of aw/auo
U Total scalar wind speed (m/s)
u Scalar wind speed, along-flow component (m/s)
u. Friction velocity (m/s)
u g Geostrophic wind speed (m/s)
u s Stack-top wind speed (m/s)
u Vector wind speed (m/s)
xu
-------
LIST OF ABBREVIATIONS AND SYMBOLS (continued)
u', u', w' Perturbation wind velocities (Boussinesq flow) in downwind, crosswind, and vertical direc-
tions, respectively (m/s)
v ' Cross-flow wind speed component (m/s)
w Vertical component of the wind speed (m/s)
w Mode of vertical wind component (m/s)
w. Convective velocity scale (m/s)
wd Mean downdraft velocity in unstable conditions (m/s)
w Stack gas exit velocity (m/s)
X Non-dimensional alongwind distance
x Alongwind coordinate (m)
xQ Alongwind coordinate at S0 , the impingement point of the approach flow on the hill (m)
x m Length parameter used in analytical expression for convolution integral /
xr,yr,zr Position in space of the receptor
x,, y,, z, Position in space of the plume release
y Distance in cross-wind direction (m)
y K Crosswind distance of receptor from plume centerline (m)
Y K' Modified crosswind distance of receptor from plume centerline due to terrain effects (m)
z Measurement height, m ; also receptor height (m)
za Surface roughness length (m)
zb zt- stack height (m)
zt Convective (daytime) mixed layer height (m)
zm Height of a streamline above the base of the hill far upwind (m)
z' Height above terrain surface (m)
a Wind speed shear (vertical component) (s~l)
a i, a 2 Parameters used in analytical expression for convolution integral /
a , Direction of the incident flow approaching a hill (deg)
(3 Strain function from theory of Hunt and Mulhearn (1973) used in computation of T y, T,,
also, rotation angle between wind flow and stagnation streamline; also the asymmetry factor
in the PDF model
6z Vertical deflection experienced by a streamline due to terrain influences (m)
A a Wind direction change within the mixed layer (degrees)
Xlll
-------
LIST OF ABBREVIATIONS AND SYMBOLS (continued)
A h Plume rise (m)
e Variable used to relate between scalar and vector wind speeds, a function of 0 9
A 9 Jump in potential temperature at the top of the mixed layer (K)
r) Height of a streamline above the surface of the hill; also denotes deflection of streamlines in
the vertical (Appendix A)
Y Constant of proportionality between the mixing length in a stable atmosphere and
au,/;V;Y = 0.52
T Constant of proportionality between the mixing length in a neutral atmosphere and
z;r-0.36
H , v Elliptical coordinates
u0 Value of [i along the boundary of an ellipse
\JL ,, v, Elliptical coordinates of the source
li T , VT Elliptical coordinates of the point where wind speed and direction are measured (tower)
4>T Wind direction measured at the tower location ( UT , VT ), (m/s)
Ym Stability correction for wind profile formulation in the surface layer
₯, Stream function through the source
p o Initial unperturbed fluid density in Boussinesq flow
a Generic representation of either a y or a z due to ambient turbulence (m)
0(6
-------
LIST OF ABBREVIATIONS AND SYMBOLS (continued)
o ya <5y due to ambient turbulence, accounting for terrain-influenced changes in turbulence
parameters (m)
0y« Effective ff y, accounting for effect of strain in the flow over terrain on horizontal diffusion
o x Vertical distance from plume centerline in a Gaussian plume at which the concentration falls
to e''/2 of the value at the plume centerline (m)
aia Value of az at s0 (m)
a * Growth in a ^ between s 0 and the receptor position of interest for flat terrain
a 'z Growth in a , between s 0 and the receptor position of interest with terrain influences con-
sidered (m)
a Ia a z due to ambient turbulence, accounting for terrain-influenced changes in turbulence
parameters (m)
a it Effective a », accounting for effect of strain in the flow over terrain on vertical diffusion
6 Potential temperature (K)
90 Surface potential temperature (at 2 = 0) (K)
9m Mean potential temperature of the mixed layer (K)
xv
-------
ACKNOWLEDGEMENTS
The authors wish to thank the following contributors to the development of CTDMPLUS: {Catherine
A. Brehme, Gary A. Briggs, William H. Snyder, Nicholas Z. Karasek, Ronald E. Stogner, and Kenneth J.
Galluppi. For their contributions to this expanded and improved user's manual, we also thank David J. Sul-
livan, Jeanne R. Eichinger, Christine D. Rafferty, Kathleen Rea, Francis Pooler, and Thomas Pierce.
Special thanks are extended to Francis A. Schiermeier, Peter L. Finkelstein, and all the members of
the EPA CTDM Technology Transfer Work Group for their steady support of the efforts to extend and
improve CTDM.
The authors of the original CTDM User's Guide wish to identify two important contributors to the
development of the CTDM code. During the formative stage of the effort, Akula Venkatram identified and
developed the two-layer approach that remains a central feature of CTDMPLUS. His guidance and insight
in the development of the theoretical basis of the model greatly contributed to its technical integrity.
Donald DiCristofaro produced much of the code for the earlier versions of CTDM, and his participation in
the many revisions to that code helped to maintain its clarity and efficiency.
xvi
-------
SECTION 1
INTRODUCTION
1.1 BACKGROUND
The Complex Terrain Model Development (CTMD) project was initiated by the U.S. Environmental
Protection Agency (U.S. EPA) to develop a practical, refined plume model for elevated point sources near
complex terrain. Accomplishments of this model development project are summarized by Strimaitis et al.
(1987). The major result of this effort was the Complex Terrain Dispersion Model (CTDM), which was
completed and released to the user community in early 1988.
Nighttime, stably-stratified conditions are generally associated with maximum short-term impacts
from elevated sources in complex terrain, as concluded by participants of an EPA workshop (Hovind et al.,
1979). CTDM was developed primarily to model these conditions.
However, since many air quality standards are written for averaging periods that generally include
daytime/unstable conditions, additional algorithms for estimating impacts from elevated sources in complex
terrain during unstable/convective conditions have been developed and integrated with CTDM to form
CTDMPLUS (CTDM PLus Algorithms for Unstable Situations). CTDMPLUS also includes several fea-
tures designed to improve user friendliness.
CTDMPLUS differs from currently available complex terrain models in the following areas:
In stable/neutral conditions, the structure of the two-layer flow (above/below the dividing-
streamline height, Hc ) is explicit in the formulation, and plume material that straddles the inter-
face remains in the respective layers (the plume is not treated as if it were all in one layer or the
other). Above Hc , the material is deflected and distorted, and the rate of dispersion is altered.
Below Hc , the stagnation streamline divides the flow, and only material that diffuses onto the
stagnation streamline is able to reach the surface of the hill. The stagnation streamline and the
concentration pattern wrap around the terrain. Plumes that lie to one side of the stagnation
streamline pass around the terrain.
The rate of plume growth depends on the turbulence and, in the case of o z, it also depends on
the degree of stratification. Sector averaging in the lateral direction is not used.
For plumes released into daytime convective layers, a probability density function (PDF)
approach is used to describe the vertical distribution of pollutants, and convective scaling con-
cepts are utilized to parameterize the lateral dispersion coefficient. CTDMPLUS considers the
effect of terrain on pollutant trajectories and on mixed-layer height deflections.
1-1
-------
Partial plume penetration into elevated stable layers is considered.
The three-dimensional nature of the terrain is used in the flow distortion calculations.
Vertical variations of windspeed, wind direction, and turbulence intensities are determined
through either interpolation of input measurements or surface layer scaling.
12 MODEL APPLICABILITY AND TECHNICAL LIMITATIONS
CTDMPLUS is a point-source Gaussian1 plume dispersion model designed to estimate hourly-
averaged concentrations of plume material at receptors on or near isolated terrain features. Although mul-
tiple sources and multiple terrain features are modeled by CTDMPLUS, the plumes impinging on a
downwind "hill" are assumed to be influenced by upwind features only in the effects found in the on-site
meteorological measurements.
The following restrictions and assumptions about CTDMPLUS should be understood:
CTDMPLUS contains no wake algorithms for simulating the mixing and recirculation found in
cavity zones in the lee of a hill. Therefore, sources within the lee of terrain features will not be
modeled with consideration given to wake effects. Estimates of concentrations at receptors in the
lee will not be reliable wiren such zones are present.
The path taken by a plume through an array of hills cannot be simulated. The model relies on
measurements of the flow obtained in the neighborhood of the source to define the incident flow
field for each of the terrain segments independently. If there is a strong channelling of the flow
due to large-scale terrain features (e.g., a valley setting), then this influence is reflected in the
modeling only insofar as it is to be contained in the meteorological measurements. Any changes
to the plume size caused by one hill are not carried forward to subsequent simulations downwind.
For situations when an individual hill is difficult to isolate from a complex terrain structure,
CTDMPLUS results may not be reliable.
Real terrain features are approximated by ideal shapes (determined by the terrain preprocessor).
CTDMPLUS considers elliptically-shaped horizontal contours and overall Gaussian shaped (ver-
tical profile) hills.
CTDMPLUS does not simulate calm meteorological periods. Therefore, under these conditions,
either no concentrations are estimated or (at the user's option) the following assumptions about
plume height variables are made in order to provide modeled concentrations:
1A skewed bi-Gaussian PDF approach is used in unstable conditions in describing the vertical pollutant
distribution.
1-2
-------
1. Minimum scalar wind speed = 1.0 m/s
2. Minimum a = 0.2 m/s
3. Minimum
-------
must also be established prior to contour digitization so that the locations specified for the stacks,
meteorological tower, hills, and receptors are consistent. A file containing ( x, y ) points along hill contours
must then be generated for each hill (see Table 1 of the terrain preprocessor guide for the format of this
file). The spacing between the points should be small enough so that the orientation and shape of each
contour is accurately represented. The contour digitization method is left to the user's discretion; a graphics
pad can be used or the map can be sent ta a vendor that specializes in digitization work.
The CTDMPLUS terrain preprocessor consists of two programs that process digitized contour data
to provide hill shape parameters in a format suitable for direct input to CTDMPLUS, and a third program
that displays the results (see the terrain preprocessor manual for complete description and user instruc-
tions). The first program, FITCON, evaluates and edits each contour and processes the data by numerical
integration to determine the following parameters for an equivalent ellipse: semi-major and semi-minor axis
lengths, contour centroid coordinates, and the orientation of the ellipse.
These parameters are input to the second terrain preprocessor program, HCRIT, which determines,
for the portion of the hill above several given elevations, the best-fit inverse-polynomial (vertical) profiles
along the hill major and minor axes. The center coordinates of the fitted hill are calculated as the mean of
the ellipse center coordinates for those contours above a given critical elevation. The orientation of the
fitted hill is calculated as a vector average of the ellipse orientations weighted by the ellipse eccentricity.
HCRIT provides an input file for CTDMPLUS that contains the following information for each critical ele-
vation:
Ellipse parameters corresponding to the contour at a critical elevation
Coordinates of the center of the fitted hill
Orientation of the major axis of the fitted hill with respect to north
The length scale and exponents for the inverse polynomial fits along the hill major and minor
axes.
CTDMPLUS uses the contour representations to provide hill shape information above the critical dividing
streamline height for each hour and each hill using interpolation between values specified at "critical" eleva-
tions.
The third program, PLOTCON, generates the following screen displays to aid in the evaluation of the
hill fitting process:
Map of the digitized input contours
Map of digitized contours that have been qualified and edited
1-4
-------
Map of the digitized contours and their associated fitted ellipses
For each critical elevation, a map showing the digitized contours and the contours of the fitted hill
at elevations corresponding to the elevations of these digitized contours that are above the critical
elevation.
13.2 Meteorological Preprocessor, METPRQ
The CTDMPLUS meteorological preprocessor, METPRO, uses on-site measurements of wind, tem-
perature, turbulence, and surface characteristics to provide predicted values of the nocturnal surface layer
height, h , or the daytime mixed-layer height, z-t and computed values of the surface friction velocity, u>, the
Monin-Obukhov length, L , and surface roughness length, z0 These values allow CTDMPLUS to estimate
profiles of wind, temperature, and turbulence above and below the range of measurements.
METPRO has four modes of operation:
Mode 0: Run for one or a few nighttime (stable) hours (surface characteristics such as surface
roughness, albedo, and Bowen ratio are assumed to be constant for all model hours); no
off-site or upper air data files are required.
Mode 1: Run for any number of stable hours that need not be contiguous, such as only nighttime
hours for several days; no off-site or upper air data are required, but site characteristics
are allowed to vary each hour.
Mode 2: Same as Mode 1, but off-site surface data are read to obtain cloud cover data.
Mode 3: Run for a series of contiguous hours that must come in blocks of complete calendar days,
although the days need not be contiguous; off-site and upper air files are required.
Input requirements are the least for Mode 0 and the most extensive for Mode 3. Mode 3 features the
determination of convective mixing heights, a computation that requires input for an entire day at a time.
The other modes do not compute this figure, but do provide all other variables.
The METPRO preprocessor output is an ASCII file that is read in by CTDMPLUS. The distinct,
modular nature of this method of processing input data for CTDMPLUS makes it possible for a user to
substitute an alternative meteorological preprocessor as long as the output format of "SURFACE" is the
same as if it were produced by METPRO. METPRO is described in detail in a separate user's guide
(Paine, 1987).
1-5
-------
133 Receptor Generator, RECGEN
The receptor generator program, RECGEN, is an interactive program written in Pascal that com-
putes receptor coordinates along selected contours. The program creates a file containing the receptor
information directly input to CTDMPLUS (file "RECEPTOR"). RECGEN also displays the contour lines
and receptor locations on the monitor screen. The user can easily add receptors, in addition to those pro-
duced by RECGEN, using a text editor.
RECGEN reads the output plot file from the terrain preprocessor FITCON. For each contour, the
user has the choice of (1) supplying a fixed receptor spacing along the contour, (2) supplying a fixed number
of receptors to be generated along the contour, or (3) having no receptors on the contour. RECGEN is
described in greater detail in Section 4.1.
13.4 Complex Terrain Dispersion Model Plus Algorithms For Unstable Situations, CTDMPLUS
The main CTDMPLUS program performs the plume transport and dispersion calculations for the
entire period of simulation. It takes the files prepared by the meteorological and terrain preprocessors,
together with the source, receptor, and (if necessary) rawinsonde and emission files, and executes under the
control of the options specified in an input file named "CTDM.IN". Modeled concentrations can be stored
in a binary or an ASCII output file, "CONC", if desired. The model lists all of the control data describing
the simulation, and pertinent source, receptor, and terrain data, to the file "CTDM.OUT". In case-study
mode, extensive tables of selected variables are also listed for computations performed for each source, hill,
and receptor. Technical aspects of CTDMPLUS are described in detail in Section 2 and Appendix A; user
instructions are described in Section 3.
1.3.5 Interactive Setup Program, SETUP
The interactive setup program allows the user to create or modify existing input files for METPRO
and CTDMPLUS. The program uses the explicit input file names required by METPRO and CTDMPLUS
in its search for existing files. Files that can be modified are "OPTIONS", "SURF1", and "SURF2" for input
to METPRO; "CTDM.IN", "SURFACE", and "TERRAIN" for input to CTDMPLUS; and "PROFILE",
which is input to both programs. Of course, all METPRO and CTDMPLUS input files can also be created
and edited by the user with any text editor that creates ASCII files.
The interactive setup program creates a batch stream to execute system commands such as renaming
and deleting files and executing the programs. Therefore, this program may not be directly portable
between operating systems. The interactive setup program is described in detail in Section 4.2.
1-6
-------
13.6 Concentration Postprocessors
There are three concentration postprocessor programs: CHIRET, CHIDIS, and CONTOUR. CHI-
RET processes the receptor and concentration information and creates a file used by CHIDIS and CON-
TOUR. CHIDIS and CONTOUR can be used to display predicted concentrations from CTDMPLUS on a
map of hill contours. Concentrations for a series of hours can be displayed sequentially, but each hill must
be done individually.
CHIDIS is a Pascal program that will draw the hill contours and provide a spatial display of concen-
trations for each hour. For a given hour, the receptor with the highest concentration is displayed as a blink-
ing, filled circle. The user can continue to display receptor locations in descending order of predicted
concentration. The current receptor blinks and all the previous receptors are displayed as non-blinking
filled circles. Information about the current receptor can be presented at the top of the plot if desired.
CONTOUR is a Pascal program that grids the concentration and receptor data and draws concentra-
tion isopleths on a background of hill contours. CONTOUR is a part of the menu driver and is not a stand-
alone program.
CHIRET and CHIDIS are described in Section 4.3; CONTOUR is described in Section 5.2.8.
13.7 Menu Driver Program, DRIVEIT
DRFVEIT is a menu driver program, written in Pascal, that displays lists of available files, prompts
the user to choose the necessary files to run the various programs in CTDMPLUS, and executes the desired
program. The user is led through the file selection process and responses are evaluated to detect potential
errors which can be viewed in the diagnostic output Gle. The menu driver requires that file names have
certain extensions and be located in specific directories. The menu driver is described hi detail in Section
4.4.
1.4 CTDMPLUS COMPUTER CONSIDERATIONS
1.4.1 Software Considerations
The complete CTDMPLUS package is designed to be run on an IBMR-compatible personal com-
puter (PC) system. However, all the essential portions can be satisfactorily run on other computer systems.
Most of the programs are written in ANSI standard FORTRAN 77 and should therefore be compiler
independent. The PC-based executables were created using the Ryan-McFarland RM/FORTRAN com-
piler. Three programs, PLOTCON, RECGEN, and CHIDIS, are written in Pascal and were compiled with
Borland International Turbo Pascal^ Finally, a text editor capable of creating ASCII input files is needed.
1-7
-------
1.4.2 Hardware Limitations
With the exception of a math co-processor chip, the programs that are included in the modeling sys-
tem require no special hardware in order to run. Attempts to run the modeling system without a math co-
processor will cause a runtime error number of 4001. The graphics programs use a procedure that
automatically detects the graphic system being used and scales the screen output accordingly. The largest
program, CTDMPLUS, requires a core space of about 360K bytes. When this program is run with the
menu driver, a core space of 640K bytes is required.
The terrain preprocessor requires input of digitized terrain information that can be obtained either
from an outside vendor or by using a graphics board.
The interactive setup program (SETUP), which is a FORTRAN 77 program, is designed for a PC
MS-DOSR operating system, in that file manipulation commands are based on MS-DOS". With some minor
modifications, SETUP could be made to run on other computer systems.
1-8
-------
SECTION 2
TECHNICAL DESCRIPTION
This section contains a description of the basic theoretical components of CTDMPLUS, and provides
the user with an overview of the subroutine structure in the code. An overview of the model is first pres-
ented in subsection 2.1. Subsection 2.2 discusses the calculation of meteorological variables at plume height.
Subsection 2.3 describes the methods of computing plume rise. Subsections 2.4 through 2.7 refer to the
modeling approach for stable/neutral conditions: 2.4 contains information on dispersion coefficients; 2.5
describes LIFT, the part of CTDMPLUS that computes concentrations at receptors in that part of flow that
is neutral or weakly stratified; 2.6 describes WRAP, the portion of CTDMPLUS that computes concentra-
tions at receptors in the strongly stratified part of the flow, and a brief note about receptors which are not
located on hills is contained in Subsection 2.7. Subsections 2.8 through 2.11 refer to model components for
computations in unstable conditions. Subsection 2.8 discusses plume penetration into elevated stable layers;
2.9 describes the probability density function (PDF) model of the vertical pollutant distribution and adjust-
ments for terrain; 2.10 describes the terrain adjustments to the daytime mixed layer height and their rela-
tionship to the PDF model; and 2.11 describes the horizontal distribution function and lateral dispersion
coefficient.
2.1 OVERVIEW OF THE COMPLEX TERRAIN DISPERSION MODEL PLUS ALGORITHMS
FOR UNSTABLE SITUATIONS (CTDMPLUS)
CTDMPLUS is a point-source steady-state plume dispersion model designed to estimate hourly-
averaged concentrations of plume material at receptors near an isolated hill or near a well-defined segment
of an array of hills. Emphasis is given to simulating situations in which the flow is toward the terrain, and in
which maximum concentrations are expected on either the windward side of the impacted terrain or near
the terrain crest. The original model (developed under contract with Environmental Research and Technol-
ogy (ERT) Inc. and named CTDM) was designed primarily to model the important situation of stable plume
impaction on the windward side of terrain. The model has since been expanded with compatible algorithms
to model all stability conditions, thus permitting the computation of concentration impacts over continuous
diurnal periods. The new model components were developed and incorporated without alterations to the
original CTDM methods. Therefore, CTDMPLUS will produce concentrations identical to those from
CTDM for releases into stable or neutral layers.
2.1.1 Pollutant Releases into Stable/Neutral Layers (An Overview)
In stable/neutral conditions, CTDMPLUS retains the Gaussian formulations of CTDM. The Gaus-
sian plume model for simulating the dispersion of pollutants from a continuous point-source describes a
2-1
-------
plume by its average properties as a function of distance along the flow downwind of the point of release.
The concentration of material in the plume is prescribed by a Gaussian distribution in a plane perpendicular
to the flow. The vertical distribution in this plane has a length scale denoted as a z , and the lateral distribu-
tion has a scale denoted as ay . Complete reflection of the plume at the ground assures that no plume
material disappears from the atmosphere. The concentration of plume material at any point downwind of
the source is determined by the size of the plume (ay and az ), the wind speed, the strength of the source,
and the distance of the sampling point from the axis of the plume. For example, if the plume is narrow in
the vertical (a, is substantially less than the height of the axis of the plume), the peak concentration is
found right at the center of the plume. Sampling points away from the center detect lower concentrations.
Over level ground, peak ground-level concentrations are found when a, is of the same order as the height
of the plume.
When a hill is present, the path of the plume changes as it flows over or around the hill causing a shift
in the relative position of a receptor and the center of the plume. But there are also mechanisms for chang-
ing the rate at which the material diffuses toward the surface, and for allowing the center of the plume to
impinge on the surface of the hill. These mechanisms are generally responsible for increasing peak
concentrations expected over terrain beyond those concentrations that would have been expected for the
same meteorological conditions over flat terrain.
In the absence of stratification, all streamlines in the flow pass over a hill. The centerline of a plume
in this flow follows the streamline that passes through the source of that plume. As the plume grows in the
vertical and horizontal directions (in the plane perpendicular to the flow), plume material diffuses across
adjacent streamlines, eventually reaching the set of streamlines that marks the surface of the terrain. Dis-
tortions in the flow which are induced by the hill change the position and relative spacing of the streamlines
from their initial distribution, and therefore change the shape of the plume as it passes over the hill. As a
result, ground-level concentrations (GLCs) of plume material change in two ways: The first and most
obvious change in the GLCs is the shift in the distribution on the surface of the hill, arising from the change
in the shape of the plume. Typically, the plume stretches in the horizontal as it passes over the crest of a
simple three-dimensional hill and this stretching produces a wider "footprint" over the hill. The second
change in the GLCs is the change in magnitude, arising from the effect of the distortion on the rate of diffu-
sion of plume material across streamlines. Typically, spacing between streamlines is reduced in the vertical
and expanded in the horizontal, while the speed of the flow increases over the crest. These changes tend to
increase the diffusion in the vertical and reduce it somewhat in the horizontal, thereby altering the magni-
tude of the GLCs over the hill. If the diffusion were not altered by the distortion in the flow, the peak con-
centration would not change from the flat terrain value.
2-2
-------
The nature of the flow changes dramatically when the flow is very stably stratified. A two-layer struc-
ture develops in which the flow in the lower layer primarily deflects around the hill, while the flow in the
upper layer travels over the top of the hill. A critical height H c defines the boundary of these two layers in
CTDMPLUS. This concept was suggested by theoretical arguments of Drazin (1961) and Sheppard (1956)
and was demonstrated through laboratory experiments by Riley et al. (1976), Brighton (1978), Hunt and
Snyder (1980), Snyder et al. (1980), and Snyder and Hunt (1984). In the layer above Hc, the approach flow
has sufficient kinetic energy to transport a fluid parcel up and over the hill against the density gradient of the
ambient stratification. In the layer below H c, the approach flow has insufficient kinetic energy to push the
parcel over the hill, so the flow below H c is restricted to lie in a nearly horizontal plane, allowing little
motion in the vertical. Consequently, plume material below H c travels along and around the terrain rather
than over it.
Above H c, the flow is similar to that just described above although the degree of distortion depends
on the stratification. Below H c, the flow is approximated as an ideal, steady, two-dimensional flow. Within
this flow, only one streamline at each elevation touches and follows the surface of a hill, and is referred to as
the stagnation streamline. Plume material reaches the surface of the hill only if it reaches the stagnation
streamline. If the plume centerline lies along the stagnation streamline and if it also lies below H c, the
center of the plume impinges on the hill. But if it lies to one side of the stagnation streamline, the centerline
will pass to one side of the hill.
The position of Hc and the stagnation streamline relative to the centerline of the plume dominates
the degree to which a hill in stratified flow is able to alter the peak ground-level concentration obtained in
the absence of the hill. Figure 2-1 illustrates this. In the model, the H c surface slices the plume into two
pieces as the hill is encountered. Plume material now residing below H c is sliced once again by the stagna-
tion streamline. Concentrations on the surface of the hill above Hc are determined by the cut made by the
H c surface because this now coincides with the bottom of the plume, which is in contact with the surface.
Concentrations on the surface of the hill below Hc are determined by the cut made by the plane of the
stagnation streamline because this cut coincides with the sides of the plume segments that are in contact
with the surface of the hill. As illustrated, receptors on the hill "see" concentrations that are much nearer
the center of the plume than do receptors in the absence of the hill. Figures 2-2 and 2-3 provide further
insight into how the plume is modeled by CTDMPLUS in stable to neutral layers.
Figure 2-2 addresses plume material above //c. The upper portion illustrates what the plume may
actually look like hi vertical cross-section as it travels along the surface. Material below Hc is removed at
s0, and the remaining material is distorted in the flow and reflected from the surface of the hill. The size
of the plume in the vertical at s depends on the amount of distortion in the shape of the plume as well as the
2-3
-------
PLUME CROSS-SECTION
LOOKING DOWNWIND
C constant
PLUME CENTERLINE
He
I
STAGNATION STREAMLINE I
I
/
PEAK GLC
(FLAT)
SURFACE UPWIND OF HILL
Figure 2-1. Idealized picture of how the dividing-streamline Hc plane and the plane of stagnation streamlines
"cut" into a plume, allowing material nearer the center of a plume to contact the surface of a hill.
2-4
-------
Hr
'ze
Figure 2-2. Depiction of plume behavior as it passes over a hill as seen from the side: upper figure shows
actual plume path, bottom figure shows treatment within CTDMPLUS with all ground-level recep-
tors placed upon a flat surface of height Hc above stack base and the plume with an effective
2-5
-------
amount of additional growth of the plume caused by changes to the rate of diffusion. The lower portion of
Figure 2-2 illustrates how the model actually treats the plume. Reflection of plume material from the sur-
face ( z - 0 ) is allowed from the source to s0, and reflection of plume material above H c is allowed from
the surface z = Hc beyond s0. Furthermore, the distortion b the flow (and the plume) beyond s0 is
scaled out, leaving only the effect of the distortion on the diffusivity in what is termed the effective sigma-z
(o ). As illustrated, a,, exceeds a x, the plume size in the absence of the hill, because the diffusion in the
vertical across streamlines is increased by the contraction in the vertical spacing of the streamlines.
Figure 2-3 addresses plume material below H c. The diagram on the left illustrates the plume in
horizontal cross-section as it splits and flows around a hill. In this case, some plume material has crossed
the stagnation streamline before the hill is encountered, so plume material is found on both sides of the hill.
Once the plume wraps around the leading edge of the hill, the stagnation streamline (which forms the
boundary of the hill) becomes a reflecting surface in addition to the plane z = 0, and material cannot dif-
fuse from the segment of the plume on one side of the hill to the segment on the opposite side. The dia-
gram on the right side of Figure 2-3 illustrates how the model treats this flow. The actual surface of the hill
is replaced by a line in the plume which corresponds to the stagnation streamline cut in Figure 2-1. All
distortion is scaled out as was done for the flow above Hc, and the effect of the distortion on lateral diffu-
sion is ignored as a second order effect. The lateral position of all receptors below Hc essentially collapses
onto the stagnation streamline, as depicted by the points labeled A and B, so that each is the same distance
from the plume centerline. However, concentrations along one side of the line differ from those on the
other side because diffusion through the line is not allowed.
Adjusting receptor positions while keeping the trajectory of the plume a straight line simplifies the
mathematics of CTDMPLUS a great deal. Rather than keeping track of the actual boundary of a hill and
the deformed trajectory of each of the segments of the plume as in the left portion of Figure 2-3, concentra-
tions are computed at receptor points A and B for a plume geometry like that of the right portion of Figure
2-3, which is only slightly more complicated than that for flat terrain. A similar adjustment of receptor
positions is employed for receptors above Hc, as illustrated in Figure 2-4. The upper portion of the figure
shows how a plume in the flow above Hc distorts over a hill as viewed from above. Three streamlines are
marked, the centerline of the plume, and streamlines passing through receptors A and B. When the deflec-
tion of each streamline is removed, and the distortion in the plume is scaled out, an equivalent plume-
receptor geometry is obtained, as illustrated in the lower portion of the figure.
2.12 Pollutant Releases into Unstable Layers (An Overview)
The methodologies employed in CTDMPLUS for modeling releases into unstable layers (where the
Monin-Obukhov length is negative) are designed to take advantage of our recent understanding of turbu-
2-6
-------
Impenetrable Boundary
Stagnation
Streamline
Plume Centerline
Figure 2-3. Depiction of plume behavior as it is deflected around a hill as seen from above: the stagnation
streamline in the left figure separates flow going around the right and left sides of the hill. In
CTDMPLUS, the hill is treated as being collapsed into an impenetrable wall (right figure) that sep-
arates the flow and plume material going on either side of the hill, with an effective crosswind dis-
tance and av.
2-7
-------
X PLUME CENTERLINE
WITH DEFLECTION * DISTORTION
PLUME CENTERLINE
B
DEFLECTION * DISTORTION REMOVED
Figure 2-4. Depiction of plume behavior in CTDMPLUS as it passes over a hill as seen from above: upper
figure shows actual deflection and distortion of the plume, bottom figure shows treatment within
CTDMPLUS with source-receptor spacing equal to spacing between streamlines upwind of hill and
appropriate adjustment of effective oy.
2-8
-------
lence and diffusion in convective boundary layers as outlined in the literature over the past two decades. At
the same time, compatibility with the original CTDM is maintained in order to utilize that existing
technology and to minimize additional user input requirements.
The CTDMPLUS algorithms for unstable conditions focus on the following areas: 1) plume rise in
unstable conditions; 2) buoyant plume penetration into an elevated stable layer; 3) application of a probabil-
ity density function (PDF) approach to define the vertical pollutant distribution with use of the original
CTDM flow module to calculate flow distortion by terrain; 4) terrain adjustments to the mixed-layer depth
(based on recent fluid modeling studies) and associated effects on the flow in the mixed layer; and 5) lateral
dispersion parameters based on convective scaling concepts.
Unstable plume rise in CTDMPLUS is calculated with methods suggested by Briggs (1975) with
modifications to account for the relationship of the risen plume to the mixed layer lid and for transitional
plume rise. CTDMPLUS uses interpolated or extrapolated vertical profiles of the meteorology to estimate
plume rise in an iterative fashion.
An important consideration in determining the amount of the pollutant that is available to impact the
surface is plume penetration into the stable layer above the mixed layer. For a steady-state condition (no
significant mixed-layer growth during the modeling period) the porjfon of the plume that penetrates aloft
can be considered unlikely to affect ground level concentrations. CTDMPLUS uses an algorithm proposed
and tested by Briggs (1984) for estimating the fraction of a buoyant plume that escapes the mixed layer.
Plume penetration in this algorithm is dependent on the strength of the stratification aloft, wind speed, the
mixed layer height, and the buoyancy flux of the plume. CTDMPLUS uses the estimated fraction of pene-
trated plume material to adjust the mass in that portion of the plume that remains in the mixed layer and
impacts on the receptors.
With the exception of cloudy or windy conditions or with very moist surface conditions, daytime
boundary layers are dominated by convection that is driven by surface heating. The strength of the convec-
tion and the depth of the mixed layer are primarily functions of sun angle and surface conditions. In con-
trast to that in stable or neutral layers, the vertical distribution of pollutants released from elevated sources
into a convectively-driven mixed layer is not Gaussian. Distinctly non-Gaussian vertical distributions have
been repeatedly observed in convective tank studies (Willis and Deardorff, 1976,1978,1981), numerical
modeling studies (Lamb, 1979,1981), and field studies (Briggs and McDonald, 1978, Briggs, 1985,1986).
Pollutant distributions reflect the distributions of vertical velocities which, due to the nature of con-
vection, are not symmetric for updrafts and downdrafts. Although the average vertical velocity may be zero,
updraft velocities tend to be "stronger" and downdrafts cover more horizontal area. This area dominance by
downdrafts results in a negative value for the mode of the velocity distribution. Consequently, the vertical
2-9
-------
pollutant distribution from an elevated source exhibits a descending mode of highest concentration (the
region of highest concentration tends toward the surface) as the plume moves downwind (as observed in the
aforementioned studies). This contrasts with the mode (centerline) of a Gaussian distribution which
remains elevated.
To model the vertical distribution, CTDMPLUS employs a probability density function (PDF)
approach which is based on a PDF suggested by Li and Briggs (1988). PDF models assume that, with con-
vective turbulence, the Lagrangian tune scale is sufficiently large and that particles released into the flow will
remain approximately on their initial trajectory until the surface or the mixed layer top are reached.
Therefore, particles initially released into a downdraft will remain in a downdraft and those released into an
updraft will remain in an updraft until a boundary is reached.
Additionally it is assumed that the initial particle trajectories are linear and have a slope of w/Lf ( w
is vertical velocity and U is mean wind speed). For steady state conditions ( U = constant), the trajectory
of the particle is controlled by the vertical velocity at the time and point of release. In other words, the
vertical distribution of pollutant material will be controlled by the probability density function of vertical
velocities at the source.
Li and Briggs have proposed an asymmetrical bi-Gaussian form of the PDF which produces pollutant
distributions very similar to those observed experimentally (Willis and Deardorff, 1981). CTDMPLUS uti-
lizes the Li and Briggs PDF; modifications to the particle trajectories for terrain influences have been added.
When terrain is present, the trajectories become curved as the flow is deflected by the terrain. The
flow field is computed by the flow module (already present in CTDM) and is used to construct the paths
Unking a receptor to the source position. An algorithm was developed which calculates the correct path by
an iterative process and determines the initial vertical velocity that will send the particle down that curved
path. This, in conjunction with the PDF, determines the probability of particles reaching the receptor and
thus the associated concentration.
Another important consideration in modeling unstable conditions is the effect of terrain on the mixed
layer depth, z f. This variable has been established as the dominant length scale for parameterizing mixed
layer variables and serves as the height that separates the layer dominated by surface effects from that which
is influenced little by the surface. Consequently, pollutant that finds itself above z f will not generally
impact surface receptors. However, the generally stable stratification above z, has a strong influence on
the adjustment of the mixed layer depth as the layer flows over terrain. This in turn has an effect on the
terrain distortion of that portion of the plume that resides in this restricted mixed layer.
2-10
-------
CTDMPLUS accounts for the effects of the elevated stable layer on restricting the lifting of the
mixed layer (and thus the plume) over the terrain. Based on results of a recent fluid modeling study (Perry
and Snyder, 1989), CTDMPLUS estimates the deflections of z t over the terrain feature in question as a
function of the hill Froude number (itself a function of flow speed, hill height, and strength of the elevated
stratified layer). It then uses this estimate to calibrate the flow model so that it yields equivalent deflections
for the z j streamline. By using a bulk calibration such as this, deflections of the pollutant trajectories cal-
culated in the PDF approach reflect the restrictions on the mixed layer created by the stable layer aloft.
Finally, there is the question of the distribution of material in the horizontal. Since significantly non-
Gaussian plume shapes have not been observed in the daytime mixed layer, analytical diffusion models have
retained the usual Gaussian form (Briggs, 1985). CTDMPLUS also adopts a Gaussian horizontal
distribution function which depends on receptor crosswind distance and a lateral dispersion coefficient.
The lateral (crosswind) distance between the streamlines passing through the source and receptor is
determined with the flow model. Because of the deflection and distortion of the flow around and over the
terrain, the effective crosswind distance can be significantly different than the geometric lateral distance
between source and receptor.
The form of the lateral dispersion coefficient, ay , is based on convective scaling concepts and analy-
ses of various data bases. Hanna (1986) correctly points out that there are about as many ay formulations
as there are data sets to analyze. The one chosen for CTDMPLUS was suggested by Briggs (1985) and
Hanna (1986) after analyses of power plant data For passive plumes, ay is dependent on the height of the
mixed layer and a non-dimensional downwind distance. For highly buoyant plumes, ay is dependent on
those two variables plus the buoyancy flux of the release.
More detailed discussions of the algorithms for unstable conditions are contained in subsections 2.8
through 2.11. Users of CTDMPLUS who are familiar with the original CTDM users manual will find that
the technical discussions of the algorithms for stable and neutral conditions have changed very little. A pri-
mary consideration in the development of the unstable algorithms was to be compatible with and compli-
mentary to the original CTDM. Therefore, users will find that concentrations estimated with CTDMPLUS
will be identical to those of CTDM for those cases where CTDM is applicable. For this reason, results of
existing evaluations of CTDM with stable and neutral hours are valid also for CTDMPLUS.
22 CALCULATION OF METEOROLOGICAL VARIABLES
CTDMPLUS reads two files of hourly meteorological data:
2-11
-------
SURFACE, a file passed from the meteorological preprocessor (METPRO) that contains
derived surface boundary layer parameters such as the Monin-Obukhov length ( L ), friction
velocity (u. ), surface roughness length (z0 ), and mixed layer height (h or z, ).
PROFILE, a user-created file of on-site measurements of wind direction, wind speed, tempera-
ture, oe (or au ), and aa taken at one or more (no more than 50) heights.
In addition, for the unstable hours, CTDMPLUS requires the daily RAWIN file (already required by
METPRO) in order to estimate the temperature profile above the mixed layer.
CTDMPLUS computes values of wind direction, wind speed, vertical potential temperature gradient
(VPTG), ou, and 0 or | L \ > 100 m), but is provided for unstable surface layers only for execu-
tion mode 3 of METPRO.
If both observed and calculated mixed layer height values are available, the CTDMPLUS user has an
option (via a model parameter switch) to use one or the other for each hour. If only one value is available,
that one is used in spite of the user's preference. If both values are missing, in stable/neutral conditions, an
"unlimited" value (99999 m) is used. Because of the importance of mixing height for scaling in convective
conditions, CTDMPLUS will not provide concentration estimates for unstable conditions when both z, val-
ues are missing.
As documented in the user's guide for METPRO (Paine, 1987), the following formula is used to cal-
culate mixing heights, h , in stable conditions (Nieuwstadt, 1981):
2-12
-------
0.3 u.
where u. = friction velocity
L - Monin-Obukov length
and / = Coriolis parameter
For unstable conditions, the Carson (1973) prognostic formulation is used, as modified by Weil and
Brower (1983) to estimate the daytime mixed layer height, z,. This formulation is also described in the
METPRO manual. Note the distinction in notation between stable mixed layer height, h , and convective
(daytime) mixed layer height, zi.
222 Wind Speed
Both scalar and vector wind speed observations can be read by CTDMPLUS. If the vector wind
speed, uu , is not available, it is calculated from the scalar wind speed, u , and o, (or av ) after Yamartino
(1984):
u. = u(l-2)'/2 (2)
where e - sino,(l -Goe) andG - 0.073864
If ou is provided instead of oe, then oe is approximated by au/u. for use in equation 2, where u is
the scalar wind speed.
The CTDMPLUS user may elect to set measured wind speeds that are less than 1.0 m/sec to 1.0
m/sec, as recommended in EPA's Guideline on Air Quality Models, Revised (1986). If the user selects this
option, proportional corrections are made to the vector wind speed (u0) and the standard deviations of the
crosswind and vertical wind speeds, ou and aa , in order to preserve the turbulence intensity values: iy =
ou/u.u and iz =OB/UU.
The assignment of the wind speed (either vector or scalar) at plume height is done by either:
interpolating between observations above and below the plume height, or
extrapolating from the nearest measurement height to plume height.
2-13
-------
When extrapolation is necessary, CTDMPLUS checks the position of plume height relative to the
surface layer. For plumes within the surface layer, the model extrapolates by using surface layer profiling
equations to be described below. For plumes above the surface layer, CTDMPLUS extrapolates the wind
speed to the top of the layer using profile equations: this speed is used aloft if the uppermost measurement
height lies within the surface layer. If the uppermost measurement height is outside the layer, then that
measurement is used for all heights above the measurement height.
The profiling equations that are used for wind speed are listed below (Lumley and Panofsky, 1964;
Businger, 1973).
where, for unstable conditions, the stability correction denoted as V m is as follows:
\
and
For stable conditions, the stability correction V m is given by:
₯m = -4.7- (6)
If the scalar wind speed is missing at ah1 observation heights, CTDMPLUS does not predict concen-
trations for that hour.
223 Wind Direction
The assignment of wind direction at plume height is handled similarly to that of wind speed. The
preferable method involves interpolation between observations both above and below the plume height. If
extrapolation beyond the uppermost measurement height to plume height is necessary, then the method
used depends upon the relationship of the plume height to the depth of the surface layer. If the plume is
within the surface layer (as is the uppermost measurement height), then the wind direction is scaled with
height if the user elects to have this calculation performed (controllable via a model parameter switch). No
2-14
-------
attempt is made to calculate the turning of the wind with height above the surface layer, as there is no defi-
nite trend that can be relied upon (Hanna et al., 1986). Therefore, the wind direction used for plumes above
the surface layer is either:
the direction calculated at the top of the layer from scaling (if the uppermost measurement height
is within the surface layer), or
the direction measured at the uppermost sensor.
The Tennekes (1981) wind direction adjustment is used for the computation of the change in direc-
tion within the surface layer:
Aa = arcfan (7)
0.4 u
where u g is the geostrophic wind speed (here assumed to be the wind speed at the top of the surface
layer) and B is obtained as a function of h/ L using graphs from Melgarjo and Deardorff (1974).
If the wind direction is missing at all heights, CTDMPLUS does not attempt to predict concentrations
for that hour. "
22.4 Potential Temperature Gradient
CTDMPLUS calculates the vertical potential temperature gradient, dQ/dz (or VPTG), from tem-
perature data contained in the "PROFILE" file. Temperature at various heights should be obtained from
measurements of temperature differences so that the temperature gradients are calculated accurately. The
height assigned to a calculated dQ/dz value is the midpoint of the interval bounded by the two tempera-
ture measurements. Linear interpolation of dQ/dz values is used for heights between these midpoint
heights.
If the plume is above the uppermost height at which a value of dQ/dz is available, the last available
value is used. No attempt is made to scale dQ/dz with height, even within the surface layer (such as by
the use of a method like Stull, 1983). The scaling methods were found during the development of CTDM to
give worse performance (values of dQ/dz that were too low) than using the uppermost measurement, so
the latter method is used.
If observed temperature gradients are not available, a default value is chosen:
a stability class-dependent default value if the plume is within the surface layer, or
2-15
-------
the greater of the stability class default value and the value for stability E (slightly stable) if the
plume is above the surface layer.
The default values are 0.0 K/m for stability classes A-D, 0.02 K/m for stability class E, and 0.035
K/m for stability class F. These default values are consistent with the complex terrain screening models.
Stability class is determined by the method of Colder (1972) using a . and L .
225 Turbulence Intensities
Values of ou (or ae ) and
-------
23 PLUME RISE CALCULATIONS
2J.1 Momentum Rise
Momentum rise is used only if the buoyancy flux is zero (stack temperature not greater than ambi-
ent). The following formulas are used (Briggs, 1975) for momentum rise:
dw,
A/i = 3 - - (lOa)
2 -i2-r N l/3
""
A/I- 1.5 [:^^ I s-° (lOb)
t / su,
where d = stack diameter
w s = stack gas exit velocity
u, = stack top wind speed
Ta = ambient temperature
T = stack gas temperature
s=(g/7a)(de/dz).
If L > 0 (stable/neutral), the minimum plume rise from Equations lOa and lOb is used; if L < 0 and
L > - 100 (unstable/neutral) Equation lOa is used.
232 Neutral/Unstable Buoyant Final Rise
These formulas from Briggs (1975) apply for plumes within the mixed layer:
Final or transitional (x f = x ) rise:
A/i = 1.6(Fx*)"3/ u (lla)
where
x/=119F0'4 for F > 55mV3 and (lib)
x, = 49F0625 for F < 55mV3 (lie)
Unstable breakup rise:
A/i = 4.3 (F/u)3/s H'2/s (12)
2-17
-------
Touchdown plume rise:
2
A/i- l.O-^ll +|£| (13)
Neutral breakup rise:
M2/3
A/i - 1.3 -4 1 + -j- (14)
uul \ Ahy/
where u = wind speed at height h, + Ah/2,
F = buoyancy flux, 0.25^^^(7, -7"0)/T,,
// = -u?(0.4Z) is the surface heat flux, and
u?tf-OAw..
Transitional plume rise is calculated using Equation Ha with x / replaced by actual downwind dis-
tance. The neutral/unstable plume rise is the minimum of the transitional rise and the final rises of Equa-
tions 11 through 14.
An iterative technique is used for calculating plume rise in any of the bouyancy equations (11 through
14), since the plume rise is assumed to be a function of the wind speed and temperature gradient at a height
halfway between the stack top and final plume height. For some plume rise formulas, this method does not
converge for certain profiles of wind speed and/or dQ/dz. Therefore, iteration is stopped after five tries
at convergence (defined by less than 1% change between successive iterations). After each iteration (and
the fifth one, if no convergence), the plume rise guess for the next iteration (or the final rise, if no conver-
gence) is the average of the plume rise estimate for the previous two iterations.
233 Neutral/Stable Buoyant Final Rise
There are several final plume rise formulas available for stable conditions, depending on whether
winds are nearly calm or not, and depending on whether conditions are close to neutrality. The final neu-
tral/stable rise that is used in CTDMPLUS is the minimum of those calculated by means of Equations 11,
14, and the following equations 15 and 16 Briggs (1975):
Neutral high wind rise:
Bent-over stable:
Ah = 1.54 [F/(uu*)]2/3 h)'3 (15)
Ah = 2.6(F/us)"3 (16a)
2-18
-------
Calm stable:
A// = 4F1/4s'3/8 (16b)
In stable conditions, the distance to final rise, x f , if given by
xf = 2.07us"1/2 (16c)
The iterative technique described in Section 2.3.2 is also used for neutral/ stable conditions.
2.4 DISPERSION PARAMETERS FOR STABLE/NEUTRAL CONDITIONS
2.4.1 Plume Spread Due to Turbulence
The formulation for a z is described in Venkatram et al. (1984). It is based on the form of Taylor's
(1921) theorem of diffusion for very short and very long times of travel:
cz =
-------
By introducing two new constants, Y and I", the length scale, I, can be expressed as
i,_L.rr-L-Y.r-^-11 (20,
where W is the Brunt-Vaisala frequency given by
and
(««>
.36
»
a » - 1.3
u.
and a and P are parameters in the surface similarity profiles (.74 and 4.7, respec-
tively).
The resulting expression for I has two distinct limiting forms. When z is very large, for non-zero
N , stratification dominates the scale of the mixing process and
When N is nearly zero and z is finite, the length scale for the mixing process is proportional to
height above the surface:
/ = Tz = / (23)
The quantities I, and / are introduced to distinguish the mixing lengths for the stable and neutral
limits.
Equation 22 states that the length scale for turbulent mixing in the stable limit is proportional to
aa/N, where Y2 is the constant of proportionality. This is consistent with the notion that a fluid element in
this limit must overcome a stable potential temperature gradient in order to be displaced vertically. Given
2-20
-------
that the velocity scale for vertical dispersion is o,,, the length scale that naturally follows is proportional to
OU/N. In the neutral limit, the size of the turbulent eddies is restricted by the height above the surface, so
that the mixing length should be proportional to z , where F is the constant of proportionality.
A simple weighting function is used to obtain a length scale between the neutral and stable iimils:
I - I + I (24)
III
1 ln ls
With this expression for / , T L is computed as
TL- (25)
o.
and a z is computed from Equation 19 as a function of the time-of-travel, including source effects
(see Section 2.4.2).
An equation similar to that used for a x is used to compute ay as a function of the time-of-travel
and the turbulence velocity scale for lateral fluctuations, a :
(26)
The departure of this expression from that for a z arises in specifying the functional form of T L' ,
which is the Lagrangian time-scale of the transverse correlogram. In this case, T L' cannot be derived from
measurements of the mean flow and its statistics. Instead, T L' is set equal to the time-of-travel required to
cover a distance of 10 km, which has the effect of reducing a y by about 18% from the value that would be
obtained using a linear growth law at a distance of 10 km.
All of the Complex Terrain Model Development field experiments included sampler locations within
10 km of the source, so that the growth in a y beyond 10 km was not documented. Nearly linear growth in
a y over the first 5 km was typical of many of the experiments. The use of the 10 km length scale is meant
to underscore the uncertainty of using linear growth beyond 10 km, while allowing nearly linear growth
within 5 km of a source. Note that the scale used by Briggs (1973) for dispersion in open country corre-
sponds to 5 km, which would reduce ay at 10 km by about 29%. Because CTDMPLUS is intended for use
in complex terrain, the use of a larger scale for TL' was judged appropriate in that lateral fluctuations in
the flow which are induced by complex topography would promote greater meandering of the plume.
2-21
-------
In summary, the following equations are used to calculate the plume spread parameters a, and a y:
1 + a,,,t\ ~rz
.540
°v = " T7^ (27b)
y [1 + ut/20,000]l/2
2.42 Allowance for Source-Induced Effects
A virtual time-of-travel is introduced to account for increases in a y and a, caused by how the
plume is released into the atmosphere. For buoyant releases in stable conditions, the initial growth during
plume rise is characteristically much greater than the growth caused by ambient turbulence alone. There-
fore, once the plume reaches its equilibrium height and ambient turbulence becomes dominant, the plume
evolves as if it had experienced a greater time-of-travel.
Write the general form for a z or a y as:
Ot
' ('*£)
and solve for t:
"2
(28)
1+1 + 4
T l-Y
TiUJ
The virtual time-of-travel ( t ) is found by setting a , equal to the plume size caused by source-
induced effects ( a i, ) and setting t equal to ts + tu , where
-------
Ambient turbulence is assumed to dominate source-induced turbulence quickly for neutral and unsta-
ble conditions, and t, is set to zero for these cases. As implemented, a,, is just the size of the plume
resulting from buoyant rise (buoyancy-enhanced dispersion, Pasquill, 1976):
AH
a" = a'b = 3~5 (31)
where A//is the plume rise. Similarly, t, is the time-of-travel to final rise. Depending on the rela-
tive values of a (ambient turbulence) and a ls for the time ts, tu can become negative in Equation 30.
Physically, this indicates that the growth of the plume due to ambient turbulence exceeds that due to source
effects. In the case of a y, this frequently happens when meandering is great. Consequently, tu is never
allowed to be less than zero, so that buoyancy enhancement is active only when it exceeds the growth rate
due to the ambient turbulence. In fact, 7" z is always greater than zero, as it has a minimum value corre-
sponding to the time it takes for the plume (a,, az) to grow to the size of the stack radius. Once t is
calculated for both the lateral and vertical scale of the plume, it is added to the actual time-of-travel in
Equation 27.
2.5 THE LIFT COMPONENT
2.5.1 Overview of LIFT
The flow above Hc is considered to be weakly stratified. That is, the stratification is strong enough
to influence the flow pattern (e.g., lee waves), but not strong enough to inhibit significant vertical motion.
To simplify the modeling task, Hc is assumed to be a level surface, and the flow above Hc is only affected
by that portion of the hill that lies above Hc.
Hc is computed for each hill from profiles of wind speed, u(z~), and temperature (Brunt-Vaisala
frequency, N ( z )) by locating the lowest height at which the kinetic energy (K.E.) of the approach flow just
balances the potential energy (P.E.) attained by elevating a fluid parcel from this height to the top of the hill,
// . The statement that defines that point of balance is:
iu2(//J = fV(z) (//-z) dz (32)
2 JHC
In practice, the value of H c for which the equality in Equation 32 is attained is found by rewriting
the integral as a series of sums over layers of constant N . For each layer, starting with the one that
includes the top of the hill, the bottom and the top of the layer are tested to see if the K.E. (left-hand side)
has decreased below the P.E. (right-hand side). The lowest layer in which the K.E. exceeds the P.E. at the
2-23
-------
top, but not at the bottom, is the layer that contains Hc. H c is then computed from Equation 32 by
assuming a linear velocity gradient within that layer, and constant N within each of the layers to the top of
the hill.
The plume is allowed to diffuse as if the terrain were perfectly flat until it reaches the point where its
trajectory crosses the height contour equal to Hc in elevation (say, at a distance s0 from the source, see
Figure 2-5). If Hc is zero, then this zone extends from the source to the base of the hill, although conceptu-
ally it could extend to any point where the hill is thought to exert a significant influence on the flow. Beyond
s0, the plume material below Hc is disregarded by the LIFT component, and the evolution of the
remaining material is modeled as if the terrain were flat, and the lower boundary were He (with full reflec-
tion). However, the rate of plume spread and the position of the plume centerline relative to the receptor
are modified to reflect the net alternation of these properties between s 0 and s (where s is the distance
from the source to the receptor) induced by the presence of the hill. The simplicity of the Gaussian plume
solution is retained in this way, while the full dilution of the plume from the source to the hill ( s 0) as well
as the effects of the hill on both flow and diffusion beyond s0 are explicitly incorporated.
Aside from the obvious distinction of incorporating the primary influence of Hc on the plume-
terrain interaction, this approach notably differs from the current regulatory modeling approach (at least as
embodied in Valley, COMPLEX I, and II) in that the terrain influence for a receptor on a hill only affects
the diffusion of the plume once it is over the terrain. The mathematical formulation of the "partial plume
height correction" approach of COMPLEX I and similar models actually "lowers" the plume at the source.
If this technique were engineered to produce the "correct" hill-influenced ground-level concentrations, the
terrain correction factor for a particular receptor would need to be a function of downwind distance, terrain
shape, and distance between source and terrain. As employed in regulatory modeling, however, the terrain
correction factor depends only on the stability class, so that its use has led to problems of interpreting "sur-
face reflection" from sloping terrain, as well as to problems in justifying values chosen for the terrain correc-
tion factor.
2.5.1.1 LIFT Structure
The terrain effect as modeled in LIFT includes re-initializing the flow at a distance s downwind of
the release. This re-initialization can be illustrated first for flat terrain and uniform flow. The concentration
at a receptor downwind of s 0 is composed of contributions from the entire concentration distribution at
2-24
-------
LIFT
WRAP
'Flat-Terrain'
Domain
LIFT Domain
WRAP Domain
'Flat-Terrain'
Domain
''mm
LIFT Domain
WRAP Domain
Figure 2-5. Idealized stratified flow about hills indicating domain of individual CTDMPLUS component algo-
rithms. The distance S0 is the distance from the source to the intersection of Hc with the hill
surface for receptors above Hc (the LIFT domain), and it is the distance from the source to the
terrain-height contour equal in height to the receptor elevation for receptors below Hc (the WRAP
domain). Note that in the upper figure the flow is into the page.
2-25
-------
Conceptually, the flux of plume material through the plane x - x , + s0 (note that the x-axis lies
along the flow direction, and the plume is released at x, , y , , zs) can be thought of as a distribution of
point sources. If we track the plume material in terms of the distance downwind of the source, s = x - x s ,
then the source strength of one of these point source elements located at the point (s0 , y , z) is given by:
dQ(s0,y,2) = C(s0,y,2)udy dz (33)
Because the flow beyond s 0 is considered to be uniform, the influence of each of these sources fol-
lows the Gaussian plume solution to the advective diffusion equation so that the contribution of the source
element at the point ( sa , y , z ) to the concentration at the point ( s , / , h ) is:
dcrs lh.s^-
dC(s,Z,/i,s0)-
where a 'y , a j denote the plume spread statistics for each point source element over the interval s - s 0
(see Figure 2-6). The total concentration at ( s . I , h ) is found by integrating Equation 34 over all point
source elements, so that
. +.
r, 7 K A f f
C(s,Z,/i;s0) =
J J
-°-5
2nava,
o _«, y 2
dydz (35)
The plume spread statistics o* and a'z for the interval s-sa are specified by the requirement that
Equation 35 for flat terrain reduces to the expression obtained for the original point source located at s - 0
(i.e., Equation 34 with s0 = 0 and a' = a(s)). Equating these two expressions for C, with h = 0, we obtain
So) a*-0*0 (36a)
5o)s<>y-ayo (36b)
Equations 35 and 36 illustrate the re-initialization technique for the limiting case of flat terrain and
uniform flow. Terrain influences are incorporated by altering the rate of diffusion within the interval
s- s0, and by changing the position of the receptor relative to the centerline of the deflected plume. Fur-
thermore, because no plume material below Hc travels over the hill and because Hc defines the lower
2-26
-------
u
C(s.z)
i )f I / t t I ) i I I
U
f ' ' ' f f f f ' ' ' ' ' '
)\ ] I ) i / f I } ) I 1 1 > t i t I I I I I
00
QS. z) = J C (s, e
where erz -
Figure 2-6. Illustration of the relationship between the crosswind-average concentration profiles at s^ and s,
and the plume from one of many point -source elements representing the flux of material across
the plane at sfl. The total concentration at a particular point C(y , 2) is constructed by summing
the contribution c(s , 2 ; s0 , z0) from each point-source element Q(s0, z0).
2-27
-------
boundary over the hill beyond s 0 (see Figure 2-7), the integration in the vertical hi Equation 35 is
performed over the domain z - H c to z-<*> and material from each point source element is reflected from
the boundary z - H c rather than z = 0.
Denote the receptor height above ground relative to the terrain-altered plume as h ' (see the lower
part of Figure 2-7). Denote its lateral position relative to the terrain-altered plume centerline as y '.
Denote the altered growth rates of the plume as a'z and a'y and define terrain-effect factors Tz and Ty
so that
-------
PHYSICAL PICTURE: material
hR above Hc rides up and over hill
in a distorted flow. Materiel
below Hc passes round the side.
LIFT CALCULATION: net effect of flow
distortion t» to increase the
effective rate of plume growth
(?n) over that in the
absence of the hill (?z)
Figure 2-7. Illustration of terrain effect on the vertical distribution of plume material above Hc as modeled in
LIFT.
2-29
-------
-o.5
« \ 2
°z \ 9
H -h''°»
fL
1 -F^F A.
)(39b)
and
(39c)
Equation 39 provides the framework for estimating concentrations due to plume material that travels
up and over a hill. It shows how the influence of the terrain affects the magnitude and the distribution of
GLCs. The most complicated part of Equation 39 is the expression for F x, the vertical distribution factor.
It contains four terms because it applies to an elevated receptor, so the image source contribution is not
equal to the contribution from the primary source (hence two terms are needed rather than one). And it
also applies to a plume segment above Hc rather than an entire plume profile, so that an image source
contribution at s0 must be explicitly maintained (hence, two more terms). The error functions also arise
from treating only the portion of plume material that lies above Hc at s0. If Hc is zero, then Equation 39
becomes:
-0.5
Qe
2nu.oyeoze
K'*,\z~\
) J40)
which, with the exception of alterations in the plume size and in the distance between the plume cent-
erline and the receptor, is the familiar Gaussian plume equation for the concentration at an elevated recep-
tor.
2-30
-------
The key quantities that need to be evaluated in order to apply Equation 39 are the effective receptor
coordinates (y x' and h,'), and the terrain-effect factors (T z and Ty). The factors Tz and Ty are
obtained from the results of Hunt and Mulhearn (1973), and require a model of the flow field over a hill as
does the receptor position relative to the deflected plume centerline. These features of LIFT are described
in Appendix A.
2.6 THE WRAP COMPONENT
2.6.1 Overview of WRAP
A particle in a steady two-dimensional flow around an obstacle will experience both accelerations and
decelerations as it passes by. The magnitude of these changes in speed depends upon how close the particle
is to the stagnation streamline of the flow. Maximum changes occur for particles on the stagnation stream-
line. Furthermore, the spacing between adjacent streamlines varies in inverse proportion to these changes
in the speed along streamlines. Figure 2-8 is a representation of a typical streamline pattern for flow around
an ellipse when the incident flow is at an angle to the axes of the ellipse.
A plume in this steady flow (with some small-scale turbulence) will follow the streamline patterns,
spreading slowly across adjacent streamlines. However, as streamlines spread apart (or contract) the plume
size in the horizontal will expand (or shrink) to the same extent. In the absence of diffusion, these kinematic
changes in the horizontal size of the plume will not alter the concentration of material within the plume.
Changes to the horizontal scale of the plume are balanced by changes in the flow speed so that the flux of
material is unchanged. With the addition of small-scale diffusion, the rate of plume growth in the horizontal
can be altered by changes in streamline spacing (Hunt and Mulhearn, 1973). However, based on the obser-
vations at CCB and Tracy field sites we choose to ignore the effects of small-scale diffusion on concentra-
tions in the WRAP component. The observations suggest that low frequency turbulence-meanders-
control crosswind plume growth over hourly averaging times.
2.6.2 WRAP Structure
To simulate ground-level concentrations due to dispersion of releases below H c in complex terrain
settings, CTDMPLUS must approximate the key features of steady two-dimensional flow around an ellipse
that were described above. Two key approximations in the WRAP component are (1) lateral diffusion is
insensitive to accelerations in the flow (i.e., the kinematic deformation of the plume has no effect on the
diffusion rate), and (2) the mean flow for the averaging period (one hour) is considered steady, while all of
the variability in the flow over the period, including that due to meandering, is considered "turbulence."
2-31
-------
Stagnation Streamline
Figure 2-8. Typical streamline patterns in two-dimensional flow around an elliptical cylinder.
2-32
-------
A primary difference between WRAP and LIFT formulations arises from the location of solid bound-
aries and the relationship between the position of these boundaries and the wind direction fluctuations. The
terrain effect is modeled in WRAP by re-initializing the flow at the distance s 0 downwind of the source
(see Figure 2-5). The concentration at a receptor downwind of s 0 is composed of concentrations from that
part of the concentration distribution at s 0 that lies below Hc, and that also lies on the same side of the
stagnation streamline as the receptor (see Figure 2-9). Reflection of plume material in the vertical is
allowed from the plane z - 0 over the entire distance s, and reflection in the horizontal is also allowed from
the hillside beyond s 0. Note that the stagnation streamline forms the boundary of the hill surface in hori-
zontal cross section.
For a receptor located on the hillside at a distance s (see Figure 2-9) and a height z,, above the
plane z - 0, the concentration due to one elemental point source located at (s0 , y , z) in the plume is given
by
,rf n , d<2(s0,y,z) -°
dC(s,0,zR\s0) = jr- 2e
2nuayaz
(41)
In Equation 41, we assume that the x-axis of the coordinate system points along the stagnation
streamline, and that the source is located at (x,, y,, z^. The total concentration at the receptor contains
contributions from those elements below H c and on the same side of the stagnation streamline as y s:
o (o)
where
2nuayoazo
(y-ys\
-0.5 -f
'
2-33
-------
Stagnation Streamline
Figure 2-9. Top view of a plume in two-dimensional flow around a hill. The shape of the hill in cross-section at
the receptor height is assumed to be invariant with height so that the deformation of the entire
plume around the hill is a function of the receptor location.
-------
These expressions are analogous to Equations 34 and 35 of the LIFT component. The integral for
dy has the limits (0) and (» ). This is meant to denote integration from 0 to + °° if the receptor lies on
the "positive" side of the stagnation streamline, and integration from - °° to 0 if the receptor lies on the
"negative" side. Note that material is allowed to diffuse upwards through H c so that receptors which lie
above Hc , and which are downwind of s0 , can receive a contribution from the elemental point-sources
below Hc at s0 . The total concentration at such receptors is the sum of both the LIFT and WRAP con-
tributions.
The integrals in Equation 42 are evaluated to obtain:
., S~\ - . *\ ^ «
ERF
(44)
Be * ^ J
Most of the notation here has already been encountered in Section 2.5. The factor 'sign (y e)' denotes
the sign of the receptor position in the coordinate system with x-axis aligned with the flow, and it results
from the choice of integrating over the "positive" or "negative" portion of the flow in Equation 42. The fac-
tors B [ and B 2 are given by
BI = ERF
Bz = ERF \ b' | * ERF [ ,' ] (45)
where
60 =
6, = Hca*
'L
The subscript R denotes the receptor location, and the subscript s denotes the source location (see
Figure 2-10). The distance from the stagnation streamline associated with the mean wind direction to the
centerline of the plume is denoted as d the total a y indicates the horizontal spread of the mean plume,
2-35
-------
Stagnation
Streamline
Figure 2-10. Sketch of the flow around an ideal cylinder of elliptical cross- section. The section shown is taken
either at the elevation of the centerhne of the plume (zs), or at Hc , depending on which is
smaller, and it indicates the relationship between the streamline through the source ( ty,), the
stagnation streamline (ip0 = 0 ), and the coordinate system with x-axis aligned with the mean
wind direction.
2-36
-------
and the total a x indicates the vertical spread. The amount of plume spread experienced over the distance
to the stagnation point is denoted as a20 and oyo , and the rate of plume growth beyond the stagnation
point to the receptor is denoted as a^ and a'y , where
o? - <*« - ol. and
a;2 = a2 - a20 (46)
The ellipse that is used to estimate the flow below Hc is taken from the horizontal cross-section of
the hill at the minimum of the following two elevations: either the elevation of the centerline of the plume,
or the elevation of Hc . In this way, the shape of the hill selected is associated with the peak concentration of
plume material found within the layer of fluid below Hc .
Concentrations at receptors located on the hillside just upwind of the stagnation point are estimated
as if the receptor sits on a pole of height equal to the receptor elevation above the base of the stack. Fur-
thermore, the lateral distance between the plume centerline and the receptor is set equal to d , so that con-
centrations at all of these receptors are controlled by the amount of material on the stagnation streamline.
In this way, plume material below Hc follows streamlines around the hill, and only material which diffuses
onto the stagnation streamline impinges on the hill. The equation for estimating these concentrations is:
c =
exp -0.5
(47)
Equations 47 and 44 provide estimates of ground-level concentrations of plume material before and
after plume material above Hc is "removed," respectively. That is, the upper and lower portions of the flow
do not become distinct until the impingement or stagnation point is reached.
2.7 RECEPTORS NOT TREATED BY LIFT OR WRAP
The LIFT and WRAP segments of CTDMPLUS compute concentrations at receptors that are
located on hills beyond the impingement point. Receptors that are upwind of the impingement point, or
those that are not associated with any of the hills in the modeling domain require separate treatment. This
treatment is described in the following two sections.
2-37
-------
2.7.1 Receptors Not Influenced by Hills
In theory (for stable/neutral conditions), the main subroutines of CTDMPLUS produce results which
are identical to flat-terrain results in the limit that the hill height goes to zero. However, the LIFT and
WRAP code is not designed to check for hill height of zero before executing statements that may require
division by the hill height; and if it were, many extensive calculations would be needlessly executed, with the
model returning no terrain-effects. Furthermore, the structure of these sections of the model requires that
receptors be associated with specific hills wherein each hill requires extensive information. A flat terrain
algorithm is included in the model to avoid such numerical problems and extra input requirements.
The flat terrain algorithm simply performs a Gaussian plume computation which assumes that there
is no mixing lid, that all receptors lie on a single ground-plane, and that plumes travel in straight lines. All
plume rise and growth algorithms match those used in the other sections of the model.
To invoke this algorithm, the user must select a hill number of zero for this receptor. Guidance on
which receptors should receive this treatment may be found in Section 3.
2.72 Receptors Upwind of the Impingement Point
All receptors on a hill that lie upwind of the impingement point sample a plume that is diffusing as if
all terrain were absent with respect to reflection from the lower boundary of the flow. That is, reflection is
allowed only from the surface z = 0 (the reference plane for the terrain features), and no changes are made
to the rates of diffusion. However, the sampling point is altered relative to the centerline of the plume.
Several possible receptor locations are sketched in Figure 2-11. The actual terrain surface can depart
significantly from the simple terrain shape used in the model, and this leads to the situation illustrated. Not
all of the receptors upwind of the impingement point may lie at elevations below H c. The surface of the
terrain beneath receptors 1 and 2 is below Hc, but because elevated receptors are allowed in the model,
receptor 2 actually lies above Hc. At receptor 3, the surface of the terrain exceeds Hc,so that this recep-
tor lies above Hc regardless of whether it is elevated or not. Each of these receptors is treated differently.
Receptor 1 lies within the two-dimensional flow below H c that is deflected around the terrain.
Therefore, the lateral position of the receptor is defined by the location of the stagnation streamline, as
described in the WRAP component. The vertical position of the receptor is given by z, in the figure.
Essentially, a standard WRAP-like computation is done for this receptor, except that plume material is
allowed to diffuse across the stagnation streamline, and plume material above Hc has not been removed.
(The removal of material above H c and the imposition of a reflecting surface along the stagnation stream-
line is initiated once the flow has reached the impingement point.)
2-38
-------
Figure 2-11. Treatment of receptors upwind of the impingement point. Receptors 1 and 2 sit above terrain
that lies below Hc . Both are modeled as "receptors-on-poles", preserving their actual height (Zj
and Zi) above the reference plane. Receptor 3 is above terrain that exceeds Hc in elevation, so
that the flow above Hc follows the surface of the hill. Therefore, the height of the "pole" Z3 is
less than the actual elevation of the receptor above the plane.
-------
Receptor 2 lies in the flow that travels over rather than around the hill. Because it lies upwind of the
impingement point, it is assumed to be upwind of the cut-off hill as well. The concentration is obtained at
receptor 2 by assuming that terrain-effects are totally absent. Hence, the flat-terrain algorithm is used in
combination with a "pole" height of z , as illustrated in the figure.
Receptor 3 is treated in nearly the same way. But because the terrain exceeds H c at this location,
the flow above H c is assumed to rise. Therefore, without accounting for terrain effects on dispersion or
streamline spacing, the streamline that passes through receptor 3 had originally been at an elevation of z 3
in the figure. The "pole" height, z3 , for this flat-terrain calculation is the sum of the actual receptor eleva-
tion above the surface of the terrain and H c .
2.8 PLUME PENETRATION INTO ELEVATED STABLE LAYERS
In modeling daytime unstable conditions, it is important to consider the interactions of elevated
plumes with the top of the mixing layer. In particular, buoyant plumes are capable of partially or wholly
exiting the mixed layer by penetrating into the stable layer aloft and thus becoming unavailable (except when
re-entrainment occurs) for impact on receptors within the mixed layer and particular those at the surface.
CTDMPLUS estimates the effects of plume penetration in the following way. Based on the integration of
the buoyancy flux equation for a bent over buoyant plume and the assumption that the vertical plume distri-
bution is rectangular ("top hat"), Briggs (1975) first defined a factor for estimating the fraction of the plume
that will penetrate into the stable layer above a well mixed layer. This formulation has received some use in
currently available models (in particular, the Rough Terrain Dispersion Model, Paine and Egan, 1987).
However, Briggs (1984) points out that his 1975 penetration factor overpredicts penetration data measured
by lidar at the Keystone power plant (Johnson and Uthe, 1969) and thus proposes a more conservative
approach where the stratification of the elevated layer is assumed constant throughout the entire layer of
rise. Then that portion of the plume above the mixing height, zit is assumed to penetrate and that below
does not. In this way, the penetration factor becomes:
where
S = (g/6)(de/dz) = the stabililty parameter,
dQ/dz = the potential temperature gradient above z\
Zb = Zj - stack height.
and Fb = plume bouyancy flux.
2-40
-------
For the majority of the Keystone data, equation 48 provides a superior fit in comparison to the Briggs (1975)
formulation and has been incorporated into CTDMPLUS for plume penetration calculations.
Note that the use of this formulation requires an estimate or measurement of vertical potential tem-
perature gradient in the stable layer aloft. A simple procedure for estimating this gradient from normally
available meteorological data will be discussed in section 2.10.1.
Once the penetration factor, P, has been determined, the remaining plume material is defined by
an adjustment of ( 1 - P ) to the original source strength and plume buoyancy flux For example, if 30% of
the plume penetrates into the stable layer then P - 0.3 and the source strength and buoyancy flux of the
remaining plume becomes 0.7Q and Q.7Fb respectively. These adjusted values are then used in any
subsequent calculations in the model.
2.9 THE VERTICAL DISTRIBUTION FUNCTION
CTDMPLUS calculates the vertical pollutant distribution separately from the horizontal distribution
(discussed in section 2.11) and then combines the two to obtain the receptor dependent concentration. Sec-
don 2.9.1 contains a discussion of the probability density function approach used to obtain the vertical dis-
tributions and section 2.9.2 contains a description of the adjustments made to the PDF for terrain induced
distortions in the flow.
2.9.1 The Probability Density Function (PDF)
The vertical distribution of pollutants (as discussed in Section 2.1.2) in a convective mixed layer has
been observed to be strongly non-Gaussian. These observations are clearly reflecting the skewed distribu-
tion of vertical velocities that are present in a convective mixed layer and which are responsible for the dis-
tribution of passive pollutants. One approach to modeling this skewed distribution is with a probability
density function (PDF) of the vertical velocities based on the observations available. Several PDF models
are described in the literature (e.g. Weil and Furth, 1981, Misra, 1982, and Venkatram, 1983, and a review
by Briggs, 1985).
CTDMPLUS employs a recently suggested (Li and Briggs, 1988) PDF approach which assumes (as
do other PDF models) that the the Lagrangian time scale is sufficiently large such that particles released
into the flow maintain their initial trajectories until a boundary is reached. With the assumption that these
particle trajectories are linear and have the slope w/U (w is vertical velocity and U is mean wind speed)
as the air flows past the release point, a steady wind, time-averaged, crosswind-integrated concentration is
given by
2-41
-------
CWIC-J"xdy-(Q/x)P1.(u»l|2f) (49)
where Q is the release rate (source strength), Pa is the probabih'ty of wz = U(z-zs)/x , zs is the
effective source height, z is the receptor height, and x is the downwind distance, w z is the convective
vertical velocity that will take the particle from source to receptor.
For the direct path between source and receptor, z, and z are the actual source and receptor
heights. For particles that reach the receptor through one or more reflections from the surface or mixing
height, z j and z must be adjusted based on source/receptor/surface/mixing lid geometry. Assuming per-
fect reflections from the surface and mixing lid, the contribution to the total CWIC due to trajectory reflec-
tions can be calculated by applying equations 49 with virtual sources that are located vertically from the real
source so that virtual trajectories intercept the receptor. Figure 2-12 illustrates the relationship between
"image" sources and the paths (including one or more reflections) from the source to the receptor. The total
CWIC as a result of direct path and reflected paths is the sum of values calculated by equation 49 for each
path.
Calculating the crosswind integrated concentration (CWIC) removes lateral diffusion as an issue for
the moment. In section 2.11 the horizontal distribution function (which is independent of the way in which
the vertical diffusion is handled in CTDMPLUS) will be discussed as will the method for coupling it to the
CWIC to obtain the concentration at a given receptor.
Li and Briggs (1988) examined a number of previously proposed PDF models and concluded that the
CWIC was most impacted by the form of P chosen; assumptions that these models make about the reflec-
tions of the particles from the surface and mixing lid were found to be of minor concern. Therefore they
elected to keep the reflection assumptions simple (perfect reflections assumed) and focus on developing
forms of Pa that best described laboratory and field observations.
The PDF formulation of Li and Briggs (1988) which is contained in CTDMPLUS is an asymmetrical
bi-Gaussian form that produces pollutant distributions very similar to those observed experimentally (Willis
and Deardorff, 1981). The bi-Gaussian form of Pu (equation 4 in Li and Briggs) is
/>- Gm exp[-uT2/(2o-2O] foru;'>0
= Gm exp[-iy'2/(2a2.)] foruT<0 (50)
2-42
-------
z.
Direct
Path
<2zi-zs> \
SSSS/S//SS/S///SSSSS/SS
z.
Ground
Reflection
SSSSSSSSS////S/////S//S
Lid-Ground
Reflection
sssssssssssssssss/sss/s
Lid
Reflection
SSSSSSSSSS/S/S/S/////SS
Ground-Lid
Reflection
ssssss/ssss/sssssssssss
Figure 2-12. Illustration of the use of image sources to model possible paths between the point (rs, zs) and the
point (rr,zr).
2-43
-------
where Gn is the maximum probability, u/'-(io- w~), and w is the mode of the vertical velocity. The
standard deviation of the velocities above and below the mode are denoted by ow. and ow. respectively.
The distribution in equation 50 is illustrated in Figure 2.13. The bi-Gaussian Pa is different from the stan-
dard Gaussian Pw in that for the bi-Gaussian w is generally less than zero ( w - w - 0 for Gaussian) and
the distribution is skewed ( o,,. - aa_ for the Gaussian).
The bi-Gaussian skewed distribution is obtained by introducing an asymmetry (skewness) factor
p-[(oa./oB)- !]-(! -OB-/O,.) where au-(au. + au.')/2 . With these definitions, equation 50 can
be evaluated with the following relationships (Appendix A of Li and Briggs, 1988):
(51)
(52)
= ow(l+p) (53)
-o^l-p) (54)
)1'2 (55)
Based on observations, Li and Briggs suggest approximating a,, (in equation 55) as 0.5 w. , where
to , is the convective scaling velocity. This approximation is used in CTDMPLUS.
The asymmetry factor, |3 , ranges from 0 to 1.0 with a value of zero yielding the standard Gaussian
distribution. A p value of 0.7 yields results very similar to the experimental distributions of Willis and
Deardorff (1981). (3 = 0.7 is used hi CTDMPLUS. Further refinement of |3 (perhaps making it a function
of release height) maybe practical after further study. For further details concerning the bi-Gaussian PDF,
the user is referred to the Li and Briggs paper.
A comparison of the vertical pollutant distribution based on the Gaussian PDF, the bi-Gaussian (BG)
PDF, and the tank experiments of Willis and Deardorff (1981) is seen in Fig. 2.14. The BG model repro-
duces the tank results very well with the maximum surface CWIC occurring with about the same magnitude
and at about the same location. The standard Gaussian results where the mode remains elevated do not
.compare as favorably.
The PDF model, as described, gives no consideration to alterations to the vertical pollutant distribu-
tion by terrain distortion to the flow. The following section describes the use of the flow model (already
present in the original CTDM code) to alter the trajectories of the pollutant particles in order to correctly
2-44
-------
Figure 2-13. Illustration of the bi-Gaussian distribution of Pw.
2-45
-------
N
1.0
" O TANK
X -<«/)»
Figure 2-14. Comparison of normalized crosswind integrated concentration contours as a function of Z = z/z\
and A' for a) Gaussian model, b) bi-Gaussian PDF model (revised from Li and Briggs, 1988), and
c) Willis and Deardorff (1981) tank data. Note the slightly expanded abscissa of figure c), but
each figure has X ranging from 0 to 2.
2-46
-------
connect the release points with the receptor points through a terrain distorted flow field. Since the convec-
tive vertical velocity producing a trajectory from the source to a receptor on a hill is different from that in
the absence of the hill, the probability (equation 50) and the crosswind integrated concentration (equation
49) will differ for the receptor on the terrain.
2.92 Terrain Influences on the PDF
A primary assumption in the PDF approach is that a particle is given an initial convective vertical
velocity at the source (or effective source) height which it "remembers" throughout its trajectory. When
complex terrain is not present and straight line flow is assumed, the flat terrain PDF described in section
2.9.1 is adequate to describe the CWIC distribution. When complex terrain is present, it is necessary to
account for the effect of perturbations in the flow field on the particle trajectories.
We know the positions of the source and the receptor but we do not know the convective vertical
velocity which, when added to the terrain-perturbed flow will result in a trajectory from source to receptor.
This must be found with an iterative approach. In each iteration a backwards trajectory is calculated from
the receptor toward the source based on the terrain-perturbed flow and the superimposed convective verti-
cal velocity. This trajectory is taken backwards until the x coordinate of the source is reached. A revised
estimate of the vertical velocity is based on how far the end point of the backwards trajectory is displaced
from the source position in the x-z plane. The iterations continue until a vertical velocity is found that con-
nects the source to the receptor through the disturbed flow.
Before beginning the iterative process, a first and second guess of the vertical velocity must be
obtained. A first guess for the vertical velocity, w,, is computed as the average of the unperturbed (no hill)
vertical velocities needed to reach receptors if they were placed at z = 0 and on a pole at z - z r (the
height of the receptor above base elevation). Defining receptor variables with subscript r and source vari-
ables with subscript s, these vertical velocities are calculated as:
to- £/z5/(x5-xr) for receptor atz = 0 (56a)
u;= i/(zj-zr)/(x,-xr) for receptor at z = zr (56b)
The first guess for the convective vertical velocity becomes the average of estimates in equations 56a and 56b
such that:
, = i/(2zs-zr)/2(xj-xr) (57)
2-47
-------
u>, is then used to estimate the second guess, w2. MI is used to calculate a back trajectory through
the perturbed flow from the receptor point ( xr, z r) to the point ( x, z ) where the flow first becomes
unperturbed (assuming that the source is upwind of the terrain perturbations). If the source is within the
perturbed flow region, the backward trajectory is calculated to the point ( xs, z) where the trajectory
crosses the y-z plane containing the source, and the point ( xs, z) should be used in place of ( x, z ) in the
following discussion.
Having determined the point ( x, z ), the second estimate of the convective vertical velocity, w z, is
made in one of two ways. The first is an intercept method, illustrated in Figure 2-15. m is defined as the
slope of the line tangent to the trajectory at the point ( x, z ) where flow perturbations become negligible.
Since both U and wt are constant in the unperturbed flow region, m = «;,/£/. This tangent line crosses
the x-axis at x, where x f = x - z /m . The second estimate of vertical velocity is based on the slope ( m')
of the line from the source to x, where m'- zs/(xs- x,)) such that:
x,-x.
u
When the slope m is small, causing numerical problems in equation 58, the intercept method is
replaced with a weighted average approach for estimating wz . This approach is illustrated in Figure 2-16.
Again, m is the slope of the line tangent to the trajectory at the point ( x, z ); m" is the slope of the line
connecting the source ( x s, z s) to the point ( x, y ). Then an average of m and m ", each weighted by
the x distance over which it applies, is used to estimate w 2:
U[m (x-xs) + m(xr-x)]
w2= (59)
xr-x,
These first two estimates of vertical velocity are used to begin an iterative process in which each sub-
sequent estimate is found, using a linear interpolation method (regula falsi technique).
Starting with w, and wz, the procedure begins by calculating backward trajectories from the recep-
tor to points ( x,, z, ) and ( x2, z2 ) respectively (refer to Figure 2-17). Using the slope of the tangent to
the trajectory at each of these points, the trajectories are continued backward until they intersect the x-y
plane of the source. (If the source is in the perturbed flow region, the points ( x,, z, ) and ( x2, z2 ) will
2-48
-------
point at which perturbations
are ignored
Figure 2-15. Illustration of the intercept method used for finding a second estimate of the vertical velocity (w)
that carries a parcel from (rs , zs) to the point (rr, zr) in a deformed flow over a hill.
2-49
-------
(xs-zs)
point at which perturbations
are ignored
m
m
(xr,zr)
Figure 2-16. Illustration of the weighted average approach used for finding a second estimate of the vertical
velocity (w) that carries a parcel from (xs , zs) to the point (xt, zr) in a deformed flow over a hill.
2-50
-------
x-plane
of source
(*S'zti)
(Xj^) point at which perturbations
(x2,z2) are ignored
(xr,zr)
Figure 2-17. Illustration of the linear interpolation method used iteratively to determine estimates of the verti-
cal velocity (w) that carries a parcel from (rs , zs) to the point (xr , zr) in a deformed flow over a
hill.
2-51
-------
already be in the x-y plane of the source). These points of intersection are denoted as ( x, , z, , ) and
( x , . 2(2 ) in Figure 2-17. The vertical difference between each of these points and the source point is calcu-
lated and defined as:
These values of f} and /2 measure how close the corresponding wl and u>2 are to the true value
producing a trajectory from source to receptor. The new estimate of vertical velocity, w 3 , is found by inter-
polating between w , and w 2 such that:
(60a)
With this new estimate of w the iterations continue by calculating a third backward trajectory and
proceeding as above, so that, in general, subsequent estimates are made using the following expression:
(60b)
As the correct value of w is approached in these iterations, the value of / is minimized; the itera-
tions continue until fN is negligible. Then WN will be the vertical velocity that will "send" the particle
along a trajectory through the disturbed flow to the receptor. It is this vertical velocity that is then used in
equation 49 to obtain the probability estimate for equation 50.
If the elevated terrain is removed or if a particular receptor is not associated with a terrain feature,
CTDMPLUS will return a flat terrain PDF result (as found in Li and Briggs, 1988).
In a convective boundary layer, there is a distance (travel time) downwind from the source where the
pollutant released from a point source will be sufficiently well mixed through the vertical layer that no signif-
icant vertical gradient of concentrations will exist. This distance will be proportional to the time scale asso-
ciated with the largest convective eddies, w,/zt and the horizontal wind speed. Beyond the point where
the pollutant is thought to be uniformly mixed, the effect of terrain is only important in terms of the
2-52
-------
horizontal distribution function. Therefore, for receptor downwind distance where X > 4, CTDMPLUS cal-
culates the CWIC as if the terrain is flat and the pollutant uniformly mixed in the vertical. Terrain effects on
the lateral separation of source and receptor are maintained at all downwind distances.
2.10 TERRAIN INFLUENCE ON MIXING HEIGHT
Current formulations for models in regulatory use provide, if any, a constant type of correction factor
which results in z i adjustments ranging from no change to values of z ( that follow the terrain. Some, for
example, apply a "half-height" correction to both the plume and to z , in unstable conditions. This approach
has merit in that it is based upon typical potential flow solutions for distortions of streamlines over a two-
dimensional ridge. However, consideration should be given to variations in controlling meteorology and ter-
rain shape. In addition, the constant correction factor approach has the effect of correcting the mixing
height over the entire distance from source to receptor even though the distortion may only occur over a
small portion of that distance.
For CTDMPLUS, we selected a somewhat different approach to modeling adjustments to the vertical
pollutant distribution caused by terrain induced distortions in the mixing height. This approach is based
upon a recent fluid modeling study (Perry and Snyder, 1989) of deflections of z t over terrain as functions of
terrain shape, upstream z i , and hill Froude number. The study results indicated that these parameters are
important for scaling the dynamic response of the mixed layer depth to changes in topography.
To simulate the mean flow in a daytime well-mixed layer beneath a deep stable layer, Perry and
Snyder used the large towing tank at the EPA Fluid Modeling Facility to measure vertical deflections of
streamlines over the crest of the model hills. They identified four significant dimensional parameters: z, ,
H (the hill height), U (the tow speed), and N (the Brunt-Vaisala frequency related to the elevated stable
layer). From these they established two dimensionless products, Fr - U/NH and z,/ H , which were
used to scale their results.
They used two terrain features to obtain some appreciation for the degree to which terrain shape
affects mixing height distortions. One feature (labeled ACCB) was an axisymmetric, polynomial-shaped hill
with a center height, H , of approximately 15% of the total depth of the tank and a half-height radial dimen-
sion of 2.5H . The second feature was an elongated (in the crossflow direction) version of the first, with a
half-height radial dimension of 5.0H .
A summary of the study is shown in figure 2-18. Figure 2-18a shows the deflection of z, as a func-
tion of both initial z f and hill Froude number for the axisymmetric hill. As the upstream height of the
mixed layer is increased, the influence of the hill is decreased and the deflections in the z i -stream line are
reduced. For example, when Fr - 1.54, a doubling of the upstream z, from one to two hill heights
2-53
-------
1.0
ACCB HILL
N
<
0.8 -
0.6 -
0.4 -
0.2 -
0.0
0.0 0.5 1.0 1.5 2.0 2.5
initial Zj/h
1.0
ELONGATED HILL
0.8 -
0.6 -
0.4 -
0.2 -
0.0
2.24
1.88
1.55
0.91
0.40
0.0 0.5 1.0 1.5 2.0 2.5
initial Zj/h
Figure 2-18. Fluid modeling results from Perry and Snyder (1989) showing deflections in the z-{ streamline as a
function of the initial z\ height and the hill Froude number for an (a) axisymmetric hill and (ft)
elongated hill.
2-54
-------
results In a 52% decrease in the deflection over the hill top. For a given mixing height, an increase in the
Froude number is associated with an increase in streamline deflections, signifying an increased ability of the
flow to overcome the elevated stratification and push up over the hill top.
Figure 2-18b shows results for the elongated terrain feature where, in contrast, the magnitude of the
deflections for any given F r and z t is greater than that for the symmetric hill. Since the elongated hill
"appears" more like a ridge to the flow than does the ACCB, there is less lateral spreading of the streamlines
over the crest and an increase in the vertical deflections.
In contrast to the half-height correction factor approach where the deflection in z t is always a con-
stant 0.5, figure 2-18 shows that, in fact, the amount of deflection is very much a function of meteorological
variables. Over the range of F r and zt/H contained in their experiments, Perry and Snyder found the
"correction factor" ranging from about 0.1 to 0.9. To account for these variations, CTDMPLUS uses the
results of Perry and Snyder (1989) to estimate the appropriate deflection to apply to the mixing height at the
hill crest. Where necessary, interpolation between values is accomplished. In the present version of
CTDMPLUS, the study results are contained in a lookup table. Analytical fits to the data may be added at a
later time.
To account for the effects of the zt deflection on the flow within the mixed layer below zt,
CTDMPLUS adjusts the flow model's stratification factor (Froude number) in order that the z t -streamline
deflection at the hill crest matches the fluid modeling results. Then, this adjusted flow model is used in the
flow distortion corrections to the PDF model of the vertical distribution function. This prevents the flow
model from calculating particle trajectories (based on streamline distortions) that might otherwise escape
the mixing layer. This type of bulk correction to the flow model calculations accounts for the influence that
the elevated stable layer has on the vertical extent of the flow in a mixed layer over terrain.
In order to calculate the hill Froude number and thus apply the corrections described above, the
Brunt-Vaisala frequency, N , for the layer above z t must be calculated. This requires an estimate of the
potential temperature gradient in that layer. The following section describes the simple method used in
CTDMPLUS to estimate that gradient.
2.10.1 Estimate of the Elevated Temperature Gradient
The gradient of potential temperature in the stable layer above the mixed layer is an important factor
in determining the potential for plume penetration and for estimating the magnitude of the deflection in z,-.
In most cases, measurements of the gradient in this elevated layer are not available for any given hour of the
2-55
-------
day since tower measurements are normally limited to heights of 100 meters or so. Therefore CTDMPLUS
is designed to store the morning temperature profile (provided in the RAWIN input file) and to use it
throughout the day for estimating the gradient above zf. The following approach is used.
Consider the morning temperature profile depicted in figure 2-19. As the sun rises and begins to heat
the surface, a convectively well mixed (uniform temperature) layer develops near the surface and "grows"
throughout the day. However, the temperature profile in a layer above z ( changes relatively little from that
in the morning sounding. Of course this assumes that there is neither significant subsidence nor cold or
warm air advection occurring in the upper layer. Having stored the morning sounding, for a given hour of
the day the model can call for the temperature gradient that is appropriate for the layer just above z,
(model currently uses a layer from zf to zt + 500m).
Field measurements (e.g. the Wangara experiment, Clark et al., 1971) of observed profiles through-
out the day lend support to this concept of using the the morning sounding throughout the day for upper
level temperature gradient estimates. These data point out the relative invariance of upper level
temperature data even during periods of intense heating at the surface.
2.11 THE HORIZONTAL DISTRIBUTION FUNCTION
Since significantly non-Gaussian plume shapes are not generally observed in relation to lateral diffu-
sion, analytical diffusion models have retained the usual Gaussian form (Briggs,1985). Therefore,
CTDMPLUS uses the standard Gaussian form for the horizontal distribution function (HDF)
HDF = \,2Q exp[-0.5(y/ay)2] (61)
where y is the effective crosswind distance from plume centerline and a y is the lateral dispersion
coefficient. Equation 61 can be combined with the crosswind integrated concentration, CWIC (discussed in
section 2.9) to extract the field of concentrations as a function of downwind and crosswind distance. The
concentration at any receptor (xr,yr,zr) is then
CI//C
C(xr,yr,zr) = exp
-0.5|
ay
2
(62)
The flow model used to calculate the trajectories for the PDF submodel (section 2.9) also provides
the lateral distance (y, - y') between the source and the trajectory passing through the receptor because
the endpoint of the trajectory returned by the submodel is at (x, y', z). This distance ( y s - y') differs
2-56
-------
2000
1800
1600
1400
1200
'1000
800
600
400
200
276
^^^ ^^ ^^"
Potential temperature (K)
292
Figure 2-19. Typical profile of mean potential temperature showing the growth of the well mixed layer
throughout the day. At each time of the day, the potential temperature gradient above z\ does not
differ significantly from the gradient of that level in the morning sounding.
2-57
-------
from the geometric source-receptor separation (y, - y r) because the deflection of the flow around the
terrain alters the relative position of the receptor within the streamlines of the flow. Therefore y of
equation 62 is set equal to ( y, - y') remembering that y' is related to y through the flow model calcula-
tions. For flat terrain (no flow distortion) y reduces to the expected ( y, - y r).
The remaining parameter to specify in equation 62 is the lateral dispersion coefficient, 0.10), is:
Oy = 1.6ztX2'3Fl'3 (63)
where
z, = mixing height
v, %
X - - non-dimensional downwind distance
(Z \l/3
^-f2 I = convective scaling velocity
L = Monin-Obukov length
u. = friction velocity
and
Ft> = plume buoyancy flux.
2-58
-------
Briggs found that equation 63 fits the data reasonably well except for very light wind conditions where
U / w. < 1.5 (these conditions may need to be handled separately but are not currently so in
CTDMPLUS). In addition, Hanna found that at certain distances downwind (where perhaps the buoyancy
is no longer important) the data better supported a linear growth with non-dimensional distance of the form
ay = 0.6ztX - (64)
Therefore, for buoyant plumes in a convective layer, the lateral dispersion coefficient is modeled in
CTDMPLUS with equation 63 except when X > 1 (or X/ F. > 17) where equation 64 is used.
For the case of passive plumes ( F, < 0.10), the choice of a y is based on the results of studies con-
cerned primarily with non-buoyant, elevated releases. Briggs (1985) accumulated data from a number of
field studies including Prairie Grass (Barad, 1958), Cabauw (Agterburg et al., 1983), and Condors (Kaimal
et al., 1986), laboratory results of Willis and Deardorff (1976,1978,1981), and numerical modeling by Lamb
(1979). Briggs suggested a reasonable fit to the bulk of the data to be
0:62:,*
* (1 + 2JQl/z
Others (e.g. Venkatram and Vet, 1981) have observed the slow departure of o y from a linear growth
with X and have proposed simple linear formulas. The scatter in the data is such that any number of possi-
ble formulas could be justified. However equation 65 is an adequate fit to a majority of the data with the
added attraction of being a particularly good fit to the convective tank data; thus it is used as the passive
plume lateral dispersion coefficient for CTDMPLUS.
2-59
-------
-------
SECTION 3
USER INSTRUCTIONS FOR RUNNING CTDMPLUS
3.1 OVERVIEW
CTDMPLUS is a refined model designed for complex terrain applications. Accordingly, a consider-
able amount of data preparation needs to be done prior to running the model, especially for terrain and
meteorological data. Figure 3-1, a schematic of the system flow, shows the input files for each of the
programs.
CTDMPLUS requires a mathematical depiction of the shape of each terrain feature (hill) that is as-
sociated with at least one model receptor. The terrain preprocessor programs must be run to generate the
terrain input to CTDMPLUS unless the required hill shape functions are already known. The terrain
preprocessor user's guide (Mills et al., 1987) gives a complete description of the procedure for generating
the required terrain input file. A brief summary is given here.
The user first obtains a topographic map of the area of interest and selects the specific features (iso-
lated hills or distinct portions of "complex" hills) that are to be modeled. Then, several contours must be
digitized; that is, (x, y ) points along the contour are determined with enough resolution so that the line
segments connecting the points are a good representation of the actual map contour. The vertical spacing
between contours that are digitized should be small enough so that detailed information about the shape of
the hill is retained. However, unnecessary detail and work could result if every contour were digitized. In
general, at least 5-10 contours for each hill is desirable. A maximum of 21 contours is imposed by the ter-
rain preprocessor. The actual digitization can be done on a graphics board or can be done by an outside
vendor, as long as the format specifications required by the terrain preprocessor are followed.
The terrain preprocessor programs FITCON and HCRIT should be run first. Both FITCON and
HCRIT require an options file, "FOPTIONS" and "HOPTIONS", respectively. The formats of these options
files are shown in Table 3-1 and Table 3-2. Once FITCON has been run to read and process the digitized
contours, the user has the terrain information necessary to run the receptor generator program, RECGEN,
if desired (RECGEN is described further in Section 4.1). The output file from HCRIT contains the mathe-
matical hill descriptions in the correct format for input to CTDMPLUS. CTDMPLUS requires that the
lowest critical elevation be at or below the common stack base (the lowest base elevation of the
meteorological tower and stacks). If the defined hill does not extend to the common base, HCRIT can be
used to extrapolate the hill downward to that level by using mode 2 (modes are discussed in the terrain pre-
processor user's guide). The mathematical fit to the actual contours can be plotted by running a third pro-
gram, PLOTCON.
3-1
-------
K>
^Digitized Contours^)
FITCON
Plot
File
HCRIT
PLOTCON
KEY
Progra/ni
------ Opllontl
RECGEN
^
Terrain )
- , ^
c
k
CDutpuTN.
UX
GRIODAT
CHIDIS
Output
File
CONTOUR
Figure 3-1. Interaction among components of the CTDMPLUS system.
-------
TABLE 3-1. FORMAT OF THE "FOPTIONS" FITCON OPTIONS FILE
Record
Group
1
2
3
4
5
6
7
8
9*
10
11**
12**
13***
14
15*
16
Format
A12
A12
12
A8
F10.4
F10.4
F10.4
Al
A4.1
11
14
14
A12
Al
A12
A12
Description
Contour master file name
Diagnostic output file name
Hill identification number
Hill name
Hill top elevation (user units)
Hill center x-coordinate (user units)
Hill center y-coordinate (user units)
Angular filtering option (y or n)
Angular filtering size
Contour selection mode
Lower bound of contour identification
Upper bound of contour identification
numbers
numbers
Contour identification number file name
Plot file option (y or n)
Plot file name
Fitted output file name
Input value required only:
*If this option is selected
**If selection mode = 2
***If selection mode = 3
3-3
-------
TABLE 3-2. FORMAT OF THE "HOPTIONS" HCRIT OPTIONS FILE
Record
Group
1
2
3
4*
5
6**
7**
Format
A12
A12
Al
A12
Al
A3
F10.4
Description
FITCON output file name
HCRIT output file name
Plot file option (y or n)
Plot file name
Selection mode for critical elevations
Number of critical elevations
Lowest critical elevation (user units)
Input value required only:
*If this option is selected
* *If selection mode = 2
3-4
-------
Meteorological data is input to CTDMPLUS from three separate files. The first of these files, TRO-
FILE", is created by the user and should represent user-validated and quality-assured meteorological data
consisting of conventionally available measurements of
wind direction;
wind speed [scalar, vector (optional)];
air temperature;
standard deviation of the horizontal wind direction (ae) or of lateral wind speed (au ); and
standard deviation of the vertical wind speed ( a ).
These measurements can be provided at a number of levels (CTDMPLUS can accept up to 50 levels)
and should extend upward into the region typically occupied by the plumes to be modeled.
The second meteorological data file, "RAWIN", consists of rawinsonde data created from a National
Climatic Data Center TD6201 file by using the READ62 preprocessor [discussed in the meteorological pre-
processor user's guide (Paine, 1987)]. The "RAWIN" file is required by CTDMPLUS only if the model is
being run for unstable hours. The "RAWIN" file contains upper air measurements of
air pressure,
height above local ground,
air temperature,
wind direction, and
wind speed.
CTDMPLUS uses this pressure and temperature information to determine the stability above the mixing
height during daytime conditions. The meteorological preprocessor, METPRO, also uses the "RAWIN"
data for calculating daytime mixed layer heights.
A third file required by CTDMPLUS, "SURFACE", consists of the following processed meteorologi-
cal variables:
Monin-Obukhov length (L),
friction velocity (u, ),
nighttime surface layer height (h ), or daytime mixed layer height (z,), and
3-5
-------
surface roughness length ( z0).
METPRO produces the "SURFACE" variables from the meteorological data in "PROFILE" and
"RAWIN" and from additional radiation information (net or total incoming radiation or cloud cover) and
surface characteristics [see METPRO user's manual (Paine,1987) for details].
The file "CTDM.IN" contains model switches, meteorological tower coordinates, and source informa-
tion. Receptor information is supplied to CTDMPLUS in a file named "RECEPTOR". This file can be
created with RECGEN, the receptor generator program, or with a text editor.
There is an optional input file to CTDMPLUS containing hourly emissions data ("EMISSION"). If
the emissions do not vary, then this file is not necessary and the constant emissions data are read from
"CTDM.IN".
CTDMPLUS executes for as many hours as are contained in the input meteorological data files. It
can be run for a period as short as a single hour. Data in multiple-hour runs must be sequential. Informa-
tion in the data files "SURFACE", "PROFILE", "EMISSION", and "RAWIN" must be chronologically con-
sistent.
When CTDMPLUS is run, it produces a concentration file, in either binary or text format (user's
choice), and a list file containing a verification of model inputs and case study output (also at the user's op-
tion). Another option is to print a table listing the contribution from each source to the total concentration
at a receptor for each hour. A summary of the top four one-hour concentrations at each receptor for the
entire run is included at the end of the list file if that display option is selected. The concentrations can be
plotted onto a display of terrain contours by the CHIDIS postprocessor (see Section 4.3). Concentration
isopleths can be displayed using the CONTOUR program available through the menu driver system (see
Section 5.2.8).
The input files required for CTDMPLUS are described in detail in Section 3.2. A menu driver pro-
gram that can be used to run the system using pre-exisiting files is described in detail in Section 5. An inter-
active program (SETUP) exists to aid users (especially new users) in constructing these files for short runs
of CTDMPLUS. The SETUP program is described in Section 4.2.
3.2 INPUT DATA REQUIREMENTS
There are five required input files and two optional input files for CTDMPLUS. The five required
files consist of:
a general file of program specifications, which consist of program switches, source data, meteo-
rological tower coordinates and hill surface roughness lengths ("CTDM.IN");
3-6
-------
a terrain data file which is obtained directly from the HCRIT terrain preprocessor ("TER-
RAIN");
a file containing receptor names, locations, and the associated hill numbers ("RECEPTOR");
a surface meteorological data file which is obtained directly from the METPRO meteorological
preprocessor program ("SURFACE");
a user-created meteorological profile data file which contains conventional meteorological data
measured at multiple levels ("PROFILE");
The optional input data files consist of:
a file of hourly emissions parameters ("EMISSION");
a file containing upper air data from rawinsonde data which is obtained from the READ62 pre-
processor ("RAWIN").
The input files are discussed in more detail in Sections 3.2.1 through 3.2.7. CTDMPLUS output files
are discussed in Section 3.3.
3.2.1 General Program Specifications
The input file "CTDM.IN" contains program options, meteorological tower coordinates, source data,
and surface roughness lengths for each hill. The inputs and file formats are listed in Table 3-3 (a sample
input file is shown in Figure B-l). Some program options control the amount of model output in the follow-
ing areas:
case-study printout (if used, a voluminous output results; use this only for short runs)
summary table for the current run (if called, a table of the top four one-hour concentrations at
each receptor is printed)
file (either text or binary) of predicted concentrations for postprocessing
hourly source contribution table.
Other program options, which control the use of input or output data, tell CTDMPLUS to:
choose which mixed layer height value to use if both measured or calculated mixed layer heights
are provided (case-specific; in general, measured values should be used if available and of good
quality)
set the minimum scalar wind speed to 1 m/s (recommended for most applications)
3-7
-------
TABLE 3-3. CONTENTS OF THE 'CTDM.IN" FILE
Line
Group #
1
2
3
4
5**
6
7
Variable
Name
Tide
ICASE
ITOPN
ICONC
EMIX
IWSI
ISIGV
IWD
icmo
ISOR
IUNSTA
HORIZ
VERT
RLAT
RLON
TZONE
IPOL
LABEL
XT
YT
ZT
SNAME
XS
YS
ZS
HS
DS
TS***
VS***
Q***
IQ
ENDS
ZOH
Columns
1-80
*
*
*
*
*
*
*
*
1-20
21-30
31-40
41-50
1-16
17-23
24-30
31-37
38-44
45-51
52-58
59-65
66-72
80
1-4
*
Format
A80
*
*
*
*
*
*
*
*
A20
F10.0
F10.0
F10.0
A16
F7.0
F7.0
F7.0
F7.0
F7.0
F7.0
F7.0
F7.0
11
A4
*
Description
80 Character header
Case study printout option. O=No, 1= Stable hours only,
2 = Unstable hours only, 3 = All hours
Create a top 4 table at the end of the run. 0 = NO,
1=YES
Concentration output file option. 0=NO, 1=BINARY,
2=TEXT, 3=TEXT with receptor information
0=use calculated mixing heights as first priority, l=use
observed mixing heights as first priority
Set minimum wind speed = 1.0 m/s. 0=NO, 1=YES
Assume a, input if 0; av input if 1
Scale wind direction with height. 0=NO, 1= YES
0= output concentrations (\i g/mi>)
1 = output chi/Q (p. s/m»)
Create a source contribution table at the end of the run. 0
= No, 1 = Yes
1 = Model unstable hours (RAWIN file required) and
stable hours
0 = Model only stable hours
Horizontal scale factor, converts user units to meters
Vertical scale factor, converts user units to meters
Site latitude (degrees)
Site longitude (degrees)
Site time zone (hours behind GMT)
Pollutant Code (1-4)
Meteorological tower name
X-coordinate of tower (user horizontal units)
Y-coordinate of tower (user horizontal units)
Z-coordinate (user vertical units) of tower base
16 character source name
X-coordinate of source, (user horizontal units)
Y-coordinate of source (user horizontal units)
Source base elevation (user vertical units)
Stack height (m)
Stack diameter (m)
Stack gas temperature (K)
Stack gas exit velocity (m/s)
Emission rate (g/s)
Variable emission rate flag for this stack (0= constant,
1-variable)
'ENDS'. Flag for end of source data.
Surface roughness length (m) for each hill in the order
that the hills appear in the TERRAIN file.
j*Free format
**One line per source, maximum of 40 sources (maximum number can be changed)
***These values are replaced if hourly emissions data are provided for this stack.
3-8
-------
indicate which variable (either a or a ) is provided in "PROFILE" for the stable/neutral
crosswind turbulence measurements
scale wind direction with height (recommended unless the user is trying to "aim" the plume at a
hill in a test run)
display concentrations in units of mass per volume (chi) or chi divided by the emission rate
(x/Q)
model all hours: stable, unstable, and neutral
The user is also asked to provide other program constants, such as:
factors by which to multiply user units to obtain meters for both vertical and horizontal coordi-
nates;
site latitude, longitude, and time zone;
pollutant code number for hourly emissions (1-4).
"CTDM.IN" includes the position (x,y,z) of the meteorological tower for the data contained in "PRO-
FILE". The horizontal position, if on or near a hill modeled by CTDMPLUS, is used by the model to deter-
mine the upwind wind direction and speed. The vertical (tower base) information is used to reference the
meteorological measurements relative to the common stack base; "PROFILE" contains the measurement
heights relative to the local ground surface at the tower site.
Point source information (in "CTDM.IN") includes stack name, horizontal and vertical coordinates,
stack height and diameter at the outlet, stack gas temperature, and exit velocity. Variable emissions for any
subset of the total number of stacks is allowed. Emission rates for up to four pollutants may be included in
the variable emissions file ("EMISSION"). The specific emission rate used for a given run is chosen via a
pollutant code number. CTDMPLUS does not require stacks to be co-located, but a common stack base is
retained for convenience; this is calculated to be the minimum of the tower base and the lowest stack base
among those input. The lowest "critical elevation" specified in the terrain preprocessor HCRIT run for each
hill must be at or below the common stack base to avoid a CTDMPLUS runtime error. The elevation of
each stack top is preserved by adjusting the stack height if the base elevation is changed. Up to 40 sources
can be input; this number can be altered by changing the code in the INCLUDE file "PARAMS.INC".
Surface roughness lengths for the local surface characteristics of each hill are given in "CTDM.IN".
The values of these roughness lengths vary according to vegetative cover and season of the year. See Table
3-4 for guidance.
3-9
-------
TABLE 3-4. EXAMPLE SURFACE ROUGHNESS LENGTHS,
FOR LAND USE TYPES AND SEASONS (METERS)
Land Use Type
1. WATER (FRESH WATER
AND SEA WATER)
2. DECIDUOUS FOREST
3. CONIFEROUS FOREST
4. SWAMP
5. CULTIVATED LAND
6. GRASSLAND
7. URBAN
8. DESERT SHRUBLAND
Spring
0.0001
1.00
130
0.20
0.03
0.05
1.00
0.30
Summer
0.0001
130
130
0.20
0.20
0.10
1.00
030
Autumn
0.0001
0.80
1.30
0.20
0.05
0.01
1.00
030
Winter
0.0001
0.50
130
0.05
0-01
0.001
1.00
0.15
DEFINITIONS OF SEASONS:
"Spring" refers to periods when vegetation is emerging or partially green. This is a transitional situation
that applies for 1-2 months after the last killing frost in spring.
"Summer" applies to the period when vegetation is lush and healthy, typical of midsummer, but also of other
seasons in locations where frost is less common.
"Autumn" refers to a period when freezing conditions are common, deciduous trees are leafless, crops are
not yet planted or are already harvested (bare soil exposed), grass surfaces are brown, and no snow is present.
"Winter" conditions apply for snow-covered surfaces and subfreezing temperatures.
3-10
-------
322 Terrain Data
The terrain data file is created by the terrain preprocessor program and is used by CTDMPLUS with-
out modification. This file contains the hill center, hill height, major axis orientation, ellipse major and
minor axis lengths, and inverse polynomial lengths and exponents which are used to define the hill used by
the WRAP and LIFT components of CTDMPLUS. The base of each hill must be at or below the common
stack base. HCRIT, when run in mode 2, allows the user to select a range of elevations between which
critical elevations will be spaced. This option allows the hill to effectively be extrapolated down to the com-
mon stack base. Table 3-5 shows the format of the "TERRAIN" input file; an example is shown in Figure
B-2.
323 Receptor Data
The file "RECEPTOR", containing receptor names, coordinates, and hill number, is read by
CTDMPLUS. This file can be used directly from the output of the receptor generator, RECGEN (see Sec-
tion 4.1) or can be altered or created using a text editor. The format of this file is shown in Table 3-6 (an
example is shown in Figure B-3). Up to 400 receptors can be input to CTDMPLUS; this number can be
altered by changing the code in "PARAMS.INC".
Each receptor must have a hill number assigned. Those receptors which are beyond a digitized hill
region (outside the "skirt" of any hill) should be assigned a hill number of 0; they will be modeled as if in flat
terrain.
32.4 Meteorological Profile Data
The meteorological data file, "PROFILE", consists of hourly averaged values of wind, temperature,
and turbulence data. This file must be provided by the user. Each record represents measurements at a
single height. The data for each hour are grouped together with the highest level for any hour indicated by
setting a flag from 0 to 1. There is no restriction on the number of hours, nor need they be contiguous;
however, the dates/hours of the data must be sequential. The input records for each hour must be in order
of increasing height. The heights do not have to be the same from hour to hour.
The data in the "PROFILE" data file are read in free format. Missing data values are indicated by
-999. There must be at least one non-missing value for each variable for each hour [although a ( a ) and
ou values are not needed for hours in which all plumes are released into unstable layers]. The order of the
variables in "PROFILE" is given in Table 3-7; a sample file is shown in Figure B-4.
3-11
-------
TABLE 3-5. FORMAT OF THE "TERRAIN" INPUT DATA FILE
(FROM THE TERRAIN PREPROCESSOR)
Record
Group
1**
1
1
1
2*
2
2
2
3*
3
3
3
3
Parameter
Name
NH***
NZ****
HTP
HNAME
ZH
XHW.YHW
MAJORW
MAJAXW
MINAXW
ZH
L
MAJORL
EXPOMA
EXPOMI
SCALMA
SCALMI
Columns
6-7
9-10
21-30
31-45
MO
11-30
31-40
41-60
1-10
10-30
31-40
41-60
61-80
Format
12
12
E10.4
A15
F10.3
2F10.3
F103
2F103
F10.3
2E10.4
F10.3
2F103
2F10.3
Description
Hill identification number
Number of critical elevations
Hill top elevation (user units)
Hill name
Critical elevations (user units)
x,y-coordinates of the ellipse centroid for the
critical elevation
Orientation (degrees) of the ellipse major
axis with respect to north
Semi-major and semi-minor axes lengths for
the ellipse at the critical elevation
Critical elevation (must match critical eleva-
tions in Record Group 2)
x,y-coordinates for the fitted cutoff hill cen-
troid
Orientation of the fitted cut off hill major
axis with respect to north (degrees)
Inverse polynominal exponent parameters
for the major and minor fitted hill axes
Inverse polynominal length scale parameters
for the major and minor fitted hill axes
*There are NZ records for group 2 followed by NZ records for group 3
**Record groups 1-3 are repeated for each hill
***There is a maximum of 25 hills
**"There is a maximum of 21 hill contours
3-12
-------
TABLE 3-6: FORMAT OF THE "RECEPTOR" INPUT DATA FILE*
Record
Group
1**
Variable
Name
RNAME
XR
YR
ZR
GE
NH***
Columns
1-16
21-30
31-40
41-50
51-60
61-65
Format
A16
F10.0
F10.0
F10.0
F10.0
15
Description
16-character receptor name
x-coordinate of receptor (user horizontal units)
y-coordinate of receptor (user horizontal units)
Height of receptor above local ground surface
(user vertical units)
Ground-level elevation (user vertical units)
Hill number of this receptor
*No special line is required to signify the end of receptor input; this is signified by the end of the file
**One line per receptor, maximum of 400 receptors
***Hill number 0 is used to indicate flat terrain algorithm to be used for this receptor
3-13
-------
TABLE 3-7. FORMAT OF USER-CREATED "PROFILE" FILE
Variable*
Name Description
JYR Year (two digits)
JMO Month (1-12)
JDY Day of month (1-31)
JHR Hour at the end of the period (1-24)
HT** Height of this observation above the tower base (m)
lEND 0 if not the highest (last) level; otherwise 1
WD*** Wind direction (deg)
WS*** Scalar wind speed (m/s)
TA*** Ambient dry bulb temperature K
SIGTH*** Sigma-theta (deg) or sigma-v (m/sec), determined by a switch in "CTDM.IN"
SW*** Sigma-w(m/s)
UV* * * Vector wind speed (m/s)
*All variables are read free format; one line per height level
**There is a maximum of 50 heights for each hour; heights must be in ascending order.
* "Missing value is indicated by -999.
3-14
-------
Meteorological Surface Data
The meteorological surface data file, "SURFACE", is created by the meteorological preprocessor and
is input directly to CTDMPLUS. It contains hourly values of mixed layer heights, surface friction velocity,
Monin-Obukhov length, and surface roughness length. The chronological sequence of the hourly data must
match the "PROFILE" data sequence. The format of the "SURFACE" file is given in Table 3-8; a sample is
shown in Figure B-5.
3.2.6 Meteorological Rawinsonde Data
The meteorological rawinsonde data file, "RAWIN", is created by the meteorological preprocessor,
READ62, from an NCDC TD-6201 file. It contains upper air measurements of pressure, temperature, wind
direction, and wind speed. If IUNSTA = 1 in "CTDM.IN", then (for unstable conditions) the temperature
values are read from the "RAWIN" file by CTDMPLUS to calculate the potential temperature profile above
the daytime mixing height. The sequence of daily 12Z soundings must match the date sequence of the
"SURFACE" file data. If no unstable hours are included in the data, the user sets IUNSTA equal to zero
and the RAWIN file is not required by CTDMPLUS. The format of the "RAWIN" file is given in Table 3-9;
a sample is shown in Figure B-6.
32.7 Hourly Emissions Data
Hourly source emission values (emission rate, stack gas exit velocity, and stack gas temperature) can
be input to the model for any subset of the modeled sources via the file "EMISSION". There is one line per
hour for each source with varying emissions. The sources with variable emissions must be in the same order
as originally input in line group 5 of the program specifications file "CTDMJN" (Table 3-3). CTDMPLUS
will check for time inconsistencies between the meteorological data and the hourly emissions data. Table
3-10 gives the format of the hourly emissions file.
33 CTDMPLUS OUTPUT FILES
33.1 CTDMPLIIS Output Listing
An output listing is always created by the model; the content is determined by the case-study printout
option. If the case-study printout option is not selected, then the listing is limited to a verification of input
data from "CTDM.IN", a line printer map showing the relative locations of sources and receptors, and the
information contained in the TERRAIN" input file.
The case-study-mode output listing also includes:
input meteorological data from "SURFACE" and "PROFILE"
3-15
-------
TABLE 3-8. FORMAT OF "SURFACE" FILE (FROM THE METEOROLOGICAL PREPROCESSOR)
Variable*
Name Description
IYR** Two-digit year
MO** Month (1-12)
IDY** Day of month (1-31)
JDAY Julian day (1-366)
IHR** Hour (at end of period, 1-24)
ZIOBS Observed mixing height (m) from on-site measurement
ZIPRE Calculated mixing height (m)
USTAR Surface friction velocity (m/s)
EL Monin-Obukhov length (m)
ZO Hourly surface roughness length (m)
*All variables free format
**Date and times must match time variables in PROFILE data file.
NOTE: Negative values are written to "SURFACE" for missing mixed layer heights.
3-16
-------
TABLE 3-9. FORMAT OF THE "RAWIN" FILE (FROM READ62 PROGRAM)
Record
Group Format
1* (A4) "6201"
A5
12
12
12
12
12
12
2** F6.1
F5.0
F5.1
D
D
Description
Series "6201" label
Station identification number
Year of sounding
Month of sounding
Day of sounding
Hour of sounding (GMT)
Number of data levels in sounding
Number of data levels extracted
Pressure (mb)
Geopotential height (m)
Temperature (K)
Wind direction (deg)
Wind speed
*Repeated each hour
"Repeated for each level
3-17
-------
TABLE 3-10. VARIABLE EMISSIONS INPUT FORMAT*
Line#
For each hour, the
user should provide
one line of data for
each source having
variable emissions.
Variable
JYR
JMO
JDY
JHR
IS
TS
VS
QS(1)
QS(2)
QS(3)
QS(4)
Description
Year
Month
Day
Hour (ending time for the hour)
Stack number
Stack gas temperature for this hour, K
Stack gas exit velocity for this hour, m/s
Emission rates for up to four pollutants
for the hour, g/s (ordering of pollutant
emission rates should be consistent with
pollutant code number, which is an input
parameter to CTDM)
*Each line is read in free format.
3-18
-------
stack data for each source
meteorological variables at plume height
geometrical relationships between the source and the hill
plume characteristics at each receptor for stable hours:
distance in along-flow and cross-flow directions
effective plume-receptor height difference
effective a y, a z values, both for flat terrain and the hill-induced case (the difference
shows the effect of the hill)
concentration components due to WRAP, LIFT, and FLAT
plume characteristics at each receptor for unstable hours:
distance in along-flow and cross-flow directions
horizontal distribution function
- «,
crosswind integrated concentration.
Two other tables may be printed in the output file depending on the selection of options in
"CTDM.IN". The user may select the ISOR option, which will print a source contribution table for every
hour. If the user selects the TOPN option, a summary table of the top four concentrations at each receptor
is given (See Figure B-6 for an example). The ISOR and TOPN switches are independent from the ICASE
(case study printout) switch. A sample of the complete CTDMPLUS output listing file is given in Figure
B-8.
332 Concentration File
A separate disk file of predicted concentrations ("CONC") is written if the user chooses this option.
Three forms of output are possible:
1) a binary file of concentrations, one value for each receptor in the hourly sequence as run
(ICONC = 1 in CTDM.IN);
2) a text file of concentrations, one value for each receptor in the hourly sequence as run
(ICONC = 2); or
3) a text file as described above, but with a listing of receptor information (names, positions, hill
number) at the beginning of the file(ICONC = 3).
3-19
-------
Hourly information provided to these files (besides the concentrations themselves) includes the year,
month, day, hour, receptor number with the highest concentration, the number of receptors, and the concen-
tration units. The file formats are listed in Table 3-11 and a sample text file is shown in Figure B-7. The
concentration file will have negative values at all receptors for hours in which CTDMPLUS cannot calculate
values. Only concentration files of the second type can be used with the CHIDIS and CONTOUR post-
processors.
3.4 ADDITIONAL COMPUTER NOTES
The size of the CTDMPLUS executable file on an IBMR PC or compatible is approximately 360K
bytes. This size is dependent upon the maximum number of sources, receptors, hills, etc., that are allowed.
These maximum values can be changed by editing the file "PARAMS.INC" and then recompiling and relink-
ing the CTDMPLUS subroutines.
Default input/output unit numbers have been assigned for access to disk files:
TERRAIN" - unit 2
"EMISSION" - unit 3
"RECEPTOR" - unit 4
"CTDM.IN" - unit 5
"CTDM.OUT" - unit 6
"SURFACE" - unit 7
"CONC" - unit 9
"PROFILE" - unit 11
"RAWIN" - unit 14.
These assignments can be changed in the CTDMPLUS main program.
Program execution will halt if any of several input problems is encountered:
Time inconsistencies among data in files "SURFACE", "PROFILE", "RAWIN", and "EMIS-
SION";
Heights in the "PROFILE" file that are negative or not in increasing order;
Read errors for various portions of the input;
Too many sources or receptors; or
Hill numbers out of sequence.
3-20
-------
TABLE 3-11. FORMAT OF TEXT CONCENTRATION OUTPUT FILE
Line
Group #
Columns
Format
Description
1* (appears only once,
before the first hour's
concentration values)
2 (each hour)
3* (one line per recep-
tor, each hour)
4** (eight values per
hne until receptor list
is exhausted, each hour)
1-4 14 Receptor number
6-13 F8.0 Receptor x-coordinate (m)
15-22 F8.0 Receptor y-coordinate (m)
24-30 F7.1 Receptor height above local ground sur-
face (m)
32-37 F6.1 Receptor relief height above common
stack base elevation (m)
39-42 14 Hill number associated with this receptor
1-5 15 Year
6-10 15 Month
11-15 15 Day
16-20 15 Hour
21-25 15 Number of receptors with maximum con-
centration
26-30 15 Number of receptors
31-46 A15 Units of concentration
1-4 15 Receptor number
6-15 E10.4 Concentration (\ng/m3 or \ns/m3,
depending on user option in "CTDM.IN")
1-10,11-20 8E103 Concentration (ng/m3 or n s / m 3 ,
71-80 depending on user option in
"CTDM.IN")
"Used only if ICONC = 3 (receptor information written at beginning of file)
**Used only if ICONC = 2 (no receptor information written at beginning of file)
3-21
-------
The messages associated with these errors are self-explanatory, and the user should be able to readily
fix the problem and re-run the model.
Another group of program halts is attributable to potential, though rare, numerical convergence
problems in the LIFT component. These problems for the most part have been solved, and thousands of
test hours have been run without error. However, future users may encounter one of two possible problems:
"ENDLESS LOOP IN PATH" (nonconvergence problem)
"NO CONVERGENCE IN MIX ROUTINE".
These program halts may, in the future, be handled by having CTDMPLUS not predict concentra-
tions for the hour affected (negative concentrations being produced). Until more experience is gained with
CTDMPLUS, these program halts are being retained; our current experience indicates that they will be
rarely encountered. If any of these program halts does occur, the case should be documented with input and
output files and forwarded to the authors for analysis (MD-80, U.S. EPA, RTP, NC 27711).
3£ CTDMPLUS SUBROUTINE STRUCTURE
A brief overview of the structure of the major subroutines in CTDMPLUS is presented in this sec-
tion. It is designed to acquaint the user with how each of the subroutines "fits in".
The main program, which is named CTDMPLUS, calls several subroutines which read in much of the
preliminary data for a run and perform several setup operations, producing much of the printed output from
the run. It then calls the major subroutine SEQMOD which actually performs all of the sequential compu-
tations. Figure 3-2 outlines the form of CTDMPLUS (MAIN).
Each of the minor subroutines referred to in the figure performs a well-defined task:
PAGE: Skips to a new page and writes a page header.
INPAR: Reads and writes the input/output switches, conversion factors, and constants for the
model. These values are passed through "PARAMS.CMN".
INPSOR: Reads and writes all source data that remain fixed (constant) over the period of simu-
lation. The source data are passed through "STACKS.CMN".
INPREC: Reads and writes receptor data. These data are passed through "RECEPT.CMN".
INPTOW: Reads and writes position of meteorological tower. The position is passed through
"TOWER.CMN".
3-22
-------
CTDMPLUS (MAIN)
initializations
- Open: CTDM.IN, CTDM.OUT
- Call PAGE
Call INPAR
Call INPTOW
Call INPSOR
Open: EMISSION (if variable)
Open RECEPTOR
Call INPTER
Close TERRAIN
Open SURFACE, PROFILE, CONG
Open RAWIN (if unstable hours)
Call MAP
Call SEQMOD
Stop
Figure 3-2. Outline of the main program CTDMPLUS.
3-23
-------
INPTER: Reads and writes the terrain information describing each hill. Data are passed through
"HILL.CMN".
MAP: Creates a line-printer map of the relative locations of the sources and the receptors.
The major subroutine SEQMOD contains all of the dispersion computations, and the overall struc-
ture that controls the sequential loops over receptors, hills, and sources for each stable/neutral time period.
SEQMOD calls another subroutine, DAYCALC, which performs these calculations and functions for the
unstable hours (HS < XMH and -100 < L < 0). Figure 3-3 outlines the form of SEQMOD, while Figure
3-4 outlines the form of DAYCALC.
Subroutines called by SEQMOD include the following:
RDSFC: Reads the SURFACE meteorological data for the current hour, computes w., and
passes these data through "SFCMET.CMN".
SUN: Computes the hours of sunrise and sunset.
INPEMS: Reads the hourly emissions data from the "EMISSIONS" file, and passes these data
through "STACKS.CMN".
HCRIT: Function that computes the dividing streamline height of the current hill. It calls the
following: KLOSE, GETWS, GETTA, GETDTH.
BULKFR: Function that computes the bulk Froude number for the flow above Hc. It calls
GETWS and GETTA.
GETWS: Function that returns the wind speed at a given height. It calls KLOSE.
GETWD: Function that returns the wind direction at a given height. It calls KLOSE.
GETTA: Function that returns the temperature (absolute) at a given height. It calls KLOSE.
GETDTH: Function that returns the vertical gradient of the potential temperature at a given
height. It calls KLOSE.
KLOSE: Function that returns the position of the data value (in a given array) that lies nearest to
a given value (but does not exceed the value).
SRISE: Computes height of final plume rise for the stable case (Monin-Obukhov length is posi-
tive). Calls GETWS, GETDTH, and PICK4.
3-24
-------
URISE: Computes height of final plume rise for non-stable cases (Monin-Obukhov length is
non-negative).
PICK4: Returns the minimum of four numbers.
PLAVG: Computes values of wind speed, wind direction, and potential temperature gradient at a
given height. Calls GETWS, GETWD, and GETDTH.
SIGB: Computes the size of a buoyant plume that results from mixing associated with plume
rise.
GETSW: Function returns a value of OM at a given height. It calls KLOSE.
GETSV: Function returns a value of a, at a given height. It also calls KLOSE.
PSRCE: Computes downwind and crosswind distances from the axis of the plume to a receptor.
XINTRP: Interpolation function.
MUNU: Computes the elliptical coordinates \i and v for the point ( x, y ), where the x-axis is
aligned with the major axis of the ellipse.
WRAPIN: Computes the streamfunction through the source, the stagnation streamline, the
impingement point, the distance between the plume centerline and the stagnation
streamline, and the orientation of the p-coordinate system for use in WRAP.
TERAX: Computes the major and minor axis lengths of an ellipse that forms the shape of the hill
above H c in horizontal cross-section at a given elevation.
LIFTIN: Computes the factors Ty and T z over the hill above Hc along a streamline that fol-
lows the surface of the cut-off hill if the plume centerline lies below Hc, or along a
streamline half of the way between Hc and the center-of-mass of plume material above
Hc. Calls TERAX, HILROT, and PATH.
HILROT: Computes a rotation factor and the length scales along the flow and perpendicular to
the flow of a Gaussian hill oriented at an angle to the flow.
PATH: Computes the position of a given streamline over a Gaussian hill, at a given point along
the flow, and returns the strain factors at that point. Calls FLOW and HILHGT.
3-25
-------
FLOW: Computes deflections and wind perturbations experienced by a streamline that passes
through a given point over a Gaussian hill, and computes the local strain factors at the
given point. Calls SPEED.
HILHGT: Computes the local hill height (scaled by the height at its peak) at a point on a Gaussian
hill.
FLAT: Computes the concentration at a given height above the ground for receptors that are
not influenced by terrain.
LIFT: Computes concentrations at receptors above Hc resulting from plume material that
lies above Hc. Calls FLAT (for receptors upwind of the cut-off hill), MIX, FLOW,
HILHGT, and LVDF.
MIX: Computes the depth of an internal mixing layer that develops at the hill surface above
Hc.
LVDF: Computes the vertical distribution factor for LIFT.
WRAP:- Computes concentrations resulting from plume material that lies below Hc. Receptors
upwind of the impingement point must be below Hc, but receptors downwind may be
above H c as well.
WRITIT: Writes hourly concentrations (all receptors) to a disk file in binary or text format.
TOPN: Initializes, updates, or prints out a top Af table; N is defined by the parameter MAX-
TOP in "PARAMS.INC".
Subroutines called by DAYCALC include the following:
DTHDZ: Returns vertical potential temperature gradient above the mixing height which is used in
calculating the penetration factor.
PENFCT: Function calculates the partial plume penetration factor according to Briggs (1984).
GETHIL: Calculates hill dimensions and transforms coordinates to system with hill as center and
x-axis parallel to the flow direction. These coordinates are used by the WPDF subrou-
tine.
TRANPR: This subroutine calculates transitional plume rise for receptors downwind of final rise.
WFLAT: Provides vertical velocity and probability values for a receptor not influenced by terrain.
3-26
-------
WPDF: Computes the vertical velocity needed to reach a receptor on a hill. Calls FLOWSP.
FLOWSP: Shortened version of FLOW, for use by WPDF. Computes the perturbation winds UP,
VP, WP, at the point (x,y,z) above the surface of a rotated Gaussian hill.
SIGMAY: Computes ay and the horizontal distribution function.
3-27
-------
SUBROUTINE SEQMOD
initializations
Hour Loop (top is line 100)
RDSFC [return if EOF is found]
SUN
PAGE (icase > 0)
read PROFILE data (loop over MAXLEV)
INPEMS (hourly emissions)
if ZQ, u*, ws, wd, sv, or sw is bad: skip to
write hourly met data (icase > 0)
zero concentration array
Preliminary Loop on Hills
' HCRIT
' BULKFR
End Preliminary Loop on Hills
Loop on Sources (do 300)
initialize variables for current source
if no emissions write source information: skip to
GETWS
SRISE (for stable)
URISE (for non-unstable)
if stack ht is in convective layer and L < 0, DAYCALC: skip to 300
SIGB
PLAVG
GETSW
GETSV
rotate coordinate system
compute virtual source, virtual time
screen out hills upwind of source
PSRCE
Loop on Hills (do 270)
Figure 3-3. Outline of the subroutine SEQMOD.
3-28
-------
if Hill # = 0 and no flat terrain receptors: skip to 270
if receptors are all upwind: skip to 270
. ellipse geometry for WRAP
' KLOSE
XINTRP
MUNU
' WRAPIN
. write WRAP info, (icase > 0)
. geometry for LIFT
' KLOSE
TERAX
r GETWS (compute wind shear)
' LIFTIN
. write LIFT info, (icase > 0)
. Loop on Receptors (do 260)
if receptor is not on current hill: skip to 260
initialize receptor
- if hill = 0:1**FLAT
I End Loop
set up LIFT
** LIFT (receptor > He)
_ set up WRAP
** MUNU
** WRAP
_ store hourly concentration
I End Loop on Receptors (260)
|- End Loop on Hills (270)
End Loop on Sources (300)
find maximum predicted concentration
write Max Concentration (icase > 0)
fill Top N arrays
TOPN (itopn = yes)
set all concentrations to -999 if no calculation
WRITIT (iconc > 0)
SOURCES (isor = yes)
print Top N Table
TOPN (itopn = yes)
End Loop on Hours (go to 100)
format statements
Return
Figure 3-3. (continued).
3-29
-------
SUBROUTINE DAYCALC
if plume height > mixing height, correct plume height to be .9 * mixing height
' GETWD (gets wind direction at 1/2 plume rise height)
- get rotation factors from wind direction
k GETWS (gets wind speed at 1/2 plume rise height)
' DTHDZ
' PENFCT
calculate source dependent variables
if icase > 2, write report
check hills to see if they have any downwind receptors
' PSRCE
Loop on Hills (do 270)
b
if all receptors are upwind: skip to 270
set up coordinate system for hill
Loop on Receptors (do 260)
get crosswind and*downwind distances for receptor
PSRCE
if receptor is upwind: skip to 260
rotate and translate receptor coordinates
calculate transitional plume rise, if final rise > rise
at receptor distance
TRANPR
get wind direction change over plume depth
GETWD
If hill # = 0
I** WFLAT: skip to 268
If hill # > 0
WPDF
calculate probabilities for all paths
calculate crosswind-integrated concentration
if icase > 2, write report
calculate concentration (268)
if icase > 2, write report
| End Loop on Receptors (260)
End Loop on Hill (270)
Return
Figure 3-4. Outline of the subroutine DAYCALC.
3-30
-------
SECTION 4
USER INSTRUCTIONS FOR AUXILIARY PROGRAMS
Three types of interactive computer programs that are associated with CTDMPLUS preprocessing or
postprocessing functions are described in this section:
a receptor generator program (RECGEN),
an interactive file setup program for CTDMPLUS (SETUP), and
graphical display programs for output concentrations from CTDMPLUS.
Each of these programs guides the user through execution and therefore does not require an exten-
sive amount of documentation, compared with a batch program like CTDMPLUS. General user instruc-
tions and examples of runs are given in the subsections that follow.
4.1 RECEPTOR GENERATOR
Sensitivity tests on the representation of a single terrain feature in CTDMPLUS have shown that
peak concentrations may change location (although the magnitude is not expected to change significantly)
for different mathematical formulations of hill orientation and shape. It is important, therefore, to blanket a
hill with receptors to assure that the peak concentration is captured. Obtaining x, y, and z input data for
receptors in terrain has always been a laborious procedure. Fortunately, the availability of digitized terrain
information allows the user to automate, to some extent, the selection of receptor input to CTDMPLUS.
RECGEN places receptors along digitized contours only. Since these contours are spaced to give an
adequate representation of the shape of the hill, the receptor coverage resulting from RECGEN may be
sufficient; if additonal receptors are desired, they can be added using a text editor.
The CTDMPLUS receptor generator program, RECGEN, is an interactive program written in Pas-
cal. RECGEN is designed to be compatible with the menu driver system and therefore uses the same nam-
ing conventions and directory setup as the menu driver. The user is first asked to select a plot file of terrain
contours (generated by FITCON) from a list of available files. As each contour is displayed, the user is
asked whether receptors on the contour are to have (1) a user-supplied spacing around the perimeter of the
contour, or (2) a spacing calculated by dividing the total perimeter distance of the contour by a user-supplied
number of points. The user then inputs the spacing or number of points, depending upon which option has
been selected. In general, the number of receptors on contours should be decreased as the contour perime-
ters get smaller toward the top of the hill. Once the mode of point generation has been specified for all
4-1
-------
contours, the locations of the generated receptors are plotted as filled circles upon a background of unedited
hill contours. RECGEN produces a file of receptors that is written to the 'C:\CTDM\FILES' directory with
an extension of '.RCT'.
Figure 4-1 shows a sample of a portion of the interactive session, while Figure 4-2 illustrates the dis-
play of receptors generated for an entire set of contours on a hill. An example of the resultant receptor file
is shown in Figure 4-3.
42 INTERACTIVE SETUP PROGRAM
An interactive program (SETUP) has been developed to acquaint new users with the operation of
CTDMPLUS. SETUP will create or modify files necessary to run the meteorological preprocessor, MET-
PRO, in mode 0 only. It will also create or modify files necessary to run CTDMPLUS for constant-
emission, stable/neutral cases. SETUP cannot be used to run METPRO for modes 1, 2, or 3; nor can it be
used to run CTDMPLUS for variable emissions or unstable hours. For these cases, a text editor should be
used to edit files.
SETUP first handles the input files to METPRO ("OPTIONS", "PROFILE", and "SURF1") if the
user wishes to run METPRO. For each file, the user may elect to modify an existing file, copy from an
existing file (by renaming it), or create a new file from scratch. Figures 4-4, 4-5, and 4-6 show the display of
values that can be modified for the files "OPTIONS", "PROFILE", and "SURF1", respectively. The SETUP
batch job then runs METPRO, creating an "OUTPUT" listing file and a "SURFACE" file.
The CTDMPLUS files that SETUP can modify (see Figure 4-7) include "CTDM.IN", "RECEPTOR",
"SURFACE", and "TERRAIN", respectively. SETUP cannot create a "CTDM.IN" or "TERRAIN" file from
scratch. Also, the "RAWIN" file must exist prior to starting SETUP for daytime hour runs. The SETUP
batch job will then execute CTDMPLUS, creating the output listing file "CTDM.OUT" and "CONC" (if the
user selects the concentration output). Additional runs of METPRO and CTDMPLUS can be easily made
with small alterations in the input files by using the SETUP facility.
43 GRAPHICAL CONCENTRATION DISPLAYS
For short runs of CTDMPLUS, two programs, CHIDIS and CONTOUR, are available for displaying
concentration values on a background of digitized terrain. The case-study printout from CTDMPLUS gives
all of the numerical information required, but a spatial representation of the concentration distribution is
also desirable. CHIRET, a FORTRAN program, reads the concentration and receptor information and
creates a sorted concentration file that is used by both display programs. The "CONC" file must be in text
format (created using ICONC = 2 in the "CTDM.IN" file) for the concentration display postprocessors to
work correctly.
4-2
-------
idatcorTi
Receptor spacing nethod?
1 - Distance Increnent
2 - tor of Receptors
Esc - Back
F10 - Exit to Dos
Press Enter
Figure 4-1. Portion of a sample interactive session from the RECGEN program.
4-3
-------
IdatconF
Press Enter
Figure 4-2. Receptor display from the RECGEN program.
4-4
-------
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
0.000
0.000
0.000
1.763
8.548
14.781
20.634
27.926
26.713
20.000
25.023
29.946
22.967
14.659
6.549
1.000
1.000
4.138
11.509
20.302
27359
18.432
22.067
21.422
10.545
2.000
2.069
8.247
17.081
25.253
16.429
21.810
12.737
3.000
3.000
6.143
11.589
16.832
4.000
12.401
20.802
28.763
26.905
23.890
29.545
27.074
21.238
16.972
10.733
4.108
2.148
2.610
0.516
7.000
18.279
28.569
21.327
27.302
24.000
19.726
9.933
5.859
5.182
9.000
21.275
22.507
23.081
25.000
18.288
8.190
8.912
10.000
18.655
23.714
17.470
11.168
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
10.000
10.000
10.000
10.000
10.000
10.000
10.000
10.000
10.000
10.000
10.000
10.000
10.000
10.000
10.000
20.000
20.000
20.000
20.000
20.000
20.000
20.000
20.000
20.000
20.000
30.000
30.000
30.000
30.000
30.000
30.000
30.000
30.000
40.000
40.000
40.000
40.000
40.000
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
Figure 4-3. Example receptor file generated by the RECGEN program.
4-5
-------
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
datcon
39
40
41
42
43
44
45
46
47
48
49
50
11.194
4.000
4.000
7.291
10.500
5.000
5.000
6.235
9.000
6.000
6.155
8.000
10.639
12.000
19.768
19.418
13.000
13.000
18.712
18.647
13.884
14.000
16.845
14.391
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
40.000
50.000
50.000
50.000
50.000
60.000
60.000
60.000
60.000
70.000
70.000
70.000
1
1
1
1
1
1
1
1
1
1
1
1
Figure 4-3. (continued).
4-6
-------
Modify OPTIONS file.
OPTIONS FILE - CASE STUDY MODE
Latitude: 57.000 [deg]
Longitude: 90.000 [deg]
Time Zone: 6
Surface Roughness Length: 0.30 [m]
Albedo: 0.20
Bowen Ratio: 1.00
Modify any of these values [Y/N]
Figure 4-4. Display from a sample interactive session of SETUP: "OPTIONS" input file for METPRO.
4-7
-------
Meteorological Profile Data
12/4/80
Hour: 19
Ambient Sigma
Temp Theta
[deg-K] [deg]
****** 15.30
273.71 *****
****** 35.00
273.76 14.70
Modify (M), Insert (I), Append (A), Continue (RETURN)
Enter "[M, I, A, RETURN]:
Level
#
1
2
3
4
Height
[m]
210.0
326.0
366.0
416.0
Wind
Dir
[deg]
292.0
*****
69.0
90.0
-Wind
Scalar
[m/s]
134
*****
0.13
0.40
Speed--
Vector
[m/s]
*****
*****
*****
*****
Sigma
W
[m/s]
0.01
*****
0.00
0.01
Figure 4-5. Display from a sample interactive session of SETUP: "PROFILE" input file for METPRO.
4-8
-------
TOTAL INCOMING NET MIXING CEILING CLOUD
YRMODYHR SOLAR RADIATION RADIATION HEIGHT HEIGHT COVER
[w/m2] [w/m2J [m] [ftxlOO]
8012419 210.0 0.0 292.0
Modify (M), Continue (RETURN):
Figure 4-6. Display from a sample interactive session of SETUP: "SURF1" input file for METPRO.
4-9
-------
Execute the CTDMPLUS meteorological preprocessor program, METPRO ? (Y/N):n
CTDMPLUS requires the following input files
- CTDM.IN
- RECEPTOR
- PROFILE
- SURFACE
- TERRAIN
The status of these files is:
CTDM.IN file does not exist.
NOTE: This program cannot create one from scratch.
RECEPTOR file does not exist.
NOTE: This program cannot create one from scratch.
PROFILE file does not exist
SURFACE file does not exist
TERRAIN file does not exist
NOTE: This program cannot create one from scratch.
Hit RETURN to continue
Figure 4-7. Display from a sample interactive session of SETUP: listing of files required by CTDMPLUS.
4-10
-------
43.1 CHIRET
The program CHIRET reads the "RECEPTOR" and "CONC" files, sorts the data, and writes an out-
put file ("CHIOUT"). The user is first asked to enter the identification number for the hill in question.
Although more than one hill can be input to CTDMPLUS for a single run, only concentrations at receptors
on a single hill can be displayed during a single run of either of the display programs. Next, the user must
specify which of the following three methods is to be used for selecting hours from the concentration file for
plotting concentrations:
Only the first hour is selected, or
All hours are selected, or
All hours between user-specified starting and ending times are selected.
For each of the qualifying hours in the concentration file, the concentrations and receptor coordinates
are sorted and written to the CHIRET output file, with only those receptors on the user-specified hill being
included. Figure 4-8 shows an interactive session involving execution of the CHIRET program.
432 CHIDIS
The concentration display program, CHIDIS, is a menu driven program written in Pascal. The user is
first asked to select the name of the plot file of digitized contours that was written by the program FITCON.
The user is also asked to select the name of the file of sorted concentrations generated by the program CHI-
RET. For a given hour, the receptor with the highest concentration is displayed as a blinking filled circle. If
the user presses the space bar, the receptor with the next highest concentration is displayed as a blinking
filled circle. The date, receptor number, rank, and concentration for the current receptor can be displayed
at the top of the plot by pressing the "C" key. Pressing the space bar then causes this text to be erased and
the receptor with the next lowest concentration to be displayed as a blinking filled circle. If the key "N" is
pressed at any time, the program begins the display of concentrations for the next hour. If all receptors have
been displayed for a given hour, then pressing the space bar causes the display or receptors to be repeated,
with the receptor having the highest concentration being displayed first. At any time, the user may press the
(ESC] key to terminate execution of the program. Sample screen displays from CHIDIS are shown in Fig-
ures 4-9 and 4-10.
4-11
-------
C:\CTDMS,EXE>CD\CTDM\FILES
C:\CTDU\FllES>COPY CC8.CON CONC
1 File(s) copied
C:\CTDM\FILES>COPY CCB.RCT RECEPTOR
1 File(s) copied
I
C:\CTDM\FILES>C:\CTDM\EXE\CHI RET
ENTER THE HILL IDENTIFICATION NUMBER(1-9?) -> 1
70 OUT OF 70 RECEPTORS WERE FOUND TO BE ON HILL 1
SPECIFY CONCENTRATION RECORD SELECTION MODE
I.) ONLY THE FIRST RECORD IS SELECTED
2.) ALL RECORDS SELECTED
3.) ALL RECORDS SELECTED BETWEEN A START AND END TIME
CHOICE?(1,2 OR 3) -> I
Figure 4-8. Sample Interactive session for the program CHIRET.
4-12
-------
Spc - Next Receptor
C - Current Cone
N - Next Hour
HO - Exit to DOS
Esc - Back
"*"<.,
Figure 4-9. Sample screen display from the CHIDIS program, showing user options.
4-13
-------
BO/ 5/26: 3 Rank 7 SMu* 23 CO.CM58
Figure 4-10. Sample screen display from the CHIDIS program, showing receptor locations and the
concentration value for the most recently displayed receptor.
4-14
-------
4.33 CONTOUR
The second concentration display program, CONTOUR, can only be run from the menu driver and is
described in greater detail in Section 5.2.8. CONTOUR allows the user to see concentration isopleths
drawn on a background of unedited hill contours. The CONTOUR program is based on the NCAR plotting
routines (McArthur, 1983) and uses the Barnes (1973) scheme for gridding the data.
4-15
-------
-------
SECTION 5
INSTRUCTIONS FOR USING THE MENU DRIVER
5.1 INTRODUCTION
The CTDMPLUS modeling system involves the use of a number of programs and input files. The
required input files are discussed in detail in Section 3. To facilitate use of the model, a menu driver was
developed that helps users select programs to be run and options and files to be used. The main compo-
nents of the menu driver are shown in Table 5-1. The menu system is designed to aid users with short
CTDMPLUS runs of a few hours or days. Although CTDMPLUS can be run on a mainframe or a PC, the
menu system is designed to work on a PC; using the menu system and CTDMPLUS together requires a PC
with 640K bytes of memory. Users must be familiar with the files required by CTDMPLUS, as well as the
fundamentals of model operation, in order to use the menu system. (Reviewing the terrain and meteorolog-
ical preprocessor manuals and previous sections of this manual will provide this information.) The require-
ments for the naming and directory location of files used by the menu driver are described in Sections 5.2.1
and 5.2.2.
52 USER INSTRUCTIONS
52.1 File Naming
A standard method for naming files has been developed for use with the menu driver. File names
consist of a descriptive name and a required extension. For example, the "PROFILE" file for use in model-
ing the CCB dataset would be named "CCB.PFL". The selection of the descriptive term (i.e., CCB) is left to
the user, but the file name extensions must be as indicated in Table 5-2 in order to use the menu driver.
During the operation of the menu driver, when a user is asked to input a filename, only the descriptive name
should be entered as the menu driver will attach the appropriate extension and pathname.
522 Directory Setup
The menu system is designed to work with a specific directory structure, which the user must create
on the hard disk. The files must be named as indicated in section 5.2.1 and in the directories as listed in
Figure 5-1. Files not properly named and located in the proper directory will not be available for use with
the menu system.
5-1
-------
TABLE 5-1. COMPONENTS IN THE CTDMPLUS MENU DRIVER
Main program
Subprograms
Terrain preprocessor
Meteorological preprocessors
CTDMPLUS
Concentration postprocessors
FITCON
HCRIT
PLOTCON
RECGEN
METPRO
READ62
CHIDIS
CONTOUR
5-2
-------
TABLE 5-2. STANDARD CTDMPLUS FILE EXTENSIONS
Description
FITCON options file
Contour master file
FITCON diagnostic output
FITCON plot file
FITCON output file
HCRIT options file
HCRIT output file
HCRIT plot file
RECGEN receptor file
METPRO OPTIONS file
SURF1 file
SURF2 file
RAWIN file
METPRO output file
RAWIN options file
TD6201 file
RAWIN diagnostic. output
CTDM.IN file
PROFILE file
SURFACE file
EMISSION file
TERRAIN file
RECEPTOR file
CONC file
CTDM.OUT file
CHIRET output file
Qridded data file
Required Name
FOPTIONS
HOPTIONS
OPTIONS
SURF1
SURF2
RAWIN
OUTPUT
OPT62
TD6201
R62OUT
CTDM.IN
PROFILE
SURFACE
EMISSION
TERRAIN
RECEPTOR
CONC
CTDM.OUT
CHIOUT
Extension
.OPT
.XY
.DAG
.PLT
.FO
.HOP
.HCO
.HPT
.RCT
.MOP
.SF1
.SF2
.RAW
.MOF
.R62
.TD
.ROF
.CIN
.PFL
.SFC
.EMS
.HCO
.RCT
.CON
.OUT
.CHI
.GRD
5-3
-------
MENU.BAT
\CTDM
- \EXE\
\TERRAIN\
\FILES\
F1TCON.EXE, HCRIT.EXE, PLOTCON.EXE, RECGEN.EXE,
METPRO.EXE, READ62.EXE, CTDMPLUS.EXE,
CHIRET.EXE, CHIDIS.EXE, DRIVEIT.EXE
*.XY, *.OPT, *.DAG, *.PLT, *.FO, *.HOP,
*.HPT
*.CIN, *.PFL, *.SFC, *.EMS, *.HCO, *.RCT,
*.OUT, *.CON, *.CHI, *.MOP, *.SF1, *.SF2,
*.RAW, *.MOF, MD, *.R62, *.ROF, *.GRD
Figure 5-1. Directory setup for menu driver use.
5-4
-------
53 Initiating the Program
The menu driver can be initiated from any directory by typing MENU. This will run the MENU
batch file located in the root directory. The batch file changes the current directory to \CTDM\EXE\ and
executes the menu driver program. DRIVEIT. The main selection menu shown in Figure 5-2 will then be
displayed.
The menu system does not allow for editing files, so all files to be used must be edited and exist in the
proper directories as listed in Figure 5-1 prior to initiating the menu system.
Certain special keys are used in the menu system:
[ESC] - back up one screen
fFiol - end the program and exit to DOS
fF9] - execute a particular FORTRAN program
(HD - end selection of files
dZ) - overwrite an existing file.
In several instances, the user is given the option to view a particular file. This is done using the DOS
command
TYPE < filename > j MORE
which enables the user to see the file one page at a time. Pressing any key will display the next page, while
Ctrl-C will return the user to the menu system. A response of "N" to the view option will also return the
user to the menu system.
5.2.4 Terrain Preprocessor Programs
The terrain preprocessor menu is displayed when the user chooses option 1 from the main menu.
The user may then choose to execute one of four programs: FTTCON, HCRIT, PLOTCON, or RECGEN.
Note that FITCON and HCRIT both require an options file. Previous versions ran interactively, but this
was not practical for multiple runs with the same hill.
The FITCON program is initiated by selecting option 1 from the terrain preprocessor menu (see Fig-
ure 5-3). The user is prompted to choose a contour master file (filenameJCY) from a list of available files.
The filename of the contour master file is used to determine if a related FITCON options file exists
(filename.OPT). If the options file exists, the parameters listed in it are used to initialize the values for the
next run, otherwise, default values are used. The options file is displayed as shown in Figure 5-4. The user
is then given the opportunity to change any of the parameters in the options file. Validation of each field
(e.g., range checking, checking for file existence) occurs when( Enter ] is pressed and errors are
5-5
-------
02/09/1989 CTDMPLUS - SELECTION MENU 07:47:10
Press a number to make a selection.
1 - Terrain Preprocessor
* 2 Meteorological Preprocessor
3 - CTDMPLUS
4 - Concentration Postprocessors
Esc - Baclc F10 - Exit to Dos
Figure 5-2. Screen display from the menu driver system, showing the main selection menu.
J
5-6
-------
01/02/1989 TERRAIN PREPROCESSOR SELECTION MENU 12:21:10
Press a number to execute the required program.
»
I - FtTCON
2 - HCRIT
*
3 - PLOTCON
4 - RECCEN
Esc - Back F10 - Exit to Dos
Figure 5-3. Screen display from the menu driver system, showing the terrain preprocessor menu.
5-7
-------
01/02/1989
FITCQN OPTIONS
12:21 :58
Use the arrow keys to hightiight the option to be changed.
CONTOUR MASTER FILE {.XY}:
DIAG. OUTPUT FILE (.DAG):
HILL ID NUMBER:
HILL NAME:
HILL TOP ELEVATION:
Hill CENTER X-COORDINATE.:
HILL CENTER Y-CQORDINA^?:
ANGULAR FILTERING (Y or N):
ANGULAR FILTER SIZE (1-22,5):
CONTOUR SELECTION MODE (1,2, OR 3)i
LOWER BOUND OF ID NUMBERS (MODE 2 ONLY):
UPPER BOUND OF ID NUMBERS (MODE 2 ONLY):
CONTOUR ID NUMBER FILE (MODE 3 ONLY):
CREATE A PLOT TILE? ? (Y or N}:
FITCON PLOT FILE(.PLT):
FITTED OUTPUT HLE(.FO):
F9 - Run FfTCON
DOWN
UP
Esc - Back
DATCON
TcO
1
DATCON
80.0000
7.0000
14.5000
n
I .0
1
I
1
DATCON
Y
tcO
tcO
FIO - Exit to Dos
Figure 5-4. Sample screen display from the menu driver system, showing the FITCON options screen.
5-8
-------
displayed at the bottom of the screen. Pressing the (FV) key causes each field to be validated a final time
and FITCON to be run. FITCON cannot be executed while errors exist in the options. After FITCON
completes execution, the user is given the option to look at the diagnostic output file (filename.DAG). The
fitted ellipses from the FITCON run can be viewed graphically using PLOTCON, or the user may elect to
continue on with HCRTT.
The HCRIT program is initiated by selecting option 2 from the terrain preprocessor menu. The user
is then prompted to choose a FITCON output file (.FO) from a list of available files. The cursor will ini-
tially be on the file that matches the .XY file chosen hi the last FITCON run, if one exists. The filename of
the .FO file is used to determine if a related HCRIT options file (.HOP) exists. From this point, the
sequence of events follows that listed above for FITCON. A sample HCRIT options file is shown in Figure
5-5. After the HCRIT run finishes, the user is prompted to press ENTER to return to the terrain menu.
The cut-off contours may be displayed graphically by using PLOTCON. Note that an HCRIT run will
create the CTDMPLUS terrain file (.HCO) for only one hill and the file will be placed in the \CTDM\FI-
LES\ directory. Files must be combined (outside of the menu driver) for CTDMPLUS runs involving mul-
tiple hills.
The PLOTCON program is activated by choosing option 3 from the terrain preprocessor menu. The
user is prompted to select a FITCON plot file (.PLT) from a list of available files. The next screen gives the
user the option to see edited or unedited contours. After the contours have been displayed, the user has the
option of seeing the fitted ellipses (from FITCON) or cut-off hill contours (from HCRIT). After displaying
the ellipses, the user is again given the option to see the cut-off contours. For viewing the cut-off contours,
the user chooses an HCRIT plot file (.HPT) and the starting elevation for displaying the cut-off contours.
After one set of cut-off contours has been displayed, the user may press [ESC) to return to the starting eleva-
tion menu, or I Enter 1 to see the next set of contours.
5.2.5 Receptor Generator
The receptor generator program, RECGEN, is initiated by choosing option 4 of the terrain prepro-
cessor menu. The operation of RECGEN is described in detail in Section 4.1.
52.6 Meteorological Preprocessors
The meteorological preprocessor menu, shown in Figure 5-6, is displayed by selecting option 2 from
the main menu. The user, at this point, may elect to run METPRO with existing files or create a new "RA-
WIN" file using the READ62 preprocessor.
5-9
-------
01/02/1989 . HCRIT OPTIONS . . 12:22:24
Use the arrow keys to highlight the option to be changed.
FITCON OUTPUT FILE (,FO): J_CO
HCRIT OUTPUT FILE {.HCO}: TOO
CREATE A PLOT FILE? (Y or N): Y
HCRIT PLOT FILE (.HPT): TCO
SELECTION MODE FOR CRIT. ELEV,: . I
NUMBER OF CRIT. ELEV. (MODE 2 ONLY): 1
LOWEST CRIT. ELEV. (MODE 2 ONLY): 0.0000
F9 - Run HCRIT DOW UP Esc - Back FIO - Exit to Dos
Figure 5-5. Sample screen display from the menu driver system, showing the HCRIT options screen.
5-10
-------
01/.02/1989 METEOROLOGICAL PREPROCESSOR MENU 12:22:44
Press a number to execute * the required program.
! - Rua METPRO
2 - Create a RAW1N File
Esc - Back FLO - Exit to Dos
y
Figure 5-6. Screen display from the menu driver system, showing the meteorological preprocessor
menu.
5-11
-------
The execution of METPRO is initiated by choosing option 1 from the meteorological preprocessor
menu. In running METPRO, the user is first prompted to choose an "OPTIONS" file (.MOP). A sample of
the screen display for selecting files is shown in Figure 5-7. The mode switch in the options file is read and
selection of files for "PROFILE", "SURF1", "SURF2", and "RAWIN" continues as needed for each mode.
After file selection has been completed, the chosen files are displayed for final approval. Files that are not
required for a particular mode have an "N/A" displayed in the area for the filename. If a filename needs to
be changed, the number corresponding to the file is pressed and the selection menu for that file is dis-
played. The user may end selection at any time by pressing iZE . The dE key is pressed to execute
METPRO with the listed files. An error message will be displayed if files needed for execution have not
been selected. The METPRO program will create a "SURFACE" file with the same filename as the options
file and an extension of .SFC, which will be placed in the \CTDM\FILES\ directory. If this file already
exists, a warning message will be displayed. The user will have the option of overwriting the existing file or
creating a backup copy of the old file (.&SF). After execution, the user is given the option to look at the
diagnostic output (.MOF) file or to return to the meteorological preprocessor menu.
File selection for the READ62 preprocessor proceeds in similar fashion to that for METPRO. The
READ62 program creates the "RAWIN" file (.RAW) needed for METPRO mode 3 and for daytime hours
in CTDMPLUS. After execution of the program, the user may look at the diagnostic output file (.R62) or
return to the meteorological preprocessor menu.
52.7 Running CTDMPLUS
Choosing option 3 from the main menu (Figure 5-2) will start the file selection process for
CTDMPLUS. The CTDM.IN file is chosen first and is read to determine which optional files are needed.
The user is prompted to select the "PROFILE", "SURFACE", "RECEPTOR", TERRAIN", and, if neces-
sary, the "EMISSION" and "RAWIN" files. The function keys, (ED and ED , operate in the same manner
as for METPRO. An example of the screen showing the selected files is given in Figure 5-8. The execution
of CTDMPLUS will produce two output files, "CONC" and "CTDM.OUT", which are then renamed to file-
name .CON and filename .OUT, respectively. Upon completion of the run, the user may view the diagnos-
tic output file or return to the main menu.
5-12
-------
01/02/1989 , yETPRO 12:37:15
Please highlight a PROFILE fite, then press ENTER:
_TC&.PFL
"TCI,PFL
* TC2.PFL
TC3.PFL
F8 - End Selection Esc - Back F10 - Exit to Dos
Figure 5-7. Sample screen display from the menu driver system, showing the selection of a file for use
with METPRO.
5-13
-------
02/09/1989 CTDMPLUS . 12:29:45 j
I
I
i
The following files have beea selected. Type a !
number to change a file name or F9 to run CTDMPLUS.
I
2
3
4
5
6
7
_Esc
V
- CTDM. IN -
- EMISSION -
- SURFACE -
- PROFILE -
- TERRAIN -
- RECEPTOR -
- RAWIN
- Back FtO -
TCO.CIN
N/A
TCO.SFC
TCO.PFL
TCO.HCO
TCO.RCT
N/A
Exit to Dos
Figure 5-8. Sample screen display from the menu driver system, showing the list of files selected
for use with CTDMPLUS.
5-14
-------
5.2.8 Concentration Postprocessors
The concentration postprocessor menu is displayed by selecting option 4 from the main menu (Figure
5-2). The user may elect to run CHIRET, CHIDIS or CONTOUR (see Figure 5-9). CHIRET is a FOR-
TRAN program that reads the concentration and receptor information and creates an output file of sorted
concentrations that is used by both CHIDIS and CONTOUR. CHIDIS, a program that displays
concentration values at the receptors in decreasing order is described in Section 4.3.1. CONTOUR, which
displays concentration isopleths, is described below.
If CHIRET is chosen from the menu, the user is asked to select a concentration file (*.CON) from a
list of available files. Next, the corresponding receptor file is chosen. Then, the user must enter the name to
be given to the output file (*.CHI). A batch job will be run to execute CHIRET. The execution of CHI-
RET continues as described in section 4.3.1. After CHIRET has been run, the user is returned to the CON-
TOUR menu. When the CONTOUR program is selected from the concentration postprocessor menu, the
user is given the option to grid the data or run CONTOUR. If the gridding routine is run, the user will be
prompted to select a ".CHI" file and provide a name for the output file of gridded data (*.GRD). Upon
completion, the gridding routine creates the output file which contains the data in gridded form, ready for
contouring.
To display the concentration isopleths, the user must select a ".GRD" file from the list of available
files. Next, the user decides whether to display the hill contours in addition to concentration contour. If the
hill is to be plotted, the user must select, from a list of available files, the FITCON plot file that contains the
hill information. Then, the user is asked to indicate the method for determining the contour interval. The
user may choose default values or enter the start, stop, and contour interval values. Finally, the isopleths are
displayed, on a background of hill contours if that option was selected. See Figure 5-10 for a sample display.
Note that CONTOUR will only provide the user with meaningful results if the receptors are widely scat-
tered over the hill.
5-15
-------
01/02/1989 CONCENTRATION POSTPROCiSSOR MENU 12:23:49
Press a aumber to execute the required program.
1 - CHIRIT
2 - CHIBIS
3 - CONTOUR
Esc - Back FIO - Exit to Dos
Figure 5-9. Screen display from the menu driver system, showing the concentration postprocessor
menu.
5-16
-------
CCB
Figure 5-10. Sample screen display from the CONTOUR program, showing the concentration iso-
pleths (solid lines) and hill contours (dashed lines).
5-17
-------
REFERENCES
Agterberg, R., F. T. M. Nieuwstadt, H. van Duuren, A. J. Hasselton, and G. D. Krijt, 1983: Dispersion
experiments with sulfur hexafluoride from the 213 m high mast at Cabauw in the Netherlands. Royal
Netherlands Meteorological Institute, DeBilt.
Barad, M. L., Ed., 1958: Project Prairie Grass: a field program in diffusion. Geophysical Research Paper
No. 59, Vols. I and H, AFCRF-TR-58-235, Air Force Cambridge Research Center, Bedford, MA.
Barnes, S.L., 1973: Mesoscale objective map analysis using weighted time series analysis. National Oceanic
and Atmospheric Administration Technical Memorandum ERL NSSL-62. National Severe Storms Labo-
ratory, Norman, OK.
Batchelor, G. K., 1970: An Introduction to Fluid Dynamics. Cambridge University Press, London NW1.
Bjorklund, J. R., and J. F. Bowers, 1982: User's instructions for the SHORTZ and LONGZ computer pro-
grams. Vol. I. EPA-903/9-82-004a, U.S. Environmental Protection Agency, Region HI, Philadelphia,
PA.
Briggs, G. A., 1973: Diffusion estimation for small emissions. ATDL Contribution File No. 79, Atmo-
spheric Turbulence and Diffusion Laboratory, U.S. Environmental Protection Agency, Research Triangle
Park, NC. -
Briggs, G. A., 1975: Plume rise predictions. In Lectures on Air Pollution and Environmental Impact Analy-
ses. American Meteorological Society, Boston, MA.
Briggs, G. A., 1984: Plume rise and buoyancy effects. Chapter 8 in Atmospheric Science and Power Produc-
tion, D. Randerson, Ed. DOE/TIC-27601, U.S. Department of Energy.
Briggs, G A., 1985: Analytical parameterization of diffusion: the convective boundary layer. /. Climate
Appl. Meteor., 24,1167-1186.
Briggs, G. A., W. L. Eberhard, J. E. Gaynor, W. R. Moninger, and T. Uttal, 1986: Convective diffusion
measurements compared with laboratory and numerical experiments. Preprint, Fifth Joint Conference
on Applications of Air Pollution Meteorology, American Meteorological Society, Chapel Hill, NC,
340-343.
Briggs, G. A., and K. R. McDonald, 1978: Prairie Grass revisited: optimum indicators of vertical spread.
Proceedings of the Ninth International Technical Meeting on Air Pollution Modeling and its Applications,
No. 103, NATO, 209-220.
Brighton, P. W. M., 1978: Strongly stratified flow past three-dimensional obstacles. Quart. J. Roy. Meteor.
Soc., 104,289-307.
Burt, E. W., 1977: Valley model user's guide. EPA-450/2-77-018, U.S. Environmental Protection Agency,
Office of Air Quality Planning and Standards, Research Triangle Park, NC.
Ref-1
-------
Businger, J. A., 1973: Chapter 2: Turbulent transfer in the atmospheric surface layer. In Workshop on
Micrometeorology, D. A. Haugen, Ed. American Meteorological Society, Boston, MA.
Carson, D. J., 1973: The development of a dry inversion-capped convectively unstable boundary layer.
Quart. J. Roy. Meteor. Soc., 99, 450-467.
Clark, R. N., D. G. Reid, and A. J. Troup, 1971: The Wangara experiment: boundary layer data. Technical
Paper 19, Division of Meteorology and Physics, Commonwealth Scientific and Industrial Research Orga-
nization, Australia.
Crapper, G. D., 1959: A three-dimensional solution for waves in the lee of mountains. J. Fluid Mech., 6,
51-76.
Csanady, G. T., 1974: Equilibrium theory of the planetary boundary layer with an inversion lid. Boundary
Layer Meteor., 6,63-79.
Drazin, P. G., 1961: On the steady flow of a fluid of variable density past an obstacle. Tellus, 13,239-251.
Colder, D., 1972: Relations among stability parameters in the surface layer. Boundary Layer Meteor., 3,
47-58.
Guldberg, P. H., J. P. Myers, K. W. Wiltsee, and P. Morgenstern, 1977: Handbook for the single source
(CRSTER) model. EPA-450/2-77-013, U.S. Environmental Protection Agency, Office of Research and
Development, Research Triangle Park, NC.
Hanna, S. R., 1983: Lateral turbulence intensity and plume meandering during stable conditions. /. Climate
Appl. Meteor., 22,1424-1430.
Hanna, S. R., 1986: Lateral dispersion from tall stacks. /. Climate Appl. Meteor., 25,1426-1433.
Hanna, S. R., and R. J. Paine, 1987: Convective scaling applied to diffusion of buoyant plumes from tall
stacks. Atmos. Environ., 21,2153-2162.
Hanna, S. R., J. C. Weil, and R. J. Paine, 1986: Plume model development and evaluation-hybrid approach.
EPRI Contract No. RP-1616-27, Electric Power Research Institute, Palo Alto, CA.
Hovind, E. L., M. W. Edelstein, and V. C. Sutherland, 1979: Workshop on atmospheric dispersion models
in complex terrain. EPA-600/9-79-041, U. S. Environmental Protection Agency, Research Triangle
Park, NC.
Hunt, J. C. R., and R. J. Mulherarn, 1973: Turbulent dispersion from sources near two-dimensional
obstacles. J. Fluid Mech., 61,245-274.
Hunt, J. C. R., and W. H. Snyder, 1980: Experiments on stably and neutrally stratified flow over a model
three-dimensional hill. /. Fluid Mech., 96,671-704.
Johnson, W. B., and E. E. Uthe, 1971: Lidar study of the Keystone stack plume. Atmos. Environ., 5,
703-724.
Ref-2
-------
Kaimai, J. C., W. L. Eberhard, W. L. Moninger, J. E. Gaynor, S. W. Troxel, T. UttaL G. A. Briggs, and G. E.
Start, 1986: Project CONDORS: convective diffusion observed by remote sensors. Report No. 7,
NOAA/ERL/Boulder Atmospheric Observatory, Boulder, CO.
Kato, H., and O. M. Phillips, 1969: On the penetration of a turbulent layer into stratified fluid. /. Fluid
Mecii., 37,643-655.
Lamb, R. G., 1979: The effect of release height on material dispersion in the convective boundary layer.
Preprints, Fourth Symposium on Turbulence, Diffusion, and Air Pollution, American Meteorological
Society, Reno, NV.
Lamb, R. G., 1981: Diffusion in the convective boundary layer. In A Short Course on Atmospheric Turbu-
lence and Air Pollution Modeling, F. T. M. Nieuwstadt and H. van Dop, Eds. D. Reidel Publishing Com-
pany, Boston, MA.
Li, Z.-K., and G. A. Briggs, 1988: Simple PDF models for convectively driven vertical diffusion. Atmos.
Environ., 22, 54-74.
Lumley, J. L., and H. A. Panofsky, 1964: The Structure of Atmospheric Turbulence. John Wiley and Sons,
New York, NY.
McArthur, G. R., 1983: An introduction to the SCO graphics system. NCAR Technical Note NCAR/TN-
161 + 1A, National Center for Atmospheric Research, Boulder, CO.
Megarejo, J. W., and J. W. Deardorff, 1974: Stability functions for the boundary-layer resistance laws based
upon observed boundary layer heights. /. Atmos. Sci., 31,1324-1333.
Mills, M. T., R. J. Paine, E. M. Insley, and B. A. Egan, 1987: The complex terrain dispersion model
(CTDM) terrain preprocessor systemuser guide and program descriptions. U.S. Environmental Protec-
tion Agency, Atmospheric Sciences Research Laboratory, Research Triangle Park, NC.
Misra, P. K, 1982: Dispersion of non-buoyant particles inside a convective boundray layer. Atmos. Emi-
ron., 16,239-243.
Nieuwstadt, F. T. M., 1981: The steady-state height and resistance laws of the nocturnal boundary layer:
theory compared with Cabauw observations. Boundary Layer Meteor., 20,3-17.
Paine, R. J., 1987: User's guide to the CTDM meteorological preprocessor (METPRO) program. U.S.
Environmental Protection Agency, Atmospheric Sciences Research Laboratory, Research Triangle Park,
NC.
Paine, R. J., and B. A. Egan, 1987: User's guide to the rough terrain diffusion model (RTDM)-revision
3.20. ERT Document PD-535-585, ERT, Inc., Concord, MA.
Paine, R. J., D. G. Strimaitis, M. G. Dennis, R. J. Yamartino, M. T. Mills, and E. M. Insley, 1987: User's
guide to the complex terrain dispersion model. Vol. 1. EPA-600/8-87-058a, U.S. Environmental Protec-
tion Agency, Research Triangle Park, NC.
Ref-3
-------
Pasquill, F., 1976: Atmospheric dispersion parameters in Gaussian plume modeling. Part IIPossible
requirements for change in the Turner workbook values. EPA-600/4-76-030b, U.S. Environmental Pro-
tection Agency, Research Triangle Park, NC.
Perry, S. G., and W. H. Snyder, 1989: Laboratory simulations of daytime mixed layer heights in complex
terrain. Accepted for presentation at the Sixth Joint Conference on Applications of Air Pollution Meteo-
rology, American Meteorological Society, Anaheim. CA.
Riley, J. J., H. T. Liu, and E. W. Geller, 1976: A numerical and experimental study of stably stratified flow
around complex terrain. EPA-600/4-76-021, U.S. Environmental Protection Agency, Research Triangle
Park,NC.
Sheppard, P. A., 1956: Airflow over mountains. Quart. J. Ray. Meteor. Soc., 82,528-529.
Smith, R. B., 1980: Linear theory of stratified hydrostatic flow past an isolated mountain. Tellus, 32,
348-364.
Snyder, W. H., and J. C. R. Hunt, 1984: Turbulent diffusion from a point source in stratified and neutral
flows around a three-dimensional hill. Part II~Laboratory measurement of surface concentrations.
Atmos. Environ., 18,1969-2002.
Snyder, W. H., R. S. Thompson, and M. S. Shipman, 1986: Streamline trajectories in neutral and stratified
flow over a three-dimensional hill. Appendix to EPA Complex Terrain Model Development: Fifth Mile-
stone Report. EPA-600/3-85-069, U.S. Environmental Protection Agency, Research Triangle Park, NC.
Taylor, G. I., 1921: Diffusion by continuous movements. Proc. London Math. Soc. Ser. 2., 20, p. 196.
Tennekes, H., 1981: Similarity relations scaling laws and spectral dynamics. In Atmospheric Turbulence and
Air Pollution Modeling, F. T. M. Nieuwstadt and H. van Dop, Eds. D. Reidel Publishing Company, Bos-
ton, MA.
U.S. Environmental Protection Agency, 1986: Guideline on air quality models (revised). EPA-450/2-78-
027R, Office of Air Quality Planning and Standards, Research Triangle Park, NC.
Venkatram, A., 1983: On dispersion in the convective boundary layer. Atmos. Environ., 17,529-533.
Venkatram, A., D. Strimaitis, and D. DiCristofaro, 1984: A semi-empirical model to estimate vertical dis-
persion of elevated releases in the stable boundary layer. Atmos. Environ., 18, 923-928.
Venkatram, A., and R. Vet, 1981: Modeling of dispersion from tall stacks. Atmos. Environ., 15,1531-1538.
Weil, J. C., 1977: Evaluation of the Gaussian plume model at Maryland power plants. Reference No.
PPSP-MP-16, Maryland Power Plant Siting Program, Martin Marietta Corporation, Baltimore, MD.
Weil, J. C., 1979: Assessment of plume rise and dispersion models using lidar data. Reference No. PPSP-
MD-24, Maryland Power Plant Siting Program, Martin Marietta Corporation, Baltimore, MD.
Weil, J. C., and R. P. Brower, 1983: Estimating convective boundary layer parameters for diffusion applica-
tions. Martin Marietta Corporation for the Maryland Department of Natural Resources, Baltimore, MD.
Ref-4
-------
Weil, J. C, and W. F. Furth, 1981: A simplified numerical model of dispersion from elevated sources in the
convective boundary layer. Extended abstract, Fifth Symposium on Turbulence, Diffusion, and Air Pollu-
tion, American Meteorological Society, Atlanta, GA, 76-77.
Willis, G. E., and J. W. Deardorff, 1976: A laboratory model of diffusion into a convective planetary bound-
ary layer. Quart. J. Roy. Meteor. Soc., 102,427-445.
Willis, G. E., and J. W. Deardorff, 1978: A laboratory study of dispersion from an elevated source in a con-
vective mixed layer. Atmos. Environ., 12,1305-1313.
Willis, G. E., and J. W. Deardorff, 1981: A laboratory study of dispersion in the middle of the convectively
mixed layer. Atmos. Environ., 15,109-117.
Yamartino, R. J., 1984: A comparison of several "single pass" estimators of the standard deviation of wind
direction. /. dim. Appl. Meteor., 23,1362-1366.
Yamartino, R. J., 1987: Appendix A~Exact solutions to the linearized equation for stratified flow over ter-
rain in multidimensional space. In EPA Complex Terrain Model Development Final Report. EPA-600/3-
88/006, U.S. Environmental Protection Agency, Research Triangle Park, NC.
Ref-5
-------
APPENDIX A
DETAILS OF LIFT AND WRAP ALGORITHMS
A.1 TERRAIN FACTORS Tz and Ty
For axi-symmetric strain, the theory of Hunt and Mulhearn (1973) shows that
2 _ 2 /"' 2
|32(O Jo
where a denotes either a y or a z, (3 is the strain function for the corresponding component, and K. is the
diffusivity for that component. In the case of a2, the strain function (3 - for small deflections is approxi-
mately
|3Z(O = e'^'*"" (A-2)
where 6z(f) is the vertical deflection experienced by a streamline at time t.
The derivation of Equation A-l depends on a local analysis of the flow field near the centerline of the
plume. The strain function is treated as a constant about the centerline, varying only with distance (time)
along the flow. In using this result, we assume that the strain function can be assigned a representative value
for the layer that contains the bulk of the plume material.
Let Th(zm, t), the factor for distortion streamlines in the vertical, be the ratio of the spacing of
streamlines in the vertical in the strained flow to that in the incident flow, evaluated at time t along a
streamline whose height far from the hill is zm . The streamline of height zm is chosen to represent the layer
in which the plume resides. Then the distortion factor for the layer at time t is defined as:
With this mean distortion, the height r\(z,f) of a streamline above the surface of the hill at time t is
approximately
A-l
-------
(A-4)
where z is the height of streamline well upwind of the hill. Because of the form of Equation A-4, 7*h is
equivalent to a "height-correction factor" for a plume.
To evaluate Tz , consider the argument of the exponential function in the vertical distribution factor
of Equation 40 (Section 2.5) when the elevation of the receptor above the surface is zero. The primary
quantity hi this argument is the ratio zr/a 2. , where a zt is given in Equation 38 (Section 2.5). The corre-
sponding ratio in the approach of Hunt and Mulhearn is nCZr.O/^CO where n ( z r , f) is given by Equa-
tion A-4, and az(() is given by Equation A-l. Tz is evaluated by equating these two quotients.
Using Equations A-2 and A-3, the strain function can be rewritten in terms of Th as:
e"(T"'U . (A-5)
For weak strain, Th is of the order 1, so that PZ(O = l/rA(zm,O. Adopting the weak-strain form of
the product P27"A, the expression for 7Z becomes:
1 Jt,
fi ^r. (A-6)
Note that the denominator is the difference of the squares of the plume size across the interval t -10 in the
absence of any effects of terrain.
This expression is implemented in CTDMPLUS by breaking the integral in 25 subintervals along the
streamline trajectory over the hill. Within each subinterval, the strain function is assumed to be constant,
equal to its value at the midpoint of the subinterval. This allows the integral for one of the subintervals to be
written as:
r
A-2
-------
where j32 is obtained from Equation A-5. The altered value of a, (aia)is given by Equation 12 (Section
2.5), except o and TL have been altered by the changes in the flow over the hill.
The Lagrangian time-scale in the absence of any strain in the flow is given by the relation
TL. i2 rz,
(A-8)
where the subscripts 'o' denote the flat-terrain quantities. The strain in the flow alters the stratification
(temperature gradient), the turbulence velocity, and the length scale of the dominant eddies, so that
_L = N"^2 + ^L T°" (A
TL V7"ft(*/n.O r~* 7"/.(z/n-0
where Taa is the ratio oa,/aa>0. The value of Taa is estimated from inner-layer theory (Britter et al.
1981) so that Tau = Tu, the speedup factor for the flow at time t. This is equivalent to saying that the
turbulence intensity (iz) is not altered by the flow over the hill.
A parallel development factor for Ty leads to a similar result. The distortion factor in the horizontal
direction is denoted by 7",( z m, 1) so that
(A-10)
This integral is also broken into subintervals, as in Equation A-7. In this case though, a ya ( O is given by
the linear growth form:
= Tav avot (A-ll)
because this form holds for distances of travel of nearly 10 km, T is set to unity.
A-3
-------
AJ EFFECTIVE RECEPTOR POSITION AND THE DISTORTION FACTORS (LIFT)
The receptor position in Equation 39 (Section 2.5) is denoted as (ft,', y, '), where the primes indi-
cate that the distance of the receptor from the centerline of the plume has been altered by the presence of
the hill. The property of the flow that is conserved as the flow is deformed over the hill is the number of
streamlines. Upwind of the hill, streamlines are straight and parallel. Over the hill, the spacing is variable
and the streamlines suffer deflections. Because concentrations are computed for a parallel flow in which the
influence of terrain is manifest in the Hc -partition of plume material, in the altered rates of diffusion, and
in the altered streamline that passes through the receptor, the chief task in computing (ft *', y *') is in
identifying the streamline that passes through the source, and the streamline that passes through the recep-
tor. The position of these two streamlines in the unperturbed flow upwind of the hill defines the effective
receptor location relative to the centerline of the plume. These streamlines are obtained from the flow
model described in Section A-4.
The design of the flow model is particularly well-suited to obtaining ft,' and yx ' because it is for-
mulated as a "backwards-looking" solution. It is designed to compute the deflection experienced by a
streamline that passes through a given point over the hill. If the given point is a receptor, then the model
will compute the deflection experienced by the streamline that passes through that receptor. Knowing the
actual position of the receptor, the deflections allow ft,' and y,' to be computed.
The deformation factors Th and T\ and the speedup factor Tu are also obtained from the model.
These are computed along the path of the streamline chosen to be representative of the layer containing the
plume material. Far from the hill, this streamline is taken to be half of the way between Hc and the eleva-
tion of the center of mass of plume material above Hc .
Because a particular streamline must be followed along the flow ( a "forward-looking" process), its
position at any time t must be found by an iterative process. Each "guess" at where the streamline may be at
time t produces a computed position in the incident flow. New guesses are dependent on the difference
between the last computed position, and the position of the targeted streamline.
A3 INTERNAL MIXING LAYER (LIFT)
A central assumption involved in the use of the LIFT equations is that the turbulence field is homoge-
neous. While this assumption is violated in nearly all applications of the Gaussian plume model, there is one
situation in modeling the dispersion of plume material over a hill where such a violation is severe.
A-4
-------
An elevated plume in a very stably-stratified flow may reside in a region of nearly laminar flow, with
very little turbulence. The turbulent boundary layer beneath this region can be very shallow. When Hc is
greater than the depth of the turbulent flow, the entire laminar region may flow over terrain, and an internal
boundary layer must form at the base of this layer. If this internal boundary layer were ignored, then the
diffusivity of plume material would remain virtually zero, and the effect of the terrain would only be seen in
the dividing streamline height partition of plume material.
To provide a more realistic treatment of the very stable limit, a simple turbulent surface layer is con-
sidered in LIFT. Within this surface layer, the vertical distribution of plume material is uniform, and the
concentration is obtained by sampling the vertical distribution in the absence of the surface layer, and taking
the average value over the depth of the layer.
The depth of the surface stable layer is estimated by focusing on the initial stages in the development
of such a layer, when turbulent mixing is very strong. Assume that the turbulence in the layer produces a
layer of constant potential temperature (see Figure A-l). As a result of the rapid mixing in the layer,
assume that the buoyancy effects are negligible in the layer, so that the wind speed takes on a logarithmic
structure in the vertical. Wind speed is continuous at the top of the layer, but the gradient is not.
Csanady (1974) discusses a theory for estimating the evolution of the depth of such a surface layer.
He uses the result of laboratory simulations by Kato and Phillips (1969) in which a layer is thickening into a
fluid with constant Brunt-Vaisala frequency N, as a result of a constant stress applied at the bottom of the
layer:
A-5
-------
8(2}
= Constant
y///////////////////////^^^
Potential Temperature Profile
'/W/W//WW////////////.
Wind Speed Profile
Figure A-l. Illustration of the structure assumed for the developing internal mixing layer over the hill above
Hc.
A-6
-------
J_ dh = 2.5u?9m -A.-12)
u, dt g A8 h ^
Here, h is the thickness of the layer, u* is the surface friction velocity, g is the acceleration due to gravity,
9m is the potential temperature of the mixed layer, and A0 is the jump in potential temperature at the top
of the layer. In essence, the entrainment depends on the turbulence generated at the surface of the layer,
and not on any assumptions about shear-generated turbulence at the interface.
The friction velocity in Equation A-12 is given by
= ka (A-13)
where u is the mean wind speed outside of the layer, k is von Karman's constant, and z0 is the roughness
length for the surface of the hill. Rewriting Equation A-12, an implicit equation for h :
(S7*)'
This implicit equation has two undesirable traits. The growth of the mixing layer is extremely rapid
for very small time-of-travel (downwind distance) from the point where the flow first encounters the hill
above Hc and the estimates of mixing depths "blow up" for weak stratification. The first of these traits is
removed by demanding that the rate-of-growth not exceed unity (slope of the interface equates 45 °). The
second trait is removed by limiting how large u/N can be in equation A-14.
For convenience, u/N is set equal to the minimum of u/N and Hc . This is done because we wish to
turn on the mixing layer primarily for those cases in which turbulence at plume height is decoupled from the
surface. As Hc becomes less than half the hill height, there is a strong likelihood that the turbulent bound-
ary layer in the approach flow encompasses most of the depth of the flow over the hill. Hence, in using the
minimum of u/N and Hc , the stability is artificially increased for Hc less than H/2 so that the depth of the
mixing layer will decrease. As Hc goes to zero, the mixing layer is completely absent, as desired.
A.4 FLOW MODEL (LIFT)
CTDMPLUS simplifies the treatment of stratified flow over a hill by separating the flow into two
regimes: a lower portion, below Hc, which is either blocked (in the case of a long ridge) or flows around the
A-7
-------
hill, and an upper portion, above Hc , which has sufficient kinetic energy to flow up and over the hill. This
simplification means that for a flow of constant speed, u , which lifts up and over the hill, one always has the
property that (H-Hc~)N/u < 1; or equivalently, that the Froude number for this portion of the flow
exceeds unity and the flow is not "strongly" stratified. This fact, coupled with the assumption that the hill is
not too steep (i.e., less than about 15°) enables one to use the linearized equations of motion for steady-
state Boussinsq flow (Smith, 1980):
PoU^L = _^L (A_15a)
dx dx
dv' dp'
pou = - (A-15b)
dx dy
dw' do'
P0u = --^- - p'gr (A-15c)
du.' dv'
dy dz
and
P' = -TI (A-15e)
These equations, in which d/dx, d/dy, and d/dz indicate derivatives with respect to downstream,
cross-stream, and vertical coordinates respectively, relate the perturbation velocities u ' , v ' , and w ' to the
perturbation density, p ', pressure p ' , and vertical fluid displacement, r| = n(x,y,z), and to the unper-
turbed initial velocity u and density p . Adding the kinematic condition for steady flow in a shear-free
approach flow
W = u
dx
Equations (A-15) can be reduced to the single partial differtial equation:
d2
dx'-
n2V^Ti = 0 (A-17a)
A-8
-------
where
n = , and N2 = -[ 1 . (A-17b)
u VP°/ a^
Yamartino (1987) discusses exact solutions to Equation (A-17), but steep terrain or strong stratifica-
tion tend to invalidate these because higher-order terms in Equations A-18 (e.g., quadratic in a perturba-
tions quantity) can be no longer neglected. Thus, the great complexity associated with the exact solution is
abandoned in favor of a more easily integrated, approximate solution. The easiest way to reach the
approximate solution is to integrate Equation (A-17) twice with respect to x ,
IT
J -n J -co
ay:
dx" dx' (A-18a)
and approximate the troublesome last term as
5 dx" dx' = n I -5 I , (A-18b)
IT
J -00 J - CO
where L x , Ly are the length scales of the hill in the along-, cross-wind directions, respectively. This
approximation, which proves to be exact for an infinite field of cosine hills, leads to the Helmholtz equation,
V2T] + m2T\ = 0 (A-19a)
where
m = n 1 + -f . (A-19b)
V L *, J
It has the simple particular solution
ac
(A-20)
A-9
-------
with G = cos(m£)/tf and R2 = x2 + y2 + z2 .
This is the Green's function solution for a delta function hill (i.e., a mathematically narrow but high
hill of unit volume). It is extended to a general hill shape h(x,y) via the convolution theorem and the lower
boundary condition that the flow at the surface follow the hill shape. Thus, for an arbitrary hill
7 =
dx' dy' h(x-x', y-y') G(x', y', z'}
(A-21)
and
a;
(A-22a)
where use of the height above terrain z', instead of absolute z reflects the fact that the lower boundary
condition has been linearized.
The other quantities needed for a complete description of the flow are obtained by using Equation
A-22a and going back to the basic equations of motion, yielding:
_
a
;(/)
«' * *,,. 2 r ^1 .
= (/) + n \ dx
u dx ay 7-oo ay
(A-22b)
(A-22c)
and
w
a
a2/
dx2
(A-22d)
It should be noted that lateral perturbation velocity, v', and the subsequent definition of lateral deflection,
6, as
u.
dy
2 C" ("' *l 1
nil dx" dx'
J-.J-. ay J
A-10
-------
make the additional assumption that streamline deflections are small enough that integrating along the
x-axis at y, z' from .*:--«> to x is equivalent to integrating along the streamline from x --<*>. Equa-
tion A-22e also provides another example of how higher-order terms are neglected, as a simple geometrical
picture of lateral deflection would indicate that the denominator factor u should actually be the true
jr-component of velocity, or approximately u. + a ' . Note that the deflections TI and 6 are those experi-
enced by a streamline that passes through the point ( x , y , z ' ) . They tell us where that streamline origi-
nates in the incident flow.
The local strain factors T\ and Th used in CTDMPLUS can also be related to / through these quan-
tities:
ay
n.2 (X (X ^-4 dx"d.x' ] (A-22g)
7-» 7-» dv I
As before, 7"h and Tj are factors that relate the spacing between adjacent streamlines in the
deformed flow to the spacing in the incident flow. The gradient of the deflection (either
dr\/d z or 36/3y ) measures the difference in the deflection experienced by adjacent streamlines and
hence, the degree to which the streamlines converge or diverge. In this backward-looking formulation, the
deflection is a function of the streamline position (coodinates) after deformation (see Figure A-2), where
the vertical position is the height above the surface of the terrain, z'. As detailed in the figure, the finite-
difference form of Equation A-22f (and A-22g) is readily obtained from the definitions of streamline posi-
tions before and after the deformation. Note that d\\/dz= 1 is a singular point. This condition would
result if streamlines that pass through z,' and z 2' were to originate at the same height in the incident
flow. Also note that Equation A-22f differs from Equation A-3 because the deflection in Equation A-3 is
that for a forward-looking formulation.
Before evaluating Equations A-21 and A-22 for a realistic hill shape, it is worthwhile to look back to
the solution, G - cos(m£)/£, given by Equation A-21, as properties of the delta function hill solution
will show up for finite hills as well. One positive aspect is that as stratification disappears (i.e., N, n, m go to
zero) the exact neutral result C = 1 / R is recovered. Thus, the approximation expressed by Equation A-18b
will not affect neutral flow but only the modifications to the neutral flow solution created by stratification.
That these stratification influences might not be too severe can be anticipated by noting that expansion of
A-ll
-------
Streamline 2
Streamline 1 -*'
0*2 A -r
" ^^ 4 » I A
.<*; i |n2 1A
h " &
z* - z1
_*2 *1
Az1
Az'-An
Az'
A-12
-------
the cosine term, as cas^rnR")*1 1 -( 1 /2) m2 R2, indicates that changes to the flow are second-order in m
(i.e., m2) rather than first-order. However, a disturbing aspect of this cos(mR) dependence is its iso-
tropic nature; that is, z has no special significance in the equation, despite the fact that the density stratifica-
tion and thus the atmosphere's "springiness" is az-oriented phenomenon. Exact solutions to Equation A-17
do not display this isotropic behaviour. Thus, one physically significant ramification of approximation
(Equation A-18b) becomes apparent. Such isotropic and thus lateral springiness does, however, occur in the
aforementioned infinite field of cosine hills problem but steps must be taken to suppress this effect in the
single isolated hill problem.
Finally, the connection to the neutral limit solution, G = 1 / R, suggests a way to inject wind speed
shear (with height) back into this shear-free solution. Crapper (1959) points out that the vertical deflection
in a neutral, shear-flow is just
n - I ^ I -, (A-233)
whereas G = 1 {R yields
(A-23b)
in the shear-free case. Unfortunately, even this simple factor is awkward to superimpose back onto G
exactly. The approximation chosel for this model yields a final Green's function of
Any unwanted residual consequences of this approximation are suppressed by enforcing the constraint that
n = /i(x,y) at z'-O.
Referring back to Equation A-22, one notes that all quantities of interest depend on integrals and
derivatives of the basic quantity 7 , given by Equation A-21. Evaluation of 7 requires a carefully chosen
shape of hill. Selection of an appropriate hill shape was governed by
i) the desire to treat as general a shape and orientation as possible, and
A-13
-------
ii) the necessity for performing the integrations in Equation A-21 analytically.
A Gaussian-shaped hill is used because it contains the overall orientation and scales of the "cut-off hill, and
because it has a tractable mathematical form. The most general Gaussian shape that is considered involves
a hill of height h , having elliptical contours with major axis.La , and a minor axisLb , oriented at an angle
ip which is counter-clockwise with respect to the flow direction. For a coordinate system with the x-axis
aligned with the direction of flow, the hill is prescribed by
h(x,y) = h exp< -
(A-25)
where
Ll
cos2ip sin2ip
COS2Tp
and
[---
US Ll
cosip sinip
For the symetric hill with La - L - L, the expression simplifies greatly, as Lx- Iy = L and y = 0 -
In fact, once the convoluted form of the hill function (i.e. h(x,y) is rewritten as h(x-x',y-y'))is
substituted into the integral for /, one finds that the exact analytic integration can be accomplished for a
receptor position ( x, y, z') at the crest (and just off the crest) of such a symetric Gaussian hill. But this
analytic result must be approximated as adequate software does not widely exist for the complex error func-
tion portion of this solution. The expression finally used is
/ =
"(0)
1/2
cos(mz')
-a,
sin(mz')
m-a.-
-I
Lx )V
. *> / \
m cos(mz') +
(A-26)
A-14
-------
where
B0mL
= RL
-| 1/2
a. is from the shear equation u(z') = u(0)
a, =
I a
2 u(0)
a, =
2
,3/2
1/2
and
The two parameters, 50 and R^, arose because there was a conflict between small and large argument,
low-order expansions of the complex error function, and ambiguity in the definition of Lz , respectively. The
parameter B0 was set to (n/2)"2, the geometric mean of the small and large argument limit values,
whereas RL was tuned to ^(2)^0.693 based on optimization studies using the tow tank data described
by Snyder et al. (1986).
The expression for 7 given above is best described as an admixture of analytic results for the symmet-
ric hill and empirical extentions for the asymmetric hill subject to the constraints that the surface boudary
condition be obeyed and that known results be covered for the 2-d ridge limit of !>.->«> and V = 0.
Higher-order terms (e.g., (y/iy)2 or other such terms arising far from the crest) are neglected. In addi-
tion, we took the liberty of switching-off the ( xm / L ) term upwind of the hill (i.e., x. < 0 ), as it seemed
to create excessive early streamline rise as the hill was approached and because appropriate damping terms
in xjj, has already been neglected in the denominator.
Evaluation of the flow variables in Equation A-22 involve straightforward integrations and differenti-
ations of the basic quantity, /, given by Equation A-26. To follow this process most easily and to bridge the
gap between mathematical notation (e.g., Equation A-26) and the computer code in SUBROUTINE
FLOW, we note that Equation A-26 may be rewritten as
A-15
-------
/ = HHXYM TZF )_^T, -AFC (A-27)
i-O
where HHXYM = h(x , y) Z./(l * 6;) contains the dominant x andy dependences in/, and
-2 r .. /-r\\ 11/2
TZF = I 1+
I- n
contains the overall multiplicative z' dependencies arising from the complex error function approximation
and wind speed shear, respectively, and is denoted TZFAC in the code. The sum over terms, T, AFt, is
organized such that all the trigonometric functions involving the vertical coordinates 2', or the 'angular'
variable, mz', are contained \nAF\, whereas multiplicative constants and remaining coordinate depen-
dences (i.e., xm/Lx} are absorbed into T,. Thus, for the quantity/, we have
T2 - -az - | _= | for xm>0 and
= 0 for xm<0 ,
whereas the AFj quantities, which along with their 2' derivatives
AD, -J|
ADD,
will be used repeatedly in integrated and differentiateed forms of I, are presented in Table A-l.
The grouping of terms by variable expressed by Equation A-27 facilitates computation of derivatives
by the well-known "chain-rule"
A-16
-------
TABLE A-l
AF{ TERMS AND THEIR DERIVATIVES
AFQ =1 + cos(mz')
ADo = ~~~V-~'> - ! i + -j-]m sin(m;0
*« .
sinC/nx') . , - . ,
^ '11 -_ m^ cos(mz'
sin(mz')
! = cos(mz')
2 1 a
AF2 = m cos(mz')+ + - sin(mz')
L n
-------
Thus, for example, the quantity dl/dz', denoted/DZ in the code, becomes
; = TZFDZ I + HHXYM TZF ^ T, ADi (A-28)
dz i_0
where TZFDZ = (d/dz" TZF^/TZF
or
2//:n 1 a
TZFDZ =
2 u(z')
Yet higher derivatives show the economy of the chosen / expression compartmentalization and its resistance
to the errors that can occur when many, many terms are generated through differation. For example,
d 2 / / d z'2, denoted TDZZ in the code, is simply
= TZFDZZ I + 2 TZFDZ
dz'2 dz'
2
+ HHXYM TZF £ T,. ADD, (A-29)
i-O
where
TZFDZZ =
or
TZFDZZ = (TZFDZy + t 2
2 + 2./L2n + l a2
+ 2 a2(z
A-18
-------
Derivatives of/ with respect tox andy are also required, and hence necessitate derivatives of the
hill shape function h fay) terra HHXYM. Determining the quantities
A (HHXYM) = -2(^0 - HHX .
ax U.J LX
ay
= -2
( Xm
1 + 2 -^
HHXYM
- -2y- HHXYM , and
L
2 (HHXYM) = -2
*r .
HHXYM
where
^x, then enables one to rapidly develop the terms:
= -2| | + HHXYM TZF
ax V /.* y /.*
(A-30)
where now TO = 7" ', = 0 , because
T" = dTi/dx in this case and TO and 7j are independent of JT , and
as all (
-------
ti
ty
HHXYM TZF
i-O
(A-33)
where 7" * are redefined now as
< ^ ,^
T , = - and yield
dy
a2
and T2 = - ;
Ly
-2
2!
(A-34)
3 + 2
dl/dy
; and
(A-35)
a2/
dxdy
= -2
Y + 2
'2/2
dl/dy
(A-36)
The above seven quantities are denoted /DA", TDAX", /ZMTZ, IDY, IDYY, ID3Y, and /DAT, respectively, in
the computer algorithm FLOW.
Computation of the integrated quantities, defined now symbolically as
I* - fdx' /(*') and
J _o,
A-20
-------
and denoted IX and IXX respectively in the code, now requires that several indefinite integrals can be com-
puted. Referring back to Equation A-27 shows that -c"s arise in HHXYM due to the Gaussian hill shape
function h(x,y) , expressed by Equation A-25, as well in the coefficient T-± . Thus, it is convenient to consider
first the integrals
exp-
\Ly
= (A-38)
for/ = 0,1 and where xm = x + ^L\y ,as before. Upon completion of the square in the exponential, the
results can be found from integral tables to be
>2 ~1
g* = Lx exp
y_
Uyy
(A-39)
where YP" 1 - Y^x^y and where the G" are presented in Table A-2. Subsequent integration in x to yield
Sf", defined as"
A-21
-------
TABLE A-2
DIMENSIONLESS GJ AND GJ* TERMS
and
definitions yield:
I
7-oo
de exp(-e2)e; =
XX
i
de
x,
L.
Gx _
n - -'
x _ _
XX _
-exp
, and
- XX _ _ _ I
These are included in flow algorithm as COX, Glx, GOXX, and G1XX, respectively.
Note: Subsequent numerical studies showed that it was preferable to substitute [l-
for [1 + ERF(xm/Lx)]. This has the major effect of allowing stream lines to return laterally to
their initial upwind positions more rapidly after passing the hill. This adjustment can be rationalized
on the basis that Equation A-26 is most valid near the crest of a hill.
A-22
-------
«" - /'
J ~\L.XJ
yields formally
r f v ^2 "1 xx
', , (A-40)
gxx = Ll exp ( ^
where the G ** are also presented in Table A-2. Computation of the G ** are made possible by use of the
infrequently tabulated inetgral of an error function:
( dx ERF(x} = x ERF(x} f 4= exp(-.x2) .
^ J n.
Armed with these quantities, one may express / * as
/- - (»\ . TZF
V *-x )
where
HHLX2F
hL.L2x
is the variable name used in the code, and where
T* - rx
1 o - ^o <
TI = -ttj GO , and
T*2 = -a2 G* .
Similarly, the doubly integrated quantity, / **, is then just
7XX = HHLX2E TZF ^ Ti ' AFi (A-42)
i-O
A-23
-------
where the T * are redefined as
7* = -OLI
, and
-a.
Achieving the goal of providing all the quantities specified in Equation A-22 now only demands that
a few differentiation operations be performed with respect toy . These are in turn facilitated by the obser-
vations that since the G * and G** are only functions of xm/Ix, but xm = yLly contains^ , then the
differentiation by y is equivalent to a differentiation by x except that d/dy generates the multiplicitive
factor yl» instead of 1/IX. Thus we obtain the useful quantities
H, = £(C*) (A-43)
= yLx exp
x,
referenced as HHJ (i.e., HHO and HH1 for j = 0,1 ) in the code, and the relations
G",' Yi,G* , and
for use in the differentiation process. Keeping track of the terms generated via the chain rule, one then
obtains:
i V
= -2
HHLX2E TZF 7 AF i (A-44)
i-o
where
A-24
-------
T
7
, and
dlx
dy
y_
L\
2(75 YP/'
HHLX2E
TZF
i-O
AFt (A-45)
and
-.2. , xx
2 T- Y
-4
+ HHLX2F TZF
i-O
(A-46)
'i AFt -
where
in both Equation A-45 and A-46; and
T\
n
H 0 ,
-a, H0 , and
= -a.
A-25
-------
4,^
\t})
(VP]
6 0
\.i-l)
6(M
6 U
(&" '
dixx r
<3y |_
(y\ d*I
(Lj dv
3 + 2yp
2v fyVl
pv^-/; J
XX
2
y
Ly
HHLX2E TZF
(A-47)
AF
where
T =
-2.H l
+2at
, and
T* = -a
The above four quantities in Equations A-44 through A-47 are denoted by IXXDY, IXDY, IXXDYY, and
IXXD3Y, respectively, in the FLOW algorithm.
One may notice that the quantities <33/dy3 in Equation A-35 and d3/"/dy3 in Equation A-47
have been defined for no apparent purpose. However, they appear when one considers derivatives £6/dy
and <326/<3y2 of 6 from Equation A-22e, that are needed for a "fully-implicit" style computation of the
lateral deflection. This compensates for the fact that lateral deflections can be quite large and therefore
badly violate the "small deflection" assumption invoked for Equation A-22e. The resulting lateral deflections
are unfortunately smaller than before the correction, but streamline crossover situations (which create a
severe problem for an iterative streamline solver) for large aspect ratio hills are eliminated.
Examples of the flow deformation resulting from the algorithm in CTDMPLUS are shown in Fig-
ures A-3 (facing downwind) and A-4 (cross-flow view).
A-26
-------
160 ^
140 -j
130 -]
120 i
0
^"s^
J
$ -
-
^x^
^*v
- . ,,
f 2
2
3
*
«
^W
^V
%
.
«3 '*
* 5
*
g
1
m
*^s_
"a.
^N
Ts
s ;
M
»
^^^,
-. « '
i
-*-»~-~.
600
Figure A-3. Behavior of an array of marked streamlines emitted in vertical columns upwind of a Gaussian hill,
as seen by an observer looking downwind in stable conditions (Froude # = I). Initial streamline
spacing is 20 m in the vertical (up to 100 m) and 69.3 m in the horizontal upwind of the hill (see
streamline column #'s 0-7 at base of diagram and above hill).
ISO
170
160
160
140
130
120
110
100
90
60
70
60
60
40
30
20
10
.X
-600
300
-100
100
300
600
X(m)
Figure A-4. Pattern of streamlines that pass over the crest of the hill as viewed from the side in stable condi-
tions. Flow is from left to right. Initial streamline spacing is 20 m in the vertical up to 100 m.
A-27
-------
A^ MODEL FOR STREAMLINES (WRAP)
The central features of WRAP are the distance between the streamline that passes through the
source and the stagnation streamline, and the relative locations of the source, the receptor, and the point of
impingment. These require a flow model to obtain streamlines. Because the flow is two-dimensional in this
strongly stratified flow, potential flow solutions are used to obtain the streamlines.
The hill below Hc is represented as a cylinder of elliptical cross-section, set on end. This shape is
chosen because it contains the overall scale of the hill, and it is simple enough that streamlin patterns can be
expressed analytically. As already stated, the ellipse used for the entire flow below Hc is the result of fitting
an ellipse to the height-contour that corresponds to the minimum of Hc and the plume height. The poten-
tial flow solution is expressed in elliptical coordinates following the notation of Batchelor (1970).
Let the x-axis be aligned with the major axis of ellipse, let a be the length of the semi-major axis, and
b be the length of the semi-minor axis. Then the elliptical coordinates (|i, v) are related to the cartesian
coordinates by the relations
x2 = (a2-b2) cosh2(n + n0) cos2(v)
= (a^-b^) cosh2(n + n0) cos2(v) (A-48)
Note that |i = n' - H0, where H0 is the value of |i' along the boundary of the ellipse ( M. is the constant
along the ellipses in the family of confocal ellispses of which the boundary of the hill is a member). Using
the elliptical coordinates, a streamline in the flow is given by
ip = 5M(a + b) sinh(pi) sinO + a^) (A-49)
where 5. is the speed of the incident flow, and a., is its direction. Note that a,,, is zero when the flow is
directed along the -x direction, and increases in the clockwise direction, as noted in Figure A-5.
The wind speed and direction measured at a tower near a hill may be influenced by the presence of
the hill, so that the observed mean speed and direction must be estimated at infinity to determine the inci-
dent flow. This is calculated from the theory of two-dimensional potential flow around an ellipse. Let
( H T v T ) be the coordinates of the point at which the speed and direction ( 5 T, $ r ) are measured, where
4»T is the direction counter-clockwise from the major axis of the ellipse. Then the wind direction is parallel
to the tangent to the streamline defined by i|> = ip T at the point ( n T , v T ) so that a , is given by
A-28
-------
Impingement Point
0
Figure A-5. Definition of modeling variables, illustrating in particular the coordinate system in which the
xB-axis is aligned with the tangent to the stagnation streamline at the impingement point (the
(3 -coordinate system). The coordinates along the xp-axis of the source are denoted by xlB,
xoB, and xrp, respectively.
A-29
-------
[, , I
^ + tan(T)tanh(|j.T)cos2(vT) U r tan((()r) [tan2(vT) + tanh2((j,T)]
cosh (p. r ) j
(tan2(vT)tanh2(nr)
v ^ /y vi-,y
f tanh((iT)
2, , +
1 cos'Cvj.)
tan(d)r)
tan(vr)
cosh
2(^r)j
(A-50)
and the speed at infinity (far away from the influence of the hill) is given by
/
V
vr 2,1
, + / ^ 1
sinhj;
sin2(vr* aw)
i + ~
sinh^CHj.)
1 ^ \
tanh(n7)
J
1/2
(A-51)
The quantity r is the ratio of the major axis length to the minor axis length, a/b .
Once the incident flow angle relative to the major axis of the ellipse and the speed S. are known,
the stream function through the source ( V,) can be calculated from
vs + cO . (A-52)
The stagnation point can also be calculated. Along the stagnation streamline, Vo = °, so that Equation
A-49 becomes
0 = -
sinh(|j.)
(A-53)
Because this must be satisfied for all n, v must be equal to -au along the stagnation streamline at the
stagnation point. This coordinate system is needed to provide a convenient Cartesian coordinate system that
allows the streamline through the source to be a single valued function of x for all aw. The rotation angle,
(3, is given by
tan(f3) = -
(A-54)
A-30
-------
The distance between the streamline through the source (ips) and the stagnation streamline
(\p0 = 0) far from the hill is related to the value of y* and the wind speed at infinity, S. . Because the
speed of the flow equals the gradient of the stream function far from the hill, we find
d = (A-55)
However, because au may be measured closer to the hill, the speed at the source is substituted for S. to
estimate d near the source.
A-31
-------
-------
APPENDIX B
TEST CASE FOR CTDMPLUS
Input and output files are shown for a test run of CTDMPLUS. Specific
figure references to these files are given below.
File Name
"CTDM.IN"
"TERRAIN"
"RECEPTOR"
"PROFILE"
"SURFACE"
"RAWIN"
"CONC"
"CTDM.OUT"
Context
Input options; source and
receptor information
Hill Parameterizations
Receptor data
Meteorological tower data
Preprocessed meteorological
variables
Preprocessed rawinsonde data
Text file of concentrations
Output listing for this run
Reference
Figure B-l
Figure B-2
Figure B-3
Figure B-4
Figure B-5
Figure B-6
Figure B-7
Figure B-8
B-l
-------
CTDMPLUS TEST RUN
3121101111
1.0 0.3048 39.5915 89.4885 6 1
TOWER 0.0 0.0 945.0
STACK-1 0.0 0.0 945.0 189.70 3.36 410.15 25.06 455.05 0
STACK-2 0.0 30.0 940.0 190.00 3.36 420.15 25.06 255.05 0
ENDS
0.76
Figure B-l. File "CTDM.IN" used in test run of CTDMPLUS (see Table 3-3).
B-2
-------
1 13
900.000
1000.000
1100.000
1200.000
1300.000
1400.000
1500.000
1600.000
1700.000
1800.000
1900.000
2000.000
2100.000
900.000
1000.000
1100.000
1200.000
1300.000
1400.000
1500.000
1600.000
1700.000
1800.000
1900.000
2000.000
2100.000
.6176E+03-.
.6176E+03-.
.6176E+03-.
.6176E+03-.
.6176E+03-.
.6176E+03-.
.6176E+03-.
.6176E+03-.
.6176E+03-.
.6176E+03-.
.6176E+03-.
.6160E+03-.
.6072E+03-.
.6176E+03-.
.6176E+03-.
.6175E+03-.
.6175E+03-.
."6175E+03-.
.6175E+03-.
.6175E+03-.
.6175E+03-.
.6175E+03-.
.6175E+03-.
.6174E+03-.
.6181E+03-.
.6290E+03-.
2240E+04
2092E+04
2092E+04
2092E+04
2092E+04
2092E+04
2092E+04
2092E+04
2092E+04
2092E+04
2092E+04
2092E+04
2095E+04
2110E+04
2093E+04
2093E+04
2093E+04
2094E+04
2094E+04
2094E+04
2094E+04
2095E+04
2095E+04
2096E+04
2097E+04
2098E+04
2086E+04
180
180
180
180
180
180
180
180
180
170
170
160
160
175
175
174
174
173
172
171
170
168
167
166
170
180
PIEDMONT HILL
.000 3072.000 1294
.000 2953.000 1195
.000 2834.000 1097
.000 2715.000 998
.000 2595.000 953
.000 2454.000 927
.000 2321.000 847
.000 2210.000 796
.000 2055.000 745
.000 1887.000 679
.000 1669.000 625
.000 1268.000 582
.000 1032
.111
.111
.666
.128
.470
.650
.589
.171
.189
.729
.960
.275
.000
3
3
3
3
3
3
3
3
3
4
4
5
2
.000
.557
.960
.949
.936
.924
.898
.863
.871
.924
.132
.898
.122
.000
405
3
3
3
3
3
3
3
3
2
2
2
2
2
.400
.700
.000
.300
.600
.100
.600
.200
.600
.500
.200
.200
.300
.076
.437
.443
.384
.311
.267
.158
.019
.853
.606
.303
.220
.000
-
2061.
1930.
1862.
1791.
1716.
1638.
1555.
1461.
1356.
1231.
1078.
966.
1115.
602
771
632
539
603
731
791
472
096
882
237
383
968
764.782
708.782
680.194
654.256
627.123
596.291
566.836
535.555
501.217
465.700
423.987
348.302
266.738
Figure B-2. File "TERRAIN" used as input in test run of CTDMPLUS (see Table
3-5).
B-3
-------
MET TOWER 710. -400. 0.0 198.2. 1
FOLLY 550. -500. 0.0 1850. 1
TREEHOUSE 550. -700. 0.0 1978. 1
SLIPPERY HILL 490. -1030. 0.0 2096. 1
RIDGE END 770. -160. 0.0 1939. 1
OVERLOOK 230. -980. 0.0 2090. 1
SOUTHSIDE 0. -1480. 0.0 2110. 1
COW PASTURE 320. -1480. 0.0 2208. 1
Figure B-3. File "RECEPTOR" used in test run of CTDMPLUS (see Table 3-6).
B-4
-------
80
80
80
80
6 26
6 26
6 26
6 26
1
1
10
10
10.
100.
10.
100.
0
1
0
1
300.
300.
328.
328.
1.2
3.9
2.2
2.5
299
299
299
299
.3
.3
.4
.4
5.0
5.0
-999
-999
.03
.03
.49
.49
-999.9
-999.9
-999.9
-999.9
Figure B-4. File 'PROFILE" used in test run of CTDMPLUS (see Table 3-7).
B-5
-------
80 6 26 178 1 92. 30. 0.057 11.2 0.150E+00
80 6 26 178 10 -999. 1242. 0.293 -7.4 0.150E+00
Figure B-5. File "SURFACE' used in test run of CTDMPLUS (see Table 3-8).
B-6
-------
a
6201 14842
990.7/ 200./302.
850.0/1538./290.
788.0/2183./289.
6201 14842
991.3/ 200./292.
900.0/1041./293,
800.0/2059./294.
80626 0
5/170/ 2
0/199/ 2
6/ 28/ 1
80 62612
O/ O/ 0
1/268/ 2
2/32S/ 3
43
978.0/ 315.
839.0/1649.
750.0/2604.
46
976.0/ 335.
880.0/1236.
750.0/2612.
/300.8/173/ 2
/289.4/194/ 2
/290.6/ ll/ 3
/297.1/235/ 1
/291.6/266/ 2
/289.S/ I/ 5
12
950.0/ 573./298.6/179/ 2
821.0/1834./288.4/187/ 1
749.0/2615./290.7/ 12/ 3
11
959.0/ 489./297.1/240/ 3
870.0/1334./293.1/262/ 2
700.0/3194./284.9/ 6/ 6
900.0/1047./294.5/195/ 2
800.0/2054./289.1/140/ 1
700.0/3187./2S5.7/ 32/ 3
950.0/ 572./2%.5/249/ 3
850.0/1535./294.2/262/ 2
Figure B-6. File "RAWIN" used in test run of CTDMPLUS (see Table 3-9).
-------
80 6 26 1 1 8 microS/M**3
0.541E+00 0.592E-03 0.619E-08 O.OOOE+00 0.214E-02 O.OOOE+00 O.OOOE+00 O.OOOE+00
80 6 26 10 6 8 microS/M**3
0.213E-01 0.152E+00 0.538E+00 0.960E+00 0.188E-05 0.155E+01 0.115E+01 0.991E+00
Figure B-7. file "CONC", concentration file from the CTDMPLUS test run (see Table 3-11).
B-8
-------
CTDMPLUS VERSION 1.0
CTDMPLUS TEST RUN
PAGE
INPUT/OUTPUT SWITCHES (1
USE THIS OPTION, 0 = DO NOT USE OPTION )
STABLE ONLY,
ICASE: INCLUDE CASE-STUDY PRINTOUT: 0 = NONE, 1
2 = UNSTABLE ONLY 3 = ALL HOURS
ITOPN: CREATE TOP A TABLE AT END OF RUN
ICONC: CONCENTRATION OUTPUT: 0 = NONE, 1 = BINARY, 2/3
TEXT
INTERNAL PROGRAM SWITCHES
IMIX: (IF 1, USE ON-SITE MIXING HEIGHT OBSERVATIONS
(OFF-SITE IF NOT AVAILABLE); IF 0, VICA VERSA) 1
IWS1: (IF 1, SET MINIMUM WIND SPEED TO 1.0 M/S) 1
ISIGV: (HORIZ. TURB. INTENSITY DATA (0=SIGMA-THETA 1=SIGMA-V) 0
IWD: (IF 1, SCALE WIND DIRECTION WITH HEIGHT) 1
ICHIQ: (IF 1, MODEL OUTPUTS CHI/Q; OTHERWISE IT OUTPUTS CHI) 1
ISOR: (IF 1, MODEL GIVES SOURCE CONTRIBUTION TABLE) 1
IUNSTA: (IF 0, MODEL WILL NOT READ RAWIN FILE OR CALCULATE
UNSTABLE HOURS) 1
"FOR HORIZONTAL SCALE, MULTIPLY USER UNITS BY 1.0000 TO GET METERS.
FOR ELEVATION, MULTIPLY USER UNITS BY 0.30A8 TO GET METERS.
SITE LATITUDE (> 0 IF NORTH) = 39.591
SITE LONGITUDE (> 0 IF WEST) = 89.A89
SITE TIME ZONE (> 0 IF WEST) = 6.
POLLUTANT # (FOR HOURLY EMISSIONS) = 1
METEOROLOGICAL TOWER COORDINATE INFORMATION:
X-COORD: 0.000 (USER UNITS) * 1.0000
Y-COORD: 0.000 (USER UNITS) * 1.0000
ELEVATION: 9A5.000 (USER UNITS) * 0.30A8
0.0 (METERS)
0.0 (METERS)
288.0 (METERS)
Figure B-8. File "CTDM.OUT", output listing from CTDMPLUS test run.
B-9
-------
CTDMPLUS VERSION 1.0 PAGE
CTDMPLUS TEST RUN
***SOURCE INFORMATION***
STK NAME
1
1 STACK-1
2 STACK- 2
EMISSION
RATE
(G/S)
455.
255.
05
05
LOCATION
X Y
(M)
0.
0.
00
00
(M)
0.00
30.00
STK
HT
(M)
191
190
.2*
.0
STK
DIA
(M)
3
3
.36
.36
GAS
TEMP
(K)
410.1
420.1
EXIT
VEL
(M/S)
25.06
25.06
COMMON BASE ELEVATION = 286.5 (METERS).
THIS BASE ELEVATION IS USED FOR ALL STACKS IN THIS RUN;
ALL STACK HEIGHTS MARKED WITH * HAVE BEEN ADJUSTED TO
RETAIN THE ACTUAL ELEVATION OF THE TOP OF THE STACKS.
MULTIPLY HORIZONTAL USER UNITS BY: l.OOOE+00 TO CONVERT TO METERS
MULTIPLY VERTICAL USER UNITS BY: 3.048E-01 TO CONVERT TO METERS
Figure B-8. (Page 2 of 14).
B-10
-------
CTDHPLUS VERSION 1.0
CTDMPLUS TEST RUN
PAGE
REC
NO.
IDENTIFICATION
RECEPTOR INFORMATION
EAST NORTH
COORD COORD
(USER UNITS)
HEIGHT ABOVE GRD LVL
LOCAL GRD LVL ELEVATION HILL
(USER UNITS) (USER UNITS) NUMBER
MET TOWER
FOLLY
TREEHOUSE
4 SLIPPERY HILL
5 RIDGE END
6 OVERLOOK
7 SOUTHSIDE
8 COW PASTURE
710.00 -400.00
550.00 -500.00
550.00 -700.00
490.00 -1030.00
770.00 -160.00
230.00 -980.00
0.00 -1480.00
320.00 -1480.00
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1982.0
1850.0
1978.
2096.
1939.
2090.
2110.0
2208.0
.0
.0
.0
.0
MULTIPLY HORIZONTAL USER UNITS BY: l.OOOE+00 TO CONVERT TO METERS
MULTIPLY VERTICAL USER UNITS BY: 3.048E-01 TO CONVERT TO METERS
Figure B-8. (Page 3 of 14).
B-ll
-------
CTDMPLUS
VERSION 1.0
PAGE
TERRAIN INFORMATION (USER UNITS FOR ALL DATA)
MULTIPLY HORIZONTAL USER UNITS BY: l.OOOE+00 TO CONVERT TO METERS
MULTIPLY VERTICAL USER UNITS BY: 3.0i3E-01 TO CONVERT TO METERS
TERRAIN INFORMATION (CONT.)
HILL t 1 PIEDMONT HILL
HILL TOP: 2240.0 (USER UNITS)
BEST FIT ELLIPSE INFORMATION FOR WRAP: PIEDMONT HILL
CONTOUR X-COORD Y-COORD MAJOR AXIS ELLIPSE AXIS LENGTHS
HEIGHT (HILL CENTER) AZIM. FROM N MAJOR MINOR
900.0
1000.0
1100.0
1200.0
1300.0
1400.0
1500.0
1600.0
1700.0
1800.0
1900.0
2000.0
2100.0
617.600 -2092.000
617.600 -2092.000
617.600 -2092; 000
617.600 -2092.000
617.600 -2092.000
617.600 -2092.000
617.600 -2092.000
617.600 -2092.000
617.600 -2092.000
617.600 -2092.000
617.600 -2092.000
616.000 -2095.000
607.200 -2110.000
180.0
180.0
180.0
180.0
180.0
180.0
180.0
180.0
180.0
170.0
170.0
160.0
160.0
3072.000
2953.000
2834.000
2715.000
2595.000
2 ASA. 000
2321.000
2210.000
2055.000
1887.000
1669.000
1268.000
1032.000
129A.AOO
1195.700
1097.000
998.300
953.600
927.100
8A7.600
796.200
7A5.600
679.500
625.200
582.200
A05.300
MULTIPLY HORIZONTAL USER UNITS BY: l.OOOE+00 TO CONVERT TO METERS
MULTIPLY VERTICAL USER UNITS BY: 3.0A8E-01 TO CONVERT TO METERS
TERRAIN INFORMATION (CONT.)
He CUT-OFF HILL INFORMATION FOR LIFT: PIEDMONT HILL
CONTOUR X-COORD Y-COORD MAJOR AXIS < INVERSE POLYNOMIAL VARIABLES >
HEIGHT (HILL CENTER) AZIM. FROM N MAJ EXP MIN EXP MAJ SCALE MIN SCALE
900.0
1000.0
1100.0
1200.0
1300.0
1AOO.O
1500.0
1600.0
1700.0
1800.0
1900.0
2000.0
2100.0
MULTIPLY
MULTIPLY
SURFACE
HILL S
ZO (M)
617.600-2093.000
617.600-2093.000
617.500-2093.000
617.500-209A.OOO
617.500-209A.OOO
617.500-209A.OOO
617.500-209A.OOO
617.500-2095.000
617.500-2095.000
617.500-2096.000
617.AOO-2097.000
618.100-2098.000
629.000-2086.000
175
175
174
17A
173
172
171
170
168
167
167
170
180
.1
.1
.7
.1
.5
.6
.6
.2
.2
.7
3
3
3
3
3
3
3
3
3
A
.0 A
.3
.0
HORIZONTAL USER UNITS BY:
VERTICAL USER UNITS
ROUGHNESS LENGTH OF
1
0.760
BY:
EACH
3.
5
2
.557
.960
.9A9
.936
.92A
.898
.863
.871
.92A
.132
.898
.122
.000
l.OOOE+00 TO
OA8E-01
3
3
3
3
3
3
3
3
2
2
2
2
2
.076
.A37
.AA3
.38A
.311
.267
.158
.019
.853
.606
.303
.220
.000
CONVERT
TO CONVERT TO
2061.
1930.
1862.
1791.
1716.
1638.
1555.
1461.
1356.
1231.
1078.
966.
1115.
602
771
632
539
603
731
791
472
096
882
237
383
968
76A
708
680
654
627
596
566
535
501
A65
A23
348
266
.782
.782
.194
.256
.123
.291
.836
.555
.217
.700
.987
.302
.738
TO METERS
METERS
HILL:
Figure B-8. (Page 4 of 14).
B-12
-------
CTDMPLUS VERSION 1.0 PAGE 5
MAP EDGES: XMIN = -521., XMAX = 1291., YMIN = -1480., YMAX = 30.
* = SOURCE, RECEPTORS SHOWN BY HILL * (0-9,A-Z)
1
1
Figure B-8. (Page 5 of 14).
B-13
-------
CTDMPLUS VERSION 1.0
CTDMPLUS TEST RUN
PAGE
INPUT MET DATA
YR MO DA HR
80 6 26 1
ADJUSTED
HEIGHT
(M)
11.5
101.5
FROM SURFACE AND PROFILE (NOTE: ****** = MISSING DATA):
MIXING
HEIGHT
(M)
92.0
WIND
DIR.
(DEG)
300.0
300.0
MONIN- SFC
SFC SFC OBUKHOV ROUGH.
TEMP U* LENGTH LENGTH
(K) (M/S) (M) (M)
299.3 0.057 11.2 0.1500
<-WIND SPEED-> AMB. SIGMA-
SCALAR VECTOR TEMP THETA
(M/S) (M/S) (K) (DEG)
1.20 1.20 299.30 5.0
3.90 3.89 299.30 5.0
SIGMA-V SIGMA-V
(M/S) (M/S)
0.20 0.03
0.34 0.03
NOTE: SCALAR WIND SPEEDS USED IN CTDM ARE SET TO A MINIMUM OF 1 M/S
NOTE: HEIGHTS ARE REFERENCED TO THE COMMON STACK BASE ELEVATION
THE ADJUSTMENT TO THE INPUT HEIGHT IS 1.5 METERS.
< SOURCE INFORMATION >
SOURCE 'OS TS VS BUOY FLUX MOM FLUX
ft (G/S) (K) (M/S) (MA/S3) (MA/S2)
1.0 A10.1 25.06
187.5
1293.A
FINAL PLUME
RISE
(M)
138.05
VARIABLES AT
PLUME HEIGHT:
HEIGHT
(M)
WDIR USCAL UVECT SIGV SIGV DTHDZ
(DEG) (M/S) (M/S) (M/S) (M/S) (DEG/M)
329.3 300. 3.90 3.89 0.34 0.0390 0.0098
Figure B-8. (Page 6 of 14).
B-14
-------
CTDMPLUS VERSION 1.0
CTDMPLUS TEST RUN
PAGE
INFORMATION FOR HILL 1: PIEDMONT HILL
HCRIT = 178.5 M; FROUDE # ABOVE HCRIT = 1.00
WRAP INFORMATION:
DISTANCE FROM SOURCE TO HILL CENTER = 2181.3 M;
ELLIPSE AXIS LENGTHS: MAJOR = 2159.4 M; MINOR =
MAJOR AXIS AZIMUTH FROM NORTH = 180.0 DEC
DISTANCE TO PRIMARY IMPINGEMENT POINT = 287.3 M
WRAP HT = 178.5 M
834.4 M
LIFT MIDPOINT HT =
554.5 M
LIFT INFORMATION:
DISTANCE FROM SOURCE TO HILL CENTER = 2183.1 M
ELLIPSE AXIS LENGTHS: MAJOR = 1528.0 MjHINOR
MAJOR AXIS AZIMUTH FROM NORTH = 171.6 DEC
DISTANCE ALONG FLOW FROM SOURCE TO HILL CENTER = 1581.8 M
CROSSFLOW DISTANCE FROM SOURCE TO HILL CENTER = 1504.7 M
DISTANCE TO PRIMARY IMPINGEMENT POINT = 287.1 M
287.4
SRC-RECP SRC-RECP RECEPTOR EFF. FLAT HIL1
L DISTANCE DISTANCE HT ABOVE SRC-RECP TERRAIN EF:
REC / ALONG FLOW CROSS FLOW STK BASE HT DIFF SIG-Y SIG-Z SIG
I V (M)
1 L 815.
1 W 814.
1 TOTAL
2 L 726.
2 W 732.
2 TOTAL
3 L 826.
3 W 838.
3 TOTAL
4 L 939.
4 W 963.
4 TOTAL
5 L 747.
5 W 737.
5 TOTAL
6 L 689.
6 W 716.
6 TOTAL
(M)
36
86
189
86
347
86
615
86
-203
86
679
86
.6
.0
.0
.0
.5
.0
.0
.0
.3
.0
.4
.0
317
317
277
277
316
316
352
352
304
304
350
350
(M)
.6
.6
.4
.4
.4
.4
.3
.3
.5
.5
.5
.5
150
11
150
51
150
12
150
-23
150
24
150
-21
(M)
.7
.7
.7
.9
.7
.9
.7
.1
.7
.8
.7
.2
71
71
63
64
72
73
81
83
65
64
60
63
(M)
.3
.2
.9
.3
.3
.2
.7
.6
.6
.7
.8
.0
(M) (M
39.
39.
39.
39.
39.
39.
39.
39.
39.
39.
39.
3?.
5
5
5
5
5
5
5
5
5
5
5
5
57.0
71.2
51.9
64.3
57.7
73.2
64.4
83.6
53.1
64.7
^9.7
63.0
EFFECTIVE TOTAL
f SIG-Z CONC
(M) (US/M**3)
39.5 3.1797E-01
39.5 O.OOOOE-01
3.1797E-01
39.5 4.7911E-04
39.5 O.OOOOE-01
4.7911E-04
39.5 5.2610E-09
39.5 O.OOOOE-01
5.2610E-09
39.5 O.OOOOE-01
39.5 O.OOOOE-01
O.OOOOE-01
39.5 2.4281E-04
39.5 O.OOOOE-01
2.4281E-04
39.5 O.OOOOE-01
39.5 O.OOOOE-01
O.OOOOE-01
Figure B-8. (Page 7 of 14).
B-15
-------
CTDMPLUS VERSION 1.0
CTDMPLUS TEST RUN
PAGE
SRC-RECP SRC-RECP RECEPTOR EFF. FLAT
L DISTANCE DISTANCE HT ABOVE SRC-RECP TERRAIN
REC / ALONG FLOU CROSS FLOW STK BASE KT DIFF SIG-Y SIG-Z
* V (M) (M) (M) (M) (M) (M)
HILL-INDUCED
EFFECTIVE TOTAL
sic-r SIG-Z CONC
(M) (M) (US/M**3)
7 L 740.
7 V 788.
7 TOTAL
8 L 1017.
8 V 1059.
8 TOTAL
1205.0
86.0
1008.5
86.0
356.6
356.6
386.5
386.5
150.7
-27.3
150.7
-57.2
65.0
69.0
88.2
91.6
39.5
39.5
39.5
39.5
52.7
69.0
69.3
91.6
39.5 O.OOOOE-01
39.5 O.OOOOE-01
O.OOOOE-01
39.6 O.OOOOE-01
39.5 O.OOOOE-01
O.OOOOE-01
< SOURCE INFORMATION >
SOURCE QS TS VS BUOY FLUX MOM FLUX
$ (G/S) (K) (M/S) (M4/S3) (M4/S2)
FINAL PLUME
RISE
(M)
1.0 420.1 25.06
VARIABLES AT
PLUME HEIGHT:
HEIGHT
(M)
199.5
VDIR
(DEC)
1262.7
USCAL UVECT
(M/S) (M/S)
140.94
SIGV
(M/S)
SIGV DTHDZ
(M/S) (DEG/M)
330.9 300. 3.90 3.89 0.34 0.0390 0.0098
Figure B-8. (Page 8 of 14).
B-16
-------
CTDMPLUS VERSION 1.0
CTDMPLUS TEST RUN
PAGE
INFORMATION FOR HILL 1: PIEDMONT HILL
HCRIT = 178.5 M; FROUDE ft ABOVE HCRIT = 1.00
WRAP INFORMATION:
DISTANCE FROM SOURCE TO HILL CENTER = 2210.0 M;
ELLIPSE AXIS LENGTHS: MAJOR = 2187.9 M; MINOR =
MAJOR AXIS AZIMUTH FROM NORTH = 180.0 DEC ,
DISTANCE TO PRIMARY IMPINGEMENT POINT = 289-3 M
WRAP HT = 178.5 M
: 834.A M
LIFT INFORMATION:
DISTANCE FROM SOURCE TO HILL CENTER = 2211.9 M; LIFT MIDPOINT HT
ELLIPSE AXIS LENGTHS: MAJOR = 1528.0 M;MINOR = 554.5 M
MAJOR AXIS AZIMUTH FROM NORTH = 171.6 DEC
DISTANCE ALONG FLOU FROM SOURCE TO HILL CENTER = 1596.8 M
CROSSFLOU DISTANCE FROM SOURCE TO HILL CENTER = 1530.7 M
DISTANCE TO PRIMARY IMPINGEMENT POINT = 289.2 M
287.4
SRC-RECP SRC-RECP RECEPTOR EFF. FLAT HIL!
L DISTANCE DISTANCE- HT ABOVE SRC-RECP TERRAIN EF!
REC / ALONG FLOW CROSS FLOW STK BASE HT DIFF SIG-Y SIG-Z SIG
# V (M)
1 L 830.
1 V 830.
1 TOTAL
2 L 741.
2 V 746.
2 TOTAL
3 L 841.
3 V 852.
3 TOTAL
4 L 954.
4 V 974.
4 TOTAL
5 L 762.
5 V 755.
5 TOTAL
6 L 704.
6V 727.
6 TOTAL
(M)
62.
94.
215.
94.
373.
94.
641.
94.
-177.
94.
705.
94.
6
1
0
1
5
1
0
1
3
1
3
1
317
317
277
277
316
316
352
352
304
304
350
350
(M)
.6
.6
.4
.4
.4
.4
.3
.3
.5
.5
.5
.5
152
13
152
53
152
14
152
-21
152
26
152
-19
(M)
.4
.3
.4
.6
.4
.6
.4
.4
.4
.4
.4
.6
72
72
65
65
73
74
83
84
66
66
62
63
(M)
.6
.6
.1
.6
.5
.4
.0
.6
.9
.3
-0
.9
(M) (M
40.
40.
40.
40.
40.
40.
40.
40.
40.
40.
40.
40.
3
3
3
3
3
3
3
3
3
3
3
3
58.6
72.6
53.3
65.6
59.2
74.4
66.1
84.6
54.5
66.3
51.1
63.9
EFFECTIVE TOTAL
y SIG-Z CONG
(M) (US/M**3)
40.3 2.2282E-01
40.3 O.OOOOE-01
2.2282E-01
40.3 1.1321E-04
40.3 O.OOOOE-01
1.1321E-04
40.3 9.2630E-10
40.3 O.OOOOE-01
9.2630E-10
40.4 O.OOOOE-01
40.3 O.OOOOE-01
O.OOOOE-01
40.3 1.8925E-03
40.3 O.OOOOE-01
1.8925E-03
40.3 O.OOOOE-01
40.3 O.OOOOE-01
O.OOOOE-01
Figure B-8. (Page 9 of 14).
B-17
-------
CTDMPLUS VERSION 1.0
CTDMPLUS TEST RUN
PAGE
10
SRC-RECP SRC-RECP RECEPTOR EFF. FLAT
L DISTANCE DISTANCE HT ABOVE SRC-RECP TERRAIN
REC / ALONG FLOU CROSS FLOW STK BASE HT DIFF SIG-'f SIG-Z
if V (M) (M) (M) (M) (M) (M)
7 L 755.
7 V 794.
7 TOTAL
8 L 1032.
8 W 1066.
8 TOTAL
HILL-INDUCED
EFFECTIVE TOTAL
SIG-Y SIG-Z COMC
(M) (M) (US/M**3)
1231.0
94.1
1034.5
94.1
356.6
356.6
386.5
386.5
152.4
-25.7
152.4
-55.5
MAXIMUM CONCENTRATION FOR THIS HOUR IS
66.3 40.3 54.1 40.3 O.OOOOE-01
69.5 40.3 69.5 40.3 O.OOOOE-01
O.OOOOE-01
89.4 40.4 71.1 40.4 O.OOOOE-01
92.2 40.4 92.2 40.4 O.OOOOE-01
O.OOOOE-01
0.54 US/M**3 AT RECEPTOR 3 1
SOURCE CONTRIBUTION TABLE
DATE 6/26/80 HOUR 1
REC I
1
2
3
4
. 5
6
7
8
STACK 1 1
0.3180E+00
0.4791E-03
0-5261E-08
O.OOOOE+00
0.2428E-03
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
STACK S 2
0.2228E+00
0.1132E-03
0.9263E-09
O.OOOOE+00
0.1893E-02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
TOTAL
0.5408E+00
0.5923E-03
0.6187E-08
O.OOOOE+00
0.2135E-02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
Figure B-8. (Page 10 of 14).
B-18
-------
CTDMPLUS VERSION 1.0
CTDMPLUS TEST RUN
PAGE 11
INPUT MET DATA FROM SURFACE AND PROFILE (NOTE: ****** = MISSING DATA):
TO MO DA HR
80 6 26 10
MIXING
HEIGHT
(M)
SFC
TEMP
(K)
SFC
U*
(M/S)
MONIN-
OBUKHOV
LENGTH
(M)
SFC
ROUGH.
LENGTH
(M)
1242.0 299.4 0.293
-7.4 0.1500
ADJUSTED WIND <-WIND SPEED-> AMB. SIGMA-
HEIGHT DIR. SCALAR VECTOR TEMP THETA SIGMA-V SIGMA-W
(M) (DEC) (M/S) (M/S) (K) (DEC) (M/S) (M/S)
11.5 328.0 2.20 -999.90 299.40 -999.0 -999.90 0.49
101.5 328.0 2.50 -999.90 299.40 -999.0 -999.90 0.49
NOTE: SCALAR WIND SPEEDS USED IN CTDM ARE SET TO A MINIMUM OF 1 M/S
NOTE: HEIGHTS ARE REFERENCED TO THE COMMON STACK BASE ELEVATION
THE ADJUSTMENT TO THE INPUT HEIGHT IS 1.5 METERS.
QS
(G/S)
1.0
TS
(K)
410.1
VS
(M/S)
25.06
j.1* f wi\i in i. j. VLI
BUOY FLUX
(M**4/S**3)
187.3
SOURCE
# '
1
SOURCE INFORMATION ADJUSTED FOR PENETRATION:
SOURCE QS BUOY FLUX
ft (G/S) (M**4/S**3)
1 1.0 187.3
PLUME LAYER AVERAGE VARIABLES:
HEIGHT WD WS
367.6 330.4 2.53
CONVECTIVE SCALING PARAMETERS:
USTAR WSTAR FSTAR
0.2930 2.1935 0.0124
MOM FLUX
(M**4/S**3)
1293.9
PENETR
FACTOR
0.0000
Figure B-8. (Page 11 of 14).
B-19
-------
CTDMPLUS VERSION 1.0
CTDMPLUS TEST RUN
PAGE
REC
4
1
2
3
4
5
6
7
DVN-WND
DIST
(M)
698.
706.
880.
1137.
519.
965.
1286.
8 1444.
SOURCE
ft
2
5
4
3-
6
4
7
9
9
CRS-VND
DIST
(M)
-419.
-231.
-132.
82.
-590.
284.
731.
452.
s
8
3
5
7
5
0
0
7
TRANS
PL HT
(M)
284
286
331
0
233
0
0
0
QS TS
(G/S) (K)
1.0 420.1
.3
.5
.8
.0
.4
.0
.0
SIGMA HDF
V
(M) (1/H)
273
275
323
384
218
346
415
.0 446
- SOURCE
VS
(M/S>
25.06
.5 1.3593E-05
.8 8.3686E-05
.9 2.9155E-04
.1 6.8837E-04
.3 2.6320E-08
.0 9.4806E-04
.3 9.1530E-04
.6 8.8439E-04
INFORMATION _
BUOY FLUX
(M**4/S**3)
199.3
CVIC
(G/M**3)
6.
8.
9.
7.
1.
8.
6.
5.
8965E-04
6034E-04
1427E-04
0466E-04
8547E-05
3413E-04
3467E-04
6739E-04
s
MOM FLUX
(M**4/S**3)
1263.1
CONC
(US/M**3)
9.3745E-03
7.1998E-02
2.6656E-01
4.8507E-01
4.8816E-07
. 7.9081E-01
5.8091E-01
5.0179E-01
PENETR
FACTOR
0.0000
SOURCE INFORMATION ADJUSTED FOR PENETRATION:
SOURCE QS BUOY FLUX
# (G/S) (M**4/S**3)
2 1.0 199.3
PLUME LAYER AVERAGE VARIABLES:
HEIGHT VD VS
373.1 330.5 2.53
CONVECTIVE SCALING PARAMETERS:
USTAR VSTAR FSTAR
0.2930 2.1935 0.0132
Figure B-8. (Page 12 of 14).
B-20
-------
CTDMPLUS VERSION 1.0
CTDMPLUS TEST RUN
PAGE 13
EC
i
1
2
3
4
5
6
7
3
DWN-UND
01 ST
(M)
724.2
732.3
906.3
1163.8
545.0
992.1
1313.6
1471.4
CRS-VND
DIST
(M)
-405.6
-217.1
-118.5
96.5
-576.2
298.0
744.7
466.3
TRANS
PL HT
(H)
297.4
299.6
345.4
0.0
246.1
0.0
0.0
0.0
SIGMA
Y
(M)
281.6
283.9
331.5
391.0
227.1
353.4
421.9
452.8
HDF
(1/M)
1.7722E-05
9.7071E-05
3.0598E-04
6.8578E-04
6.3662E-08
9.3428E-04
9.0485E-04
8.7274E-04
CWIC
(G/M**3)
6.7111E-04
8.2709E-04
8.8877E-04
6.9192E-04
2.1805E-05
8.1337E-04
6.2484E-04
5.6083E-04
CONC
(US/M**3)
1.1894E-02
8.0287E-02
2.7194E-01
4.7450E-01
1.3882E-06
7.5992E-01
5.6539E-01
4.8946E-01
MAXIMUM CONCENTRATION FOR THIS HOUR IS
1.55 US/M**3 AT RECEPTOR #
SOURCE CONTRIBUTION TABLE
DATE 6/26/80 HOUR 10
REC #
1
2
3
4
5
6
7
8
STACK # 1
0.9374E-02
0.7200E-01
0.2666E+00
0.4851E+00
0.4882E-06
0.7908E+00
0.5809E+00
0.5018E+00
STACK # 2
0.1189E-01
0.8029E-01
0.2719E+00
0.4745E+00
0.1388E-05
0.7599E+00
0.5654E+00
0.4895E+00
TOTAL
0.2127E-01
0.1523E+00
0.5385E+00
0.9596E+00
0.1876E-05
0.1551E+01
0.1146E+01
0.9913E+00
Figure B-8. (Page 13 of 14).
B-21
-------
CTDMPLUS VERSION 1.0
CTDMPLUS TEST RUN
SEC TOP 4 CONCENTRATIONS [US/M**3] AT EACH RECEPTOR
* HIGHEST (JCD.HR) SECOND (JCD.HR) THIRD (JCD.HR)
PAGE
FOURTH (JCD.HR)
i
2
3
4
5
6
7
8
5.41E-01 (178, 1)
1.52E-01 (178,10)
5.38E-01 (178,10)
9.60E-01 (178,10)
2.UE-03 (178, 1)
1.55E+00*(178,10)
1.15E+00 (178,10)
9.91E-01 (178,10)
2.
5,
6.
0.
1.
0.
0.
0.
132-02* (173,
92E-04
19E-09
OOE-01
88E-06
OOE-01
OOE-01
OOE-01
(178,
(178,
(178,
(178,
(178,
(178,
(178,
10)
1)
1)
1)
10)
1)
1)
1)
*********(
******** (
******** (
******** (
******** (
******** (
******** (
******** (
o,
0,
0,
0,
0,
0,
0,
0,
0)
0)
0)
0)
0)
0)
0)
0)
*********(
******** (
******** (
******** (
******** (
******** (
******** (
******** (
0,
0.
0,
0,
0,
0,
0,
0.
0)
0)
0)
0)
0)
0)
0)
0)
Figure B-8. (Page 14 of 14).
B-22
-------
-------
United States
Environmental Protection
Agency
Center for Environmental Research
Information
Cincinnati OH 45268
Official Business
Penalty for Private Use, $300
Please make all necessary changes on the above label,
detach or copy, and return to the address in the upper
left-hand corner.
II you do not wish to receive these reports CHECK HERE D,
detach, or copy this cover, and return to the address in the
upper left-hand corner.
EPA/600/8-89/041
------- |