I 600385064
(7
GREEN RIVER AIR QUALITY MODEL DEVELOPMENT
VALMET - A Valley Air Pollution Model
ATMOSPHERIC SCIENCES RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
RESEARCH TRIANGLE PARK, NORTH CAROLINA 27711
-------
GREEN RIVER AIR QUALITY MODEL DEVELOPMENT
VALMET - A Valley Air Pollution Model
by
C. D. Whiteman, and K. J. Allwine
Pacific Northwest Laboratory
Rich!and, Washington 99352
Interagency Agreement AD89F20970
with the U.S. Department of Energy
Project Officer
Alan H. Huber
Meteorology and Assessment Division
Atmospheric Sciences Research Laboratory
U.S. Environmental Protection Agency
Research Triangle Park, NC 27711
ATMOSPHERIC SCIENCES RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
RESEARCH TRIANGLE PARK, NORTH CAROLINA 27711
-------
DISCLAIMER
The research described in this report has been funded wholly or in part
by the United States Environmental Protection Agency through Interagency
Agreement AD-89-F-2-097-0 to the Pacific Northwest Laboratory. This docu-
ment has been reviewed in accordance with U.S. Environmental Protection
Agency policy and approved for publication. Mention of trade names or
commercial products does not constitute endorsement or recommendation for
use.
-------
PREFACE
This final report is submitted as part of the Green River Ambient Model
Assessment (GRAMA) project conducted at the U.S. Department of Energy's
Pacific Northwest Laboratory for the U.S. Environmental Protection Agency.
The GRAMA Program has, as its ultimate goal, the development of validated
air quality models that can be applied to the complex terrain of the Green
River Formation of western Colorado, eastern Utah and southern Wyoming. The
Green River Formation is a geologic formation containing large reserves of
oil shale, coal, and other natural resources. Development of these
resources may lead to a degradation of the air quality of the region. Air
quality models are needed immediately for planning and regulatory purposes
to assess the magnitude of these regional impacts. This report documents
one of the models being developed for this purpose within GRAMA—
specifically a model developed to predict worst case air pollutant
concentrations caused by elevated continuous point sources of pollution
located in deep valley terrain. The model shows promise for use as a
planning tool and eventually as a regulatory tool. Further testing and
evaluations of the model are needed to gain a measure of confidence in the
model's performance.
-------
ABSTRACT
An air quality model is described for predicting air pollution con-
centrations in deep mountain valleys arising from nocturnal down-valley
transport and diffusion of an elevated pollutant plume, and the fumigation
of the plume on the valley floor and sidewalls after sunrise. Included is a
technical description of the model, a discussion of the model's applica-
tions, the required model inputs, sample calculations and model outputs, and
a full listing of the FORTRAN computer program.
IV
-------
EXECUTIVE SUMMARY
Pacific Northwest Laboratory has developed an air quality model for
application in valleys for the U.S. Environmental Protection Agency (EPA)
under their Green River Ambient Model Assessment (GRAMA) program. This
program was initiated in response to the need for air quality assessment
tools applicable in the Green River Oil Shale Formation region of western
Colorado, eastern Utah, and southern Wyoming. This region has the potential
for large-scale growth because vast energy resources, especially oil shale,
are located in the region.
Following a thorough analysis of meteorological data obtained from deep
valleys of western Colorado, a modular air pollution model has been
developed to simulate the transport and diffusion of pollutants released
from an elevated.point source in a well-defined mountain valley during the
nighttime and morning transition periods. This initial version of the
model, named VALMET, operates on a valley cross section at an arbitrary
distance down-valley from a continuous point source. The model has been
constructed to include parameterizations of the major physical processes
that act to disperse pollution during these time periods.
The nighttime simulation uses a modified Gaussian plume to calculate
concentrations on the valley floor and sidewalls. The Gaussian plume model
uses Briggs1 plume rise equations and Pasquill-Gifford diffusion coeffi-
cients, as modified for complex terrain diffusion enhancement. Dilution of
the plume due to tributary flows can be handled in the model, and a new
method has been incorporated to ensure conservation of plume mass as the
plume is transported down the valley axis. The purpose of the nighttime
simulation is to provide air pollution concentrations on a cross-valley
section at sunrise, as an initial condition for the daytime simulation.
The daytime simulation uses numerical techniques and a valley energy
budget to predict how concentrations will vary in time as pollutants from
the elevated plume are fumigated onto the valley floor and sidewalls during
the post-sunrise temperature inversion breakup period. The effects of CBL
growth, temperature inversion subsidence, transport and diffusion in uoslope
flows, warm air advection above the valley, albedo, sensible heat flux, and
other physical processes are incoroorated in the model. The daytime
simulation is driven by solar radiation and a simplified surface energy
budget, taking account of valley topography.
Valley meteorological observations are required in order to obtain the
necessary input parameters to the model. The required observations can be
obtained using standard meteorological instrumentation, and the costs for
-------
the observations necessary to perform a simple screening analysis are
expected to be within the means of a small industrial applicant wanting a
permit for operation in a mountain valley. The VALMET model could form the
basis for a more comprehensive site-specific model if further meteorological
and tracer data were available in a given valley. The model outputs are in
a form suitable for regulatory decision making. The model predicts the
maximum 1- and 3-hr-average concentrations for locations on the valley floor
and sidewalls.
The model shows promise for use as a planning tool and, eventually as a
regulatory tool. Further development, testing, and tracer evaluations of
the model will be necessary before sufficient confidence can be gained to
justify the model's use in a regulatory setting. The priorities for further
development and testing have been provided in the body of the report. Test-
ing of the model's sensitivity to input parameters, and an initial evalua-
tion of the model with tracer experiment data are high priority tasks.
These tests will, no doubt, result in future modifications to the initial
version of the model.
The authors stress that the model's ultimate utility in addressing and
providing solutions to potential air pollution problems in mountain valleys
will depend on the further evaluation of the model. In order to have
confidence in model predictions it is necessary to test the model against
actual air pollution data. Several parameters in the model (AQ, k, a , and
a ) are, at present, poorly understood for mountain valleys due to a aearth
of experimental data and theoretical research. We hope, by pointing out
these deficiencies, that attention will be focused on the need for
information on both turbulent diffusion and valley energy budget studies.
The use of full physics models may help in providing some of the answers
necessary to improve the present model.
-------
CONTENTS
Preface i i i
Abstract 1v
Executi ve Summary v
Figures ix
Tab! es xi i
Acknowledgments xiii
1. Introduction 1
2. Background 4
Nocturnal Val 1 ey Meteorol ogy 4
Valley Meteorology During the Temperature Inversion
Breakup Period 5
3. Technical Discussion 14
Coordinate System 14
Topographic Cross Section 14
Nocturnal Model 17
Pollutant Source 17
PIume Rise 18
Gaussi an PI ume Model 21
Dispersion Coefficients 22
Channel ing 22
PI ume Di 1 uti on 24
Reflections 26
Calculation of Steady-State Nighttime
Concentrations 31
Daytime Model 32
Extraterrestrial Solar Radiation 33
Parameterization of Sensible Heat Flux 36
Model Grid 36
Thermodynamic Equations for CBL Ascent and
Inversion Descent 38
Advection in the Slope Flows 42
Pollution Concentration Calculation Method 44
Exponential Decay of Concentration 47
Maximum 1- and 3-Hour Average Concentrations 48
4. Overvi ew of Modul ar VALMET Model 49
Features of the Computer Code 50
Model Inputs 50
Specification of the Model Inputs by the User 53
Valley Characteristics 53
Date 55
Site Location 56
-------
Model Characteristics.. 56
Val 1 ey Atmosphere 57
Stack Characteristics. 57
Inversion Characteristics... 58
Gaussian Plume Parameters 58
Sensible Heat Flux 60
Error Messages 62
5. Technical Description of Individual Modules.. 64
VALMET-Mai n Program 64
INPUT-Input Module 69
JULIAN-Julian Day Module 71
PRISE-PIume Rise Module.... 72
DILUTE-PIume Dilution Module 74
INGRAT-Valley Plume Reflection Module 75
NORMAL-Normal or Gaussian Curve Integration Module 77
GAUSS-Gaussi an PI ume Modul e 78
SOLAR-Extraterrestrial Solar Radiation Module 79
EBDGT-Surf ace Energy Budget Modul e 81
DESCNT-Inversion Descent and CBL Growth Module 82
PROFIL-Concentrati on Prof i 1 e Modul e 84
VELOCY-Upsl ope Fl ow Vel oci ty Modul e 86
BRKUP-Pol 1 utant Mass Budget Modul e 88
PSTPRC-Post Processor Modul e 90
6. Samp!e Model Runs 92
Simulation Number 1 92
Simulation Number 2 97
7. Further Work 100
Guidance from the 1982 Experiments 100
Suggested Modifications to the VALMET Model 103
Deposi ti on 103
Emission Above or Below Stable Core 103
Energy Budget 104
Cross Vail ey Fl ows 104
Turbulent Erosion of the Valley Inversion by
Overlyi ng Fl ows 104
Effect of Tributary Flows on the Enhancement of
Di f f usi on 105
Diffusion Coefficients 105
Time-Varying Wind Speeds in the Stable Core 105
Temperature Inversion Buildup 105
Differential Heating 105
Stacks Located on SidewalIs 106
References 109
Appendices
A. FORTRAN Listing of VALMET 115
B. FORTRAN Listing of VALMET Output Plotting Program 147
C. Research Paper Entitled "Breakup of Temperature Inversions
in Deep Mountain Valleys: Part II. Thermodynamic Model".. 155
D. Summary of Modifications to VALMET 173
vn
-------
FIGURES
Number Page
1 Three patterns of temperature structure evolution during the
i nversi on breakup peri od 7
2 Typical mid-morning wind structure over and within a deep
valley on the western slope of the Rockies, illustrating
the five interrelated wind systems identified in field
studi es 8
3 Relationship between temperature structure layers and wind
systems 9
4 Dual soundings from a valley floor and a valley sidewall
illustrating the up-slope flow found within the CBL over
the si dewal 1 9
5 Illustration of the hypothesis of inversion destruction.
In the center of the diagram cross sections of a valley
are shown at times t, , t^, t,, t^, and tg 10
6 Air pollution implications of CBL growth and inversion
top descent 12
7 Illustration of the nocturnal down-valley transport and
dispersion of a pollutant plume 15
8 Illustration of the local coordinate system on a cross
val 1 ey secti on 16
9 Parameters used in the VALMET model to approximate a
valley topographic cross section 16
10 Illustration of the definitions of secondary topographic
variables used in the model 16
11 Cross section of pollutant plume and valley topography
illustrating the integral method of calculation for
plume reflection and diffusion out the top of the
val 1 ey temperature i nversi on 28
ix
-------
Number Page
12 VALMET grid configuration on a valley cross section
illustrating the nomenclature used in the model 32
13 Extraterrestrial solar heat flux as a function of time,
showing solar model nomenclature 37
14 Cross section of the valley floor and sidewalls illustrating
the grid elements whose height corresponds to the CBL
heights 38
15 Illustration of the effect of topography in controlling the
heating rates of the air within a valley temperature
inversion versus the air within an inversion over the
pi ai ns 39
16 Diagram showing the changes in CBL and inversion depth at a
given time step above a valley floor grid element 43
17 Representation of the volumetric element of mass incorporated
into the growing CBL above a valley floor grid element at
each model time step 43
18 Schematic diagram of an individual model grid element
illustrating the pollutant mass balance 45
19 Flow diagram of VALMET model, showing the modular structure
of the model and the physical processes parameterized
in the modules 51
20 Illustration of the VALMET model input table as it appears
on the user's interactive screen 54
21 Listing of summary output file generated by Sample
Simulation 1 93
22 Plots of Sample Simulation 1: (a) nocturnal vertical
concentration profile through plume centerline, (b) CBL
height (long dashes) and inversion height (short dashes)
as a function of time, and (c) nocturnal cross-valley
concentration profile through plume centerline 95
23 Pollutant concentration versus time for selected grid
elements for Sample Simulation 1 96
24 Same as Figure 22, for Sample Simulation 2 98
25 Same as Figure 23, for Sample Simulation 2 99
-------
Number Page
26 Illustration of differential solar flux on opposing sidewalls
for Brush Creek Colorado on August 4, 1982 107
27 Brehm's (1981) conceptual model of temperature inversion
destruction in Austria's Inn Valley, showing differential
CBL growth over the opposite sidewalls and continued
down-valley flow in the elevated stable core 108
XI
-------
TABLES
Number
1 Values of the Constants I, J, and K, for ay as a Function
of Downwind Distance, for Six Stability Conditions 23
2 Values of the Constants I, J, and K, for az as a Function
of Downwind Distance, for Six Stability Conditions 23
3 Relationship Between Weather Conditions and Stability
Categori es 24
4 Equation of Time Correction 35
5 Default values of VALMET Input Parameters 55
6 Limits on the Values of Input Parameters. 56
-------
ACKNOWLEDGMENTS
The work reported here had its origins in research conducted at the
Department of Atmospheric Science at Colorado State University while one of
the authors (COW) was a student there. Dr. Thomas B. McKee's contributions
to the work at that time, as a thesis advisor, are greatly appreciated.
The authors' colleagues at PNL contributed suggestions and criticisms
as the work progressed that led to significant improvements in the model's
design and in the clarity of the text. Drs. Thomas W. Horst and
J. Christopher Doran are especially thanked for their help in that regard.
Mr. Alan Huber, the EPA Project Officer, and Mr. Rich Fisher, EPA
Region VIII meteorologist, are thanked for their guidance, support, and
encouragement of the work. Mr. Fisher tested early versions of the model
and provided a number of very useful suggestions from the viewpoint of a
model user that resulted in improvements of the model code. The final set
of technical revisions to this report was made following suggestions by
Mr. Bill Bernardo of EPA's Region VIII office, arising from sensitivity
tests he conducted on the model in late 1983. We wish to thank Mr. Bernardo
for his help.
xi n
-------
SECTION 1
INTRODUCTION
This report documents an air quality model that was developed to pre-
dict concentrations of nonreactive pollutants arising from elevated con-
tinuous point sources that emit pollutants within well-defined deep mountain
valleys. The model, termed VALMET, is intended to simulate the effects on
pollutant transport and diffusion of various meteorological processes that
are thought to result in worst-case pollutant concentrations. The model is
run for situations when pollutants are carried in locally developed circula-
tions within a valley when these circulations are "decoupled" from prevail-
ing circulations above the valley. The primary physical processes included
in the model are
Nocturnal Simulation:
Transport by down-valley drainage flows
Plume channeling within the valley
Enhanced horizontal and vertical diffusion due to topography
Plume reflections off valley floor and sidewalls
Pollutant diffusion out the top of the valley
o Dilution of the plume due to clean air inflow from tributaries.
Post-Sunrise Simulation During Temperature Inversion Breakup Period:
• Convective boundary layer growth
o Plume subsidence in the valley inversion
• Fumigation into growing convective boundary layers
• Transport and diffusion in upslope flows over the sidewalls.
The model, while including a variety of meteorological processes, is
highly parameterized so that it is simple in concept and easy to run. The
model is composed of 13 modules, or subroutines, arranged in such a way that
an improved understanding of individual valley meteorological phenomena can
be easily incorporated in future versions of the model. The modules within
the model can be replaced by data if they are available. Thus, the model
can be used in one of two modes. It can be used in a "screening" mode to
calculate pollutant concentrations within a valley when little site specific
data is available, or the model can be "calibrated" with site-specific data
so that it can be used as a site-specific model.
-------
The two-dimensional model was developed primarily to predict pollutant
concentrations on the valley floor and sidewalls on a valley cross section
an arbitrary distance down-valley from a pollutant source during the post-
sunrise temperature inversion breakup period. It is necessary, however, to
know the air pollution concentration within the valley cross section at sun-
rise, as an initial condition for the post-sunrise simulation. The model is
therefore comprised of two parts—a nighttime part whose purpose is to pre-
dict concentrations on the valley cross section at sunrise, and the daytime
part which predicts concentrations on the valley floor and sidewalls during
the post-sunrise temperature inversion breakup period. The temperature
inversion breakup period has been identified by previous investigators [1]
as a period when diurnal fumigations [2] can produce high pollutant con-
centrations in valleys.
The nighttime simulation, which applies during the steady-state period
after valley temperature inversions and drainage wind systems have become
established, uses a modified valley-following Gaussian plume algorithm to
calculate air pollution concentrations for points on the valley floor and
sidewalls. A plume rise formulation is used to simulate the initial rise of
a pollutant plume at the stack due to momentum and buoyancy of the efflu-
ent. Pasqui11-Gifford diffusion coefficients are modified to account for
enhanced nocturnal diffusion caused by rough terrain. The Gaussian plume is
also modified to allow for dilution of the plume during its down-valley
transport caused by clean air flowing into the plume from valley tributaries
or by converging downslope drainage flows. An integral constraint on pol-
lutant mass is applied to ensure that pollutant mass is conserved during the
plume's transport down the valley and within any valley cross section down-
valley from the emission source, except for pollution diffusion out the top
of the valley.
The daytime simulation uses numerical techniques that simulate the
fumigation of the nocturnal plume onto the valley floor and sidewalls as a
convective boundary layer grows upwards from the heated valley surfaces and
as subsiding motions occur over the valley center after sunrise. The rate
of growth of convective boundary layers and subsidence within the valley
temperature inversion are simulated using the bulk thermodynamic model of
Whiteman and McKee [3]. This model is driven by sensible heat flux, esti-
mated as a fraction of the solar radiation using a highly parameterized
surface energy budget. The effects of such factors as snow cover, soil
moisture, cloud cover, or surface albedo are not explicitly included in the
model but can be incorporated into the model in the future through an
expanded energy budget module. The shape of the topographic cross section
of the valley is explicitly included in the model through the specification
of valley floor width and sidewall inclination angles at the valley cross
section of interest. The retarding effect on temperature inversion breakup
and pollution dispersion due to warm air advection above the inversion is
also included in the model. Fumigated pollutants are transported from the
valley cross section in upslope flows that develop within the convective
boundary layers over the slopes. Pollutants are diffused through model grid
elements during this transport up the slopes in the growing convective
-------
boundary layer. Pollutant concentrations decay exponentially within indi-
vidual grid elements high on the sidewall as they are dropped from the simu-
lation as the inversion top subsides below them.
The output from the nighttime simulation includes the steady-state
pollutant concentration at valley floor and sidewall grid elements on the
valley cross section of interest. The fraction of plume mass that has dif-
fused out the top of the valley during the plume's travel is also an output
of the model. Since an analytical formula describes the concentrations
within a valley cross section, cross valley and vertical profiles of pollu-
tant concentration can be calculated and plotted. The plume centerline con-
centration is an output of the model.
The primary output of the daytime simulation is the maximum 1- and 3-hr
average pollutant concentrations in each of the model grid elements on the
valley floor and sidewalls. The time varying 5-min-average concentrations
for each of the grid elements between sunrise and the time of inversion
destruction is also an output of the model. In addition to these primary
outputs, intermediate model outputs come from individual modules in the
program. The local standard time of sunrise, the duration of the daylight
period, and the solar flux on a horizontal surface at solar noon come from
the solar module. The convective boundary layer height and inversion top
height as a function of time come from the temperature inversion breakup
modu1e.
Twenty-seven input parameters are necessary to drive the model. These
input parameters include the date, site location, topographic characteris-
tics of the valley cross section, temperature inversion characteristics at
sunrise, emission and stack characteristics, down-valley wind speed(s),
atmospheric stability, grid element length, and sensible heat flux para-
meters. If known, the rate of warm air advection above the valley can be
input. The necessary model inputs can be obtained from topographic maps,
engineering information on the pollutant source, and one or more seasonal
meteorological data collection campaigns in the valley of interest using
tethered balloon data collection systems and/or doppler acoustic sounders.
This report is written following U.S. Environmental Protection Agency
(EPA) reporting guidelines [4,5] and is organized in the following way.
First, meteorological observations taken in a number of deep valleys in
western Colorado are summarized in order to provide a brief description of
the salient meteorological features that must be included in a realistic air
quality model. Second, a technical discussion is given of the mathematical
representation of these processes used in VALMET. Third, a brief overview
is given of the air quality model's structure. Fourth, the individual
modules of the model are described in detail, including the FORTRAN code,
inputs, outputs, etc. Fifth, sample runs of the model are presented.
Finally, the need for further data and verification studies is indicated and
an overall summary of the modeling approach is presented. Appendices
include a complete listing of the VALMET source code, a listing of a sample
plotting program for depicting model outputs, and a reprint of a scientific
paper which describes the theoretical basis of one of the major components
of the modeling approach.
-------
SECTION 2
BACKGROUND
The VALMET model is a further development of ideas that were published
in research papers and reports over the last 5 years, based on a comprehen-
sive set of valley meteorological field experiments conducted in western
Colorado. These earlier papers summarize a case study of temperature inver-
sion buildup [6]; describe the data analysis, hypothesis, and mathematical
foundation of the temperature inversion destruction part of the model
[3,7-10]; describe the initial idea of an air pollution model framework able
to incorporate upslope winds and fumigations arising after sunrise [11]; and
present a general summary of Colorado valley meteorology based on data
analyses [12]. Progress in the development of VALMET has been described in
two progress reports submitted to EPA in 1981 and 1982 [13,14].
The present version of VALMET, as described in this report, must be
considered as an initial or preliminary version of the model. This report
is written to facilitate the limited distribution of the model to users who
will provide comments and suggestions regarding the further development of
the model. The present version has not yet undergone a full set of sensi-
tivity tests, has not been fully evaluated with respect to its performance
in simulating air pollution concentrations in actual valleys, and contains
modules that have not yet been fully developed.
The user is cautioned against applying the air quality model for
meteorological or topographical conditions that are contrary to assumptions
made in the model. In order for the user to understand fully the assump-
tions and hypotheses used in the model, it is necessary to summarize briefly
the results and hypotheses that arose from the observational program and the
work of previous investigators. In this section we present a condensed sum-
mary of this material for both the nocturnal and temperature inversion
breakup periods, emphasizing the physical mechanisms in the valley that are
responsible for the transport and diffusion of pollutants. It is these
individual mechanisms that are parameterized in the air pollution model.
NOCTURNAL VALLEY METEOROLOGY
The most universal features of nighttime valley meteorology, reported
in valleys all over the globe [15-19, and others], is the drainage of cooled
air down the valley sidewalls and valley axes. The initiation of these
flows, which is closely related to changes in the temperature structure
above valley surfaces, has received inadequate study to date and cannot be
adequately parameterized for inclusion in the model. After the short (2 to
-------
4 hr) initiation period, the down-valley flows become well-established
through the valley depth, and attain a speed that is characteristic of the
valley and is reasonably steady-state. Observations [9,12] show that the
strength of the down-valley flows is quite variable from valley to valley,
and no model is yet available which can accurately predict their strengths
in individual valleys. In Colorado's valleys, however, observations show
that windspeeds on a clear undisturbed night in one season are typical of
windspeeds on similar nights in other seasons. Observations thus seem to be
the best way of specifying down-valley wind system strength.
Down-slope flows on the valley sidewalls become weaker and less
organized after the valley inversion develops, since air flowing down the
slopes must work against a very stable ambient stratification. Slope flows
can more easily become detached from the sidewalls during this period,
providing an important stirring mechanism which causes enhanced horizontal
diffusion of pollutant plumes within the valley inversion. The effects of
terrain roughness and mechanically generated turbulence also play an impor-
tant role in enhancing dispersion in the down-valley flows.
Pollutant concentrations can be strongly affected by dilution of the
plume due to clean air flowing in from valley tributaries. This feature of
meteorology, unique to valleys, can make substantial changes in pollutant
concentrations, especially where high volume fluxes of clean air are
involved, as when a pollutant plume originates in a small tributary and
emerges into a major valley.
VALLEY METEOROLOGY DURING THE TEMPERATURE INVERSION BREAKUP PERIOD
Observation in Colorado's deep valleys [9,10] have provided detail on
the changes in atmospheric structure during the morning transition period
when the nocturnal down-valley flows are reversed to daytime up-valley
flows. The transition period is much longer than expected in these deep
valleys and the physical processes leading to the transition are expected to
be of great importance for air pollution transport and dispersion. In this
section, we will summarize the observations of inversion destruction, point-
ing out typical characteristics of the meteorology of these valleys, and the
physical processes that must be included in realistic air pollution models
simulating the inversion destruction period.
At sunrise nocturnal inversions in most of the valleys were built up to
about the level of the surrounding ridgetops. The average depth, based on
21 cases studies, was 604 m. Vertical potential temperature gradients
within the inversions averaged 0.0295°K m" , but ranged from 0.0187 to
0.0566°K m" . The strength and direction of prevailing winds aloft, as
determined from the Grand Junction, Colorado, morning rawinsonde sounding,
had no demonstrable effect on the characteristics of valley inversions. The
valley inversions showed much less variability from day to day and from
season to season than inversions over Grand Junction. This suggests that at
least for a narrow range of synoptic conditions, valley topography produces
more consistent inversions, perhaps by protecting them from winds aloft.
-------
When the along-valley winds were weak, the valley vertical potential
temperature profiles frequently were hyperbolic, especially near the
ground. In most valleys where along-valley winds were of at least moderate
strength, the inversion profiles were nearly linear with height.
Temperature inversions in all of the valleys investigated were
destroyed after sunrise following one of three patterns of temperature
structure evolution (Figure 1). The first pattern, observed in the widest
valley studied, approximates inversion destruction over flat terrain, in
which the nocturnal inversion is destroyed after sunrise by the upward
growth from the ground of a warming convective boundary layer (CBL). The
second pattern, observed in snow covered valleys, differs significantly from
the first. Here the growth of the CBL, which begins after sunrise, is
arrested once the CBL has attained a depth of 25 to 50 m. The inversion is
then destroyed as the top of the nocturnal inversion descends into the
valley. Successive profiles of the valley atmosphere show a warming
consistent with a simple subsidence of the previous profiles. The third
pattern of temperature structure evolution was observed in all the valleys
when snow cover was not present and describes the majority of case studies
observed in field experiments. Following this pattern, inversions are
destroyed by two processes: the continuous upward growth from the valley
floor of a warming CBL and the continuous descent of the top of the
nocturnal temperature inversion. Warming of the elevated inversion layer
above the CBL is consistent with a simple subsidence of the previous
profiles. In the Colorado valleys studied, the time required to break an
inversion and establish a neutral atmosphere within the valley was typically
3-1/2 to 5 hr after sunrise. Temperature structure evolution during clear,
undisturbed weather was surprisingly uniform from day to day and season to
season. Thus, in future work, one may be fairly confident of observing
typical inversion breakup during short field studies in undisturbed weather.
The common element of all three patterns of temperature structure
evolution is the development of a CBL over the valley floor after it is
illuminated by direct sunlight. Observations taken from the sidewall of one
valley also show the development of a CBL after direct sunlight illuminates
the sidewall. Due to the shading effects of surrounding topography, the
different valley surfaces can be illuminated at significantly different
times, thus affecting the initiation of CBL growth. The temperature struc-
ture of the sidewall CBL is similar to that of the CBL over the valley
floor, but winds blow up the sidewall CBL at speeds of up to 3 m sec .
Five different temperature structure layers have been observed during
the inversion destruction period. Above the valley floor CBL and the
sidewall CBL just mentioned is the stable core of the potential temperature
inversion. A neutral stability layer above the stable core appears to be
part of a large-scale convective boundary layer that forms over the western
slope of the Rocky Mountains. Above this layer is the free atmosphere.
-------
zt
(a)
(b)
H(t)
e
t
Pattern 1. Growth of CBL.
zt
(0
(d)
\
h(t)
\
\
\
H(t)
8 tj t2 t3 t0 t
Pattern 2. Descent of inversion top and arrested growth of CBL.
Zt
to
(f)
h(t)
\
\
.H(t)
j i
9 t| tz t3 t0 t
Pattern 3. Descent of inversion top and continuous growth of CBL.
Figure 1. Three patterns of temperature structure evolution during the
inversion breakup period. Potential temperature profiles
are on the left, and time-height analyses of CBL height (H)
and inversion top height (h) are on the right.
Each of the five temperature structure layers, identified primarily by
their potential temperature structure, can also be identified by the winds
that prevail within them (see Figures 2, 3 and 4). During inversion
destruction, the CBLs over the valley floor and sidewalls contain winds
which blow up the floor of the valley and up the slopes. The CBL, or
-------
UP-FLOOR
A. DOWN-VALUEY
Figure 2. Typical mid-morning wind structure over and within a deep valley
on the western slope of the Rockies, illustrating the five
interrelated wind systems identified in field studies.
neutral layer, above the valley inversion has winds which blow up the
inclined Western Slope of the Rocky Mountains during the day. Winds in the
stable core typically continue to blow down-valley after sunrise until the
stable core is nearly destroyed. Winds in the stable free atmosphere may
blow from any direction with speeds determined by snyoptic-scale pressure
gradients. Despite variability in the strength and timing of reversal of
the winds, the temperature structure evolved uniformly from day to day in
individual valleys.
On the basis of the wind and temperature observations summarized above,
an hypothesis (Figure 5) has been developed to explain the temperature
structure evolution. Since energy is required to change the temperature
structure, and the changes begin at sunrise, it is reasonable to hypothesize
that solar radiation is the driving force. A fraction of the solar radia-
tion, received on the valley floor and sidewalls, is converted to the sen-
sible heat flux that provides energy to the valley atmosphere. Sensible
heat flux from a surface, as over flat terrain, causes a convective boundary
-------
t
FREE
ATMOSPHERE
NEUTRAL LAYER
STABLE
CORE
CBL
GRADIENT
DP-INCLINE
DOWN-VALLEY
UP-FLOOR
e
Figure 3. Relationship between temperature structure layers
and wind systems. The temperature profile is a
typical mid-morning sounding from the floor of a
deep valley.
Z f NEUTRAL
LAYER
STABLE
CORE
CBL
UP-SLOPE
ff
Figure 4. Dual soundings from a valley floor and a
valley sidewall illustrating the up-
slope flow found within the CBL over
the sidewal1.
-------
9
Figure 5. Illustration of the hypothesis of inversion destruction. In the
center of the diagram cross sections of a valley are shown at
h.
times t-i , t~, t.,, t., and tn. On the left are corresponding
potential temperature profiles as taken from the valley center.
On the right are corresponding up-valley wind components as a
function of height. At sunrise, tn-, an inversion is present in
the valley. At t?, a time after sunlight has illuminated the
valley floor and slopes, a growing CBL is present over the val-
ley surfaces. Mass and heat are entrained into the CBLs from
the stable core above and carried up the sidewalls in the
upslope flows. This results in a sinking of the stable core and
growth of the CBLs (t3 and t^) until the inversion is broken
(t,,) and a turbulent well-mixed neutral atmosphere prevails
through the valley depth. Down-valley winds continue to blow in
the stable core during the inversion breakup period. Winds in
the CBL below, and in the region above the stable core often
blow up-valley during this same period.
10
-------
layer to develop over the surface. Mass and heat are entrained into the CBL
from the stable core above. Mass entrained into the valley floor and side-
wall CBLs, however, is carried from the valley in the upslope flows that
develop in the convective boundary layers over the sidewalls. This removal
of mass from the base and sides of the stable core causes the elevated
inversion to sink deeper into the valley and to warm adiabatically due to
subsidence, and decreases the rate of growth of the CBLs. Following this
hypothesis, the rate of warming depends directly on the rate of energy input
into the valley atmosphere. This energy may be used to deepen the CBLs or
to move mass up the sidewalls, allowing the stable core to sink. From this
hypothesis a thermodynamic model of temperature inversion destruction has
been developed [3]. This model forms the basis for parameterizations of
inversion breakup in the VALMET model.
The valley inversion destruction mechanism has important implications
for the dispersion of valley air pollution. These implications can be
investigated by assuming that the nocturnal inversion at sunrise contains
pollutants emitted from a continuous elevated source within the valley.
During the night such pollutants would be carried down the valley in the
along-valley wind system, undergoing both vertical and horizontal disper-
sion. The horizontal dispersion in these flows is known to be much greater
than over flat terrain [20]. After sunrise a convective boundary layer
forms over the valley floor and sidewalls (Figure 6a). The subsequent dis-
position of pollutants in the stable core will be affected by two competing
processes—the sinking of the stable core and the growth of the CBL. The
air pollution implications of the two processes can be considered by study-
ing a valley cross section down-valley from the source at a later time. The
three inversion breakup patterns discussed above have the following air
pollution implications [8]:
1. Pattern 1 - Growth of convective boundary layer (Figure 6b): Pure
growth of a CBL will result in fumigation [2] of pollutants at the
valley floor as the CBL grows upward into the polluted stable core.
This process is favored when the slope flows are ineffective in
removing mass from a valley, and thus will occur in very wide or
shallow valleys.
2. Pattern 2 - Sinking of stable core (Figure 6c): Failure of the CBL to
grow once it has formed over the valley floor and sidewalls results in
inversion destruction by sinking of the stable core. Thus, the pol-
lutants sink into the top of a shallow mixed layer, producing high con-
centrations at the ground. The pollutant plume, once entrained into
the CBL, is advected up the sidewalls and dispersed into the neutral
layer aloft. This process is favored for narrow-to-wide valleys when
sensible heat flux is weak.
3. Pattern 3 - Combination (Figure 6d): A combination of CBL growth and
stable core descent results in the sinking of pollutants into the top
of a growing mixed layer. Pollution concentrations should be interme-
diate between the previous two cases. This pattern is the most common
one in Colorado mountain valleys.
11
-------
Figure 6. Air pollution implications of CBL growth and inversion
top descent.
Another process, that of cross-valley advection of the pollutant plume
due to differential heating of the valley sidewalls, may be important in
determining the cross-valley position of maximum pollutant concentration.
Unfortunately, relatively little observational evidence is presently avail-
able on cross-valley advection or the effects of cross valley advection on
pollutant dispersion.
Observations in snow-free valleys of Colorado show that inversions are
typically destroyed within 3-1/2 to 5 hr after sunrise, regardless of
season. Thus, high pollutant concentrations at the ground may be expected
following sunrise in polluted Colorado valleys. But inversions are typ-
ically broken every day and, once the inversion is destroyed, good mixing
will prevail in an up-valley flow regime that is coupled to flows above the
valley.
Moist surface conditions or high albedo due to snow cover may change
the surface energy budget components so that sensible heat flux is
reduced. In these conditions inversion destruction will be delayed or an
inversion may persist all day. In Colorado's snow-covered Yampa Valley on
February 23, 1978, the CBL failed to grow deeper than 35 m and the observa-
tions were characterized by a slow descent of the stable core. The inver-
sion failed to break during this day, although the top of the stable core
descended to within 150 m of the ground.
To summarize, the primary factors affecting a valley pollutant plume
after sunrise are
» continued transport of the pollutant plume down the valley in the
stable core
• sinking motions in the stable core as inversion destruction progresses
after sunrise, bringing the plume closer to the valley floor
12
-------
• fumigation of the plume into convective boundary layers growing over
the valley floor and sidewalls
• transport and diffusion of pollution in upslope flows that develop over
the sidewalls.
These processes continue from sunrise until the valley temperature inversion
is destroyed and the valley atmosphere becomes coupled with the overlying
circulations. The amount of time required for temperature inversion
destruction in Colorado Mountain Valleys is typically 3-1/2 to 5 hr, but it
depends on many factors including the depth and strength of the initial
temperature inversion, the incoming solar energy, the sensible heat flux,
the pattern of inversion destruction followed, etc.
13
-------
SECTION 3
TECHNICAL DISCUSSION
In this section we present a technical discussion of the VALMET valley
air pollution model. We begin with a short discussion of the coordinate
system and the means used in the model to incorporate topography. Discus-
sions of the mathematical equations used in both the nighttime and daytime
portions of the model then follow along with a statement of any assumptions
used in the model development. Unless otherwise noted, numerical values
of the parameters in the equations are assumed to be in the MKS (meter-
kilogram-second) system of units, except for pollutant concentrations
expressed in micrograms per cubic meter.
COORDINATE SYSTEM
VALMET is a two-dimensional, air pollution transport and diffusion
model in which pollutant concentrations are determined on a valley cross
section located an arbitrary distance down-valley from a pollutant source.
The distance x is measured following the axis of the valley down the center-
line of the valley floor (Figure 7). The local coordinate system utilized
on the valley cross section is shown in Figure 8. The x-coordinate system
specifies the distance of the cross section down-valley from the pollutant
source, as measured on the valley floor centerline following any turns in
the valley. The y-coordinate system has its origin at the valley floor
centerline, is orthogonal to the x-axis, and increases to the right when
viewed from the mouth of the valley. The z-coordinate system is a vertical
coordinate system with origin at the valley floor centerline.
TOPOGRAPHIC CROSS SECTION
For a particular valley cross section, the model assumes that the val-
ley topography can be simply represented by a horizontal valley floor of
width, w, and two inclined sidewalls of inclination angle a^ and o^, as
shown in Figure 9, below.
To simulate a real valley, the modeler must approximate the topography
of any cross section of interest using a topographic map, from which w, a t
and <*2 are estimated. Given the simplified valley cross section and the
local coordinate system, one may calculate several secondary topographic
parameters that will be used in the model. For example, using definitions
shown in Figure 10, the area of the cross section of a valley containing a
temperature inversion of depth h is
14
-------
Figure 7. Illustration of the nocturnal down-valley transport and
dispersion of a pollutant plume. The terrain following
coordinate system, with cross valley section an
arbitrary down-valley distance x, is illustrated.
15
-------
Figure 8. Illustration of the local coordinate system on a
cross valley section.
Figure 9. Parameters used in the VALMET model to approximate a valley
topographic cross section. The three parameters include
the valley width and the two sidewall inclination angles.
t
Figure 10. Illustration of the definitions of secondary topographic
variables used in the model.
16
-------
h2
A = h w +• j- (cot c^ + cot <*2) (1)
The y-distance to the left and right sidewalls from the vertical coordinate
axis is
yL(z) = -(£+ z cot a2J (2)
and
w
z) = j + Z cot a, . (3)
The topographic cross section may be considered to have a thickness of
unity. In this case the area AT of the top of the inversion (z = h) is
= yR(h) - y (h) = w + h(cot a + cot a ) . (4)
NOCTURNAL MODEL
Pollutant Source
The model begins with an elevated pollutant source on the valley floor
at the topographic cross section where x=0. The pollutant source may be
located anywhere on the valley floor. The y-position of the source is
specified relative to the coordinate origin. For example, the y-coordinate
of a pollutant source on the valley floor 200 m off the valley floor center-
line would be specified as yQ=+200 m if the source were in the positive y
direction from the centerline, or as y =-200 m if in the negative direc-
tion. The present version of the mode? cannot handle sources on the valley
sidewalls. The source is assumed to emit a continuous elevated pollutant
plume into a nocturnal valley temperature inversion of depth h=hi. This
pollutant plume is approximated in the model by a modified Gaussian plume
algorithm. The potential temperature gradient within the inversion is
T=86/3z, and the down-valley wind speed at plume carrying height is u$.
Wind speed is assumed not to vary with height. The model, being run for
worst-case dispersion under nocturnal temperature inversion conditions,
would be run for Pasquill-Gifford stability categories D, E or F. In later
sections, further details of the Gaussian plume formulation will be given
along with details concerning dispersion coefficents.
17
-------
Plume Rise
Plume rise is treated in the model according to formulas described by
Briggs [21-24] and used in the MPTER model [25].
The "2/3 law" of plume rise is used for all stability categories up to
a certain downwind distance Xf to calculate the plume rise using the
equation
1 6 F1/3 x 2/3
Ah = U6 F x (5)
where F is the plume's buoyancy flux, x is the downwind distance (m), and u,
is the wind speed at stack height (m/s). The parameter F is calculated as
-) (6)
_p
where g is the acceleration due to gravity (9.8 ms ), Ve is the exit
velocity of gases emitted from the stack (m/s), R is the stack inside radius
(m), and T and T are the exit and ambient temperatures (°K), respectively.
The ultimate rise Ah-r and the downwind distance x^ beyond which it applies
are calculated as shown below.
Unstable or Neutral Atmospheric Conditions—Calculations of Xf and Ahf
depend on whether the plume rise is "buoyancy-dominated" or "momentum-
dominated." This is determined by a crossover temperature difference,
calculated as
1/3 2/3
(AT). = 0.0297 T V_ ' /(2R) ' for F < 55
\* j C
or (AT) = 0.00575 T V 2/3/(2R)1/3 for F >_ 55.
v« j C
(7)
The actual temperature difference AT = T -T is compared to the crossover
temperature difference (AT)C and if AT _<_ (AT) then the plume rise is
momentum-dominated. If AT > (AT) then the plume is buoyancy-dominated.
18
-------
Momentum-Dominated Plume Rise—For a momentum-dominated plume rise the
distance to final rise is zero, such that
xf = 0. (8)
The final plume rise is given by
Ahf = 6RVe/us. (9)
Buoyancy-Dominated Plume Rise—For a buoyancy-dominated plume rise the
downwind distance of final plume rise is x.f - 3.5 x*, where x* is the
distance at which atmospheric turbulence begins to dominate entrainment,
given by
x* =
14 F when F < 55 m sec
9/c; A T (10)
34 F ' when F > 55 nT sec .
The final plume rise under these conditions is
1 C. C1/3 /O C *>
Ah 1.6 F (3.5 x*)
-
Mlx
(11)
Stable Atmospheric Conditions—Calculations of x^ and Ahf, as for unstable
conditions, requires a determination of whether the plume rise is buoyancy-
dominated or momentum-dominated. This is determined by a crossover
temperature difference, calculated as
(AT)C = 0.01958 T Vg s1/2 (12)
where s is a stability parameter given by
ae_
s = ~^- , (13)
19
-------
and ae/az is the vertical gradient of potential temperature within the
valley temperature inversion at sunrise. The actual temperature difference
AT = T$-T is compared to the crossover temperature difference and if AT <_
(AT)C the plume rise is momentum-dominated. If AT > (AT)C the plume rise" is
buoyancy-dominated.
Momentum-Dominated Plume Rise—For a momentum dominated plume rise the
distance to final rise is zero, such that
xf = 0.
(14)
The final plume rise is given by
=1.5
2 2
ve R T
VsJ
1/3
.-1/6
(15)
Equation (9) for unstable-neutral plume rise is also evaluated and the
smaller of the two plume rise calculations is used as the final plume rise.
Buoyancy-Dominated Plume Rise—When the plume rise is buoyancy-
dominated, the distance to final rise depends on wind speed at stack
height. The actual wind speed is compared to a certain critical wind speed,
uc, and the distance to final rise is calculated using one of two formulas
depending on whether the actual wind speed is greater than or less than the
critical speed.
The critical wind speed is calculated as follows:
uc = 0.2746 f1/4 s1/8
(16)
The distance to final rise is given by
,3/2
F-i/8,-9/1S
or
xf =
1.6.
Us s
for u >
20
-------
The final plume rise is given by
1.6 F1/3 xf2/3
Ahf = 1 . (18)
Gaussian Plume Model
Steady-state nocturnal pollutant concentrations within a valley cross
section are determined by means of a Gaussian plume model which is modified
to account for channeling and reflections off the sidewalls within the val-
ley, dilution due to clean air volume flux convergence within the valley,
and diffusion of pollutants out the top of the valley inversion.
We begin with the basic Gaussian plume equations having no reflection
terms:
XG(x,y,z;5;y0) = 2u a u
y z s
where
Gt(y) = exp [-^
(20)
G2(z) - exp {-(-}} (21)
and where XQ = (x,y,z;C;y0) is the pollutant concentration at a receptor
located at point P(x,y,z) due to a plume with center!ine height C and offset
distance from the valley floor centerline yQ. The factor 10 is the conver-
sion factor from kilograms to micrograms. a and a are the dispersion
coefficients, Q is the source strength and u is the transporting down-
valley wind speed at the stack. The pollutant stack is constrained to be on
the valley floor, so that
21
-------
Dispersion Coefficients
Nocturnal dispersion coefficients a and az for the Gaussian plume
equations are determined from the empirically derived Pasquill-Gifford
curves [26,27], as developed from atmospheric tracer experiments conducted
over homogeneous terrain. The empirical Pasquill-Gifford curves relate the
values of a and az to downwind distance and stability category. A numer-
ical approximation to the curves has been provided by McMullen [28] and is
given as:
a = exp [I + J (In x) + K (In x)2] (23)
where
a = standard deviation of the concentration in the horizontal (a )
or vertical (az) in meters, y
x = downwind distance (km), and
I,J,K = empirical constants for a given stability condition for a and
V
Tables 1 and 2 provide McMullen's [28] values of the constants I, J, and K
for determination of o and az> respectively. Table 3 is provided to deter-
mine the proper stability, category from meteorological observations.
The above method of determining a and az is appropriate for homoge-
neous terrain. Previous work [20,29-31] shows that dispersion is enhanced
in valleys over that experienced over homogeneous terrain where the original
empirical Pasquill-Gifford dispersion coefficients were determined. This
enhancement is particularly marked for stable conditions [20]. Enhanced
dispersion in valleys is due to a number of physical mechanisms that have,
to date, been inadequately studied or characterized. In the VALMET model we
have taken the approach, which we do not consider entirely satisfactory, of
calculating valley dispersion coefficients from the original empirical
Pasquill-Gifford curves [26] by simply adjusting the stability categories.
From existing data it appears that the enhancement of the horizontal disper-
sion can be approximated by dropping the stability two categories (e.g., F
becomes D). The vertical enhancement of dispersion is approximated by drop-
ping the stability one category (e.g., F becomes E). The use of these sim-
ple adjustments to the stability categories should be considered a stopgap
measure. Further work will be required to evaluate its appropriateness and
seek more appropriate means of handling nocturnal dispersion in deep valley
terrain.
Channeling
Calculations with the Gaussian plume model are made assuming the
Gaussian plume follows the valley axis, despite meanders in the valley's
course. The model user must therefore use a topographic map to determine
the down-valley distance of a given cross section.
22
-------
TABLE 1. VALUES OF THE CONSTANTS I, J, AND K, FOR a AS A
FUNCTION OF DOWNWIND DISTANCE, FOR SIX STABILITY
CONDITIONS [28]
Stability
Condition*
A
B
C
D
E
F
*As defined
I
5.357
5.058
4.651
4.230
3.922
3.533
in Table 3.
J
0.8828
0.9024
0.9181
0.9222
0.9222
0.9181
TABLE 2. VALUES OF THE CONSTANTS I, J, AND
FUNCTION OF DOWNWIND DISTANCE, FOR
CONDITIONS [28]
K
-0.0076
-0.0096
-0.0076
-0.0087
-0.0064
-0.0070
K, FOR CF_ AS A
SIX STABILITY
Stability
Condition*
A
B
C
D
E
F
I
6.035
4.694
4.110
3.414
3.057
2.621
J
2.1097
1.0629
0.9201
0.7371
0.6794
0.6564
K
0.2007
0.0136
-0.0020
-0.0316
-0.0450
-0.0540
*As defined in Table 3.
23
-------
TABLE 3. RELATIONSHIP BETWEEN WEATHER CONDITIONS AND STABILITY
CATEGORIES [26]
Day
Surface Wind
Speed
m
(at 10 m),
sec'1
<2
2-3
3-5
5-6
6
Incomi
Strong
A
A-B
B
C
C
ng Solar Radi
Moderate
A-B
B
B-C
C-D
D
ation
Slight
B
C
C
D
D
Night
Thinly Overcast
or >4/8
Low Cloud
E
D
D
D
<3/8
Cloud
F
£
D
D
The neutral class, D, should be assumed for overcast conditions during
day or night.
Plume Pilution
An important factor affecting plume concentrations is the dilution of
the plume during its travel down a valley due to tributary flows that bring
clean air into the plume. This process could be extremely important in
decreasing plume concentrations, especially when the volume of diluent air
is large relative to the plume volume, as when a pollutant's source is
within a small tributary to a larger valley. Other processes occur in val-
leys that can affect plume dilution, including converging downslope drainage
flows and entrainment of clean air into a valley from above. Any process
capable of producing volume flux divergence between two vertical valley
cross sections is capable of diluting a plume. At present, the physical
understanding of these processes is insufficient to incorporate them
explicitly into an air pollution model. The extent to which volume flux
divergence is a feature of a particular valley's meteorology can be deter-
mined, however, by wind observations.
The approach used in the present version of the model is to modify the
Gaussian plume equation in a simple way that makes an initial correction for
the dilution of a pollutant plume during its down-valley travel. The
correction is applied only when data are available to show that a plume
would actually be diluted in the valley being modeled. The correction is
applied to the plume transporting wind speed at the pollution source cross
section to obtain a modified or virtual wind speed uv- The virtual wind
speed is a function of the ratio of the volume fluxes across the source and
receptor cross sections, such that
24
-------
ucAc
usAs
(24)
where uc is the measured plume-transporting wind speed at the receptor cross
section and
source cross
Ac and AS are
sections,
the cross sectional areas at the receptor and
._ , respectively. The value of uv calculated using this
equation is simply substituted for u? in the denominator of the Gaussian
plume equation to calculate pollutant concentrations at the receptor cross
section. To avoid physically unrealistic results (i.e., plume concentra-
tions increasing with transport distance) uv must always be equal to or
greater than u . This is accomplished by specifying
ucAc
usAs
> 1
(25)
Using this plume dilution approach, pollutant concentrations at the receptor
cross section are calculated using the equation
y =
xv
126)
The plume dilution method described above is an initial approach that will
need improvement as the physics of the dilution process are more fully
understood. The advantages of the method are that it should produce more
realistic concentrations than when plume dilution is not considered, and the
method can be supported by simple measurements. The above method suffers
from several drawbacks, however. First, the separate effects of individual
physical processes that can lead to dilution are not clearly separated using
the current method. Second, the dilution method does not consider the
effects on plume dilution of the enhanced lateral and vertical mixing where
tributary flows merge into a main valley or of the drift of a olume center-
line towards one sidewall caused by incoming tributary flows. Third, the
method does not consider stratification effects that may arise in merging
flows of different temperature. Finally, dilution air is assumed in the
model to be clean with respect to the pollutant species.
method seems appropriate, however, as an initial means of
dilution until a better understanding is obtained for the
tion (and dispersion) processes. The wind observations necessary to use the
dilution method could be obtained from a doppler acoustic sounder or special
pibal and/or tethersonde campaigns. A user should use data to show that the
plume would be diluted due to volume flux divergence.
The dilution
handling plume
individual dilu-
25
-------
Reflections
Dispersion of pollutants during the plume's nocturnal transport down
the valley will result in a spreading of the plume in both the vertical and
horizontal, resulting in the plume's eventual contact with the valley sur-
faces (valley floor and sidewalls). Additional spreading of the plume,
using the basic Gaussian equations, would give the physically unrealistic
result of appreciable concentrations occurring "beyond" the valley sidewalls
or "below" the valley floor. Concentrations calculated for receptors within
the valley cross section would be too low and the total pollutant mass of
the plume would not be conserved within the valley. Several methods have
been proposed to remedy this situation [26,32]. The valley may, for
example, be treated as a channel having vertical sidewalls. Then, single or
multiple reflections of the plume from the sidewalls can be simulated by the
method of "image plumes", in exact analogy to the way that reflections off
an elevated inversion surface are now handled in many Gaussian models
[26]. Other similar means of handling such plume reflections have been
tried by the authors giving consideration to the inclined sidewalls of real
valleys, but all such methods thus far investigated have suffered from draw-
backs of one type or another. The most satisfactory solution came from a
consideration of pollutant mass conservation within the plume, to be
discussed in the following paragraphs.
The Gaussian plume equations are a solution to the classical diffusion
equation obtained by means of a number of restrictive assumptions (steady-
state motion, au/9z=0, etc.) and application of certain boundary condi-
tions. Among the boundary conditions is an expression of plume mass conser-
vation, given by
Q = / / u x dy dz for z > o . (27)
This integral constraint, which applies to the basic mathematical form of
the Gaussian equation [Equation (9)], has a simple physical meaning.
Suppose that a source emits a pollutant at the rate Q of 1 kg/s into a
horizontal wind of 1 m/s. This will result in a kilogram of pollutant mass
being carried in a meter of transported air. Following the Gaussian plume
equation this kilogram of pollutant is simply distributed in the 1-m slice
of air following a bi-normal distribution. The total plume mass within the
1-m slice, however, does not change. It is 1 kg regardless of whether the
kilogram is relatively concentrated into a tightly confined plume or widely
dispersed about the plume center. Thus, when one integrates the pollutant
concentration over the volume of the slice from y=-o> to y=+°= and from z=-<=°
to z=+°°, and uv is constant, we obtain
26
-------
£-- / / X.. dydz=i^. (28)
"y -co _co
In terms of a worst-case model of pollutant dispersal in a valley, one must
consider ways of preserving or conserving the total pollutant mass so that
this relation still holds.
The total amount of pollutant mass in a unit cross section perpendicu-
lar to the plume path is determined by Q and uv, but the distribution of
pollutant mass about the plume centerline is determined by the bi-normal
probability density function. The characteristics of this function are such
that
/ / ~T^exp U --— exp {- % j dYdH (29)
a
z
where
V - and „ .
A valley cross section can be superimposed on the theoretical plume cross
section (Figure 11) and, through integrations, one can determine the frac-
tion of total plume mass within the cross section within the valley and
below the inversion top, as well as the mass that has diffused above the
inversion top and the mass that has diffused beyond the valley sidewalls and
valley floor. For example, the fraction of the total mass in the cross sec-
tion that is within the valley cross section below the inversion top hi is
^ /R 1 / A 1 / H2\
f, = / / ^r"— exp - 7 J • ^^- exp (- f J dYdH (30)
1 o Y, z - '
27
-------
GAUSSIAN PLUME
TEMPERATURE INVERSION TOP
^>;^;lj POLLUTION DIFFUSING OUT THE INVERSION TOP
|| | | || POLLUTION REFLECTED BACK INTO VALLEY CROSS SECTION
Figure 11. Cross section of pollutant plume and valley topography
illustrating the integral method of calculation for
plume reflection and diffusion out the top of the
valley temperature inversion.
where
h '
y " y
y " y
The fraction of plume mass that has diffused out the top of the inversion is
fo =
Hi -
—— exp
exp
/2-n: a
- / dYdH
(31)
The fraction of plume mass which has diffused "beyond" the valley sidewalls
(below the inversion top) and "below" the valley floor is
= 1-frf2 •
(32)
28
-------
In our simple model of plume reflection we must find a means of folding the
mass fraction f3 back into the valley cross section so that pollutant mass
will be conserved within the valley. Although the pollutant mass could be
distributed into the valley cross section any number of ways, we choose the
simple expedient of mixing it uniformly within the valley unit cross section
volume, so that concentrations within the cross section include an offset
concentration and are given by
x(x,y,z;r;yo) = fXy + (33)
where the factor PQT/pT adjusts the concentrations to standard atmospheric
conditions, with pQ = 1013 mb, T = 293.16°K and T and P are the ambient
temperature (°K) and pressure (mB). The offset concentration is given as
fol 1 ows:
v - f
yoff T3
Equation (33) applies within the spatial domain where
i) x > o,
ii) yLfz) £ y _£ yR(z), and
iii) o <_ z <_ h.j.
By mixing the lost mass uniformly we are assured that, in the limit of long
travel distances x, plume concentrations within the valley cross section
will approach uniformity.
The above method provides a simplified means of allowing some plume mass to
escape upwards out of a valley during its down-valley transport, as well as
providing for "reflections" of the plume from sidewalls and the valley
floor. An objection to the volumetric mixing of the reflected plume mass
through the entire cross-sectional volume seems necessary, however. This
mixing of reflected plume mass through the whole valley cross section will
result in an underestimate of concentrations near the sidewalls and valley
floor. The amount of underestimation is not easily calculated due to a lack
of understanding of the physics of the reflection processes at the bounda-
ries during nighttime conditions. The effects of deposition on valley sur-
faces may reduce ambient air concentrations near the valley surfaces,
however, and will tend to counteract these errors. Further work is deemed
necessary to evaluate both the effects of deposition and the presence of
concentration underestimates on the valley surfaces.
29
-------
The integrations in Equations (30) and (31) are accomplished numeri-
cally in VALMET using a polynomial equation given by Abramowitz and Stegun
[33] attributed to Hastings [34], The Gaussian or Normal probability
density function is given by
Z(x) =--exp (- f-) . (35)
where variable x is the standard deviation.
The area under the Gaussian function from -°° to x is given by
x 2 x
P(x) = — / exp (- f-] dt = / Z(t) dt . (36)
/2lt -°° -05
Using these definitions, we note that Z(-x) = Z(x) and P(-x) = l-P(x).
Following Hastings [34] we may approximate P(x), for 0 <_ x _<_ °°, by the
polynomial equation
P(x) = 1 - Z(x) [bxt + b2t2 + b3t3 + b4t4 + b5t5] + e(x) (37)
where
1
t =
1 + 0.2316419 x '
-8
|e(x)| < 7.5 x 10
and
b-, = 0.319381530
bi = -0.356563782
b3 = 1.781477937
b4 = -1.821255978
b5 = 1.330274429.
30
-------
The method for the numerical integration is illustrated below for
Equation (30).
(38)
exp L- 7
Az
where
n = h/Az
and z.j = iAz.
The numerical integration for f^ is accomplished in an analogous manner from
Equation (31). The approximation is made in VALMET using +5 standard devia-
tions as the limits of integration for y as this makes the calculations
tractable and produces no appreciable error in the results.
Calculation of Steady-State Nighttime Concentrations
Calculation of nighttime pollutant concentrations at receptor point
P(x,y,z) can be made in a straightforward way with Equation (27) assuming
steady-state conditions. The receptor must be located within a valley cross
section, so that Equation (33) can be applied for
i)
11
any
0 <
<_ h.j, and
iii) yLTz) <_y <_yR(z)
In practice nighttime calculations are performed for a limited number of
points on the valley surfaces (floor and sidewalls), choosing points that
correspond to locations where daytime pollutant concentrations will be
calculated. The fixed daytime grid configuration is shown below in
Figure 12. Since the model is run on a half-valley cross section, the cross
section sidewall angle a is chosen to be the average of the observed side-
wall angles a, and ou- Points for which nighttime pollutant concentrations
are calculated are snown on the figure with x's.
If the grid element numbers are given by n = 1,2, ...NBOX, locations at
points P(y,z) on the valley floor are at
y = BOX,LEN + (n - 1) BOXLEN and z = 0 .
31
-------
w/2
Figure 12. VALMET grid configuration on a valley cross section illustrating
the nomenclature used in the model. The cross section is
representative of the model configuration at sunrise, when the
initial CBL-height (l-h), inversion top-height (h^
center-line height (?) are as shown. The x's on th
illustrate the locations where nocturnal pollutant
concentrations are determined.
and plume
e figure
Calculation points on the valley sidewalls (n = NBOX + 1, NBOX + 2,
NTOTI) are at points P(y,z) where
y = 2" + z cot a ,
BOXLEN tan a
z = - + (n - NBOX - 1) BOXLEN tan a ,
and a =
DAYTIME MODEL
The daytime portion of VALMET begins with solar radiation calculations,
a parameterization of the surface energy budget, and estimation of the time
variation of sensible heat flux which destroys the temperature inversion
after sunrise. The valley energy budget model of Whiteman and McKee [3] is
invoked to simulate CBL growth and inversion descent. The nighttime
32
-------
Gaussian plume on a valley cross section is used as the initial condition
from which the numerical model proceeds to calculate concentrations in grid
elements fixed to the valley floor and sidewalls. The depth of these verti-
cally growing grid elements is identified with the CBL height. Calculations
are made of concentrations within the grid elements taking account of fumi-
gation of the Gaussian plume into the elements and upslope advection from
grid element to grid element. Pollutants are mixed uniformly within each
grid element, roughly simulating atmospheric diffusion processes in the con-
vective (actually, advective) boundary layer.
Extraterrestrial Solar Radiation
Solar radiation on a plane surface above the earth's atmosphere can be
calculated using basic principles of spherical trigonometry. The basic
formula for solar flux on a horizontal surface is
Q=S(j cosZ (39)
_
where SQ is the solar constant (1367 Wm ), defined as the solar flux on a
plane surface perpendicular to the sun's rays at the mean earth-sun dis-
tance, d. The factor (d/d ) modifies the solar constant to account for the
fact that the earth-sun distance, d, varies during the year as the earth
travels in an elliptical orbit about the sun. The cosine term accounts for
the projection of solar flux onto the horizontal plane. This term is
clearly a function of earth-sun geometry, depending on the time of day, day
of year, and latitude of the surface. Z is the angle between the normal to
the surface and the direction to the sun, known as the zenith angle. From
spherical trigonometry [35,36] this is given as
cos Z = sin sin & + cos = latitude
6 = sun's declination
X. = sun's hour angle (i.e., an angular measure of time reckoned from
solar noon where 15° represents one hour). For example, 30°
represents 2 hours after solar noon and -30° represents 2 hours
before solar noon.
— 2
In these equations both the declination and the factor (d/d ) are a function
of the day of year (D, from 1 to 365). McCullough [37] gives some approxi-
mate analytical formulas that allow the calculation of these terms, utiliz-
ing the longitude of the earth (A) in its orbit around the sun as reckoned
from the earth-sun radius vector at vernal equinox (D = DQ). The formulas
are
33
-------
<5 = sin" (sin £ sin \) ,
= (D - DQ) + 2e(sin cuD - sinuDQ),
and
-T 2
(±)
- e cos
(41)
(42)
(43)
where
E = 23°26' = maximum solar declination,
01 = 2ir/365,
DQ = 80 = day of vernal equinox, and
e = 0.0167 = eccentricity of earth's orbit.
Using these equations, the extraterrestrial solar flux on an arbitrarily
oriented plane surface can be calculated as a function of hour angle, n.
Input to the equations must include the latitude and the day of year for
which calculations are desired.
The equations are used to calculate the solar flux at any time from the
hour angle of sunrise, £_R, to the hour angle of sunset, £<-(-. For a
horizontal plane, the times of sunrise and sunset are generally specified in
terms of half-day length, R, calculated from Equation (40) for cos Z = 0, or
Then
R= cos"1 (-tan tan 6)
= -R and
(44)
(45)
Calculations using the above set of equations are made in terms of hour
angle and are thus referenced to the time of solar noon at a particular
site. To convert the hour angles to local standard time it is necessary to
determine the local standard time of solar noon, and to adjust all the times
accordingly. This is accomplished by making two corrections to the time of
solar noon (12h 00m OOs) in the local solar time coordinate system [38].
Then the standard time of solar noon is given by
tNOON = 12h 00m OOs + Equation of + Longitude
time correction correction
= 12h 00m OOs
ET
(46)
First, the equation of time correction is applied by adding the appropriate
number of minutes and seconds from Table 4. Second, a longitude correction
accounting for the difference in geographical longitude between the site and
the reference meridian of its time zone is applied. The correction is
+4 minutes for each degree of longitude west of the meridian. For example,
34
-------
TABLE 4. EQUATION OF TIME CORRECTION
0* CFT
1 3m
16 9
31 13
46 14
61 12
76 8
**
12s
32
24
16
23
39
D
91
106
121
136
151
166
C
4m
-0
-2
-3
-2
0
ET
08s
01
51
44
33
10
0
181
196
211
226
241
256
C
3m
5
6
4
1
-3
ET
21s
46
21
44
07
50
D
271
286
301
316
331
346
C
-9m
-13
-16
-15
-12
-6
ET
06s
33
06
53
34
34
D CET
361 Om 47s
* D = day of year
** Cry = equation of time correction (minutes and seconds).
the reference meridian of the Mountain Time Zone is 105°W longitude. For a
site at 106°W longitude, the correction is CL = 4(106-105) = +4 minutes.
The above equations can be used to determine extraterrestrial solar
flux at any site on a given day for any given time (i.e., hour angle).
Calculations using these equations for a latitude of 40°N and longitude of
109°W show that the solar flux QSb as a function of time on a given day
closely follows a sine curve. This simple analytical form proves so accu-
rate at the latitudes and longitudes of interest and is so easy to calculate
that it is used in the model instead of the full formulation. The sine
approximation is given by
Qsh = Ax sin ~ (t - tSR) (47)
where
Ai = solar flux at solar noon
T = length of daylight period
tSR = time °^ Sunnse» and
t = time.
One can calculate A-^ by evaluating Equation (40) with 2 = 0, such that
_ 2
Al = So tfj (sin sin S + cos cos <5) .
The length of the daylight period is given by
T = 2R = 2 cos" (-tan £ tan 6) . (49)
35
-------
The time of sunrise is evaluated in local standard time using
If the model equations are used far from the latitudes of the central Rocky
Mountains, the model user should verify that the sine approximation is a
valid representation of solar flux, by comparing its calculations with the
full analytical equations.
Parameterization of Sensible Heat Flux
The sensible heat flux which constitutes the basic driving term of the
daytime portion of the model is parameterized from the extraterrestrial
solar flux using the formula
= A0[Al
sin
I (t - tSR)] » AQ Qsh (51)
where A0 is a fraction between 0 and 1. This formula should be recognized
as a crude parameterization for sensible heat flux, simply stating that the
sensible heat flux is a constant fraction of the time-dependent extrater-
restrial solar flux, as illustrated below in Figure 13.
AQ depends on a great many factors which are not yet included
explicitly in the model. In general, A should be a time-varying function
dependent on the
transmissi vity of the earth's atmosphere,
cloud cover,
surface albedo,
soil moisture,
surface cover or vegetation,
longwave radiation budget, and
other factors.
A future version of VALMET should include some of these factors explicitly, at
least the ones which may affect worst-case pollutant concentration, using as a
basis the extensive work already conducted on these topics by other
investigators.
Model Grid
The daytime portion of VALMET uses a numerical technique to simulate time
varying pollutant concentrations in grid elements that are fixed to the valley
floor and sidewalls. Calculations are made for individual time steps during
the temperature inversion breakup period beginning at sunrise and using the
nocturnal steady-state concentrations within the valley cross section as an
initial condition.
36
-------
HEAT FLUX
Q
Sh
A0Ai _
Figure 13. Extraterrestrial solar heat flux as a function of time, showing
solar model nomenclature. The figure illustrates how the
sensible heat flux F is parameterized as a fraction AQ of the
extraterrestrial solar heat flux on a horizontal surface above
the site. The times of sunrise and sunset
daylength T is defined in the model as the
between sunrise and sunset.
are shown. The
time difference
The geometry of the valley cross section, illustrating the numerical
grid, is shown in Figure 14. Calculations in the model are performed on a
half-valley cross section.
Grid elements are arranged from the valley center across the valley floor
and up the sidewall, with grid element height representing the height of the
CBL, assumed not to vary from grid element to grid element. At sunrise the
model is initiated with a 25 m CBL height, and an inversion top height that is
obtained from observations in the modeled valley. These initial heights
change after sunrise as sensible heat flux drives the destruction of the
valley inversion. Specifically the height of the CBL, H(t), increases and the
height of the inversion top, h(t), decreases. The simulation is accomplished
by treating the Gaussian plume within the cross section as being "frozen"
within the inversion. After sunrise, concentrations within individual grid
elements change as the grid elements grow upwards into the frozen plume, as
the inversion sinks causing pollutants to sink into the tops of the grid ele-
ments, and as upslope flow develops within the CBL. These upslope flows
develop as air parcels are heated by sensible heat flux over the sidewalls and
rise up the slope. The speed of the slope flows is calculated under a con-
tinuity of mass constraint within the cross section below the (sinking) inver-
sion top. That is, sinking of the top of the inversion with time implies that
mass is removed from the cross section below the inversion top level. In the
37
-------
Figure 14. Cross section of the valley floor and sidewalls illustrating the
grid elements whose height corresponds to the C8L heights. The
rate of sinking of the top of the inversion (arrow) is related
to the upslope transport of mass (other arrow) by the equation
of mass continuity.
model, the assumption is made that the upslope flows carry the mass required
to allow the inversion to sink at the rate predicted from Whiteman and McKee's
[3] bulk thermodynamic model. This model is driven by sensible heat flux
considering the thermodynamic energy budget of the valley cross section.
One of the assumptions in the model is that the top of the valley inver-
sion at sunrise extends horizontally across the valley cross section. This
assumption is supported by data from an observational program conducted in
Colorado mountain valleys [9], In that study several experiments were con-
ducted by dual tethersondes located on a cross valley section. These obser-
vations showed that the inversion top was nearly horizontal. This is not
unexpected, since hydrostatic forces should produce settling of cold air
masses in topographic depressions or valleys and would result in horizontal
isentropes in the absence of major cross valley inflows. This effect would be
more pronounced for strong inversions where the hydrostatic forces would be
stronger. Gravity waves are known to propagate frequently on the upper
boundary of such inversion layers and are frequently noted in acoustic sounder
records. These waves are generally of small amplitude and should not greatly
affect the validity of the assumption of horizontal homgeneity.
Thermodynamic Equations for CBL Ascent and Inversion Descent
For details of Whiteman and McKee's thermodynamic model of temperature
inversion breakup in mountain valleys the reader Is referred to the original
article (Appendix C). Here a brief explanation is given.
38
-------
One can view a valley temperature inversion as a pool of cold air trapped
within a valley. This pool has an energy deficit relative to the air above
the valley, and the temperature inversion can be destroyed only by warming the
pool to the same potential temperature as the air above the valley. The
amount of energy required to warm the pool can be calculated from the First
Law of Thermodynamics knowing the depth of the pool, the vertical temperature
gradient within the pool (assumed constant) and the cross valley dimensions of
the pool, by means of a volumetric integration through the valley cross sec-
tion. This energy deficit is overcome by sensible heat flux, which can be
estimated as a fraction of the downward solar radiation coming across the area
of the top of the valley temperature inversion. Conversion from solar energy
to sensible heat occurs at the ground. Downwards solar heat flux across the
inversion top will be intercepted on a valley surface below the inversion top
and will be available as sensible heat flux to warm the inversion airmass.
The topography of the valley strongly affects the course of the inversion
destruction by determining which of the patterns of inversion destruction will
be followed and by controlling the overall rate of heating of the valley atmo-
sphere. This control on the rate of heating can be explained by Figure 15.
Consider, for example, a valley that is completely filled by a temperature
inversion. Solar energy coming downwards across the top of the valley
inversion will be converted to sensible heat flux at valley surfaces and will
be used to heat the inversion volume. The ratio of the area at the inversion
top through which this radiation comes to the inversion volume governs the
rate of warming of the valley atmosphere and, consequently, the time required
to destroy the inversion. Given an equally strong and deep inversion over the
plains, the valley inversion will be destroyed more rapidly since the same
energy input goes to heat a smaller volume of air in the mountain valley. The
thermodynamic model includes this important factor, which arises from the
volumetric integration in the derivation of the model equations.
I J_J LI
L_L_L_L_LLJ
1
VALLEY
PLAIN
Figure 15. Illustration of the effect of topography in controlling the
heating rates of the air within a valley temperature
inversion versus the air within an inversion over the
plains. The same incoming energy heats a smaller volume of
air in the valley case.
Whiteman and McKee's model is a "bulk" thermodynamic model in that it
does not differentiate between sensible heat flux over different valley sur-
faces, dealing only with a bulk heat flux. The heat flux which drives the
valley inversion destruction is partitioned or distributed in a fundamentally
different way from that for an inversion over homogeneous terrain. There the
39
-------
sensible heat flux destroys the inversion by driving an upward growth from the
ground of a convective boundary layer, which warms the inversion air mass from
below until the temperature deficit is overcome and the inversion is
destroyed. In contrast, in a valley the upward heat flux over valley surfaces
develops a convective boundary layer but also, due to sensible heat flux
convergence over the slopes, causes warmed air parcels to flow up the slopes
in an upslope flow that develops within the convective boundary layer. These
upslope flows remove mass from the base of the temperature inversion in the
shallow slope flows, and through mass continuity, results in a general subsid-
ing motion over the valley center. The atmospheric energy budget approach
used by Whiteman and McKee is capable of partitioning energy between these two
different processes to produce inversion destruction solely by CBL growth,
solely by inversion descent (assuming a non-growing CBL is present initially
in a simulation), or by a combination of the two processes. The partitioning
is controlled by a single parameter, k, defined as the fraction of sensible
heat flux going to CBL growth. The remaining fraction, 1-k, is assumed to be
responsible for mass transport up the CBL which results in inversion
descent. The proper partitioning of energy (i.e., the most appropriate value
of k) in a particular valley can be determined by fitting the thermodynamic
model to observations in a particular valley. This was done for a Colorado
valley by Whiteman and McKee [9], although they cautioned that the value of k
will depend on other factors besides valley geometry, including the magnitude
of sensible heat flux.
Further research on the functional value of k is suggested as a high
priority in this modeling approach, although the model is functional using the
curve fitting approach to determining k.
The thermodynamic model is composed of two coupled equations. The first
equation is a prediction equation for CBL height where, in accordance with the
bulk nature of the model, the CBL depth H is assumed not to differ over the
valley floor and sidewalls. The first equation is
= HJ + AHJ (52)
where j is a time step index and
n , w + H .C A A,
f
Note that a non-zero initial CBL height is necessary to make this prognostic
numerical equation tractable. In the model, this requirement is met by using
an initial CBL height at sunrise of 25 m. The second equation, describing the
inversion top height h is
40
-------
(54)
where
h.+n .
h. g
h,Y(w + f- C) -j (t-t.)(w+h.C)
J £• ^- ' J
At' (55)
T"
1000, .286
p = atmospheric pressure,
c = specific heat of air at constant pressure,
C = cot ctj + cot a2,
Ho E Hi * °»
ho = hi'
t = (j+1) At + t-j, and
j = 0,1,2,3,....n.
The terms in the numerator and denominator of Equation.. (54) involving the
warming coefficient for the air above the inversion, 6(°Ks~ ), allow the model
to incorporate the retarding effect on temperature inversion breakup caused by
warm air advection above the valley temperature inversion. Extra energy is
required to destroy the valley temperature inversion if this warming occurs
during the temperature inversion breakup period since the inversion cannot be
broken until the entire valley atmosphere is warmed to the temperature of the
air above the valley.
The numerical simulation using these coupled equations proceeds with
discrete time steps and is completed when the inversion is destroyed at the
first time step n at which the CBL height becomes greater than the inversion
top height, such that
(56)
In the limit of an infinitely wide valley (i.e., a plain) and k=l (i.e., all
sensible heat flux going into CBL growth) the equations reduce to a single
equation describing CBL growth over homogeneous terrain. In the limit as k
goes to 0, the equations describe inversion destruction where the CBL fails to
41
-------
grow after attaining an initial height, and inversion destruction occurs due
to sinking of the inversion and adiabatic warming of the valley atmosphere.
Whiteman and McKee [3] have shown that this approximates inversion destruction
in snow covered valleys. Values of k between 0 and 1 allow the coupled
equations to describe intermediate methods of inversion breakup in which the
sensible heat flux is partitioned between CBL growth and inversion sinking. A
typical inversion destruction in the fall in the Eagle Valley of Colorado was
simulated well using a value of k=0.14. The appropriate value of k for other
valleys would be selected based on careful analysis of valley atmospheric
data, as discussed further in Section 4.
Sinking of the inversion layer will result in a sinking of the frozen
plume within the inversion. The pollution mass can be calculated by consider-
ing the vertical Gaussian pollution profile above each of the grid elements.
As shown by Equations (52) and (53) for k*0 the CBL depth (and hence grid
element height) increases with each time step. This CBL growth entrains
pollutant mass from the inversion's stable core into the top of each grid
element. These two processes, sinking of the frozen plume into the top of
each grid element and growth of the grid element upwards into the frozen
plume, combine to determine the amount of pollutant mass introduced into each
grid element during each time step.
As illustrated in Figure 16 for a grid element on the valley floor, these
two processes result in a volume of air of depth AH. + Ah- and length y, - y~
being entrained into the growing grid element during time step j. The average
pollutant concentration in this air volume can be calculated from the Gaussian
plume equation [Equation (33)] by averaging the pollutant concentrations at
the centerpoints of the four sides of this volumetric element (Figure 17). A
z-coordinate transformation is useful in these calculations so that one may
think of the coordinate origin on the valley floor rising into the fixed
plume, rather than the frozen plume sinking towards the valley floor.
Advection in the Slope Flows
Conservation of total air mass on the two-dimensional valley cross sec-
tion of interest implies that there must be a relationship between the sinking
of the inversion top and removal of mass from below the inversion top due to
upslope transport of mass in the convective boundary layer over the slope.
Thus, at the level of the inversion top, a slow subsidence of the broad stable
core must be balanced by stronger upward motions through this level due to
flow up the narrower sidewalls CBLs.
The continuity equation for air of constant density on the two-dimen-
sional cross section may be written as
_ _ . }
3y ' • 3l • (57)
42
-------
Ah;
I
'AH?
CBL
J+1
//////////^^^^^
Figure 16. Diagram showing the changes in CBL and inversion depth
at a given time step above a valley floor grid
element. The amount of air volume incorporated into
the growing CBL at each time step is seen to be repre-
sented by (y2 - y-^'Uh. + AH.) cubic meters.
Ah:
AH:
v v
M '2
Figure 17. Representation of the volumetric element of mass
incorporated into the growing CBL above a valley
floor grid element at each model time step. The
x's in the figure represent the positions where
pollutant calculations are performed at each time
step to estimate the average pollutant concentra-
tion carried within the volumetric element during
the time step.
43
-------
This equation may be used to calculate wind velocities at the boundaries or
septa of the grid elements in the slope flow layer. For the first grid
element on the valley floor, we assume that the wind velocity on the left
septum at the valley center is zero (i.e., there is no mass transferred across
the valley center within the CBL). At a given time step a certain amount of
mass sinks into the grid element of length L and depth D. Integrating,
, ,
/ §7 dy dz = - / / |f dy dz
\ * r\ r\
(58)
or
L -
(59)
Thus, an inversion sinking at the rate of .01 m/s into the top of a 100-m-long
grid element of depth 50 m will produce a wind speed increase from the left
septum of the grid element to the right septum of the grid element of
VL - Vo =
100
.01 = .02 m/s .
(60)
The velocity, v, normal to a given septum of any grid element can thus be
determined by summing the velocity changes calculated across each grid element
beginning at the valley center and ending at the septum of interest.
Pollution Concentration Calculation Method
Pollution concentration calculations are made for each time step for each
grid element following the general outlines of a numerical method described by
Whiteman and McKee [11]. Following this method, calculations at each time
step are performed sequentially for all of the grid elements, beginning with
the grid element at the valley center and progressing up the valley sidewall.
Calculations within each grid element (Figure 18) are made under a pollutant
mass conservation assumption. The pollutant mass sinking into the top of each
growing grid element is calculated using the method described above. Simi-
larly, the advection of pollutant mass into and out of each grid element is
accomplished using the velocity calculations at grid element septa as deter-
mined using a total air mass conservation constraint as described above. The
concentration Cn t of pollutant in grid element n at the end of the t-th time
step of length At is given by a simple pollutant mass balance:
44
-------
»»M5
M3
Figure 18. Schematic diagram of an individual model grid element
illustrating the pollutant mass balance.
M. .„. . + M. - M ,.
_ initial in out
'n,t
(61)
M1n,t + M2n,t + M3n,t + M4n,t " M5n,t
where
Vn is the volume of the nth grid element
Ml is the mass of pollutant coming into the growing grid element n during
time step At
M2n is the mass of pollutant advected into box n from the adjacent upstream
grid element (n-1) during time step At
M3n is the net mass of pollutant introduced into grid element n during time
step At from pollutant sources and sinks within the grid element
45
-------
M4n is the initial mass of pollutant in grid element n at the beginning of
time step At
M5n is the mass of pollutant advected out of grid element n into the adjacent
downstream grid element (n+1) at the end of time step At after complete
mixing of masses Ml through M4 has taken place within the grid element.
The separate terms are calculated as follows:
Ml
n,t
AHJ
(62)
where xn is the average pollutant concentration within the volume element
(Ah, + AH.) that is incorporated into the grid element n at time step t (or,
following our previous nomenclature, j). Ai , the area of the grid element t
(Figure 18), should not be confused with the A, used previously in
Equation (40).
top
n .. = C .
n,t n-l,
(63)
where
A2 = area of box septum
vi = wind velocity component into grid element normal to box septum.
M3
n,t
- S) AlAt
(64)
where
O 1
Q = pollutant emission rate [ML T ]
_p I
S = pollutant removal rate [ML" T ].
M4n,t = Cn,t-l Vn
(65)
M5
n,t
Ml + M2 + M3 + MA
n,t n .t n ,t n.t
V
n
(66)
46
-------
where
V£ = wind velocity component out of box normal to box septum.
Exponential Decay of Concentration
Calculation of concentrations in individual grid elements follows the
procedure outlined. A modification of this procedure is necessary, however,
for individual grid elements on the sidewalls when the top of the inversion
sinks below the level of the grid element, as the depth of the CBL suddenly
increases when this occurs. Concentrations would then be expected to decrease
rapidly. This is handled in the model by simply allowing the last concentra-
tion calculated in the grid element to decay exponentially with time after the
inversion top descends below the grid element. In view of the lack of
observational evidence, the time constant of the exponential decay is chosen
arbitrarily so that the concentration decreases by 2% at each 10-s time step,
such that
r -go At/10 r (67]
cn,t >yB Ln,t-l ' lb/;
As with any numerical model, one must be cognizant of the possibility of
numerical instabilities affecting model results. The criteria for maintaining
stability [39,40] constrains the model user to maintain a certain relationship
between the model time step and the grid element length, depending on the
upslope wind speed encountered in the simulation. This relationship is
vmax At < y2 - YI . (68)
Thus, the maximum upslope wind speed simulated in the model, when multiplied
by the time step length must be less than the grid element length in order to
maintain computational stability. In the VALMET model, the user may choose
(yo-y-i) but has no control over At and v . Calculations are made auto-
matically within VALMET to ensure computational stability. Since the maximum
upslope wind velocities should never exceed 10 m/s, we may ensure computa-
tional stability by setting
y9 - y,
At = *1Q X . (69)
Thus, for example, a time step of 10 s is sufficiently small to ensure com-
putational stability in a model simulation with 100 m grid elements.
Following the numerical method outlined above, concentration may be
calculated for every time step for every grid element. Rather than storing
each of these concentrations, which would take a great deal of computer
47
-------
memory, 5-min-average concentrations are stored. This greatly reduces
computer storage requirements. The 5-min-average concentrations constitute
one of the main outputs of the model.
Maximum 1- and 3-Hour-Average Concentrations
Another ouptut of the model is the maximum 1- and 3-hr-average concentra-
tions and their times of occurrence in each of the grid elements. These
calculations are made from the basic 5-min-average concentration array for
each grid element by means of a moving average method. Initial tests have
been conducted to confirm that the maximum 1- and 3-hr-average concentrtions
calculated from the 5-min-average array do not differ significantly from
concentrations calculated from a 10-s-average concentration array.
48
-------
SECTION 4
OVERVIEW OF MODULAR VALMET MODEL
The initial valley air quality model called VALMET has been designed and
constructed based on hypotheses presented in Section 2 that arose from an
observational study of Colorado valley meteorology. The model uses the tech-
nical approach and equations outlined in Section 3.
The VALMET model can be used to simulate the transport and diffusion of
pollutants released from an elevated source in a well-defined mountain valley
during the nighttime and morning transition periods. The model operates on a
valley cross section an arbitrary distance down-valley from an air pollution
source and has been constructed to include parameterizations of the major
physical processes that act to disperse pollution during these time periods.
Before a modeler attempts to use VALMET to simulate dispersion in a
particular valley, he should critically review Sections 1 through 3 of this
report. Since the modeling approach is phenomenological, and individual
physical processes affecting pollutant dispersion are parameterized, the
modeler should be wary of applying the model to a meteorological or topo-
graphical situation where model assumptions are invalidated or physical pro-
cesses parameterized in the model are clearly not occurring. The modeler
should carefully review existing meteorological data for the valley of
interest to see whether the model applies.
Even if the model's assumptions seem valid, the present state of the art
in such modeling results in the use of rather arbitrarily specified parameters
within the model. Further work is necessary to refine our estimates of the
values of these parameters (e.g., k, a and a ). This work would benefit
from comparisons of model results witlr actual diffusion trials in real val-
leys. In view of our rather tentative understanding of valley meteorology,
the model has been constructed in a modular fashion. This modularization of
the code is a major design feature of the model and should allow the code to
be modified easily as we learn more about the meteorology of valleys and find
better parameterizations for individual physical processes. The modules have
been designed, where possible, so that the analytical or numerical calcula-
tions in the modules can be replaced by observational data, when available,
with minimal modifications in the computer code structure. For example, if
acoustic sounder data are available to measure CBL height and inversion top
height, these data can replace the analytical scheme for predicting these
heights.
49
-------
The modular structure of the model is illustrated in Figure 19, where the
modules are named and the functions of the modules are indicated. In this
section we give the model user an overview of the VALMET model, explaining its
structure and its input requirements. A full technical description of the
modules, including variable name definitions and module inputs and outputs is
given in Section 5 for the interested computer programmer. In Section 6 the
model outputs are illustrated for two sample simulations. The full VALMET
code is given in Appendix A. A brief historical summary of modifications to
the code is given in Appendix 0.
FEATURES OF THE COMPUTER CODE
The VALMET model is coded following a FORTRAN 77 standard [41] and should
be useable without modification on any computer system having an up-to-date
FORTRAN compiler. The model is documented internally through the liberal use
of comment statements. These comments explain the purpose of sections of code
or individual FORTRAN statements, or serve as variable name definition tables.
Several special features or protocols are included in the code. The
dimensions of arrays in the code are generally set using PARAMETER state-
ments. Thus, if the user finds that the array sizes are too small for the
length of the simulation or too large to fit in a smaller computer he may
easily change the dimensions by modifying PARAMETER statements in the main
program and in selected subroutines. COMMON blocks are fully utilized in the
code to reduce the memory requirements of the operating progam. Most of the
parameters passed to the subroutines are generally passed through the sub-
routine argument lists. Arrays, however, are usually passed through COMMON
blocks.
The model was developed on a VAX/VMS 11/780 computer. After compilation,
a single model run generally takes 1 minute of central processor time on the
VAX, loads in 95,000 bytes (25K words) of core, and costs approximately $1.00.
Costs of optional plotting will vary from installation to installation.
MODEL INPUTS
The VALMET program runs in an interactive mode in which the user controls
program execution from a remote interactive terminal. The user enters model
inputs by following directions given to the user on the terminal screen.
As with any model, model performance will be a function of the suitabil-
ity of the input data entered by the user. Where worst-case values of the
input parameters are estimated on the basis of few data, the user should
consider the model results to be a screening analysis only. On the other
hand, an industrial installation having a great deal of pollutant concentra-
tion and meteorological data could use the model as a site-specific model by
making specific modifications to the model to include observed plume dilution,
dispersion coefficients, etc. This would require considerable field experi-
mentation and data analysis.
50
-------
PROCESSES
PLUME RISE
DILUTION
POLLUTANT MASS
CONSERVATION
STEADY STATE TRANSPORT
AND DIFFUSION
I
EBDGT
DESCNT
SOLAR FORCING
SURFACE ENERGY BUDGET
INVERSION DESCENT AND
CBL GROWTH
Figure 19.
Flow diagram of VALMET model, showing the modular structure of
the model and the physical processes parameterized in the
modules. The modules indicated between the two-horizontal
dashed lines constitute the nocturnal portion of the model.
The daytime portion of the model follows the second dashed
line.
51
-------
y
N=0
"t ,N=N+1
NO -X"
f iv
PROFIL
\
VELOCY
I
BRKUP
A
l-IMCTPDC
__ GAUSS
FU
> ^s
FUMIGATION
WRITE
SUMMARY/
*\ OUTPUT,
FILE
WRITE
OPTIONAL
•^.PLOTTING/
FILE
\
POLLUTANT MASS FUMIGATION
TOTAL MASS CONSERVATION
AND ADVECTION
POLLUTANT MASS CONSERVATION
MULTI-HOUR-AVERAGE
CONCENTRATIONS
Figure 19. (contd)
52
-------
When the user runs VALMET from an interactive terminal, a data table
appears on the user's screen. This table gives sample input values for a
model run. The user may follow directions given on the screen to change any
of the input values to simulate air pollution dispersion in any valley of
interest, or he may simply specify that the sample (or default) simulation be
run.
The input table (Figure 20) is arranged according to categories of input
data. Table input values for the sample simulation are generally specified in
the MKS system of units, although pressures are given in millibars and
parameters involving potential temperatures are given in degrees Kelvin in
conformance with general meteorological practice.
Table 5 below gives a full listing of the appropriate units of all input
parameters used in VALMET as well as the default values of these parameters.
VALMET is programmed so that the user can change input values as many times as
necessary until the input table is properly completed. Once the input table
is satisfactory to the user, he instructs VALMET to proceed with the model
run.
SPECIFICATION OF THE MODEL INPUTS BY THE USER
As a convenience to the user, the VALMET program checks to see that any
parameter values input by the user are within a normal atmospheric range of
values for that parameter in the units required by VALMET. See Table 6 for a
listing of the ranges of values assigned to each parameter. If the input
value is outside normal atmospheric ranges, a message is sent to the user to
respecify the input value. The intent of this feature is to help the user
identify gross input errors caused, for example, by using the wrong units for
input variables.
In this section of the report we will instruct the model user how to
determine the input parameters that are required by the model. The parameters
are entered into the model by the user as he follows instructions given on the
interactive terminal screen. The input parameters are discussed in groups, as
follows:
Valley Characteristics
The characteristics of the valley topographic cross section must be
obtained from topographic maps. The user should plot the terrain cross sec-
tion at the down-valley distance of interest. From this cross section the
user should estimate the valley floor width (estimates to the nearest 100 m
should be sufficient) and the two sidewall inclination angles. The user must
use his judgment when making these estimates, depending on the representative-
ness of the topographic cross section plotted. Since the model will deal with
the valley sidewalls only to the initial depth of the temperature inversion,
the user should ignore the sidewalls that extend above this height when esti-
mating representative sidewall angles.
53
-------
• O CM O O O CO
in o co o o o Ed
r- O . • • C5
. ,=r .=r on z
o <
x
o
EM 0>
O -=
CO JEdXCJCOCC CC c
Ed «a!CQMO3eO Ed o
a CQ
S ^
2: a;
Q.
• a.
II
CM
<
X
CU
J
«a!
in
in
n
«s
H
Ed
ca
<*
in
CM
n
CC
>H
M
*
O «-
O CM
II II
CJ CO
a a
cr> CM
<- CM
o •
o o
n
Q
<4
CC
CO
in
CM
O •
O O
TO THE FOLLOWINi
EH
Ed
Ed
CC
CO
CC
Ed
EH
Ed
S
CC
cu
z
o
M
N
M
t •)
•a!
H
EH
M
M
Z
•a!
ce
o
o
cc
a,
Ed
33
EH
in
o
n
CM
CM
e
O
II
o
ai
in
n
33
CM
<
.=r
0
0
VD
II
33
EH
Q
H
3
in
CM
o
o
n
•
CC
Ed
EH
• O
CC *aj
Ed CC
H «J
CO 35
«! O
CC
•a! Z
33 0
CJ H
CO
J CC
Ed Ed
Q >
0 Z
S M
1 URLONG=105.00
I IDAY = 21
«- on
O O>
0
«
0
jy
II II
EH
•a! O
J S
O CM
»- t-
z
0
H
EH
«aj
O
O
J
Ed Ed
H £H
H <
CO Q
O •
O O
O CM
in
c-
II II II
CO O
CO CC CQ
Ed Ed <
CC N EH
CM >M CO
vO CO ^
T- T- f\j
O « O
O 0 i-
• O O
000
t- o •
T- O
II II II
cu
s
Cd 0
EH X O
in c— o
r~ r- CM
Ed
CC
Ed
33 Ed
*** i«
CO 0
O J
S CU
EH
aj z
<$
SH M
Cd CO
J CO
J »
«< <£
> O
O •
O O
0
o
II II
2
EH CO
CO «B
sr i>-
CM CM
O 0
o o
• »
0 ^T
in CM
CM
II II
J
Ed
33 CO
ro vo
CM CM
•
(X
Ed
E-i
CJ
•a
OS
•a
33
O
M
CJ
•a!
EH
CO
cn
a:
Ed
EH
Z
Ed
1
1
1
Ed
•1
>— 1
ca
H
CO
Ed
p-i
J
•a;
{>
Sg
M
Ed
O
^
CQ
O
g-t
Cd
CC
*3j
CO
Ed
O
y%
«3J
X
O
o
z
EM
M
co
CC
Ed
EH
y*
Ed
1
1
1
1
0
X
CJ
Ed
CQ
0
EH
CO
M
S
CQ
c-
**—
O
c
o
4-1
tO
•!->
I/I
3
r—
'""
1—1
•
O
CM
CL)
^
C71
>r~
lj-
.
CJ
0)
{_
0
in
OJ
>
•i—
-^J
O
2
O)
c
•r—
to
c_
d)
to
3
1
54
-------
TABLE 5. DEFAULT VALUES OF VALMET INPUT PARAMETERS
Variable
Name
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.
AO
K
WIDTH
ALPHA1
ALPHA2
NBOX
HZERO
GAMMA
BETA
LAT
URLONG
MO
IDA
IYR
TEMP
PRESS
X
YZERO
UC
Q
STAB
US
H
STMP
SRAD
SVEL
AS
Definition
Fraction of solar flux converted
to sensible heat flux
Fraction of sensible heat flux used
to grow CBLs
Val ley floor width
Sidewall #1 elevation angle
Sidewall #2 elevation angle
Number of grid elements on valley
floor half-width
Initial inversion depth
Inversion vertical potential
temperature gradient
Rate of potential temperature
increase above valley
Latitude
Longitude
Month
Day of month
Year
Ai r temperature
Air Pressure
Down-valley distance from stack
Stack offset distance from valley
floor centerline
Plume-carrying wind speed at
receptor cross section
Stack emission rate
Atmospheric stability
Plume-carrying wind speed at
source cross section
Plume centerline height
Stack emission temperature
Stack radius
Stack emission velocity
Cross-sectional area at source
Default
Units Value
0-1
0-1
m
o
o
-
m
°K/m
°K/s
°N
°W
1-12
1-31
00-99
°C
mb
m
m
m/s
kg/s
D-F
m/s
m
°C
m
m/s ,,
10J nr
0.24
0.15
600
15
15
3
500
0.025
0
40
105
9
21
82
10
750
10000
0
4
0.001
F
4
250
100
3
0
0
Date
The user must specify the date of the simulation by giving the month, day
and year. The date is necessary in the solar model to calculate day length,
time of sunrise, etc. The year is necessary to account for the minor effect
of leap years on the simulation.
55
-------
TABLE 6. LIMITS ON THE VALUES OF INPUT PARAMETERS
Variable
Name
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.
AO
K
WIDTH
ALPHA1
ALPHA2
NBOX
HZERO
GAMMA
BETA
LAT
URLONG
MO
IDA
IYR
TEMP
PRESS
X
YZERO
UC
Q
STAB
US
H
STMP
SRAD
SVEL
AS
Units
_
-
m
o
o
-
m
°K/m
°K/s
°N
°W
-
-
-
°C
mb
m
m
m/s
kg/s
D-F
m/s
m
°C
m
m/s 0
10V
Lower Limit
0
0
0
0
0
1
0
0
0
-90
-180
1
1
0
-50
600
100
-50000
0
0
D
0
0
0
0
0
0
Upper Limit
1
1
50000
90
90
30
2000
0.050
0.0005
90
180
12
31
99
50
1050
100000
50000
15
1
F
15
1000
600
5
50
150000
Site Location
Both latitude and longitude of the site must be
order to support solar radiation calculations within
value is necessary to account for corrections to the
Model Characteristics
specified by the user in
the model. The longitude
time of sunrise.
The user must specify the number of grid elements on the valley floor
half width in order to define the model's numerical
considerations that the user should keep in mind when
number. The implications of the specification can be
to Figure 12.
grid. There are several
specifying this
considered by referring
56
-------
The number of grid elements on the valley floor half width (NBOX) is used
in the model to figure the length of the grid elements. The length is given
by
BOXLEN=w/(NBOX*2)
where w is the valley floor width. Once BOXLEN is calculated, the model then
calculates how many grid elements of this length (NTOTI) will be necessary in
order to reach the top of the inversion when the grid elements are arrayed end
to end. NTOTI must not exceed 30, as currently dimensioned in the model. On
the other hand, NTOTI should not be too small, since this parameter affects
the spatial resolution of concentrations calculated within the model. As a
rule of thumb, NBOX may be chosen to produce grid elements of 100 or 200 m
length.
Valley Atmosphere
The mean temperature and pressure within the valley temperature inversion
at sunrise should be estimated using the altitude of the site and knowledge of
seasonal changes in atmospheric temperature. Observational data would be use-
ful in estimating both of these parameters, but the accuracy of these
estimates is not critical to the model results. The values are used to
correct the model concentration values to standard pressure and temperature
and to calculate the atmospheric density for use in the inversion breakup or
energy budget calculations.
Stack Characteristics
The stack characteristics should be obtained from engineering calcula-
tions for the specific air pollution source. The required parameters include
physical stack height, effluent temperature, stack inside radius, and exit
velocity. If the user does not wish to engage the plume rise module, the
stack radius or exit velocity should be specified as zero. The plume center-
line height over the entire length of the valley will then be equal to the
value specified as the physical stack height. The final parameter that must
be specified by the user is the area of the valley cross section at the
stack. This area, bounded by the valley floor, the two sidewalls and the
temperature inversion top at sunrise, is used in the plume dilution module to
calculate clean air volume flux across the cross section at the pollutant
source. The area should be calculated from Equation (1), but is input in
thousands of square meters (i.e., the results of the calculation with
Equation (1) should be divided by 1000 before being input into the model).
The model user will need to know the inversion depth at sunrise in order to
calculate the area. If the user does not wish to engage the plume dilution
module or does not have the wind speed observations on the cross valley
section where pollution calculations are to be made, he should specify the
cross-sectional area at the stack to be zero. The wind speed specification at
57
-------
the downwind cross section is then deactivated (i.e., is not used in the
model) and the plume is advected down the valley with the wind speed specified
at the pollutant source with no plume dilution during transport.
Inversion Characteristics
The user must specify inversion depth and inversion potential temperature
gradient at sunrise, as well as the rate of potential temperature increase at
the inversion top during the period of inversion breakup. The user should
specify these parameters on the basis of observations of temperature inversion
development within the valley of interest. These observations are critical to
the model results and should be based on vertical soundings taken from the
valley floor, if available. The modeler should be cautioned against using
observations taken from surface-based instruments, unless these have been
shown to be clearly related to vertical profiles over the valley center. A
qualified analyst may be able to determine the inversion depth from acoustic
sounder records. The modeler must use his judgment in specifying a single
number for the potential temperature gradient at sunrise. Actual soundings
are often complicated, especially in the lower levels. The potential
temperature gradient obtained by subtracting the surface value from the value
at the inversion top and dividing by the inversion depth may not be a
representative number for the inversion as a whole when a shallow and very
strong inversion layer is present near the ground. The temperature gradient
is used in the model primarily to calculate the energy deficit within the
valley represented by the presence of the temperature inversion. The total
energy deficit is obtained under an assumption of cross valley homogeneity by
weighting the temperature deficits at various levels by the cross-sectional
area of the valley. Since the valley is much wider at its top than near the
valley floor, errors made in estimating temperature deficits (i.e., deter-
mining a representative potential temperature gradient) at the upper levels
are more critical to model results than at the lower levels. A better esti-
mate of the valley energy deficit will be obtained by ignoring shallow, but
very strong, temperature inversion layers present on the valley floor if this
means that the potential temperature gradient estimate is more accurate for
upper levels.
In lieu of observations one may use the observations of Whiteman [9] to
obtain first estimates of inversion depths, strengths, and inversion top
warming rates in Colorado valleys.
When the model is used to estimate worst-case air pollution concentra-
tions, one should naturally attempt to arrive at worst-case estimates of
inversion characteristics parameters.
Gaussian Plume Parameters
The model user must specify Gaussian plume parameters including the dis-
tance to the cross section of interest, the stack distance offset from the
valley floor centerline, the plume-carrying wind speed, the source emission
rate, and the atmospheric stability.
58
-------
The down-valley distance is determined from a topographic map by
measuring the distance from the source to the cross valley section along the
valley floor center!ine following any curves or turns in the valley's course.
The stack distance offset is determined as the shortest distance from the
valley floor center!ine to the pollutant source on the valley floor. The sign
of the distance follows from the orientation of the valley coordinate system,
as discussed in Section 3.
The wind speeds used in VALMET are critical parameters that must be
carefully specified from available data. Down-valley wind speeds in Colorado
valleys have been shown to vary significantly from valley to valley [9], so
that wind speeds in one valley are not necessarily indicative of wind speeds
in another valley. In short, observations are necessary for the valley being
modeled. Furthermore, the wind speeds vary on a cross valley section due to
vertical and horizontal position within the cross section [42]. On a vertical
profile through the nocturnal temperature inversion, the peak winds are often
located about midway through the inversion depth with lower wind speeds near
the ground and at the inversion top. Near-surface wind speeds, therefore, are
not generally indicative of wind speeds at a typical plume carrying level.
Near-surface wind direction observations, in fact, may be 180 degrees in error
shortly after sunrise when wind continues to flow down-valley in the elevated
stable core while the winds reverse to up-valley in the C8L which develops
over the valley floor. Near-surface wind observations should, therefore, not
be used to drive the VALMET model.
Two wind speeds are required inputs for the VALMET model. Both are
required in order to allow for plume dilution due to clean air inflow between
the source and the cross section where air pollution calculations are made.
The user must specify the average nocturnal wind speed at plume-carrying level
at both the source and the receptor characteristic of the nocturnal down-
valley flow. They may be obtained from a doppler acoustic sounder or from
successive tethersonde or pilot balloon ascents.
The source emission rate may be obtained from engineering calculations
for the stack of interest. "Worst-case" pollutant calculations should be made
on the basis of worst-case emission rates. The model can not handle emission
rates that vary in time, so a realistic time average emission rate is most
appropriate for driving the model.
The model uses Pasquill-Gifford diffusion coefficients [26,27,43] to
characterize diffusion in the nighttime, high-stability, elevated stable core
of the valley temperature inversion. The stability class in this stable
region should be D, E, or F following the original definitions of atmospheric
stability given in Table 3. The user-input stability class is automatically
modified internally within VALMET to account for enhanced diffusion in complex
terrain.
59
-------
Sensible Heat Flux
The final two input parameters in the VALMET model are the two sensible
heat flux parameters. The parameters, which are both fractions between 0 and
1, are based on the valley surface energy budget. Little research has yet
been done on valley energy budgets and, since the energy budget is known to
strongly influence the local diurnal meteorology of valleys, the initiation of
such research should be considered an important priority in future studies.
The first sensible heat flux parameter is A0, the fraction of extra-
terrestrial solar heat flux converted to sensible heat flux within the
valley. The conversion, of course, occurs at the valley surfaces (valley
floor and sidewalls), but the bulk thermodynamic model used as the basis of
VALMET cannot distinguish between energy budgets over the different sur-
faces. The fraction required in the model, then, is the fraction of the
(extraterrestrial) solar heat flux coming downward across the upper surface of
the temperature inversion which is converted to sensible heat flux within the
valley. This fraction is assumed to be constant during the temperature
inversion breakup period. AQ depends on the transmissivity of the earth's
atmosphere to the solar beam, the albedo of the valley surfaces and vegetative
cover, and the presence of clouds, as well as the partitioning of the solar
beam upon striking the valley surfaces. The long wave radiative balance also
plays a role at these surfaces which affects sensible heat flux. Whiteman [9]
has analyzed the individual terms in the surface energy budget to determine
that the maximum value of AQ is unlikely to be greater than 0.6. A valley
having a low reflectivity would approach this value under clear skies only if
the valley had no soil moisture. Whiteman and McKee [3] were able to simulate
temperature inversion breakup well in a dry Colorado valley in October using
an AQ value of 0.45. Maximum values of heat flux observed in the 1968 Kansas
experiments [44], if assumed to be representative of solar noon values, corre-
spond to AQ values of about 0.25. The value of A may be smaller in winter
when snow cover is present in a valley and the albedo is large. The overall;
albedo in winter may depend strongly on the presence of forest cover within
the valley or the relative proportions of other lower albedo surfaces.
Whiteman and McKee [3] were able to simulate inversion breakup in Colorado's
snow-covered Yampa Valley using AQ = 0.19.
The modeler may estimate a value of AQ from the discussion in the para-
graph above or, if sufficient meteorological data are available for the valley
in question, he may base his estimate on one of two alternate quantitative
methods. The first method is described below. Then the second, and more com-
prehensive method, which also provides a means for estimating the second
sensible heat flux parameter, k, is described. The first method of determin-
ing A involves the solution of Whiteman and McKee's [3] Equation (21), Equa-
tion (28), or a combination of the two for AQ knowing, for a particular day::
a) the sunrise temperature inversion potential temperature gradient,
b) the inversion depth at sunrise,
c) the valley floor width and sidewall inclination angles,
d) the solar parameters on the day of interest, and
e) the observed time of temperature inversion breakup.
60
-------
By using such meteorological information for a variety of days in dif-
ferent seasons, the modeler may determine a climatology of A0 values that
would be useful in air pollution modeling. For the convenience of the reader,
a reprint of the Whiteman and McKee article is included in this report as
Appendix C.
The value of the second sensible heat flux parameter, k, is the fraction
(0-1) of the sensible heat flux used to cause CBL growth. The remaining frac-
tion (1-k) represents the fraction of sensible heat flux used to transport
mass up the valley sidewalls causing the inversion top to sink. The value of
k is 1 for flat terrain since temperature inversions are destroyed there
entirely by CBL growth, but approaches zero in snow-covered valleys where CBL
growth is arrested after reaching a certain height [9], Whiteman and McKee
[3] were able to simulate inversion breakup in a dry Colorado valley using k =
0.14. Inversion breakup in the snow-covered Yampa Valley was simulated well
using k = 0.
The search for a mathematical relationship between k and parameters that
are more easily measured is a high priority for future work. The value of k
is expected to depend strongly on valley width and, perhaps, on sensible heat
flux itself. Given the present lack of an appropriate mathematical relation-
ship, the best approach for determining k is to simply determine the value of
k which best fits observations of inversion breakup in the particular valley
being modeled. The means by which this can be accomplished was described by
Whiteman and McKee [3]. They described how the value of AQ determines the
temperature inversion breakup time, while the value of k determines the height
at which the CBL and inversion top meet to cause inversion breakup, and have
presented plots showing the effect of varying k and AQ on temperature inver-
sion breakup. Both the time and height can be observed using an acoustic
sounder or multiple serial tethersonde ascents in the valley of interest.
Observed time variations of CBL depth and inversion top height can be fit with
a trial and error procedure using Whiteman and McKee's [3] Equations (29) and
(30). The reader should refer, again, to Appendix C. These equations are
programmed in VALMET module "OESCNT" where they may be solved using a
numerical method. The modeler who wishes to use the procedure should convert
a copy of subroutine DESCNT to a main program where the subroutine arguments
are explicitly specified. The number of time steps to inversion destruction,
and the CBL and inversion top heights as a function of time would be required
outputs of such a program. The effect of changing AQ and k could then be
easily investigated and fit to actual observations. A time step of
600 seconds would be sufficient for the computations desired. The numerical
method requires a CBL height and inversion depth at sunrise as mathematical
initial conditions. The valley meteorologial data necessary to use this
approach to determine k and AQ include:
a) inversion depth at sunrise,
b) inversion potential temperature gradient at sunrise,
c) valley floor width and sidewall angles,
d) solar model parameters A]_ and T,
61
-------
e) the time of inversion destruction, and
f) the height at which the rising C8L and descending inversion top meet
at the time of temperature inversion destruction.
Solar parameters can be determined from model equations provided earlier in
this report. Topographic parameters can be determined from a topographic
map. The potential temperature gradient at sunrise requires a vertical tem-
perature sounding. The other required parameters can also be obtained from
such soundings, if they are made frequently through the valley depth during
the temperature inversion breakup period. They may also be determined from
monostatic acoustic sounder records by a qualified analyst.
ERROR MESSAGES
The VALMET model has several built-in checking and correction routines
which will produce error messages when the model input parameters have been
improperly specified by the user, or when the user attempts to use the model
in a way unintended by the model developers. These messages are written on
the screen of the user's interactive terminal. The error messages and their
explanations are given below.
1. PARAMETER NO. nn OUT OF RANGE, PLEASE RESPECIFY
This error message is written from Subroutine INPUT when the user has
tried to change the default value of an input parameter (number nn=l to 27),
but has selected a new value that is outside normal ranges of that parameter.
The user should check that he has used the right units for the input specifi-
cation. Execution of the program is not terminated by this error. The user
may respecify the input parameters as many times as he wishes before executing
the VALMET calculations.
2. STACK MUST BE ON THE VALLEY FLOOR - RESPECIFY YZERO
This error message is written from Subroutine INPUT when the user
specifies a stack location that is not on the valley floor. VALMET was not
designed to handle this situation, so the user must respecify the stack
location. To do this, respecify YZERO to be smaller in absolute value than
half the valley width W. Execution of the program is not terminated by this
error.
3. DOWNVALLEY WINDSPEED AT STACK HAS BEEN SET TO 1 M/S
This error message is written from Subroutine INPUT when the user has
specified a very low (0-1 m/s) plume-carrying wind speed at stack height. The
nighttime VALMET simulation has, as its basis, the Gaussian plume formulation,
which fails to predict concentrations realistically as winds become calm
(us-K)). To obviate this difficulty, VALMET automatically respecifies the wind
speed to 1 m/s and continues execution. Mile the basis for this modification
to the wind speed is primarily mathematical, one must recognize that perfectly
calm conditions are unusual in the valley atmosphere due to the generation of
62
-------
the prevalent locally developed circulations. The frequency of occurrence of
calm conditions has been overestimated in the past, and is still being over-
estimated by the performance characteristics of available mechanical wind
sensors.
4. TOO MANY GRID ELEMENTS — PLEASE REDUCE NBOX OR RE-DIMENSION THE MODEL
TO ALLOW MORE GRID ELEMENTS
This error message is written from Subroutine INPUT when the input
parameters specified by the user result in the model requiring more grid
elements than are currently dimensioned in the model (i.e., 30). The total
number of grid elements, NTOTI, may not exceed 30 unless the user increases
the dimensions in the model by changing the PARAMETER statements so that NB is
the new number desired (don't forget to also change NB1, since it should be
equal to NB+1). If your computer memory is limited, you must decrease NBOX
and live with the poorer spatial resolution which results, as documented at
other places in this report. Execution of the model is not terminated by this
error. If you reduce NBOX the program will continue to execute, but if you
decide to increase the model dimensions you must stop execution and make some
basic changes in the VALMET code.
5. YOU HAVE SPECIFIED A WRONG STABILITY CLASS - F WILL BE USED INSTEAD
This error message is written from Subroutine INPUT when the user enters
a stability class other than D, E, or F. «VALMET automatically uses F
stability when this occurs. Execution of the program is not terminated by
this error.
6. THE PLUME, AFTER PLUMERISE, IS NOT WITHIN THE STABLE CORE
This error message is written from the main program. If the plume rises
above the stable core, or never exceeds the initial CBL height, it is not
transported down the valley in stable nighttime conditions and it makes no
physical sense to run VALMET under these conditions. Model execution is
terminated by this error.
63
-------
SECTION 5
TECHNICAL DESCRIPTION OF INDIVIDUAL MODULES
VALMET-Main Program
VALMET is the main program that controls the modules or subroutines
which, taken-together, form the valley air quality model. The function of the
main program is to provide the basic structure of the air quality model and to
call the specific modules as required. In addition, several computing house-
keeping functions are performed in the main program including the establish-
ment of PARAMETER statements and COMMON blocks, and the opening and closing of
output files. Model outputs are written from the main program, as well.
Inputs.
Inputs to the main program come through Subroutine INPUT, a subroutine
that is run interactively by the user and which features a default input
table. Data coming into the main program from subroutines usually come
through arguments listed in the subroutine call statement. Large arrays of
data, however, are frequently transferred from the main program to the sub-
routines (and vice versa) through COMMON blocks.
Outputs
The results of a model run are written to two output disk files. The
first output file, named VALMET.OUT, is a formatted file which contains the
results of the model run. After the model is run, the contents of this file
may be printed to obtain a summary record of the model run, including model
input parameters and outputs. A second file, named VALMET.PLT, is also auto-
matically generated with every computer run. This formatted file contains
time series of pollutant concentrations, temperature inversion depths and con-
vective boundary layer depths, as well as many of the basic parameters used in
the model run. The file is created for the model user who wishes to develop
his own plotting programs to plot the results of a model run. An example pro-
gram that the authors used to generate the plots used in this report is shown
in Appendix B. Since every user's computer installation will have different
plotting software and hardware, this approach of creating a basic data file
but not including a specific plotting program in the air pollution model
seemed the most appropriate way to proceed. Interested users can develop
their own specific plotting software using our example.
64
-------
Conventions
The following indexing conventions are used in VALMET where it has proven
convenient:
I=grid element index
Variable Name Definitions
INPUT/ VARIABLE
OUTPUT NAME
N or MN=time step index
DEFINITION
AO
Al
A2
AC
ALPHA1
ALPHA2
AS
BETA
BOXLEN
CCHI
CHIOFF
-
-
I
-
I
I
-
I
-
-
-
I
-
CLCONC
DELTAT
DLH
DF
GAMMA
H
HCBLI
HZERO
I
11
IAVG
IDA
IFIX2
I NO
ISTAB
ISTABY
ISTABZ
IT
ITIME
Fraction of extraterrestrial solar flux converted to
sensible heat flux (0-1)
Solar flux on horizontal surface at solar noon (Wm )
Same as AO
Valley cross-sectional area below inversion top
at cross-sectional distance XC (m )
Sidewall #1 elevation angle (rad)
Sidewall #2 elevation angle (rad)
Valley cross-sectional area below inversion top
at pollution source cross section (m )
Rate of increase of potential temperature at
inversion top (°K/s)
Length of grid elements (m)
Nocturnal pollutant concentration in a grid element
Concentration offset value due to reflection off
valley floor and sidewalls (ng/m )
Center!ine concentration (ng/m )
Time step size (s)
Plume rise (m)
Inverse of dilution factor
Vertical potential temperature gradient (°K/m)
Plume center!ine height (m)
Initial CBL height (m)
Initial depth of inversion (m)
Grid element index
Intermediate variable
Averaging period (s)
Day of month (1-31)
Number of grid elements on sidewall at given time
step
Loop index
Stability index (1-6)
Stability index for determination of a (1-6)
Stability index for determination of o^ (1-6)
Intermediate value used in calculating time
Time unit=LST
65
-------
IYR
JULDAY
K
LAT
MN
MO
N
HAS
NB
MB1
NBOX
ND
NINDX
MS
NSTEPS
NTOT
NTOTI
NTSA
NTSR
NY
PRESS
PSTD
Q
RHO
SIGMAY
SIGMAZ
SRAD
STMP
STNDCOM
SUM1
SUM2
SVEL
Tl
T2
T3
T4
TAU
TEMP
TMAXC1
averaging intervals in NINDX
giving the maximum number of model grid
for use in dimensioning arrays
Year (00-99)
Julian date (1-366)
Fraction of sensible heat flux used to grow CBL
(0-1)
Latitude (°M)
N-l
Month (1-12)
Time step counter
Number of
Parameter
elements,
NB+1
Number of grid elements on the valley floor
half-width
Loop index
Number of time steps in simulation, including
exponential decay after breakup
Parameter giving the maximum number of model time
steps, for use in dimensioning arrays
Number of time steps required to destroy initial
temperature inversion
Total number of grid elements left in the model at
a given time step
Total number of grid elements in the model at
initiation
Number of time steps in averaging interval
Time of sunrise (LST)
Intermediate variable used to count grid elements
Atmospheric pressure at center of temperature
inversion (mb)
Standard pressure (1013 mb)
Stack emission rate (kg/s)
Air density (kg/m )
Sigma y (m)
Sigma z (m)
Stack radius (m)
Stack temperature (°K)
Factor to adjust concentration to standard conditions
Fraction of plume mass within valley cross section
(0-1)
Fraction of plume mass above valley cross section
(0-1)
Stack exit velocity (m/s)
Intermediate variable used to
Intermediate variable used to
Intermediate variable used to
Intermediate variable used to
Length of daylight period (s)
Air temperature (°K)
Ending time of the 1-hr average
(h LST)
calculate
calculate
calculate
calculate
time
time
time
time
concentration max
66
-------
TMAXC2 Ending time of the 3-hr average concentration max
(h LSI)
TSTD Standard temperature (293.16°K)
TSTP Time in averaging period (s)
I UC Mean down-valley plume transport velocity at cross
section distance XC (m/s)
I URLONG Longitude (°W)
I US Mean down-valley plume
pollutant source cross
I W Valley floor width (m)
I XC Distance from stack to model cross section (km)
Y Y-coordinate (m)
I YZERO Stack offset from valley floor centerline (m)
Z Z-coordinate (m)
transport velocity at
section (m/s)
Array Name Definitions
INPUT/ ARRAY
OUTPUT NAME
DEFINITION
AVG(NB,2)
CCH(2,NB)
CHI(NB)
CHIBAR(NB)
CONC(NS,NB)
HC(NS)
HITEL(NB)
HITEU
HT(NS)
NDXTIM(NB,2)
NTS(2)
SUM(NB)
TIME(NS)
V(NB1)
Maximum 1-hr and 3-hr average concentration array
Ug/m3)
Short-term storage location for concentrations
at two adjacent time steps d-ig/m )
Nocturnal steady state concentration array
array (p-g/m )
grid element on sidewall
grid element on sidewall
Average pollutant concentration injected into
the top of grid elements at each time step
Ug/m3)
Pollution concentration
Height of CBL top (m)
Height of lower side of
(m)
Height of upper side of
(m)
Height of inversion top (m)
Index number of TMAXC1 and TMAXC2
Number of time steps in 1-hr and 3-hr
Summing array to accumulate pollutant
concentrations during averaging interval
Midpoint time array for averaging intervals
(h MST)
Slope flow velocity array (m/s)
Subroutines Called
BRKUP
DESCNT
DILUTE
67
-------
EBDGT
GAUSS
INGRAT
INPUT
JULIAN
PRISE
PROFIL
PSTPRC
SOLAR
VELOCY
Common Blocks
BLK1
BLK2 Not used in main program
BLK3 Not used in main program
BLK4
BLK5
BLK6
68
-------
INPUT-Input Module
Purpose
The VALMET program runs in an interactive mode in which the user controls
program execution from a remote interactive terminal. The user enters model
inputs from the remote terminal by following directions given to the user on
the terminal screen by subroutine INPUT. Once the input values are specified
by the user he can direct the program to begin the air pollution simulation.
The functions of the INPUT module are to obtain the proper input data from the
model user, to convert the input data to the proper units for processing in
later modules of VALMET, to check the input data for errors, and to notify the
user of any errors.
Inputs
Inputs to this module come entirely through user input from an inter-
active terminal, unless default values of the input parameters are selected by
the user.
Outputs
Outputs from the module go to the main program through the subroutine
argument list.
Variable Name Definitions
INPUT/ VARIABLE
OUTPUT NAME
DEFINITIONS
I
0
I
I
I
I
I
0
0
I
AO
AC
ALPHA1
ALPHA2
AS
BETA
BOXLEN
DELTAT
GAMMA
H
HCBLI
HZERO
IDA
IFIX
ISTAB
IYR
J
J6
Fraction of solar flux converted to sensible heat
flux (0-1)
Valley cross-sectional area at distance XC (m )
Sidewall #1 elevation angle (°)
Sidewall #2 elevation angle (°)
Cross-sectional area at source (m )
Rate of potential temperature increase above valley
temperature gradient
Grid element length (m)
Time step (s)
Inversion vertical potential
Plume centerline height (m)
CBL height at sunrise (m)
Initial inversion depth (m)
Day of month (1-31)
Number of grid elements on sidewall
Atmospheric stability (1-6)
Year (00-99)
Index
Index
69
-------
J12 Index
J13 Index
J14 Index
I K Fraction of sensible heat flux used to grow C8L (0-1)
I LAT Latitude (°N)
LU Logical unit number
I MO Month (1-12)
I NBOX Number of grid elements on valley floor half width
NTOTI Total number of grid elements
NU Number of input parameters to be changed from default
values
I PRESS Air Pressure (mb)
I Q Stack emission rate (kg/s)
RHO Air density (kg/nv3)
RNO NBOX
I SRAD Stack radius (m)
I STMP Stack emission temperature (°C)
I SVEL Stack emission velocity (m/s)
I TEMP Air temperature (°C)
I UC Plume-carrying wind speed at cross section (m/s)
I URLONG Longitude (°W)
I US Plume-carrying wind speed at source (m/s)
I w Valley floor width (m)
I XC Down-valley distance to cross section from stack (m)
I YZERO Stack offset distance from valley floor centerline
(m)
Array Name Definitions
INPUT/
OUTPUT
_
-
-
-
-
-
-
-
ARRAY NAME
ASTAB(6)
CAT(IO)
10(27)
NAM(27)
VAL(27)
VALLO (27)
VALHI (27)
VALNEW(27)
DEFINITIONS
Atmospheric stability array
Input parameter category array
Array of identification numbers for
parameters
Input parameter name array
Default input value array
Minimum values of input parameters
Maximum values of input parameters
New values of input parameters
input
Subroutine Called
None
Common Blocks
None
70
-------
JULIAN-Julian Day Module
Purpose
JULIAN calculates the Julian date, given the month, day and year of the
air pollution simulation.
Inputs
The month, day, and year are input to the subroutine from the main pro-
gram through the subroutine argument list.
Outputs
The Julian date is sent to the main program through the subroutine
argument list.
Variable Name Definitions
INPUT/ VARIABLE
OUTPUT NAME
DEFINITION
I IDA
I IYR
0 JULDAY
I MO
Leap year indicator
Day of month (1-31)
Year (00-99)
Julian date (1-366)
Month (1-12)
Array Name Definitions
INPUT/
OUTPUT ARRAY NAME
DEFINITION
NDAY(12)
Cumulative days of year, by month
Subroutines Called
None
Common Blocks
None
71
-------
PRISE-Plume Rise Module
Purpose
PRISE calculates plume rise using the Briggs plume rise algorithms docu-
mented in the MPTER User's Manual, taking account of both plume momentum and
buoyancy effects.
Inputs
The inputs to PRISE come through the subroutine argument list and include
stack characteristics (radius, exit velocity, exit temperature), ambient atmo-
spheric characteristics (air temperature, wind speed, stability class, poten-
tial temperature gradient) and down-valley distance to the cross section of
interest.
Outputs
The output of PRISE is the plume centerline rise above the physical stack
height. This value is passed to the main program through the subroutine
argument list.
FORTRAN Library Subroutines
Name
Function
AMINl(xl,x2,...xn)
Determine the minimum value of a list
of real numbers
Variable Name Definitions
INPUT/ VARIABLE
OUTPUT NAME
DEFINITION
DELTT Plume temperature excess (°K)
0 DLH Plume centerline rise above stack height (m)
DLH1 Intermediate plume rise variable
DLH2 Intermediate plume rise variable
DTC Crossover temperature difference (°K)
F Buoyancy flux (ms ) „
G Acceleration due to gravity (ms )
I GAMMA Potential temperature gradient (°K/m)
I ISTAB Vertical atmospheric stability index (3-6)
S Stability parameter (s~ )
I SRAD Stack radius (m)
I STMP Stack gas exit temperature (°K)
I SVEL Stack gas exit velocity (m/s)
I TEMP Ambient air temperature (°K)
UCUT Critical wind speed (m/s)
I US Ambient plume-carrying wind speed (m/s)
X Down-valley distance (m)
72
-------
I XC Down-valley distance (km)
XF Distance to final plume rise (m)
XX Intermediate distance variable
Array Name Definitions
None
Subroutines Called
None
Common Blocks
None
Special Note
The user might wish to specify an "effective stack height" directly, rather
than specifying a physical stack height and allowing this subroutine to
calculate the additive plume rise. This can be done by inputting the effec-
tive stack height in place of the physical stack height, and inputting a zero
value for either stack radius or stack gas exit velocity. PRISE returns a
zero value for plume rise when a zero value is input for either of these
parameters.
73
-------
DILUTE-Plume Dilution Module
Purpose
DILUTE calculates the dilution factor necessary to account for clean air
dilution of the nocturnal pollutant plume during its transport down the valley
due to clean air inflow from valley tributaries, slope flows, or entrainment
at the top of the down-valley flow layer.
Inputs
The inputs to DILUTE come from the main program through the subroutine
argument list. They include the cross-sectional areas and wind speeds at the
pollutant source and the receptor cro'ss sections.
Outputs
The output from DILUTE is the dilution factor calculated from the ratio
of volume fluxes at the two cross sections. It is passed to the main program
through the subroutine argument list.
Variable Name Definitions
INPUT/ VARIABLE
OUTPUT NAME DEFINITION
O
I AC Area of valley cross section at distance x (nr)
I AS Area of valley cross section at pollutant source (m2)
0 OF Inverse of dilution factor
I UC Wind speed at cross section at distance x (m/s)
I US Wind speed at pollutant source cross section (m/s)
Array Name Definitions
None
Subroutines Called
None
Common Blocks
None
74
-------
INGRAT-Valley Plume Reflection Module
Purpose
The purpose of Subroutine INGRAT is to calculate Gaussian diffusion
coefficients and to accomplish cross wind integrations of plume mass on a
valley cross section in order to determine the amount of plume mass within the
valley cross section and the amount of plume mass which has escaped out the
top of the valley inversion. These mass calculations are necessary to simu-
late valley plume channeling so that nocturnal plume reflection off the valley
floor and sidewalls can be handled in the main program.
Inputs
Inputs to Subroutine INGRAT come into the subroutine through the sub-
routine argument list. Outputs include the fraction of plume mass within the
valley inversion, the fraction which has diffused out the top of the inver-
sion, and the two dispersion coefficients, sigma y and sigma z.
Variable Name Definitions
INPUT/ VARIABLE
OUTPUT NAME
DEFINITION
I ALPHA1 Sidewall #1 elevation angle (rad)
I ALPHA2 Sidewall #2 elevation angle (rad)
F Analytical formula for Gaussian distribution
I H Plume center!ine height (m)
I HZERO Initial inversion height (m)
I Height increment counter
I ISTABY Horizontal atmospheric stability index (1-6)
I ISTABZ Vertical atmospheric stability index (1-6)
NP Number of height increments to inversion top
P Number of height increments to inversion top
PHIY Area under Gaussian curve from yl to y2
PHIY1 Area under Gaussian curve from minus infinity to yl
PHIY2 Area under Gaussian curve from minus infinity to y2
PI Trigonometric constant
R Height increment counter
0 SIGMAY Standard deviation of plume concentration in y-direction (m)
0 SIGMAZ Standard deviation of plume concentration in z-direction (m)
0 SUM1 Fraction of pollution within valley inversion (0-1)
0 SUM2 Fraction of pollution diffusing out inversion top (0-1)
I W Width of valley floor (m)
I XC Down-valley distance from stack to cross section (km)
I Yl Y-coordinate of sidewall #1 at height z (m)
Y2 Y-coordinate of sidewall #2 at height z (m)
I YZERO Off-centerline displacement of stack (m)
I Vertical coordinate (m)
ZINC Height increment for integration (m)
75
-------
Array Name Definitions
INPUT/
OUTPUT ARRAY NAME DEFINITION
SIGGY(18) McMullen's coefficients for P-G a
SIGGZ(18) McMullen's coefficients for P-G a
z
Subroutines Called
NORMAL
Common Blocks
None
76
-------
NORMAL-Normal or Gaussian Curve Integration Module
Purpose
Subroutine NORMAL is used to calculate the area under a Gaussian or
Normal distribution curve from minus infinity to X standard deviations. The
polynomial approximation technique used in this subroutine comes from
Abramowitz and Stegun [33, p. 932] following a formulation attributed to
Hastings [34].
Inputs
The sole input is the standard deviation X where X is between 0 and
positive infinity. This input comes into the subroutine through the
subroutine argument list.
Outputs
The sole output is the area under the Gaussian curve, PHI. This output
is carried back to the calling program (Subroutine INGRAT) through the
argument list.
Variable Name Definitions
INPUT VARIABLE
OUTPUT NAME DEFINITION
0
-
-
-
I
PHI
PI
SIGH
T
X
Area under curve (0-1)
Trigonometric constant
Intermediate variable
Intermediate variable
Standard deviation
Array Name Definitions
INPUT/ VARIABLE
OUTPUT NAME DEFINITION
C(6) Coefficients of polynomial
Subroutines Called
None
Common Blocks
None
77
-------
GAUSS-Gaussian Plume Module
Purpose
This subroutine uses a Gaussian plume algorithm to calculate pollutant
concentration at an arbitrary receptor location P(x,y,z). The present version
uses Pasquill-Gifford dispersion coefficients, sigma y and sigma z, calculated
using McMullen's [28] method.
Inputs
The subroutine requires information on the coordinates of the receptor
and source, the plume centerline height, source strength, wind velocity, dis-
persion coefficients, atmospheric temperature, and atmospheric pressure. The
inputs come into the subroutine through the subroutine argument list.
Outputs
/
The subroutine output is the pollutant concentration at the receptor in
micrograms per cubic meter. The output is sent to the calling program (MAIN
or PROFIL) through the subroutine argument list.
Variable Name Definitions
INPUT/ VARIABLE
OUTPUT NAME DEFINITION
0
-
I
-
-
I
I
I
I
I
I
I
I
CHI
CHIOQ
OF
Gl
G4
H
Q
SIGMAY
SIGMAZ
US
Y
YZERO
Z
Pollutant concentration (ng/m )
Concentration/source strength (ng s m" kg" )
Inverse of plume dilution factor
Of f-centerl ine term in y direction
Off-centerline term in z direction
Plume centerline height (m)
Source strength (kg/s)
Sigma y (m)
Sigma z (m)
Wind speed (m/s)
Receptor coordinate (m)
Stack offset distance from valley floor centerline
Receptor coordinate (m)
(m)
Array Name Definitions
None
Subroutines Called
None
Common Blocks
None
78
-------
SOLAR-Extraterrestrial Solar Radiation Module
Purpose
The purpose of Subroutine SOLAR is to calculate the time of sunrise,
length of day, and extraterrestrial solar flux on a horizontal surface at
solar noon given the day of year and the latitude and longitude of the site.
The outputs of SOLAR are necessary to drive the daytime portion of VALMET.
Inputs
Necessary subroutine inputs include the latitude and longitude of the
site and the Julian date of the simulation. These three inputs come into
SOLAR through the subroutine argument list.
Outputs
Outputs from Subroutine SOLAR include the local standard time of sunrise,
the length of the daytime period, and the extraterrestrial solar flux on a
horizontal surface at solar noon.
Variable Name Definitions
INPUT/ VARIABLE
OUTPUT NAME
0 Al Solar
DEFINITION
flux on horizontal surface
at solar noon
Conversion factor-degrees to radians (rad/deg)
Cosine of zenith angle
Julian date (1-366)
Date interpolation variable for equation of time
Declination (rad)
Maximum declination (rad)
Julian date of vernal equinox
Eccentricity of earth's orbit
Date interpolation variable for equation of time
Time unit conversion variable
Julian date (1-365)
Latitude (°N)
Longitude correction (h)
Longitude of earth in its orbit around sun (rad)
Time of sunrise (hhmm)
Earth's position in orbit around sun on date of
interest (rad)
OMDZRO Earth's position in orbit around sun at date of
vernal equinox (rad)
OMEGA Revolution rate of earth (rad/day)
ONEHR Conversion factor (rad/h)
PI Trigonometric constant
CONV
COSZ
D
02
DECLIN
DECMAX
DZERO
ECCENT
ID
IT
JULDAY
LAT
LONCOR
LONG
NTSR
OMD
79
-------
RDVCSQ Earth-sun distance factor
SC Solar constant (W/m2)
SR Hour angle of sunrise (rad)
STOLON Longitude of standard meridian of time zone (°W)
Tl Time unit conversion variable
T2 Time unit conversion variable
TAU Length of daylight period (h)
TIMCOR Correction due to equation of time (h)
TMNOON Local standard time of solar noon (h)
URLONG Longitude of site (°W)
Array Name Definitions
INPUT/
OUTPUT ARRAY NAME DEFINITION
EQNTIM(25) Equation of time correction at 15-day intervals
(h)
Subroutines Called
None
Common Blocks
None
80
-------
EBDGT-Surface Energy Budget Module
Purpose
The purpose of EBDGT is to determine the fraction of extraterrestrial
solar heat flux that is converted to sensible heat flux in the valley of
interest. It is this sensible heat flux that drives the post-sunrise inver-
sion breakup that leads to air pollution fumigations on the valley floor and
sidewalls. The fraction depends on factors associated with the surface energy
budget, including surface albedo, soil moisture, surface cover, cloud cover,
atmospheric transmissivity, and long-wave radiative transfer within the val-
ley. The present version of the subroutine is not fully developed and does
not yet explicitly include these factors. The present skeletal version of the
subroutine merely specifies a fraction that the user has previously input
using guidance given in earlier sections of this report. Further development
of this subroutine is a suggested priority for future work.
Inputs
The fraction of extraterrestrial solar flux converted to sensible heat
flux in the valley is the sole input to the present version of EBDGT, and is
input to the subroutine through the subroutine argument list.
Outputs
The present version of EBDGT simply outputs to the main program the user-
specified input value. The output is transferred to the main program through
the subroutine argument list.
Variable Name Definitions
INPUT/ VARIABLE
OUTPUT NAME DEFINITIONS
I AO Fraction of extraterrestrial solar flux converted
to sensible heat flux (0-1)
0 A2 Same as AO (0-1)
Array Name Definitions
None
Subroutines Called
None
Common Blocks
None
81
-------
DESCNT-Inversion Descent and CBL Growth Module
Purpose
DESCNT uses Whiteman and McKee's [3] bulk thermodynamic model of tem-
perature inversion breakup to calculate post-sunrise changes in valley inver-
sion and CBL depths. The method considers that the valley temperature
inversion represents an energy deficit relative to the warmer air above the
inversion. The energy deficit is destroyed after sunrise by solar energy
input into the valley as it is converted to sensible heat flux at valley sur-
faces. The energy is used in two ways: (1) it is used to grow CBLs over
valley surfaces and (2) it is used to cause air to flow upslope in the valley
sidewall CBLs. Air flowing up the sidewalls in the CBLs causes corresponding
descending motions over the valley center, resulting in a descending inversion
top after sunrise. The module uses a numerical method for the calculations,,
and incorporates the effects of warming above the valley in retarding the
inversion breakup.
Inputs
The inputs to DESCNT include solar flux parameters, valley characteris-
tics parameters, model time step, initial CBL and inversion top heights, an
above valley warming rate parameter, inversion potential temperature gradient,
and mean valley air temperature and pressure. These inputs are passed to the
subroutine through the subroutine argument list.
Outputs
The main outputs of DESCNT are CBL heights and inversion top heights as a
function of time (actually, time step). These outputs are passed back to the
main program in arrays located in common block BLK1. The number of time steps
required to destroy the valley inversion is also an output of the subroutine,
passed to the main program through the subroutine argument list.
Variable Name Definitions
INPUT/ VARIABLE
OUTPUT NAME DEFINITIONS
I
I
I
I
I
_
-
I
-
-
-
AO
Al
ALPHA1
ALPHA2
BETA
C
CP
DELTAT
DHC
DHDT
DHT
Fraction of Al converted to sensible heat flux (0-1)
Solar flux on horizontal surface at solar noon (Wm )
Sidewall #1 elevation angle (rad)
Sidewall #2 elevation angle (rad)
Rate of temperature change at inversion top (°K/s)
Intermediate factor ' . ,
Specific heat at constant pressure (J kg" K~ )
Time step (s)
Change in CBL height (m)
Change of CBL height with time (m/s)
Change in inversion top height (m)
82
-------
DHTDT Change in inversion top height with time (m/s)
FACT1 Multiplicative factor
FACT2 Multiplicative factor
FDENOM Intermediate variable-denominator
FNUM Intermediate variable-numerator
I GAMMA Vertical potential temperature gradient (°K/m)
HCBL Height of CBL top (m)
I HCBLI Initial CBL height (m)
HTOP Height of inversion top (m)
I HZERO Initial height of inversion top (m)
I K Fraction of sensible heat flux going into CBL growth
(0-1)
N Time step counter
0 NSTEPS Number of time steps required to destroy inversion
PI Trigonometric constant
I PRESS Average pressure at inversion center at sunrise (mb)
I RHO Air density (kg/nr)
T Elapsed time since sunrise (s)
I TAD Length of daylight period (s)
THTOT Potential temperature/ambient temperature
TI Time of sunrise (s)
I W Valley floor width (m)
Array Name Definitions
INPUT/
OUTPUT ARRAY NAME DEFINITIONS
0 HC(NS) Array of CBL heights (m)
0 HT(NS) Array of inversion top heights (m)
Subroutines Called
None
Common Blocks
BLK1
83
-------
PROFIL-Concentration Profile Module
Purpose
The purpose of Subroutine PROFIL is to calculate the pollutant concentra-
tion injected into the top of each of the model grid elements at each time
step.
Inputs
Subroutine PROFIL requires a number of inputs, dealing primarily with the
valley characteristics, grid element locations and specifications, and
Gaussian plume characteristics. Subroutine PROFIL is called for each model
time step and does calculations for each of the grid elements. All inputs
come into the subroutine through the subroutine argument list.
Outputs
For each model time step, Subroutine PROFIL calculates the pollutant
concentration injected into the top of each of the growing grid elements. To
maintain accuracy in the calculations, the concentration for each grid element
is calculated as the average of four concentration determinations made at the
various sides of the quadrilateral Gaussian plume element which sinks into the
grid element.
Variable Name Definitions
INPUT/ VARIABLE
OUTPUT NAME
DEFINITIONS
_
I
I
I
I
mf
~
-
-
I
I
I
I
I
I
I
-
I
I
ALPHA
ALPHA1
ALPHA2
BOXLEN
CHIOFF
CHIY1
CHIY2
CHIZ1
CHIZ2
DF
H
HCL
HCU
HTL
HTU
HZERO
I
NBOX
NTOT
Average of ALPHA1 and ALPHA2
Sidewall #1 elevation angle (rad)
Sidewall #2 elevation angle (rad)
Grid element length (m)
Offset air pollution concentration
Concentration at P(X,Y1,Z) (ng/m;:
Concentration at P(X,Y2,Z) (p-g/m^
Concentration at P(X,Y,Z1) (ng/m;r
Concentration at P(X,Y,Z2) (M-g/m'3
Inverse of plume dilution factor
Plume centerline height (m)
CBL height at lower time step (m)
CBL height at upper time step (m)
Inversion top height at lower time
Inversion top height at upper time
Inversion depth at sunrise (m)
Grid element index
Number of grid elements on valley
Number of model grid elements at a
due to reflections
}
)
)
)
step (m)
step (m)
floor half-width
given time step
84
-------
I
I
I
I
I
I
I
-
-
-
I
-
-
NTOTI
PRESS
Q
SIGMAY
SIGMAZ
TEMP
US
Y
Yl
Y2
YZERO
Z
Zl
Z2
Initial number of model grid elements
Ambient pressure (mb)
Source strength (kg/s)
Sigma y (m)
Sigma z (m)
Ambient temperature (°K)
Down-valley wind speed at stack cross section (m/s)
Y-coordinate of receptor (m)
Y-coordinate at left side of sinking mass element (m)
Y-coordinate at right side of sinking mass element (m)
Y-distance of stack from valley floor centerline (m)
Z-coordinate of receptor (m)
Effective z-coordinate at bottom of sinking mass
element (m)
Effective z-coordinate at top of sinking mass element
(m)
Array Name Definitions
INPUT/
OUTPUT ARRAY NAME
DEFINITIONS
0 CHIBAR(NB) Average concentration injected into top of each
grid element
Subroutines Called
GAUSS
Common Blocks
BLK2
85
-------
VELOCY-Upslope Flow Velocity Module
Purpose
Subroutine VELOCY uses the mass continuity equation for air to calculate
the slope flow velocity profile, under the assumption that the convective
boundary layer depth does not vary as a function of distance up the slope. An
implicit assumption is that all mass lost from the valley inversion cross sec-
tion as the inversion top sinks is transported from the valley in the upslope
flows which develop in the convective boundary layers over the slopes. The
upslope velocities are calculated at the boundaries between grid elements.
Inputs
The inputs to this subroutine come through the subroutine argument list
and through the common block labeled BLK1. BLK1 provides the inversion top
and CBL heights for every time step. Inputs from the argument list include
the time step index number, the time step length, and information about the
size and location of the grid elements.
Output
The output of VELOCY is an array of upslope wind velocities calculated at
the NTOTI+1 boundaries between the NTOT grid elements. The subroutine is
called from the main program at each model time step, and the calculated
velocities are passed back to the main program through the common block
labeled BLK3.
Variable Name Definitions
INPUT/ VARIABLE
OUTPUT NAME
DEFINITIONS
I BOXLEN Length of grid elements (m)
CIH Inversion top displacement during time step (m)
I OELTAT Time step (s)
I Grid element index
I MN N-l
I N Time step index
I NBOX Number of grid elements on valley floor half-width
I NTOT Total number of grid elements in simulation at a
given time step
I NTOTI Number of grid elements in simulation at beginning
of simulation
86
-------
Array Name Definitions
INPUT/
OUTPUT ARRAY NAME DEFINITIONS
I
I
0
HC(NS)
HT(NS)
V(NB1)
CBL height array (m)
Inversion top height array (m)
Upslope wind velocity array (m/s)
Subroutines Called
None
Common Blocks
BLK1
BLK3
87
-------
BRKUP-Pollutant Mass Budget Module
Purpose
The purpose of Subroutine BRKUP is to calculate pollutant concentrations
in each of the grid elements using a pollutant mass balance. The calculations
are made for each model time step, taking account of pollution sinking into
the top of each grid element, CBL growth, upslope advection of pollution,
mixing within each grid element and carryover of pollution in each grid ele-
ment from the previous time step. A currently unused feature of the pollution
mass budget in each grid element is the provision for sources and sinks of
pollution within the element. A special feature of BRKUP is the exponential
decay of air pollution concentrations in grid elements that are dropped from
the simulation as the top of the inversion sinks below the grid element.
Inputs
Inputs to the subroutine coming through the subroutine argument list
include time step index, number of grid elements on the valley floor half-
width, total number of grid elements at model initiation, the grid element
length, and model time step length. Other inputs are available through
labeled common blocks BLK1, BLK2, and BLK3. BLK1 provides information on CBL
and inversion depths for each time step. BLK2 provides information generated
in Subroutine PROFIL on the amount of air pollution mass sinking into the top
of each grid element at the given time step. BLK3 provides the upslope wind
velocities used to advect pollutant mass up the valley sidewalls at the given
time step. These velocities were generated in Subroutine VELOCY.
Outputs
BRKUP calculates air pollution concentrations in the model grid elements
for each time step. These calculations constitute the main output of the
VALMET model, and are output to the main program through the subroutine argu-
ment list. The main program performs further calculations on these concentra-
tions to determine 5-min-average concentrations which are stored in the main
program for further processing into maximum 1- and 3-h averages.
Variable Name Definitions
INPUT/ VARIABLE
OUTPUT NAME DEFINITIONS
A Pollutant mass within the grid element at previous
time step (ng)
B Pollutant mass coming from adjacent downhill grid
element (|ig)
I BOXLEN Grid element length (m)
C Pollutant mass sinking through growing top of CBL
(ng)
D Pollutant mass transported into adjacent uphill grid
element (,ug)
-------
DELTAH Entrainment at top of grid element (m)
I DELTA! Time step (s)
I Grid element index
I MN Index of previous time step
I N Current time step index
I NBOX Number of grid elements on valley half-width
I NTOTI Total number of grid elements at model initiation
S Exponential decay time constant (s"1)
Array Name Definitions
INPUT/
OUTPUT ARRAY NAME DEFINITIONS
0 CCH(2,N8) Grid element concentration array (ng/m)
CHIBAR(NB) Average pollutant concentration injected into
top of grid element (ug/m )
HC(MS) CBL height array (m)
HT(NS) Inversion top height array (m)
V(NB1) Upslope velocity array (m/s)
Subroutines Called
None
Common Blocks
BLK1
BLK2
BLK3
89
-------
PSTPRC-Post Processor Module
Purpose
This subroutine is a post-processor. It processes the time series of
calculated 5-min-average concentrations for each grid element to obtain
maximum short-term (1- and 3-hr) average concentrations which can be compared
to regulatory standards. The maximum 1- and 3-hr-average concentrations are
determined by means of a moving average of the 5-min-average concentration
time series, determining maximum values for each of the grid elements. A time
index value is saved for each of the maximum values so that the time of
occurrence of the maximum can be calculated in the main program. The
subroutine is written so that the averaging intervals can be easily changed,
if necessary.
Procedure
The moving average is calculated by stepping backwards through the
5-min-average concentration array using a 1- and 3-hr averaging "window",
beginning at the last 5-min average calculated for each grid element. As they
travel backwards through the 5-min concentration array the 1- and 3-hr
averaging intervals will eventually pass the time of sunrise. The nocturnal
steady state concentrations are added into the average at that point.
Inputs
The inputs to PSTPRC include the series of 5-min-average concentrations
for each of the grid elements, the steady-state nocturnal concentrations in
each grid element, the averaging interval (300 s), total number of grid ele-
ments and index number of the last 5-min concentration values. The time
series of 5-min averaged values is input through common block BLK2. The
nocturnal concentrations are input through common block BLK. The other
parameters are input through the subroutine argument list.
Outputs
The outputs of Subroutine PSTPRC are passed to the main program through
common block BLK6. They include the maximum 1- and 3-hr average concentra-
tions for each grid element, the time index numbers for each of these
concentrations, and the number of 5-min time steps in the regulatory averaging
intervals (1- and 3-hr).
Variable Name Definitions
INPUT/ VARIABLE
OUTPUT NAME DEFINITIONS
I Index for grid elements
K Index for averaging intervals
=1 for 1-hr averages
=2 for 3-hr averages
90
-------
N Index for 5-min-average time series
MAS Number of elements in 5-min series
NSTART NAS+1
NTIMES Time index intermediate variable
NTOTI Total number of grid elements in model at sunrise
ONEHR Number of seconds in 1 hr (s)
Ql Intermediate summation variable
SUM Summation variable for concentrations
THREE Number of seconds in 3 hr (s)
TS Number of 5-min-averages in averaging interval
TSTP Averaging interval of basic model concentration
array =5 min= 300 s
Array Name Definitions
INPUT/
OUTPUT ARRAY'NAME DEFINITIONS
0 AVG(NB,2) Maximum 1- and 3-hr-averages for each grid
element (|j.g/m )
I CHI(MB) Nocturnal concentrations in grid elements
Ug/rr)
I CONC(NA,NB) Array,,of 5-min-average concentrations
(ug/m-3)
0 NDXTIM(NB,2) Array of time index values for maximum 1- and
3-hr averages
0 NTS(2) Number of 5-min intervals in 1 and 3 hr
SUMM(NB,2) Summing variable for 5-min-average
concentration
Subroutines Called
None
Common Blocks
BLK4
BLK5
BLK6
91
-------
SECTION 6
SAMPLE MODEL RUNS
In this section we present sample model runs to illustrate model perfor-
mance and model outputs. In the first model run we present the model results
for a run using default input parameters. The first-time model user may
compare his results with the results from this run to see that the model is
operating properly at his installation. In the second sample run we make
calculations for a downwind distance of 30 km rather than the default 10 km,
but maintain all other default values for model inputs.
SIMULATION NUMBER 1
The first simulation to be illustrated was obtained using default input
values obtained from Table 5. A full listing of the summary output file
(VALMET.OUT) generated from this simulation is given below in Figure 21. This
output file would be printed by the model user at the completion of the model
run and consists of three pages of printer output. Outputs shown in the
figure include:
1. a model input table,
2. a listing of numerical simulation parameters,
3. a summary listing of nocturnal model output values,
4. a listing of the output values from the solar model,
5. an output table of maximum 1- and 3-hr-average concentrations in each
grid element, and
6. a table of 5-min-average (300 s) concentrations as a function of time for
each model grid element.
The user who wishes to get a graphical depiction of model results may manually
plot the model results from the output listing. Alternatively, he may develop
a plotting program using as input the special file of model output data,
VALMET.PLT, that is generated with each model run. Figures 22 and 23
illustrate the types of plots that may be generated using such a plotting
program. A listing of the plotting program that generated these figures is
given in Appendix B.
92
-------
I
»3
U
S
g
H
u a a u
a
i~t m rjoc4oo
o • in o fe
3 O *n
.
i. t -3 < U W < S
_; < =; - a: « e- i- '/j
< 'J— pH&i>i05U5<
•* —« ^ ** p* EN r«
^ .en > O ff) O • S Q O
CM -3 S> 3 OO^OO
. c :> • • o o • •
o« in o ooooo
t .^ r* ^ fM C4 *N
<
S
z
H
X
0*
=
X
3
Eu
S
U
U
2
yj
z
v
'•&
a
u
ARACT
=
U
>*
u
_:
1
u
-•§
s§
2s
as
< z
S2
yj
j s
w u
a >
i5
z
a
S
<
8
>J
U
u
e
u
X U
a. s
co a
i*
^S
U U]
Id- U
s- C-* « —
w « <
« 3 > a
CB
RACTE
<
u
a«
U
<
JH
U
a
'ARAHE'J'tH
^
<
O
at
i
r*
an
W
fH
V
1
5
1
in
3
Z
O
i
»•»
as
u
-£
r*i <
-« >•
S
zg
Sj
§S ~
«z in
6-1 U
M z
S2 1
z282
-« WO «
W £d O -3
Sous
a tax
Z U O «!
M z O
1-1 m u
l*;.g
M u o a z
u to **<'•«
C SS S •« OS W
Off* 3 36- <
CQ fS U CC £ Z >
a c -H ^ o:
z u a to u
o too o .« fr«
U Z O fa f» IM H.
ua w o ^
3C *-* Cfi CQ < :?
•H < M W < <
a u u _ a
U 7-
U ^ -3 H
,&* as w ae ts, z cs
03 CJ £d U U 03
Siilill
an
j§
S8
U2
uu
M a; ^
§gj
as> to
<~a3
sSg£
o3&* Z
irt < O O *
M > -* «
a* 03 tii
Mp a: e-
2PSSS
3 M Z
asp-si
^g^-ss
a- o c-o •
SsS*«*
u O C 3 rd IN
£ U U Q & ^
ZJJ =
23£B2
ca H &< a < w
3Z« K BW
uaava
u &* o w
an ens
S2S5yU
gssie*
< = = u.^o
x u ^ « a «
u Q* a* 2 •
C 3 O fci
X U U O
u ssa
-i £-• &» cd &• ua
*i.u.5°i
OOO X Z <
!N O W
M . of •• Cfl — 3
d-^-i-Si
> in o W a e-t
< ot ai &• o
K A* < S Z O
r- oir 2 u^o
ni M U
>j z u yi z *-*
< o o o * q o
K &4.9Q U .-
u f* cd o w
O S £«. i* Z fa. Z
;J O tO pt
lsiz°i§£
r? tu^t t» ec ia z z
Z -4 Q Q £ 5" c-
o a o < —
3 z > > •-» £ a
Ci. U < < Z p —
_ ji§ 5 ^ w S
a ca z z c* <
I ca M a o <
M O Z 3 CS M —
l
ca
S: M S W
< S — -3 I 1 1 t
3£u~33SS
SxstU ^^
Ji-S
H GO ooor-
u3£'~
U £ ( I t l
u 3 H
1 §*
a u a»
3£5S~^,.
1 1 1 1 1 1 I 1 ( l 1 1 1 I i l
00030090000<=0<303
t 1 , 1 1 ' 1 1 1 1 1 I ' I L 1
JgS*"^£'5'°S^5'n2a"n
t i i i i j i i i i ! : , i i
r--»o^f-r--.-,«^-.c=^
0)
01
c
01
93
-------
i-*^nr4(N»«M'v>.3p*r-aio>ff*a>oor-^ir4f«N-^ooooQoaaooocjooooooooooi3OOQOooooooooo
•-» CQ 0000000000000000000000000000000000000000000000000020000=90 o
III * °IIIII°ir°°IIo!IIIIIIoIII°°II°III^^ *
£O CH^r*fflflOCDflOflQWOOWff»W
-------
600
o
200
b)
\
\
0 1 2
CONC (ug/m3)
789
TIME (LSI)
10
11
n
E
o>
O
z
o
o
c)
0-4-3-2-10 1 2 3 4
CROSS-VAUH DISTANCE (km)
Figure 22. Plots of Sample Simulation 1: (a) nocturnal vertical concentra-
tion profile through plume centerline, (b) CBL height (I°n9
dashes) and inversion height (short dashes) as a function of
time, and (c) nocturnal cross-valley concentration profile
through plume centerline.
Figure 22 illustrates some of the characteristics of the simulation. It
presents plots of the nocturnal vertical air pollution concentration profile
at the valley center (Figure 22a), the nocturnal cross-valley concentration
profile through the plume centerline (Figure 22c), and plots of CBL height
(long dashes) and inversion top height (short dashes) as a function of time
(Figure 22b).
Figure 23 presents plots of the 5-min-average pollution concentrations
for selected model grid elements. Twenty grid elements are necessary for the
95
-------
1.0
0.8
X= 10km
CL CONC= 1.20jjg/m3
% OUT TOP= 0.1 %
E
^
O
o
0.4
0.2
0.0
8
10
11
TIME (1ST)
Figure 23. Pollutant conentration versus time for selected grid elements for
Sample Simulation 1. The down-valley distance, centerline con-
centration and fraction of mass diffusing out the top of the
valley inversion during the plume's nocturnal transport are
indicated on the figure.
default simulation—three on the valley floor and 17 on the valley sidewall.
Plots are shown for grid elements 3, 6, 9, 12, and 15. Sunrise on the date of
the simulation was at 0550 am LSI. Low steady-state concentrations are seen
in the figure before sunrise, with time-varying concentrations arising from
fumigations and upslope advection after sunrise. Nocturnal concentrations on
the valley floor and sidewalls are low since the plume, traveling down the
valley centerline at 250 m height, has diffused insufficiently both vertically
and horizontally to produce high concentrations on the valley surfaces. The
predominant feature of the concentration curves at the 10 km travel distance
is the post-sunrise fumigation of the pollutant plume on the valley floor and
sidewalls. The highest concentrations occur on the valley floor as the high
concentrations at plume centerline subside into the growing CBL. The highest
concentration shown occurred in grid element number 3 and is 0.77 u.g/m . This
may be compared to the nocturnal plume centerline concentration of 1.20 ng/m .
Grid elements on the valley sidewalls are strongly affected by upslope advec-
tion. For example, the concentration in grid element 15 increases rapidly
just following sunrise as the higher concentrations in downslope grid elements
96
-------
(e.g., see grid element 12) are advected up the sidewall. Another feature of
the simulations, apparent in the figure, is the exponential decay of concen-
trations in individual grid elements as the temperature inversion top subsides
below their elevation. This occurs earliest at the highest grid elements.
The exponential decay in concentrations simulates the effect of a sudden
increase in mixing depth that occurs when the temperature inversion descends
below the grid element. The concentration in grid element 3 begins to drop
rapidly to zero at the time of temperature inversion destruction (0933 1ST).
SIMULATION NUMBER 2
The results of sample simulation number 2 are shown below in Figures 24
and 25, in the same way that the results of simulation number 1 were illus-
trated. Input parameters for this simulation, 30-km down-valley from the
pollutant source, differ from the previous simulation only in down-valley
distance.
At the 30 km cross section, the nocturnal plume had diffused sufficiently
during travel that substantial reflections had occurred from the valley sur-
faces. The vertical and horizontal concentration profiles through the plume
center!ine (Figure 24) show that the concentrations were more uniform within
the valley than at the 10-km section. At the 30-km section, 2.7% of the plume
had diffused out the top of.,the valley inversion, and the plume centerline
concentration was 0.39 ng/m . The CBL growth, inversion top descent, and
inversion breakup occur exactly as shown in the first simulation. Concentra-
tions (Figure 25) differ during both the pre-sunrise and post-sunrise periods,
however, due to the different nocturnal concentration profiles within the
cross section at the longer travel distance. Nocturnal concentrations are
higher at the valley surfaces due to the broader distribution of the plume
about its centerline. The highest concentration in sidewall grid element 12
occurs during the night. Since the pollutant is more evenly distributed
across the cross section, the effects of fumigations are less pronounced than
at the 10-km cross section. Again, as at the 10-km cross section, the
greatest effect of fumigation is on the valley floor.
The example shown here as simulation number 2 is used merely to illu-
strate the effect of changing down-valley distance, other parameters remaining
constant. Other parameters could have been changed along with the down-valley
distance (such as valley width, for example) to simulate the physical charac-
teristics of a real valley more accurately.
97
-------
0 400
<
£
§200
X
0
a)
\
r\
\
f 1
i.i.
600
/•^
< ***
J,
O
| 200
b)
h
"
h
L 3
—••**"*'
1 i 1 1 1
0 1 2
CONC (ug/m3)
789
T1ME(USIO
10
11
/*"V A
lO
E
3
^ 1
0
z
O
0
0
•
. _ _ a
_^ — »-^""— "~~"'^™*>—- ^.
"T^ T . .
0-4-3-2-10 1 2
CROSS-VALLEY DISTANCE (km)
Figure 24. Same as Figure 22, for Sample Simulation 2.
98
-------
,-N
fO
E
1?
V— X
o
z
o
o
1.0
0.8
0.6
0.4
0 7
u.z.
/N A\
. X=30km
. CL CONC= 0.39jjg/m3
. % OUT TOP= 2.7%
i
i
i
^1 -^^=3«~, in — . —a -»*>3f*^-. __ "*^^>^ ft"I
"" ' ^""Z"*^^? \ ^ \ \
"~ — — •" **^^«-^^*^' \ « ^nQ*rtc
- — \15 \'2 \ , "\
\ \ \% \\
1 , L , i ^-^v X>— 1 ^ — ^--S^a—
7 8
TIME (1ST)
10
11
Figure 25. Same as Figure 23, for Sample Simulation 2.
99
-------
SECTION 7
FURTHER WORK
In this section, we will discuss the need for the further development
and testing of the VALMET model. VALMET, as described in this report, is a
working model able to simulate many of the observed features of the
meteorology of Colorado's deep valleys. The model has been developed in a
modular fashion so that the modules can be upgraded as we learn more about
valley meteorology. The model, while incorporating what are thought to be
the relevant physical processes that lead to air pollution transport and
dispersion, has not yet been fully tested against actual air quality data.
Nevertheless, some initial results [52,53,54] are now available from
meteorological and tracer experiments conducted in the Brush Creek Valley of
Colorado in 1982. The tracer experiments, supported by the EPA in coopera-
tion with the U.S. Department of Energy's (DOE) Atmospheric Studies in
Complex Terrain (ASCOT) Program, were designed to provide the data necessary
to perform an initial evaluation of the VALMET model. The evaluation is not
yet complete, but is sufficient to provide us with initial guidance on the
further development of VALMET. In addition to the 1982 tracer experiments,
a more comprehensive set of tracer and meteorological experiments was
conducted in the same valley in September and October of 1984 as part of
DOE's ASCOT program. The 1984 ASCOT data are not yet available, but they
will constitute an important resource to further test the model. A portion
of the tracer experiment, supported by the EPA, was designed specifically to
test the VALMET and MELSAR [55] air quality models developed in its GRAMA
program. In the following, we will briefly summarize the results from the
1982 experiments, pointing out the implications for future modifications to
the VALMET model. Suggested modifications resulting from these experiments
and from other sources will then be individually discussed.
GUIDANCE FROM THE 1982 EXPERIMENTS
The 1982 tracer experiments were designed to provide the data required
to evaluate the initial version of VALMET [52]. We did not consider it
sufficient to simply collect tracer concentration data on a cross-valley arc
and compare this with model calculations. Rather, the approach taken was to
collect meteorological and tracer data to test the full range of
meteorological assumptions and parameterizations used in modules within the
model. For example, the model predicts that convective boundary layers will
grow over heated surfaces after sunrise, that upslope flows will develop
within these boundary layers, that pollutants from the elevated nocturnal
plume will fumigate into the convective boundary layers, and that they will
be transported out of the valley by the upslope flows. Thus, within the
100
-------
restraints of the resources available, it was necessary to observe the
development of convective boundary layers over the slopes, the upslope wind
systems, fumigation of pollutants, and transport of pollutants up the slope.
In addition, it was necessary to simulate an elevated release of pollutants
and to observe the characteristics of the nocturnal plume.
The EPA tracer experiments were conducted in a valley chosen by DOE
using criteria unrelated to the testing of the VALMET model. The Brush
Creek Valley was a useful "target of opportunity" for the initial evaluation
of VALMET, but, as is usual with such opportunities, there were advantages
and disadvantages to the choice of this particular valley.
There were several advantages to choosing the Brush Creek Valley for
the initial evaluation of VALMET. First, the valley has a rather simple
topography. The narrow, 25-km-long valley has no major changes in valley
orientation along its length. It has nearly equal sidewall inclinations.
The valley drains a plateau, so that the ridges are at a constant altitude
regardless of location along the valley axis. The valley has no major
tributaries. Second, the valley axis is oriented from NW to SE so that the
sidewalls would be exposed to quite different insolation during the post-
sunrise temperature inversion breakup period. The effect of this unequal
heating was a major uncertainty in the model formulation. On the basis of
meteorological data collected in wider Colorado valleys, and numerical model
results, the VALMET model was developed under an assumption of horizontal
homogeneity of atmospheric structure on a valley cross section. This
assumption could be readily tested in the Brush Creek Valley, where the
narrowness of the valley and the NW-SE orientation of the valley would
clearly maximize any horizontal gradients in atmospheric structure between
the sidewalls. Third, the Brush Creek Valley would be heavily instrumented
with meteorological sensors by the ASCOT program. Access to their
meteorological data would be a great benefit to the model evaluation effort.
Along with the above advantages, there was a major disadvantage to
conducting an initial evaluation of VALMET using data from the Brush Creek
Valley. This disadvantage was related to the short segment of the valley
that was accessible for tracer instrumentation. VALMET is a two-dimensional
model, predicting concentrations on a cross section oriented perpendicular
to the valley axis some distance down-valley from a source. Restrictive
assumptions are present in VALMET regarding a required homogeneity of the
temperature and wind structure in the along-valley direction. The Brush
Creek Valley, however, is a short tributary valley that flows into the Roan
Valley a few kilometers below the valley cross section where most measure-
ments were made. Consequently, tracer plume carried down the Brush Creek
Valley during the night would be carried into Roan Creek. Reversal of the
down-valley winds (to up-valley) after sunrise would result in a large part
of the tracer plume being carried up the Roan Creek Valley, rather than
being carried back up the Brush Creek Valley as assumed in the model.
Evaluation of VALMET would be complicated by this violation of a major
assumption in the model, which had been designed for longer valleys.
101
-------
The full evaluation of the VALMET model will be the subject of future
work. It is appropriate here, however, to make some initial qualitative
statements concerning the evaluation of the model. First, with respect to
the nocturnal portion of the model, the nocturnal plume was carried down the
valley, as expected. The nocturnal plume, although released above the
valley center, was found to be displaced somewhat towards one sidewall as it
was transported down the valley. No completely satisfactory physical
explanation for this phenomenon has yet been offered, and no mechanism is
available in the VALMET model to model it. The valley is not strictly
linear, but turns slightly with down-valley distance. Because the plume was
displaced towards the "outside" of the turn, it is conceivable that inertia!
effects are responsible for the displacement of the plume from the valley
centerline. Further information on this phenomenon will be available from
the 1984 experiments.
The nocturnal plume was carried down the valley in a rather strong
"jet" of down-valley winds, with the level of maximum winds at about release
height. The nocturnal model, based on the Gaussian formulation, is
incapable of treating vertical shears in transport winds but, when winds at
release height are used for transport, the model approximates transport and
diffusion along the valley direction fairly well.
Assumptions in the daytime portion of the model were verified with
actual meteorological and tracer data [53,54]. The post-sunrise period was
characterized by the growth of convective boundary layers over the sunlit
valley surfaces. The tracer plume fumigated the valley sidewalls as
convective boundary layers grew upwards into the remnants of the nocturnal
temperature inversion containing the elevated tracer plume. Tracer was
carried from the valley by upslope flows, which developed within the growing
convective boundary layers. Corresponding subsiding motions over the valley
center were noted in the temperature profiles at several of the tethered
balloon sites, but the limited vertical resolution of the tracer plume did
not allow this feature to be seen in the tracer concentration analyses.
Due to the northwest-southeast orientation of the deep, steep-walled
valley, very significant differences occurred in the timing and rates of
convective boundary layer growth on the opposing sidewalls following sun-
rise. As a result of the unequal heating of the different sidewalls, a
cross-valley flow developed, carrying the elevated plume towards the warmer
sidewall [53]. Due to the cross valley advection, tracer concentrations
were higher on this sidewall than predicted by the model. A future
modification of the VALMET model will be required to handle this situation
properly in narrow valleys where post-sunrise insolation on the opposing
sidewalls is quite different. The Brush Creek tracer experiments were the
first direct experimental confirmation of the importance of this physical
effect on tracer plume dispersion.
The short length of the Brush Creek Valley, as expected, affected the
results of the tracer experiments. The primary effect, from initial
analyses, seems to be that the tracer concentrations in the valley fell more
rapidly than expected after the post-sunrise wind reversal. This is thought
102
-------
to be due to the nocturnal plume being carried largely up Roan Creek after
the wind reversal rather than reversing direction to come back up Brush
Creek.
SUGGESTED MODIFICATIONS TO THE VALMET MODEL
The basic modeling approach in VALMET is to explicitly parameterize a
number of physical processes acting to disperse pollutants in the valley
atmosphere. Based on the parameter!zations now included in the model and
initial results from the 1982 tracer experiments, a number of suggestions
can be made for model improvement. Some of these suggestions are presented
below. The simpler modifications are listed first.
Deposition
Deposition of pollutants during transport could be incorporated into
the model at two points. Deposition during along-valley transport of the
plume could be included in modules INGRAT and NORMAL, while deposition
occurring in the slope flows could be rather easily incorporated into module
BRKUP using a surface sink term in the pollutant mass budget calculations.
Emission Above or Below Stable Core
The present model assumes implicitly that pollutants are emitted into
the stable core, which is then advected down-valley. Concentrations are
produced by fumigation of pollutants from this elevated plume as th*e CBL
grows upward into it. Modifications to the model will be required if the
emissions fail to enter the stable core. This will happen at the stack
location (x=0) once the CBL grows to the effective stack height. Pollutants
would then be emitted into the CBL. Similarly, once the inversion top sinks
below the effective stack height, the pollutants are dispersed above the
valley inversion and do not affect the stable core.
Calculations of air pollution concentrations down-valley of the stack
should take account of these possibilities by modifying pollutant
calculations when changes in the stable core concentration profiles are
advected over the cross section of interest. The advection time depends on
the distance x down-valley of the stack and the wind speed in the stable
core, such that the advection time is
tx - t0 - x/us
where t~ is the time of sunrise and t, is the time when the inversion top
descends below the effective stack height or the CBL grows above the
effective stack height. When calculations are made far downwind of the
stack (x»0), or down-valley flows are weak (ug small) the travel time will
be so long that the inversion will be broken before any change in
concentration calculations is necessary. It may be necessary, however, to
modify the model calculations in short, well-drained valleys. Similarly,
the model, when applied to ground level sources, will need modification to
account for the post-sunrise emission directly into the CBL.
103
-------
Energy Budget
The energy budget model needs a great deal of further work. It
presently includes a simple bulk parameterization of the fraction of
extraterrestrial solar radiation converted to sensible heat flux at the
valley surfaces. This fraction, in general, should vary throughout the day
and should depend on cloudiness, soil moisture, albedo, ground heat flux,
vegetation, atmospheric absorption of the solar beam, etc. It may be
worthwhile to adapt existing surface energy budget parameterization schemes
for use in the energy budget module. This will, of course, increase the
input requirements of the overall model. A different approach is
possible. One could look at existing data sets and try to determine the
value of this fraction leading to the "worst-case" pollution episodes in
different parts of the country. The model could then be used with these
fractions to simulate such conditions. Such an empirical approach should
also be focused on determining the other free parameter in the model, the
fraction k of sensible heat flux causing CBLs to grow. A major effort is
required to do a good job on either approach.
Cross Valley Flows
The model should be modified to incorporate cross valley flows in the
stable core. Such flows have been postulated by others [45] and were
clearly observed in the 1982 Brush Creek tracer experiments. They are
expected to be strongest in valleys where differential heating occurs on the
sidewalls in the morning. An example would be a north-south valley where
the morning sun would shine on one sidewall while the other was still in the
shade. The existence of such flows may strongly affect post-sunrise concen-
trations on the sidewall since the plume, containing high concentrations, is
advected towards the most-strongly sunlit sidewall. Several methods seem
possible for incorporating the modification into VALMET, but little
observational information is available to test the different methods.
Guidance may be available from the more-sophisticated primitive equation
models [46-49].
Turbulent Erosion of the Valley Inversion by Overlying Flows
Too little is presently known about postulated turbulent erosion
[15,50] of the top of a valley temperature inversion to include it in the
present version of the model. It is not clear that it should be
incorporated anyway, since the model is more conservative without it.
Lenschow et al. [50] have developed an approach that could be rather easily
incorporated into the present model if a means could be found to relate a
local Richardson number at the top of the inversion to the rate of turbulent
erosion. No way is yet available to do this. Additionally, recent
primitive equations modeling indicates that turbulent erosion is not a major
factor affecting temperature inversion breakup under conditions of light to
moderate winds aloft when inversions are strong.^ ' Fluid modeling
experiments may be useful here in the future.
(a) Personal communication with T. B. McKee, Colorado State University,
Fort Collins, Colorado, 1983.
104
-------
Effect of Tributary Flows on the Enhancement of Diffusion
The present version of the model incorporates a simple means of
handling dilution due to valley tributary flows. The model does not treat
the effect of these tributary flows on enhancing diffusion within the main
valley. At present there seems to be too few observational data to
incorporate this effect or the related thermodynamic effects caused by
converging air masses of different temperature or stratification. The DOE
ASCOT program may look into these effects in their future Colorado field
programs.
Piffusion Coefficients
The diffusion coefficients in the model need improvement. At present
we are modifying Pasquill-Gifford [43] flat terrain coefficients on the
basis of the Huntington canyon data of Start et al. [20] and other data
summarized by others [29-31]. We need to look at more data or conduct
specific diffusion experiments to determine the most appropriate values of
these coefficients. Alternatives to the Gaussian plume equation should also
be seriously considered.
Time-Varying Wind Speeds in the Stable Core
The model could be modified to handle time-varying wind speeds in the
stable core. There is no doubt that the down-valley wind speeds in the
stable core decrease after sunrise during the temperature inversion
destruction period, since this has been frequently observed [9]. It is not
entirely clear what the effect of this decrease will be on calculations,
since the wind speed, in addition to modifying plume concentration
calculations, will also affect plume rise calculations.
Temperature Inversion Buildup
The present version of the model simulates the nocturnal steady-state
period and the inversion breakup or morning transition period. The model
could be extended to simulate the temperature inversion formation and growth
period, when down-valley winds first become established within a valley. A
large quantity of data on this period was collected in Whiteman's National
Science Foundation-supported dissertation program in 1977 and 1978. The
data have been processed, but not analyzed. A case study [6] indicated some
promising research approaches that could lead to an air pollution model for
the inversion buildup period.
Differential Heating
One might expect CBLs to grow at different rates over the different
sidewalls of a valley, especially when the energy budgets of the sidewalls
are substantially different. This might occur frequently in north-south
oriented valleys where the differential insolation on the sidewalls varies
strongly with time (see Figure 26 for an example). The present version of
VALMET has as its basis a bulk energy budget model of the valley inversion
that does not treat the sidewall energy budgets separately. A major
105
-------
modification to the model would be required to simulate CBLs growing at
different rates on the individual sidewalls. This would require an addi-
tional free parameter in the model, specifying the fraction of sensible heat.
used for growing CBLs over one sidewall, as opposed to the other. At pre-
sent there seem to be no vertical structure data available taken on a single
day on two opposite sidewalls. Without this kind of data it is difficult to
come up with accurate parameter!zations of the differential CBL growth.
Brehm [51], based on the Innsbruck, Austria, Slopewind Experiment of 1978,
has presented a conceptual model snowing differential CBL growth on opposite
sidewalls (Figure 27). His conceptual model also shows subsidence in the
stable core. Bader and McKee [47], however, in their two-dimensional simu-
lation of valley inversion breakup, noted no significant differential
boundary layer growth over opposite valley sidewalls even when substantially
different sensible heat flux functions were used on the different valley
sidewalls. They explained this as due to the effects of horizontally
propagating gravity waves that act to keep isotherms horizontal within the
stable core. They also noted the existence of multiple closed cross valley
circulations within the stable core early in the post-sunrise period that
acted to warm the valley temperature inversion differentially with height,
destabilizing the valley atmosphere. This process acts in consort with
stable core subsidence to warm the valley atmosphere. As shown by the
different conclusions reached in these two studies, it seems appropriate to
verify differential CBL growth using specially designed field studies.
Stacks Located on Sidewalls
The model could be modified to handle emissions from stacks located
above the valley floor on the sidewalls. In view of the results from the
1982 tracer experiments, this modification should be combined with a new
treatment of differential heating and cross valley flows, since the
development of cross valley flows due to differential heating may strongly
affect concentrations on the sidewalls in the vicinity of the elevated
stack.
106
-------
BRUSH CR.-AUGUST 1, 1982
1400 r
E-FAC1NG SIDEWALL
W-FACING SIDEWALL
10 12 14 16 18 20
Time (MDT)
Figure 26. Illustration of differential solar flux on opposing sidewalls
for Brush Creek Colorado on August 4, 1982. The Brush Creek
Valley drains from NW to SE. The curves represent extra-
terrestrial solar flux (i.e., they assume no atmosphere) and
were determined using a solar radiation model.
107
-------
®
MOUNTAIN
WIND
VALLEY
WIND
SLOPE
WIND
STABLE
LAYER
Figure 27. Brehm's [51] conceptual model of temperature inversion
destruction in Austria's Inn Valley, showing
differential CBL growth over the opposite sidewalls and
continued down-valley flow in the elevated stable core.
108
-------
REFERENCES
1. Hewson, E. W., and G. C. Gill. 1944.
ver Valley Near Trail , B.C.
Columbia Ri
of Interior,
453, 23-228.
"Meteorological Investigations in
" Bur. Mines Bull., U.S. Department
2. Bierly, E. W., and E. W. Hewson. 1962. "Some Restrictive Meteorological
Conditions to be Considered in the Design of Stacks." J. Applied Meteor.,
1, 383-390.
3. Whiteman, C. D., and T. B. McKee. 1982. "Breakup of Temperature
Inversions in Deep Mountain Valleys: Part II. Thermodynamic Model." j_._
Applied Meteor., 21, 290-302.
4. U.S. Environmental
Office of Research
Protection Agency. 1978. Handbook for Preparing
and Development Reports, Revised.
EPA-600/9-78-032,
Cincinnati, Ohio.
December 1978, U.S. Environmental
Report
Protection
Agency,
5. Peterson, W.
Handbook for
B., J. S. Irwin,
Preparing Users'
D. B. Turner and J. A.
Guides for Air Quality
U.S. Environmental
North Carolina.
Protection Agency, Research Triangle Park,
Catalano. 1982.
_Mod_e_ls. August 1982,
6.
7.
8.
Whiteman, C. D.
Rocky Mountains.
Nov. 9-12, 1981,
1981. "Temperature Inversion Buildup in Valleys of the
1 Preprints, Second Conference on Mountain Meteorology.
Steamboat Springs, Colorado, pp. 276-272.
Whiteman, C. D. and T. B. McKee. 1977. "Observations of Vertical
Atmospheric Structure in a Deep Mountain Valley." Arch. Meteor. Geophys.
Bioklimat., Ser. A, 26, 29-50.
Whiteman, C. D., and T. B. McKee. 1979. Temperature Inversion
Destruction in Mountain Valleys - Implications for Air Pollution
Dispersion. Proceedings of the Tenth International Technical Meeting
Air Pollution Modeling and Its Application, NATO/CCMS Pilot Study,
Oct. 23-26, 1979, Rome, Italy, 217-224.
on
9.
Whiteman
Mountain
C. D.
Valleys.
ey;
Til
1980. Breakup
Ph.D.
of Temperature Inversions in Colorado
Dissertation, Colorado State University,
pp. 250.Available from University Microfilms, Ann Arbor, Michigan, or as
Atmospheric Science Paper No. 328 from Department of Atmospheric Science,
Colorado State University, Fort Collins, Colorado.
109
-------
10. Whiteman, C. D. 1982. "Breakup of Temperature Inversions in Deep
Mountain Valleys: Part I. Observations." J. Applied Meteorology, 21,
270-289.
11. Whiteman, C. 0., and T. 8. McKee. 1978. "Air Pollution Implications of
Inversion Descent in Mountain Valleys." Atmospheric Environment, 12,
2151-2158.
12. Whiteman, C. D., and T. 8. McKee. 1981. "Valley Wind Systems Under
Temperature Inversion Conditions." Proceedings, Symposium on Intermediate
Range Atmospheric Transport Processes and Technology Assessment.
October 1-3, 1980, Gatlingburg, Tennessee, pp. 351-358.
13. Whiteman, C. D., and R. L. Drake. 1981. The Green River Ambient Model
Assessment Program September 1981 Progress Report on Local Scale Modeling
Activities. PNL-4058. September 1981. Pacific Northwest Laboratory,
Richland, Washington, 249 pp.
14. Whiteman, C. D., and K. J. Allwine. 1982. Green River Ambient Model
Assessment Program FY-1982 Progress Report. PNL-4520. October 1982.
Pacific Northwest Laboratory, Richland, Washington, 187 pp.
15. Davidson, 8., and P. K.. Rao. 1958. Preliminary Report on Valley Wind
Studies in Vermont, 1957. Final Report AFCRC-RT-58-29, Contract AF
19(604)-1971, College of Engineering, New York University, 54 pp.
16. Beuttner, J. K., and N. Thyer. 1966. Valley Winds in the Mount Rainier
Area. Arch. Meteor. Geophys. Bioklimat., Ser. B., 14(2):125-147.
17. Tyson, P. 0., and R. A. Preston-Whyte. 1972. "Observations of Regional
Topographically-Induced Wind Systems in Natal." J. Applied Meteor., 11,
643-650.
18. Sterten, A. K. 1963. A Further Investigation of the Mountain and Valley
Wind System in South-Eastern Norway. Intern Rapport K-254, Forsvarets
Forskningsinstitutt, Norwegian Defence Research Establishment, Kjeller,
Norway, 51 pp.
19. Urfer-Henneberger, C. 1970. Neuere Beobachtungen uber die Entwicklung
des Schonwetterwindsystems in einem V-formigen Alpental (Dischmatal bei
Davos). Arch. Meteor. Geophys. Bioklimat., Ser. B, 18, 21-42.
20. Start, G. E., C. R. Dickson and L. L. Wendell., 1975. "Diffusion in a
Canyon Within Rough Mountainous Terrain." J. Applied Meteor., 14(3):333-
346.
21. Briggs, G. A. 1969. Plume Rise. AEC Critical Review Series, TID-25075.
22. Briggs, G. A. 1971. "Some Recent Analyses of Plume Rise Observations."
In: Proceedings of the Second International Clean Air Congress,
H. M. Englund and W. T. Beery, eds. Academic Press, New York,
pp. 1029-1032.
110
-------
23. Briggs, G. A. 1973. "Diffusion Estimation for Small Emissions." Atmos.
Turbulence and Diffusion Laboratory, Contribution File No. (Draft) 79.
Oak Ridge, Tennessee, 59 pp.
24. Briggs, G. A. 1975. "Plume Rise Predictions." Lectures on Air Pollution
and Environmental Impact Analysis, D. A. Haugen, ed. American
Meteorological Society, Boston, Massachusetts, pp. 59-111.
25. Pierce, T. E., and D. B. Turner. 1980. User's Guide^for MPTER: A
Multiple Point Gaussian Dispersion AlgoriThm with Optional Terrain
Adjustment. EPA-600/8-80-016. April 1980, U.S. Environmental Protection
Agency, Research Triangle Park, North Carolina.
26. Turner, D. B. 1969. Workbook on Atmospheric Dispersion Estimates.
Public Health Service Publications No. 999-AP-26.
27. Slade, D. H. (ed.). 1968. Meteorology and Atomic Energy 1968. U.S.
Atomic Energy Commission, Division of TechnicalInformation, Oak Ridge,
Tennessee, 445 pp.
28. McMullen, R. W. 1975. "The Change of Concentration Standard Deviations
with Distance." J. Air Poll. Control Assoc., 25(10):1057-1058.
29. Minnott, D. H., D. L. Shearer and Y. Gotaas. 1979. "Measurements of the
Vertical Dispersion Rate in Deep-Valley Terrain." Preprint Volume:
Fourth Symposium on Turbulence, Diffusion and Air Pollution. January 15-
18, 1979, Reno, Nevada. Published by American Meteorological Society,
Boston, Massachusetts.
30. Minnott, D. H., D. L. Shearer and and R. S. Marker. 1977. Development of
Vertical Dispersion Coefficients for Deep-Valley Terrain. Report
COO-4026-3 prepared for U.S. Energy Research and Development Administra-
tion by TRC Environmental Consultants, Inc.
31. Shearer, D. L., D. H. Minnott and G. R. Hilst. 1977. Development of
Vertical Dispersion Coefficients for Rolling Terrain Environments. The
Research Corporation of New England (TRC) Report COO-4026-1, January 1977.
32. Gotaas, Y. 1972. A Model of Diffusion in a Valley from a Continuous
Point Source." Arch. Meteor. Geophys. Bioklimat., Ser. A, 21. 13-26.
33. Abramowitz, M., and I. A. Stegun. 1965. Handbook of Mathematical
Functions with Formulas, Graphs, and Mathematical Tables. Dover
Publications, Inc., New York, 1046 pp.
34. Hastings, C., Jr. 1955. Approximations for Digital Computers. Princeton
University Press, Princeton, New Jersey.
35. Sellers, W. D. 1965. Physical Climatology. The University of Chicago
Press. Chicago, Illinois, 272 pp.
Ill
-------
36. Kreith, F., and J. F. Kreider. 1978. Principles of Solar Engineering.
Hemisphere Publishing Corp., Washington, D.C., 778 pp.
37. McCullough, E. C. 1968. Total Daily Radiant Energy Available Extrater-
restrial ly as a Harmonic Series in the Day of the Year. Arch. Meteor.
Geophys. Bioklimat.. Ser. B. 16. 129-243.
38. Waugh, A. E. 1973. Sundials: Their Theory and Construction. New York,
Dover Publications, Inc., 228 pp.
39. Thompson, P. D. 1961. Numerical Weather Analysis and Predictions. The
MacMillan Co., 170 pp.
40. Haltiner, G. J. 1971. Numerical Weather Prediction. John Wiley and
Sons, Inc., New York, 317 pp.
41. Friedman, F. L., and E. B. Koffman. 1981. Problem Solving and Structured
Programming in FORTRAN, Second Edition. Addison-Wesley Publishing Co.,
Reading, Massachusetts.
42. MUller, H., R. Reiter, R. Sladkovic, and K. Munzert. 1982. "Aerologische
Untersuchungen des Tagesperiodischen Windsystems im Loisachtal." XVII
Internationale Tagung fur Alpine Meteorologie. Berchesgaden, 21-25
September 1982. Annalen der Meteorologie, (Neue Folge), Nr. 19, 186-188.
43.. Hanna, S. R., G. A. Briggs and R. P. Hosker, Jr. 1982. Handbook on
Atmospheric Diffusion. DOE/TIC-11223. Technical Information Center, U.S.
Department of Energy. 102 pp.
44. Sundararajan, A., and S. A. Macklin. 1976. "Comments on the Heat Flux
and Friction Velocity in Free Convection Near the Ground." J. Applied
Meteor., 33, 715-718.
45. Gleeson, T. A. 1951. "On the Theory of Cross-Valley Winds Arising From
Differential Heating of the Slopes." J. Meteor., 8(6), 398-405.
46. Bader, D. C. 1981. Simulation the Daytime Boundary Layer Evolution in
Deep Mountain Valleys. M.S. Thesis, Atmospheric Science Department,
Colorado State University, Fort Collins, Colorado.
47. Bader, D. and T. B. McKee. 1983. "Dynamical Model Simulation of the
Morning Boundary Layer Development in Deep Mountain Valleys." J. Climate
and Appl. Meteor, 22(3), 341-351.
48. McNider, R. T. 1981. "Investigation of the Impact of Topographic
Circulations on the Transport and Dispersion of Air Pollutants." Ph.D.
Thesis. Department of Environmental Science, University of Virginia,
Charlottesville, Virginia (available from University Microfilm, Ann Arbor,
Michigan).
112
-------
49. McNider, R. T., and R. A. Pielke. 1981. "Diurnal Boundary-Layer
Development Over Sloping Terrain." J. Atmos. Sci.. 38(10), 2198-2212.
50. Lenschow, D. H., and B. B. Stankov. 1979. "The Rapid Morning Boundary-
Layer Transition." J. Atmos. Sci., 36, 2108-2124.
51. Brehm, M. 1981. Hangwindexperiment Innsbruck 1978 - Gebirgswindsystem
und Inversionsauflosung.Ph.D. Dissertation, Ludwig-Maximinans-
Universitat, Munich, West Germany.
52. Whiteman, C. D., R. N. Lee, M. M. Orgill and B. D. Zak. 1984. Green
River Air Quality Model Development. Meteorological and Tracer Data -
July/August 1982 Field Study in Brush Valley, Colorado. PNL-5163.
Pacific Northwest Laboratory, Richland, Washington, 139 pp.
53. Whiteman, C. D., A. H. Huber, R. W. Fisher, and B. D. Zak. 1984.
"Atmospheric Tracer Experiments in a Deep Narrow Valley." Conference
Volumes: (Joint Session) Fourth Joint Conference on Applications of Air
Pollution Meteorology and Third Conference on Mountain Meteorology,
October 16-19, 1984, Portland, Oregon, pp. J1-J4.
54. Orgill, M. M., R. N. Lee, R. I. Schreck, K. J. Allwine and C. D. Whiteman.
1984. "Early Morning Ventilation of an SF§ Tracer from a Mountain
Valley." Conference Volumes: (Joint Session) Fourth Joint Conference on
Applications of Air Pollution Meteorology and Third Conference on Mountain
Meteorology, October 16-19, 1984, Portland, Oregon, pp. J36 J-39.
55. Allwine, K. J., and C. D. Whiteman. 1984. Technical Description of
MELSAR: A Mesoscale Air Quality Model for Complex Terrain. PNL-5048,
Pacific Northwest Laboratory, Richland, Washington.
113
-------
APPENDIX A
FORTRAN LISTING OF VALMET
115
-------
1 PROGRAM VALHET
2 c****************************************************************************
3 C**** THIS COMPUTER MODEL PREDICTS AIR POLLUTION CONCENTRATIONS ON A *
4 c**** VALLEY FLOOR AND SIDEWALLS ARISING FROM AN ELEVATED CONTINUOUS *
5 C**** POINT SOURCE OF POLLUTION DURING NOCTURNAL STEADY-STATE *
g c**** DOWN-VALLEY DRAINAGE FLOWS AND DURING THE POST-SUNRISE *
7 C**** TEMPERATURE INVERSION BREAKUP PERIOD. A MODIFIED GAUSSIAN PLUME *
8 C**** ALGORITHM IS USED DURING THE NOCTURNAL PERIOD. A NUMERICAL SCHEME *
9 C**** IS USED FOR THE CALCULATIONS AFTER SUNRISE WHEN FUMIGATIONS OCCUR. *
10 C**** SIMPLE PAP.AMETERIZATIONS ARK USED IN THE POST-SUNRISE SIMULATIONS +
11 C**«* TO ACCOUNT FOR SOLAR FORCING, SENSIBLE HEAT FLUX, CBL GROWTH,
12 C**** INVERSION DESCENT, AND UPSLOPE TRANSPORT AND DIFFUSION. A USER'S
13 C**** MANUAL IS NOW AVAILABLE FOR THIS MODEL. CONTACT C.D. WHITEMAN
14 c**** OR K.J. ALLWINE, GEOSCIENCES RESEARCH AND ENGINEERING DEPARTMENT,
15 c**** BATTELLE PACIFIC NORTHWEST LABORATORIES, RICHLAND.WA 99352.
16 C**** VERSION 1.1 DECEMBER 1, 1984
^7 £****************************************************************************
18
19 PARAMETER(NS=6000,NA=100,NB=30,NB1-31)
20
21 COMMON/BLK1/HC(NS),HT(NS)
22 COMMON/BLK2/CHIBAR(NB)
23 COMMON/BLK3/V(NB1)
24 COMMON/BLK4/CONC(NA,NB)
25 COMMON/BLK5/CHKNB)
26 COMMON/BLK6/AVG(NB,2),NDXTIM(NB,2),NTS(2)
27
28 REAL HITEL(NB),HITEU(NB),TIME(NS)
29 REAL K,LAT
30 REAL SUM(NB),CCH(2,NB)
31 CIIARACTER*3 ITIME
32
33 c**** OPEN FILES
34 OPEK(UNIT=2,NAME='VALMET.OUT',TYPE='NEW',FORM='FORMATTED1)
35 OPEN(UNIT=3,NAME='VALMET.PLT',TYPE='NEW',FORM='FORMATTED')
36
37 c**** INITIALIZE PARAMETERS
38 HCBLI=25.
39 ITIME='LST'
40 C***« STANDARD TEMPERATURE AND PRESSURE
41 PSTD=1013.
42 TSTD=20.+273.16
43
44 c**** READ INPUT DATA FROM SUBROUTINE
45 CALL INPUT(HCBLI,H,W,ALPHA1,ALPHA2,DELTAT,NBOX,HZERO,
46 + GAMMA,BETA,LAT,URLONG,MO,IDA,IYR,K,RHO,PRESS,
47 + Q,ISTAB,XC,UC,STMP,SRAD,SVEL,TEMP,AO,YZERO,US,
48 + AS,BOXLEN,NTOTI,AC)
49 C NUMBER OF TIMESTEPS IN THE CONCENTARTION AVERAGING PERIOD
50 NTSA=30
51 C CONCENTRATION AVERAGING INTERVAL IN SECONDS
52 IAVG=NTSA*INT(DELTAT)
53 C FACTOR TO CONVERT TO STANDARD CONDITIONS
54 STNDCON=PSTD/PRESS*TEMP/TSTD
55
56
57 c**** CALCULATE JULIAN DATE
58 CALL JULIAN(MO,IDA,IYR,JULDAY)
59
60 c**** PLUME RISE
61 CALL PRISE(TEMP,SVEL,SRAD,STMP,XC,US,ISTAB,GAMMA,DLH)
62 H=H+DLH
63 IF (H. GT. HZERO. OR. H. I.E. HCBLI) THEN
64 WRITE(6,115)
117
-------
65 115 FORMAT(1HO, "THE PLUME, AFTER PLOMERISE, IS NOT WITHIN',
66 +' THE STABLE CORE1)
67 STOP 4
68 ENDIF
69
70 c**** ADJUST STABILITY CATEGORIES TO ACCOUNT FOR ENHANCED
71 C DIFFUSION IN COMPLEX TERRAIN, AFTER START ET AL. (1975) .
72 ISTABY=ISTAB-2
73 ISTABZ-ISTAB-1
74
75 C**** VOLUMETRIC DILUTION DUE TO SLOPE FLOW CONVERGENCE, TRIBUTARY
76 C**** FLOWS, OR ENTRAINMENT.
77 CALL DILUTE(US,UC,AS,AC,DF)
78
79 C**** DETERMINE DIFFUSION COEFFICIENTS AND CALCULATE THE FRACTION
80 c**** OF POLLUTANT MASS IN THE VALLEY CROSS-SECTION THAT HAS
81 c**** DIFFUSED "BEYOND" THE VALLEY WALLS AND "BELOW" THE GROUND.
82 CALL INGRATJXC,YZERO,ISTABY,ISTABZ,W,ALPHA1,ALPHA2,
83 + HZERO,H,SUM1,SUM2,SIGMAY,SIGMAZ)
84
85 c**** DISTRIBUTE THE LOST MASS FRACTION BACK INTO THE CROSS-
86 c**** SECTION. USE WELL-MIXED ASSUMPTION.
87 CHIOFF=(1.-SUM1-SUM2)*Q*1.E09/(US*DF*AC)
88
89 c**** USE GAUSSIAN PLUME EQUATION TO CALCULATE THE MAXIMUM
90 c**** CENTERLINE CONCENTRATION WITHIN THE VALLEY CROSS-SECTION
91 CALL GAUSS(YZERO,H,YZERO,H,Q,US,DF,SIGMAY,SIGMAZ,
92 + CLCONC)
93 CLCONC=(CLCONC+CHIOFF)*STNDCON
94
95 c**** USING GAUSSIAN PLUME EQUATION, CALCULATE STEADY-STATE
96 c**** NOCTURNAL CONCENTRATIONS AT POINTS ON THE VALLEY FLOOR AND
97 c**** SIDEWALLS CORRESPONDING TO THE GRID ELEMENTS USED FOR THE
98 c**** AFTER-SUNRISE CALCULATIONS.
99 DO 1 I=1,NTOTI
100 IF(I.LE.NBOX) THEN
101 C VALLEY FLOOR
102 Z»0.
103 Y=(BOXLEN/2.)+BOXLEN*(I-D
104 CALL GAUSS(Y,Z,YZERO,H,Q,US,DF,SIGMAY,SIGMAZ,CCHI)
105 CCHI=CCHI+CHIOFF
106 C SIDEWALLS
107 ELSE
108 Z=(BOXLEN*TAN(ALPHA1))/2.+(I-NBOX-1)*BOXLEN*TAN(ALPHA1)
109 Y=W/2.+Z/TAN(ALPHAl)
110 CALL GAUSS(Y,Z,YZERO,H,Q,US,DF,SIGMAY,SIGMAZ,CCHI)
111 CCHI=CCHI+CHIOFF
112 ENDIF
113 CHI(I)=CCHI
114 1 CONTINUE
115
116 C**** THE FIRST ELEMENT IN THE MODEL OUTPUT CONCENTRATION ARRAY IS
117 C**** DEFINED AS THE STEADY-STATE NOCTURNAL CONCENTRATION.
118 DO 2 I=1,NTOTI
119 CONC(1,I)=CHI(I)
120 CCH(1,I)=CHI(I)
121 2 CONTINUE
122
123 C**** CALCULATE SOLAR PARAMETERS FOR THE SITE AND DATE.
124 CALL SOLAR(LAT,URLONG,JULDAY,NTSR,TAU,A1)
125 TAU=TAU*3600.
126 T4=NTSR
127 IT=NTSR/100
128 T2=IT
118
-------
129 T1=T4/100
130 T3=(T1-T2)*100./60.+T2
131
132 c**** ENERGY BUDGET
133 CALL EBDGT(AO,A2)
134 AO=A2
135
136 C**** USE BULK THERMODYNAMIC MODEL TO DETERMINE INVERSION TOP
137 c**** DESCENT AND CBL RISE AS FUNCTION OF TIME. THE NUMBER OF
138 C**** TIME STEPS REQUIRED TO DESTROY SHE INVERSION IS ONE OF
139 c**** THE OUTPUTS OF THE THERMODYNAMIC MODEL.
140 CALL DESCNT(AO,A1,ALPHA1,ALPHA2,BETA,DELTAT,GAMMA,
141 + HCBLI,HZERO,K,W,PRESS,RHO,TAU,NSTEPS)
142 NINDX=NSTEPS+NINT(3600./DELTAT)
143 NAS=NINDX/NTSA-H
144
145 C**** CALCULATE AIR POLLUTION CONCENTRATION AS A FUNCTION OF TIME
146 C**** AFTER SUNRISE.
147 N=l
148 DO 9 IND=1,NAS
149
150 C INITIALIZE CONCENTRATION SOMMING ARRAY TO ZERO
151 DO 3 I=1,NTOTI
152 SUM(I)=0.
153 3 CONTINUE
154
155 C LOOP ON NUMBER OF TIMESTEPS PER AVERAGING PERIOD
156 DO 7 ND«1,NTSA
157 N-N+1
158 MN«N-1
159
160 IF(N.GT.NSTEPS)THEN
161 DO 4 I=1,NTOTI
162 V(I)=999.9
163 4 CONTINUE
164 GO TO 5
165 ENDIF
166
167 IFIX2=(HT(MN)-HC(MN))/(BOXLEN*TAN(ALPHA1))
168 NTOT=NBOX+IFIX2
169
170 C**** CALCULATE POLLUTANT MASS COMING INTO GROWING BOX AT ITS TOP.
171 CALL PROFIL(MBOX,NTOTI,NTOT,BOXLEN,HC(MN),ALPHA1,
172 + ALPHA2,HZERO,H,HC(N),HT(MN),HT(N),YZERO,Qr
173 + US,DF,SIGMAY,SIGMAZ,CHIOFF)
174
175 c**** CALCULATE VELOCITIES AT BOX SEPTA.
176 CALL VELOCY(N,MN,NTOTI,NBOX,NTOT,BOXLEN,DELTAT)
177
178 c**** CALCULATE POLLUTANT MASS BUDGET FOR EACH BOX.
179 5 CALL BRKUP(N,MN,NTOTIJ.BOXLEN,DELTAT,CCH)
180
181 DO 6 I»1,NTOTI
182 SUM(I)=CCH(2,I)+SUM(I)
183 CCH(1,I)=CCH(2,I)
184 6 CONTINUE
185
186 7 CONTINUE
187
188 DO 8 I=1,NTOTI
189 CONC(IND,I)=SUM(I)/REAL(NTSA)*STNDCON
190 8 CONTINUE
191
119
-------
192 IF(IND.EQ.1)THEN
193 TIME(IND)=T3 +REAL(NTSA)*DELTAT/2./3600.
194 ELSE
195 TIME(IND)=TIME(IND-1)+REAL(NTSA)*DELTAT/36 00.
196 ENDIF
197
198 9 CONTINUE
199
200 C**** EMPLOY THE POST-PROCESSOR TO CALCULATE MAXIMUM 1- AND 3-HOUR
201 C**** AVERAGE CONCENTRATIONS FOR REGULATORY APPLICATIONS.
202 TSTP=REAL(NTSA)*DELTAT
203 CALL PSTPRC(TSTP,NTOTI,NAS)
204
205 c**** WRITE OUT RESULTS TO UNIT 2.
206 WRITE(2,100) DELTAT,NTOTI,BOXLEN,NSTEPS,NINDX
207 100 FORMAT(1HO,'NUMERICAL SIMULATION PARAMETERS:',/,IX,
208 +'TIMESTEPS ARE',F5.1,' SECONDS LONG',/,IX,'NUMBER OF GRID',
209 +' ELEMENTS"1,13,/,IX,'GRID ELEMENTS ARE',F6.1,' METERS',
210 -I-1 LONG',/, IX,'NUMBER OF TIMESTEPS FROM SUNRISE TO INVERSION',
211 +' DESTRUCTION"1,I6,/,IX,'NUMBER OF TIMESTEPS IN SIMULATION1,
212 +' (INCLUDES EXPONENTIAL DECAY AFTER INVERSION DESTRUCTION)=',
213 +16)
214 WRITE(2,101) IAVG,NAS
215 101 FORMATdH ,'CONCENTRATION AVERAGING INTERVAL" ' , 15 , ' SECONDS',
216 +/flH ,'NUMBER OF AVERAGING INTERVALS IN THE SIMULATION=',14)
217 WRITE(2,102) CLCONC,H,XC,(SUM1*100.)
218 102 FORMAT(1HO,'PLUME DIFFUSION DURING NOCTURNAL TRAVEL:',/,lH ,
219 + 'A PLUME CENTERLINE CONCENTRATION OF',F7.2,
220 + ' MICROGRAMS PER CUBIC METER OCCURS',F6.0,' M ABOVE ',
221 •*• 'THE VALLEY FLOOR CENTER. ' ,/,lH , ' AT THE TRAVEL DISTANCE OF1
222 + ,F6.1,' KM',F5.1,'% OF THE PLUME IS CONTAINED WITHIN THE'
223 + ' VALLEY INVERSION CROSS SECTION.')
224 WRITE (2,103) XC,(SUM2*100.),((1.-SUM1-SUM2)*100.),CHIOFF
225 103 FORMATdH , ' AT THE TRAVEL ',
226 + 'DISTANCE OF',F6.1,' KM',F5.1,'% OF THE PLUME HAS DIFFUSED',
227 + ' OUT THE TOP OF THE VALLEY.',/,lH ,'THE REMAINING',F5.1,
228 + '% OF THE PLUME MASS IS MIXED UNIFORMLY THROUGHOUT THE ',
229 + 'INVERSION CROSS',/,1H ,'SECTION TO PRODUCE AN OFFSET COHC',
230 + 'ENTRATION OF',F7.3,' MICROGRAMS PER CUBIC METER.1)
231 WRITE (2,104) (AC/1000.),(1./DF)
232 104 FORMATdH ,'AREA AT SIMULATION CROSS SECTION IN THOUSANDS OF ',
233 + 'SQ METERS IS',F8.0,/,lH ,'CLEAN AIR DILUTION FACTOR IS',F8.3)
234 WRITE(2,105) NTSR,ITIME,(TAU/3600.),Al
235 105 FORMATdHO,'SOLAR MODEL RESULTS:',/,lH ,'TIME OF SUNRISE=',
236 + I6,1X,A3,/,IX,'LENGTH OF DAYLIGHT PERIOD=',F6.2,' HRS',/,
237 + 1H ,'EXTRATERRESTRIAL SOLAR FLUX ON HORIZONTAL SURFACE AT',
238 + ' SOLAR NOON=',F7.1,' WATTS PER SQUARE METER')
239 WRITE(2,106)
240 106 FORMAT(1HO,'MODEL OUTPUT: CONCENTRATIONS IN MICROGRAMS/CUBIC',
241 + ' METER:1,/,1H ,2X,'GRID',5X,'HEIGHT ABOVE',3X,'MAXIMUM1,4X,
242 + 'TIME OF',5X,'MAXIMUM1,4X,'TIME OF',/,1H ,IX,'ELEMENT1,
243 + 3X,'VALLEY FLOOR',4X,'1-HOUR',2X,'OCCURRENCE1,5X,'3-HOUR',
244 -I- 2X,'OCCURRENCE',/,lH ,IX, ' NUMBER' ,9X, ' (M) ' ,7X,'AVERAGE',
245 + 5X, ' (LST) ' ,6X,'AVERAGE',5X, ' (LET)')
246
247 DO 10 I=1,NTOTI
248 IF(I.LE.NBOX)THEN
249 HITEL(I)=0.
250 HITEU(I)=0.
251 ELSE
252 HITEL(I) = (I-NBOX-1) *BOXLEN*TA£J(ALPHA1)
253 HITEU(I)=(I-NBOX)*BOXLEN*TAN(ALPHA1)
254 ENDIF
255 10 CONTINUE
120
-------
256
257 DO 11 1=1,NTOTI
258 TMAXC1=TIME(NDXTIM(I,1))
259 TMAXC2=TIME(NDXTIH(I,2))
260 WRITE(2,107) I,HITEL(I),HITEU(I),AVC{I,1),TMAXC1-1.,
261 + TMAXC1,AVG(I,2),T«AXC2-3.,TMAXC2
262 107 FORMATdH , 14 ,9X,F4 .0 , ' - ' ,F4 .0 ,2 (2X,F8.2 ,2X, F5 .2 ,'-' , F5 .2 ) )
263 11 CONTINUE
264
265 WRITE(2,108) IAVG,IAVG
266 108 FORMAT(1H1,//,15X,I5,' SECOND AVERAGE POLLUTANT CONCENTRATIONS',
267 + ' (ug/m3) IN MODEL GRID ELEMENTS AS A FUNCTION OF TIME',//
268 + 1H ,32X,'INDICATED TIME IS AT THE MIDPOINT OF',15,' SECOND',
269 + ' AVERAGING PERIOD',/,1H ,32X,'PRE-SUNRISE STEADY STATE',
270 + ' CONCENTRATIONS ARE GIVEN IN FIRST LINE')
271 WRITE(2,109)
272 109 FORMAT (1 HO, 2X, 'TIME' ,2X, (20CBOX1 ,3X) ) ,/,lH ,2X, ' LST',4X,
273 +'01',4X,'02',4X,'03',4X,'04',4X,'05',4X,'06',4X,'07',4X,
274 +108',4X,'09',4X,110I,4X,'11I,4X,'121,4X,'13',4X,'14',4X,
275 +'15',4X,'16',4X,'17',4X,'18',4X,'19',4X,'20')
276
277 IF (NTOTI. LE. 20) THEN
278 NY=NTOTI
279 ELSE
280 NY=20
281 ENDIF
282
283 WRITE(2,110) T3,(CHI(I),1=1,NY)
284 110 FORMATdH ,1X, ' < ' ,F4 .2 ,20F6 .2)
285 DO 12 MN-1,NAS
286 WRITE(2,111) TIME(MN),(CONC(MN,I),1=1,NY)
287 111 FORHAT(1H ,21F6.2)
288 12 CONTINUE
289
290 IF(NTOTI.GE.21)THEN
291 WRITE(2,112)
292 112 FORMAT(1H1,//////,2X,'TIME',2X, (20CBOX1 ,3X)) ,/,lH ,2X,
293 + 'LST',4X,'21r,4X,'22',4X,'23',4X,'24',4X,'25',4X,'26',
294 +4X,'27',4X,'28',4X,'29',4X,'30',4X,'31',4X,'32',4X,'33',
295 +4X,'34',4X,'35',4X,'36',4X,'37',4X,'38',4X,'39',4X,'40')
296 WRITE(2,113) T3,(CHI(I),1=21,NTOTI)
297 113 FORMATdH ,1X, ' < ' , F4 .2 ,20F6 .2)
298 DO 13 MN=1,NAS
299 WRITE(2,114) TIME(MN),(CONCfMN,I),1=21,NTOTI)
300 114 FORMATdH ,21F6.2)
301 13 CONTINUE
302 ENDIF
303
304 CLOSE(UNIT=2)
305
306 C**** WRITE OUT RESULTS TO UNIT 3 FOR PLOTTING.
307 11=1
308 TIME(1)=T3
309 WRITE(3,201) NSTEPS,NTSR,NINDX,NTOTI,NAS,ITIME
310 201 FOSMATdH ,518,A3)
311 WRITE(3,202) CLCONC,DELTAT,SUM1,SUM2,XC,CHIOFF,
312 + SIGMAZ,SIGMAY,ALPHA1,ALPHA2,W,HZERO,H
313 202 FORMATdH ,13F10.5)
314 DO 14 1=1,NTOTI
315 WRITE(3,203) I,II,T2,CHI(I)
316 DO 14 IKD=1,NAS
317 WRITE(3,203) I,IND+1,TIME(IND),CONC(IND,I)
318 203 FORMATdH , 12 ,15 , F12 .7 , F12 .7)
319 14 CONTINUE
121
-------
320 DO 15 I=1,NSTEPS
321 WRITE(3,204) HC(I),HT(I)
322 204 FORMAT(1H ,6(2F9.3))
323 15 CONTINUE
324 CLOSE(UNIT=3)
325
326 STOP 'NORMAL EXIT1
327 END
122
-------
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
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
SUBROUTINE INFUT(HCBLI, H, W, ALPHA1 , ALPHA2 ,DELTAT,
1 NBOX, HZ ERO, GAMMA, BETA, LAT, URLONG, MO, IDA, IYR, K, RHO,
2 PRESS, Q, ISTAB, XC, UC, STMP, SRAD, SVEL, TEMP, AO , YZERO, US,
3 AS,BOXLEN,NTOTI, AC)
PARAMETER (NS=6000 , NA=100 , NB=30 ,NB1=31 )
C**** THIS SUBROUTINE CONTROLS AN INTERACTIVE TERMINAL THROUGH WHICH
C**** INPUT DATA ARE RECEIVED BY THE MODEL.
REAL K, LAT, VAL(27), VALNEW(27), VALLO(27), VALHI(27)
REAL ASTAB(6)
INTEGER ID (27)
CHARACTER*? NAM (27)
CHARACTER*20 CAT (10)
*********** Jjjpprp pgfJfljrpjQjqg***** ****************************
VALLEY SIDEWALL SLOPE (DEC ELEVATION) *
OTHER VALLEY SIDEWALL SLOPE (DEG) *
CROSS SECTIONAL AREA AT SOURCE (THOUSANDS OF M2) *
RATE OF TEMP RISE ABOVE VALLEY (DEG K/SEC)
TIME STEP (SEC)
VERTICAL POTENTIAL TEMPERATURE GRADIENT (DEG K/M)
INITIAL HT OF POLLUTION MAX CONCENTRATION(M) .
INITIAL INVERSION HEIGHT(M)
DAY OF MONTH (1-31)
ATMOSPHERIC STABILITY CLASS (1-6), E.G., 6=F
YEAR (00-99)
FRACTION OF SENSIBLE HEAT GOING INTO CBL GROWTH (0-1)
LATITUDE (DEG N)
MONTH OF YEAR (1-12)
NUMBER OF BOXES ON HALF OF THE VALLEY FLOOR.
PRESSURE (MB)
POLLUTANT SOURCE STRENGTH (KG/SEC) *
DENSITY (KG/M**3) *
STACK RADIUS («) *
STACK EFFLUENT TEMPERATURE (DEG C) *
STACK EXIT VELOCITY (M/S) *
AMBIENT AIR TEMPERATURE (DEG C) *
DOWNVALLEY WIND VELOCITY AT CROSS SECTION (M/SEC) *
LONGITUDE (DEG W) *
PLUME-CARRYING DOWNVALLEY WINDSPEED AT STACK (M/S) *
FULL WIDTH OF VALLEY FLOOR (M) *
OFF-CRNTERLINE DISPLACEMENT OF STACK (M) *
DOWNVALLEY DISTANCE FROM SOURCE TO CROSS SECTION (M) *
C************************************************************************
C**** SPECIFY NAMES OF CATEGORIES AND PARAMETERS.
c***********
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
ALPHA1
ALPHA2
AS
BETA
DELTAT
GAMMA
H
HZ ERO
IDA
ISTAB
IYR
K
LAT
MO
NBOX
PRESS
Q
RHO
SRAD
STMP
SVEL
TEMP
UC
URLONG
US
W
YZERO
XC
DATA
1
2
3
4
DATA
1
2
3
4
5
6
CAT /'SENSIBLE HEAT FLUX
'MODEL CHARACTER.
'SITE LOCATION
'VALLEY ATMOSPHERE
•STACK CHARACTER.
NAM /'AO
'ALPHA2=
1 BETA =
1 IYR
'UC
'Q
1 SRAD =
'K
'NBOX =
'LAT
1 TEMP =
' STAB =
1 SVEL =
'VALLEY CHARACTER.
'INVERSION CHARACTER.
'DATE
'GAUSSIAN PLUME
'
•WIDTH =' , 'ALPHA1=' ,
•HZERO =' , 'GAMMA =' ,
1 URLONG= ' , ' MO = ' , ' IDA!
•PRESS =' , 'XC =', 'YZEF
'US = ' , ' H = ' , ' STMI
'AS ='/
123
-------
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
^78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
C**** SET DEFAULT VALUES.
DATA VAL /. 24, .15, 600
.,15. ,15. ,3. ,500. ,.025,0.,
1 40. 00, 105. 00, 9., 21. ,82., 10., 7 50. ,10000. ,
2 0.,4.,.001,
'F',4.,250.,100.,3.,0.,0./
C**** SET LOWER LIMITS OF PARAMETERS
DATA VALLO /5*0.,1. ,3
1 600. ,100.
*0.,-90.,-180.,1.,1.,0.,-50.,
,-50000., 2*0. ,'D' ,6*0./
c**** SET UPPER LIMITS OF PARAMETERS
DATA VALHI /I . ,1 . ,50000 . ,90 . ,90 . ,30 . ,2000 . , .050 , .0005 ,
1 90. ,180.,
12., 31., 99., 50., 1050., 100000., 50000.,
2 15., 1.,'P1, 15., 1000., 600., 5., 50., 150000. /
C**** NAME STABILITY CATEGORIES.
DATA ASTAB /'A' , 'B' , '
C','D' ,'E' ,'F'/
C**** WRITE OUT TABLE OF INPUT VALUES.
LU=6
10 WRITE(LU,1CO)
WRITE(LU,110) CAT(l),
WRITE(LU,120) CAT(2),
J6=6
NBOX=INT(VAL(6))
WRITE(LU,130) CAT(3),
WRITE(LU,140) CAT(4),
WRITE (LU, 150) CAT (5),
J12=12
J13=13
J14=14
MO=INT(VAL(12) )
IDA=INT(VAL(13))
IYR=INT(VAL(14»
(J,NAM(J) ,VAL(J) ,J=1,2)
(J,NAM(J) ,VAL(J) ,J=3,5)
J6,NAM(6) ,NBOX
(J,NAM(J) ,VAL(J) ,J=7,9)
(J,NAM(J) ,VAL(J) , J=10,ll)
WRITE(LU,160)CAT(6) ,J12,NAM(12) , MO, J13 ,NAM(13) , IDA, J14 ,NAM(14) , IYR
WRITE(LU,170) CAT(7)
WRITE(LU,180) CAT(8)
WRITE (LU, 190)
WRITE(LU,192) CAT(9),
WRITE(LU,194)
IF(LU.EQ.2) GO TO 70
(J,NAM(J) ,VAL(J) , J=15,16)
(J,NAM(J) ,VAL(J) ,J=17,19)
(J,NPM(J) ,VAL(J) ,J=20,22)
(J,NAM(J) ,VAL(J) , J=23,25)
(0,HAM(,1) ,VAL(J) ,J=26,27)
100 FORMAT(/,1H , ' THE PROGRAM INITIALIZATION PARAMETERS ARE SET TO THE
1 FOLLOWING VALUES:',/)
110 FORMAT(1H ,A20,2(2X,
120 FORMATflH ,A20,3(2X,
130 FORMATdH ,A20,2X,'(
140 FORMATdH ,A20,2X,'(
1 F6.3))
150 FORMATdH ,A20,2(2X,
160 FORMATdH ,A20,3(2X,
170 FORMATdH ,A20,2(2X,
180 FORMATdH ,A20,2X,'(
1 2X,'(
2 2X,'(
190 FORMATdH ,20X,2X,'(
1 2X,'(
2 2X,'(
192 FORMATdH ,A20,3(2X,
194 FORMATdH ,20X,2X,'(
1 2X,'(
(',12, ) ',1X,A7,F6.2))
(',12, ) MX,A7,F6.0))
,12,') ,1X,A7,I6)
,12,') ,1X,A7,F6.0,2(2X,'(',I2,') ',1X,A7,
(| ,12, ) |,1X,A7,F6.2))
\ t i2> t } fJ-XfA/fXa)/
(',12, ) ',1X,A7,F6.2))
,12,') ,1X,A7,F6.0,
,12,') ,1X,A7,F6.0,
,12,') ,1X,A7,F6.2)
,12,') ,1X,A7,F6.4,
,12,') ,1X,A7,A6,
,12,') ,1X,A7,F6.2)
(',12, ) ',1X,A7,F6.2)>
,12,') ,1X,A7,F6.2,
,12,') ,1X,A7,F6.0)
124
-------
124 C**** DETERMINE IF CHANGES TO INPUT VALUES ARE TO BE HADE.
125 WRITE(6,200)
126 200 FORMAT)/,1H ,'IF NO CHANGES ARE TO BE MADE IN VALUES IN TABLE',
127 1 ' ENTER 99',/,lH ,'IF "STAB" IS TO BE CHANGED ENTER1,
128 2 ' 98',/rlH ,'IF CHANGES OTHER THAN TO "STAB" ARE TO BE MADE',
129 3 ' ENTER TOT. NUMBER OF CHANGES')
130 READ(5,*) NU
131 IF(NU.EQ.99) GO TO 40
132 IF(NU.EQ.98) GO TO 30
133
134 C**** CHANGE VALUES IN TABLE.
135 WRITE(6,210)
136 210 FORMAT(/,1H ,'SPECIFY ID. NUMBER OF PARAMETER TO BE CHANGED1,
137 1 ' FOLLOWED BY THE NEW VALUE.',/,lH ,'ALL CHANGES ARE ENTERED',
138 2 ' SEPARATED BY COMMAS (E.G., 1,200.,2,50.,4,600., ETC).1)
139 READ(5,*) (ID( J) , VALNEW( J) , J=l, NU)
140
141 DO 20 J=1,NU
142
143 IF(VALNEW(J).LT.VALLO(ID(J)).OR.VALNEW(J).GT.VALHI(ID(J)))THEN
144 WRITE(6,240) ID(J)
145 240 FORMAT(/,1H ,'PARAMETER NO. ',12,' OUT OF RANGE, PLEASE',
146 1 ' RESPECIFY')
147 GO TO 20
148 ENDIF
149
150 VAL(ID(J))=VALNEW(J)
151
152 IF(ID(J).EQ.18)THEN
153 IF(ABS(VAL(18)>.GT.VAL(3)/2.)THEN
154 WRITE(6,250)
155 250 FORMATf/,1H ,'STACK MUST BE ON VALLEY FLOOR-RESPECIFY YZERO1)
156 ENDIF
157 ENDIF
158
159 IF(ID(J).EQ.22)THEH
160 IF(VAL(22).LE.1.1THEN
161 VAL(22)=1.
162 WRITE(6,260)
163 260 FORMAT(/,1H ,'DOWNVALLEY WINDSPEED AT STACK HAS BEEN SET',
164 + ' TO 1 M/S')
165 ENDIF
166 ENDIF
167
168 20 CONTINUE
169 GO TO 10
170
171 C**** CHANGE STABILITY DESIGNATION.
172 30 CONTINUE
173 WRITE(6,220)
174 220 FORMAT(/,1H ,'SPECIFY NEW STABILITY: D, E OR F ')
175 READ(5,230) VALI21)
176 230 FORMAT(A4)
177 GO TO 10
178
179 C**** SET VALUES FOR EACH PARAMETER.
180 40 CONTINUE
181 C SENSIBLE HEAT FLUX
182 AO=VAL(1)
183 K=VAL(2)
184 C VALLEY PHYSICAL CHARACTERISTICS
185 W=VAL(3)
186 ALPHA1=(VAL(4)/J60.)*6.28
187 ALPHA2=(VAL(5)/360.)*6.28
125
-------
188 C MODEL CHARACTERISTICS
189 NBOX=INT(VAL(6))
190 RNO=NBOX
191 C LENGTH OF BOXES.
192 BOXLEN=.5*W/RNO
193 C TIME STEP
194 DELTAT=((W/2.)/RNO)/10.
195 C INVERSION CHARACTERISTICS
196 HZERO=VAL(7)
197 GAMMA=VAL(8)
198 BETA=VAL(9)
199 C CROSS-SECTIONAL AREA OF VALLEY
200 AC= HZERO*W+(.5*(HZERO)**2)*(1-/TAN(ALPHAl)+
201 + l./TAN(ALPHA2))
202 C TOTAL NUMBER OF WHOLE BOXES IN MODEL AT INITIATION.
203 IFIX=(HZERO-HCBLI)/(BOXLEN*TAN(ALPHA1))
204 NTCTI=NBOX+IFIX
205 IF(NTOTI.GT.NB)THEN
206 WR1TE(6,280)
207 280 FORMAT(/,1H ,'TOO MANY GRID ELEMENTS— PLEASE REDUCE NBOX1
208 + ' OR RE-DIMENSION THE MODEL TO ALLOW MORE GRID ELEMENTS')
209 GO TO 10
210 ENDIF
211 C SITE LOCATION
212 LAT=VAL(10)
213 URLONG=VAL(11)
214 C DATE
215 MO»INT(VAL(12))
216 IDA=INT(VAL(13))
217 IYR=INT(VAL(14))
218 C VALLEY ATMOSPHERE
219 TEMP=VAL(15)
220 PRESS=VAL(16)
221 RHO=PRESS/(TEMP+273.16)/2.8704
222 C GAUSSIAN PLUME
223 XC=VAL(17)
224 YZERO=VAI,(18)
225 UOVAL(19)
226 Q=VAL(20)
227 DO 50 J=4,6
228 IF(VAL(21).EQ.ASTAB(J)) GO TO 60
229 50 CONTINUE
230 WRITE(6,290)
231 290 FORMAT(/,1H ,'YOU HAVE SPECIFIED A WRONG STABILITY CLASS-
232 1 F WILL BE USED INSTEAD')
233 J=6
234 60 ISTAB=J
235 C STACK CHARACTERISTICS
236 US=VAL(22)
237 H=VAL(23)
238 STMP=VAL(24)
239 SRAD=VAL(25)
240 SVEL=VAL(26)
241 AS=VAL(27)*1000.
242 LU=2
243 GO TO 10
244
245 c**** CONVERT DOWN-VALLEY DISTANCES TO KILOMETERS
246 70 XC=XC/1000.
247
126
-------
248 c**** CONVERT TEMPERATURES TO DEGREES
249 TEMP"=TEMP+273.16
250 STMP=STMP+273.16
251
252 RETURN
253 END
127
-------
2 SUBROUTINE JULIAN(MO,IDA,IYR,JULDAY)
3 C THIS SUBROUTINE CALCULATES THE JULIAN DATE GIVEN THE
4 C MONTH, DAY, AND YEAR.
5
6 DIMENSION NDAY(12)
7
8 DATA NDAY/0,31,59,90,120,151,181,212,243,273,304,334/
9
10 JULDAY=IDA+NDAY(MO)
11 C ADJUST FOR LEAP YEAP
12 A=FLOAT(IYR)/4-IYR/4
13 IF(A.EQ.O. .AND. MO.GE.3) JULDAY=JULDAY+1
14
15 RETURN
16 END
128
-------
1
2 SUBROUTINE PRISE(TEMP,SVEL,SPAD,STMP,XC,US,ISTAB,GAMMA,DLH)
3 c**** THIS SUBROUTINE IS USED TO CALCULATE PLUME RISE. IT IS
4 C**** DERIVED FROM THE PLUME RISE ALGORITHM IN THE E.P.A. MODEL
5 c**** "MPTER".
6
7 C**** DEFINITIONS **********************************************
8 C THE INPUT VARIABLES ARE: *
9 C GAMMA - POTENTIAL TEMPERATURE GRADIENT (DEC K/M). *
10 C ISTAB - P-G STABILITY (1,2,3,4,5 OR C). *
11 C SVEL - STACK EXIT VELOCITY (M/S). *
12 C SRAD - INSIDE STACK RADIUS (M). *
13 C STMP - STACK EXIT TEMPERATURE (DEC K). *
14 C TEMP - AMBIENT AIR TEMPERATURE (DEG K). *
15 C US - AMBIENT WIND SPEED (M/S). *
16 C XC - DISTANCE TO DOWNVALLEY CROSS SECTION (KM). *
17 C THE OUTPUT IS: *
18 C DLH - PLUME CENTER-LINE RISE ABOVE STACK HEIGHT (M). *
19 C***************************************************************
20
21 DATA G/9.8/
22
23 IF(SVEL.EQ.O..OR.SRAD.EQ.O.) THEN
24 DLH=0.
25 RETURN
26 ENDIF
27
28 X=XC*1000.
29 DELTT=STMP-TEMP
30 F«G*SVEL*SRAD*SRAD*DELTT/STMP
31
32 C **COMPUTE FOR NEUTRAL-UNSTABLE
33 IF(ISTAB.LE.4) THEN
34 IF(F.GE.55.) THEN
35 DTC=0.00575*STMP*SVEL**.6667/((2.*SRAD)**.3333)
36 ELSE
37 DTC=0.0297*STMF*SVEL**.3333/((2.*SRAD)**.6667)
38 ENDIF
39 IF(DELTT.LE.DTC.OR.STMP.LE.TEMP) THEN
40 C *MOMENTUM DOMINATED
41 DLH=6.*SRAD*SVEL/US
42 ELSE
43 C *BCUYANCY DOMINATED
44 IF(F.GE,55.) THEN
45 XF=3.5*34.*F**(.4)
46 ELSE
47 XF=3.5*14.*F**(.625)
48 ENDIF
49 XX=AMIN1(X,XF)
50 DLH=1.6*(F**.3333)*(XX**.6667)/US
51 ENDIF
52
53 C "COMPUTE FOR STABLE CONDITIONS
54 ELSE
55 S=G*GAMKA/TEMP
56 DTC=0.01958*TEMP*SVEL*SQRT(S)
57 IF(DELTT.LE.DTC.OR.STMP.LE.TEMP) THEN
58 C *MOMENTUM DOMINATED
59 DLH1=6.*SRAD*SVEL/US
60 DLH2=(1.5*((SVEL**2*SRAD**2*TEMP)/(STMP*US))**,3333)
61 + /(S**.1667)
62 DLH=AMIN1(DLH1,DLH2)
129
-------
63 ELSE
64 C *BOOYANCY DOMINATED
65 UCUT=((2.6/4.)**3)*(F**.25)*(S**.125)
66 IF(OS.GT.UCUT) THEN
67 XF=((2.6/1.6)**1.5)*US/SQRT(S)
68 ELSE
69 XF-((4./1.6)**1.5)*(US**1.5)/((S**(9./16.))*(F**.125)]
70 ENDIF
71 XX=AMIN1(X,XF)
72 DLH-1.6*(F**.3333)*(XX**.6667)/US
73 ENDIF
74 ENDIF
75
76 RETURN
77 END
130
-------
1 C+++++-H-+++++++++++++++++++-M
2 SUBROUTINE DILUTE(US,UC,AS,AC,DF)
3 C**** CALCULATE DILUTION FACTOR TO ACCOUNT FOR CLEAN AIR DILUTION OF THE
4 c**** NOCTURNAL POLLUTANT PLUME DURING ITS TRANSPORT DOWN THE VALLEY AS
5 c**** A RESULT OF CLEAN AIR FLUX INTO THE VALLEY FROM TRIBUTARIES, SLOPE
6 C**** FLOWS, AND/OR ENTRAINMEST AT THE TOP OF THE DOWN VALLEY FLOW
7 C**** LAYER.
8
9 IF(AS.EQ.O.)THEN
10 DF=1.
11 ELSE
12 DF=(UC*AC)/(US*AS)
13 IF(DF.LT.1.)DF=1.
14 ENDIF
15 RETURN
16 END
131
-------
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
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
DATA SIGGS/ 6.035, 2.1097,
SUBROUTINE INGRAT(XC,YZ ERO,ISTABY,ISTABZ,W,ALPHA1,
1 ALPHA2,HZERO,H,SUM1,SUM2,SIGMAY,SIGMAZ)
C**** CALCULATE GAUSSIAN DIFFUSION COEFFICIENTS AND DO A
C**** CROSSWIND INTEGRATION WITHIN THE VALLEY CROSS-SECTION
C**** TO CALCULATE THE MASS OF POLLUTION FLOWING THROUGH
C**** THE SECTION (KG/SEC). SIMILARLY, DO A CROSS-WIND INTEGRATION
C**** ABOVE THE VALLEY TO CALCULATE THE MASS OF POLLUTION FLOWING
C**** ABOVE THE VALLEY. THE REMAINING MASS (CALCULATED FROM THE
C**** EQUATION OF CONTINUITY AS Q MINUS THESE TWO QUANTITIES) CAN
C**** BE FOLDED BACK INTO THE CROSS-SECTION. THIS IS DONE IN THE
C**** MAIN PROGRAM UNDER AN INSTANTANEOUS PERFECT-MIXING
C**** ASSUMPTION.
REAL SIGGY(18),SIGGZ(18)
DATA SIGGY/ 5.357, 0.8828, -0.0076, 5.058, 0.9024, -0.0096,
1 4.651, 0.9181, -0.0076, 4.230, 0.9222, -0.0087,
3.922, 0.9222, -0.0064, 3.533, 0.9181, -0.0070/
6.035, 2.1097, 0.2770, 4.694, 1.0629, 0.0136,
1 4.110, 0.9201, -0.0020, 3.414, 0.7371, -0.0316,
2 3.057, 0.6794, -0.0450, 2.621, 0.6564, -0.0540/
C**** DEFINITIONS *************************************************
C H PLUME CENTERLINE HEIGHT (M) *
C HZERO INITIAL INVERSION HEIGHT (M) *
C ISTABY HORIZONTAL STABILITY INDEX (1-6) *
C ISTABZ VERTICAL STABILITY INDEX (1-6) *
C SIGMAY STD DEVIATION OF PLUME CONC IN Y DIRECTION (M) *
C SIGMAZ STD DEVIATION OF PLUME CONC IN Z DIRECTION (M) *
C ALPHA1 VALLEY SIDEWALL SLOPE (RAD) *
C ALPHA2 OTHER VALLEY SIDEWALL SLOPE (RAD) *
C W FULL WIDTH OF VALLEY FLOOR (M) *
C XC DOWNVALLEY DISTANCE OF CROSS SECTION (KM) *
C YZERO OFF-CENTERLINE DISPLACEMENT OF STACK (M) *
C Y,Z COORDINATE AXES (M) *
C ZINC VERTICAL HEIGHT INCREMENT FOP INTEGRATION (M) *
£******************************************************************
C**** CALCULATE THE SIGMAS USING MCMULLEN'S METHOD (.1975) **********
SIGMAY-EXP(SIGGY(3*ISTABY-2)+SIGGY(3*ISTABY-l)*LOG(XC)+
1 SIGGY(3*ISTABY)*(LOG(XC))**2)
SIGMAZ = EXP(SIGGZ(3 *ISTABZ-2)+SIGGZ(3*ISTABZ-1)*LOG(XC) +
1 SIGGZ(3*ISTABZ)*(LOG(XC))**2)
C**** INITIAL VALUES
PI=3.14159
P=100.
NP=P
ZINC=HZERO/P
SUM1=0.
C**** RATE OF MASS FLOW ACROSS THE VALLEY SECTION
DO 1 1=1,NP
R=I
Z=(R-1.)*ZINC
F=(1./(SQRT(2.*PI)*SIGMAZ))*EXP(-.5*((Z-H)/SIGMAZ)**2)
Yl=(W/2.+Z/TAN(ALPHAl)-YZERO)/SIGMAY
Y2=(-W/2.-Z/TAN(ALPHA2)-YZ ERO)/SIGMAY
132
-------
60
61 IF(Y1.LT.O.)THEN
62 Y1=-Y1
63 CALL NORMAL(Y1.PHIY1)
64 PHIY1-1.-PHIY1
65 ELSE
66 CALL NORMAL(Y1,PHIY1)
67 ENDIF
68
69 IF(Y2.LT.O.JTHEN
70 Y2=-Y2
71 CALL NORMAL(Y2,PHIY2)
72 PHIY2=1.-PHIY2
73 ELSE
74 CALL NORMAL(Y2,PHIY2)
75 ENDIF
76
77 PHIY=PHIY1-PHIY2
78 SUM1=SUM1+PHIY*F*ZINC
79 1 CONTINUE
80
81 C**** RATE OF MASS FLOW ABOVE THE CROSS SECTION.
82 SUM2=0.
83 DO 2 1=1,NP
84 R=I
85 Z=(R-1.)*ZINC*2.+HZERO
86 ARG=-.5*((Z-H)/SIGMAZ)**2
87 IFfARG.LT.-80.JTHEN
88 F=0.
89 ELSE
90 F=(1./(SQRT(2.*PI)*SIGMAZ))*EXP(ARG)
91 ENDIF
92 Yl=5.
93 Y2=-5.
94
95 IF(Yl.LT.O.JTHEN
96 Y1=-Y1
97 CALL NORMAL)Yl,PHIYl)
98 PHIY1=1.-PHIY1
99 ELSE
100 CALL NORMAL(Y1,PHIY1)
101 ENDIF
102
103 IF(Y2.LT.O.JTHEN
104 Y2=-Y2
105 CALL NORMAL(Y2,PHIY2)
106 PHIY2=1.-PHIY2
107 ELSE
108 CALL NORMAL(Y2,PHIY2)
109 ENDIF
110
111 PHIY=PHIY1-PHIY2
112 SUM2=SUM2+PHIY*F*ZINC*2.
113 2 CONTINUE
114
115 RETURN
116 END
133
-------
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
SUBROUTINE NORMALfX,PHI)
C**** CALCULATE THE INTEGRAL OF THE AREA UNDER THE GAUSSIAN CURVE
C**** FROM MINUS INFINITY TO X. THE POLYNOMIAL APPROXIMATION
C**** USED HERE COMES FROM ABRAMOWITZ AND STEGUN, 1965,
C**** HANDBOOK OF MATHEMATICAL FUNCTIONS, DOVER PRESS, NINTH
C**** PRINTING, P. 932. THEY ATTRIBUTE THE FORMULA TO HASTINGS (1975)
REAL C(6)
DATA C/0.2316419, 0.31938153, -0.356563782,
1 1.781477937, -1.821255978, 1.330274429/
PI=3.14159
T-1./(1.+C(1)*X)
SIGH-(l./SQRT(2.*PI))*EXP(-X**2/2.)
PHI=1.-SIGH*(C(2)*T+C(3)*T**2+C(4)*T**3+C(5)*T**4+C(6)*T**5)
RETURN
END
134
-------
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
SUBROUTINE GAUSS ( Y, Z,YZERO,H,Q, OS, DF,SIGMAY,SIGMAZ , CHI)
C**** THIS PROGRAM USES THE WELL-KNOWN GAUSSIAN PLUME EQUATIONS
C**** (TURNER, 1969) TO CALCULATE VALLEY POLLUTION CONCENTRATION
C**** DOWNVALLEY OF A POLLUTANT SOURCE.
DEFINITIONS *************************************************
CHI POLLUTANT CONCENTRATION (MICROGRAMS/M**3) *
CHIOQ CHI OVER Q (S/M**3) *
CHIOQN NORMALIZED AND STP CORRECTED .CONCENTRATION (S/M**3)*
DF DILUTION FACTOR *
H PLDME CENTERLINE HEIGHT (M) *
Q SOURCE STRENGTH (KG/S) *
SIGMAY SIGMA Y (M) *
SIGMAZ SIGMA Z (M) *
US WINDSPEED AT SOURCE (M/SEC) *
YZERO OFF-CENTERLINE DISPLACEMENT OF STACK (M) *
XC,Y,Z COORDINATES OF RECEPTOR (KMFM,M) *
C
C
C
C
C
C
C
C
C
C
C
C**** CALCULATE INDIVIDUAL TERMS IN THE GAUSSIAN PLUME EQUATION ***
C Gl : OFF CENTERLINE IN Y DIRECTION *
C G4 : OFF CENTERLINE IN Z DIRECTION *
£******************************************************************
C**** GAUSSIAN PLUME EQUATION
G1=EXP(-0.5*((Y-YZERO)/SIGMAY)**2)
G4=EXP(-0.5*((Z-H)/SIGMAZ)**2)
CHIOQ=1./(2.*3.14*SIGMAY*SIGMAZ*US*DF)*(G1)*(G4)
C**** CONCENTRATION IN MICROGRAMS PER CUBIC METER
CHI-CHIOQ*Q*1.E09
RETURN
END
135
-------
2 SUBROUTINE SOLAR(LAT,URLONG,JULDAY,NTSR,TAU,Al)
3 C**** CALCULATE TIME OF SUNRISE, LENGTH OF DAY, AND SOLAR FLUX ON
4 c**** HORIZONTAL SURFACE AT SOLAR NOON FOR ANY SITE GIVEN LATITUDE,
5 C**** LONGITUDE, AND JULIAN DATE. THIS SUBROUTINE USES KCCULLOUGH'S 1968
6 C**** APPROXIMATIONS AS GIVEN IN ARCH. METEOR. GEOPHYS. BIOKLIMAT.,
7 C**** SER. B, 16, 129-143. REFER TO WHITEMAN, 1980, ATMOSPHERIC
8 C**** SCIENCE PAPER NO. 328, COLORADO STATE UNIVERSITY, FORT
9 C**** COLLINS, COLORADO, 250 PP.
10
11 REAL EQNTIM(25)
12 REAL LAT,LONG,LONCOR
13
14 C**** EQUATION OF TIME IN HRS FAST OR SLOW. .HHHH VALUES GIVEN FOR
15 c**** EVERY 15 JULIAN DAYS, FOR EXAMPLE, 1,16,31,46 ..360.
16 DATA EQNTIM/-0.0533,-0.1589,-0.2233,-0.2378,-0.2064,-0.1442,
17 1 -0.0689,+0.0003,+0.0475,+0.0622,+0.0425,-0.0028,
18 2 -0.0558,-0.0961,-0.1058,-0.0789,-0.0186,+0.0639,
19 3 +0.1517,+0.2258,+0.2683,+0.2647,+0.2094,+0.1094,
20 4 -0.0131/
21
22 C**** CONSTANTS
23 C JULIAN DATE OF VERNAL EQUINOX.
24 DZERO=80.
25 PI-3.14159
26 C CONVERSION-DEGREES TO RADIANS.
27 CONV=PI/180.
28 C ECCENTRICITY OF EARTHS ORBIT
29 ECCENT=.0167
30 C MAX DECLINATION.
31 DECMAX=(23.+27./60.)*PI/180.
32 C REVOLUTION RATE OF EARTH.
33 OMEGA=2.*PI/365.
34 C SOLAR CONSTANT IN WATTS/M**2
35 SC=1367.
36 C 15 DEGREES OF LONGITUDE IS EQUIVALENT TO 1 HOUR.
37 ONEHR=15.*PI/180.
38 C LONGITUDE OF STANDARD MERIDIAN FOR YOUR TIME ZONE.
39 STDLON=105.
40 LAT=LAT*CONV
41
42 C JULIAN DATE
43 D=JULDAY
44
45 c**** MAKE COMPUTATIONS.
46 OMD=OHEGA*D
47 OMDZRO=OMEGA*DZERO
48 C RATIO OF RADIUS VECTORS SQUARED. ACCOUNTS FOR VARYING DISTANCE
49 C BETWEEN EARTH AND SUN.
50 RDVCSQ=1./(1.-ECCENT*COS(OMD))**2
51 LONG=OMEGA*(D-DZERO)+2.*ECCENT*(SIN(OMD)-SIN(OMDZRO))
52 DECLIN=A SIN(SIN(DECMAX)*SIN(LONG))
53 SR=-ABS(ACOS(-TAN(LAT)*TAN(DECLIN)))
54 C LENGTH OF DAYLIGHT PERIOD
55 TAU=-SR*2./ONEHR
56
57 c**** CALCULATE TIME OF SOLAR NOON AT SITE AS A PRELIMINARY TO CALCULATING
58 c**** THE LOCAL TIME FOR ANY GIVEN HOUR ANGLE.
59 LONCOR=(STOLON-URLONG)*1./15.
60 IF(D.GT.361.) GO TO 1
61 ID=((D-1.)/15.)+1.
62 D2=(ID-1)*15+1
63 TIMCOR=(EQNTIM(ID+1)-EQNTIM(ID))*(D-D2)/15.+EQNTIM(ID)
136
-------
64 1 IF(D.EQ.362.) TIMCOR=-0.0211
65 IF(D.EQ.363.) TIMCOR=-0.0292
66 IF(D.EQ.364.) TIMCOR=-0.0372
67 IF(D.EQ.365.) TIMCOR=-0.0453
68 TMNOON-12.00-LONCOR-TIMCOR
69
70 C**** CALCULATE SOLAR FLUX AT SOLAR NOON
71 COSZ=SIN(LAT)*SIN(DECLIN)+COS(LAT)*COS(DECLIN)*1
72 C SOLAR FLUX AT SOLAR NOON
73 A1»SC*RDVCSQ*COSZ
74 IF(A1.LT.0.) Al=0.0
75 T1-(12.+(SR/ONEHR)+(TMNOON-12.00))
76 IT=T1
77 T2 = IT
78 C TIME OF SUNRISE.
79 NTSR-(T1-T2)*60.+(T2*100.)
80
81 RETURN
82 END
137
-------
1
2 SUBROUTINE EBDGT(AO,A2)
3 c**** DETERMINE FRACTION OF EXTRATERRESTRIAL SOLAR FLUX THAT PRODUCES
4 C**** SURFACE SENSIBLE HEAT FLUX.
5 c**** AT PRESENT, THIS SUBROUTINE IS NOT BEING USED, SINCE THE USER NOW
6 c**** SIMPLY INPUTS FRACTION AO. IN FUTURE, THIS SUBROUTINE WOULD BE
7 c**** DEVELOPED TO EXPLICITLY INCORPORATE FEATURES OF THE SURFACE ENERGY
8 C**** BUDGET SUCH AS ALBEDO, SOIL MOISTURE AND CLOUD COVER.
9
10 A2=AO
11
12 RETURN
13 END
138
-------
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
39
j t\
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
C++-f++++-H"t-+++++++++ +++++++++++++++++-f +++++++++++ +++++++++++++++++
SUBROUTINE DBS CUT (AO , Al , ALPHA1 , ALPHA2 , BETA, DELTAT, GAMMA,
1 HCBLI,HZERO,K,W, PRESS, RHO, TAU, NSTEPS)
C**** THIS PROGRAM CALCULATES THE HEIGHTS OF THE INVERSION TOP
C**** AND CBL TOP AS A FUNCTION OF TIME AFTER SUNRISE AS VALLEY
C**** INVERSIONS ARE DESTROYED, USING WHITEMAN'S (1980)
C**** NUMERICAL METHOD AS GIVEN IN HIS EQUATIONS 41-50.
PARAMETER (NS=6000)
COMMON/BLK1/HC(NS) ,HT(NS)
REAL K
C AO FRACTION (0-1) OF Al CONVERTED TO SENSIBLE HEAT
C Al SOLAR FLUX ON HORIZ SFC AT SOLAR NOON (W/M2)
C ALPHA1 SIDEWALL NO. 1 INCLINATION ANGLE (RAD)
C ALPHA2 SIDEWALL NO. 2 INCLINATION ANGLE (RAD)
C BETA RATE OF TEMP CHANGE AT INVERSION TOP (K/SEC)
C CP SPECIFIC HEAT AT CONSTANT PRESSURE (J/KG/K)
C DELTAT TIME STEP (SEC)
C GAMMA VERT POT TEMPERATURE GRADIENT IN INVERSION (K/M)
C HC CBL HEIGHT (M)
C HCBL HEIGHT OF CBL TOP (M)
C HCBLI INITIAL CBL HEIGHT (M)
C HT INVERSION TOP HEIGHT (M)
C HTOP HEIGHT OF INVERSION TOP (M)
C HZERO INITIAL HEIGHT OF INVERSION TOP (M)
•++++•
FLUX
*
*
*
Ik
*
*
*
*
*
It
*
*
*
C K FRACTION (0-1) OF SENS HT FLUX GOING INTO CBL GROWTH
C N TIMESTEP COUNTER
C NSTEPS NUMBER OF TIME STEPS RQD TO DESTROY INVERSION
C PRESS AVG PRESSURE (MB) AT INV CENTER AT SUNRISE
C PI TRIGONOMETRIC CONSTANT
C RHO AIR DENSITY (KG/M3)
C TAU LENGTH OF DAYLIGHT PERIOD (SEC). 43200SEC=12HR
C T ELAPSED TIME SINCE SUNRISE (SEC)
C TT TIME OF SUNRISE (SEC)
*
*
it
*
*
*
*
*
C W VALLEY FLOOR WIDTH (M) *
°
C**** INPUT PARAMETERS
CP=1005.
PI=3.14
N=0
TI=0.
T=0.
C**** CALCULATE SECONDARY PARAMETERS
FACT1=AO*A1/(RHO*CP)
THTOT=(1000./PRESS)**.286
C=(1./TAN(ALPHA1))+(1./TAN(ALPHA2))
c**** D0 CALCULATIONS
HTOP= HZERO
HCBL=HCBLI
GO TO 2
1 CONTINUE
T=T+ DELTAT
2 CONTINUE
N=N+1
FACT2=FACT1*SIN(PI*(T-TI)/TAU)
139
-------
64 IF(K.NE.O.) THEN
65 DHDT=FACT2*THTOT*K*(W+HCBL*C)/
66 1 (GAMMA*HCBL*(W+(HCBL*C)/2.))
67 ELSE
68 DHDT=0.
69 ENDIF
70
71 FNUM=(W+HTOP*C-K*(W+HCBL*C))*FACT2-((BETA*
72 l(HZERO-HTOP)*(W+( (HZERO+HTOP)/2.)*C) )/{2.*THTOT) )
73 FDENOM=HTOP*GAMMA*(W+(HTOP*C)/2.)
74 l+(BETA*(T-TI)*(W+HTOP*C))/2.
75 DHTDT=-THTOT*FNUM/FDENOM
76 DHT=DHTDT*DELTAT
77 HTOP=HTOP+DHT
78 HT(N)=HTOP
79 DHC=DHDT*DELTAT
80 HCBL=HCBL+DHC
81 HC(N)=HCBL
82
83 c**** STOP CALCULATIONS WHEN INVERSION TOP SINKS BELOW CBL TOP
84 IF(HTOP.LE.HCBL) THEN
85 NSTEPS=N
86 RETURN
87 ELSE
88 GO TO 1
89 ENDIF
90
91 END
140
-------
1
2 SUBROUTINE PROFIL(NBOX,NTOTI,NTOT,BOXLEN,HCL,ALPHA1,
3 1ALPHA2,HZERO,H,HCU,HTL,HTU,YZ ERO,Q,OS,DF,
4 2SIGMAY,SIGMAZ,CHIOFF)
5 C**** THIS SUBROUTINE DETERMINES THE AVERAGE POLLUTANT CONCENTRATION
g c**** INJECTED INTO THE TOP OF EACH BOX.
7
8 PARAMETER (NB=30)
9
10 COMMON/BLK2/CHIBAR(NB)
11
12 c**** DEFINITIONS ***************************************************
13 C HCL CBL HEIGHT AT LOWER TIME STEP (M)
14 C HCU CBL HEIGHT AT UPPER TIME STEP (M)
15 C HTL INVERSION TOP HEIGHT AT LOWER TIME STEP (M)
16 C HTU INVERSION TOP HEIGHT AT UPPER TIME STEP (M)
17 C Yl Y COORD AT LEFT SIDE OF SINKING MASS ELEMENT (M)
18 C Y2 Y COORD AT RIGHT SIDE OF SINKING MASS ELEMENT (M)
19 C Zl EFFECTIVE Z COORD AT BOTTOM OP SINKING MASS ELEMENT (M)
20 C Z2 EFFECTIVE Z COORD AT TOP OF SINKING MASS ELEMENT (M)
2i £********************************************************************
22
23 ALPHA=(ALPHAl+ALPHA2)/2.
24
25 C**** CALCULATE CHIBAR FOR THE BOXES ON THE VALLEY FLOOR.
26 DO 1 I=1,NTOTI
27 Y1=(I-1)*BOXLEN
28 Y2-I*BOXLEN
29 Y»(Yl+Y2)/2.
30
31 IF(I.GT.NTOT) THEN
32 CHIBAR(I)-999.9
33 GO TO 1
34 ELSE IFd.GT.NBOX .AND. I.LE.NTOT) THEN
35 Zl=HCL-l"HZERO-nTL-l-(FLOAT(I-(NBOX+l))+.5)*BOXLEN*TAN(ALPHA)
36 Z2=HCU+HZERO-HTU+(FLOAT(I-(NBOX+1))+ .5)*BOXLEN*TAN(ALPHA)
37 ELSE
38 Z1=HZERO-HTL+HCL
39 Z2=HZERO-HTU+HCU
40 ENDIF
41
42 Z=(Zl+Z2)/2.
43 CALL GAUSS(Y,Z1,YZERO,H,Q,US,DF,SIGMAY,SIGMAZ,CHIZ1)
44 CALL GAUSS(Y,Z2,YZERO,H,Q,US,DF,SIGMAY,SIGMAZ,CHIZ2)
45 CALL GAUSS(Y1,Z,YZERO,H,Q,OS,DF,SIGMAY,SIGMAZ,CHIY1)
46 CALL GAUSS(Y2,Z,YZERO,H,Q,US,DF,SIGMAY,SIGMAZ,CHIY2)
47 CHIBAR(I)=(CHIZH-CHIZ2+CHIYl+CHIY2)/4. + CHIOFF
48
49 1 CONTINUE
50 RETURN
51 END
141
-------
1
2 SUBROUTINE VELOCY(N,MN,NTOTI,NBOX,NTOT,BOXLEN,DELTAT)
3 c**** THIS SUBROUTINE CALCULATES VELOCITY AT EACH SEPTUM BASED ON
4 C**** MASS CONTINUITY.
5
6 PARAMETER (NS=6000,NB1=31)
7
8 COMKON/BLK1/HC(NS),HT(NS)
9 COMMON/BLK3/V(NB1)
10
11 V(l)=0.
12 CIH=HT(MN)-HT(N)
13 IF(NTOT.EQ. NBOX. AND. HT(N) .LE.HC(N)) GO TO 3
14
15 DO 2 I=1,NTOTI
16 IF(I.GT. KTOT) GO TO 1
17 V( I-H) = (I *BOXLEN*CIH/DELTAT) /HC (N)
18 GO TO 2
19 1 V(I+1)=999.9
20 2 CONTINUE
21 GO TO 5
22
23 3 DO 4 I=1,NTOTI
24 V(I+1)=999.9
25 4 CONTINUE
26
27 5 RETURN
28 END
142
-------
1
2 SUBROUTINE BRKUP(N,UN,NTOTI,BOXLEN,DELTAT,CCH)
3 c**** THIS SUBROUTINE CALCULATES POLLUTANT CONCENTRATIONS IN EACH
4 c**»* op THE GRID ELEMENTS (BOXES) USING A POLLUTANT MASS BALANCE.
5
6 PARAMETER (NS=6000,NB=30,NB1=J1)
7
8 REAL CCH(2,NB)
9
10 COMHON/BLK1/HC(MS),HT(NS)
11 COHMON/BLK2/CHIBAR(NB)
12 COMMON/BLK3/V(NB1)
13
14 C*DEFINITIONS **********************************************************
15 C A IS THE POLLUTANT MASS WITHIN THE BOX AT THE PREVIOUS TIME STEP.
16 C B IS THE MASS COMING FROM THE ADJACENT DOWNHILL BOX.
17 C C IS THE MASS COMING INTO TOP OF BOX DUE TO SINKING INVERSION
18 C AND/OR GROWING CBL
19 C D IS THE MASS GOING OUT UPHILL SIDE OF THE BOX.
20 C E IS THE MASS SOURCE OR SINK FROM SURFACE OF BOX. *
21 c************************************************************************
22
23 E=0.0
24 DELTAH=HC(N)-HC(MN)+HT(MN)-HT(N)
25 DO 4 1=1,NTOTI
26 IF(V(I+1).EQ.999.9) GO TO 3
27 A=CCH(1,I)*BOXLEN*HC(MN)
28 IF(I.EQ.l) GO TO 1
29 B=D
30 GO TO 2
31 1 B=0.0
32 2 C=CHIBAR(I)*DELTAH*BOXLEN
33 D=(DELTAT/BOXLEN*V(I+1)+(HC(N)-HC(MN))/HC(N))*A
34 1 -HDELTAT/BOXLEN*V(I+1))*(B+C)
35 CCH(2,I)=(A+B+C+E-D)/(BOXLEN*HC(N))
36 GO TO 4
37 C**** EXPONENTIAL DECAY
38 3 S=ALOG(.90)
39 CCH(2,I)=CCH(1,I)*EXP(S*DELTAT/10.)
40 4 CONTINUE
41
42 RETURN
43 END
143
-------
1
2 SUBROUTINE PSTPRC{TSTP,NTOTI,NAS)
3 c**** THIS SUBROUTINE PROCESSES THE TIME SERIES OF CALCULATED
4 c**** 5-MINUTE (300 S) AVERAGE CONCENTRATIONS TO OBTAIN MAXIMUM
5 c**** SHORT TERM (1 AND 3 HOUR) AVERAGE CONCENTRATIONS WHICH CAN
g c**** BE COMPARED TO REGULATORY STANDARDS.
7
8 PARAMETER (NA=100,NB=30,NB1=31)
9
10 REAL SUMM(NB,2)
11
12 COMMON/BLK4/CONC(NA,NB)
13 COMMON/BLK5/CHI(NB)
14 COMMON/BLK6/AVG(NB,2),NDXTIM(NB,2),NTS(2)
15
16 c**** TIME IN SECONDS
17 ONEHR-3600.
18 THREE«10800.
19
20 c**** NUMBER OF TIMESTEPS
21 NTS(1)=OHEHR/TSTP
22 NTS(2)=THREE/TSTP
23 NSTART=NAS+1
24
25 c**** CALCULATE MAX 1-HR AND 3-HR AVERAGE CONCENTRATIONS
26 DO 7 K=l,2
27 TS=NTS(K)
28 DO 7 1=1,NTOTI
29 SUM=0.
30
31 c**** AVERAGING INTERVAL LESS THAN LENGTH OF DAYTIME CONC SERIES.
32 IFfNSTART.GT.NTS(K)) THEN
33 DO 1 N=HAS,NSTART-NTS(K),-1
34 SUM=SUM+CONC(N,I)
35 1 CONTINUE
36 SUMM(I,K)=SDM
37 Q1=SUM
38 NTIMES=NAS
39 DO 2 N=NAS,NTS(K)-H,-1
40 SUMM(I,K)=SUMM(I,K)+CONC(N-NTS(K),I)-CONC(N,I)
41 IF(SUMM(I,K).GT.Ql)THEN
42 Q1=SOMM(I,K)
43 NTIMES=N-1
44 ENDIF
45 2 CONTINUE
46 DO 3 N=NTS
-------
64 Q1=SOMM(I,K)
65 NTI«ES=N-1
66 ENDIF
67 6 CONTINUE
68 ENDIF
69
70 NDXTIM(I,K)=NTIHES
71 AVG(I,K)=Q1/TS
72 7 CONTINUE
73
74 HETURN
75 END
145
-------
APPENDIX B
FORTRAN LISTING OF VALMET OUTPUT PLOTTING PROGRAM
147
-------
1 C PROGRAM PLOTT.FOR
2
3 C THIS PR06RAM IS DESIGNED TO PRODUCE PLOTS OF CONCENTRATION
4 C VERSUS TIME FOR ARBITRARY 6RID ELEMENTS IN MODEL VALMET.
5 C THIS PROSRAfl ALSO PRODUCES A SINGLE PA6E CONTAININ6 THREE OTHER PLOTS:
6 C CONC VS HEIBHT, CONC VS X-VALLEY DISTANCE, AND HEI6HT VS TIME.
7
3 PARAMETER ,CONCEN(NA),)C80RD(5),
11 1 YBORD(5),NBOX(10),l(MAn(10),YMAn!10)
12 DIMENSION X1(4),K2(4),Y1(4),Y2(4),ICNT!3,2)
13 1 ,T(3,4!,Z(1000),TIC(3,2},
14 2 YU000),CHIZUOOO>,CHIY(1000>,HC(NS),HT(NS),
15 3 TIME!NS),SZ1(3},SZ2(3),TIK(3,2)
16
17 CHARACTER120 INFILE
18 CHARACTERS ITIME
19 CHARACTERI2 ICHARI30)
20 CHARACTER126 ILBL(3,2)
21 CHARACTER}! ANSi,ANS2
22 CHARACTERI2 JLBL1
23 CHARACTERJ5 JLBL2
24 CHARACTERI5 JLBL3
25
26 DATA Xl/3.0,2.5,5.5,5.5/
27 DATA X2/10.0,4.0,10.0,10.0/
28 DATA Yl/3.0,5.0,5.0,2.0/
29 DATA Y2/7.5,8.0,8.0,3.5/
30
31 DATA ICHAR/'Or,'02','03','04','05','06','07','08',
32 1 '09','10','11','12','13','14','15','16',
33 2 '17','18','19','20','21','22','23','24',
34 3 '25','26','27','28','29','30'/
35
36 SZH1K35
37 SZ1(2)=.35
38 SZ1(3)=.70
39 SZ2(1)=.40
40 SZ2(2)=.40
41 SZ2(3)=.80
42
43 1CNT(1,1)=12
44 ICNT(2,1)=10
45 ICNT(3,l)=2i
46 ICNT(1,2)=14
47 ICNT!2,2!=14
48 ICNT!3,2)=12
49
50 ILBLU,1!='CONC (ug/»3)'
51 ILBL!2,l)='TIt1E (LOT)'
52 ILBL(3,1)='CROSS-VALLEY DISTANCE !k«)'
53 ILBL(1,2)='HEIGHT !• A6L)'
149
-------
54 ILBL(2,2)='H£IBHT !a A6L)'
55 ILBL(3,2!='CQNC (ug/§3)'
54
57 TIC(l,l)=i.O
58 TIC!1,2)=200.
59 TIC(2,1)=1.0
60 TIC(2,2)=200.
61 TIC(3,l)=i.O
62 TIC<3,2)=1.0
63
44 TIK(1,1)=0.5
65 TIKU,2)=100.
66 TIK(2,1)=1.0
67 TIK(2,2)=iOO.
68 TIK(3,1)=0.5
69 TIK(3,2)=0.5
70
7! HR!TE!6,1QO)
72 100 FORHATtlH ,'UHAT IS THE NAHE OF THE DATA FILE?')
73 READ(5,'(A20)')INFILE
74 QPEN
-------
107 501 FORMAT(Al)
108 IF(ANS2.NE.'Y') 60 TO 502
109 WR!TE(6,1Q3>
110 103 FORMAT(1H ,'THE FIRST PLOT IS CONCENTRATION VS TIME')
111 XMAX=TTIME(1)+NINDXtDELTAT/3600.
112 HRITE!6,104) NTSR,XHAX,CLCONC
113 104 FORMATUH ,'NIN AND MAX VALUES ARE: TIME ',I5,2X,
114 1F8.2,5K,'CONCENTRATION 0.0',2X,F4.2)
115 WRITE(6,105)
Hi 105 FORMATdH ,'ENTER T«IN,TMAX,CONCMIN,CONCMAX')
117 READ(5,I> XMIN,XMAX,YHIN,YHAX
118 XBQRD!1)=XHIN
119 XBQRD(2)=XNAX
120 IBORDttMHAX
121 XBORD(4>=XNIN
122 XBORD(5)=X«IN
123 YBORD(1)=YMIN
124 YBORD(2)=YMIN
125 YBORD(3)=YMAX
126 YBORD!4)=YMAX
127 YBORD(5)=YMIH
128 WRITE(6,106) HTOTI
129 106 FQRMATilH ,'THERE ARE ',12,' GRID ELEMENTS')
130 WRITE (6,107)
131 107 FQRHATUH ,'HO« MANY CURVES DO YOU HAMT TO PLOT?')
132 READ(5,t) NC
133 WRITE!6,108)
134 108 FORHATIJH ,'CURVES ARE TO BE DRAWN FOR WHICH 6RID ELEMENT
135 1 NUMBERS? EX:3,6,9')
136 READ(5,I) (NBOX(N),N=1,NC)
137
138 Cttll DRAW THE PLOTS
139 CALL KTERH(IMODE)
140 CALL KBE6IN
141 CALL KPA6E(X1(1),X2(1),Y1(1),Y2(1),0)
142 CALL KSCALE(XMIN,XMAX,0,YMIN,YHAX,0)
143 CALL KAXIS(XMIN,YMIN,XMAX,1.0,0.5,0,0.0,1,0.35,0.,
144 1ILBLI2,!),10,1,0.40,0.)
145 CALL KAXISUMIN,YNIN,YHAX,0.2,0.1,1,0., 1,0.35,0.,
146 1ILBL(3,2),12,1,0.40,0.)
147 CALL KCURVEC(BDRD,Y80RD,5)
148 IMAXI(N)=0.
149 YMAXI(N)=0.
150 DO 3 N=1,NC
151 DO 2 MN=1,NAS+1
152 CONCEN(MN)=CONC(«N,NBQX!N»
153 IF'MN.EQ.l) SO TO 14
154 IF(CONCEN(«N).ST.CONCEN(MN-D) THEN
155 XMAXI!N)=TTIME!MN)
156 YMAXI!N)=CONCEN(MN)
157 ENDIF
!53 14 CONTINUE
15" 2 CONTINUE
151
-------
160 HRITE(6,109) XHAXI(N),YMAXI!N)
161 109 FORMAT UH ,2F8.2)
162 CALL KLINE(N,0,0.,0.)
163 CALL KCURVE(TTIME,CONCEN,NAS+i)
164 3 CONTINUE
165 DO 4 N=1,NC
166 YNAXI(N}=YNAmN}+.10l
167 HRITE«6,HO) Y«AXI(N)
168 110 FORMAT(F8.2)
169 4 CONTINUE
170 DO 5 N=1,NC
171 !=N80X(N)
172 »RITE(6,111) I.ICHARd)
173 111 FORMAT(1H ,I3,A2)
174 CALL KUALFA(XHAXI(N),Y?1AXI(N),ICHAR(I),
175 12,1,0.35,0.,0.)
176 5 CONTINUE
177 IX=X
178 XA=XH!N+.74
179 XB=XMIN+1.45
180 XC=XHIN+l.i5
181 XPOS1=XNIN+0.4
182 YPOSl=.9tYf!AX
183 YPOS2=.8»YMAK
184 YPQS3=.7tYHAX
185 HTL=0.4
186 HRITE(JLBL1,300) IX
187 300 FORHAT(I2)
138 WITE(JLBL2,30i) CLCONC
139 301 FOR«AT(F5.2)
190 URITE
-------
213 TINE(HN)=TTI«E(1)
214 ELSE IF!!1N.E9.2)THEN
215 TIME(nN)=TTIKE(2!
216 ELSE
2i7 TINE(HN)=TTIME(2)HMN-2)tDELTAT/3600.
218 EHDIF
219 i CONTINUE
220
221 C WRITE
-------
266 CALL K8E6IN
267 DO 12 H=l,3
268 CALL KPASEai(Htl),X2(«tl),YHH+l),Y2«N+l)tO)
269 CALL KSCALE(T(H,l),T(n,2),0,T(l1,3),T!N,4},0)
270 mjRD(l)=T(N,l)
271 XBOR»(2)=T(«,2)
272 JBORD(3)=T(«,2)
273 HBORD(4)=T(H,1)
274 XBQRD(5)=T
230 CALL KLINE(1,0,0.0,0.0)
281 CALL KCURVEUBORD,YBORD,5)
282 CALL KAXIS,0.)
286 IF(n.EQ.l) THEN
287 CALL KLINEI4,0,0.0,0.0)
288 CALL KCURVE!CHIZ,Z,K)
289 ELSE IF(H.E9.2) THEN
290 CALL KLINEI4,0,0.0,0.0)
291 CALL KCURVE(TIHE,HC,NSTEPS)
292 CALL KLINE(2,Q,0.0,0.0)
293 CALL KCURVE(TIHE,HT,NSTEPS)
294 ELSE IFIN.EQ.3) THEN
295 CALL KLINE!4,0,0.0,0.0)
296 CALL KCURVE!Y,CHIY,J>
297 ELSE
298 STOP 4
299 ENDIF
300 12 CONTINUE
301 13 CONTINUE
302 STOP 1
303 END
154
-------
APPENDIX C
RESEARCH PAPER ENTITLED "BREAKUP OF TEMPERATURE
INVERSIONS IN DEEP MOUNTAIN VALLEYS:
PART II. THERMODYNAMIC MODEL"
155
-------
Reprinted from JOURNAL OF APPLIED METEOROLOGY, Vol. 21, No. 3, March 1982
Amman Mcuorotofical Society
Pnmod m U. S. A.
Breakup of Temperature Inversions in Deep Mountain Valleys:
Part II. Thermodynamic Model
C. DAVID WHITEMAN
THOMAS B. McKEE
157
-------
290
JOURNAL OF APPLIED METEOROLOGY
VOLUME 21
Breakup of Temperature Inversions in Deep Mountain Valleys:
Part II. Thermodynamic Model
C. DAVID WHITEMAN
Pacific Northwest Laboratory, Richland, WA 99352
THOMAS B. McKEE
Department of Atmospheric Science. Colorado State University. Fort Collins 80523
(Manuscript received 31 March 1981, in final form 7 December 1981)
ABSTRACT
A thermodynamic model is developed to simulate the evolution of vertical temperature structure during
the breakup of nocturnal temperature inversions in mountain valleys. The primary inputs to the model are
the valley floor width, sidewail inclination angles, characteristics of the valley inversion at sunrise, and an
estimate of sensible heat flux obtained from solar radiation calculations. The outputs, obtained by a numerical
integration of the model equations, are the time-dependent height of a convective boundary layer that grows
upward from the valley floor after sunrise, the height of the inversion top, and vertical potential temperature
profiles of the valley atmosphere. The model can simulate the three patterns of temperature structure
evolution observed in deep valleys of western Colorado. The well-known inversion breakup over flat terrain
is a special case of the model, for which valley floor width becomes infinite. The characteristics of the model
equations are investigated for several limiting conditions using the topography of a reference valley and
typical inversion and solar radiation characteristics. The model is applied to simulate observations of inversion
breakup taken in Colorado's Eagle and Yampa Valleys in different seasons. Simulations are obtained by
fitting two constants in the model, relating to the surface energy budget and energy partitioning, to the
data. The model accurately simulates the evolution of vertical potential temperature profiles and predicts
the time of inversion destruction.
1. Introduction
In Part I (Whiteman, 1982) observations of ver-
tical temperature structure and wind evolution were
summarized for 21 case studies of nocturnal tem-
perature inversion breakup in deep Colorado moun-
tain valleys. Breakup occurred following one of three
patterns of vertical temperature structure evolution.
A hypothesis was offered to explain the observations
in which sensible heat flux was the driving force and
one of the three patterns occurred, depending on
whether the sensible heat flux was used primarily to
cause convective boundary layers (CBLs) to grow
over the valley floor and sidewalls (Pattern 1), to
remove mass from the inversion in the upslope flows
(Pattern 2), or to accomplish both (Pattern 3). In
Part II, a bulk thermodynamic model of temperature
inversion destruction is developed based on the hy-
pothesis, using several simplifying assumptions.
2. Mathematical model of inversion destruction
a. General equations
Two approaches can be used to develop a math-
ematical model able to simulate temperature changes
in the valley atmosphere. In the first approach de-
tailed mathematical equations can be developed for
the individual components of the overall system, in-
cluding the various boundary layers and stable core.
However, the valley atmosphere consists of many
interrelated layers, and the coupling of the equations
for the different layers to simulate potential tem-
perature changes in the valley atmosphere as a whole
would be difficult due to geometrical considerations
and lack of detailed information on physical char-
acteristics of the various layers. Consequently, a sec-
ond approach is taken in which a bulk thermody-
namic model is developed for the valley inversion.
As more is learned about the individual components
of the system, the model can be refined to introduce
greater detail into the simulation of the individual
components.
The thermodynamic model of valley temperature
structure evolution developed here is based on the
hypothesis of Part I. Fig. 1 shows a unit thick cross
section of a mountain valley having sidewalls of in-
clination or, and a2 and valley floor width /. At sunrise
(t,) the valley is assumed to have an inversion of
depth hi and constant vertical potential temperature
gradient 7. The variable width of the valley at the
0021-8952/82/030290-I3S07 25
9 1982 American Meteorological Society
159
-------
MARCH 1982
C. DAVID WHITEMAN AND THOMAS B. McKEE
291
INVERSION TOP AT t«ti
H ••
FIG. 1. Valley geometry and potential temperature profiles used to formulate a mathematical model of inversion destruction.
top of the inversion is designated by L, and the origin
of a y-z coordinate system is placed at the center of
the valley floor at point 0. After sunrise a typical
Pattern 3 potential temperature evolution ensues, in
which a CBL develops over the valley floor and
sidewalk. Removal of mass from the valley in these
CBLs allows the stable core to sink so that a sounding
taken at an arbitrary later time r will show a lower
inversion top height h(t ) and a shallow CBL of height
H(t) near the ground. A later sounding, taken at a
time t0 when the inversion has just been destroyed,
will show a neutral atmosphere having a potential
temperature 8 = 9h.
From the first law of thermodynamics, the incre-
ment of energy required to increase the potential
temperature of a mass of air m by the potential tem-
perature increment A# is
(1)
where T/9 = (P/IOOQ)*'' = 1, p is the density (as-
sumed constant) and V is volume. Using (1), the
energy required to change the valley potential tem-
perature profile at time t, to the profile at time t can
be obtained by an integration over the valley volume
below the height of the inversion top. The total en-
ergy requirement is composed of two pans: the en-
ergy increment g2 that removes mass from the valley
and allows the top of the inversion to sink, and the
energy increment Q3 that causes a CBL to grow.
These energies are represented by the areas desig-
nated in Fig. 1 and are given by
T r c*i f?* f>
- MJxdydz
9 [_Jo JV£ Jo
f (•>•<< /-i
- &6tdxdydz
Jo Jxt. Jo
where
(2)
(3)
(4)
yL = - - H
v2 tana2/'
Vt> = 1" .
y 2 tana,
~ l . + . l
tana, tana2
= f(h - Z),
and
J fH «-« (M
- Mydxdydz,
0 JO Jyi JO
where
y(H - z).
(5)
(6)
(7)
(8)
(9)
(10)
(11)
In order to simplify the integration, it was assumed
that the valley temperature structure is horizontally
homogeneous across the valley, that mass is removed
from the valley in such a way that the potential tem-
perature gradient in the inversion layer does not
change with time, and that p and cf are constant. By
differentiating the individual energies Q2 and Q3 with
respect to time, the rates of change of the height of
the top of the inversion and the height of the CBL
are obtained, such that
dt
dt
. dH,
HC\
(13)
The total rate of energy input into the valley to
accomplish these changes is the fraction A0 of solar
irradiance F coming across the area L of the top of
the inversion that is converted to sensible heat. The
solar irradiance may be approximated by a sine func-
tion having a certain amplitude A\ and period T, so
that the total rate of energy input becomes
dO,
-^ = AoLF = A0(l
at
hC)A, sin - (/ - t,). (14)
160
-------
292
JOURNAL OF APPLIED METEOROLOGY
VOLUME 21
An energy balance for the valley inversion is obtained
by equating (14) to the sum of (12) and (13). Al-
ternatively, a fraction of the energy input is available
to dnve the growth of the CBL while the rest of the
incoming energy is used to remove mass from the
valley. The fraction of energy input used to drive the
CBL growth is assumed to be of the form
where k is a number between 0 and 1. This form is
chosen in order to simplify later equations. Equating
this fraction of the energy input to (12) and the re-
mainder to (13) and solving for dH/dt and dh/dt
results in the final model equations
dt
dt
1JL
Tpc.
hC~
1 +
HCY1
—• sinT- (r - r()] . (16)
T« LT J
These equations specify the dependence of the rate
of ascent of the CBL and the rate of descent of the
inversion top on inversion characteristics, incoming
energy and valley topography. An integration of the
coupled equations allows the simulation of the time-
dependent behavior of the heights of the CBL and
inversion top. If the potential temperature #/, at the
top of the inversion is known and is independent of
time, and 7 is constant, knowledge of the variation
of h and H with time is sufficient to specify how
vertical profiles of potential temperature change with
time. When k = 0, the equations provide an ap-
proximate simulation of Pattern 2 inversion destruc-
tion in which destruction occurs solely due to the
removal of mass from a valley in the slope flows,
resulting in a descent of the inversion top. When k
= 1, the equations provide a simulation of Pattern
1 inversion destruction, in which destruction occurs
mainly due to the growth of a CBL over the valley
floor. When k = 1 and the valley floor becomes very
wide, the simulation approaches that of inversion
destruction over the plains. When k is between 0 and
1, the equations provide a simulation of Pattern 3
inversion destruction in which the inversion is de-
stroyed by the combined effect of a growing CBL
and a descending inversion top. A more complete
description of the characteristics of the model equa-
tions for Pattern 1, 2 and 3 temperature structure
evolution follows.
b. Pattern 2 inversion destruction
The physical hypothesis of Pattern 2 inversion de-
struction requires a shallow CBL to form over the
sidewalls so that additional energy can be used to
cause mass to flow up them. Pattern 2 destruction
may be approximated by assuming that all the energy
available to destroy the inversion goes solely to move
mass up the sidewalls, causing the top of the inversion
to descend. This can be accomplished by setting k
equal to zero in (15) and (16) resulting in the two
equations
f-a
dh
dt
'I + hC
pc.yh
sin -(t~t,)\.
LT J
(18)
Following these equations, the CBL does not grow
as a function of time, and the inversion is destroyed
as the top of the inversion sinks. The rate of descent
of the top of the inversion increases as the inversion
descends. The descent rate is faster when more en-
ergy is available and when the potential temperature
gradient of the inversion is weaker. The factor in
large parentheses in ( 1 8 ) is a topographic factor that
varies from 1 to 2 depending on the shape of the
valley cross section. This accounts for the reduced
volume of air within the mountain valley relative to
that over the plains for the same energy flux on a
horizontal surface. Since the valley has less volume
to be heated by the same incoming energy, it warms
more rapidly. By separating variables h and / and
integrating from the initial conditions h = h, and H
= 0 at / = r, to h - h and H = H at r = t, analytical
expressions are obtained which describe how H and
h vary with time, i.e.,
H=Q,
pc,y ir
(19)
(20)
Fig. 2 illustrates the shapes of the curves of h vs
t for a reference simulation and for two cases where
one parameter in the reference simulation is changed.
The reference simulation uses representative values
of valley parameters observed in western Colorado
including / = 1000 m, a, = a2 = 15°, h, = 500 m,
7 = 0.025 K m-', r = 12 h = 43200 s, T/B = 1 and
p = 1 kg m~!. The reference inversion takes nearly
4'/> h to be destroyed when A^A\/pc,, = 0.25 K m s"'.
Fig. 2 was obtained by a numerical integration of
(18) using a forward finite difference scheme and a
time step of 10 min.
An analytical expression for the time required to
161
-------
MARCH 1982
C. DAVID WHITEMAN AND THOMAS B. McKEE
500 r
293
400
300
200
100
REFERENCE SIMULATION
012345
t-tj (hours)
FIG. 2. Descent of inversion top as a function of time for the reference inversion
simulation and for two simulations for which the single parameters indicated were
changed. Pattern 2 destruction.
destroy an inversion can be obtained by integrating
(18) from the initial conditions to the final conditions
of h = 0 at i = to. This expression,
tD — t, = - cos"
IT
8
«.»
enumerates the factors affecting inversion breakup
time. The sensitivity of inversion breakup time (mea-
sured from sunrise) to the various parameters is il-
lustrated in Fig. 3 using the reference simulation
above. The reference inversion takes 4.4 h to break
(vertical line in Fig. 3). The effect on the time re-
quired to destroy the reference inversion by varying
the individual parameters is obtained by following
the labeled curves. Thus, varying the initial height
of the inversion from 400 to 600 m, other parameters
being equal, changes the time required to destroy the
inversion from 3.5 to 5.4 h. For the reference inver-
sion, the most sensitive parameters affecting the time
required to break an inversion are the available en-
ergy and the initial potential temperature gradient
X
04O
030
020
010
h;
-800
SOO
(m)
4OO
200
3 4 ~S
tD-tj (hours)
FIG. 3. Sensitivity of inversion destruction time to various model parameters for Pattern 2 destruction.
162
-------
294
JOURNAL OF APPLIED METEOROLOGY
VOLUME 21
500
t-t|( hours)
FIG. 4. Ascent of CBL and descent of inversion top as a function
of time for Pattern 1 inversion destruction in a valley for the
reference simulation and for two simulations in which single pa-
rameters were changed to the values indicated.
and inversion height. In the normally dry Colorado
valleys the most important factors affecting the avail-
able energy are albedo (snow versus no snow) and
latent heat flux. The effect of valley shape on the
breakup time is relatively small for normal ranges
of a and / encountered in valleys of western Colorado.
Nevertheless, it is apparent that the valley width and
sidewall angles may affect the mode of inversion de-
struction since they control, to a certain extent, the
divergence of mass in the CBL's and thus determine
whether inversion destruction more nearly follows
Pattern 1 or Pattern 2. Overall, predictions of the
time required to destroy inversions, given typical val-
ues of the parameters observed in field experiments,
are consistent with the observed range of 3.5-5 h.
c. Pattern I inversion destruction—Valley case
A useful approximation to Pattern 1 inversion de-
struction can be obtained from the model equations
by setting k = \ in (15) and (16). The general equa-
tions then reduce to the two equations
dH=9_ f I + HC
dt ~ T\l
dh = _9_ \~(h-H)C~\
dt~ Tll
AoA, . fir .
— —• sin -(r-
pCpyh ]_r
, (22)
, .
• (23)
Eq. (22) is formulated so that the entire fraction
AQ of the energy coming across the area (/ -t- HC~)
of the top of the CBL is used to cause the CBL to
grow. The energy used to cause the top of the in-
version to descend is the fraction A<> of the difference
between the energy coming across the top of the in-
version and the energy coming across the top of the
CBL. The time-dependent behavior of the height of
the CBL can be obtained by an integration of (22)
from the initial condition of H = 0 at t = t, to the
final condition of H = H at / » t, such that
(24)
The integration of (22) and (23) can be accomplished
numerically to determine how h and H change with
time in a Pattern 1 inversion destruction in a moun-
tain valley. This is done for the reference simulation
and for two simulations in which a single parameter
of the reference simulation is changed. The resulting
plots are shown in Fig. 4. The characteristics of the
plots include a near-linear growth of the CBL with
time, a slow descent of the inversion top, and a more
rapid breakup than for Pattern 2 destruction. Pattern
1 destruction takes 3.7 h versus the 4.4 h required
for Pattern 2 breakup. The same total amount of
energy is required to destroy the reference inversion,
whether it is destroyed following Pattern 1 or Pattern
2. However, since the energy available to destroy the
inversion comes across the area of the top of the
inversion, and this area is larger when the inversion
top sinks more slowly, the total amount of energy
required to destroy the inversion is attained earlier
in the day, resulting in an earlier inversion breakup.
d. Pattern I inversion destruction—Flat plains case
Application of (22) and (23) to a valley that is
very wide or approaches a plain (/ — oo) results in
the equations
dH 8
dt
(25)
dt
(26)
Thus, over flat terrain where no topographically in-
duced mass divergence occurs from the CBL, the
inversion is destroyed solely by the growth of a CBL.
Performing integrations on (25) as described in the
previous section results in the analytical expression
J2£l4>dl i-cos-U-r,) , (27)
TT —
and an expression for the breakup time,
, (28)
163
-------
MARCH 1982 C. DAVID WHITEMAN AND THOMAS B. McKEE
500 r
REFERENCE SIMULATION x
500
400
300
I
200
100
295
01 234567
t-t, (hours)
FIG. 5. Growth of CBL over Sat terrain as a function of time for the modified reference
inversion (a " 0, / —• vo) and for two simulations in which single parameters were changed
to the values indicated. Pattern ! destruction.
Eq. (27) is nearly identical to an equation for CBL
growth over homogeneous terrain as derived by
Leahey and Friend (1971). Following (27), Fig. 5
presents plots of H vs time for the reference simu-
lation and for two cases where one parameter in the
reference simulation has been changed. The refer-
ence inversion is destroyed in 5.65 h by a near-linear
increase in the depth of the CBL, If the incoming
energy is doubled, the inversion takes 3.8 h to break,
and if the potential temperature gradient is increased
to 0.035 K. m"', the inversion is broken in ~7.1 h.
Fig. 6 indicates the sensitivity of the time required
to break an inversion on the different parameters of
(28). The effect on the breakup time of changing
individual parameters in the reference simulation is
obtained by following the individual curves.
e. Pattern 3 inversion destruction
A simulation of Pattern 3 inversion destruction
uses the general model [ Eqs. (15) and (16) ], in which
a partitioning of energy is required to allow both
CBL growth and inversion top descent. In order to
use the general model equations, the fraction of sen-
•340
456
t0-ti (hours)
30
(•Km/sec.)
1 20
-lo
FIG. 6. Sensitivity of inversion destruction time to various model parameters
for Pattern 1 inversion destruction over flat terrain.
164
-------
296
JOURNAL OF APPLIED METEOROLOGY
VOLUME 21
500
400
300
r
o
aoo •
100 •
012345
t-f, (hours)
Fie. 7 Ascent of CBL and descent of inversion top as a function
of tune for Pattern 3 destruction of the reference inversion for
different values of k.
sible heat flux K = k[(l + HC)/(l + he)] that drives
the growth of the CBL must be determined. It is
apparent that the fraction K is a function of time,
since the initial energy input must be used primarily
to develop the CBL's before appreciable mass can
be carried up them. Factor K also depends on the
topographic characteristics of the valley, since K
must approach 1 as the valley width approaches in-
finity. It seems probable that K may also be a func-
tion of sensible heat flux. Since the functional de-
pendencies of K are not yet known, it is assumed that
k is a constant and, by comparing model simulations
to actual data, the constant value of k that results
in the best fit to data is determined. This approach
allows an investigation of the effect of k on the sim-
ulation. Further research is necessary to determine
the actual functional form of K, so that a better un-
derstanding of the energy partitioning phenomenon
can be obtained, resulting in more accurate simu-
lations.
The equations used to simulate Pattern 3 destruc-
tion are thus
dH
dt ''
8 k f l + HC
X rinjj (t - t,) \ , (29)
^L = _ 1 _L f/ + hC-k(l +
J ~fh
x sinTj (f - r,)l, (30)
0 < (k = constant) « 1. (31)
These equations can be integrated numerically to
determine how H and h vary with time from a given
where
initial state. Fig. 7 shows several numerical integra-
tions using this method, time steps of 10 min, the
reference simulation, and several values of k. Re-
plotted on the same figure are some of the limiting
cases of the general equations for destruction of the
reference inversion discussed earlier. The time re-
quired to break the reference inversion decreases
from 4.4 to 3.7 h as the value of k is increased from
0 to 1. For all values of k the growth of the CBL is
nearly linear.
Fig. 8 shows the effect of varying individual pa-
rameters in the reference simulation with k fixed.
The inversion is destroyed when the ascending CBL
and the descending inversion top meet at a height
of HD = h0 = 205 m. The fact that the fraction k
uniquely determines the height at which the ascend-
ing CBL meets the descending inversion top at the
time of inversion destruction, suggests a means of
fitting the model results to actual data that will be
used in a later section.
/ Model modification to account for warming of the
neutral layer
Field observations show that the potential tem-
perature 8), at the top of the inversion usually in-
creases slowly with time. The warming rate varies
from valley to valley and from day to day with the
average warming rate being about 0.4 K h"1. From
the physical hypothesis, this warming requires that
more energy be spent to move mass up the sidewalls,
since the parcels must be warmed to a higher tem-
perature 9h(t) = 8k(t,) -I- (d8h/dt)6t to be removed
from the valley. The energy requirement can be cal-
500
400
300
200
100
035
REFERENCE SIMULA!,ON
01 2345
f-t, (hours)
FIG. 8. Ascent of CBL and descent of inversion lop for Pattern
3 destruction (k = 0.2) of the reference inversion, and for two
simulations in which single parameters were changed to the values
indicated.
165
-------
MARCH 1982
C. DAVID WHITEMAN AND THOMAS B. McKEE
297
culated by considering that the mass of air removed
from the valley to allow the top of the inversion to
sink from h, to h is given by
dQ±
—r-
T Q\~
c!,--\(
mass removed
pV
p(h, - h)(l
~(t- O (/ +
. (35)
Vi(A, + h)C], (32)
and that the energy required to move this mass across
a potential temperature jump at the top of the in-
version, 0»(r) — #*('() can be approximated by
(33)
If the wanning in the neutral layer occurs linearly
in time, the potential temperature jump is given by
j3(/ — r,) where ft = d8hjdt is the rate of wanning.
Then
X [/ + H(h, + h)C], (34)
The model equations are modified to account for this
extra energy requirement by specifying that the frac-
tion k((l + HC)/(l + hC)] of the inversion energy
input dQt/dt drives the CBL growth, or
*(
I + HC
\l + hCj dt
dQ2
dt '
(36)
while the remainder of the energy input drives the
descent of the inversion and carries parcels across
the potential temperature jump, such that
hCj dt
dt
dt
(37)
Substituting Eqs. ( 1 2 ), ( 1 3 ), ( 1 4) and ( 3 5 ), the mod-
ified model equations are
dt
/
dh^ _^J_
dt ' Tpcf
sin -
r
-O-pc,-^(Af-A){/
1 sinj - (t - t,]
h)C] ]
I
- f,)(/ + AC)
(38)
(39)
These equations can be used to simulate inversion
breakup when significant wanning occurs in the neu-
tral layer above. Eq. (39) reduces to (30) when /3 is
zero. It is important to note that (38) has the same
form as before but, since more energy is required to
move mass up the sidewalls, the partitioning of en-
ergy may be affected. If this actually occurs in na-
ture, then the functional dependency of K is more
complicated than previously discussed, since it will
depend not only on time, energy input and valley
width, but also on the rate of warming of the neutral
layer above the inversion. As before, Eqs. (38) and
(39) may be integrated numerically to simulate val-
ley inversion breakup if Bh(t) is known. If 9* does not
vary with time, Eqs. (29) and (30) can be used for
the simulation.
3. Comparison of model results with data
In this section a finite difference form of the model
equations (38) and (39) will be applied to simulate
actual data collected in the valleys of western
Colorado. The input parameters needed to solve the
equations are given in Table 1 along with a summary
of how they may be obtained. The output of the
model is h(t) and H(t). From these outputs and the
assumptions that 7 is constant for altitudes between
H and h, that potential temperature is independent
of height in the CBL and neutral layer, and that
B*(t,) and ddhjdt are known, 8(t, z) for the CBL.
stable core and neutral layer can be determined.
Unfortunately, since two of the input parameters to
the thermodynamic model (k and A0) were not ob-
served in the field programs, the equations cannot
be applied directly. Instead, arbitrary values for k
and A0 are chosen until the best simulation of the
data is obtained with the model. It is then determined
whether the values of k and A0 are reasonable for
the situation at hand. Both k and A0 are bounded,
since they are fractions between 0 and 1. The value
of k specifies the constant fraction of sensible heat
flux that is used to cause the CBL to deepen. The
model results are presented below for a winter Pat-
tern 2 inversion destruction in the wide Yampa Val-
ley and for a fall Pattern 3 inversion destruction in
the Eagie Valley.
a. Pattern 2 simulation—Yampa Valley, 23 Feb-
ruary 1978
The model input parameters required to simulate
the Pattern 2 inversion destruction observed in the
snow covered Yampa Valley on 23 February 1978
were obtained as follows. First, 9T~{ = 1.07 and
pc? = 1040 J m"J K"1 were calculated using the ap-
166
-------
298
JOURNAL OF APPLIED METEOROLOGY
TABLE 1. Model input parameters.
VOLUME 21
Model input
External conditions (neutral layer warming)
Energy partition
Surface energy balance
Numencal
(viii)
(ix)
(x)
(xi)
(xii)
(xiii)
(xiv)
Source
Constants
Valley topography
Initial inversion characteristics
Solar irradiance
(i)
(ii)
(iii)
(iv)
(v)
(vi)
(vii)
8 IT
PC,
I
C
1
h,
^
from average P and T of
from topographic maps
from sunrise sounding
from extraterrestrial solar
valley atmosphere
i [radiance model
or field notes for t,
from sequential soundings taken during
observation of inversion breakup or from
climatic data
comparisons of theory and data
measurements, if available
arbitrary
proximate mean pressure (780 mb) of the early
morning inversion and an average temperature of the
valley atmosphere (—10°C) from the defining equa-
tion for potential temperature and from the equation
of state, respectively. Second, the valley topographic
parameters (/ = 2580 m, a, = 9°, a2 = 16°, C
= 9.80) were obtained from topographic maps.
Third, the initial inversion parameters (h, = 530 m,
1 = 0.0345 K m"1) were estimated from a straight
line fit to the top of the 0714 MST sounding on this
date. The hyperbolic lower region of the sounding
could not be adequately fit with a straight line, so
it is ignored in the analysis. Fourth, a sinusoidal fit
to the extraterrestrial solar flux curve obtained from
a standard solar irradiance model (e.g., Sellers,
1965) provides solar irradiance parameters (At
= 878 W m'2, t, = 0655 MST, r = 10.9 h). The
output of the solar irradiance model for the hori-
zontal surface of interest is presented in Fig. 9. For
reference, the solar fluxes on extraterrestrial surfaces
with the same aspect and inclination angles as the
valley sidewalls are indicated on the figure. The input
required to run the solar irradiance model includes
the date and the latitude and longitude of the Yampa
Valley site. The remaining parameters necessary to
run the inversion destruction model are fractions k
and A& and the neutral layer wanning rate 0. To
simulate a Pattern 2 destruction, k must be zero. The
neutral layer warmed only 1.1 K during the inversion
destruction, for a warming rate $ of 2.8 X 10~5 K
' SOM8HCRO RANCH
YAMPA VALLEY
23 FEBRUARY 1978
II 13
TIME (MST)
17
19
FIG. 9. Extraterrestrial solar flux calculated for valley floor and sidewall surfaces
of the Yampa Valley, 23 February 1978.
167
-------
MARCH 1982
C DAVID WHITEMAN AND THOMAS B. McKEE
299
SCO
500
J 400
t-
| 300
I 200
100
0
YAMPA VALLEY NEAR STEAMBOAT SPRINGS, COLORADO
FEBRUARY Z5, 1978
THEORETICAL SUNRISE 0655
SUNRISE ON VALLEY FLOOR 0802
0 DATA POINTS
m MODEL
0700
0800
0900
1200
1300
1400
1500
1000 1100
TIME (MST)
FIG. 10. Comparison of model simulation -of h(t) with actual data for the Yarapa Valley, 23 February 1978.
s~'. The best fit of the model output h(t) to the in-
version top data of Fig. 10 was obtained with A9
= 0.19 by a trial-and-error procedure. Using this
value, the simulation of the height of the top of the
inversion agrees to within 25 m of the data over a
7 h period. In mid-afternoon, however, the simulation
continues to call for the descent of the inversion top,
when the actual data indicate that the descent
stopped. Thts is probably due to afternoon shading
of the valley by a mountain peak southwest of the
site. The model, using a simple sine function to sim-
ulate solar flux, does not account for this shading.
In Fig. 11 the potential temperature profile simula-
tions are compared to actual sounding data. The
sounding data consist of eight consecutive potential
temperature soundings taken at •—1 h intervals
throughout the day. The characteristics of the initial
inversion were obtained from the 0714 MST sound-
ing and are indicated by the straight line fit to the
sounding in the figure. The slow neutral layer warm-
ing is also apparent in the figure. The excellent fit
of the model simulations to the data in the stable
core above the 50 m deep surface layer is shown for
soundings 4 and 6. The value of A0 required to obtain
these results seems reasonable for the snow-covered
valley with evergreen forests covering much of the
valley sidewall above the observation site. It is par-
ticularly interesting that the model assumption of
constant potential temperature gradient was satisfied
well during the long period of inversion destruction
in the complicated topography of the Yampa Valley.
The case study is an excellent example of the effect
of enhanced albedo due to snow cover in retarding
the normal breakup of an inversion. Despite the fact
that the temperature inversion was not destroyed on
this clear day, the diurnal range of temperature at
the ground was quite large.
b. Pattern 3 simulation—Eagle Valley, 16 October
1977
Pattern 3 inversion breakup in the Eagle Valley
is remarkablv consistent in all seasons when snow
cover is not present in the valley. To test the math-
ematical model, the Pattern 3 breakup on the clear
day of 16 October 1977 is chosen for simulation. The
input parameters to the model are obtained as for
the Pattern 2 simulation above. Thus, values of the
constants 07*"' (1.08) and pcf (990 J m'3 K~') are
determined from the approximate mean pressure
(768 mb) and temperature (0°C) of the valley in-
version. The valley topographic parameters (/ = 1450
m, a, = 21°, «2 = 10°, C = 8.28) are determined
from topographic maps. The initial values of the in-
version parameters (h, = 650 m and 7 = 0.0269 K
m~') were taken from the 0650-0719 MST sounding
since the pre-sunrise sounding was of insufficient
height to determine h,. A solar irradiance model was
used to determine the parameters A\ = 906 W irT2,
t, = 0621 MST and r = 11 h (Fig. 12). The value
of /3 (8.3 X 10"5 K s"1) was taken from valley tem-
perature sounding data.
The fitting of the model output to the data on Fig.
13 was accomplished first by choosing the value of
|- YAMPA VALLEY
FEBRUARY 23, 1978
300 -
200 -
iooL
SIMULATION
ze:
270 275 23O 285 290 295
POTENTIAL TEMPERATURE ( K!
Fio. 11. Comparison of model simulation of potential temper-
ature structure with data for the Yampa Valley, 23 February 1978
Upsoundings 1-3 were initiated at 0714, 0905. 0959, 1100, 1202.
1259, 1359 and 1508 MST, respectively.
168
-------
300
JOURNAL OF APPLIED METEOROLOGY
VOLUME 21
u.
IT
Id
CC
o:
uj
I4OO
1200
IOOO
Z~ aoo
6
J 500
400
200
STEVE MILLER RESIDENCE
EASLE VALLEY
. 16 OCTOBER 1977
jO'SLOPE. 194-ASPECT
S-FACIN6 SLOPE
0600 0800
IOOO 1200 .1400
TIME (MST)
Fic. 1 2. Extraterrestrial solar flux calculated for valley floor and sidewall
surfaces of Eagle Valley, 16 October 1977.
k so that the ascending CBL and descending inver-
sion top met at the proper observed height at the
time of inversion destruction. The value of A<> was
then varied until the model outputs h(t) and ff(t)
fit the data well. The fit in Fig. 13 was obtained with
k = 0.14 and A0 » 0.45. The value of A0 seems
realistic considering the dry nature of the valley sur-
face on this date. The equations were integrated from
the initial conditions of h = 650 m and H = 10 m
at i = 0705 MST using a time step of 10 min. From
Fig. 13 it is clear that a good fit to the data is obtained
with the chosen values of the two parameters Aa and
k. The simulation of CBL height and inversion top
height agrees with the data within 50 m over most
of the period of inversion destruction. The CBL
height was overpredicted early in the inversion pe-
riod, due to the bulk nature of the model. Following
sunrise, the CBL develops first over the illuminated
sidewall or sidewalls, and somewhat later over the
valley floor when it is sunlit. The bulk model, how-
ever, does not differentiate between CBL growth over
the three different valley surfaces. All energy going
into CBL growth is attributed, in the model, to
growth of the valley floor CBL. Thus, the initial
700i
600^
400
200
100
STEVE MILLER RESIDENCE
EAGLE VALLEY
OCTOBER 16,1977
0700
0800 09OO
TIME (MST)
IOOO
overprediction of CBL growth over the valley floor
is a characteristic feature of the model equations.
The behavior of the simulation at mid-levels of the
valley atmosphere near the time of inversion destruc-
tion should also be mentioned. The data typically
show a more sudden inversion breakup is to be ex-
pected in nature since the final remnants of the stable
core will break up in convective overturning, once
the stable core becomes thin enough and the con-
vective plumes rising from the valley floor become
vigorous enough. Due to the chaotic nature of the
breakup, actual soundings taken during this time will
often show deformations in the vertical potential
temperature profiles. Estimates of heights h and H
from soundings are difficult to make during this time.
A continued flux of energy into this region is nec-
EAGLE VALLEY
FIG. 13 Comparison of model simulation of H(t) and hd)
with actual data for the Eagle Valley, 16 October 1977
285 290 293 300 305 310
POTENTIAL TEMPERATURE («K)
FIG. 14. Comparison of model simulation of potential temper-
ature structure with data from the Eagle Valley, 16 October
1977
169
-------
MARCH 1982
C. DAVID WHITEMAN AND THOMAS B. McKEE
301
essary before the deformations in the profiles are
destroyed and well-organized convection through the
entire depth of the valley results in smooth neutral
profiles. When this occurs, the breakup can be con-
sidered as finished.
The potential temperature profiles corresponding
to the data and simulation of Fig. 13 are presented
in Fig. 14. The potential temperature data are from
the tethersonde up-soundings taken at 0650-0719,
0829-0848, 0924-0948, 1013-1033 and 1106-1121
MST. The corresponding potential temperature sim-
ulations are plotted for the approximate midtimes
of the tethersonde up-soundings at 0840, 0935, 1025
and 1115. Again, the model provides a good simu-
lation of the actual data, reproducing the wanning
and growth of the CBL, the warming of the stable
core, and the descent of the inversion top.
4. Summary and conclusions
A thermodynamic model of valley temperature
structure evolution has been developed. By differ-
entiating the first law of thermodynamics, the rate
of energy input into the valley atmosphere is equated
to the rate of descent of the inversion top and the
rate of ascent of the CBL under the following as-
sumptions:
I) The potential temperature gradient in the sta-
ble core is constant during the period of inversion
breakup.
2) The temperature structure is horizontally ho-
mogeneous in the cross-valley direction.
3) The valley topography can be adequately rep-
resented by a horizontal valley floor of arbitrary
width and two linear sidewalls of arbitrary slope.
4) The initial inversion at sunrise can be ade-
quately represented by a constant potential temper-
ature gradient layer of arbitrary depth.
5) Convective boundary layers can be adequately
represented as constant potential temperature layers
of arbitrary height.
6) The rate of energy input into the valley at-
mosphere is a constant fraction of the solar energy
flux (assumed to be a sinusoidal function of time
from sunrise to sunset) coming across the horizontal
upper surface of the inversion.
To complete the model, the partitioning of the
energy input into CBL growth and mass transport
must be estimated. As a first approximation, it is
assumed that a constant fraction of the energy input
is used to cause the CBLs to grow. The remaining
energy is used for mass transport in the upslope flows.
Energy going into the growth of the CBL causes the
CBL depth to increase with time. Energy causing air
parcels to flow up the sidewall CBLs results in the
descent of the top of the inversion. When most of the
available energy drives the growth of the CBLs, a
temperature structure evolves in which the inversion
is destroyed predominantly by the upward growth
of a CBL from the ground. When most of the avail-
able energy drives mass up the sidewalls, a temper-
ature structure evolves in which the inversion is de-
stroyed by the descent of the inversion top.
The inputs to the model are 1) the initial inversion
characteristics (depth and average potential temper-
ature gradient); 2) valley topography (floor width
and inclination angles of the two sidewalls); 3) the
rate of energy input; and 4) the fraction of energy
available to increase the depth of the CBL. The effect
on inversion evolution of warm air advection above
the inversion layer can be investigated by using a
modified version of the model in which the wanning
rate is input.
The model simulates the changes with time of the
height of the inversion top and the depth of the CBL
during the inversion breakup period. From these sim-
ulations potential temperature profiles of the valley
atmosphere can be constructed for any time during
the period.
Sensitivity analyses were conducted for the lim-
iting cases of the model. The results indicate that the
time required to destroy an inversion depends pri-
marily on the initial height of the inversion, on its
potential temperature gradient, and on the amount
of energy available to destroy it. Using a reference
simulation in which model parameters were given
values typical of valley inversions, inversion destruc-
tion took approximately 3'/i to 4Vi h after sunrise.
These times correspond well with actual observa-
tions. Less time is required to destroy a valley in-
version than an inversion of like dimensions over the
plains, because the available energy is used to warm
a smaller volume of air. For the dry valleys of western
Colorado, the amount of energy available depends
to a large extent on the presence or absence of snow
cover or surface moisture in the valley. Valley in-
versions were destroyed sooner by the growth of a
CBL than by the descent of the inversion top. Valley
width and inclination angles of the sidewalls had only
a limited effect on the time required to destroy an
inversion. Increased valley width and steeper side-
walls both increased slightly the time required.
The thermodynamic model was used to simulate
two specific sets of inversion breakup data for Pattern
2 and 3 temperature structure evolution in the to-
pographically diverse Eagle and Yampa Valleys.
Simulations were obtained by fitting two constants
in the model (relating to the surface energy budget
and energy partitioning) to the data. The model out-
put fit the Pattern 2 inversion breakup in the snow-
covered Yampa Valley very well using an energy
input equal to 19% of the extraterrestrial solar flux
on a horizontal surface, and assuming that all of this
energy was used to drive the slope flows. A good fit
to the Eagle Valley data was obtained using an en-
ergy input equal to 45% of the extraterrestrial solar
170
-------
302
JOURNAL OF APPLIED METEOROLOGY
VOLUME 21
flux and assuming that 14% of this energy was used
to cause the valley floor CBL to grow. The remaining
energy was used to remove mass from the valley in
the slope flows.
Results of the present study provide new insights
into the evolution of valley temperature structure and
quantify the influence of the various parameters af-
fecting temperature inversion breakup. The model
explains the importance of the initial sunrise inver-
sion characteristics; the observed timing of the be-
ginning of inversion destruction; the mean time re-
quired to destroy typical inversions in the deep
valleys of western Colorado; the weak seasonal de-
pendence of the time period required to destroy the
inversions; the effects of snow cover and ground mois-
ture and of valley topography; the patterns of warm-
ing observed in the various layers of the temperature
structure; the typical observed inversion top descent
rates of 40-150 m h~'; and the retarded growth of
the valley CBL's relative to the flat plains case. The
thermodynamic model, while implicitly incorporat-
ing up-slope mass transport, is able to simulate tem-
perature structure evolution in a wide range of valley
topography without taking account of along-valley
wind systems.
Acknowledgments. The research was accom-
plished at the Department of Atmospheric Science,
Colorado State University, under funding from
Grant ATM76-84405, Atmospheric Sciences Sec-
tion, National Science Foundation. The manuscript
was prepared under funding from the Environmental
Protection Agency through Interagency Agreement
AD-89-F-0-097-0, with the U.S. Department of En-
ergy.
REFERENCES
Leahey, D. M., and J. P Fnend, 1971: A model for predicting
the depth of the mixing layer over an urban heat island with
applications to New York City. / Appl. Meteor., 10, 1162-
1173.
Sellers, W. D., 1965: Physical Climatology. University of Chicago
Press, 272 pp.
Whiteman, C. D., 1982: Breakup of temperature inversions in deep
mountain valleys: Part I. Observations. J. Appl. Meteor., 21,
270-289.
171
-------
APPENDIX D
SUMMARY OF MODIFICATIONS TO VALMET
173
-------
APPENDIX D
SUMMARY OF MODIFICATIONS TO VALMET
VERSION 1.0, DATED APRIL 30, 1983
Original version of VALMET.
VERSION 1.1, DATED DECEMBER 1. 1984
Major code changes:
1. Main Program
- The call to subroutine PRISE was fully incorporated .into
the code (it had been "commented out" in the original
version).
- The code which corrects concentrations to standard
conditions was moved from subroutine GAUSS to the main
program.
- The call to subroutine GAUSS for the calculation of the
plume center!ine concentration was corrected to allow the
plume to be offset from the valley centerline. Correction
to standard conditions is now accomplished following the
call to GAUSS.
2. Subroutine INPUT
New code was developed to check the model parameter values
that are input by the user. An error message is written
when the user inputs values outside normal atmospheric
ranges.
- Calculations of NBOX, NTOTI and AC were moved from the main
program into the INPUT module to facilitate the new error
checking code.
- A new error message scheme was developed for the INPUT
module.
175
-------
Subroutine PRISE
- A new plume rise module was developed from the MPTER code
to replace the old module derived from the CRSTER code.
- A new code was developed in order to bypass the plume rise
calculations when SRAD and/or SVEL are set to zero.
Subroutine GAUSS
- The section of code that corrects concentrations to
standard atmospheric conditions was removed and placed in
the main program.
Subroutine INGRAT
- A statement was added to this subroutine.
Subroutine BRKUP
- A new exponential decay algorithm was added to handle
variable time step lengths.
176
-------
TECHNICAL REPORT DATA
ffletue read Imtructions OH tfie rtvtrte be fort comptellngj
REPORT NO.
2.
3. RECIPIENTS ACCESSION NO.
4. TITLE AND SUBTITLE
GREEN RIVER AIR QUALITY MODEL DEVELOPMENT
VALMET - A Valley Air Pollution Model
6. REPORT DATE
6. PERFORMING ORGANIZATION CODE
7. AUTMOR(S)
. PERFORMING ORGANIZATION REPORT NO.
C. D. Whiteman and K. J. All wine
PNL-4728
. PERFORMING ORGANIZATION NAME AND ADDRESS
Battelle, Pacific Northwest Laboratory
Richland, Washington 99532
10. PROGRAM ELEMENT NO.
CDTA1D/09-0726 (FY-85)
11. CONTRACT/GRANT NO.
IAG AD89F20970
12. SPONSORING AGENCY NAME AND ADDRESS
Atmospheric Sciences Research Laboratory, RTF, NC
Office of Research and Development
Environmental Protection Agency
Research Triangle Park, North Carolina 27711
13. TYPE OF REPORT AND PERIOD COVERED
Final FY-85
14. SPONSORING AGENCY CODE
EPA/600/09
18. SUPPLEMENTARY NOTES
16. ABSTRACT
Following a thorough analysis of meteorological data obtained from deep valleys
of western Colorado, a modular air pollution model has been developed to simulate
the transport and diffusion of pollutants released from an elevated point source in
a well-defined mountain valley during the nighttime and morning transition periods.
This initial version of the model, named VALMET, operates on a valley cross section
at an arbitrary distance down-valley from a continuous point source. The model has
been constructed to include parameterizations of the major physical processes that
act to disperse pollution during these time periods. The model has not been fully
evaluated. Further testing, evaluations, and development of the model are needed.
Priorities for further development and testing are provided.
17.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.IDENTIFIERS/OPEN ENDED TERMS C. COSATI Field/Croup
18. DISTRIBUTION STATEMENT
19. SECURITY CLASS (Tilts Report)
UNCLASSIFIED
21. NO. OF PAGES
RELEASE TO PUBLIC
20. SECURITY CLASS (Tills page I
22. PRICE
UNCLASSIFIED
EPA F*rm 2220.1
. 4-77) ^««VIOU1 BDITION i« OBIOLCTC
------- |