May 1984
v-.'.--»- > A
H > C ^
-1 3 l« „! V
^i r* * /%, -
i^v;-.^'
,^
A REGIONAL-SCALE (1000 KM) MODEL OF PHOTOCHEMICAL AIR POLLUTION
Part 2. Input Processor Network Design
ENVIRONMENTAL SCIENCES RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
RESEARCH TRIANGLE PARK, NC 27711
-------
A REGIONAL-SCALE (1000 KM) MODEL OF PHOTOCHEMICAL AIR POLLUTION
Part 2. Input Processor Network Design
Robert G. Lamb
Meteorology and Assessment Division
Environmental Sciences Research Laboratory
Research Triangle Park, North Carolina 27711
ENVIRONMENTAL SCIENCES RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
RESEARCH TRIANGLE PARK, NC 27711
-------
NOTICE
This document 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
After the model described in Part I of this report was formulated, a
draft of an instruction manual was rather hastily prepared to guide computer
programmers in the task of transforming the theory into an operational model.
The present report, Part 2, evolved from that manual. The purpose of this
document goes beyond that of an instruction manual, however. The broader
objective is to provide in conjunction with Part 1 of detailed description of
what we regard as EPA's first-generation regional oxidant model.
In attempting to use science as a tool for treating the types of applied
problems that are of concern to the EPA, one is not allowed the luxury of
simplifying assumptions that reduce problems to forms that possess concise
elegant solutions. Instead, one must face the harsh realities of the physical
world and search for approximate descriptions of phenomena that strike an
acceptable compromise between scientific rigor and practicability. Just what
constitutes an "acceptable" compromise in this case is a subjective judgement
that each person must make for himself. In my view, several of the techniques
presented in this report represent compromises that are not wholly acceptable,
but they must suffice for now because time constraints dictate that we move on
to the task of model testing. Hopefully, the flexibility that we have built
into the basic framework of the model will foster efforts by others to expand
and improve upon the work we have done;" and a second generation model will
emerge significantly better than the model presented here.
R.G. Lamb
November 1983
m
-------
ABSTRACT
Detailed specifications are given for a network of data
processors and submodels that can generate the parameter fields
required by the regional oxidant model formulated in Part 1 of
this report. Operations performed by the processor network in-
clude simulation of the motion and depth of the nighttime rad-
iation inversion layer; simulation of the depth of the convec-
tive mixed and cloud layers; estimation of the synoptic-scale
vertical motion fields; generation of ensembles of layer-aver--
aged horizontal winds; calculation of vertical turbulence fluxes,
pollutant deposition velocities, parameters for a subgrid-scale
concentration fluctuation parameterization scheme; and many
other functions. This network of processors and submodels, in
combination with the core model developed in Part 1, represent
the EPA's first-generation regional oxidant model.
IV
-------
CONTENTS
Preface i i i
Abstract i v
Figures viii
Tab 1 es xi
Acknowledgements xi i i
1. Introduction 1
General Discussion 1
Summary of Processor Functions 5
Processor PI 5
Processor P2 5
Processor P3 7
Processor P4 7
Processor P7 7
Processor P8 7
Processor P9 7
Processor P10 8
Processor Pll 8
Processor P12 8
Processor P15 8
BMC 8
Summary of Model Equations 9
Layer 1 9
Layer 2 12
Layer 3 14
Layer 0 15
2. Processor PI 18
General Discussion 18
Step 1 18
Step 2 19
Step 3 20
Step 4 20
Step 5 : 21
Step 6 21
Step 7 22
Step 8 23
Step 9 23
Step 10 24
Stage INT 24
3. Processor P2 30
General Discussion 30
Stage LBC: Estimating lateral boundary conditions
from monitori ng data 34
-------
CONTENTS (continued)
Stage 1C: Estimating initial conditions from
monitoring data 34
Stage UBC: Upper boundary conditions 36
4. Processor P3 39
General Discussion 39
Step 1 39
Step 2 40
Step 3 41
Step 4 42
Stage INT 42
5. Processor P4 46
General Discussion 46
Step 1 47
Step 2 48
Step 3 50
Step 4 51
6. Processor P7 56
Motion of a viscous, hydrostatic fluid of constant
density over irregular terrain 56
The pressure force 58
The fri cti on force 61
The Coriolis force 65
The momentum equations 65
The fluid depth equation 67
Simplified model equations 76
Solution of the u, v, and z equation 79
Stage ZT Y? 82
Stage DELRO 83
Stage ETA 86
Stage PCD 89
Stage IBC 91
Stage H1HO 94
Stage SIG 95
Stage FLOMOD 96
Appendix A to Section 7 101
Appendi x B to Secti on 7 110
Calculation of the BAR variables 120
Calculation of the PRIME variables 122
Calculation of the dependent variables uX and h1 124
7. Processor 8 131
Introduction 131
Derivation of basic equations 131
Stage ZQ 147
Stage PATH 161
Stage WEWC 165
Stage W2 172
VI
-------
CONTENTS (continued)
Stage WV 182
Mode 0: Apt=0 182
Mode 1: Api=l 185
8. Processor P9 197
Devel opment 197
Stage DEN 197
Stage ZEN 203
Appendix to Section 8 195
9. Processor P10 207
Stage DELH 208
Stage S 211
Stage ZTA 214
10. Processor Pll 219
Summary 219
Introduction 220
Average horizontal winds in a layer bounded by two
arbitrary pressure surfaces 222
Stage UV11 247
11. Processor P12 254
Devel opment 254
Stage K 254
Stage WWO 261
Stage WW1 264
12. Processor P15 271
Devel opment 271
Step 1 271
Step 2 275
13. The B-Matrix Compiler 279
Introduction 279
The B-Matrix elements 280
Preparation of terms in the F-equation 285
Step 1 287
Step 2 288
Step 3 289
Step 4 290
Step 5 290
References 296
VI1
-------
FIGURES
Number Page
1-1 Schematic illustration of the regional model and the
network of processors that supply it information 2
1-2 Modeling region in the NEROS study. Each dot represents
a grid cell 6
1-3 Illustration of Processor PI and its input and output
interfaces with the processor network 29
2-1 Illustration of Processor P2 and its input and output
interfaces with the processor network 38
3-1 Schematic illustration of Processor P3 and its input and
output interfaces with the processor network 45
4-1 Schematic illustration of Processor P4 and its input and
output interfaces with the processor network 55
7-1 Illustration of the variables used in the flow model 57
7-2 Air parcel considered in the force balance analysis 58
7-3 Friction forces on a fluid parcel of horizontal dimensions
(Ax,Ay) bounded by z and z. 62
7-4 Projection of the horizontal rectangle (Ax,Ay) onto the
surface H centered at point Q 69
7-2 Illustration of the 54 5 min x 5 min cells that are used in
the calculation of z.(I,J) 83
7A-1 Illustration of the points used in the numerical solution
of Eq. 7A-1. Circled grid points are those from which
values f and S are taken to derive a biquintic expansion
of F(x',tJ about the point (IST,JST) (see Eq. 7-49) 104
*^* n
7B-1 Flow chart of FIOMOD operations 119
7B-2 Grid network on which p , i , $.> S , S , and S, are computed.
Different spatial derivative operators Ax and Ay are required
in these calculations as indicated 122
7B-3 Illustration of the southwest corner zone and the west and
south boundary zones of the mode! domain 129
-------
FIGURES (Continued)
Number Page
7-6 Schematic illustration of Processor P7 and its input and
output interfaces with the processor network 130
8-1 Illustration of surfaces H2 and H3 during situations in
which convective clouds cover only a portion of the
modeli ng regi on 132
8-2 Illustration of the air parcel trajectory that arrives at
point x at hour ti. Point x' denotes the parcel's
location at time t0 < tx .. .7. 141
8-3 (a) Idealized profiles of mixing ratio q and potential
temperature 0 in dry, convective conditions.
(b) Second derivatives of the profiles illustrated in
panel a 150
8-4 Illustration of th points (dots) at which measurements of £
are available; and the point z" at which a measure of d2£/dz2
is desired. Values of £ measured at the points z.,i=l,...4
centered at z" are used to approximate £(z) in a polynomial
about z" 151
8-6 Schematic illustration of Processor P8 and its input and
output interfaces with the processor network 196
9-1 Schematic illustration of Processor P9 and its interfaces
with the processor network 206
10-1 Illustration of the influence of terrain on model layers 0,
1 and 2 for given values of the penetration fractions crT1
and aTQ .. 211
10-2 Schematic illustration of Processor P10 and its interfaces
with the processor network 218
11-1 Surfaces bounding layer p in which M soundings of the
horizontal wind are available on the bottom surface p
of Layer p 224
11-2 Illustration of Processor Pll and its input and output
interfaces with the processor network 253
12-1 Illustration of the superposition of five hypothetical
realization of an ensemble of point source plumes. The
width I of the ensemble of plumes is controlled by the
character of the flow field ensemble. The width a of
the plumes in the ensemble is controlled by the turbulent,
eddy diffusivity K 256
-------
FIGURES (continued)
Number Page
12-2 Illustration of Processor P12 and its input and output
interfaces with the processor network 270
15-1 Illustration of Processor P15 and its input and output
interfaces with the processor network 278
BMC-1 Schematic illustration of the B-Matrix Compiler (BMC).
The input interface is the Model Input File (MIF). The
output of the BMC is the "b-matrix tape" which is read
by the model code CORE (see Figure 1-1) 293
-------
TABLES
Number Page
1-1 Summary of input and output variables of each step of
Processor PI 26
2-1 Pollutant species concentration (ppm) thaken to be
representative of "clean" atmospheric conditions 31
2-2 Summary of input and output variables of each stage of
Processor P2 37
3-1 Summary of input and output parameters of each step of
Processor P3 44
4-1 Summary of input and output parameters in each step of
Processor P4 53
7-1 Summary of the input and output requirements of each stage
of Processor P7 97
8-1 Summary of the input and output requirements of each stage
of Processor P8 188
9-1 Summary of the input and output variables of each stage
of Processor P9 and thei r sources 205
10-1 Summary of input and output parameters for Processor P10
and thei r sources 215
11-1 Input requirements of stage UV11 and their sources 251
12-1 Summary of the expressions used to estimate the horizontal
eddy diffusivity Kn in each of the model's three layers 259
12-2 Input and output variables of each stage of Processor P12 267
15-1 Deposition resistances (sec/m) for S02 as a function of
land use type n and stability L 272
15-2 Deposition resistances (sec/m) for ozone as a function of
1 and use type n and stabi 1 i ty L 273
15-3 Deposition resistances of several pollutants relative to that
of ozone over agricultural land. (Deduced from data of Hill
and Chamberlain, 1971) 274
XI
-------
TABLES (Continued)
Number Page
15-1 Summary of the input and output parameters of each step
of Processor P15 277
BMC-1 Definitions of parameters in the Model Input File (MIF) 294
-------
ACKNOWLEDGEMENTS
I am indebted to Dr. Don Relies of Computer Data Systems, Inc. (CDSI) for
his invaluable assistance in developing the ideas on the flow field ensemble
presented in Section 10, and to Ms. Achamma Philip, also of CDSI, for her
steadfast assistance in developing the numerical boundary condition scheme
presented in Section 6. I also want to acknowledge the excellent work of
Gayle Webster and Sherry McCoy of Systems Research and Development Corporation
(SRD), in preparing the manuscript.
XI11
-------
SECTION 1
INTRODUCTION
GENERAL DISCUSSION
In Part 1 of this report, (Lamb 1983d) we developed a theoretical basis
for a regional-scale model of photochemical air pollutants. Realizing that
the operational model would be very complicated and would require considerable
time and effort to develop, we proposed that, rather than integrate all of the
various components of the model into a single unit, we construct it by
partitioning the mathematical descriptions of small groups of individual
physical and chemical processes into discrete modules interconnected by fixed
communication channels. Such modular design would facilitate troubleshooting
operations, provide a natural division of labor for the tasks required to
implement the model, and permit continual incorporation of state-of-the-art
techniques without the need to overhaul the model code each time a new
technique was introduced.
The overall structure of the proposed model system is illustrated in
Figure 1-1. The box labeled CORE represents the computer language analogue of
the differential equations that describe all the governing processes
considered in the model development in Part 1. (These equations are listed at
the end of this section.) The CORE module is expressed in a very primitive
mathematical form in the sense that its inputs are matrices and vectors whose
------- An error occurred while trying to OCR this image.
-------
elements are composites of meteorological parameters, chemical rate constants,
etc. For example, the link between CORE and the output of the module labeled
CHEM, which contains the analogue of the chemical kinetics scheme, consists of
two vectors P and Cj, each of length N, where N is the total number of chemical
species simulated. The n-th element of P is the net rate of production of
species n due to source emissions and chemical reactions among all other
species, and the n-th element of Q is the net rate of destruction of species n
due to its chemical interaction with all other species. Thus, any chemical
kinetics mechanism can be incorporated into the model as long as it is
expressed in a form that is compatible with the vector interfaces that link
CORE with the chemistry module CHEM.
The remainder of the inputs required by CORE are prepared by the module
designated BMC (b-matrix compiler) in Figure 1-1.The BMC performs essentially
the same task that language compilers perform in computers. It translates the
parameter fields in the model input file (MIF) into the matrix and vector
elements that are required to operate the algorithms in CORE. These
parameters include layer thicknesses, horizontal winds in each layer,
interfacial volume fluxes, and deposition velocities.
The variables in the model input file (MIF) are supplied in turn by a
series of interconnected processors, labeled P7, P8, etc. in Figure 1-1,
several of which are rather complex models in themselves. These processors
generate the wind fields, the interfacial surfaces that separate the layers,
turbulence parameters, source emissions, and many other variables. Their
inputs consist of information generated by other processors in the network
-------
and partially processed raw data that are transferred through the processor
input file (PIF).
The purpose of this report is to provide detailed specifications of the
processor network illustrated in Figure 1-1. This network will consist of
both permanent and interchangeable components. The permanent components are
the CORE, the BMC, the MIF and PIF, the communication channels that link the
processors and the PIF and MIF, and the interfaces between the processors and
the communication channels. All processors, i.e., PI, P2, etc. and CHEM are
interchangeable components of the network. Any or all of the interchangeable
elements can be replaced by other modules as long as the replacements are
compatible with the communication channel interfaces. By "interface" we mean
the set of input and output variables assigned to each processor. We can
think of each processor as being analogous to an electronic device that plugs
into a multireceptacle socket (the interface). Each receptacle that provides
an input to the processor is connected to a fixed signal source, and each
receptacle that receives an output from the processor acts as a signal source
on the network of communication lines. Considering the interchangeability of
the processors and chemistry module CHEM, one should view the processor
designs that we develop in this report as "first generation" versions that may
be replaced in the course of future tests and refinements of the model.
Neither the chemistry module nor the results of any test simulations are
discussed in this report. These topics will be discussed in later parts of
this series of reports. In the following sections we present designs of each
processor in the network including the BMC. Theoretical descriptions of the
-------
mathematics contained within the module 'CORE were given in Section 9 of Part
1. We will not elaborate further on this part of the model system in this
report other than to describe in more detail the numerical scheme that we
apply to the transport and diffusion portion of the model equations. These
details are given in Appendix A of Section 6 which describes Processor P7.
Figure 1-2 shows the region in the Northeastern United States to which
the regional model will be applied. Each point in the figure represents the
center of a grid cell in which values of all pollutant species concentrations
are computed by the model.
SUMMARY OF PROCESSOR FUNCTIONS
The functions performed by each processor in the model network,
illustrated earlier in Figure 1-1 are summarized below.
Processor PI
Prepares upper air data for use by other processors in the network.
Processor P2
Uses surface air monitoring data to estimate initial, lateral and upper
boundary conditions on pollutant species concentrations.
-------
000'69 A
01
CO
-------
Processor P3
Prepares surface meteorological data for use in other processors.
Processor P4
Estimates surface roughness, Obukhov length, surface heat flux, and
friction velocity.
Processor P7
Simulates the depth and motion of the nighttime surface inversion layer
and the depth of the daytime shear layer and provides smoothed terrain
elevations.
Processor P8
Computes depths of the convective mixed layer and cloud layer, cloud and
turbulent entrainment velocities at the mixed layer top, synoptic-scale mean
vertical motion on the top surfaces of model Layers 2 and 3, and layer
averaged horizontal wind divergences.
Processor P9
Calculates factors for correcting the chemical rate constants for
temperature, density, and sunlight variations.
-------
Processor P10
Transforms the source emissions inventory into source strength functions,
and estimates the plume volume fraction parameter required in the subgrid
scale chemistry parameterization.
Processor Pll
Computes ensembles of volume averaged horizontal winds for each model
layer (except the nighttime surface inversion layer which is treated in P7).
Processor P12
Calculates layer interface turbulent volume fluxes; horizontal eddy
diffusivities; and cumulus cloud flux partition parameter.
Processor 15
Computes pollutant species deposition velocities.
BMC
Compiles processor network outputs for input into the model core
algorithms.
-------
SUMMARY OF MODEL. EQUATIONS
The equations upon which the regional-scale mode is based were derived in
Part 1 of this study. The basic forms of the governing equations in each of
the model's four layers are summarized below.
Layer 1
where
\ a\
+ u ai
Fo,i
a cos(|)
a = earth radius at MSL (in meters)
\ = longitude
0 = latitude
A = a2 (A(j)AX) cost))
A({),AA. = latitude, longitude grid cell dimensions
= constants (A<|) = AA = ^°)
-------
A+AA/2 <)>+A<|>/2
2
acos<|) MAX {0,
A.-AA/2 <|>-A/2
- MAX [zT(A',f ), zo
VU L -, /
> subgrid scale fluxes of c (see Part 1, Section 7)
1J
-, = all chemical, rainout and washout processes.
<$>-, = all emissions of c in Layer 1 (includes stacks and surface
sources above nighttime radiation inversion.
-,, , = layer averaged horizontal wind components
(u = east-west component, positive eastward),
(v = north-south component, positive northward).
Fl,l = (1-aTl> K1 - 2)Wlm + 1 V
P = deposition velocity of species c
Oy-,(\,<|>,t) = fraction of surface H^ penetrated by terrain
n w - w
T
C -
- ^ wm " ZT, neutral and unstable conditions;
wR1 - | ui l
Hi> ( = given inversion layer depth growth rate), stable
10
-------
wn, = mean vertical velocity on H-, (terrain induced component
Ui excluded) x
a = rms vertical turbulence on H, ( = 0 in stable cases)
"1 •*•
w = threshold of cumulus "root" updraft velocity on H,
»l>, o~ , w (see Layer 2 equations)
c c
9 I -f-2
erf(x) = — e L dt
Fo,l = (orTo
T = fraction of surface H in given grid cell penetrated by
terrain
Ozone:
F0(03) = <03>1 w.A.
- (1 - a) •
Nitric Oxide:
FQ(NO) = 1 w_X_ - NO w
- (1 - a)
, otherwise
- a)
, otherwise
-------
Nitrogen Dioxide:
FQ(N02) = lW_\_ - N02
- a)
- (1 - a)
2 , otherwise
all other species, x-
F0(X) = i*_\_ - XW+A+(1 - |)(1 - a) - (1 - a)(v£x + S )(1 + p
A. A
a =
1 - a
- - fe/A6
See Layer 0 equations for specification of 03, NO, N02> x, w+, A+, etc.
Layer 2
af
= +
V2(\,<)),t) = a2cos<(>
X-AA/2
{z2(X',<))l,t) - MAX [zT(Xl,
-------
2
subgrid scale fluxes of c (see Part 1, Section 7)
2
2 = all chemical, rainout and washout processes in Layer 2
2 = source emissions in Layer 2 (if any)
2, 2 = Layer 2 averaged winds.
F2 2 = c, c (c - ) + - 9 [.
t. 31. ACT C (L c. O\» tit/
l,2 = aTlP2 + (1 " CTT1)C1W1 + (1 " 2)Wlm + (2 "
W-, ,w, , see Layer 1; WQ, , see Part 1, Eq. 4-44c)
a = fractional area coverage of cumulus clouds
w = mean upward velocity in cumulus clouds
fe
7^ = turbulent entrainment rate of inversion air into mixed layer
= mean vertical velocity (mean of both terrain induced component w-™
and divergent component w™)
-2
= fraction of cloud air from surface layer
0 < 1)1 < 1 and
"""~~
.. Tn
— T T | Q C
13
-------
- w
Vc = -
c, c', I, w+, \+ = See Layer 0 equations
Layer 3
a3 .A rr: r i _
5*\7T Lro o ("o oJ - v-^3
X+AX/2 ()>+A3 = all chemical (wet and dry) rainout and washout processes
3, 3 = layer averaged winds
W ~ W
F _ _ r !2 *£ 4
« - - o L ^ _ _.. - iu/tiujv-_ ^'—'s^ ' ^'3 aF"
-------
3 ^ , if dH3/dt < 0 ;
I
3,3
- 3) dHT/dt + , ^ , otherwise
«3 O O t
c» = concentration of species c above z3
dH3/dt = given volume flux through model top surface
c = (1 - <|02 + i|/[|c' + (1 - 4)c]
where c1 and c are defined with the Layer 0 variables.
Layer 0
<03>0 = (1 - C)
0 = w_\_<03>1
3 "
» otherwise
0 = (1 - O NO +
N° =
w.\.1
"
w+\+ +(1 -
15
-------
NO1
I , ,Mn
, otherwise
0 =
N02
•N0£v -H 03v * 3
N02v + (NQ + NQ
z NU NU2 , otherwise
Species x other than 03, NO, and N02:
o =
X =
V^Ti^T^
vv + ?
X """ ft 4- \)
Parameters:
X. = %[1 - erf
16
-------
\ =
(j •«
W Z^
- — exp ( -
w
+ [l - erf (
v = u* (plume entrapment velocity)
= plume volume fraction
a = rms vertical turbulent velocity on H
17
-------
SECTION 2
PROCESSOR PI
GENERAL DISCUSSION
This Processor performs a number of standard operations on the raw rawin
data to put them into the forms required by the higher level processors in the
network. It operates on the rawin data only. We assume that the "raw"
rawin record at the m-th station consists of a sequence of I "vectors"
[SPD.DIR.T.TD.p]^, i=lf...IB (1-1)
where the i-th vector represents the wind speed (m/sec), direction (degrees),
temperature (°C), and dew point depression (°C) at pressure level p
(millibars) at a given station m. The number of observation levels Im in each
record can vary from station to station and from hour to hour. The steps
necessary to convert the raw data into the desired forms follow.
Step I
Convert the speed and direction (SPO.DIR)^ into cartesian components
(u,v) at each level i and at each station m by
18
-------
where
eim = (DIRim)'(2*/360). (1-3)
These values must be processed further (in step 10) before recording in the
PIF.
Step 2
Convert the temperature and dew point depression TD into the water vapor
mixing ratio q (mass of water per total mass of air) at each level i and
station m by
0.622e.
*• = TV£
where
R = 0.461 joule g" °K~
L = 2500 joule g (latent heat of vaporization)
e = 40 mb (saturation vapor pressure at T )
TQ = 29 + 273 = 302°K (reference temperature)
T... = T. - TD. + 273 (dew point temperature at level i, station m)
uim im im
p. = pressure (mb) at level i, station m.
The mixing ratio q. and dew point T,,. are outputs of this stage and
require further processing in Step 6.
19
-------
Step 3
Compute the relative humidity RH. for each level 1 and station m by
where e is found from (1-5) by replacing the dew point Tn. with the
n.
U I Ul
temperature T. .
Step 4
Convert the pressure p. at the levels 1=1,2,...I of the upper air
observation at station m to altitudes z. (AGL) by first converting the
temperature measurements T. to virtual temperatures
1 HI
Tvim •
-------
Step 5
Let z be the known elevation (MSL in meters) of the m-th rawin station,
and let pQm and TVQ(n be the station pressure and virtual temperature. Using
(1-7) we now calculate the sea-level pressure at station m:
p,-^ = sea-level pressure at station m (in millibars)
zm 9
= Pomexp[JU - ]; (1-8)
om Kd 'vom
and we calculate the geopotential height Qm (m2/sec2) of the 1000 mb surface:
*om '
+ Vvo.1" ITOF •
Using the value from the previous observation time, say, t-At, compute
Record om and 4>om in the PIF.
Step 6
At this stage the virtual temperature, relative humidity, dew point and
mixing ratio are known at each elevation z. in the sounding (from step 3).
These data should now be interpolated to give semi-continuous soundings, e.g.,
T (z) and q (z) of virtual temperature and mixing ratio at each rawin station
m. That is, we convert
21
-------
where
z = z + nAz , n=0,...
Az=50 meters, and
5000-z.
m bO VJ
The third-order interpolation scheme described in Processor P8 (see Eqs. 8-56,
8-60, ...69, and 8-70) should be used in performing the transformations
indicated by Eqs. (1-10). Record the Tfflf qm> TDm, and RHm values for each of
the z given by (l-10d) in the PIF.
Step 7
The pressure-height pairs (p^, z,-m) should be transformed into a semi-
continuous sounding as in Step 6 above. In this task the sequence begins with
the sea-level pressure, followed by the station pressure, and the observation
levels aloft, i.e.,
(P-lm'0)'(Pon,'y'(Plm'zlm)---(Plmm'zImm) => Pm(z)> (1'lla)
where
z=nAz, n=0,l,...100; (1-llb)
and Az=50 meters. In performing the interpolation (1-lla), an exponential
formula should be employed based on the hydrostatic condition
(1-12)
22
-------
Record the p values at each z given by (1-llb) in the PIF.
Step 8
Using the Pm(z) profile obtained in Step 7, determine the geopotential
_2
(m2sec ) of the 850, 700, and 500 millibar surfaces at each station m. For
example,
*(850)m = 9Z850
where g is gravity and zQ50 is the elevation (MSL) that satisfies
= 85° mb"
By the same process compute
-------
-2 _1 _3
where pm is in millibars and [^=287 m2sec °K~ , Tvm is in °K, and p=kg m~ .
From 9 and p compute the profile a (z) of static stability from
din
_ -[em(z+Az)-em(z-Az)]-io~2
CTsm(z) ~
-3 ..2
With p in millibars and p in kg m , a has units of m4 sec2 kg . Record the
profiles 0m(z),pm(z) and (Jsm(z) in the PIF for each rawin station m.
Step 10
Interpolate the wind components u- and v. of Eq. (1-2) above to obtain
the profile
where z takes the values given by (lOd).
Stage INT
The raw rawin data are generally available only at 12-hour intervals,
but the output variables produced by this processor, PI, are required each
hour by processors further along in the network. Therefore, an interpolation
of all output variables must be performed to provide values at hourly intervals.
The specific interpolation formula that is used for this purpose is left to the
discretion of the user.
24
-------
Table 1-1 summarizes the inputs and outputs of each step of Processor PI
and Figure 1-3 illustrates the processor and its data interfaces.
25
-------
Table 1-1 Summary of input and output variables
of each step of Processor PI.
Input
Variable
SP01m
DIRfm
Description Source Step
wind speed (ms~ ) at RAW 1
level i at rawin
station m.
Wind direction RAW
(compass degrees) at
level i at rawin
station m
Output
Variable
uim
-------
Table 1-1 Summary of input and output variable
of each step of Processor PI. (Continued)
Input
Variable
zm
in
p
om
vnm
V Ulll
vim
V 1 III
q.
RHim
i in
zim
TDim
*"i m
1 ill
zim
1 III
pm«
Pm(z)
in
T (z)
vm
Description
elevation (meters,
msl) of rawin
station m
station pressure
(mb) at station m.
virtual temperature
(°C) at station m.
see above
see above
see above
see above
see above
see above
see above
pressure at elevation
z(MSL) over station m.
see above
see above
Source Step
RAW 5
RAW
Step 4
Step 4 6
Step 2
Step 3
Step 4
Step 2
RAW 7
Step 4
Step 7 8
Step 7 9
Step 6
Output
Variable Description
'OV geopotential (m2sr2)
om K of 1000 mb surface
at station m, hour t^.
*««,(*!,) time rate °f change
om k (m2s-3) of geo-
potential of 1000 mb
surface at station m.
T (z,t. ) Temperature, mixing
ratio, relative humidity,
q (z,tk) and dew point profiles
at station m resolved
RH_(z,t. ) to Az=50m as prescribed
by Eqs.
TDm(z,tk) (l-10d,e) at hour tk.
pm(z,t.) pressure (mb) at
elevation z(msl) above
station m, resolved to
Az=50m as prescribed
by Eqs. (l-10d,e) at
hour t^.
•fr/ocmm geopotential (m2sec~2 )
*/7ftnX of the 850, 700, and
*r(nitt 500mb surfaces at
^ouujm _4.,4.i~_ _.
station m.
9 (z) potential temperature
m (°K) at elevation z
over station m, resolved
to Az=50m as prescribed
by Eqs. (l-10d,e).
p (z,t. ) air density (kgm~ ) at
elevation z (resolved
as above) over station m,
hour t, .
acm(z,t,,) static stability
Sm K /n,4c2Un-2>\ a-f olowatirtn
27
z at station m,
hour t. .
-------
Table 1-1. Summary of input and output variables
of each step of Processor PI. (Concluded).
Input Output
Variable Description Source Step Variable Description
u. see above Step 1 10 "(x^.z.t. )• east-west and north-
lra ~" k south wind
v. see above ' Step 1 ^Xm'2'^ at elevat"'on z (as
[also given by Eqs. l-10d,e
"m>vm] at location x of
m m station m, atmhour
28
-------
T
YYYYYYYYYYYYYYY
-L
29
-------
SECTION 3
PROCESSOR P2
GENERAL DISCUSSION
This processor determines initial, upper and lateral boundary conditions
on all of the pollutant species simulated by the regional model. Presently,
it is impossible to estimate these conditions with an.accuracy anywhere
near that of available emissions estimates, because: (1) the pollutant
monitoring data from which the initial and boundary conditions must be derived
are few and nonuniformly distributed, (2) no measurements are routinely
available aloft, and (3) most of the intermediate pollutant species that must
be treated to simulate the chemistry properly are not measured at all. In
addition, the model requires conditions on the cell averaged concentrations
but only point measurements are made at the monitoring sites.
The problems caused by the paucity of data can be mitigated by
initializing the model on a "clean" day, such as a day immediately following
the passage of a cold front, and by choosing a model domain that is large
enough that the quantities of pollutants emitted by sources just outside the
model boundaries are small compared to those generated within the simulation
area itself. In such a case, we can use "clean atmosphere" conditions for the
initial, upper and lateral boundary values of each species. Table 2-1 lists
these values for each of the 23 pollutants modeled by the Demerjian-Schere
(1979) kinetics scheme that we have been using during the development phase of
the regional model.
30
-------
Table 2-1. Pollutant species concentrations (ppm) taken to be
representative of "clean" atmospheric conditions*.
(Notation: 1.000-01 = 10" .)
NO
HC2
HN02
H202
H02
R20
4.
3.
3.
8.
2.
2.
499-04
390-03
473-05
784-05
078-05
452-09
N02
HC3
HN03
0
H04N
R102
1.
7.
7.
1.
7.
1.
404-03
679-03
215-04
284-10
583-06
550-05
CO
HC4
PAN
N03
RO
R202
1.
6.
3.
5.
5.
5.
010-01
911-04
808-04
434-07
848-08
319-06
HC1
03
RN03
HO
R02
1.
3.
1.
3.
5.
084-03
522-02
004-06
734-07
578-05
*Values reported here were obtained by initializing the Demerjian-Schere
(1979) 23-species kinetics mechanism with the following species concentra-
tions and allowing reactions for 90 minutes in full sun conditions:
NO = .OOlppm, N02 = .002ppm, total non-methane hydrocarbon = .05ppm,
CO = .Ippm, 03 = .02ppm, all other species = Oppm.
Hydrocarbon classes are as follows:
HC1 = olefin, HC2 = paraffin, HC3 = aldehydes, HC4 = aromatics. Initial
values of each of these lumped species were obtained from the total
nonmethane hydrocarbon concentration using the speciation method
suggested by Demerjian (1983).
31
-------
The problems caused by the unavailability of volume-averaged
concentration data cannot be eliminated because there is no unique
relationship between the concentration values measured at one or several
discrete points within a given volume of space and the concentration averaged
over that volume. This is the so-called subgrid-scale closure problem
encountered in all modeling studies in which it is impossible to choose a grid
size small enough to resolve the smallest spatial variations in the simulated
variables. In Part 1, Section 5, we developed a crude subgrid-scale closure
scheme for use in treating the pollution chemistry. Proceeding along lines
similar to those described there we could formulate an approximate way of
extracting volume-averaged concentration estimates from point measurements.
However, we will not attempt this here because the improvement in accuracy
gained in the initial concentration fields would probably not be significant
enough to warrant the development and implementation efforts. Perhaps future
modeling studies can investigate this problem in detail if the need is great.
In the remainder of this section we outline a procedure for obtaining
rough estimates of initial and boundary conditions on pollutant concentrations
in situations where the "clean atmosphere" assumption does not apply. An
important point to note in the preparation of initial and boundary conditions
is that, due to inaccuracies and uncertainties in the methods used, the set of
concentrations obtained for any given grid cell or boundary point will
generally not be consistent with chemical equilibrium conditions. For
example, if one deduces concentration values for 03, NO, N02 and olefin from
the monitoring data and assigns "nominal" or zero values to all the other
species included in the chemical kinetics scheme, one would find upon
initializing the model with these values that a period of rapid chemical
32
-------
transformations immediately ensues. These rapid changes indicate that the
concentration conditions selected for the initial state do not represent a
state near equilibrium. These spurious reactions are an artifact of the
chosen set of concentrations and are not representative of the chemical
activity that occurs just after the initial moment. This is analogous to the
initialization problem in meteorological models where failure to prescribe an
initial state that is in geostrophic balance results in the generation of
spurious gravity waves.
The transient concentration variations excited by errors in the initial
state diminish the accuracy of the model's predictions within some finite
period following the initial instant tQ. They also exact a significant
penalty in computer time. Because when the chemical state is far from
equilibrium, the mathematical algorithm in the model that handles the chemical
rate equation must utilize many small time steps to track the approach of the
state of the system to equilibrium. Since this operation must be performed
initially at every grid point in the model and at all boundary points at all
times where the boundary condition specification is inexact, a considerable
portion of the computation time required by the model can be wasted on this
fictitious phenomenon. The remedy is to use the initial concentrations
deduced from the monitoring data as the initial state in a batch reactor
model; to run that model for a time long enough for the chemical state to
settle down to a point where changes are occurring relatively slowly (say,
time scales ~10 min); and then to use the concentrations of each species given
by the batch reactor model at that point as the initial conditions in the
regional model. The same procedure should be used to obtain the upper and
lateral boundary conditions.
33
-------
Stage LBC: Estimating lateral boundary conditions from monitoring data
(1) Collect hourly surface monitoring data for each species x-(ppro) i =
1....IMAX from all stations within 20 km either side of each of the
regional model's four lateral boundaries. Each station within this 40 km
wide boundary zone is treated as though it lies on the model boundary at
the point closest to the station location.
(2) Use a cubic spline or other acceptable interpolation formula to estimate
concentrations x^ U-g, <]>g, t)(ppm) at each grid point (Ag, <|>g) on the
boundary. Here (\,<|>) denotes the longitude and latitude of a cell center
on the boundary of the regional model domain.
(3) Using the functions x -j obtained in step 2, estimate layer averaged
concentrations along the boundaries as follows:
= B x,- (AR, V t), n=l,2,3 (2-1)
1 B b n n i B B i=l,...IMAX
where the B are empirical constants to be estimated from the NEROS field
n
experiment data.
(4) For each hour t , each boundary point (An, B), and each layer n perform
the batch reactor equilibration process, described in the introduction to
this section, to the set of concentrations g. t)>n, i=l,...IMAX
given by (2-1). Record the results in the ICBC portion of the model
input file MIF (units = ppm for each species).
STAGE 1C: Estimating initial conditions from monitor-ing data
(1) Collect surface monitoring data for all pollutant species at all stations
within the regional model domain at the initial hour t (=1200 EST) of
34
-------
the period to be simulated. Noon is chosen as the initialization hour
because, at this time, pollutants are usually distributed nearly
uniformly in the mixed layer and initial values for Layers I and 2 can be
equated with minimum error.
(2) In places where more than one monitoring station lies within a grid cell,
compute a weighted average of the reported values taking into account the
proximity of each station to sources and the distribution of land use
types within that cell. For example, if one monitoring site is in a
rural area and 70% of that cell is in rural land use, the rural station
would be given a weight of 0.7 in computing the cell average
concentration.
(3) Fit a two-dimensional function to the finite set of cell averaged
concentrations obtained in step 2 and from this function derive values of
the concentration of each measured pollutant species at all grid cells in
the model region.
(4) Subject each set of concentrations i=l,...IMAX obtained in step 3 to the
batch reactor equilibration process. Call the results of this operation
X,-(X> > t ) where (X, ) ranges over all grid points in the model
region. Now record the following initial layer averaged concentrations
(ppm) in the ICBC portion of MIF:
, t0)>! = X^X, , to) (2-2)
2 = Xjtt, , t0) . (2-3)
%(A, *, t0)>3 = ^Xjfr, 4>, t0) (2-4)
where |. is an empirical constant to be derived from the NEROS field
experiment data.
35
-------
Stage UBC: Upper boundary conditions
(1) In the first generation model we will use the "clean" atmosphere species
concentrations listed in Table 2-1 for the upper boundary conditions XQ,,
for all cells and all hours. Record these in the ICBC portion of
MIF(units = ppm for all species), i.e.
X* ,(X, «», t) = XCLEAN , 1XL....IMAX; (2-5)
jl LLtAN' 1 all (\, 4, t)
Table 2-2 summarizes the inputs and outputs of Processor P2, and Figure 2-1
illustrates the processor and its interfaces with the processor network.
36
-------
Table 2-2 Summary of input and output variables of each
stage of Processor P2.
Input
Variable
Description
Source Stage
Output
Variable
Description
X..|,(t )
1K m
concentration (ppm)
of i-th pollutant
at hour t measured
at surface monitoring
station k=l,...K.
RAW
LBC
n
average concen-
tration (ppm) of
i-th pollutant
in layer n=l,2,3
at boundary cell
(\D, d»D) at hour
tf
,-i,
IK
(t _)
Ml
see Stage LBC input. RAW
ic
average concentration
(ppm) of i-th pollut-
ant in layer n=l,2,3
in grid cell (A,$) at
the initial instant
t_.
X.j(tm)
see Stage LBC input. RAW
UBC
tn)
concentration (ppm)
of i-th pollutant
at top surface of
model over grid
cell (\,<)>) at hour
V
37
-------
H
D
CN a.
Q. h-
YYY
c
A
+*
CD
^
CO
C
A
o
+*
*
•e-
8*
3
a
*^
o
TJ
C
CO
**
3 .
a..*
IS
•— tt
s °
C/) b.
®
8^
a- j=
§ «
3 «-•
CN
£
3
a>
u.
38
-------
SECTION 4
PROCESSOR P3
GENERAL DISCUSSION
This processor performs standard operations on the surface meteorological
data to put them into forms required by the higher level processors in the
network. The surface and rawin data are treated in separate processors to
facilitate future alterations in the data analyses and to permit easier
incorporation of additional data.
We assume that at given time intervals (not necessarily hourly intervals)
the surface observations consist of the set
[SPD,DIR,TJD,p]n
where n denotes the surface station, whose location is x • and the other
*>Tl
variables represent wind speed (ms ), direction (compass degrees),
temperature (°C), dew point depression (°C) and station (not sea-level)
pressure (mb).
Step 1
Convert the wind speed and direction into north-south and east-west
components as follows:
Qn = Qfen> = " SPVsin 8n (
vn = v(xn) = - SPDn-cos 6n (
39
-------
where
= DIRn-(2n/360).
The components (u ,v ) are outputs of this processor for each time interval
Step 2
Convert the temperature and dew point depression to mixing ratio, dew
point temperature, relative humidity and virtual temperature as follows:
0.622e
qn = e (=mixing ratio) (3-2a)
n n
where
L 1 1
R = 0.461 joule g"1 °K~1
L = 2500 joule g"
e = 40mb (saturation vapor pressure at temperature T )
TQ = 302 °K
p = station pressure (mb)
and
TDn = Tn " TDn + 273 (=dew point temperature) (3-2b)
p ~e
RHn = qn On622e (-relative humidity) (3-2c)
where e is obtained from (3-2a') by replacing the dew point temperature Tg
in the formula with the dry bulb temperature T (expressed in °K);
Tvn " Tn ^ + 0l61
-------
where T (and hence T ) is in °K.
The mixing ratio q (dimensionless) and the relative humidity RH
(dimensionless) can be recorded in PIF (they are outputs of this processor).
Before recording the dew point and virtual temperature in PIF, they should be
converted to degrees Celsius as follows:
TDn
where TDn (°K) is from Eq. 3-2b and Tyn (°K) is from (3-2d).
Step 3
Let zn be the elevation (meters, MSL) of surface station n. Compute the
geopotential of the 1000 mb surface
*ons«VRdTvnlnl5ro (3'3)
where Tyn is from (3-2d) above (in °K) and Rrf = 287 m2s" °K~ . Using the
value of from the previous time interval, estimate the time rate of change
of * by
*on = t*on(t) " ^n^^0^ (3'4)
Compute the sea-level pressure p , at each station site x ,n=l,2,...N as
follows.
gz
Rsl n = Pnexp [^p-S- ] (3-5)
sl,n n RT
41
-------
where pn is the station pressure (mb) and all other variables have the same
definitions as above. If the surface temperature T undergoes large
variations between day and night, a 12-hour moving averaged temperature must
be used in (3-5) to avoid fictitious variations in the estimated sealevel
pressure.
Step 4
Compute the surface potential temperature and air density at each station
en = Tn (1M)°-286 (=potential temperature) (3-6)
"n
Pn = Tpp- -ICf2 (=density, kg m"3) (3"7)
d vn
where Typ is from Eq. 3-2d and is in degrees Kelvin (°K), pn is the station
pressure in millibars and
Rd = 287 m2sec" °K~1.
Record en(°K) and pn (kg m" ) in PIF.
Stage INT
Values of each of the parameters produced by this processor are required
at hourly intervals by higher level processors in the network. However, most
of the inputs to this processor, P3, may only be available at 3-hour
intervals. Thus, it is necessary to interpolate each of the parameter values
produced in this processor by some means adequate to produce reasonable
42
-------
estimates of parameter values each hour. Stage INT is intended to perform
this task. We leave the detailed specification of the interpolation algorithm
to the discretion of the user.
The inputs and outputs of each step of Processor P3 are summarized in
Table 3-1 and illustrated schematically in Figure 3-1.
43
-------
Table 3-1 Summary of input and output parameters
of each step of Processor P3.
Input
Variable
Description
Source
Step
Output
Variable
Description
SPD.
DIR.
surface wind speed RAW
(m/s) at station n
surface wind direction RAW
(compass degrees) at
station n.
east-west surface wind
component at station n,
hour t .
m
north-south surface wim
component at station n,
hour t .
TD
surface temperature RAW
(°C) at station n
surface dew point RAW
depression (°C)
at station n.
station pressure RAW
(mb) at station n.
surface mixing ratio
(dimensionless) at
station n, hour t .
surface relative humidf
(dimensionless) at
station n, hour t .
m
surface dew point (°C)
at station n, hour t .
surface virtual
temperature (°C) at
station n, hour t .
vn
elevation (meters MSL) RAW
of surface station n.
virtual temperature Step 2
geopotential (m2s~2) of
1000mb surface at locat'
x of station n, hour t
~n n
station pressure
(mb)
RAW
on
time rate of change
(m2s~3) of geopotential
of 1000mb surface at
station n location,
,
n
at hour t
.
m
p , (t ) sea-level pressure
Sl'n m (mb) at station n,
hour t
.
vn
see above
see above
RAW
Step 2
44
0 (t )
n m
p (t )
n m
surface potential
temperature (°K) at
station n, hour t
m
surface air density
(kg m'3) at station n,
hour t
ffl.
-------
YYYYYYYYYYY
H
3
CO O.
Q.H
D
^
4
^
<:
••»
*
•^
c
3
^
4-
^
<2
^
•»
*
^>
C
»
JM
<>
««i
c
^
r
^
c
r
\
•mv
1
I
GC
^
••
•
\
*|
^
£
h
«
^
c
D
•
«•
^
H
/
t 4>
^ •••
C
>
• •€
(
*
••»
c
0
>
y
(
\
(
*
•^
~€
1
>
^
•*
c
o
^
^
•
f
M>
•
I
*••••*
4-» —
"c i
3? c
j
V
CO >
Q.
LU
«•.
1^
»•
C
1
>
^
*•
*
^
G
"K
•
•*,
7
^
c
)
*m
7
•M
(
^
l_
;
CO
c
G
Q.
CO
C
oc
c
H-
C
Q
h-
c
a
c
_
o
V)
o
- o
•B I
5 ®
.2 c
(0 o
i3 (A
II
m
•£- 4-1
o -3
(/) ^
CO
03
3
a>
iZ
V V V V v V
45
-------
SECTION 5
PROCESSOR P4
GENERAL DISCUSSION
This processor estimates the surface roughness z , the Obukhov length L,
the surface heat flux Q, and the friction velocity u* in each cell of the
NEROS grid. The last three parameters are treated in a single processor,
rather than distributed among several, because they are interrelated variables
and the value estimated for one can be altered by a change in the method used
to estimate another. The method we use in this first generation processor
network to estimate Q is based on the scheme proposed by Golder (1972). More
refined methods have recently been reported, e.g., Holtslag and van Ulden
(1983), and these should be considered in the development of the second
generation model.
We should also point out that the estimates of the friction velocity u*
that we outline here are derived from the raw surface wind observations, i.e.,
(u,v), rather than from the flow fields that are finally used in the model,
i.e., the output (ls j) of Pll. To utilize data from the latter source
would result in complex interconnections of Pll with other processors that
require L, Q and u* values. The difficulty of operating such a system cannot
be justified considering that the method of estimating u* is itself quite
crude, and that the flow fields that Pll generates are constrained locally in
space-time by the observed winds. That is, the observed surface winds that we
will use here to estimate u* are explicitly incorporated into the flow fields
generated by Pll.
46
-------
Each of the steps below must be performed for each cell in the NEROS
grid.
Step I
Determine the wind speed and exposure classes C and C for use in
W c
estimating the Obukhov length L. First, using the locations XR of the N
surface meteorological stations (n=l,2,...N) and the observed winds (u ,v ) at
each station and at each hour tm, estimate the wind speed u in the given
grid cell at hour t by performing a weighted average of the observations
_2
(un,vn) at the nearest sites xn. The r weighting formula (given by Eq.
11-88) may be used. The wind speed class for this grid cell is now defined to
be
Cw = '-2-, if 0 < u < 8 m s"1; (4-1)
4, otherwise
Next, using the fraction o>T (x,t_) of local sky coverage by all cloud types
u i ~ m ~—
and the local land use distribution T(x,j), where x refers to points on the
NEROS grid and j=l,...10 refers to the 10 land use types (see P15), compute
the exposure class C :
3, aCT < 0.2
2, 0.2 < aCT < 0.7
1, 0.7 < aCT < 1.0 (4-2)
0, aCT = 1.0
-1, afT > 0.5
\ nighttime hours only
.-2, arT < 0.5
47
-------
In order to account for the nighttime heat fluxes from cities, we will limit
the minimum value of C to 1 for those cells with "urban" land usage (j=l)
greater than 30%, i.e.,
if T(x,l) > 0.30, Ce(x) > 1. (4-3a)
Similarly, to account approximately for the effects of large bodies of water
on the surface heat flux, we assume
-I, daylight hours
if T(x,7)>0.6, C (x) = \ (4-3b)
e '0, night.
Now compute the Obukhov length L in the cell at x at hour tm by
L = I (anS + a <53) • zn 1 2 3 I , (4-4)
where
al = 0.004349,
a2 = 0.003724,
bl= 0.5034,
b2 = 0.2310,
b3 = 0.0325.
The stability function, S (a digital version of the Pasquill stability category),
is expressed as
- C + CQ ) • Sign(CJ ,
w
e
48
-------
where
Sign(Ce)
1, if Ce > 0;
0, if Ce = 0;
l, if C < 0,
Step 2
Estimate the effective surface roughness z in each cell of the NEROS
grid using the following expression
10
10
x) = 2 T(x,n)z (n)/X T(x,n)
n=l
n=l
(4-5)
where T(x,n) is the fraction of the NEROS cell centered at x that is in land
use category n, with n (=1,2,,.. .10) referring to the following use types:
1. Urban Land 6.
2. Agricultural Land 7.
3. Rangeland 8.
4. Deciduous Forest Land 9.
5. Coniferous Forest Land 10.
Mixed Forest Land (including Forested Wetland)
Water
Land Falling Outside the Study Area
Non-Forested Wetland
Mixed Agricultural Land and Rangeland
In Eq. 4-5 the surface roughness zQ(n) associated with each land use type are
the following (based on the values suggested by Sheih et al., 1979 with
modifications as noted below):
(2)
(3)
(4)
(5)
= 0.7 nT
= 0.2
= 0.1
= 1.0
= 1.0
49
-------
(6) = 0.7 (*2)
(7) = 0.05
(8) = 0.00(*3)
(9) = 0.3
(10) = 0.15(*4)
Notes:
*1. Sheib, et al. (1979) recommended a value of z =lm for "metropolitan
city." In our case the "urban" category includes all "built-up"
lands, including residential, industrial, commercial, and other
areas characterized by building heights much lower than those
characterized as metropolitan. Our choice of 0.7 m for urban areas
is an estimated mean value.
*2. The value for Category 6 represents a rough average between the
values suggested.
by Shieh, et al. for marshland and ungrazed forests.
*3. None of the cells in the NEROS grid fall in this category, i.e.,
T(x,8)=0 everywhere, and hence the value assigned to zQ(8) is
immaterial.
*4. Category 10 is a rough average of cropland and rangeland values.
The values of z computed from Eq. (4-5) for each NEROS cell should be
recorded in the PIF.
50
-------
Step 3
Compute the local friction velocity u*(x, t ) using
u* =
G(z,L,z0)
(4-6)
where k = 0.4 is the von Karman constant; z is the elevation of the surface
wind observations above ground—use
2=10m;
(4-7)
and G is the similarity function
G = In
z+z.
if 1/L = 0
(4-8a)
G = In
z+z
MIN(5.2r,5.2) , if z < 6L and L > 0
(4-8b)
(4-8c)
The Obukhov length L is from (4-4), ZQ is the local surface roughness from
(4-5), and
-1 -l (M)(40+l)
= 2(tan x 4 - tan 4Q) + In [/fe+i\/fe -isL if L < 0.
4 =
Z+Z
40 = U-15-)
51
-------
Step 4
Estimate the effective surface kinematic heat flux in each grid cell at
hour t :
where T is the surface virtual temperature (°K) in the given cell obtained by
weighting the nearest observations T of virtual temperature (provided by P3)
_2
as in Step 1 above; and g = 9.8 ms is the acceleration due to gravity.
(Note that the T data from P3 are in degrees Celsius.)
Table 4-1 summarizes the input and output variables in each of the three
steps that comprise this processor, and Figure 4-1 illustrates the processor
and its data interfaces.
52
-------
Table 4-1 Summary of input and output parameters
in each step of Processor P4.
Input Output
Variable Description Source Step Variable Description
"n^m) east-west surface , P3 1 L(x,t ) Obukhov Ten
gth (meters)
wind component (ms )
at surface weather
station n, hour t .
north-south surface
wind component at
station n, hour t .
location of surface
weather station n
land use fraction
of category j in
NEROS cell centered
at x.
CFp,-(x,t ) fractional sky
coverage, total of
all cloud types over
cell centered at x
at hour t .
P3
PIF
PIF
in NEROS cell centered at
x, at hour t .
~ m
|LI (x,tm)| wind speed (m/sec) in NERI
~ cell centered at x, at hoi
t (for use in thTs Proce:
sor only).
nr
PIF
T(x,j) see above
PIF 2 z (x) effective
roughness
NEROS cell
at x.
surface
(m) of
centered
u (x,t )Iwind speed
r*t *** 111
m
L(x,t ) Obukhov length (m)
~ m
z (x) surface roughness
~
Step 1 3
Step 1
Step 2
friction velocity
(m/s) in NEROS
cell centered at
x, at hour t.
T,,~(x,t.J surface virtual
\/ ri />^ ill . . ~ ** _.
vn
m
temperature (°C)
at surface weather
station n, hour t .
P3
surf ace-, heat flux
(m °Ks"1) in NEROS
cell at x, at hour
^•w *
V
53
-------
Table 4-1 Summary of input and output parameters
in each step of Processor P4. (Concluded)
Input Output
Variable Description Source Step Variable Description
u.,t(x,t ) friction velocity Step 3
~ m (m/s) in cell at
x at hour t .
~ m
L(x,t) Obukhov length Step 1
~ ra (m) in cell at x
at hour t .
54
-------
YYYY
T3
CO
<**
3
a
CO
CL w
i> O
O CO
« %
CO V
80
2 S.
Q.
-------
SECTION 6
PROCESSOR P7
MOTION OF A VISCOUS, HYDROSTATIC FLUID OF CONSTANT DENSITY OVER IRREGULAR
TERRAIN.
We are concerned here with nighttime flow regimes, where winds in the
cold air adjacent to the ground, i.e., the radiation inversion layer, are
influenced by buoyancy, terrain, warm cities, geostrophic forcing and
friction.
We will treat the cold, radiation inversion layer as one of constant
density pQ whose upper surface is described by
Hvs(x,y,z,t) = zvs(x,y,t) - z = 0 (7-1)
The subscript "vs" designates that this is also the "virtual surface" of the
atmosphere above. We adopt this point of view later in formulating processor
Pll which describes the flow field above the inversion layer at night. Figure
7-1 illustrates H and other terms that we shall use in the following
analyses.
Terrain elevation (MSL) is represented by zt(x,y) and z is the vertical
coordinate whose origin is at seal eve!. In keeping with the level of
simplicity adopted in formulating the regional diffusion model, we shall treat
the cold air layer as a two-dimensional fluid. That is, the horizontal
velocity in the cold layer u=(u,v) is assumed to be invariant with
56
-------
respect to z. Keep in mind that during nighttime hours the cold, inversion
layer considered here is by definition Layer I of the regional model. The
flow speed in the atmosphere above the cold layer will be represented by
U =(U ,V ) (see Pll): and the density in the layer above the inversion will be
~m m m
represented by p , where P and Fxf represent
the x-components of the pressure, Coriolis and friction forces, respectively,
exerted on the parcel. Consider first the form of pressure force.
57
-------
The pressure force
We consider the pressure forces acting on the fluid parcel shown in
Figure 7-2. Since we assume the flow speed to be uniform in z, we take the
parcel to be a vertical column of fluid extending from z. to the top surface
u
zys of the cold layer and having horizontal dimensions (Ax,Ay).
Since the cold fluid and the air above are assumed to be in hydrostatic
balance, within the cold layer we have
3p/3z = -pQg.
zvs
(7-3)
Figure 7-2. Air parcel considered in the force balance analysis.
58
-------
And since p is assumed to be constant, at any level z within the parcel the
pressure is
p = pvs + pog(zvs " z) (7"4)
where p is the pressure at the top of the parcel, that is on H .
The total force on the left face of the parcel due to the hydrostatic
pressure is
zvs
fXL = Ayfpdz.
zt
Substituting (7-4) into this expression and making use of the definition
h(x,y,t) = zvs(x,y,t) - zt(x,y) (7-5)
we get
(7-6)
Keeping in mind that the pressure force is exerted inward on all faces, the
force on the right face of the parcel is
9p
fXR = " (pvs + ~^ AzVS)(h+Ah)Ay " P09Ay(h+Ah) /2 (7-7)
dZ
The force exerted by the ground on the parcel is directed normal to the
lower face which is inclined at an angle 0. with the horizontal plane and it
has a magnitude equal to the total fluid pressure force exerted by the fluid
on the ground, i.e. ,
59
-------
fxt = -ft sin 0t
= - (P09h + Pvs)(Ax2 + AZj.2)3* Ay sin 9t (7-8)
= - (p0gn + PVS)AV Azt
In a similar manner we find that the force on the top face of the parcel is
fxvs = PVS^VS' (7~9)
On combining (7-6) - (7-9) and noting that 9p /8z = -p g we obtain the total
x-component of the pressure force on the parcel:
Fxp = fXL + fXR + fxt + fxvs
= Ay(Azvs - Azt - Ah)pvs + pmghAyAzvs (7-10)
- pQghAy(Ah + Azt)
We can obtain Ah from (7-5) whereupon we can reduce (7-10) to the final form
(for small Ax, Ay)
3z
Fxp = - (Ao)gh -^ AxAy (7-11)
where
*P = po"pm (7-12)
Now the mass m of the fluid parcel is simply
m = pQghAxay,
60
-------
hence
37
1 c _ Ap __ ys ,-, 10_v
i FxP - " f g IT • (7"13a)
By analogy the y-component of the pressure induced acceleration is
I AP Zvs
Later (see Eqs. 7-34 and 7-35) we will add force components resulting from
synoptic scale pressure gradients.
The friction force
The friction force on the fluid parcel is the result of momentum fluxes
across the parcel's boundaries. These fluxes are caused by molecular
diffusion and by sub-parcel scale velocity fluctuations, or turbulence. For
example, at the earth's surface the velocity must be zero, and thus momentum
is drained from the fluid in much the same way as heat is removed from a fluid
at a cold, constant temperature surface. In the case of our flow where there
is no heat transfer into the fluid, the no-slip ground surface acts to
transform the bulk kinetic energy of the moving fluid into internal heat
energy. The fluid can also be accelerated by influxes of momentum from the
atmosphere above.
Referring to Figure 7-3 we note that the viscous force on the left faco
of the parcel is
fvl = Ayht (7-14)
AL XX
61
-------
where TXX is the net flux of x-momentum in the x-direction across the left
face of the parcel. The force on the right-hand face is
fXR ~
- (t
xx 8x
3tvv 3z
Ax)(h + -R— Ax)Ay
r -i^YX Aw
TY«+-#" Ay
(7-15)
Figure 7-3. Friction forces on a fluid parcel of horizontal
dimensions (Ax,Ay) bounded by z and zt.
V 5 w
Continuing this analysis for each of the other faces and assuming that the
slopes of z. and z are small enough that the areas of the bottom and top
t V 5
62
-------
faces of the parcel are approximately AxAy, we get
xf
hAx
t Ax
Ax2)
(7-16)
In the limit as Ax and Ay become very small, the mass m of the parcel is
m = p hAxAy
(7-17)
and hence
F = - -
m xf p
3x 8y
(7-18)
t + t
xx h 3x yx h 3y
For the lateral stress components T and t we will adopt the gradient
XX jrX
transfer forms
P0Txx
P0V
K §H
xx 8x
V ay
(7-19)
(7-20)
where K and K are elements of the eddy vicosity tensor, to be defined
A/\ y X
later. We will treat the stress T^ on the upp^r surface of the fluid as a
£^
prescribed variable. The surface stress T will be given later.
63
-------
By analogy with (7-18) we have
IP =. _1
m yf p L 8x " 8y T h
(7-21)
Txy h 8x Tyy H 8y -"
with
-1 r - -K §¥ /-7-o-j\
p^ yy ~ yy §y' ^ '
The stress T^ is considered to be a prescribed variable; t . and T are given
by
Tzx = " P0U*COS 8 (7-25a)
Tzy = " Pou*s1n 8 (7"25b)
where U* is the friction velocity which we shall express in the form,
following Melgarejo and Deardorff (1974),
u* = CQ(u2 + v2)15; (7-26)
where Cn is the surface drag coefficient and 8 is the wind direction,
0 = tan"1 ^. (7-27)
From (7-25) - (7-27) we obtain
5* 4 = - Cj (u2 * v2)^ (7-28.)
64
-------
It 22
-± T;W = - c£ (u^ + v
p zy U
According to Melgarejo and Deardorff,
C§ = k2 [(ln(r*) - b)2 + a2]"1 (7-29)
o
where z is the local surface roughness, and a and b are functions of h/L,
where L is the Obukhov length. (Approximate forms for a and b are given later,
7-89.)
The Coriolis force
The Coriolis force is given simply by
Fxc = *fv (7-30)
i Fyc ' -fu (
where f is the Coriolis parameter (=2fi sin ((», where = latitude and fi = earth
angular speed of rotation = 2IT/24 hr ).
The momentum equations
On combining (7-2), (7-12), (7-18) - (7-20), (7-28a) and (7-30) we obtain
the equation governing the u component of the flow speed:
65
-------
au + ,£H + »/au - -
at + % * vay ' ax~ xx 2 yx
+ "• 3x + h 3x -1 3x + "• 3y +
If S7 VS
J^ vs 3v
h 3y 3x
3y 3x
For convenience we will express the shear stress at z in the form
where avs is a vertical exchange velocity scale at 2-2 and UM is the
j is a vci i> i i_a i CAV.IICIIIUC VCIULIOV o«_aic at, L. — L anu u»j
W VS M
x-component of the "observed" flow in the layer above z (see processor Pll,
Stage UV).
Substituting (7-33) into (7-32) we obtain after rearranging terms
§" + r,, 8K** ^ azvs-, 3u . , 9Kyx ^x ^vs, 3u
3t L 3x h 3x J 3x L 3y h 3y J 3y
- *xx ^ - Kvx ^ * -R CCD <"2 + ^ + C^ •
3x y 3y
(7-34)
|! + . ,
w UM + 3y 3^ + Ky
x i
zvs avn . i apsi
h 3y 3xJ p 3x
66
-------
Likewise we have
-. + , . - -.
at ax h ax J ax L ay h ay J ay
a? + *
Q ^p ' 3y u h aw M 3x 8y xy
ffvs §u, . 1
h 3x 3yJ p
In Eqs. (7-34) and (7-35) the last terms on the right-hand side represent the
accelerations caused by synoptic scale pressure gradients in the boundary
layer, and p , is the sea-level pressure.
Next we derive an equation for the virtual surface elevation z (x,y,t).
The fluid depth equation
Let Q be a point on the virtual surface (cf 7-1)
Hvs(x,y,z,t) = zvs(x,y,t) - z = 0 (7-36)
and let its coordinates at time t be (XQ, y , ZQ). If the surface is moving
with a velocity v = (u , vvgj wyg) at the point Q, then at the later time
t + At Q will have the coordinates
o ^
(xr ylf Zl) = (XQ -H uvsAt, yo + v^At, ZQ + w^At). (7-37)
67
-------
Since by definition Q lies on the surface HVS at all times we must have
Hvs
-------
AS =
(7-44)
Hvs = zvs(x,y.t)-z =
Ax"
Figure 7-4. Projection of the horizontal rectangle (Ax,Ay)
onto the surface H centered at point Q.
If the slope of z is small, say
vs
8x
3z
vs
«i
(7-45)
then (7-44) can be approximated by
ft 7
As * AxAy
, if 7-45 holds.
(7-46)
Let v, be the velocity of the fluid at the point Q on H . Then the
normal downward component of the fluid velocity relative to H at point Q is
_ -1
(7-47)
69
-------
Making use of (7-42) we can write (7-47) in the equivalent form
(7'48)
where
The total downward volume flux of fluid crossing As is vAs. Using
(7-46), (7-48) and
we have
dH
~dt
vAs = AxAy -~ (7-51)
Thus, the net downward flux or fluid per unit horizontal area is
"vs • TF <7-52>
or
= V (7-53)
where (u, v, w) is the fluid velocity on H .
~~ V S
An expression for nvs can be derived from the first law of
thermodynamics. Assuming that all sources and sinks of heat are on the
boundaries of the fluid, e.g., ra'diative cooling on HVS and heat transfer to
the terrain surface H., we can write the first law in the form
70
-------
de _ ae ae . ae . .ae _
- + u- + v + w-
(7-54)
where 6 is the potential temperature defined as
p R/cp
e = T(£) p
p
(7-55)
where P is a reference pressure (usually 1000 mb) and p is local pressure.
Following the procedure described in Part 1, Section 2 for obtaining layer
averaged variables, we get from (7-54)
+ V
3y h
•ar
.t "^
t cit~
vs
(7-56)
V-VH f
~ ~ vs
vs
= 0
vs
where TT denotes averaging over the terrain surface Ht, T~7 denotes averaging
over the virtual surface H ,
h(x,y,t) = zvs(x,y,t) - zt(x,y)
(7-57)
and
i ( rvs
<6> H r I edzda
Az,
and A is the area of a grid cell in the model
(7-58)
Look at the first term in brackets in (7-56). In analogy with (7-52) we
have
dH
(7-59)
71
-------
where n. is the downward fluid volume flux at the ground. Thus
t
where w'e1 is the kinematic heat flux at the ground.
Consider now the second term in brackets in (7-56). Since 3zt/3t = 0
_t
and since n. = 0 (there is no net flux of air through the terrain), we have
t
dH. t
= F = 0. (7-61)
The third term in brackets in (7-56) can be written with the aid of
(7-52) as
vs
~dH^ vs
*!T =envs (7-62)
and the fourth term becomes, using (7-52) and (7-53),
vs 32
[nvs - -§£*] (7-63)
where n and 2 are taken to be averages over an area A of a grid cell
surrounding any given point.
Substituting (7-60) - (7-63) into (7-56) and collecting terms we get
dH -, vs
•£ + ± [-wT0To - 0nvs + <6>HVS] = 0 (7-64)
where
dt = 3t "ax V3y
72
-------
Following Zeman (1979) we shall assume that within the cold fluid, i.e.,
the nighttime stable boundary layer, the potential temperature has a linear
variation with height, namely
z-zt
8 = 8h " (8h " 8o)(1 " ~fi ^ (7"66)
where 0. is the value of 6 at the elevation h above ground and 0 is the
potential temperature at the surface. At elevation h the turbulent heat and
momentum fluxes are negligibly small; so if we take this to be the elevation
of the virtual surface z , then the kinematic heat flux on H is simply that
V & V5
due to the motion of the surface itself. In this case we have
_ vs
envs =ehnvs. if nvs > o (7-6?)
Also, on integrating (7-66) in the manner of (7-58) we get
<0> = h (0h + 0Q) (7-68)
and hence
d,, d,, d,,
at h
Substituting (7-67) and (7-68) into (7-64) we obtain
Making use of (7-69) in this equation we obtain
Hvs = B(he-h) (7-71)
where
-. d
K (eh + V
73
-------
2wTern
d -
It is instructive to compare (7-71) - (7-72) with models of the stable
boundary layer developed by others. Earlier efforts have considered only
homogeneous surfaces with no mean vertical motion (w=0). In this case (7-71)
reduces to (see 7-53 and 7-57)
f£ = B(he - h). (7-73)
and d^/dt •* 8/3t. Nieuwstadt and Tennekes (1981) (NT) recently proposed
f£=B'(h;-h) (7-74)
where
(7-75.)
2
h; = c4 fG sina cosa (7.?5b)
where f is the Coriolis parameter, G is the geostrophic wind speed, a is the
angle between the geostrophic and the surface winds, and c. is a constant
whose value is estimated by Nieuwstadt and Tennekes to be about 0.15. Good
agreement was found between the predictions of this model and observations of
the boundary layer depth h.
On comparing (7-71) - (7-72) with the NT model (7-74) - (7-75) we find
that except for the constant 4/3 in B1, the two models are the same, provided
74
-------
that the surface heat flux satisfies
|wTgr ( = 0^5 fG2sina cosa (
To determine whether this is a plausible relationship, we note first that
for large values of h/L, where L is the Obukhov length
(7-77)
gk vF0J
Brost and Wyngaard (1978) found
G2sinacosa = u£ ()2 (*) (7-78)
where k = 0.4 in the von Karman constant. Substituting (7-77) and (7-78) into
(7-76) we find that (7-76) implies
h = 0.8 (^)% (7-79)
But this is just the formula that Zilitinkevich (1972) derived for h from
similarity theory [in particular h = Cg (%=)]; and thus we conclude that
(7-76) is a plausible expression, at least for large h/L.
Having an expression now for q , let us return to (7-53) and consider
the vertical velocity w. Since the cold layer is shallow we can approximate
the continuity equation by
Integrating from z. to z and using the assumption that u and v are invariant
with respect to z, we get
75
-------
Since there is no flux F of fluid through the ground, we see from (7-61)
that
9zt 8z. 3zf
8T + "55T * V - «t • ° V-™
Solving (7-82) for w. and using the result in (7-81) we obtain
u, - .8zt . t,8zt . rau . 8vn ,, 0,v
wvs - "alT + V8y- ' *% + 3y] (7"83)
where h is given by (7-57). Substituting (7-83) into (7-53) we obtain finally
where
- zt (7-84b)
This is the model equation that we shall use to obtain z .
Simplified model equations
Going back to the momentum equations (7-34) and (7-35) and recalling that
z is the elevation where vertical heat and momentum fluxes due to turbulence
are negligible, we assume that
avs~ 0. (7-85)
Let us also assume in this "first generation" model that the lateral eddy
76
-------
viscosity tensor K is also zero. With these assumptions (7-34) and (7-35)
reduce to
1 8psl
5* T*l - fv (7-86)
"ap
sl
where p , is the sea-level pressure and h and z are given by (7-84), which
we rewrite here for completeness:
h = zvs - ZT. (7-88b)
The set of equations (7-86) - (7-88) is a closed system that we can solve for
u, v, and h given CQ> Ap, p , p , , z,. and q . Earlier we described how nvs
is related to the surface heat flux w'e' and the temperature distribution 0
within the cold layer (see 7-71 and 7-72); and we presented formulas that one
might use to derive these variables (see 7-66 and 7-76). Below we summarize
the expressions that we propose to use for these and the other parameters
listed above in the first generation version of processor P7.
77
-------
(1) nvs = Eq. 7-71,72 with Q^ and 0Q from surface and upper air mete-
orological data (Processors PI and P3); vTe7 = Q = kinematic
heat flux at the ground (from Processor P4);
(2) Zy(x,y) = 30 x 45 min (1 at- Ion) smoothed terrain heights (meters,
v
MSL, from P7)
(3) CD = k{[ln(^) -b]2 + a2}'*5 (7-89a)
o
where k = 0.4 is the von Karman constant; h is the local
solution of Eq. (7-88); z is the surface roughness (m), from
P4;
T5, if L>0,
a =i
i'l, otherwise; (7-89b)
f5, if L>0,
b (7-89c)
4, otherwise;
and L is the Obukhov length (from P4). Note that L, z , h, and
hence a, b, and CQ are defined for each grid point in the model
domain. (Eq. 89 is based on values given by Melgarejo and
Deardorff, 1974.)
(4) : -. (7-90.)
where 6, and 6 are from Processors PI and P3;
(5) p0-fl (7-90b)
where p , is the sea-level pressure (from P3); Ty is the virtual
temperature (from P3); and R is the gas constant.
78
-------
_2
(6) g and f are gravity (9.8 m sec ) and the Coriolis parameter (=
2ftsin), respectively. Here Q = 2jt/(24-60*60) sec'1 and $ is the
local latitude in radians.
Later we outline the stages of calculations that are needed to compute the
parameters listed above and the additional quantities, not used in the flow
simulation, that P7 must provide to the model and to other processors in the
network.
Solution of the u, v, and z equations
Each of the Eqs. (7-86) - (7-88), which govern u, v and zys,
respectively, has the form
v + br = s <
where U and V are functions of F. We will approach the task of solving (7-91)
numerically using the technique developed in Section 9 of Part 1. That is, we
assume that within each small time interval At, the coefficients U and V in
(7-91) can be treated as independent variables whose values are determined
using the value of r at the beginning of the interval, t say. In this case
the solution of (7-91) can be expressed in the closed form
r(x,to+At) = fp(x,t0+At x',t0)r(xit0)dx'
(7-92)
VAt
p(x,t +At|xU')S(x',t«)dt'dx'
79
-------
where p is the Green's function of (7-91). In the present instance it has the
form
P(x,y,t|xjyit0) = 6[x-(x'+x)]6[y-(y'+y)]e"b(t"to) (7"93)
where 6[x] is the delta function and
t
:|x',y',t0) = U[x' + x(t'), y1 + y(t'),t']dt' (7-94a)
x =
t
= y(t|x',y',tn) =fv[x' + x(t'),y' + y(t'),t']dt'. (7-94b)
oj
•
t.
The assumptions that U and V are approximately constant during the time
step At should be valid provided that At is small enough to satisfy
IT] -5
At <
-------
exceeds unity, flow discontinuities such as hydraulic jumps will occur, and
these will cause considerable problems in the numerical model. We do not
expect the nighttime flows that we will simulate to become supercritical
often. When it happens, it will occur in isolated portions of the model
region and we will be able to anticipate it by monitoring the temporal
behavior of the flow. In those grid cells where we predict that the flow is
about to become supercritical, we propose to prevent it by applying enhanced
eddy viscosity.
A detailed description of the numerical scheme used to solve the flow
model equations, specifically equations reducible to the form (7-91), is given
in Appendix A of this section. The scheme described there is the same one
that is used to solve the transport and diffusion component of the regional
pollution model equations. The scheme is an explicit, 5th order space, 1st
order time algorithm that permits the model domain to be treated in piece wise
fashion. This is a particularly valuable feature in models such as ours that
are too large to load entirely into the computer memory.
One of the principal problems associated with the numerical solution of
equations like (7-86) - (7-88) is the treatment of the lateral boundary
conditions. Since no generally valid method exists for specifying the
boundary values required by the difference equations, a common practice in
mesoscale flow models is to place the boundaries far from the edges of the
spatial domain of interest. Another approach is to imbed the flow model
within another one of coarser resolution which provides boundary values.
81
-------
Neither of these methods of circumventing the boundary problem can be used in
our studies because the additional computer storage and time that they would
require would make the overall simulation effort impractical. Consequently,
we have formulated an approximation of the boundary conditions, described in
Appendix B of this section, that is sufficient for treating the limited
periods of concern to us in modeling the nighttime boundary layer flow.
Below we outline the various stages of computation that are necessary
within processor P7.
StaQe ZT
The "raw" topography data available in the PIF will be denoted here by
zt(\,). These represent terrain elevations averaged over 5 min x 5 min
latitude - longitude sectors. The regional model and the flow model developed
in this section require the elevation z.(A.,) averaged over 30 min x 45 min
latitude - longitude sectors.
Let zt(I,J) denote the value of the 30 x 45 min smoothed terrain at grid
point (1,0) (column I, row 0) of the NEROS region. Then
ztd,0)=£ z,(i,j) (7-97)
1 r
where the summation is over the 54, 5 min x 5 min cells surrounding grid point
(1,0) (see figure 7-5). Note that the 54 cell smoothing area used in the
definition of zt(I,0) overlaps the smoothing areas associated with the 8 NEROS
U
grid points nearest the point (1,0).
82
-------
The input and output requirements of Stage ZT are summarized in Table 7-1
later in this section.
®_.,
(I,J)
©_.
— ,j
M.J+1
®...
-
(1+1. J)
®, -.-
Figure 7-5 . Illustration of the 54 5 min x 5 min cells
that are used in the calculation of zt(I,J).
Stage DELRO
The parameter Ap, is used in this and other processors as an indicator of
the presence of an inversion layer at the ground. This parameter will be a
scalar and a function of time only in the first generation model. We define
1, if surface inversion is present
Ap,(t ) = ) over most of the model region
.0, otherwise.
(7-98)
Under this definition the magnitude of Ap1 is an indicator of the density of
83
-------
air in Layer 1 relative to that in Layer 2. In the first generation model we
will compute Ap,(t ) using the following procedure.
where
(1) Define
6Pl(tn) =
+1, if Ap
-1, if Ap-j^C^
0, otherwise
) = 0 and Q(tn) < 0;
= 1 and TD~EF(tn) <0;
Q(tn)=IIQ(i,j,tn)
(7-99)
is the surface kinematic heat flux (°Kmsec ) summed over all grid points (i,j)
of the model domain, and TDEF is the mean "temperature deficit" of the cold layer
which we approximate by
TDEF(t ) = allhL. d.j.t ) - 112 q(i,j,n)At.
n ... max ID ijn=0
(7-100)
In this expression h is the depth of the cold layer at the time the surface
heat flux reverses, i.e.,
h(i,j,tm), if Q(i,j,tm) > 0
and. QOJ.Vp < 0;
0, otherwise
(7-101)
The value assigned to h_,v before the surface heat flux reversal, in this case
uiaX
0, is immaterial since the value is never used. In (7-100), the variable q is
given by
!(U,tm), if Q(i,j,y > 0;
I, otherwise
(7-102)
84
-------
The constant a that appears In (7-100) is a temperature gradient that we
define in Stage ETA.
(2) With SpjCt^) determined by (7-99) a functional form for ApjCt ) that
is consistent with (7-98) is the following:
Apl(tn) = Apl(tn-l) + 6pl(tn} (7-103a)
with initial value
Ap-^) = 0 (7-103b)
and
t0 = 1200 EST. (7-103c)
Under definition (7-103) we assume in effect that an inversion layer forms
over the entire model domain at the hour the average surface heat flux over
the whole domain becomes negative; and it disappears everywhere at the hour t
that TDEF, defined by (7-100), first becomes negative. This is clearly a
crude approximation in a model such as ours which spans 15° of longitude, but
to relax it would require an escalation in the complexity of the flow field
description that would put the overall modeling effort beyond the scope of the
"first generation" effort.
In summary, stage DELRO computes the single scalar variable Ap,(t ),
n=0,...N using Eqs. (7-99) - (7-103). This variable is used throughout
processor P7 as an indication of when specific functions are to be performed,
and it is an output of P7 that guides the use of the fields generated here in
other processors in the model network, A list of all the input and outputs of
this stage is given later in Table 7-1.
85
-------
Stage ETA
This stage operates only during those time steps t when Ap-,(t )=1 and it
operates in unison with Stage FLOMOD, described below, which solves the
equation set (7-86) - (7-88). In this stage we compute nvs> using Eqs. (7-71)
and (7-72), and Ap/p , using (7-90a). Both of these calculations require
values of the temperature 6. at the top of the cold layer. Rather than
attempt to estimate these from upper air data, which would be a difficult task
given the shallow depth of the cold layer and the limited frequency of the
upper air measurements, we propose instead to estimate 0. from the observed
ground-level temperatures 6 assuming a constant temperature lapse
|f = a (7-104)
in the cold layer. The observations reported by Godowitch and Ching (1980) of
the nighttime surface inversion in rural areas around St. Louis indicate that
a typical value is
a = 1.2 °C/100m. (7-105)
Using (7-104) we have
0. = 0 + ah (7-106)
n o
with a given by (7-105); and hence from (7-90a)
_ ah ~ ah
" 26? (
We will estimate the surface temperature QQ at each grid point (I,J) using the
virtual temperature observations Tyn, n=l,...N made at the M surface weather
stations. We assume here that the T"vn(°C) are available from Processor P3 for
each hour. In this case (M.t) is computed as follows:
86
-------
N
T -2
N -2
I r
n=l n.
(7-108)
where
(7"109)
is the distance between grid cell (I,J) and the site (xn,y ) of surface
station n. Equation (7-108) should be evaluated at each time step t that
The inversion layer growth rate parameter n can now be computed from
(7-71) and (7-72) for each grid cell and time step. Since the physical
assumptions on which the governing equations (7-86) - (7-88) are based do not
hold in situations where the fluid is being heated from below, we must limit
the minimum value of n to zero. The net effect of this is to assume that
over warm cities at night, the surface heat flux is never positive. Thus, we
have
nvsd,J,tm) = max (0, B[he - h(I,J,y]} (7-110)
where
B = -
87
-------
and
2Q(I,J,t)
In these expressions we can use the approximations
8T2 - 9°(I>J>tm)At0(I>J>Vl) <7-112>
3h n(I.J,tm) - h(ItJ>Vl)
at At (7"113)
where 6 is given by (7-108) and h is the solution of Eq. 7-88 provided by
Stage FLOMOD.
Note in (7-lllb) that we use estimates of the surface heat flux Q in each
cell rather than the expression (7-76) discussed earlier in the estimate of
h . The latter is based on the assumption of homogeneous, flat terrain and
homogeneous heat flux and hence estimates of w'6 ' derived from it using
observations of the ageostrophic wind angle a acquired from wind observations
would probably be erroneous.
In summary, Stage ETA provides values of Ap/P0 (using (7-107)) at each
grid point and time step that Ap-, = 1; values of nvs at the same locations and
times using (7-110); and the surface temperature 0Q at all hours. A list of
the input and output parameters of this Stage is provided in Table 7-1.
-------
Stage PCD
This stage computes the pressure and drag coefficient terms that are
required in the momentum equations (7-86) and (7-87). Therefore, this stage
operates only at times t when Ap-iCO = 1.
We assume that the sea- level pressure data p , measured at each of the
s I ,n
n=l,2,...N surface weather stations is available for each hour t that PCD
operates. At each of these hours sea-level pressure values should be
interpolated at each grid point in the model domain using the following
formula:
1 P,l.n>si»'J>V • —^
I r"1
n=l"
where
rn = [(IAX - xn)2 + (JAy - y^2]55 (7-115)
and N is the number of surface stations at which sea- level pressure values are
available. Using the same interpolation formula, derive estimates of the
surface air density at each grid point at hour t :
N
S
n
n=l n
where r is given by (7-115) and p is the density at surface station n (from
P3).
89
-------
At each grid point we require values of the functions
1 apsl
P £ - (7-117a)
1 9psl
py ' ^ TT (7-117b)
We will approximate the spatial derivatives that appear in (7-117) using the
fourth-order finite difference scheme discussed in the appendix to Section 8.
With those expressions we get
Px (I.J.y = A[p0(I,J,tm)]"1{2/3[psl(I-H,J,tm)
- ps1(I-l,J,tm)] - l/12[psl(H-2,J,tm) (7-118)
- Ps1(I-2,J,tm)]}(Ax)"1
Py(I,J,tm) = A[p0(I,J,tm)]"1{2/3[Psl(I>J-H,tm)
- psl(I,J-l,tm)] - l/12Cpsl(I,J+2,tm) (7-119)
- Psl(I,J-2,tm)]}(Ay)"1
where Ax and Ay are the x and y separation distances (m) of the grid points
nearest (I,J); and A = 100 is the conversion factor required to transform the
-1 ~2
sea-level pressure values from units of millibars to kg m sec '. (Note that
_o
p has units of kg m and Ax and Ay have units of meters.) The gridded
functions P and P are outputs of Stage PCD each hour that Ap-, = 1; and the
x y i
sea-level pressure p , given by (7-114) is an output for all hours.
The computation of the drag coefficient value Cg in each grid cell and
hour requires a straightforward implementation of Eq. (7-89). We have
90
-------
hd.J.t ) 2 2 _h
CD d,J,tm) =k{[ln ( z (I>J^ )-b]2 +a2} * (7-120)
where k=0.4; z is the surface roughness (m) of cell (I,J); h is the depth (m)
of the surface inversion layer in (I,J) at t (from Stage FLOMOD); and
5, if L(I,J,t ) > 0;
a = i m
-1, otherwise
-5, if L(I,J,t ) > 0;
b = { . m
4, otherwise.
Here L is the Obukhov length (m) in cell (I,J) at hour t (from P4). The
coefficient C.,(I,J,t ) is the third and final output of Stage PCD. (Refer to
Table 7-1 for a summary of the inputs and outputs of this Stage.)
Stage IBC
This stage computes boundary and initial values of u, v, and h for use in
solving the system of equations (7-86) - (7-88). For this purpose we assume
that initially and along the boundaries of the simulation region the velocity
field is in equilibrium with the friction, pressure and Coriolis forces. That
is, we assume (cf 7-86, 87)
2r2
—2 = -Pv + fs (sine) (7-121)
9 9
^CoSin 6
= -P - fs (cose) (7-122)
91
-------
2 2 ^
where s = (u + v )^ is the flow speed and 0 is its direction. Cross-
multiplying (7-121) by (7-122) we get
P¥ sine - Pw cose = fs (7-123)
* y
and on squaring (7-121) and (7-122), adding the results and Baking use of
(7-123) we obtain
s4 4
+ fV - (P/ + Pv/) = 0. (7-124)
x y
Since all parameters except s in (7-124) are known, we can solve this equation
for the flow speed s at each grid point and then substitute the results into
(7-123) to obtain the corresponding flow directions 0.
Next, we define the time tKQ:
tKQ = time during a given 24- hour period when
Ap.^ changes from 0 to 1. (7-125)
This marks the initial instant at which the equations (7-86) - (7-88) apply,
and it is the function of Stage IBC to provide the necessary values of u, v,
and h at this time. Thus,
if t = t,,n, Stage IBC computes initial values as follows:
Ml 1\U
h(i,j,tKO) = h0 = 30m, (7-126)
u(1,j,tK(J) = s(i,j,tKO)cos0(i,j,tKQ) (7-127)
v(i,j,tKO) = s(i,j,tKO)sin0(i,j,tKO) (7-128)
92
-------
where s(i,j,t.,Q) is the solution of (7-124) at grid point (i,j) based on
values of P , P , Cn and h (from 7-126) at that point; and 6(i,i,t,.,n) is the
X y U I\U
corresponding solution of Eq. (7-123).
if t > tj,Q and Ap-,(t ) = 1, Stage IBC computes
boundary values of u, v, and h as follows:
°Bc(i'J'V = AT MU.V - "(u.Vi)] (7-130)
In these equations (i,j) are boundary points only. Also, n is given by
V 5
Stage ETA and u and v are computed on the boundaries using the same method
employed to obtain the initial conditions, (7-127) and (7-128) above. We
should point out that in our treatment of the boundary conditions (see
Appendix B to this chapter), we give the values of u, v, and h at inflow
boundary points in terms of their initial values, i.e., h(i,j,tKQ) etc., and
we require that the subsequent time derivatives of these variables be
specified, namely
H, U, and V.
Stage IBC provides the boundary conditions (7-129) - (7-131) at each time
step t that Ap, = 1. The input and output parameters for this stage are
summarized in Table 7-1.
93
-------
Stage H1HO
This stage performs the final operations in the calculation of the
elevation z,, of the top surface of Layer 1 of the regional model; and it
converts this and other surface elevations into pressure coordinates. The
specific operations are defined below.
rn(i,j,t ), if Ap,(t ) = 1;
\ m I m (7-132)
VU'VH
(_max{100, min{500, L(i,j,tm) }},
if Ap,(t ) = 0.
Here h is the solution of (7-88), as provided by stage FLOMOD,, and L is the
Obukhov length, used earlier. The value assigned h, by Eq. (7-132) when Ap-, =
0 is constrained to prevent the Layer 1 depth from approaching zero
thickness — which could happen in extremely unstable conditions when L •*
0 — and also to prevent Layer 1 from becoming as deep as the mixed layer —
which would force the thickness of Layer 2 to zero. Although these
constraints are not necessary in principle, they are applied to prevent
numerical problems in the code which might arise as cell thicknesses approach
arbitrarily small values.
The virtual surface elevation zvs is computed as follows:
f zT(i,j) + MlJ.t ), if Apx (t ) = 1;
zvs(i,j,tj= T l m l m (7-133)
VS m I zT(i,j), if AP-L (tm) = 0.
Since ZT is given in meters above sea-level, the values of z1 and zys obtained
from (7-132) and (7-133) are also in units of m (MSL).
94
-------
Several other processors require the surface elevations in pressure
coordinates. These are computed as follows:
Zi(1,j,t )g
- 1 (7"134)
RdT(U,tm)
where g = 9.8m sec"2, Rd = 287m2 sec"2 °K"1, z^ = ZT + h-^ and
i,j,tm) + 0.0065zt(i,j)] . (7-135)
In this expression 0.0065 °K m is the temperature lapse rate of the Standard
Atmosphere which, for altitude calculation purposes, we assume holds between
ground elevation and sea-level.
Similarly
Pvs(i,j,tm) = Eq. (7-134) with ^ replaced by zvs (7-136)
Finally, in accordance with the analyses presented in Appendix B of Part
1 of this report, we shall prescribe the depth h of Layer 0 to be
no(i,j,tm) = h^iJ.y/lO. (7-137)
The inputs and outputs of Stage H1HO are summarized in Table 7-1.
Stage SIG
This stage estimates the fractions O--Q and a-j-g of surfaces zi(x«y>tm)
[Zj(x,y) + h (x,y,t )] that are penetrated by terrain.
Referring to Fig. 7-5 which illustrates the process of computing Zj in
Stage ZT, we define
95
-------
o-T1(I,J,t ) = -i 2-\(1,j,t ) (7-138)
11 m b4 i,jeD(I,J)m
where \ is defined as
l,if Mi.j) ^ Zi^.J.tn,); (7-139)
C -L III
0, otherwise
and the summation in (7-138) is over all 54 of the 5 min x 5 min subcells
contained within the 3 x 3 grid cell area in which Zj (I,J) is defined.
For simplicity in this first generation model we will assume
o-TOd,J,tm) =aT1(I,J,tm), (7-140)
where ov., is given by (7-138). The inputs and outputs of Stage SIG are
summarized in Table. 7-1.
Stage FLOMOD
This stage solves the system of equations (7-86) - (7-88) using the
numerical method discussed earlier, beginning on page 79. A detailed
description of the numerical procedures used to solve the equations is given
in the appendices to this chapter; and the input and output parameters of this
stage are summarized in Table 7-1.
A schematic view of the interrelationships among the various stages that
comprise Processor P7 is given in Figure 7-6, page 130.
96
-------
Table 7-1. Summary of the input and output requirements
of each stage of Processor P7.
Input
Variable
Description
Source Stage
Output
Variable
Description
QdJ.t )
surface kinematic -.x
heat flux (°K m s"i;
at grid cell (i,j)
at hour t_
P4 DELRO Ap1(t|n)
surface inversion
indicator: = 1 if
an inversion is
present over the
entire model domain
at hour t • = 0
otherwise.
depth of cold layer Stage
(m) at grid point FLOMOD
(i,j) at hour t
5 min x 5 min
smoothed terrain
elevation (m MSL)
centered at longitude
A, latitude .
RAW
ZT
30 min x 45 min
smoothed terrain
centered at NEROS
grid cell (i,j).
virtual temperature P3
(°C) at surface
weather station
n=l,... N at hour
V
depth of cold layer Stage
(m) at grid point FLOMOD
(i,j) at hour t
ETA Ap/p (i,j,t ) buoyancy parameter at
(operates grid point (i,j) at
only when time t
nvs(i,j,tn)
e0(U,tn)
cold layer growth
rate (m s"1) at grid
point (i ,j) at time
V
Surface temperature,
Q(i,j,t )
m
surface heat flux
(°K m sec"1)
P4
"sl,n
-------
Table 7-1. (Continued)
Input
Variable
Description
Source Stage
Output
Variable
Description
KI.J.V
z0(I,J)
depth (meter) of
cold surface layer
in cell (I,J) at
hour m.
surface roughness
(m) in cell (I,J)
Obukhov length (m)
in cell (I,J) at
hour m.
Stage PCD Cn(I,J,tJ
FLOMOD (cont.) u m
P4
P4
Psl(i,j,tm)
drag coefficient
(dimensionless)
in cell (I,J) at
hour m.
sea-level pressure i
cell (i,j) at hour t
Pv(i>j»t ) horizontal pressure Stage IBC
m
force term (see 7-118) PCD
horizontal pressure Stage
force term (see 7-119) PCD
h(i,j,tKO)
'KO
initial depth (m)
of the cold layer in
cell (i,j)
u(i,j,t,,n) initial east-west
flow speed (m sec )
component in cell
CD(i,j,tm) drag coefficient
nwcO',J>t_) cold layer growth
vs m rate (m/sec)
z-rO'.J) elevation (m, MSI)
of smoothed terrain
Stage
PCD
Stage
ETA
Stage H1HO
ZT
v(i,j,tKO) initial north-south
flow speed component
in cell (1,j)
Hor-O'jj.t ) time derivative (m/s
BL m of cold layer depth
at boundary cell
(i,j) at time tffl
URr(i,j,t ) time derivative of
BC m the east-west flow
component on the
boundary at time t
VRr(i,j,t ) time derivative of
BU m the north-south flow
component at boundar
point_2(i,j) (units:
m sec )
h-,(i,j,t ) thickness (m) of
1 m Layer 1 in cell (i.j!
in grid cell (i,j)
at time t
98
-------
Table 7-1. (Continued)
Input
Variable
Description
Source Stage
Output
Variable
Description
hd.j.V
L(i,j,tm)
depth (m) of cold
surface layer
Obukhov length (m)
sea-level pressure
(mb) in cell (i,j)
at hour t
surface temperature
(°K) in cell (i,j),
hour t
m
Stage H1HO
FLOMOD (cont.)
P4
Stage
PCD
Stage
ETA
,
m
elevation (m, MSL)
of virtual surface
in cell (i,j), hour
m
elevation of surface
z, in pressure (mb)
coordinates.
elevation of the
virtual surface z
in pressure (mb)
coordi nates
indicator of
presence of cold
surface layer
Stage
DELRO
depth (m) of Layer
0 in cell (i,j) at
hour t
.
zf(i,j) smoothed terrain Stage SIG
elevations (m, MSL) ZT
z+.0',j) 5 min x 5 min averaged RAW
terrain elevation
(m, MSL) of cell (i,j)
elevation (m, MSL) Stage
of top surface of H1HO
Layer 1.
aT1(i,j,t ) fraction (0 - cell averaged east
west flow speed
component (m/sec)
at time t in the
cold layer1
99
-------
Table 7.1. (Concluded)
Input
Variable
Description
Output
Source Stage Variable
Description
cn(i,j,tm) drag coefficient
Ap/p (i,j,t )buoyancy parameter
Stage
PCD
Stage
ETA
same as ,
except nortn-sou
component
Ap,(t )
cold layer growth Stage
rate ETA
indicator of presence Stage
of cold surface layer DELRO
u(i,j,t,,n) \ initial values of
v(i,j,t™) ) h, u, and v
Stage
IBC
) time rate of ,
m change (m sec"1)
of cold layer
depth at boundary
point (iij1) at
time t
m
m -2
) acceleration (m sec )
of east-west flow
component at boundary
point (iij1) at time
V
Stage
IBC
Stage
IBC
m'
same as URr except Stage
north-south flow IBC
component.
100
-------
Appendix A to Section 7.
Here we describe a numerical procedure for deriving solutions of
differential equations of the form
or* sir* 5ii"*
°L + LpL + y°L + t>f =
where all coefficients are functions of time and all except K and b are
functions of (x,y,t). We pointed out in Eq. 7-92 that the solution of (7-A1)
can be expressed in the form
Cx,t, I xiOrCxlOdx1
10 ° (7-A2)
1
pCx.t-JxJt^SCxlt^dt'dx'
t
where
p(x,t|x',t') =^?exp [-^^; - ^ ^> -) b(t")dt"] (7-A3)
t1
x = x(t | x'.yit1) = UCx'+xCt'^.y'+yCt'^.t'^dt" (7-A4)
t
Y = y(t|xiy',t') = rv(x'-Hx(t"),y'+y(t"),t")dt" (7-A5)
t1
101
-------
and
a2 = 2K(t-t'). (7-A6)
In both our flow model and in the regional model we are interested in the
values of the dependent variable r only at the grid points
and only at the discrete time intervals t =nAt, n=l ..... Furthermore, in the
situations of concern to us the spatial variations in S are of a scale much
larger than U...wAt or V^/wAt and the temporal variations are generally
slower than the time step At. Under these conditions we can express (7-A2) in
the approximate form
where
p(lAx,JAy,tn+1|x;tn)F(x;tn)dx' (7-A7)
r(lAx,JAy,tn) (7-A8)
and
FQ$;tn) = F(xjtn) + S(xJtn)At. (7-A9)
Eq. 7-A7 expresses the value of f at grid point (I,J) at the future time t -,
in terms of its known values at the present time t . Since the kernel p has
the form (7-A3), we can evaluate (7-A7) analytically if we express F(x' ,t ) in
polynomial form (or in a Fourier series).
To do this we note first that the kernel p in (7-A7) has a maximum value
at the point
x' = x* = lAx-x 1 (7_A10)
JAy-y j
102
-------
and it decreases to zero rather rapidly away from this point. In fact, if
K=0, p has the delta function form given earlier by (7-93). Thus, the
polynomial that we use to represent F in (7-A7) should have maximum accuracy
in the vicinity of the point x' = (x*,y*), which we can find by solving the
transcendental equations (7-A10). In the simple case where the spatial
variations in U and V are much larger than Ax and Ay and the temporal changes
in U and V are slow compared to At, (7-A10) yields the approximate solutions
(cf 7-94)
x* = IAX - U(lAx,JAy,tn+1)At 1 (7.An)
y* = JAy - V(lAx,JAy,tn+1)At j
The points (x*,y*) and (lAx.JAy) are illustrated in Figure 7-A1.
Using computer programming notation (to facilitate comparisons of the
theory presented here with the actual computer code), we define
1ST = [IFIX(xVAx) + 0.5]Ax (7-A12)
JST = [IFIX(yVAy) + 0.5]Ay (7-A13)
As illustrated in Figure 7-A1, (1ST, JST) is the grid cell center closest to
(x*,y*) (taking the grid points (I,J) to lie at the corners of each grid
cell). In preparation for the expansion of F(x' ,t ) in polynomial form, we
i«v p|
define the new coordinates
(7'A14)
C'-AIS)
Note that the origin of the (n.,4) system is the cell center (1ST, JST). We can
103
-------
now express F(x',t ) by the biquintic (5th order) Lagrange polynomial
Ffx't ^ = 2J z TF
C~J n; XL j=i L ij
n -
n
k=l
1=1
(7-A16)
O O O O
O O O O,
0 0 •
/a o •
o o o o/ o o
o o o
o o o
o o o o o o
/ Back trajectory starting
from (I,J) at time t +-,,
going back in time
totn.
GRID POINT (I,J)
x*,y;(c)
IST.JST)
Figure 7-A1.
Illustration of the points used in the numerical
solution of Eq. 7-A1. Circled grid points are those
from which values of F and S are taken to derive a
biquintic expansion of F(x',t_) about the point
(IST.JST) (see Eq. 7-A9).~ n
where
a2 =
a3 =
a4 =
a5 =
a6 =
(7-A17)
-n _
F(iAx/2+IST,jAy/2+JST,tn)
(7-A18)
104
-------
and
6 ,
J[ ( ) = product over al1 k except k=i,
k=l
6 ,
]~[ " ( ) = product over all 1 except l=j.
1=1
Note that (a.,b.), 1,j=l...6 are the coordinates in (n,4) space of the grid
points at which the F.. are evaluated. These are the circled grid points
shown in Figure 7-A1.
To simplify (7-A16) let
6 ,
X,- =F1 (n-ak) (7-A19)
1 k=l *
6 ,
Y, =IT (4-b,) (7-A20)
J 1=1 '
6 ,
L.. =H (a,-a.) (7-A21)
1 k=l n K
6 ,
M. =f[ (b.-b,). (7-A22)
J 1=1 J '
The last two parameters can be evaluated directly using (7-A17). We get
Lj = MX = -3840
L2 = M2 = 768
L3 = M3 = - 384 { (7-A23)
L4 = M4 = 384
[_5 = M5 = - 768
L, = Mc = 3840
b b
105
-------
Now (7-A16) can be written
6 6 n X,Y
Now look at the expansion of X.:
Xx = n5 - 5q4 - 10n3 + 50n2 + 9n - 45 (7-A25)
X2 = n5 - 3n4 - 26n,3 + 78n2 + 25n - 75 (7-A26)
•
•
•
Xg = n5 + 5q4 - 10n3 - 50q2 + 9n + 45. (7-A27)
Thus, we can express the X. in the alternate form
X. = Z a. na (7-A28)
1 a=0 la
where the coefficients a^ are the known values given in (7-A25 - A27).
Similarly
Y. = 1 a.ft4P. (7-A29)
J «=n JP
Substituting (7-A28) and (7-A29) into (7-A24) and we obtain
5 5
F(x'ftn)=I Z A"na4P (7-A30)
n a=0 p=0 afi
106
-------
where
6 6
n n a. a.
...
We now have the polynomial expansion of F(x' ,t ) that we need to evaluate
(7-A7), but first we must change the integration variables in (7-A7) from x1 to
(n,C). We get
00
-00
F(n,4,tn)dnd4
where
2a2
(7-A34)
(see Eq. 7-A3); and b* = b(x*,y*,tn).
Let
e = x* - 1ST (7-A35)
/\
e = y* - JST. (7-A36)
Then, representing x1 in (7-A33) in terms of n we obtain (from 7-A14 and
7-A10)
>x = -=~— exp C- -352*- ] (7-A37)
107
-------
where
sx = ex/(Ax/2) (7-A38)
CT = a/(Ax/2) (7-A39)
An expression similar to (7-A37) describes A.
Substituting (7-A30), (7-A37) and the analogous expression for 6 into
(7-A32) we obtain
n+1 ^ ^ n _h*At
FT, =2- Z A. e b At X v (7-A40)
10 n n_n UP UP
where
oo
1 a (n+ ^x )2
[- -s2- ] dp (7-A41)
oo
00
1 8
V = - P
2cr
-00
These last two integrals can be evaluated analytically and we get
= 1
304
A A
=
(7-A42)
x (7-A43)
= + 2
108
-------
Similar expressions give v , 0=1,,..6, except £ is replaced by e .
In summary, the value of the dependent variable f governed by (7-A1) is
obtained at time level n+1 at grid point (I,J) from (7-A40) where A"B is
derived from (7-A31) using known values of r at time level n; and where the \
and v parameters are given by (7-A43).
109
-------
Appendix B to Section 7
Here we describe the scheme we have developed to provide the boundary
values that the numerical analogues of the differential equations (7-86) -
(7-88) require but which the differential equations themselves do not need.
The essence of our method is that we divide the u, v, and h fields into a
base state and a perturbation component, and then in the governing equations
we set all sources of the perturbation fields to zero outside the region in
which the prognostic equations are applied. In addition we assume that
perturbations do not exist outside inflow points of the modeling domain. At
outflow points, we extrapolate interior values to estimate u, v, and h outside
the boundary, but these estimates are only included in the advection and
diffusion terms of the equation and are excluded from the forcing terms.
Rewriting the governing equations (7-86, 87, and 88) for future reference
and manipulation, we have
(7-Bl)
(7-B2)
110
-------
(7'B3)
where
P = 1 §E p = 1 §E a'
^x pQ3x ' Fy pQ3y ' g
Let
u = u + u1 (7-B5a)
v = v + v1 (7-B5b)
h = h + h1 (7-B5c)
where barred variables refer to the base state and primed variables denote
fluctuations from this state. Our aim is to absorb the effects of the given
large scale forcing terms P and P as well as initial and boundary conditions
x y
in the base state variables (u, v, h), and to let (u' , v1, h1) represent the
perturbations from this state that arise as a result of forcing by the terrain
z. and cooling n. within the simulated domain. Substituting (7-B5a) into
(7-B1) we get
* * '• + '• (7"B6)
-P + fv + fV- + KF^ + |!H +1!§' +1^']
" ^uavZ 3v^ &x. 3v
Now let u and v be solutions of
= -Pv + fv (7-B7a)
x
111
-------
= -P - fu. (7-B7b)
h y
and assume
8u _ 9v ~ n ~ J"a2u . a2u~| (7-B8)
at ~ at " u " "
Subtract (7-B7a) from (7-B6) and use (7-B8):
-v'-f) (7-B9)
a2u'
From (7-B2) and (7-B5b) we have
§*'*=£
. a2v'
T i) I
Z J
112
(7-B10)
-------
Using (7-B7b) and (7-B8) we can reduce (7-B10) to
-g.fh' - u.(fv+f)
y 8y V3x '
Now combine (7-B5c) and (7-B3):
,-, D10x
(7"B12)
where
_ , r
(7-B13)
and A is the model domain. We define
(7-B14)
113
-------
Subtracting this equation from (7-B12) we have
where we have assumed
* ' 0.
(7-B16)
Now we write the basic equations (7-B9), (7-B11) and (7-B15) in the following forms:
(7-B17)
= + S
(7-818)
where
= 3u/8x
= 8v/8y
= 8u/8x + 3v/8y
(7-B19)
(7-B20)
(7-B21)
(7-B22)
114
-------
C
C2
" i
- r(u2+v2)'5]) if u1 t 0;
, otherwise
(7-B23)
] , if v1 t 0;
-(u2+v2P , otherwise.
8x 8y
(7-B24)
(7-B25)
s
(7-B26)
(7-B27)
sh =
9v
8y
(7-B28)
(7-B29)
(7-B30)
(7-B31)
and h satisfies (7-B14); u and v satisfy (7-B7a, B7b) and (7-B8); and u, v and
h are defined by (7-B5).
115
-------
We require only first-order derivatives in evaluation of the p's and S's.
We use the notation
x=!Ax
y=JAy
A r „ \ - 9g
V9l-j) ' 7y
x=lAx
/=JAy
We define several different A and A operators, each of different orders of
x y
accuracy.
(7-B33)
= (9
I+lfJ
(7-B34)
+ ^I,J+3 ' fli,j-3)]/(606y)
(7-B35)
= C8(fl
I,J+l
(7-B37)
where 6 and 6 are the grid mesh dimensions.
x y
The following operators are used to compute first derivatives on
boundaries:
116
-------
1059i+2,J ' 1609l+3,J
) (7-B38)
" 16°9l,J+3
EAx(gI}J) =
) (7-B40)
^j) = [171fllfJ - j >
•659l,J-4 + 99l,J-5 - 9l,J-6^120V (7"B41)
The operators defined above are used to compute the variables jj , Pv, ... S'
defined by Eqs. (7-B20) - (7-B31) in the sequence illustrated by the flow
chart in Figure 7B-1. Following are descriptions of the operations indicated
in the flow chart.
Although the model domain is a grid of 60 x 42 cells, we use an array of
66 x 48 cells in solving the governing equations (7B-1, 2 and 3) in FLOMOD.
The extra cells comprise a "frame" 3 cells wide around the 60 x 42 modeling
domain. Values of all parameters within this frame of cells are specified
whereas those within the modeling domain are predicted using the governing
equations. Specification of the BAR variables Pu, J3V, Ph, $u, S"v and $h
within the frame is straightforward since the variables u, v, h, etc. on which
they depend are given at all points in the 66 x 48 cell region. Specification
of the PRIME variables p^, p^, p^, S^, S^, and S^ within the frame region is
guided by the desire to avoid the spurious generation of disturbances just
outside the modeling region that can subsequently enter the area of interest.
117
-------
The prime variables u1, v1, and h1 represent perturbations from the "base"
state represented by (u, v, h). We assume that all sources of perturbation
energy, such as terrain z.t cooling n'» etc. are within the 60 x 42 model
region. Therefore, we assign zero value to all prime variables outside this
region, namely in the boundary frame area. This is a key feature of our
boundary scheme. Finally, specification of u1, v1, and h1 in the frame zone
must be done arbitrarily. The prognostic equations (7-B1, B2, B3) cannot be
applied in the frame region because that would require values of all
parameters outside the 66 x 48 domain. A common method of estimating the
values of the dependent variables outside the modeling area is to extrapolate
values from the interior of the simulation area. We have found that even
crude extrapolation techniques give acceptable results when used with the
advection-diffusion equation. But the same methods generally fail when
utilized with systems of equations like (7-B1, etc.) because the errors in the
extrapolated values outside the simulation region give rise to disturbances
that spoil the accuracy of the solutions obtained within the model region. We
attempt to alleviate this problem by setting all source terms that involve the
dependent variables equal to zero outside the model domain (the 60 x 42
region) but we extrapolate values of these variables for use in the advection
and diffusion terms of the equation. In particular, if a given point on the
edge of the 60 x 42 region is a point of inflow, then we assume u' = v1 = h1 =
0 at all 3 cells of the frame zone adjacent to this point. For example,
u'(l,J) = u'(2,J) = u'(3,J) = 01
v'U.J) = v'(2,J) = v'(3,J) = 0 if [u(4,J) f u'(4,J)]> 0
h'(l,J) = h'(2,J) = h'(3,J) = OJ
At points of outflow, we use the following extrapolation, illustrated for the
case of point (4,J):
118
-------
READ FIXED
FIELDS
INITIALIZE
u'=v'=h'=0
READ
u. v, h, CD. g'. n. n'
COMPUTE
BAR VARIABLES
(0
c
o
+•«
CO
O
O
O
COMPUTE
PRIME VARIABLES
/8'u. P'v. (8'h,
S'u. S'v. S'h.
CALL SOLVER
COMPUTE NEW
u',v", h'
N
CD
SI
U
1
u.
00
r-
£
O3
119
-------
u'(l.J) = u'(2,J) = u'(3,J)
v'(l,J) = v'(2,J) = v'(3,J)
h'(l,J) = h'(2,J) = h'(3,J)
if [u(4,J) + u'(4,J)]< 0
where
= - 2u'(5,J)
(5/2)u'(4,J)
etc.
Calculation of the BAR variables.
Figure 7B-2 shows the grid on which the BAR variables J3 , S , etc. are to
be evaluated and the derivative operators that are to be employed at each
point.
From (7-B20) - (7-B22)
Pud,J)=Ax(uI>J)
Ph(I,J) = Pu
(7-B42)
(7-B43)
(7-B44)
and from (7-B26) - (7-B28)
120
-------
(7-B45)
(7-B46)
"I,JLVUI,J' V I,JJJ 'I,J (7-B47)
In evaluating the six expressions above, the operators A and A should be
x y
selected as follows (see Fig. 7B-2):
A • if I=1«
EAx , if 1=66;
2AX , if 1=2 or 1=65; (7-B48)
4AX , if 1=3 or 1=64;
6AV , if 1=4-63.
NAy , if J=48;
2A , if J=2 or J=47
4A , if J=3 or J=46
6A, , if J=4-45
where 2AY, 2A , etc. are defined by (7-B32) - (7-B41).
x y
(7-B49)
121
-------
:rrrcz:::::.:
46
45
J
f
1 1
1
1
1 I
6-
5.
i
ll
2^
f
'
j i
•^ —
\x .
! :•: :•:•
\ l;:i
> • «<• •;
| is:
1 f<: •
' * "::'•
. . I
\ \ t
^;
?
^;
;!|i
. .c/\
2
,-.-'J
6^y
1 = 1
i. >
- •™p*^.(7f.f!VTf^*'?!'!**.T**T.T'*"~
:':-x:::•:•:• ':•:x:: ::::.i
>;:i
i
fe& u ilitfu'-.'^'fe'-l ™*.A»—
i
\ \ \
1
60 61 62 63 64 65 66
Figure 7B-2.
Grid network on which pu, py, ph, §u, Sy, and §h
are computed. Different spatial derivative operators
AX and A are required in these calculations as indicated.
Calculation of the PRIME variables
The PRIME variables p^, p^,...S^ are defined on the same 66 x 48 grid
system as the BAR variables but the calculation procedure is different. As we
noted earlier all PRIME variables are set to zero outside the 60 x 42 cell
simulation area. Thus, we assume <
122
-------
p^(I,J) = P^CI.J) = P/,(I,J) = S^(I,J) = S^(I,J) = S^(I,J) = 0, (7-B50)
if I = 1,2,3,64,65, or 66; or J=l,2,3,46,47, or 48.
At all other grid points, namely 1=4-63, J=4-45 compute the PRIME variables as
follows:
0 , if u'Ij0 = 0;
(7-B51)
C2 u i u
UT i nT -i *,J 1 ,J r 1 ,J 1»*J
I,J I,J I J
if ui , ± 0
where
where
UI,J = UI,J + UI,J '
0, if
p;d,J) = J (7-B52)
[ (uf j + vf , - (uf , + vf
i h i,J i , J - i,J 1,0
VI,J hI,J hI?j
(7-B53)
(7-B54)
123
-------
sh(I>J) = ui,A + vi,jV"hi,J) + Ri,jPh(I'J) (7'B56)
In equations (7-B53,. . ,7-B56) the derivative operators A and A should be
selected as follows when they operate on PRIME variables:
EAX if 1=63;
Ax(4JsJ) =< 2AX if 1=5 or 1=62; (7-B57)
4AX if 1=6 or 1=61;
( _ 6AX if 1=7-60.
NAy if J=
/ " 2Ay if J=5 or J=44; (7-B58)
4A if J=6 or J=43;
if J=7-42.
In Eqs. (7-B53,.. ,7-B56) use 6AX and 6A in a11 applications to the BAR
variables u, v, and h. In this regard it should be repeated that Eqs.
(7-B51,..,7-B56) are applied only to columns 1=4-63 and to rows J=4-45, and at
these points the 6-th order operators can be applied to BAR variables (see
7-B48 and 7-B49). At all other points the PRIME variables are set to zero
(see 7-B50).
Calculation of the dependent variables u', v;. and h'.
The three dependent variables are computed at each time step using the
prognostic equations (7-B17,.. ,7-B19) within the 60 x 42 cell' modeling region
bounded by columns 1=4 and 63 and by rows J=4 and 45. These calculations are
done in the following steps.
124
-------
First we define the function
A(I,J) = SOLVER(I,J,IST,JST,£6x6,K) (7-B59)
which represents the biquintic algorithm described in Appendix 7A that is used
to solve the differential equations numerically. In Eq. (7-B59), 1ST and JST are
given functions of I,J,u' and v1 (see 7-A12, A13); £gxc is an array of 36
variables at grid points surrounding (1ST,JST) which we specify below (see
also Fig. 7-A1); and K is the diffusivity. The function SOLVER can be used to
predict the next values of u1, v1, and h1 at each of the 60 x 42 points
defined above. Consider each of the 3 variables in turn.
1. u'(I,J,N+l) in 60 x 42 model domain.
Step 1.
B*(I,J) = SOLVER(I,J,IST,JST,46x6,0)
where
^6x6 = PU + Pu
Step 2.
C*(I,J) = SOLVER(I,J,IST,JST,46x6,K)
where
= u'(N) +
Step 3.
= C*(I,J)*exp(-B*(I,J))
2. v'(!,J,N+l) in 60 x 42 model domain.
Same steps as in u1 calculation except replace Pu, p^, u'(N), §u and
Su by PV' Pv' V'(N)> ^v and Sv'
125
-------
3. h'(I,J,N+l) in 60 x 42 model domain.
Same steps as u1, except replace pu, 3^, u'(N), §u, and S^ by ph,
P^' n'(N)» $n and S^, respectively. Replace K by K^. In order to
avoid the problems that zero or negative fluid depth predictions
would cause, perform the following operation after each time step:
IF((h'(N+l) + h(N+l)).LT.1.0) h1 (N+l) = 1.0-h(N+l)
As we noted earlier, values of the dependent variables in the boundary
frame region are assigned zero values at all inflow points and are predicted
by simple extrapolation of the interior values at points of outflow. Consider
first the western boundary zone.
4. u'(N+l), v'(N+l) and h'(N+l) in western boundary zone 1=1-3, J=4-45.
If (u(4,J,N) + u'(4,J,N)) > 0, then
u'(l,J,N+l) = ul(2,0,N+l) = u'(3,J,N+l) = 0
Same holds for v1 and h1.
If (u(4,J,N) + u'(4,J,N)) <0), then
u'(l,J,N+l) = u'(2,J,N+l) = u'(3,J,N+l)
= -2u'(5,J,N) + (l/2)u'(6,J,N) + (5/2)u'(4,J,N)
Similar expressions are used for v1 and h1.
5. u'(N+l), v'(N+l) and h'(N+l) in the eastern boundary zone 1=64-66,
J = 4-45.
If (u(63,J,N) + u'(63,J,N)) < 0, then
u'(64,J,N+l) = u'(65,J,N+l) = u'(66,J,N+l) = 0
126
-------
Same holds for v1 and h1.
If (u(63,J,N) + u'(63,J,N)) < 0, then
u'(66,J,N+l) = u'(65,J,N+l) = u'(64,J,N+l)
= (l/2)u'(61,J,N) - 2u'(62,J,N) + (5/2)u'(63,J,N)
Similar expressions hold for v1 and h'.
6. u'(N+l), v'(N+l) and h'(N+l) in the southern zone 1=4-63, J=l-3.
If(v(I,4,N) + v'(I,4,N))> 0, then
u'(I,l,N+l) = u'(I,2,N+l) = u'(I,3,N+l) = 0.
Same holds for v1 and h1.
If (v(I,4,N) + v'(I,4,N)) < 0, then
-2u'(I,5,N)
Similar expressions hold for v1 and h1.
7. u'(N+l), v'(N+l) and h'(N+l) in the northern boundary zone
I = 4-63, J = 46-48.
If (v(I,45,N) + v'(I,45,N)) < 0, then
u'(I,46,N+l) = u'(I,47,N+l) = u'(I,48,N+l) = 0
Same holds for v1 and h1.
127
-------
If (v(I,45,N) + v'(I,45,N)) > 0, then
'(I,43,N) - 2u'(I,44,N) + (5/2)u'(I,45,N)
Similar expressions hold for v1 and h1.
The specifications of u1, v1 and h1 given in steps 4-7 above give the
dependent variables at time step N+l at all boundary zone areas except the 4
corner zones. The southwest corner zone is illustrated in Figure 7B-3. In
each of the four corner zones we will assign the dependent variables the
average of the values computed on the edges of these zones, as illustrated in
the Figure. For example, referring to Figure 7B-3 and keeping in mind that
u'(4,l,N) = u'(4,2,N) = u'(4,3,N),
and similarly for v1 and h1 ; and that
u'(l,4,N) = u'(2,4,N) = u'(3,4,N),
and similary for v1 and h', we assume in the southwest corner zone that
u'(I,JtN+l) = l/2[u'(4,l,N+l) + u'(l,4,N-KL)]
1=1,2,3 and J=l,2,3.
and similarly for v1 and h'. We use a similar method to compute the dependent
variables in the northwest, northeast, and southeast corner zones.
128
-------
1
!•
i
f
6*
i
I
5|
i
"WEST*
BOUNDARY
• 1 •
ZONE
'I'
.*.
lllltiilillP^
||;CpRiN!ER^|
io'ftBinDniiUTci
t: • :•-:•: i'xx^x-x-x'x-x::- .
I: : •••,:': ::^^x.::x:;xv:'x:x:::x:.Vx-:
j. x : •.•: ' :•: :•:-.•:•: :•:-.-. •••:•:•••:•:•:•:•:•:•:•:• •.•'
i;: xAMobEL ooMAiN ^•••^
1
1
I-* SOUTH BOUNDARY Z(
i
1
i
J = 1
1=1 234
Figure 7B-3. Illustration of the southwest corner zone and the
west and south boundary zones of the model domain.
129
-------
SYYYYYYYYYYYY
V V V V V V V
130
-------
SECTION 7
PROCESSOR P8
INTRODUCTION
This processor determines the top surfaces H2 and HS of layers 2 and 3 of
the regional model, it computes mean vertical velocities on these surfaces,
and it estimates convective cloud updraft speeds and other parameters required
in the specification of pollutant fluxes across surfaces H2 and HS. All of
these quantities play important roles in the regional model, but unfortunately
none of them is directly measurable. In this section we outline a procedure
for deriving estimates of these parameters that are consistent both with
observational data and physical principles.
DERIVATION OF BASIC EQUATIONS
Recall from Part 1 that
H2(x,y,z,t) = Z2(x,y,t) - z (8-1)
where z2 is nominally the elevation of the top of the mixed layer. During
clear daylight hours, z2 is the highest elevation that nonbuoyant pollutants
can reach. When convective clouds are present, we take z2 to be the elevation
of cloud b?ses, specifically, the so-called lifting condensation level (LCL).
In this case pollutants entering the updrafts that feed individual cumulus
clouds can rise above elevation z2 and proceed as far as the elevation of the
cloud tops, which we define to be the elevation zs of the top of Layer 3.
131
-------
When convective clouds are absent, z3 is defined to be z2 + h3, where h3 is
some constant depth of order 100 meters. Figure 8-1 illustrates the surfaces
H2 and H3 in a regional domain in the typical condition where cumulus clouds
are present over only a portion of the modeling region.
We showed in Part 1, Eq. (4-20) that z2 satisfies the differential
equation
(8-2)
Here w is the entrainment velocity (represented in Part 1 by f6/A6); w2 is
the cell averaged vertical air speed at elevation z2; w is the effective,
cell averaged cumulus updraft velocity at cloud base (i.e., at z2); V2u is the
horizontal wind velocity at elevation z2; and QC is the fractional area
covered by convective (cumulus) clouds. The last parameter is measurable from
Figure 8-1. Illustration of surfaces H2 and H3 during situations in which
convective clouds cover only a portion of the modeling region.
132
-------
satellite photographs and it is the only variable in (8-2) that can be
estimated reliably. We will assume that ac(x,y,t) is an input to this
processor.
All the other variables in (8-2) must be inferred through comparisons of
this equation with certain measured data and other physical principles. This
is the task we undertake in this section.
An auxiliary relationship that will aid in the estimation of w is the
water vapor mass conservation equation. Using Eq. 2-29 of Part 1 we can
express the layer averaged water vapor mixing ratio q in the form
a 31 nV. .
^ j + j -inr1 + vV^j + V FMJ ' FJ,J] = ° (8"3)
j
where . is the cell averaged mixing ratio in layer j, A is the horizontal
J
area of a cell, V. is the volume of a cell lying in layer j and F. ,, is the
J J > *
flux of material across surface H. to or from layer k. To arrive at (8-3), we
J
dropped the horizontal flux terms in (2-29) of Part 1 that represent the
effects of subgrid scale variations in the horizontal wind v.u. Our interest
~jn
is in 3 during periods when cumulus clouds are present (i.e., 0 2 0). If
we assume that the advection term in (8-3) is negligible in these conditions
compared to the other terms, we obtain with the aid of Eqs. (4-3) and (4-21b)
of Part I
3 + i3-i2] + _!_ [CT (i^±£
H 3 z3-z2 3 2 z3-z2 L cv l-ac
we)(qc-3) +
-3Zs] = 0
133
-------
where z2 = 3z2/8t, Z3 = 8z3/3t and
H3 = z3 + vs'Yz3 - w3. (8-5)
Also, qc represents the mixing ratio in cumulus clouds and q^ is the mixing
ratio just above the elevation of cloud tops, i.e., just above z3. If we
assume that the horizontal advection term in (8-5) is negligible, we can
reduce (8-4) to
at 3 + zz Cac + we)(qc-3)
(8-6)
and after multiplying by (z3-z2) and rearranging terms we obtain
[w2 * (l-ac)we][3 + Js!J£j
c (8-7)
+ w3(qoo - 3)
where
M3(x,t) = q(x,zft)dz (8-8)
z2(x,t)
and
-? (8-9)
At each of the sites x of rawin stations we have available measurements
~ro
of q(z), v2u(z) and 6(z) at each observation time. We also know a and z3
(where a ^ 0) everywhere from satellite data. Using all this information we
can derive estimates of all the parameters in (8-7) except WG and wg, and
134
-------
thereby we can derive a relationship that w and w must satisfy at each
rawin site x . This relationship can subsequently be employed in conjunction
with (8-2) to obtain estimates of z2, w and w throughout the modeling
domain.
In order to use the rawin data and Eq. 8-7 in this way, it is convenient
to integrate (8-7) between each set of measurement times. Let t0 and ^
denote two successive times when rawin observations are gathered at a given
station. Ordinarily, tt and t0 will be 12 hours apart, but shorter intervals
are possible. Integrating (8-7) we get
Cp~ wc - we(l-oc)3 - weqcoc]dt = M3(ti) - M3(t0)
to
(8-10)
x q a
[w2(3 + i~ ) - Va + w3(qoo-3)]dt
where all variables are evaluated at a given rawin site x . We assume that
/V|||
only w and w are unknowns in this equation. For all other parameters we can
use the following expressions.
First, we adopt the approximation
qc ~ q(z2). (8-11)
Then
qc(t) ~ q(z2,t0) + qc(t-t0) (8-12)
135
-------
where
and
qc = Cq(z2,ti) - q(z2,to)](t1-t0)"1 (8-13)
q(z2,t0) = q(xnj>z=z2(to))to).
Next, we approximate z2 and z3 by
z2(tn) = FCqCx^z.t^, e(xm,z,tn)) n=0,l (8-14)
where F is a function that is defined later; and
fz2(t ) + 100 meters, if a (x t )=0;
z3(tj = c m n n=0,l (8-15)
(_ from satellite data, otherwise
z3(t) = [zs(ti) - Z3(t0)] (ti-to)"1. (8-16)
Since (8-4) bears the implicit assumption that fls > 0, which we adopted
in anticipation of application to situations where convective clouds are
present and growing, we must ensure that the value of z3 given by (8-16)
satisfies
z3 > 0. (8-17)
This constraint is most likely to be violated when a = 0 either at
observation time to or at ti, but not both.
In this instance we can adjust the estimate of z2 at the hour that QC = 0
to achieve z3 > 0. Adjusting z2 is consistent with the position adopted later
that only a range of z2 values can be estimated with confidence from the
measured data at any rawin site and time.
136
-------
We approximate the temporal variations in qw by
(8-18)
and from this we assume in analogy with the form we adopted for q , i.e., Eq.
8-12,
qjt) = q(z3,t0) + qoo(t-t0) (8-19a)
where
cia, = CqUs.ti) - q(z3,t0)] (tx-to)" . (8-19b)
From (8-8) we get
M3(xm,t1) = 1 q(xm,z,t1)d2 (8-20)
and
Z3(xm,t0)
q(xn),z,t0)dz (8-21)
Similarly,
3 = 3 + 3(t-t0) (8-22)
where
3 = [3 - 33/(ti-t0) (8-23)
and
The vertical air speeds w2 and w3 that enter into (8-10) will be
approximated by
137
-------
*2(xm,t) = w(xB,z2(t),t) , to < t < ti; (8-25)
ws(xm,t) = w(xm,z3(t),t) , to < t < ti; (8-26)
where w(x,z,t) is the vertical air speed at site x, elevation z(MSL) at hour
t. There are various methods of determining w from the set of rawin data but
we will not consider any of them here. We will simply assume that some means
of estimating w exists and we will leave it to the model user to choose a
feasible scheme. [In implementating the present regional model, we used the
method developed by Bullock (1983)]. The elevations Za(t) and z3(t) at which
w(t) is evaluated in (8-25) and (8-26) are taken to be
z2(t) = Z2(xm,t0) + Z2(t-t0) t0 < t < ti (8-27)
where
Z2 = [Z2(xm,ti) -Z2(xm,to)]/(ti-t0);
z3(t) = Z3(xra>t0) + Z3(t-t0) (8-28)
with z3 given by (8-16).
Now that we have formulated approximate expressions for each of the
parameters in (8-10) except WG and we, we can define a constant wc which
satisfies
ti
T^T dt = \ ^^3 + dr(qc-s)]wpdt
J. U_ J I. v. c
to to
A(xm; to.
(8-29)
138
-------
where WG = wc(xm; t0sti) and
r'1
- M3(t0) + [w2(3 (8-30)
to
W3(qco-3)]dt.
The parameter w is the effective cumulus cloud updraft speed w character!' s-
L* C*
tic of the air mass in the vicinity of the rawin site xm during the interval
to < t < ti between observations. In other words, in the vicinity of x we
assume
wc(x>t) = w^; to.ti) to < t < ti, and
i i (8~31)
x-xm < 6
I ~ ~m I
where wc satisfies (8-29). Note that (8-29) relates w to the unknown
entrainment speed w , which we consider next.
Continuing our examination of the region around x , let us look at the
air parcel trajectory that ends at xm at time ti. Figure 8-2 depicts it
beginning from its position at time to. Writing Eq. (8-2) in the form
dz w a
and integrating it from to to ti along the trajectory ending at xm at time
we get
139
-------
dt e,ti), or more precisely a range of values in which we expect z2 to
lie, from (8-14). But x' is generally not the location of meteorological
measurements so z2(x(J),to) can only be estimated from interpolation procedures.
Let us suppose that through a combination of measurements and interpolation we
are able to say with confidence that
(8-34)
Z2(xffl,ti) - z2(x^,t0) = Az2(xm;t1,t0) ±
6z2(xm;ti,t0)
where Az2 and 6z2 are known, partly from (8-14). We can now write (8-33) in
the form
ti tj ti
r w . r o- c
Az2 ± 6z2 =(t) J-g-dt - w (D TT£- dt +(K w dt. (8-35)
J uc J c T
to to to
The cloud cover fraction a (x,t) is known everywhere and we assume that w2 is
C* ***
also available at all points and hours. Thus, (8-35) contains only 2
unknowns: WG and wg.
140
-------
Figure 8-2. Illustration of the air parcel trajectory
that arrives at point x at hour ti. Point
x' denotes the parcel's location at time
?m <- t
to < ti.
Following a number of previous investigators (see, for example, the
review by Artoz and Andre, 1980) we shall assume that the kinematic heat flux
at z2 is proportional to that at the ground, namely
= = aQ
(8-36)
where Q is the known surface heat flux and a is a constant. Recall that
we =
A8
(8-37)
where A8 is the effective jump in potential temperature across the elevation
Z2- By comparing simple mixed-layer growth rate formulas with measurements,
141
-------
Artoz and Andre (1980) found that a simple, reasonably accurate expression is
(8.38)
where Y = d6/dz evaluated at z=Z2+£ and h=Z2-zf, where zt is the local terrain
\f U
elevation. Eq. (8-38) suggests that within the vicinity of rawin station x
we can approximate w by
wp(x,t) = G(xm;t0,t1)Q(x,t) , t0 < t < ti; (8-39)
x-xm < 6
~ ~m
where G is an unknown function that we suspect from (8-38) is of the order
G(xm;t0jt1)"1 ~ [3! e(xm,z=z2+e,T)][z2(xm,T) -
(8-40)
where t is some instant in the interval to to
Substituting (8-39) into (8-35) we obtain
W 10 L
AZ2 ± 6Z2 =Ci)^- dt " ^r Q-FF- dt + G(P
c J c J
t(j to to
where the integrals are all along the trajectory shown in Fig. 8-2. Making
use of (8-39) in (8-29) and substituting the result into (8-41) we get
142
-------
tl tt
Az2 ± 6z2 = G> ^- dt - i- [G j [3 + ac(qc-3)]
to t0
(8-42)
Q(xm,t)dt + A] + G© Qdt
~m '
Jt0
where
5 = kll£ . (8-43)
qcac
dt
Fa dt
to °
We repeat that in our notation^ £dt denotes integration of £ evaluated at the
& 11 i
f" " i. *• o
fixed point x while© ?dt denotes integration of t(£,t) along the space-time
» tg
coordinates illustrated in Figure 8-2.
In expression (8-42), G is the only unknown parameter. For each number
AZ2 in the interval indicated on the left-hand side of (8-42), i.e.,
Az2 - 622 < AZ2 < AZ2 + 6z2 (8-44)
there corresponds a G given by
[AZ2 -0- dt + A]{(1> Qdt - - [3
^
' 4- ^* *C * 4- ^*
to ^0
(8-45)
Q(xm.t)dt}-1
143
-------
Let us assume that the set of G associated through (8-45) with the set of AZ2
defined by (8-44) lies on the interval
Gmin < G i G«x <8-46>
Only positive values of G are physically meaningful because both w and Q are
positive during convective conditions. Therefore, if the interval1 defined by
(8-46) contains only negative values, one or more of the parameters entered in
(8-45) have erroneous values. When this situation arises in practice, we
propose to alter the values of the various terms in (8-45) until the resulting
G interval at least contains positive values. At the most we would want the
largest positive values of G to be consistent with (8-40). Alterations of the
parameters in (8-45) should proceed according to a fixed rule in which the
parameter suspected of being the most inaccurate is altered first, and
succeeding terms in the hierarchy are modified only after adjustments to the
least accurate terms have failed to produce the desired values of G. The
alterations made in the value of any one parameter should be confined to the
smallest interval in which one could reasonably assume that the correct value
lies. Regarding the order in which the parameters in (8-45) should be
altered, we propose the following hierarchy based solely on intuition:
qc, qc, 3, Q, w2, AZ2, CTC. (8-47)
Suppose that through a process like that described above we have managed
to extract from (8-45) ah interval of G that contains some positive values.
Our interest is only in the positive values, i.e., the G in the interval
i G i Gmax (8'48)
144
-------
To each G in this interval there corresponds a WG through relationship
(8-29), namely
gc (qc "
(8-49)
Q(xm,t)dt}
Let the interval defined by (8-48) and (8-49) be represented by
A A
w - 6w < w < w + 6w . (8-50)
Like G, WG is intrinsically a positive quantity. Therefore, if (8-50) does
not lie at least partially on the positive, real axis, we must perform
alterations on parameter values in (8-49) until the interval (8-50) contains
some positive values. The procedure should be similar to that used to obtain
positive G values, except in the case of w we should alter only those
I*
parameters in (8-49) that do not appear in the expression for G. Otherwise,
we would cause changes in the G values as well. Thus, the suggested hierarchy
of parameters that should be modified, if necessary, to achieve positive w
values is the following:
q^, M3, w3, Z3. (8-51)
At the conclusion of this operation we obtain a positive set of w values
and a positive set of G values, both of which apply only to the time interval
to < t < ti and within the vicinity of rawin station m. The next step is to
apply the same procedures to the rawin observations made at station m+1 at
hours to and ti to obtain the corresponding sets of WG and G values for that
145
-------
site. Once values have been obtained for all rawin stations for the interval
to < t < ti, we can interpolate values for w and G at all grid points in the
model domain. With these fields and the known wind fields v2u(x,t), w2(x,t),
and the cloud cover distribution a , we can solve (8-2) for z2 at each grid
point and each hour t in the interval t0 < t < ti. By repeating the entire
process for the next observation interval ti < t < t2, we can obtain z2 and
all the other necessary fields during this period.
Before outlining the specific sequence of steps necessary to implement
the procedure described above, let us comment briefly on the philosophy of
this approach. Many prognostic models have been developed in recent years for
predicting the mixed-layer depth z2 given the surface heat flux Q, the mean
vertical velocity w2 and other physical quantities. These models are not well
suited to our needs for several reasons.
First, nearly all these models are one-dimensional and therefore they do
not take into account advection and horizontal variations in z2, w2, Q and the
other governing fields.
Second, our interest in regional modeling is with historical situations
where available observations exist from which Z2 can be inferred, at least
approximately, not only at the initial moment t0 of the simulation period but
also at discrete intervals throughout it. In general, the z2 predictions of a
model initialized at to will become increasingly inconsistent with later
observations due to deficiencies in the model and errors in the input data.
Our position is that the values of z2 inferred from meteorological observa-
tions made during the simulation period are more credible than the predictions
146
-------
made by a model. Therefore, in our approach the physical principles on which
models are based are employed in the role of interpolating and extrapolating
the discrete observations. This is the essence of the roles performed by Eqs.
(8-2), (8-45) and (8-49) in our scheme above.
At this time, the proposed procedure has not been tested. Therefore, it
should be viewed as the starting point in the development of a scheme capable
of providing the various required parameter fields.
In the remainder of this section we present the detailed steps needed to
produce an operational processor.
Stage ZQ
As noted in the introduction, when cumulus clouds are not present during
daylight hours, z2 is the highest elevation that dry thermals produced by
surface heating can reach; however, once cumulus clouds form, z2 is defined to
be the lifting condensation level. In either case, z2 can be inferred, at
least approximately, from radiosonde data and we expressed this in the form
(8-14) of a function F that relates z2 to the potential temperature and mixing
ratio vertical profiles. Our first task here is to develop an explicit form
for F.
i
Our approach stems from the realization that turbulence acts to destroy
spatial variations in scalar quantities. According to the empirical K-theory
147
-------
description of turbulent mixing, the rate at which spatial fluctuations are
eliminated is inversely proportional to the square of the size A. of the
fluctuation. This is evident from the Fourier transform of the classical
diffusion equation
H =
3z
From this equation one finds that the Fourier amplitude A of spatial variation
in the scalar c of wave number k = 2n/\ decays in time at the rate
f£ = -Kk2A. (8-53)
Sources of c can generate small scale variations but turbulence always
destroys them.
These observations suggest that if second and higher order derivatives of
a conservative, scalar quantity remain large for an extended period of time at
points in the fluid that are far from sources, the intensity of turbulent
mixing at these points, as manifested in the diffusivity K in Eqs. (8-52) and
(8-53), must be very small. Otherwise, the small-scale variations that cause
large values of these derivatives would be eradicated.
Thus, our basic premise is that in cloud free conditions we can estimate
the elevation at which the vertical diffusivity becomes vanishingly small by
examining the vertical profiles of the second, and perhaps higher, order
derivatives of the mixing ratio and potential temperature, both of which are
(approximately) conservative quantities that can be derived from the
radiosonde measurements.
148
-------
Consider for example the classical profiles of mixing ratio q and
potential temperature 6 in dry, convective conditions illustrated in Figure
8-3a. In Figure 8-3b we show the corresponding profiles of d2q/dz2 and
d28/dz2. In this idealized example, the product of these derivatives has a
large negative value at the top of the mixed layer, i.e.,
<« o at z = z2. (8-54)
dz2 dz2
We can estimate the derivatives of a given parameter £ measured at
discrete points in space by representing the parameter values measured in the
vicinity of the point of interest by a polynomial. Suppose that we want an
estimate of d2£/dz2 at the point z=z" but that £ is known only at discrete
points ZL z2, ... Zj that are not necessarily equally spaced. Let us denote
the £ values at these points by
&i =C(V zi'to) • i=1.2.---I (8'55)
where we assume that £ is a parameter measured at rawin station m at time to-
We can estimate the second derivative of t, as well as other properties of
interest at the desired elevation z" by expanding £ in a polynomial about z".
Thus, let
t(z) = a0 + ain + a2q2 + a3r}3 (8-56)
where
n = z-z". (8-57)
149
-------
d20
(a)
(b)
Figure 8-3. (a) Idealized profiles of mixing ratio
q and potential temperature 0 in dry,
convective conditions.
(b) Second derivatives of the profiles
illustrated in panel a.
We can obtain the four constants a0, ... a3 in (8-56) from four values of £
and their corresponding measurement locations z.. The most accurate
150
-------
representation is obtained by using the four successive measurements of £,
shown in Figure 8-4, that straddle the point z". Two of these are from
elevations below z" and two are from higher elevations. For convenience, let
us denote these four points by zl9 z2, z3 and z4, as indicated in the Figure;
and let the associated £ measurements be designated £x ... £4. In this
notation we have from (8-56)
Figure 8-4. Illustration of the points (dots) at which
measurements of £ are available; and the point
z" at which a measure of d*£/dz2 is desired.
Values of £ measured at the points z., i=l,...4
centered at z" are used to approximate £(z) in
a polynomial about z".
= a0 + atni + a2n.? + asnf
*
(8-58)
where
n. = z.-z". (8-59)
I i
Solving the system of equations (8-58) for the a's is straightforward and we
obtain
151
-------
a = BiAtg - 2An
2 A2A12-A11A2
2112-1122
(8-60)
(8-61)
(nf-ni)
(8-62)
(8-63)
(8-64)
(8-65)
(8-66)
(8-67)
(8-68)
(8-69)
The coefficients a0, ... a3 provide the following useful properties of £:
A22 =
- (nl-n§)(n?-ni)
- (n§-nl)(n|-nl)
t(2=z") - a0
(8-70)
(8-71)
2=2"
152
-------
= 2a2 (8-72)
(8'73)
z"+6z
C(z)dz = 2(6z)a0 + |(6z)3a2 (8-74)
z"-6z
Now that we have formulated a method of estimating the derivatives of
parameters measured at discrete points, we can proceed to estimate Z2.
Let
ezz = . (8-76)
The steps for Stage 2Q are as follows:
(1) Applying the procedure outlined above for estimating derivatives.
[specifically (8-60) - (8-69) and (8-72)] to the mixing ratio and
potential temperature soundings q (z,to) and 0 (z,to), respectively,
compute
k = 1, ... 60
at rawin station m (=1 on the first pass through this stage) at hour to, where
153
-------
zk= zt(v + 25m + ^ (8"77)
and Az = 50 m.
(2) Next form the products
Pk = ^zk''W . ""I. - W (8-78)
and from this find
P = max(-P.). (8-79)
k k
This value of P pertains to the site x of rawin station m and to hour to of
the sounding.
(3) Next we estimate z2 to be (units of m MSL).
= F(q(xm,t0),
(8-80)
= zt(xm) + 25 + 50k*
where k* is the smallest integer for which
Pk* = - P. (8-81)
(4) We also need for use in Eq. 8-34 an estimate of the range ±6z2 of
elevations in which the actual elevation z2 of the mixed layer top
is likely to lie.
154
-------
We anticipate that a measure of 622 is the width of the interval
centered at z=zk* in which Pk has a magnitude of, say, 60% of its
peak value . Examination of the third order derivatives of q and 0
might also provide a useable measure of 622- At this time, no
information is available on which to base a quantitative rule for
estimating the range of 22 values. Therefore, we will assume for
now that
622 = 100m (interim assumption) (8-82)
and after we have gained experience with actual soundings, we will
attempt to formulate an empirical rule for estimating 622.
(5) Compute the lifting condensation level (LCL) at each surface weather
station n=l,...N at the hour to of the upper air observations used
above in step 3 to estimate 22- The elevation of the LCL is found
first in pressure coordinates as follows.
Let q and 0 be the mixing ratio (dimensionless) and potential
temperature (°K), respectively, at surface weather station n.
These data are available from Processor P3. If a parcel of air
originally at ground level in the vicinity of the station is lifted
without mixing with ambient air, both the mixing ratio and potential
temperature will be conserved. Thus, at any altitude p(mb) the
vapor pressure in the air parcel will be
pq
(cf Eq. 3-2a), and its temperature will be
155
-------
(8-83b)
n p
(cf Eq. 3-6). The elevation p.p. of the lifting condensation level
is defined as the altitude where the parcel's vapor pressure e is
the saturation vapor pressure e . The latter is a function of
temperature alone, namely
es = 8oexf4(T ' T)] (8"83c)
where eQ = 40mb, TQ = 302°K, L = 2500 joule g"1, and R=0.461 joule
g °K . Hence, p.p. is the solution of the equation
0.622"+ q = eoexp[S (3U2 " T)] (8"83d)
where T is given by(8-83b). Eq. (8-83d) can be solved approximately
by substituting successively smaller values of p into the equation,
beginning with the surface pressure; and considering the solution
PLCL to be the Pressure at wnich tne left Slde of (8-83d) first
exceeds the right side.
Using the pressure-height functions PmC2'^) available from
processor PI, convert p.p. into elevation Zj^, (m MSL) as follows
zLCL(xn,t0) = z* (8-84)
where z* is the elevation for which
p(xn,z*,t0) = PLCL(Vt0). (8"85)
156
-------
The function on the left side of (8-85) can be obtained by applying
an inverse r weighting interpolation to the set of Pm(z,t0)
profiles, namely
M
P(Xn,z,to) =
At the end of this step we have values of Z|Q.(to) at all N surface
weather stations. The M rawin stations whose soundings we are using
to determine Z2 are a subset of the N surface stations.
(6) At this step we must decide whether z2 is to be the elevation given
by (8-80) or the elevation z,p. given by (8-84). As we noted
earlier, the decision rests solely on whether a (x ,to) is greater
than or equal to zero. In particular
Z2(xm,t0) if ac(xm,t0) = 0;
z2(xm,t0) =
zi_CL(~m'to) (Eq> 8~84)> otherwise.
where
ZjsCX.to) = value given by Eq. 8-80. (8-88)
There are two situations here that signal the presence of an error
in our estimates of z£ and/or Z. The first is
(xm,t0) = 0 and z^LCx^to) <[zKxm,t0) - 6z2]; (8-89)
157
-------
and the second is
ac(Vto) * ° 204 z|_CL(Vto)
The first condition indicates that clouds are not forming even
though our estimates of z^ and the mixed layer depth indicate
that they should; and (8-90) implies the opposite, namely, that
clouds are forming even though our estimates indicate that they
should not. If either (8-89) or (8-90) is true at site xm at hour t0
this should be recorded for output from Stage ZQ. Frequent
occurrences of these conditions would indicate some systematic error
in the calculation procedures.
(7) Next we compute 23. As we noted earlier (see Eq. 8-15)
Z2(xn,to) + 100 , if ac(xm,t0) = 0;
(8-91)
m'to) * Otnerw1se
where ZjClJ is the average elevation (m MSL) of the tops of cumulus
clouds, an input to P8 derived from satellite data.
(8) From (8-11) and (8-70)
qc(xm,t0) = q(xm,z2,t0) = a0>q(z2,t0) (8-92)
where a is the coefficient aQ of the expansion of q about the
point z"=z2. This coefficient is found from the mixing ratio
sounding and Eqs. (8-60) - (8-69).
158
-------
(9) From (8-18) and (8-70)
q00(xm,t0) = q(xm,z3,t0) = a0jq(z3,to) (8-93)
where a is obtained as in step (8-8) except the expansion is
about the elevation z"=z3.
(10) From (8-8) and (8-74) we have
M3(xm,t0) = /[Az a0)q(zk,to) + ^ (Az)3a2)q(zk,t0)] (8-94)
where k2 is the altitude interval given by (8-77) that is nearest
Z2(xm,t0), i.e.,
Iz2(xm,t0) - (zt(xm) + 25 + k2Az)| = minimum (8-95)
and similarly k3 is the integer that minimizes
I23(xm>t0) - (zt(xm) + 25 + k3Az)| = minimum; (8-96)
and Az= 50m as in (8-77). The coefficients a and a2 in (8-94)
are derived from the q sounding data using formulas (8-60) - (8-69).
(11) From (8-9)
- z2(xmJt0)
(8.g7)
(12) Repeat steps 1-8 at rawin station m for the next observation hour
to obtain Z2(xm,t1), Z3(xn),t1), qc(xm,t1), q^C
,t!) and 3. (In general, ti = t0 + Atm, where Atm is
159
-------
the interval between observations at station m; so we should
actually write tx rather than tj to designate the second
observation time. However, this distinction is not important
because in the analyses that follow, we treat the observations on a
station-by-station basis and interpolate hourly values of the
desired quantities at each site.)
(13) Construct linear functions of qc, q^ and 3 for use in integrating
these values with respect to time from t0 to ti. That is
t0 < t < ti (8-98)
3(t-t0)
where
% = tQc(4.t!) - qc(xm,t0)](trt0)"1 (8-99)
with similar expressions for q^ and 3 [see (8-19b) and (8-23)].
(14) Repeat steps 1-10 above for each of the M rawin stations, and then
repeat all steps for each observation interval in the period for
which the regional model is to be operated. The observation
interval is usually not the same at each upper air station. Figure
8-5 illustrates a hypothetical situation in which the model
simulation period is of length T beginning at hour t0 and the
stations make soundings at a variety of intervals.
160
-------
Station m
Station 3
Station 2
Station 1
t t t
om im 2m
Model Simulation Period
to t0+ T.
Figure 8-5. Illustration of possible relationships among
the intervals at which upper air soundings are
made at a set of rawin stations.
At the end of the operations in Stage ZQ there should be values of
Z2, 6z2, 23, qc> q^, 3 and M3 for each station m=l,...M and for
each observation interval in the model simulation period T.
illustrated in Figure 8-5.
Stage PATH
In Stage ZQ we derived estimates of the parameters required to evaluate
the fixed point time integrals that enter in Eqs. (8-45) and (8-49). In this
stage we estimate the trajectories necessary to evaluate the path integrals §
that enter into these equations.
161
-------
We are still considering the time interval to < t < tj. and we require the
backward trajectories that begin at each of the rawin stations at hour ti and
go back in time to to- That is, we need to know the location x(t;x jtO
during the interval to < t < tj of the air parcel that arrives at station x
at hour tj. By definition
(8-100)
and by previous declaration (see Figure 8-2)
i). (8-101)
We will compute the trajectories using the horizontal velocities v2u
measured in the vicinity of z2. These velocities control the rate of
advection of the mixed layer height (see Eq. 8-2); but the effective path
along which the surface heat flux Q is determined is more likely fixed by the
vertically averaged horizontal flow beneath z2. In actuality the temperature
of a vertical column of air of depth z2-zt at x at time tj is affected by the
surface heat flux in a plume-shaped area whose width a increases from zero at
xm at t=ti to a value of the order
~m x
c-y(x;) ~ (t1-t0)Av (8-102)
at the upstream starting point of the trajectory that ends at x at t^ In
(8-102) Av is the magnitude of the vertically integrated horizontal wind shear
in the mixed layer. Since we are representing the heating rate of the column
as an unknown function G multiplied by the surface heat flux variations along
the trajectory between x1 and x , we can assume that the difference in the
J J ~m ~m
integrated heat flux along the trajectory extracted from the winds v2^ at the
level z2 and that along the "correct" trajectory is absorbed in the function
G.
162
-------
To compute the trajectory x1 to x we will assume that there exists a
~m ~m
routine, or stage, which we will call WV, which returns an estimate of the
vertical wind speed w and the horizontal wind vector VM at any given
space-time point (x,z,t) in the NEROS domain. Thus, when we write w(x,z,t) or
vu(x,z,t) it will be understood that these values are available from the
~n ~
routine WV which we will not specify here. (In the first generation model we
will employ the scheme developed by Bullock (1983) in the role of stage WV.)
For later reference let us signify this in equation form:
w(x,z,t)
from Stage WV (8-103)
The trajectories can now be generated using the following recursive
formula:
-AtvH(x(t1-lAt;xm,t1),z2,t1-lAt)
where At is a time step of order 30 minutes. Execution of (8-104) is the
first step in the operation of Stage PATH, i.e.,
(1) Solve (8-104) for the trajectories x(t;x .tj), t0 < t < tj at each
of the M upper air stations.
In the next steps we use these trajectories to evaluate the path
integrals in (8-45) and (8-49).
163
-------
(2) Determine the variations in w2, a., and Q along each of the M
trajectories generated in step 1 for the period t0 < t < tj. and
express the results in the following functional forms:
to
W2m(t) = *<*(*: V*^' V^ <
acm(t) = (^(xCtjx^tO.t) (8-106)
Qm(t) = QCxCtjx^tO.t) (8-107)
where
*2m = ^Z2(xm,to)^2(xm,t1)]. (8-108)
(3) Now compute the following path integrals:
r » — *• i~n "" *"rrnx "•*• o*"*~' *• iU*/y
t0 c J u cm
(8-111)
to
where
J = (t!-t0)/At (8-112)
(4) Repeat the three steps above for each of the M rawin stations and
164
-------
for each of the observation intervals within the model simulation
period (see Figure 8-5). At the end of Stage PATH, there should be
values of the path integrals I , I and In (Eqs. 8-109 - 8-111) for
W U w
each rawin station for each observation interval in the simulation
period T.
Stage WEWC
Here we attempt to solve (8-45) for the function G(xm;t1,t0), which will
provide the entrainment velocity field w (x,t) through (8-39); and equation
(8-49) for wc(xm;ti,to) which will provide the cumulus updraft speed w (x,t).
(1) Using the estimates (8-98) of the temporal variations in q , and
3 at each upper air site x^ in the period t0 < t < ti, and the
surface heat flux Q(x_>t) in this period, compute the following
integral:
J
Ux_;ti,t0) = At2[3 + a_(x .ti-jAt)'
~m ^o m c ~®
(qc(xm,t1-jAt)-3)]Q(xm,t1-jAt) (8-113)
(2) Compute the parameter (see 8-43):
-i JL qr(xm,ti-jAt)a (x ,ti-jAt)
= [I.-CtO] AtI-£-S $-2
am J=o i-arx^trjAt)
where xcn)(ti) is from Stage PATH, Eq. 8-110.
165
-------
In Eq. 8-45, AZ2 is a measure of the change in z2 between the starting
point x^ and the end point x^ of the trajectory that arrives at x at time tj.
As we noted earlier, we have an estimate of z2 at xm at tt but no measurement
of z2 at x' at t0.
Let us assume that t0 is the initial instant of the time period in which
data for the regional model are required and that t0 is also an hour at which
rawin data are routinely collected. At this initial moment we must use some
objective analysis method to determine z2 in the model domain. Subsequently,
we can use the prognostic equation 8-2 to obtain z2. Thus, to estimate z2 at
x^ at time t0, we first estimate z2 at each of the m rawin stations using
(8-87) (Stage ZQ), and we estimate the likely error bounds ±5z2(xm,t0) on this
estimate at each x^ (see step 4 of Stage ZQ). Next, we apply the r"
interpolation formula to obtain Z2(x||),t0). To do this we must take into
account the possibility that cumulus clouds are present between rawin
stations. In this instance the z2 estimates obtained at the station sites x
will not represent the lifting condensation level, which is the altitude of z2
wherever ac(x,t0) 1 0. Thus, we assign Z2(x(|1,t0) values according to the
following rule:
M
I tr,.
I j r. I z2(x.,to)
i=l' ~1m ]
' 1f ac(*m'to) = 0;
(8-115a)
, otherwise
(8-115b)
166
-------
where
= E(x -x')2 + (y-yi)2]11, (8-116)
~nm u""n "m' wn 'm
ZLCL 1S 91ven by (8-84), Z2(x.:,to) is given by (8-87) and the summations in
(8-115a) are over all M rawin stations while in (8-1155) they are over all N
surface weather stations.
i -1
We apply the same |r interpolation formula to estimate 6z2 at x'. At
- ~m
all times t = t1} t2, etc. after t0, we will have estimates of z2(x,t) from
the prognostic equation 8-2 (output of Stage W2 described below). In this
case z2(xm,t) must be assumed to be exact (i.e., 6z2(x' t)=0). Errors are
allowed only at the end points of the time intervals that we treat. For
example, in tx < t < t2> we assume 6z2(tx)=0 but we allow finite values of
6z2(t2) at each rawin station.
Thus, step 3 of this stage (WEWC) is as follows:
(3) a. If to is the initial instant of the regional model simulation
period, determine Z2(xm,t0),m=l,...M from (8-87) (Stage ZQ) and
Z2(xl,t0) from (8-115). Then estimate 6z2(x ,t0),m=l,...M (see
step 4, Stage ZQ) and subsequently 6z2(x' t0) using these
values and (8-115).
b. If t0 is not the initial instant of the model simulation, then
Z2(x,t0) is available at all grid points x from Stage Z2,
described below; and 6z2(x^,t0) = 0 everywhere.
(4) We can now estimate the range of values (8-44) in which the
parameter AZ2 that enters in (8-45) lies. Recall from (8-34)
167
-------
that AZ2 is the difference in the z2 values measured at (x',t0)
and (x .ti). We find after some thought that
v~ffl i' a
- [Z2(x;,t0)
(8-117)
~m*
- I[z2(x' t0)
(8-118)
(5)
We now solve (8-45) for the minimum and maximum values of G
(see 8-46) in the vicinity of xm during the interval t0 < t <
~ — —
IOT(ti) and ^(x^) and
and IQm ^ 0;
, if
0;
0, if I(ti) 1 0.
Qm
= (8"119) with
replaced by
(8-119)
(8-120)
is from (8-117), AZ2C )max is from
In Eq. (8-119), AZ2(
(8-118), Iwm(ti) is from (8-109), IQm(ti) is from (8-111),
J (tx) is from (8-110), I( ) is from (8-113), qr( ) is from
Oui ^»
(8-114), and A is given by (8-30).
168
-------
(6) Next we must check whether the upper bound G on G is
iiiciX
positive. If it is not, we must return to Stage ZQ and alter
the value of q and, if necessary, the other parameters listed
in (8-47) until the resulting integrals I, q , etc. in (8-120)
\*
yield a value for G__v that is positive, and preferably of an
maX
order of magnitude consistent with (8-40). The procedure for
altering the variables in (8-47) must be developed through
experimentation (refer again to the paragraphs preceding
(8-47)) and therefore we will not attempt to define it here.
Thus, if G(x ,tt) > 0» 9° to steP ?'» otherwise return to
Stage ZQ and begin modification of q , etc. as described above.
(7) We can arrive at this step only after the previous analyses
have produced a range of G values that lies at least partially
on the positive, real axis. We assume, therefore, that
Go < GCt.) < 6M (8-121)
where
Go = •nax[0,G(xm,t1)m.n]. (8-122)
We now compute limiting values for the cloud updraft velocity
parameter w .
tu
FA(x ,ti) + GQ!(X iti,to)JLQ,.(\.,»ti,to)*
~m ~m c ~m
cm ^ , if q and * ^ 0; (8-123)
169
-------
where
t1-jAt)[3
qc(xm,t1-jAt)ac(xm,t1-jAt)
- 3] } (8-124)
In this expression M3 is from (8-94); z3 is from (8-91); 3,
q , and q^ are from (8-98); w(x,z,t) is the vertical velocity
function solution or stage WV defined earlier just before
(8-103); and z is defined by (8-108) with a similar
definition for z . .
sm
^(xm,t1)mav = Eq. (8-123) with G0 replaced by
c ~m max (8-125)
(8) If w (x ,ti) < 0 and LmC^i) and qc(x ;ti,t
then return to Stage ZQ and alter q^ and, if necessary, the
other parameters in list (8-51) until Eq. (8-125) yields a
positive value for w c(xin>ti)max- The alteration process is
similar to that discussed earlier in connection with the
calculation of G (see the paragraph preceding (8-51) and step 6
above. When a range WG( )min < wc < WG( )max has been found
that lies at least partially on the positive real axis, proceed
to step 9.
170
-------
(9) At this point a range of values of GCx^ti) has been
established (in step 7, Eqs. 8-121 and 8-122); and a range of
w(.(x[n,t1) has been computed in steps 7 and 8, namely
where the upper bound is positive. We now select from these
ranges the following values for site xm for the interval t0 to
G(x ,tj) = value in the range (8-121) that
~m satisfies (8-40) closest; (8-127)
wr(x _,ti) = solution of (8-49) with G
c ~m given by (8-127) and all (8-128)
other parameters with values
determined in step 8.
We assume that these values are constant at xm during the
entire interval t0 < t < ti, specifically, we assume
G(x t) = G(x ,
ra ^ t0 < t < ^ (8-129)
These values should be recorded for each hour in the interval
t0 to tt for site x . Note that tlt which is the hour of the
rawin observation following that at t0, is not necessarily the
same hour for all rawin stations. (See Figure 8-5). This does
not pose a problem because we apply (8-129) to each station
separately to obtain hourly estimates of G and w .
171
-------
(10) Repeat step 9 for each interval in the model simulation period.
This will produce hourly values of G and w at rawin station m
throughout the simulation period t0 to to + T.
(11) Repeat steps 1-10 for each rawin station in=l.,...M. At the
conclusion of Stage WEWC, there will be hourly values of G and
w at every rawin station throughout the simulation period.
Stage W2
We now begin to assemble the information gathered in the previous stages
to solve the prognostic equation (8-2) for z2 over the model space-time
domain. First, we construct the entrainment velocity field.
(1) Starting at the initial instant t0 of the simulation period, collect
the G values at this hour generated in Stage WEWC at each of the M
rawin stations and interpolate them onto the NEROS grid using the
.1
r weight:
G(x,t0) = - - - (8-130)
* |
m=l
~
~m
where x ranges over all grid points.
(2) Convert the G field into the w field as follows (see 8-39):
we(x,t0) = G(x,t0)Q(x,t0) (=f6/A0) (8-131)
172
-------
where Q is the surface heat flux in the grid cell centered at x at
time t0. Note that wg(x,t) (designated f0/A8 in Part 1), t0 < t <
t0 + T is an output of this processor, P8, for each grid point in
the regional model domain.
(3) We construct the initial z2 field at each grid point using the
interpolation scheme (8-115) employed earlier in Stage PATH.
Specifically,
(" Eq. 8-115a, if o (x,t0) = 0;
Z2(x,t0) = \ c (8-132)
(_ Eq. 8-115b, otherwise
where x ranges over all grid points in the model region.
(4) Interpolation of the cloud updraft velocities w estimated at the
rawin stations in Stage WEWC requires caution because it may often
happen that cumulus clouds are present in isolated areas that do not
contain an upper air station, or are present at station locations
between observation times. To handle these situations we propose to
compile a semi-empirical relationship between w and cloud depth
using all w , z2 and ZJQM estimated in the earlier analyses during
the simulation period T. This function can be based on field data.
For now compute the average of the w values obtained in Stage WEWC
for various values of the cloud thickness h3 = Zjp., - z2, say values
of h3 at intervals of 200m; and call the resulting function W(h3).
Then we assume
wr(x,t0) = (l-ff(x)Mz (x,t0)-z2(x,t0))
c ~ ~ ICU - - (8-133)
173
-------
where
W(ha) = w for clouds of depth ha (an empirical (8-134)
function);
a(x) = distance weighting function
(8-135)
=exP[-|x-xmi/SIG]
x is the rawin site closest to x where w 2 0; SIG is a distance
constant, e. g., 50km; and
wcLOCAL(*'to) = wc(~m'to)' (8-136)
Formula (8-133) is an heuristic expression that assigns w a value
V*
that is determined partly by any w estimates that have been derived
for that location in the earlier stages and partly by the assumed
empirical relationship between cloud depth and updraft velocity.
Note that wc(x,t), with x ranging over all grid points and t over
all hours in the simulation period T, is an output of this
processor, P8.
(5) For the advection velocity field v2H required in Eq. (8-2) at time
to, we will use
V2H(x,t0) = [22(x)to)-zt(x)f1 vH(x,z,t0)dz (8-137)
where v(x,z,t) is the horizontal wind at level z from the routine WV
described earlier. The horizontal velocities that WV provides
174
-------
should be those that are used to compute the vertical velocity
w(x,z,t) and they should be interpolated from the wind observation
stations to each grid point x using an interpolation scheme that
maintains a consistent relationship between w(x,z,t) and the
horizontal winds vu(x,zit) at elevations z. < z'< z.
~n ~ t — —
(6) The vertical air speed at each grid point x at the initial moment to
is obtained from the function routine WV as before, namely
w2(x,t0) = w(x,z2(x,t0),t0) (8-138)
where the function on the right side is a part of routine WV.
(7) We now solve (8-2) for z2(x,t0+At) using the difference scheme
described in Appendix A to Processor P7. With this scheme we have
Z2(x,t0+At) = Jp(x,to+At|x;t0)[z2(x;t0)
(8-139)
At S(x;t0)]dx'
where
cc
(l-ac(xU0)) + we(xlt0)
where a (x,t0) is the known fractional coverage of cumulus clouds at
time t0 in the grid cell centered at x;w2 is from (8-138); wc is
from (8-133); wg is from (8-131); and
p(x,tix',t') = known function of y2H- (see Chapter 9, Part 1)
(8-141)
175
-------
(8) At this point we have computed z2(x,to+At). From this field we
construct z3 in the manner of (8-91), namely
Z2(x,t0+At) + 100, if o-c(x,t0+At) = 0;
Z3(x,t0+At) = { (8-142)
z-j-pu(x i t0+At) , otherwi se
where z-p~y is the known elevation of cumulus tops. This calculation
should be performed at all grid points x and the results retained
for output from P8.
(9) The material surface flux H3 across surface H3 must be computed at
each grid point at time t0 + At as follows:
H3(x,t0+At) = z3(x,t0+At) + y3H-VHz3
-W(x,z3(x,t0+At),t0+At) (8-143)
where
Z3(x,t0+ At)-z3(x,t0)
Z3(x,t0+At) = - - - (8-144)
ysH(x,t0+At) = y(x,z3(x,t0+At),t0+At) (8-145)
and
v3u4VHz3 = calculated from v3H and z3 in the (8-146)
manner of (8-A2)7 See the Appendix
to this section.
In (8-143), w( ) denotes the function routine WV. The same applies
to v( ) in (8-145). The flux H3 is an output of processor P8 but it
is not used in the calculations performed within this processor. At
the initial moment t0, assume
176
-------
H3(x,t0) = HaCx.to+At) (8-147)
(10) We must also output the local time derivative z2 at each grid point
and each hour:
z2(x,t0+At) - z2(x,t0)
Z2(x,t0+At) = . (8-148)
assume
0) = z2(x,t0+At). (8-149)
These fields are outputs of P8 at each hour.
(11) Compute the thickness of Layer 3 at each grid point at time t0+At:
h3(x,t0+At) = Z3(x,t0+At) -Z2(x,t0+At) (8-150)
where z3 is from (8-142) and z2 is from (8-139). Evaluate (8-150)
for output only when t0+At is an integral number of hours.
(12) Compute the thickness of Layer 2 for output:
h2(x,t0+At) = Z2(x,t0+At) - ZT(X) - Mx.to+At) (8-151)
\
where z2 is from (8-139) and ZT and hi are inputs to Processor P8.
Compute (8-151) only at integral hours.
(13) Compute the elevation of surface z3 in pressure coordinates.
P3(x,t0+At) = p(x,z3(x,t0+At),t0+At) (8-152)
177
-------
where p(x,z,t), which is given by (8-86), is the pressure (mb) at
elevation z (m MSL) over site x at hour t. This function is
interpolated from the pressure-height functions Pm(x,t) that are
inputs from Processor PI. The field p3 is an output of this
processor, P8, at each grid point and hour.
(14) Compute the elevation of surface z2 in pressure coordinates:
P2(x,t0+At) = p(x,z2(x,t0+At),t0+At) (8-153)
where p is given by (8-86). The field p2 is an output of this
processor at each hour of the simulation period.
(15) In Part 1, Eq. 5-44, we defined the function >H(x,t) to be the
fraction of the volume flux entering cumulus clouds in the grid cell
centered at x at time t that originates in Layer 0. Here we compute
an interim estimate V* of this parameter based strictly on heuristic
notions. In Processor P12 we test whether the I'1 values derived
here satisfy criterion (5-47) of Part 1. If they meet this
condition they become the final estimate of the parameter 4*.
Otherwise y(x,t) is given the largest value that satisfies (5-47).
As a rough, heuristic approximation we shall assume
2h3(x,t0+At) 2
¥' (x,t0+At) = 0.5 exp -[- 3
Under this assumption, ¥ -» 0.1 as the depth of cumulus clouds
becomes large compared to the depth (z2-Zy) of the subcloud layer.
For shallow cumulus, t -» 0.6. The H*1 field is an output of
Processor P8 each hour.
178
-------
Stage WV
This stage is in effect a function routine that provides vertical and
horizontal wind speeds at any grid cell and elevation, and horizontal
divergence between any two surfaces in any grid cell.
As with nearly all of the input parameters required in the regional
model, the variables just cited can be estimated by a number of methods.
Solutions of the omega equation, isentropic analyses, and differentiation of
functions fit to wind data are several approaches. Our intent here is not to
prescribe the use of one particular technique but rather to specify which
quantities are required for subsequent use in this and other processors in the
network, and in the regional model itself. Selection of specific techniques
is left to the "user". Indeed, one of the reasons for structuring the model
system in a modular form was to facilitate the use of various techniques
interchangeably. In the "first generation" processor network we plan to
employ in this stage, WV, the method developed by Bullock (1983).
In order to insure compatibility between the horizontal divergences <6>
required in Processor Pll and the vertical velocities w required in this and
other processors, it is advantageous, to compute all these fields in this one
stage. Following is a list of the divergences that must be computed each hour
for output to Processor Pll:
P3(x,y,t)
<6(x,y,t)>3 = p^- (VH-v)dp (8-155)
3 2
179
-------
P2(x,y,t)
1 ^(V^dp , if APl(t)=l;
<6(x,y,t)>2 = "
Pvs(x,y,t)
(8-156)
P2(x,y,t)
(V^)dp ' if ^i^0
Px(x,y,t)
<6(x,y,t)>1 = p .* f (VH-v)dp , defined only v/heri Ap^O (8_157)
Pvs(x,y,t)
Since the input wind data are in z coordinates, it will be necessary to
convert them to p coordinates using the function p (z,t) available from PI.
Note also that since Po^p. P?0 and the reciprocals
of the pressure differences that appear outside these integrals are also
negative. The net result is that <6> has the same sign as VM-V. We leave the
method of computing the divergences (V,,'V) to the user.
The layer averaged divergences <6> given by (8-155, 156 and 157) are
related to the vertical velocity wn on the pressure surfaces p that bound
each layer by (see Eqs. 11-lb and 11-40)
180
-------
<6(x,y,t)>n = ji- [u»n(x,y,t) - sT^x.y.t)] (8-158)
n
where
Apn = |Pn(x'y't} " Pp-i^y^l (8-159)
and ui denotes the vertical velocity in pressure coordinates (mb sec )
averaged over the surface p within the grid cell centered at (x,y). We shall
assume that the spatial variations in u> have such large scales that
wn(x,y,t) - u>n(x,y,t) (8-160)
where u> denotes the local value of u> on pressure surface p at the cell
n rn
center (x,y).
By definition
where p is pressure and w is the vertical velocity in z coordinates at the
point where p is measured. Making use of the hydrostatic approximation and
(8-161), we can write
% = 5? Pn+ Wn - Vn* (8-162)
Two of the fields that we need in this processor, P8, are w^ and u^, the
vertical velocities in z coordinates (m sec ) on the top surfaces of Layers 3
and 2. Eqs. (8-158), (8-160), and (8-162) provide a means of estimating these
velocities in terms of the measured horizontal winds and their divergence.
181
-------
We must distinguish between the two cases Ap,=l and Ap,=0. Recall from
(7-98) that Ap-^t^l if a surface inversion is present over the modeling
region at time t and it is zero otherwise. When an inversion is present the
effective ground surface for the flow aloft is the "virtual surface" p
defined as the top of the inversion layer.
Mode 0: Ap, = 0.
In this case the virtual surface coincides with the terrain. Using
(8-158) and (8-159) we can express u in the form
+ u>Q (8-163)
where uu is the value of u> at ground level.
To estimate ui we note first that the vertical velocity w at ground
level is
W0 = WT (8"164)
where ZT is the terrain height and v is the horizontal wind velocity near the
ground. Since horizontal gradients in terrain level are generally much larger
than those of the synoptic scale pressure field p , we expect that the last
term on the right side of (8-162) will have a much larger magnitude than the
term involving Vu*p . We also expect the last term to exceed the magnitude of
~ti o
the local time rate of change of p . Thus, to good approximation we can
assume
u> - -p gv -VHZT. (8-165)
182
-------
W3 = p [-H + *3-V3 " Ap3<5>3 '
Combining (8-162), (8-163), and (8-165) we have
9p,
'• + ^3*^HP3 " AP3<6:
mode 0. (8-166)
J. J. U ~U ~~ll I
where
Apj = jp1(x,y,t) - Pvs(x,y,t)| , mode 0. (8-166a)
We assume the velocity and divergence fields required to evaluate this
equation have been extracted from the wind data that are inputs to this
processor. The pressure surfaces p are specified by those stages within P8
that call stage WV for iu estimates.
To estimate the densities p3 and p , we use
P(x,t)
n
"n' RTTxTt) <8-167>
n "*
where R is the gas constant and T is the temperature at pressure level p and
grid point x = (x,y). The temperature fields will be an input to this
processor. To evaluate VHp3 and V^ in Eq. (8-166), it will be necessary to
fit a polynominal of at least 3-rd order to the pressure andi terrain data (see
Appendix to this chapter).
The expression for w« is obtained in the same manner that (8-166) was
derived. We get
183
-------
1 ™2
wo = TTZ ["of + vo*vuPo ~ Ap0<6>0 -
2 P/>g 3t ~2 ~Hr2 r2 2
(8-168)
+ p gv *V,,zT], mode 0
Also in mode 0 we require for direct use in the regional model the
component of the vertical velocity w, on surface p., that is due to horizontal
divergence of the flow in Layer 1 (see Eq. 3.1b of Part 1). We denote this
velocity component by wn, where
wl = WT1 + "Dl = pTg [8T + *l-Vl • Apl<6>!
1 (8-169)
Since by definition w ^ is the divergence induced vertical motion on p,,
the terrain induced component is just the right side of (8-169) with <6>,=0.
Thus,
AP-,
wni = ' rr <6>i » mode 0- (8-170)
UJ- P^9 -1-
and Ap, is given by (8-166a). In order to emphasize that this value applies
in mode 0 only, we designate it wi, at the output of Pll. Processor P12 will
take this field as input and generate the final, complete function w,,,
184
-------
Mode 1: Ap., = 1
In this state an inversion layer is present at the ground and the
effective bottom of the flow aloft is the virtual surface p , which coincides
with terrain features that extend up through the inversion layer but which
otherwise is the top of the inversion layer itself. In this situation the
expression for u>_ becomes
+ u)vs. (8-171)
In order to obtain an estimate of uu that is consistent with the model used
in P7 to simulate flow in the ground level inversion, we must impose the
constraint that the volume flux across p is continuous.
Just below the top z of the cold inversion layer the downward volume
flux is
dt
where (u ,y ,w ) is the fluid velocity in the cold layer; and nvs is the
cooling rate, a known function of space and time. Just above zys, i.e., at
the base of Layer 2, the downward volume flux is
Continuity of the flux across HVS requires that (8-172) and (8-173) be equal,
which is true if
185
-------
3zvs azvs azvs
If we assume that the spatial gradients in z are also much larger than those
in the synoptic pressure field, we can make use of the earlier approximation
(8-165) to obtain
wvs " "Pvs^vs (8-175)
where w is given by (8-174).
Combining (8-162), (8-171), (8-174), and (8-175) we get
i 3P.
P3g
3
Csr5 + v~-Vupo - Ap,<6>- - Ap9<6>«]
o t» ^o ^n o o o £ £
--n]. model (8-176)
and similarly
W2 =
(8-177)
+ [57 + % *£lJZ ~ ^ ^ » mOC^e 1
In both (8-176) and (8-177)
Ap2 = |p2(x,y,t) - pvs(x,y,t) | , mode 1 (8-178)
In mode 1 no estimate of w at the top of layer 1 is required because in this
situation the volume flux is given by n.vs-
186
-------
The input requirements and outputs of Stage WV are summarized in Table
8-1. Figure 8-6 illustrates processor P8 and its data interfaces
schematically.
187
-------
Table 8-1. Summary of the input and output requirements
of each stage of Processor P8.
Input
Variable
Description
Source Stage
Output
Variable
Description
qm(z,t)
em(z,t)
qn(t)
en(t)
mixing ratio at PI
elevation z (m MSL)
at hour t over rawin
station m.
potential temperature PI
(°K) at elevation z
(m MSL) at hour t over
rawin station m.
ZQ
mean terrain ele-
vation (m MSL) in
grid cell centered
at x .
P7
mixing ratio (dimension- P3
less) at surface weather
station n at time t.
potential temperature P3
(°K) at surface weather
station n at time t.
elevation (m MSL)
of top surface of
model layer 2 over
rawin site x at tim
t of rawin soundings
expected range of
values, centered at
at which actual mixe
layer top lies over
rawin site x at
times t of rawin
soundings.
elevation (m MSL) of
top surface of model
Layer 3 over rawin
station m at times t
of rawin soundings.
estimated water mixir
ratio entering cumuli
clouds over rawin
station m at time t
(available at 30 mini
intervals).
Pm(z,t)
pressure (mb) at PI
elevation z(m MSL)
at time t over
rawin station m.
fraction of sky RAW
covered by cumulus
clouds in grid cell
centered at x at
time t.
H.QS.,1)
estimated water
mixing ratio above
Layer 3 at rawin
station m at time t
(available at 30
minute intervals).
vertically integrated
water mixing ratio in
Layer 3 over rawin
site m (units=m, see
Eq. 8-20) at times t
of rawin soundings.
188
-------
Table 8-1. (continued)
Input
Variable
Description
Source Stage
Output
Variable
Description
zTCU(x,t)
elevation (m MSL)
of highest cumulus
cloud tops in grid
cell centered at x
at time t. ~
RAW
(Cont.)
average mixing ratio
in Layer 3 over rawin
site m at time t
(available at 30 min.
intervals).
elevation (m MSL) of
the lifting condensa-
tion level over sur-
face station n at
time t (required
only at the hours
t of rawin observa-
tions)
w(x,z,t)
QCx.t)
horizontal wind Stage
vector (m/sec) WV
at site x, ele-
vation z (m MSL)
at time t.
vertical air speed Stage
(m/sec) at site WV
x, elevation z (m
MSL) at time t.
surface kinematic i P4
heat flux (m°K sec" )
in cell at x at time
sv
t.
fraction of sky RAW
covered by cumulus
clouds in grid cell
centered at x at
time t.
elevation of top Stage
surface of Layer 2 ZQ
at rawin station m.
PATH
path integral of w
leading to rawin
site m at time t
(see 8-109) units
= m.
path integral of
cumulus cloud cover
fraction leading to
rawin site m at hour
t (see 8-110) units =
sec.
path integral of
kinematic surface
heat flux leading
to site x at time t
(see 8-1H) units =
m °K.
189
-------
Table 8-1. (continued)
Input Output
Variable Description Source Stage Variable Description
average mixing ratio Stage WEWC G(x ,t) entrainment velocit;
in Layer 3 over ZQ ~ scale factor (units
rawin station m at °K~1) (see 8-39) at
time t (available at rawin site x at
30 min. intervals). hourly intervals t.
mixing ratio Stage wc^~mjt^ cumulus updraft
entering cumulus ZQ ~ velocity scale (m/
clouds over site sec) at rawin
x at time t site x at hourly
TSvaiTable at 30 min. intervals t.
intervals).
surface kinematic P4
heat flux (m °K
sec- ) in grid
cell at x at time t
(at hourTy intervals).
fractional coverage RAW
of cumulus clouds
cumulus cover path Stage
integral PATH
Z2(>c,t) elevation of top Stage
surface of Layer 2 ZQ
at rawin site x
at times t of ~m
rawind soundings.
z, r,(x ,t) lifting condensation Stage
LLL ~n ievel at surface ZQ
weather station n
at times t of rawin
soundings.
T. (t) w path integral Stage
*"" leading to rawind PATH
site x at time t
of rawrn sounding
In (t) surface heat flux Stage
gm path integral PATH
leading to site
x at hour t of
rawin sounding.
190
-------
Table 8-1. (continued)
Input
Variable
Ms(xfn,t)
Description
vertically inte-
grated water mixing
Source
Stage
ZQ
Stage
WEWC
(Cont.)
Output
Variable
Description
w(x,z,t)
ratio in Layer 3
over rawin site x
(at sounding hours t
only).
vertical air speed
(m/sec) at (x,z,t)
Stage
WV
water mixing ratio Stage
above Layer 3 at ZQ
rawin site x (at
t = 30 min Tfltervals).
elevation of top of Stage
Layer 3 at rawin ZQ
station m at sounding
times t.
expected range of Stage
z2 at rawin site ZQ
x at observation
K8ur t.
zTCU(x,t)
hourly values of Stage
the entrainment WEWC
velocity scale factor
at rawin site xm.
~m
hourly values of Stage
the cumulus updraft WEWC
speed (m/sec) at
site xm.
~m
cumulus top RAW
elevation at hour t
in grid cell
centered at x.
estimated top of Stage
Layer 2 at rawin ZQ
site x at sounding
time t.
W2
w (x,t) = entrainment velocity
f§/A6(x,t) at hour t at grid eel
centered at x (m/sec)
w (x,t) cumulus updraft speed
(m/sec) at hour t in
grid cell centered at
x.
vertical air speed
(m/sec) on surface
H2 at grid cell
centered at x at hour
t.
Z2(x,t) elevation (m MSL) of
~ surface H2 (top of
Layer 2) in grid cell
centered at x at hour
t.
191
-------
Table 8-1. (continued)
Input
Variable
hW
Description
mean terrain
elevation (m MSL)
in grid cell
centered at x.
Source Stage
P7 W2
(Cont.)
Output
Variable
z3(x,t)
Description
elevation (m MSL) o-
surface HS (top of
Layer 3) in grid ce'
centered at x at hoi
t.
w(x,z,t)
P.(z,t)
hi(x,t)
horizontal wind Stage
vector (m/sec) WV
at site x, elevation
z, time t.
vertical air speed Stage
(m/sec) at site WV
x, elevation z (MSL)
at time t.
pressure (mb) at PI
elevation z (m MSL)
at hour t over rawin
station m.
depth (m) of Layer 1 P7
in grid cell centered
at x at time t.
H3(x,t) volume flux (m/sec)
through top surface
Layer 3 in grid cell
at x at hour t.
/N/
Z2(x,t) local time derivativ
of elevation Z2
(m/sec) at grid
cell centered at x
at hour t.
hs(x,t) thickness (m) of
Layer 3 at hour t
in grid cell centere
at x.
**
h2(x,t) thickness (m) of Lay
2 at hour t in grid
cell centered at x.
cv
Ps(x,t) pressure (mb) at
elevation of top
of Layer 3 at hour
t in grid cell
centered at x.
P2(x,t) pressure (mb) at
elevation of top of
Layer 2 at hour t
in grid cell centerei
at x.
^'(x,t) interim estimate of
the cumulus flux
partition function
(see Eq. 5-44 of
Part 1) at hour t
in grid cell at x
(non-dimensionalT
192
-------
Table 8-1. (continued)
Input
Variable
Description
Source Stage
Output
Variable
Description
Qn(t)
On(t)
Ap,(t)
P3(x,t)
terrain elevation
(m MSL) in grid cell
centered at x.
P7
WV
observed east-west PI
wind component (m/
sec) at elevation
z (m MSL) at observa-
tion hour t at rawind
station m.
same as u except PI
north-south wind
component
observed east-west P3
wind component (m/
sec) at observa-
tion hour t at sur-
face weather station
n.
same as u (t) except P3
north-south wind
component.
surface inversion P7
indicator (see 7-98).
elevation (m MSL) at PI
pressure level p (mb)
at rawin station m
at hour t.
elevation in pres- Stage
sure coordinates (mb) W2
of top surface of
Layer 3 in grid cell
at x at time t.
v,,(x,z,t) horizontal wind
vector (m/sec)
[VM=(U,V)] at (arbitrc
elevation z (m MSL),
time t at site x.
w(x,z,t)
W[),(x,t)
<6(x,t)>3
<6(x,t)>2
<6(x,t)>1
vertical air speed
(m/sec) at (arbitrary.'
elevation z (m MSL)
at time t at site x.
divergence induced
vertical air speed
(m/sec) on top surfaa
of Layer 1 (Defined
for daytime hours
average horizontal wii
divergence (sec-1) in
Layer 3 in grid cell
centered at (x) at
hour t.
same as <6>3 except
applies to Layer 2.
Same as <6>3 except
applies to Layer I
(values of this
quantity are computed
only for daytime hour
same as <5> except
applies to Layer 2.
same as <5> except
applies to Layer I
(values of this
quantity are computed
only for daytime hour
P2(x,t) same as p3 except Stage
elevation of top W2
of Layer 2.
193
-------
Table 8-1 (Concluded)
Input
Variable
P^x.t)
Pvs(x.,t)
Description
same as p3 except
top of Layer 1.
same as p3 except
Source Stage
P7
P7
Output
Variable
Description
-vsv~'
nvs(x,t)
inversion
elevation (m MSL)
of virtual surface
in cell at (x) at
hour t. Note: z
when Ap-,
= 0.
vs
growth rate (m/sec)
of the radiation
inversion layer
depth.
P7
P7
194
-------
Appendix to Section 8
Some of the equations in this Processor, such as (8-166), contain terms
of the form
v-VHp (8-A1)
that must be evaluated at each grid point (I,J) of the regional model domain.
In (8-A1) both v and p are variables represented in discrete form at each grid
point of the model region. To fourth-order accuracy we can represent (8-A1)
at grid point (I,J) by
vI,JAy(pI,J) (8"A2)
Ax(pI,J} = 3 (PI+I,J " PI-1,J} ' 12 (PI+2,J " PI-2,J}
where
and
\i (PT i) = T (PT 1+1 " PT 1-1) " T? (PT 1+7
y -Ij" -5 1,J*J. 1,J 1 LL i,J+t
Here PT , is the value of p at point (I,J), i.e., column I, row J, with rows
i,j
parallel to the x axis and columns parallel to y.
195
-------
YYYYYYYYYYY YYY
Q.
•*•*
3
O
•o
03
+*
3
Q.
C
tn
+*
•a
c
to
00
o
O v.
o o
o a
.
*± 0)
ra o
ECO
„ H-
o ^
OJ -E
CO
00
(D
3
0)
196
-------
SECTION 8
PROCESSOR P9
DEVELOPMENT
This processor prepares the information necessary to correct the chemical
rate constants for variations in atmospheric density, temperature, cloud
cover, and solar zenith angle. Often the top of the regional model will lie
near the middle of the troposphere and therefore significant variations in air
density and temperature can exist between each of the model's layers. All
rate constants for the intermolecular reactions are affected by density and
many are strongly sensitive to temperature as well. The photolytic rate
constants are affected by the variations in solar radiation induced by cloud
scattering and absorption, and by the variations in radiation that accompany
changes in the solar zenith angle. We treat the density and temperature
correction terms in Stage DEN, and the cloud and zenith corrections in Stage
ZEN.
Stage DEN
The easiest way to handle the effects of atmospheric density variation on
the pollutant concentrations is to work with the species mass continuity
equation in its mixing ratio form. Thus let
Y(r,t) = c(r,t)/p(r,t) (9-1)
197
-------
be the mixing ratio of a given species, where c is the species concentration,
_3
say, moles m , and p is the local air density in the same units. The mixing
ratio y is dimensionless and is often expressed in terms of parts per million
(ppm). Making use of (9-1) in the species mass continuity equation (Eq. 2-1
of Part 1) we have
§| (w) + V(WY) + 9z (Ypw) = s + R" w- (9"2)
The mass continuity equation for the atmosphere has a similar form, namely
at P + YH-(pv) + §1 (Pw> = °- <9'3)
Multiplying (9-3) by -y and subtracting the result from (9-2) we obtain
&+ *„•(«)+ st 103m) that air density variations are important.
198
-------
Normalization of the emissions function S by the air density is performed
by Processor P10. Here we are concerned with the normalization of the
chemical reaction term R, which has the general form
R = kcd + kjc (9-5)
where d is the concentration of a species with which the given pollutant
-1 0 -1
reacts, k is the rate constant (units of mole nrsec ) of this second-order
reaction, and k^ is the rate constant for a first-order reaction that consumes
species c. We can write (9-5) in the form
(9.6)
where
Yd = d/p
is the mixing ratio of species d. Substituting (9-6) into (9-4) we have the
equation governing y- In it the second-order rate constants have the form
k* = kp. (9-7)
but the first-order rate constants are unchanged. In order for the model
equations to predict mixing ratio, we must supply the kinetics algorithm with
effective air densities for each layer so that the rate constants k can be
modified in the manner of (9-7). For this purpose, and also for normalizing
the emission strengths S in Processor P10, we compute in this stage average
air densities for each of the model's four layers. Concurrently, we compute
average temperatures for each cell and layer for us? within the model itself
to make temperature corrections on the rate constants. The rationale for
supplying temperature data to the model rather than temperature corrected rate
constants is to avoid any procedure that would "hard wire" a particular
chemical kinetics scheme into the model or the input processor network.
-------
Let Tvm(z,t) and pn(z,t) be the virtual temperature and density,
respectively, at elevation z(MSL) at hour t over rawin station m, and let z
be the elevation (m MSL) of that station. Now define
Sm om
P.todz (9-8)
A <
zm
zm + hon, + hlm
hom
hon, + hlm + h2n,
hom
where
- 1000
~ 2O7
_3
is the factor necessary to change the units of density from (kg m ), the
_3
units of Pm(z) as supplied by PI, to (moles m ), the units in which the
chemical rate constants that we will modify are expressed. In (9-8) - (9-11)
h- is the thickness of model layer j in the grid cell that contains rawin
station m.
200
-------
The integrals on the right-hand side of Eqs. (9-8) - (9-11) can be
evaluated by subdividing the integration interval into 50m subintervals that
coincide with the intervals at which TVI|)(Z) and PmU) are available (from PI).
Similarly, we define
(m om 1m
Tvm(z)dz (9-12)
hom
(9-13)
•
m "'* 3m
T™(z)d* (9"14)
m + hom +hlm *h2m
The temperature and density profiles T and p are available from PI
each hour but they give values only at the measurement station sites xm.
~m
Therefore, the layer averaged density and temperature values derived from
_i
(9-8) - (9-14) must be interpolated to each grid cell location. Use the r
weighting scheme as follows:
_ _
10 I r nm
n = S^j - n = 0,1,2,3 (9-15)
I r"1
m=lm
= -^ n = 1,2,3 (9-16)
IN n JLJ
™ -1
I rm
m=lm
201
-------
where
]35. (9-17)
-6
The factor 10 in the density equation (9-15) is necessary when the mixing
ratio Y is expressed as parts per million (ppm). The factors generated
by this stage will pass through the B-matrix compiler (BMC) unchanged and will
be used in the model to modify the second-order rate constants k in the manner
described above in Eq. (9-7), namely
kj = k n (9-18)
where k* is the modified rate constant in Layer n. (Keep in mind that the
first-order constants are not modified by the density.) As we noted above, we
assume that the second-order rate constants k are expressed in units of (moles
-3 .1. _1
m ) sec . (Consistency of the concentration units throughout the model
should be confirmed by comparing the parameters generated in this processor,
the source strength functions provided by Processor P10, the initial
concentration fields produced by Processor P2, and the rate constants
contained in the chemical kinetics subroutine CHEM that operates in unison
with the model.)
The layer averaged temperature values given by (9-16) above also
pass unmodified into the model at operation time. These data are made
available there for temperature corrections to those rate constants that
require it.
The input-output summary of stage DEN is given in Table 9-1, and the
processor is illustrated schematically in Figure 9-1.
202
-------
Stage ZEN
All of the photolytic rate constants require adjustments for the local
solar zenith angle and cloud cover. We assume following Jones, et al
(1981) that these "constants" can be expressed in the form
ka= ys)H(cc) (9"19)
where k ($ ) is the "clear sky" photolytic rate constant for reaction a and E
Of o
(cc) is an empirical function of the cloud cover cc = cc(h), which is the
fraction (0 < cc < 1) of the sky covered by clouds of height h. We assume
further that the clear sky constants k are contained in the chemical kinetics
subroutine of the regional model code in the functional forms k ( ) and that
the solar zenith angle <)>s must be supplied for each grid cell and hour to
evaluate them (note that is virtually independent of altitude). Thus,
Processor, P9, particularly Stage ZEN, must supply the following fields to the
model by way of the model input file MIF:
(I,J,t ) = solar zenith angle in cell (I,J) at hour t
s m (units of degrees of arc); m (9-20)
and
E(I,J,t ) = cloud cover correction factor for photolytic
m rate constants in cell (I,J) at hour t (9-21)
(dimensionless).
The zenith angle <|>s can be obtained from standard astronomical formulas given
the latitude (JA) and longitude (IAX) of cell (I,J), and the date and hour tm
of the period of interest. The factor E is obtained from the formulas given
in Jones, et al (1981) using the observed cloud cover cc(h) in cell (I,J) at
hour t . We will not elaborate on the computation of s and E here.
203
-------
In summary, the solar zenith angle <|> is used in each grid cell and hour
to determine the clear sky rate constant k for each photolytic reaction a.
These values are multiplied in turn by the cloud cover factor H for that cell
and hour to arrive at the corrected rate constant
"Stf.J.V = W^'V^^'V- (9~22)
The photolytic constants are not modified by the density correction terms
-•- n
generated in Stage DEN.
The inputs and outputs of Stage ZEN are summarized in Table 9-1.
204
-------
Table 9-1. Summary of the input and output variables
of each stage of Processor P9 and their sources.
Input
Variable
Description
Source Stage
Output
Variable
Description
virtual temperature
(°c) at elevation z
(m MSL) at hour t.
over rawin station
PI
DEN
0
Pm(z,tk) same as T except
m K density (Kgm-3)
m
elevation (m MSL)
of rawin station m
h (I,J,t.) thickness (m) of
Layer 0 at time
t. in grid cell
PI
RAW
P7
2
3
density correction
for rate constants
in Layer 0, cell
(I,J) hour t|< (units
10~6 moles nr3)
same as except
applies to Layer 1
same as except
applies to Layer 2
same as except
applies to Layer 3
average temperature
(°c) in Layer 1, ce"
(1,0), hour tk (for
temperature correct'
of rate constants ii
model)
hn(I
h2(I
h,(I
*3
,J,t
,J,t
,J,t
\('
l\
k)
i<)
N
same
Layer
same
Layer
same
Layer
as
1
a-s
2
as
3
h0
U
"o
h0
U
except
except
except
P7
P8
P8
9 same
K ^ . Layer
, same
K d Layer
as
2
as
3
-, except
J.
-, except
cc(x ,h,t ) fractional sky
coverage of clouds
of height h (low,
middle, and high)
at surface weather
station n at hour
t
RAW ZEN
-------
h-
D
a.
YYYYYYYY
A
OJ
A
«x
V
CO
A
Q.
V
A
V
A
V
CO
A
V
T3
C
«3
S- O
o s
CO -M
co
-------
SECTION 9
PROCESSOR P10
DEVELOPMENT
Processor P10 transforms the emissions inventory into the source strength
fv
functions S, S^ and $2- It is assumed that the emissions data have already
been structured to provide the following information: (1) total emission rate
(moles hr ) of each primary pollutant from all mobile and minor stationary
sources in each of the model's grid cells, at each hour of the simulation
period; (2) emission rate (moles hr" ) of each primary pollutant each hour
for all major point sources; and (3) the physical parameters necessary to
determine in which grid cell each major point source lies and what its
effective source height will be under given meteorological conditions. Here
"major" point source refers to any point source whose discharge rate of a
given pollutant exceeds some prespecified threshold. This processor also
computes the plume volume fraction £ which is used in the Layer 0 equations to
parameterize subgrid scale chemistry effects.
The three basic tasks involved in computing the source strengths are
estimating the effective heights of the major point sources; partitioning the
point and area source emissions among the three layers 0, 1, and 2; and
converting the results to units of ppm. The first of these operations is
performed in Stage DELH — the last two are done in Stage S. We should add
here that the units conversion is necessary because concentration must be
207
-------
expressed in ppm in the governing equations in order to account for the
variation of atmospheric density with elevation (see the discussion presented
with the description of Processor P9). The calculation of the parameter £ in
this processor is performed in Stage ZTA.
Stage DELH
Suppose that there are a total of K major point sources in the entire
modeling region, and let the subscript "k" designate any one of them. We
compute the plume rise Ahu(t ) of the k-th source at hour t as follows.
Ix III In
Let [I(k),J(k)] be the grid cell (column I, row J) in which the k-th
source lies. Then
0, if L[I(k),J(k),tm] < 0;
Ahk(V = i F. 1/3 (10-1)
3 (-|) , otherwise.
Here L[I,J,t ] is the Obukhov length (meters) in cell (I,J) at hour tm> and
VTak 2
F,, = g[ 3 (%)0kw. . (10-2)
K Tk + 273 K K
In this expression T. is the exhaust temperature (°C) of the stack, Dk is the
-1 -2
stack diameter (m), w. its exhaust velocity (m sec ), g = 9.8 m sec is
gravity, and
Z_/nk Tvi/V
T = Jt± (10-3)
K| — 0
208
-------
is the estimated surface air temperature (°C) at the location of source k. In
(10-3)
(xn,yn) is the location of surface weather station n and T is the virtual
temperature (°C) measured at that station at hour t .
In Eq. (10-1)
2
u = [(,)
where -. and , are the Layer 1 averaged wind components (m sec ).
Finally, the parameter s in (10-1) is approximated by
(10-5)
VTak + 273)
2 -2 n -1
where c = 1003 m sec K . This expression for the stability parameter s
assumes an isothermal temperature lapse rate throughout the depth of the layer
that the buoyant source emissions traverse.
The effective height of source k at hour t is now estimated to be
where zsk is the stack height (m) of source k.
By virtue of the assumption embodied in (10-1) for the case of negative
L, expression (10-6) yields effective source heights equal to the actual stack
209
-------
height z under neutral and unstable atmospheric conditions. Although plume
rise does occur under these conditions, it is not particularly important in
our regional scale grid model because under these conditions vertical mixing
is so strong that pollutants are nearly uniformly spread through the mixed
layer before horizontal transport has moved them out of the grid cell in which
they were released.
For example, with a horizontal grid size of 18,000 m arid a horizontal
wind speed of 4 m sec , material has a residence time of 4500 sec in a grid
cell. The time scale of vertical mixing during unstable conditions is on the
order of 2h/w*. Typical values of the mixed layer depth h and velocity scale
w* are 1500 and 1.5 m sec , respectively. Hence, vertical mixing is
completed within the time material is resident in the cell surrounding the
source and consequently the actual height of the emission is not significant.
(A much more important factor in this instance is the lack of complete
horizontal mixing of point source plumes within a grid cell. Although the
material is well mixed vertically, it may occupy a volume only 1 km wide
whereas the grid model treats the emissions as though they fill the entire
cell uniformly. This discrepancy is a subgrid scale phenomenon that the
current regional model does not treat in Layers 1, 2 and 3. It is potentially
a source of considerable error in the photochemical reaction simulations that
should be considered in future modeling applications.)
In stable conditions, vertical mixing is very weak or nonexistent and in
that case point source emissions must be placed in the proper layer.
Table 10-1 summarizes the input and output parameters associated with
Stage DELH.
210
-------
Stage S
Having estimated the effective discharge heights of the major point
*M
sources in Stage DELH, we can now compute the source strength functions S, S^
and $2-
In this task we must take into account that the ground surface can extend
into Layers 1 and 2 where hills penetrate the model layer surfaces H and H-,.
Figure 10-1 illustrates the relationships between the layer interfaces and the
terrain. (See Fig. 4-7 of Part 1. Recall that a^ is the fraction of
surface H-^ that is penetrated by terrain, and that a-™ is the corresponding
fraction for surface H .)
Let
be the emission rate (moles h ) of a given primary pollutant
from the k-th major point source at hour tm, and let E(I,J,t ) be the total
emission rate (moles hour ) of that pollutant from all other sources in grid
cell (I,J) at time t . The source strength functions for each primary
pollutant have the following forms.
H2
LAYER 2
Figure 10-1. Illustration of the influence of terrain on model layers 0, 1 and 2
for given values of the penetration fractions aT1 and ovn.
211 IJ- IU
-------
K
2 = [crT1E + £ Ek Uk2]/(A2h2) (10-7)
k=l
K
1 = C(aTO-aT1)E + £ Ek ^/(A^) (10-8)
k=l
K
S(I,J,tm) = [(l-aTO)E + £ Ek Uko]/AQ (10-9)
k=l
where aTQ, aT1, E, h-^, and h2 are functions of (I,J,tm); and Ek and II.
depends on t . The latter is defined by
U
where
fl, if (z. i - ZT) < H < (Z. - ZT)
,(tm) = J 1 I - * - J l (10-10)
J (_0, otherwise.
.j = zT(I(k),J(k)).
zo =
zT(I(k),J(k))
= h1(I(k),J(k),y
= h2(I(k),J(k),tm)
and h is the thickness of Layer n. Also
212
-------
A2 = A (10-11)
A! = (1 - an)A (10-12)
AQ = (1 - aTO)A (10-13)
A = a2(cosih)(6X)(6$) (10-14)
o
a = 6,367,333 m = earth radius
6A = l/4() = longitude grid interval
6 = 1/6(3!^) = latitude grid interval
and <|>j is the latitude of row J. In the evaluation of Eq. 10-9, the total
number of major point sources lying in Layer 0 should be determined in each
grid cell for later use in Stage ZTA. Let us denote this variable
The source strength functions derived from (10-7) and (10-8) have units
-3 -1 "" -? -1
of (moles m hr ) and (10-9) gives S in units of (moles m hr ). These
must be converted to units of (ppm sec ) and (ppm m sec ), respectively.
This is done as follows:
=
^>'tm)>n 1S tne average air density (units = 10 moles m ) in
Layer n, time t in cell (I,J). Eqs. (10-15) - (10-17) should be evaluated
213
-------
for each primary pollutant at each hour, level and grid cell; the results
should be recorded in the MIF. The inputs and outputs of Stage S and their
sources and destinations are summarized in Table 10-1.
Stage ZTA
The plume volume fraction parameter £(I,J) is used in the Layer 0
equations to parameterize subgrid scale chemistry effects (see Chapters 5 and
8 of Part 1). Here we use a very simple approximation to estimate it.
Let L(I,J) be the total length (m) of major line sources in cell (I,J),
let N (I,J) be the total number of minor point sources in cell (I,J), and let
N .(I,J) be the total number of major point sources that lie in Layer 0 in
cell (I,J> (from Stage S). Then we assume
J-V = tow V'V ^"o37* (1(M8)
where h is the depth of Layer 0 in cell (I,J) at hour t and A is given by
(10-14). (This is Eq. 8-7 in Part 1.)
Input and output information for this stage are summarized in Table 10-1,
and a schematic illustration of the relationship among the stages and I/O is
provided in Figure 10-2.
214
-------
Table 10-1. Summary of input and output parameters for
Processor P10 and their sources.
Input
Variable
Description
Source Stage
Output
Variable
Description
J,t ) Obukhov length
m (m) in cell (I,J)
at hour t
m
T (t ) virtual temperature
weather station n
at time t .
, Layer 1 averaged
east-west wind
(m/sec) in cell
(I,J), hour tffl.
, same as t except
1 north-south
component
Tt^m) temperature (°C)
K m of exhaust gas
of major point
source k at hour
V
D. diameter (m) of
K stack k.
w. (t ) exhaust velocity
(m/sec) of stack
k.
P4
P3
DELH
w
Pll
Pll
RAW
RAW
RAW
effective source heighl
(m) of k-th major poinl
source at hour t.
.
[I(k),J(k)] grid cell coordinates RAW
(row J, column I) of
major point source k.
•sk
stack height (m) of
source k
RAW
E,(t ;a) emission rate (moles/ RAW
hr) of species a at
time t from major
point source k
S(I,J,tm;a)
emission rate (ppm
m/sec) of species a
at hour t from all
sources in Layer 0,
cell (I,J).
215
-------
Table 10-1. (Continued)
Input
Variable
E(I,J,tm;a)
Description
emission rate
(moles/hr) of
Source Stage
RAW S
(Cont.)
Output
Variable
Sl(I'J'tm;a)
Description
emission rate (ppm
sec"1) of species <
species a from all
except major point
sources in grid cell
(I,J,) at hour t
m
o-jQ(I,J,tm) fraction (non-
m
dimensional) of
model surface H
penetrated by
terrain.
t ) same as a-™ except
applies to surface
z.(I,J) average terrain
elevation (m.MSL)
in cell (I,J).
h (I,J,t ) thickness (m) of
Layer 0 at hour
tm in cell (I,J)
h-i(I,J,t ) thickness (m) of
Layer I at time
tm, cell (I,J).
h?(I,J,t ) same as hn except
alies to Laer 2.
(t )
m
effective source
height of major
point source k
•2
mean air density
(moles m-3) in
Layer 2 at hour
t , cell (I,J)
same as «
except applies to
Layer 1.
P7
P7
P7
P7
P7
P8
DELH
P9
P9
S2(I,J,tm;cO
at hour t from all
sources in Layer 1,
cell (I,J).
emission rate (ppm
sec~l) of species a
at hour t from all
sources in Layer 2,
cell (1,0).
total number of major
point sources in Laye
0, cell (I,J) at hour
V
216
-------
Table 10-1. (Completed)
Input
Variable
Description
Source Stage
Output
Variable
Description
same as 2
except applies to
Layer 0.
P9
L(I,J) total length (m)
of major line
sources in cell
d,J)
N (I,J) number of minor
p point sources in
cell (I,J).
N .(I,J,t ) total number of
J major point sources
in Layer 0, cell
(I,J) at hour t
h (I,J,tm) depth (m) of Layer
0 m 0 in cell (I,J)
at time t .
RAW ZTA
CCI.J.V
fraction (0<£<1) of
Layer 0 in cell
(I,J.) filled by line
and point source plume
at hour t
m
RAW
Stage S
P7
217
-------
YYYY
OD
218
-------
SECTION 10
PROCESSOR Pll
SUMMARY
Processor Pll generates families of vertically integrated horizontal
winds for each of the model's layers. During daytime hours, namely times when
the surface heat flux is positive, this processor produces wind fields for
each of the model's three layers. However, at night when surface inversions
exist, it provides winds for Layers 2 and 3 only. In this case the flow field
in Layer 1 is generated in Processor P7, and the results are passed into
Processor Pll for amalgamation with the wind fields derived here.
219
-------
INTRODUCTION
Previous air pollution models have represented the wind fields either by
continuous functions fit to discrete meteorological data, or by flow fields
derived from mesoscale meteorological models. In this study we use neither
of these approaches. Following the consideration presented in Chapters 6 and
7 of Part 1 of this report and in Lamb (1983a), we use wind observations and
physical principles jointly to define a set or family of flow fields each
member of which is a possible description of the flow that: existed during the
time that the observations were made. In the aforementioned papers it is
argued that in the case of the atmosphere, observations and physical laws
together are inadequate to specify flow patterns to within better than a set
of functions. The ability of so-called objective analysis schemes to produce
a single functional description of a given discrete set of data stems from
the imposition of additional constraints that are not founded in physical
laws and which therefore lack universal validity. In the more sophisticated
of these schemes, use is made of empirical data such as spatial auto-correlation
functions of the flow in the given region. Our position is that this type
of empirical information is of great value, but its proper role is in estimating
the probabilities of the members of the set of possible flows that specific
observations and physical principles together define. These members of the
wind field set that have finite probability constitute an ensemble of wind
fields. For a given distribution of pollutant sources, there is a one-to-one
correspondence, through the equation of species mass conservation and chemical
reaction, between each member of the flow ensemble and each member of a
concentration ensemble. In other words, having defined the ensemble of flow
fields one can generate the ensemble of concentration fields associated with a
220
-------
given emissions distribution by "driving" the air pollution model with each
member of the flow ensemble and recording the outcomes. Ensemble statistics
of the concentration can then be derived by performing averages of the desired
quantity over the ensemble set, weighting each member of the set by the
probability of the flow field from which that particular concentration
distribution was derived.
This approach to pollution modeling is much more costly than the
conventional method of formulating differential equations that yield the
desired statistics directly. However, our approach has at least two
advantages that in our judgement compensate several fold for the added cost.
The first is that by deriving ensemble properties from a subset of the
ensemble itself, we avoid major sources of error associated with the
assumptions that one must invoke to formulate a set of differential equations
that describe a particular statistical property. In the case of modeling the
long range transport of photochemical pollutants, the nonlinear character and
interaction of the processes involved are so complex that it is very unlikely
that a single set of equations exist that would yield accurate estimates of
even the simplest statistical properties under all conditions.
A second advantage of our approach is that it provides direct access to
all of the statistical properties of concentration, such as the mean,
variance, frequency distribution, spectrum, etc., whereas the conventional
method gives only those limited properties, usually only the mean, for which
equations have been hypothesized. Following we describe the basic
mathematical steps implemented in this processor for deriving the ensemble of
flow fields in each of the model's layers.
221
-------
AVERAGE HORIZONTAL WINDS IN A LAYER BOUNDED BY TWO ARBITRARY PRESSURE SURFACES
In this section we develop a general method for deriving vertically
averaged winds in a given layer. Later we will apply it to each of the
model's three layers to obtain the necessary horizontal flow fields.
As in Part I we define the cell averaged value between two pressure
surfaces p and pfi (p > pR, i.e., surface a has the lower elevation)
u pup
Pp(x',y')
> = 1 ff(t(xjy'p)dpdx'4y' (11-1.)
P Vcrp J J J
A Pa(x',y')
where A denotes the rectangular domain x-|^ pfi, both V „ and the integrals in (11-la) are negative
(assuming 0).
Suppose that observations of v are available at N arbitrary sites on
surface pa, i.e.,
222
-------
,yn,P (xn,yn,t),t) (11-2)
n=l,2,...N.
where (x ,y ) denotes the location of observation n. (Throughout this
section, we shall use the caret O to signify observed values of a
parameter.) Suppose further that we have measurements of v and other
meteorological parameters along vertical lines between surfaces p and pg at a
total of M locations, i.e., we have
These surface observations and sounding sites are illustrated in Figure 11-1.
In general the number N of surface stations must equal or exceed the number of
soundings, i.e.,
N > M. (11-4)
Later we would like the option of supplementing the rather sparse upper
air data with estimates of the flow aloft extrapolated from the numerous
surface measuring stations. Toward this end we define the function
9a(x,y,p,t) to be the ratio of the wind velocity at (x,y,p,t) to the value at
(x,y,p ,t) where a is any pressure level. Specifically, we have gfl =
(gua,gva) where
n fv v n
gua(x,y,p,
n fv u n tl - V(X»y.P,t)
gvaU,y,p,tj = V(x,y,p' t).
223
-------
VERTICAL SOUNDINGS
AT M SITES
SURFACE OBSERVATIONS1
AT N SITES
Figure
Aoo B^dof the Horiz
Available on the Bottom Surface °' ] of Lajer p °0 "" dls°
Assunlng that g and v vary by only sma11 fractions of twelves as (x y)
ranges over the area A of any grid cell, we get from (11-1) and (11-5)
or
where
0 s
Pa -
ga(x,y,p,t)dp
p - vcr(xiy,t)(
= v(x,yspa,t),
(11-6)
(11-7)
(11-8)
224
-------
and ~ is defined as in (11-1).
At the M sites of the measurements of the vertical wind profile, we can
evaluate Q. We have
~ P
« g at each point (x,y) in space by
interpolation of the "observed" values _. That is, we assume
p
I Wm(r)
m=l m ~
_2
where W (r)is a weighting function, e.g., r , to be selected later, and r is
Ml ^ *** r^
the vector separation of (x,y) and the observation site (xm.ym)- Using
(11-10) and the N (>M) observations v of the wind on surface p , we can ap-
— ""Of 0(
proximate the layer average winds at all those points (x ,y ), n^m, where wind
observations are available only on the surface p . In other words, we use the
horizontally interpolated 0 function to approximate layer averaged winds at
- p
points where only surface wind observations are available. Thus, we define
225
-------
where
6 is given by (11-10), and the vector product on the right side of (11-11)
is defined by
5 (axbx>ayy
where a = (ax,ay) and b = (bx,by).
We will work later with the stream function and velocity potential which
are functions of the vorticity and divergence, respectively, of the wind
field. The vertical component of vorticity £ is defined by
s k-Vxv =
i a 5
8x 3y 8p
u v w
3u
3y »
(11-12)
and the horizontal divergence 6 is
3u . 8v .
3x 3y
(11-13)
Since we are concerned with layer averaged winds , we must determine
how and are related and how <6> and are related. From (11-1) and
(11-12)
.
3y
(11-14)
226
-------
From Part I, page 21, we have
3n An"
53T - v T# ) + •£ in (Vap) (H-15a)
i n .. (n-i5b)
where A is the horizontal area on which < > is defined; V Q is defined by
ap
(11-lb);
and
s -i-J J —A
dx'dy'-
227
A
with a similar definition of the surface average
Assumption: On pfl, y(x,y,pa) - and on Pp,
v(x,y,pp) - . (11-17)
-------
Note that
£- ln(V ) = 1H _§_ dpdx'dy'
8x |IUV V 8x
A
3||r f i i\ / i i 11 -i i .1 i
"^VR^ lUPaWy'J-p^xly')] dx'dy'
ap
A
where the domain A is a function of x inasmuch as A = x-Ax/2 ~ a
3x 8x
Thus, from this result and (11-14) we conclude that
(given (11-17)) (11-21)
and similarly
' (g1ven
228
-------
Using Helmholz theorem that any vector can be written as the sum of a
rotational and a divergent component we can write
\f — V +
from which we have
(11-23)
= • (11-24)
Since by definition v^ is nondivergent, we can express it in terms of the
stream function ¥:
= JS x S* (11-25)
and since v is irrotational we can express it in terms of the velocity
A
potential x:
= Vx-
~X ~
(11-26)
From (11-25) we get
< V =
and from (11-26)
! j 5
001
= |X i + |X
~X 9x ~ 3y
Combining (11-24), (11-27), and (11-28) we get
(11-27)
(11-28)
(ll-29a)
229
-------
= §3! + §* + v (ll-29b)
8x 3y
where u and v are the average values of and over the entire model
region and are known functions of time only.
Differentiating these expressions in the manner of (11-21), (11-22), we
get
(ll-30a)
= V2x- (ll-30b)
We seek the functional forms of ¥ and x that are consistent both with physical
laws and with our observations OB,n=l...N (see Equation 11-11).
(We assume that the wind observations are error free.)
Our observation set might represent a time period in the past for which
we would like to determine the origin of pollutants responsible for an
episode; or it might represent so-called worst-case meteorological conditions
for which we wish to test a control strategy; or so on.
We will consider only two of the physical laws as constraints on the flow
fields that we consider, the momentum law and the mass conservation principle.
In pressure coordinates these are expressed, respectively, by
9~ + (v-Vu)v + u>8~ + fk x v = ~ V
-------
where u> is the "vertical velocity" in pressure coordinates, i.e.,
= If + 2* V + w If (11-33)
* is the geopotential height, f the Coriolis parameter, K the "eddy
viscosity", FM represents the viscous source of negative momentum, and v =
(u,v) as before. Expanding (11-31) into its two components and using (11-13)
we get
2 2
9u * ..9u j. w9u * 9u f< - 3* . ,/r3 u . 3 un . c
9t + ^x + V3y + % ' fV - " 3* + KC2 + + FHx
6 = - 9ui/3p (11-36)
It is useful to convert (11-34) and (11-35) into the vorticity equation,
partly as a means of eliminating explicit dependence of the velocity on *.
Guided by (11-12) we differentiate (11-35) with respect to x and (11-34) with
respect to y and then subtract the two equations to get
r3uj 3v 3uj 8u, ,, r32C . 32£n A 9FHy - 8FHx,
L3x 3p 3y 3pJ " l^2 9y2J ' 3x 3y
(11-37)
231
-------
The equation governing <6> is found by averaging (11-36) in the manner
of (11-1):
<5> = - gdpdx'dy' (11-38)
K
> (x;y') - u>s(x;y')]dx'dy' (11-39)
where u>a denotes u)(x;y',pa). Using (11-16) we can write (11-39) in the final
form
<6> = § 0*- J#]. (11-40)
Later when we apply (11-40) we will describe approximations of io*, etc.
Turning next to the vorticity equation (11-37) and the derivation of the
equation governing , we see that under assumption (11-17) we can write
^, ^ - 0. (under assumption 11-17) (11-41)
op dp
This permits the term in brackets in (11-37) to be neglected. (This term
represents the generation of £ by the rotation of horizontal vorticity.) We
3f /
will also omit the "beta" term 9y in (11-37). With these approximations we
can write (11-37) in the simplier form
232
-------
= K _
O \i 5iv/
\J *\ O y
(11-42)
3y ^
where we have used the approximation (see Haltiner 1971, page 151)
FH = - g (11-43)
and i is the shear stress.
Averaging (11-42) in the manner of (11-1) we generate the following
terms, expressed following the analyses of Part I:
§1 *> + <> si ln(V
(11-44)
3 3
3x 3x (Kp
(11-45)
. ~ap
v
ap
V
233
-------
< =
32t
The terms like < — *> are quite complex and we will approximate them simply by
2 2
> - -2<£> (11-48)
9x 3x
Combining all the above and assuming <£6> = <£><6>, we get
= - f<6> + u
Vap t dt ~H ~ (11-49)
aB
where
3^
+ KC
(11-50)
(11-51)
Pa = P^
and
al 5 si
In deriving (11-49) we assumed that i_ = 0 on pQ in anticipation of the
~z p
eventual use of this equation in applications where the upper surface pg is
above the friction layer. We will use Haltiner's approximation (page 152) for
i , namely
234
-------
where pn is air density in layer ot.p, in particular the mixed layer; CD is the
drag coefficient and is the layer averaged wind vector given by (11-51).
It is convenient to group the terms in (11-49) according to whether they
represent sources or sinks of <£>. The terms involving £' do not fall into
either category, but for now we will combine them with the source terms and
write (11-49) in the form
af
where
af
3x
3y
(11-55)
3
H s - f<6> +
"
06
3^
(11-57)
(11-58)
We are now ready to collect the equations that we will use to generate
the ¥ and x fields. First, let the superscript denote the time variable,
i.e.,
235
-------
+ —<£> + <£> G
3x 9y »
(11-60)
a2 a2 11
f_ + d i <•/->" - ur
L—o^ —oJ ^s,^ ~ n
Let 4 be the complex Fourier transform of <£> , i.e.,
00
i i
<>J = -
dkdk
471
-co
oo
f (* -ik*x
J J <^(x)>Je ~ ~ dxdy
-03
where k = (k ,k) and x = (x,y). We will represent transform pairs by
4J ~ °, (11-62)
and will make use of the convolution theorem
00
dl-63)
(2n)
-00
where x and k are 2-D and
f(x) «. F(k).
236
-------
Note from (ll-61a) that
00
i rr i ^'X
i rr,. ,j,,.^ dk dk (11_64)
-co
7 1 I *Y* vt'c dMkv> (11-65)
4^ ii * x y
-00
and in general
~<£>J ^ dk ) 4 (11-66)
3xn x
Now, replacing each term in (11-60) by its Fourier transform and making use of
(11-63) and (11-66) we get
09
-oo
oo
-CD
(11-67)
00
J(k')GJ(k-k')dk'dk' = 4/tV
^ XA^ ' xys/ ** ' XV
-oo
237
-------
where
GJ *-» GJ (ll-68a)
HJ ~ HJ (ll-68b)
UJ «-» J (ll-69a)
l/J «-» J (ll-69b)
We want equations for 4* and x- Thus, we obtain from (ll-30a) and (11-66)
(11-70)
where
H
-------
The transforms U and V of the velocity components and can now be
expressed in the form (see 11-29)
UJ = -ik AJ + ikxBJ + uJ6(k) (ll-73a)
where
= ikAJ + ik S° + vJ6(k) (ll-73b)
x y
X°(x) «-> BJ(k) (11-73C)
and 6(k) is the delta function of the wave number vector k. Then from
(ll-30b), (11-40), and (11-66) we find that 8J must satisfy
where
kx+ky
<6> = - - [^"(x.JAt) - ui (x.JAt)] (11-75)
~
and ID is considered to be a known variable.
Finally, the observations QB given by (11-11) must be satisfied.
Hence, from (11-29), (ll-61a), (11-66), and (11-11) we have
00
-00
n=l,2...N;
239
-------
00
—2 CkxAJ(k) + kvBJ(jc)]e~ n dk = „ - v(JAt)
(2n) JJ y
y (11-77)
-OB
n=l,2...N
where xn = (xn,y ) are the observation sites described above, and u and v are
the averaged observed winds in the entire model domain at time t = JAt.
Equations (11-72, 74, 76 and 77) comprise a set of 2N + 2 relationships
that the transforms A and B, or equivalently and , must satisfy. A
fundamental difference between this system of equations and the system
currently used in meteorological modeling is the presence of the 2N equations
(11-76, 77) associated with the wind observations made at time JAt. In
conventional meteorological models, the observed data are transformed through
ad hoc means into a description of an initial velocity field; and using it,
the two equations (11-72 and 74) are solved as an initial value problem. As
we shall discuss below, this practice is not supportable on scientific
grounds, and when it is employed, the velocity field that one obtains is
largely an artifact of the ad hoc procedure that was used in the formulation
of the initial state.
The fact is that the N observations of the wind velocity at any instant
JAt, say, do not uniquely determine the velocity distribution at that moment.
Rather, these measurements, specifically (11-76, 77), together with the
continuity equation (11-74) define a hyperpiane in the phase space of wind
velocity. Each point on this plane represents an entire vector function
[u(x,JAt), v(x,JAt)], xeD. We pointed out in Section 6 of Part 1 that each of
240
-------
these functions in a possible description of the actual velocity field that
existed at the moment t=JAt that the observations contained in Eqs. (11-76)
and (11-77) were made. If we had to choose from this set of possible velocity
functions the one that we thought described the flow that actually existed at
the moment the observations were made, our choice would be guided by whatever
previous wind data we had seen for the region in question and on our knowledge
of the nature of atmospheric motion generally. On this basis we would declare
many of the possible functions to be "unlikely" descriptions because they have
characteristics that are incompatible with our knowledge. For example, even
though all functions in the set satisfy the observations and the continuity
equation (11-74) exactly, some contain hurricane force winds between the
observation stations; some contain intense vortices and jets; and so on.
However, even after all these are dismissed, there remains a virtually
infinite set of functions no one of which can be ruled out as a plausible
description of the actual flow.
This suggests that we should assign to each point on the hyperplane of
possible flow descriptions a weight p, say, whose magnitude is a measure of
our conviction that that particular function is the one that describes the
flow field that existed when the observations were made (at time JAt). To
those bizarre functions mentioned above, we would assign the weight p=0; and
to the remainder of the points on the hyperplane of possible descriptions we
would assign weights l>p>0. Without going into details, we point out that the
weight p is synonymous with the "probability" of occurrence of the flow field
to which the p is assigned (see Lamb, 1983c).
241
-------
The problem is to formulate a quantitative rule for assigning the
weights. Clearly, we cannot examine every function on the hyperplane of
possible flows and assign it a weight subjectively. We need a mathematical
procedure.
In Lamb (1983c) several approaches to this problem are proposed. It is
shown in that paper that in effect the method that has been employed to date
in both diffusion and meteorological modeling studies has been to assign a
zero value of p to all of the possible flows except one and to assign the
value p=l to that single function. This is equivalent to declaring that we
know the correct description of the flow field at t=JAt without any doubt.
This "correct" description is obtained using one of the so-called objective
analysis formulas. But many such formulas exist — no formula of this type
can be universally valid ~ and therefore it is illogical to label any one
description of the flow the "correct" description. The point is that there is
no scientific basis whatever for assigning unit weight to only one function in
the set of possible flows and zero weights to all others.
One rational method of assigning p is the following. It has been
theorized that in 2-D fluids kinetic energy is partitioned among the spatial
_o
fluctuations in the flow in proportion to | k| , where k is the wave number
vector of the fluctuation. To a large extent this theory is supported by
measurements of kinetic energy in the free atmosphere. Thus, let Q represent
the manifold in velocity phase space of functions whose Fourier transform is
_ ^
consistent with an energy spectrum of the form E(k) ~ | k | . And let Q
represent the manifold formed by the intersection of Q and the hyperplane of
possible flow descriptions. We now advance the hypothesis:
242
-------
q(v), if v sQ;
p(v(x,JAt),xeD) = (ll-78a)
~ 0 , otherwise
where l>q(v)>0 is a function of the entire velocity field v(x,JAt), xeD that
we have yet to specify. Hypothesis (ll-78a) is simply a statement of our a
priori belief that the description of the true velocity field at time JAt is a
member of the subset ft of the hyperplane of possible flows. Observations of
the flow field at times following JAt and observations of the rate of
dispersion of material in the flow might force us a posteriori to reject this
hypothesis.
To supplement (ll-78a) we might advance the following hypothesis:
q(v) ~ E(v)"1 (ll-78b)
where E is the kinetic energy of the flow field v integrated over the entire
model domain D. In practice one can treat only a finite number of the
velocity fields contained in Q. For example, in accordance with (ll-78b) we
might select, say, the 20 members of Q that contain the least kinetic energy,
namely the 20 points in fi closest to the origin of the phase space; and assign
p values to each using (ll-78b) and the constraint that the sum of the 20 p's
be unity.
In hypothesis (ll-78a,b), and in the conventional approach discussed
earlier, consideration is given only to the observations made at the single
time JAt. However, observations made after JAt provide valuable information
that can be used by way of the momentum equation (11-72) to obtain a better
approximation of the probabilities p of the flow fields at JAt than we can
243
-------
obtain from empirical-theoretical considerations alone. To implement this
approach requires rather complicated dynamic programming procedures. We plan
to utilize this method in the "second generation" regional model.
For the present we shall replace Eq. (11-72) with the much simpler, and
much weaker, diagnostic constraint
dx = C(t) (11-79)
D
where C(t) is the circulation around the perimeter of the modeling region D.
This equation is a statement of Gauss' theorem. We will derive estimates of
C(t)from observed winds and treat it as a known function of time.
We see from (ll-30a) that Eq. 11-79 is a constraint on the stream
function, and hence on A (See 11-71). Substituting (ll-30a) into (11-79) and
making use of (11-71) and Stoke1s theorem, we can express (11-79) in terms of
the transform of the stream function:
f f ik L ik L k2 k2
-M A(k)[(e x x-l)(e y y-D( XkV )ldk dk = c(t) (11-80)
4^JJ KxKy x y
Thus, the initial version of processor Pll will be based on the equation set
=>AJ(k),BJ(k)
244
-------
Each solution set [A(k), o(k)] derived from (11-81) yields layer averaged
winds [, ] through the following equations, whose origin
is (11-29):
oo
ik-x
= _i-^\[-kyAJ(k) + kxBJ(jk)]e "dk + u(JAt) (11.82a)
-08
00
= - [k A(k) + k vB(k)]e~ ~dk + v(JAt)
(2n)z x
-oo
In the first generation model we will derive between ten and twenty flow
fields from (11-81, 82) for each hour of the model simulation period and
assign weights p to each using hypothesis (11-78) above. Let us call this
finite ensemble of flow fields at hour t=JAt W . (We should add that even
though the members of W are not explicitly related to those of W , K^J, due
to our replacing the momentum equation (11-72) with the simple diagnostic
expression (11-79), there is implicit coupling of these ensembles by virtue of
the fact that each of them is defined in terms of the actual winds observed at
each hour.) We generate a corresponding finite ensemble C of concentrations
(containing M=10 to 20 members) associated with a given emissions distribution
S(x,t), xsD, 0 < t< T, by "driving" the dispersion model with M wind fields
[, ] , xeD, J=0,l,. . . JMAV, u=l,...M. Each of the M flow
~ ~ (J ~ nnft
fields is created by selecting a and a at random from
the ensembles ^ for each hour J=1,...JMAX in the period of interest. Either
ensemble mean values of the concentration or the frequency distribution of
concentration or any other statistical quantity of interest can be computed
245
-------
directly from the ensemble C for any desired receptor sites and times.
Details will be presented in a future report. (See also Lamb 1983 a,b,c).
Below we outline the basic structure of Processor Pll which utilizes the
techniques described above to compute the layer averaged winds . for each
of the regional model's three layers. Winds are determined for Layers 2 and 3
at all hours and in Layer 1 during the day only. The nighttime air flow in
Layer 1 is simulated in Processor P7 and passed into Processor Pll in the form
(VL,VL) for amalgamation with the flow fields computed here.
246
-------
STAGE UV11
This stage computes layer averages of the measured winds aloft that are
used as the observations in solving Eq. (11-76) and (11-77). Consider first
the calculations required for Layer 3.
(1) At each of the M upper air weather stations compute the layer averaged
horizontal wind components as follows:
P2(xm,t) ~~} m=l,2,...M
mode 0 and
(11-83)
mode 1
^pl
3 2
where u and v are the measured wind profiles. Since the input winds (u, v)
are in z coordinates, it will be necessary to transform them to p coordinates,
e.g., u(x ,p,t), using the pressure-height function p (z,t).
~ro m
(2) Since the bottom surface of Layer 2 is p in mode 1 and p^ in mode 0,
the corresponding expression for the Layer 2 observed winds is
247
-------
pl J
PB(xm,t)
8 m
(11-84)
p0(x ,t) V mode 0 and
t. ~m f
v(xm,P,t)dp
m
where
~Pvs(x,t) , mode 1;
PB(x,t) =| (11-85)
P1(x,t) , mode 0.
(Recall that mode 1 applies when Ap,=l, and mode 0 applies when Ap-,=0.)
(3) Processor Pll is not used in Layer 1 during mode 1 conditions. In that
case the flow field in Layer 1 is generated in P7. During mode 0
conditions, the depth of Layer 1 is set in P7 to a value that is
approximately the top of the shear layer. Assuming, then, that the flow
speed and direction in Layer 1 are roughly uniform we adopt the following
expressions for OB -,:
OB , = u(x t) -v
n OB»1 n n=l,2,...N
OB,l = KZn'
mode 0 only
(11-86)
where x , n=l,...N are the sites of the N surface weather stations.
~n
(4) Optional. In order to supplement the sparse upper air data on which the
observations in Layers 2 and 3 are based, we define the following "shear
functions" for use in estimating the flow in Layer 2 over surface wind
stations.
248
-------
2 u(x p (x_,t;,t;
~m uc~m ' BFI....M
(11-87)
mode 0 only
In mode 0, pvg coincides with the ground surface and hence the function
is the ratio of the measured vertically averaged wind in layer 2 to the
observed wind at ground level.
Next we interpolate values of 2 and 2 at each of the N surface
stations from the values given by (11-87) for the rawin station sites as
follows:
M
' ""- 2 — (11-88)
where
J = C(x-o2 + (y.-yJ2]"1- (ii-89)
A similar operation yields o- We can now define pseudo Layer 2 winds
over each surface station as follows
= 20(xn,t) ^^
mode'6'only. (11-90)
,t) optional
This optional method of supplementing the upper air data can be implemented as
desired.
249
-------
Stage PHLOB
This portion of Processor Pll solves the equation set (11-81) to obtain
members of the flow field ensemble. When a surface inversion is not present
(mode 0), this stage computes averaged winds for all three layers of the
model. Otherwise (Mode 1) simulations are performed only for layers 2 and 3,
and the corresponding Layer 1 flow is assigned the values
= v, 1
1 VL I mode 1
1 = VL J
where v^ and VL are the wind components in the cold inversion layer
generated by Processor P7 at each grid cell and at each hour that the
inversion (mode 1) exists.
The input requirements of Processor Pll are summarized in Table 11-1 and
its outputs are listed in Table 11-2. Processor Pll and its interfaces with
the processor network are illustrated in Figure 11-2.
250
-------
Table 11-1 Input Requirements of Processor Pll and
Their Sources
Variable Description Source
p^(x,t) top surface of Layer 3 in pressure P8
(mb) coordinates
P2(x,t) same as above, except top of Layer 2 P8
PycQ^t) pressure (mb) at the virtual P7
surface
Pi(x,t) top surface of Layer 1 in pressure P7
(mb) coordinates
u(>e,z,t) observed east-west wind speed component PI
(m/sec) at rawin station m(=l,...M)
at hour t at elevation z(m MSL)
v(x ,z,t) same as above, except north-south PI
wind component
u(x ,t) observed east-west wind component P3
(m/sec) at surface station
n(=l,...N) at hour t
v(x ,t) same as above, except north-south P3
wind component
p (z,t) pressure (mb) at elevation z PI
m
(m MSL) over rawin station
m(=l,...M) at hour t.
<6(x,t)>3 average-,horizontal wind divergence P8
(sec ) in Layer 3 in grid cell
centered at x at hour t.
<6(x,t)>2- Same as <6>3 except applies to Layer 2 P8
<6(x,t)>-. Same as above except applies to Layer 1 P8
y. Average east-west wind component in P7
Layer 1 during mode 0 (generally night-
time) hours (m/sec)
y. Same as above except north-south wind P7
component
251
-------
Table 11-2 Outputs of Processor Pll
Variable Description
3 Vertically averaged wind (m/sec) in Layer 3 at
grid cell centered at x at hour t (east-west
component)
3 same as above, except north-south wind component
2 Vertically averaged wind (m/sec) in Layer 2 at grid
cell centered at x at hour t (east-west component)
2 same as above, except north-south wind component
. Layer 1 averaged wind (m/sec) in grid cell centered
at x at hour t (east-west component)
, same as above, except north-south wind component
252
-------
T
YYYYYY
(B
< *?
? 0)
w
9
i 8 c
^ O (D
c 2 w
o a £
*•• ® fe
S? JE o
•fc ** w
w j= c
3 ** S
0)
3
O)
V v
253
-------
SECTION 11
PROCESSOR P12
DEVELOPMENT
Processor P12 calculates parameter fields required in the description of
interfacial material fluxes between Layer 0 and Layer 1 and between Layer 1
and Layer 2. It also provides estimates of the horizontal eddy diffusivities
in all layers. The mathematical expressions for the parameters we require
were presented in Part 1 of this report. Below we provide detailed
descriptions of the implementation of these expressions and further commentary
on their meanings and origin. The discussion is divided by stages just as the
calculations within P12 are divided.
Stage K
This stage computes estimates of the horizontal eddy diffusivity
K (I,J,t ) at each grid point (I,J) at each hour t , for each of the model's
three layers n = 1, 2, and 3. Here eddy diffusivity refers to the action of
the small scale wind fluctuations generally associated with turbulence. These
include wind fluctuations generated by wind shear at the earth's surface and
across the top of the mixed layer, and lateral fluctuations induced by
convective thermals, principally at the top of the mixed layer and at the
ground. (The effects of small scale vertical fluctuations are handled by
parameters like ft'-, that we consider later.)
254
-------
It is important to keep in mind the role that the eddy diffusivity will
play in the regional model. In conventional studies, the magnitude of the
horizontal diffusivity controls the rate at which a plume expands about its
"centerline" and this in turn affects the magnitude of the "mean"
concentration that the model predicts. In the present model the eddy
diffusivity plays a lesser role. As we discussed in Part 1, Chapter 7, and in
the description of Processor Pll in this report, the regional model will be
applied to simulation of the ensemb1e averaged concentration. This requires
the model to be run a number of times, each time using a different flow field
[n, n] from the ensemble of flow fields described in Processor Pll
following (11-77) but using the same eddy diffusivity K in each case. The
results of all the individual simulations are eventually superposed to arrive
at the ensemble averaged concentration properties. In the case of a point
source plume, the superposition of the individual realizations within the
ensemble might appear as shown in Figure 12-1. Note that K controls the rate
of expansion of the plume in each realization, but the spatial variability of
the wind fields within the flow ensemble are what govern the envelope of the
superposed results and hence it is these larger scale variations in the flow
that dominate the ensemble mean concentration. In this case it can be shown
that K affects the mean square concentration. This in turn affects the rates
of second-order chemical reactions and also the expected deviation of the mean
concentration the model predicts from the actual concentration values that one
might measure.
Another way of viewing the role of the eddy diffusivity in the regional
model is to consider that between the largest turbulent fluctuations
represented by the diffusivity K and the smallest perturbations in the
255
-------
subsynoptic flow that our network of wind sensors can resolve, there lies a
range of "mesoscale" wind fluctuations that are not accounted for either in
the K's or in the "transport flow". And it is these fluctuations that give
the ensemble averaged plume the width I illustrated in Figure 12-1.
ENVELOPE OF ENSEMBLE
AVERAGE PLUME (2)
CENTERLINE OF ENSEMBLE PLUME WIDTH (o).
AVERAGE PLUME _ CONTROLLED BY
s
CENTERLINE OF PLUME
IN ONE REALIZATION
Figure 12-1. Illustration of the superposition of five hypothetical
realizations of an ensemble of point source plumes. The
width I of the ensemble of plumes is controlled by the
character of the flow field ensemble. The width a of the
plumes in the ensemble is controlled by the turbulent, eddy
diffusivity K.
On the basis of the points raised above, we expect that the values
assigned to Kn in the regional model will directly affect the mean square
concentration values more than the mean. However, they may have a significant
indirect effect on the mean concentration through the influence that they will
have on photochemical reaction rates. We can only determine through test
256
-------
calcuations just how large these effects can be. For now we will postulate
simple forms for the lateral diffusivities K that we can use in the first
generation model.
Lacking a firm basis on which to estimate the eddy diffusivity, we will
use heuristic arguments to formulate expressions for K . Our first assumption
is that the lateral eddy diffusivity is effectively zero during nighttime
hours when convection is absent and/or stable stratification damps shear
generated velocity fluctuations. Our second assumption is that during the
day, convection is the primary source of horizontal wind fluctuations. Let D
be the horizontal scale of convective cells, or rolls, and let H be the depth
of air in which the convection is confined (usually by a stable layer an
elevation H above the ground). We will base our parameterization on the
concept that convective circulations have a toroidal shape with rising motion
along the axis of the toroid, outward horizontal motion along the top, sinking
motion along the outer surface of the toroid, and horizontal motion at its
bottom directed in toward the axis. In this picture the outer diameter of the
toroid is D and its height is H. If the circulation were indeed this simple,
then a particle would receive a horizontal displacement of order D in the time
At required for the particle to traverse the depth H of the circulation.
We imagine that each convection toroid is surrounded by others and that a
particle in any one toroid can escape into an adjacent one during the time the
particle is moving downward along the outer edges of the toroid or during
times when the convective cell is dying. If these migrations of particles
from one cell to another occur at random, then as a rough approximation we can
treat them as a classical random walk process. In this case the effective
257
-------
diffusivity K is
K ~ D2/At. (12-1)
(Incidentally, the random walk and the gradient transfer, or K theory, models
of diffusion are essentially equivalent.)
It has been found in empirical studies of dry convection in the planetary
boundary layer that D ~ 1.5H. If we represent the transit time At by
At = H/W, (12-2)
where W is a characteristic velocity scale of the convective motion, then we
get from (12-1) and (12-2)
K ~ WH (12-3)
which we could have guessed at the outset on purely dimensional grounds. In
dry convection W = w*, where
w* = (HQg/e)1/3 (12-4)
where Q is the kinematic heat flux at the ground and 6 is the mean temperature
in the mixed layer. When cumulus clouds are present (moist convection), the
appropriate measures of W and H are unknown. In this initial work we shall
assume that when convective clouds are present
W = wc
where w is the average upward air speed within clouds and
H = H2 + H3 (12-5)
258
-------
where H^ is the depth of the convective layer below cloud base and HO is the
depth of the layer occupied by clouds. These estimates will be tested in
future studies.
Based on the considerations presented above, we propose the expressions
summarized in Table 12-1 for the diffusivities K in each of the model's
layers.
Horizontal
Diffusivity
Ki(i»J>V
K2(i,j,tm)
K3(i,j,tm)
uu,y
>0 <0 and <0 and
ac = 0 ac * 0
0 .lw*H2 .lw*H2
0 .lw*H2 .lwc(H2+H3)
0 0 .lwc(H2+H3)
Table 12-1 Summary of the expressions used to estimate the horizontal
eddy diffusivity K in each of the model's three layers.
259
-------
The relationships between H2 and H3 and the thickness h of the model layers,
which are inputs to P12, are as follows:
H2 = hl * h2 (12-6a)
H3 = h3 . (12-6b)
In Table 12-1, a is the fraction of the sky covered by cumulus clouds in a
given grid cell at a given hour. The velocity wc(i,j,tm) required in the
expressions for K is an input parameter to processor P12; and w* should be
computed by
w* = (H2Qg/00)1/3 (12-7)
where hL is given by (12-6a); Q is the surface heat flux, an input parameter;
g = 9.8m sec ; and 0 is the surface temperature interpolated from the
surface virtual temperature measurements T as follows (compare with Eq.
7-108)
*r »— *—TT /i "\ _»
n=l "
where
rn = [(iAx-xn)2 + (jAy-yn)2]% (12-9)
and (x ,y ) are the coordinates of surface weather station n. The values of
w* obtained from (12-7) will be needed again in stage WW1 below.
The input and output requirements of Stage K are summarzied in Table
12-2.
260
-------
Stage WWO
Here we calculate parameters A.+ , w+, w_ and u associated with turbulence
phenomena in Layer 0 and on the surface H that separates Layers 0 and 1. The
definitions of these parameters were developed in Part 1.
First we compute the vertical velocity variance on surface H :
where
Of,
mj
2/3
4 =
(12-10)
(12-11)
-164)"4 , if 4 < 0;
, if 4 > 0
(12-12)
"0.4 [0.4 + 0.6exp(4|)], if | < 0;
_0.4 [1.0 + 3.44 - 0.254*], if 4 > 0.
(12-13)
and u* is the friction velocity. In Eqs. (12-10) - (12-12), the variables
a , u*, and | are all functions of (i,j,t ).
Next we approximate the time derivative of h by
-i
(12-14)
where At = 3600 sec is the time interval at which hQ and other variables in
the processor network are available.
261
-------
With these preliminary variables computed at each grid point and hour t ,
we can now calculate the following output parameters (see Eqs. 8-10, 8-12, and
8-13 of Part 1).
where erf is defind by
X ,2
erf(x) = -£ fe"X d\. (12-16)
-Jn J
n (12-17,
where a^Q) given by (12-10), and h , given by (12-14), are functions of
w.(i,j,tm) = n0(i,j,y + w+(1fj,tB) (12-18)
v(i,j,tm) = u*(i,j,tm) (12-19)
The parameter v is the entrainment velocity of ambient air into plumes from
surface emissions and is discussed in detail in Chapters 5 and 8 of Part 1.
At this point it is advantageous to determine the values of the cumulus
flux partition parameter ¥ defined in Part 1, Eq. 5-44. Interim values ¥' of
this parameter were developed in Processor P8, but we must test whether they
satisfy criterion (5-47) of Part 1, namely
262
-------
w+A+(l-aTO)
V < j v— (12-20)
v* c
where w+ and \+ are the parameters that we just computed in (12-17) and 12-15)
above, and where
w9-w
Vc = - -jr^ - fQ/A6 (12-21)
(see Eq. 5-45 of Part 1). All of the variables in (12-20, 21) are functions
of space and time and each of them, except w+ and \+, is an input to this
processor. Consult Table 12-2 for definitions.
To obtain the final values of 4* for output to the model, proceed as
follows. First, at each grid point (i,j) and at each hour t compute the
variable V(i,j,tm) defined by (12-21) using the input fields w2(i,j,t ),
W-O'.j.t ), etc. Next, compute
f (i j t ) - (12.22)
*(1J'V oc(l,J.t,)vc(i,J,tB) (1222)
where w+ and X+ are from (12-17) and (12-15) and VjQ and a are inputs to this
processor. Finally, we have
ni,j,tm) = min{«{"(i,j,tm)r'(i,j,tm)} (12-23)
where V is the interim estimate of f that is input from processor P8. The
values of Y should be recorded in the model input file MIF. This parameter is
also required in the next stage of this processor.
263
-------
Parameter summaries for Stage WWO are given in Table 12-2.
Stage WW1
This stage computes the parameters W w, , and w^, (defined in Part 1,
Chapter 4) that are associated with material fluxes across surface H,. To
compute these quantities we must first define several auxiliary parameters.
m
) = vertical velocity at level z, during neutral and
unstable conditions (input from P8)
0.6w*(i,j,y if L(i,j,tm) < 0;
(12-24)
0 , if L(1,j,tJ > 0
where (12-25)
where At = 3600 sec is the time interval at which h-, and hQ values are
available.
Let w be the solution of the following transcendental equation:
2 — —
W«T W - W.
(12-26)
264
-------
where erf(x) is defined by (12-16), V is from (12-23), etc. (Note that , cr , w , 0,., and hence w are all functions of (i,j) and tm). Now let
w(i,j,t)=w'(i,j,t) - O'J''
R1
mD1
m
At each grid point (i,j) and each hour t compute
(12-27)
erf(-
if
wl
iJ. , if L(i,j,tm) > 0;
(12-28)
Jwl
exp t
'wl
W
R1
, if L(i,j,tm)<0
,
wl
nvs(i,J,tm) , if L(i,j,
(12-29)
'WniOJ.tJ , if L(1,j,tJ < 0;
0.
(12-30)
265
-------
To evaluate the function erf in (12-28 and 12-29) in the limit as a ,
0, the following procedure should be used.
if a , = 0, erf (— — ) = SIGN(A)-1 (12-31)
This condition should be encountered only rarely since a , =; 0 when L > 0 (see
Eq. 12-24) and in this case W-, and W-,M have the forms given in (12-28) and
(12-29) that do not involve the error function.
The input and output requirements of Stage WW1 are summarized in Table
12-2. The processor and its interfaces with the processor network are
illustrated in Figure 12-2.
266
-------
Table 12-2 Input and output variables of each
stage of Processor P12.
Input
Variable
V'-J'V
h2(i,j,y
Description
thickness (m) of
Layer 1 in cell
(i,j) at hour tm
thickness (m) of
Layer 2
Source
P7
P8
Output
Stage Variable
K ICjd .J.V
K2(i,j,tm)
Description
horizontal eddy
diffusivity (m/sec)
in Layer 1, cell (i,
hour t
same as K, except
Layer 2
h-a(i»J»O
J m
thickness (m) of
Layer 3
P8
,_
m
same as K, except
Layer 3 i
fraction (0
-------
Table 12-2 (Continued)
Input
Variable
Description
Source Stage
Output
Variable
Description
Obukhov length (m)
in cell (i,j) at
hour t
P4
WWO
(cont.)
average speed (m/s
of fluid moving uf
relative to surfa(
thickness (m)
of Layer 0 in
cell (i,j) at
hour t
P7
same as w+ except
speed of downward
moving fluid
entrainment veloc
(m/sec) of plumes
in Layer 0
fraction (0
-------
Table 12.2 (Concluded)
Input
Variable
Description
Source
Stage
Output
Variable
Description
vertical velocity
(m/sec) on
surface H-, in
neutral and
unstable conditions
P8
WW1
w-i(i>j»t )
vertical velocity
parameter (m/sec)
on surface H,
vertical velocity
parameter (m/sec)
on surface H,
w*(i,j,t ) convective velocity Stage K
m
scale (m/sec)
vertical velocity
(m/sec) on surface
H, due to horizontal
divergence in the
flow field.
thickness (m) of P7
Layer 1 at (U,tffl)
thickness (m) of P7
Layer 0
L(iJ,t )
«Ki»j»t )
(i,j,t )
Obukhov length
scale (m) in cell
(i,j), hour t
ffl
cumulus flux
partition function
(dimensionless)
fraction (0<0 <
of sky covered by
cumulus clouds in
cell (i ,j) at hour
P4
Stage
WWO
RAW
w (i,j,t ) cumulus updraft
velocity scale
(m/sec)
crT1(i,j,t ) fraction (0J> V
cold layer
rate (m/sec)
P7
269
-------
YYYYYYYYYYY
VV v V
\J7
270
-------
SECTION 12
PROCESSOR P15
DEVELOPMENT
Processor P15 computes deposition resistances r and deposition velocities
P for those pollutant species for which empirical data are available. These
parameters are related through the friction velocity u*, which is a function
of local wind speed. To avoid coupling this processor to processor Pll, which
computes the flow fields, we will use friction velocity estimates based on the
"raw" surface wind observations, namely those generated by P3, which are
available in the PIF. These "raw" values constrain the magnitudes of the wind
velocities generated by Processor Pll; consequently, any inconsistencies
between deposition velocity estimates based on the u* values derived from the
raw data and those implied by the output of Pll are most likely much smaller
than the level of uncertainty in the empirical relationship between (3 and u*
that we employ here.
Step 1
Compute deposition resistances for each pollutant species. In this task
we must adopt a numbering convention for the pollutants. In Part 1 of this
report, we initiated the convention that species NO is pollutant "1"; N02 is
pollutant "2"; and 03 is pollutant "3". Beyond this the numbering is
unspecified and one is free to select whatever system is convenient. In this
271
-------
Processor description we will refer to all pollutants other than NO, NCL, and
Oo by their names and leave the assignment of numbers to them as a task to be
performed during the implementation phase of this work.
The deposition resistance r(x,t;a) (sec/m) for species a averaged over
the NEROS grid cell centered at x at hour t is assumed to be
10 10
r(x,t;a) = I T(x,n)rn(a,L)/I T(x,n) (15-1)
n=l n=l
where rn(a,L) is the deposition resistance of species a over land use type n
during atmospheric conditions characterized by Obukhov length L at hour t. In
Tables 15-1, 2, and 3 we list the r values currently available for several
pollutants.
Table 15-1. Deposition resistances (sec/m) for S02 as a function
of land use type n and stability L
(From Sheih, et al., 1979)
n= L<0 L = « L>0
1
2
3
4
5
6
7
8
9
10
900
100
100
100
100
100
0
0
70
100
900
300
300
300
300
300
0
0
250
300
272
103
3 — >
J
103
3 I
3 (
103 |
103
X
0
0
800
103
!
)» use 0.0 if
L > 0 and
RUM> 0.9
1 **
\
-------
The surface relative humidity RH(x,t) in each cell can be interpolated from
the surface station values RH available from P3 using a weighting scheme like
that employed in Eq. 12-8.
Table 15-2. Deposition resistances (sec/m) for ozone as a function
of land use type n and stability L.
(From Wesely, 1981)
n= L<0 L=» L>0
1
2
3
4
5
6
7 lake
ocean
8
9
10
300
70
150
60
150
70
3
0.0
100
100
400
200
200
300
400
300
2.103
0.0
250
200
400
400
300
1.5-
1.5-
1.5-
2.1
0.
300
350
103
103
103
.3
0
Resistances for other pollutant species are not available in the current
literature. Deposition velocities for several pollutants other than S02 and
03 onto "vegetation" were reported by Hill and Chamberlain in 1975 (see Table
6 of McMahon and Dem*son (1979)), but these pertain to only one land use type
and moreover, the associated atmospheric stability conditions were not
specified. Lacking any other source of data, we propose to use the Hill and
Chamberlain data to estimate the resistances of each species relative to 03,
273
-------
and subsequently to use that estimate and the values given in Table 12-2 to
deduce crude estimates of resistances. From the Hill and Chamberlain data we
estimate that the difference A^Ca) between the resistance of species a and
that of ozone over agricultural land (n=2) (under unknown conditions of L) are
as follows.
Table 15-3 Deposition resistances of several pollutants relative to that of
ozone over agricultural land. (Deduced from data of Hill and
Chamberlain, 1971)
Species Ar2 [=r2(03)-r2(a)](sec/m)
°3
so2
CO
NO
PAN
N02
0
+25
<-105
-940
-100
+ 7
If we assume that these values apply for L<0, then we can estimate
r2(PAN,L<0), for example, to be
r2(PAN,L<0) = r2(3,L<0)-Ar2(PAN)
= 70 + 100
= 170 sec/m.
One conflict that is readily apparent with this method is that the value of
r2(S02, L<0) deduced from Table 12-3 does not agree with the value given in
Table 12-1. This discrepancy is indicative of the wide scatter in the
274
-------
reported deposition rates of the various species. In the interim, we propose
that approximate values of rn(a,L) be derived from the Ar2 method above. That
is, we assume
r2(a,L<0) = r2(3,L<0)-Ar2(a) (15-2)
for species a = CO, NO, N02 and PAN. The values for the other land use
categories n=l,3,4...10 and other L values can be taken as having the same
ratio to the n=2, L<0 values as those given above for ozone.
The output of Step 2 should be deposition resistances r(x,a) (sec/m) for
species a=03, S02, CO, NO, N02 and PAN for each cell x of the NEROS grid for
each hour of the simulation period. These values are used in the next Step to
compute deposition velocities.
Step 2
Compute the average deposition velocity p(x,t;a) (m s ) of each species
a in each NEROS grid cell at hour t using the following expression from Sheih
et al. (1979):
P(x,t;cr) = 0.4u*[ln(z/zJ+ 2.6 + 0.4u*r(x,t;a) - ^f1 (15-3)
where u* is the friction velocity (m s ) in the cell centered at x at hour t
(from P4); ZQ is the effective surface roughness (m) for this cell (from P4);
z is the elevation at which the concentration of species a is taken for the
purpose of estimating the deposition flux from p (we will use
z = h Az - 10m (15-4)
275
-------
see Part I, Eq. 5-6b); and the function 41 is given by (Sheih et al., 1979)
{0.598 + 0.3901n (- /L) - 0.090[ln(-VL)r}, if L<0; (15-5)
- 5z/L, if L>0;
where L is the Obukhov length (m) for the cell centered at x at hour t (from
P4).
The inputs and outputs of each step of Processor P15 are summarized in
Table 15-1; and the processor and its interfaces with the processor network
are illustrated in Figure 15-1.
276
-------
Table 15-1 Summary of the input and output parameters
of each Step of Processor P15.
Input
Parameter
Description
Source Step
Output
Parameter
Description
T(x,n)
RHn(t)
Ux,t)
fraction of NEROS cell RAW
centered at x in land
use type n, n=l,2,...10.
(dimensionless)
relative humidity P3
(dimensionless) at
surface weather station
n at hour t.
Obukhov length (m) in P4
cell centered at x at
time t.
r(x,t;a) effective deposition
resistance (sec/m) of
cell centered at x
to deposition of species
a=0,, NO, N0~, etc. at
hour t. *
surface roughness P4
(m) of cell centered
at x.
u*(x,t) friction velocity P4
(m/sec) in cell at x
at hour t.
L(x,t) Obukhov length (m) in P4
cell at x at hour t.
r(x,t;«) deposition resistance Step 1
(see Step 1 output).
P(x,t;a) effective deposition
velocity (m/sec) of
species a=03, NO, N02,
etc., in ceil centered
at x at hour t. Reference
level z=10m)
277
-------
I
X
w £
.ti -C
(0 -|
Ifl >
r—
(O 3 O
111
U)
£
3
OJ
V V V V
278
-------
SECTION 13
THE B-MATRIX COMPILER
INTRODUCTION
The b-matrix compiler (BMC) translates the variables generated by the
network of processors into the elements of the "b-matrix" defined in Part 1,
Section 9, and into the variables required in the transport and diffusion
terms of the "r-equation," also defined in Part 1, Section 9. The interface
between the BMC and the processor network is the model input file (MIF). This
is illustrated schematically in Figure 1-1.
The mathematical relationships among the MIF variables and the b-matrix
elements were summarized in Appendix C of Part 1. A complicating aspect of
these relationships is that the values of a small group of the b-matrix
elements are dependent on the local ozone concentration, which is not known
until during execution of the model. Fortunately, all of the matrix elements
in this group can have only one of two predeterminable values, depending upon
the local value of the function IL which is defined as follows:
, if ^' (16-1)
0
O, otherwise
where ?N« is the local emission rate of NO in Layer 0 and
(16-2)
279
-------
Thus, the BMC must provide to the model both of the possible values of each
b-matrix element that is concentration dependent. The appropriate value can
then be selected during execution of the model based on the local value of IL.
In the next section we present the expressions for each of the b-matrix
elements in terms of the MIF variables. (Definitions of each variable can be
found in Table 13-1 at the end of this section.)
THE B-MATRIX ELEMENTS
Let
" (1'aTi)(ft/i+V
* =
6 = w_A. - w+\+(l-|)(l-a)Q
(16-4)
Ozone (03)
bn=
-------
Nitric oxide (NO)
.} 9lNO ' % <°3>1' 1f U0 =
, if UQ = 0.
where
„* - cNO
Nitrogen dioxide
. , NO, N00
b!2 = K
b13 = 0 (16-11)
N0 (16-13)
**- = 1 ' 1f Uo = l>
g = i 1NU2 N02 d 1 o
1 a** U = 0
91N02 ' Uo U'
281
-------
where
N02 (l-aTO)(l-cQ
91NO, = Sl r h, (pwn + u) UJNO (16-18)
C -L llUn
N02 (l-aTn)(l-o)u _
** = c ^
*i kin "™ ^T
'Dl;
g1NO s Lii «%Q0 (16-20)
1NU2 h1ONQ + u) U3
species x (excluding 03, NO and
b = + fi + 1-
b!2 = E^-OlAm
b,, =~0 (16-23)
all species
(16-25)
b22 = H [
u n2
h _ 1 /, > f8 (16-27)
D23 ~ " ~h7U"CTcJ A8
^ (16-28)
UI/T fw —w ^
~|J/,V."0 ™rJ ffl
+ CiJ C + S d-a. + a.f)] (16-26)
282
-------
Ozone (Q3)
31
(16-29)
b,9 = r (l-f)M
*<• n3
1 ,
1% — / •• I f » LI
DOO — 7- I M*n.
jo rU
i* = i H U C°3 if U
- K "-J0-?^13 > ' ' u,
3 3
= 0.
(16-30)
(16-31)
(16-32)
Nitric oxide (NO)
h -
b31 -
b33 - y-
' if
' if
where
NO
3NO
(16-33)
(16-34)
(16-35)
(16-36)
(16-37)
283
-------
Nitrogen dioxide (N09)
where
NO
284
13 (16-39)
3NO " h(3
b = —£- + l~^ (16"40)
b3? = K (I'^M (16-41)
Ot. ''0 x
b33 = H (" M+ |!|3U3) (16-42)
cO^,, if U =1;
J -1 ° (16-43)
9
2 ] (16-44)
= U- -- MO NOO] (16-45)
-------
species x (excluding 0-, NO and N02)
where
«.-!
0, if H, < 0.
b31 =
b32 = H(1~'1')M (16-48)
b33 = J [-M+ H3U3] (16-49)
1, if H3 > 0;
U =
-------
Consequently, the basic equation 2-29 of Part 1 must be cast in a form
compatible with (16-54). The only terms in (2-29) that are involved in this
transformation are
_3
3t
3
31 nV
3.
Vv>j
81 nV
(16-55)
= a are metric factors and a is the radius of
where u, = (a cos .. In Section 7 of Part 1, we
J
proposed to express the fluxes ., etc. associated with the subgrid scale
J
fluctations in the gradient transfer form (see Eq. 7-2, Part 1)
. = -K.
8.
8\
(16-56a)
3 .
(16-56b)
where the K. are diffusivities specified in Stage K of processor P12. Making
J
use of (16-56) in (16-55) we get
3.
3K, 3,
31 nV.
31 nV. 3K. 3.
J _ ..2 JT J
55>K "<*> 3A J JJA
286
-------
2 2
3 . 3 • •
— " — =
This equation is of the form (16-54) required by the algorithm in the model
that handles the advection and diffusion processes (see Appendix A of P7).
Two basic operations are now necessary to complete the preparation of the
equation for numerical solution: (1) to convert the units of all
parameters — the model treats the equations in (A.,) space rather than
physical space (x,y); (2) to convert the effective advection velocities,
i.e., the terms in brackets in (16-57), into coordinates (\*,*) of the
upstream trajectories associated with each grid point in the model domain (see
Figure 7-A1 of P7). Below we outline the details of specific operations that
are required.
Step 1
. Four parameter fields are involved in these operations: the two velocity
components . and ., the diffusivites K.; and the cell volumes V..
J J J J
Values of each of these are available in MIF for each layer j=l,2,3 and for
each grid cell in the model domain. All data in the MIF are in mks units.
-1 2 ~1
This means that the units of and are m*s ; K has units of m s ; and
3
V has units of m . The first task is to convert all length units from meters
to radians.
This is accomplished using the metric factors u, and u^ which are
measures of the arc angles in radians of longitude and latitude, respectively,
per unit length on the earth's surface. Thus, for example, I-V^.- is the
east-west component of wind speed in radians (of longitude) per second.
-------
Let (I,J) denote (column, row) of any cell in the model domain. The
metric factor u,. varies with J, namely
ux(J) = (a cos ^j)"1 (16-58)
where 4>j is the latitude in radians of row J, but u, is a constant, namely
Mf = a"1 (16-59)
where a is the earth's radius:
a = 6,367,333. meters (16-60)
Now convert the velocity and diffusivity fields obtained from the MIF into
their corresponding values in (A.,(j>) space for each layer j=l,2,3 and grid
point (I,J) as follows:
J)>.j, (16-61)
j, (16-62)
= [HX(I,J)J2K..(I,J), (16-63)
(16~64)
Kyj(I'J)
The K*. and K* . fields should be recorded directly in the output file of the
^j yj
BMC for input into the model.
Step 2
Take the natural log of each cell volume V. in each layer at each grid
\j
288
-------
point to produce arrays of (InV.). For each layer j compute estimates of
j
g
g-^ InV. at every grid point (I,J) as follows
^ lnV.(I,J) - DXLV(I,J,j) =
(16-65)
where
6\ = k (^=fi) = model grid interval in \.
On the western boundary where 1=1, use
DXLV(I,J,j) = [ln(Vj(2,J))-ln(VJ.(l,J))]/6\ . (16-66)
and on the eastern boundary
DXLV(IMAX,J,j) = [ln(VJ.(IMAX,J))-ln(VjdMAX-l,J))]/6\ (16-67)
In the same manner define
at InV-(I.J) - DYLV(I,J,j) = [ln(V.(I,J+l)) n(. cfi>
O = 1/6 (•o^) = grid interval in <}>. (16-69)
On the south and north boundaries, where J=l and JMAX, respectively, use
approximations similar to (16-66) and (16-67) for DYLV(I,J,j).
Step3
Form approximations of the derivatives of K*. and K*. as follows:
DXKX(I,J,j) =[K*.(I+1,J)-K*.(I-1,
(16-70)
289
-------
9KJi
gjU- = DYKY(I,J,j) = [KJj(I,J+l)-K*j(I,J-l)]/(26(t))
Use approximations similar to (16-66) and (16-67) to treat grid points on the
boundaries.
Step 4
Form the effective advection velocity components, i.e., the terms in
brackets in (16-57) as follows.
ueff(I,J,j) = * - KJj.d.
- DXKX(I,J,j) (16-71)
veff(I,J,j) = J - K*.jd,J)*DYLVd,J,j) (16-72)
- DYKYd.J.j)
Values of u ., and v - , must be computed for each layer j=l,2,3; each grid
point in the model domain; and each time step At (= 30 min.).
Step 5
Compute the "back track" points associated with each grid point (I,J),
each layer j=l,2,3, and each time step. This is done as follows:
Let ueff (I.J.j.N), veff(I,J,j,N) denote the
effective velocity components , (16-71) and (16-72), at time step N, i.e., at
time t=t +NAt. Consider a given grid cell (I,J,j) at a given time step N.
Compute
290
-------
Ml = "ueff (I.J.J.N)AT (16-73)
A4>1 = "veff d.J.J.N)AT (16-74)
where AT = At/n, say AT = At/3 = 10 minutes. These increments define a point
(A.,) in (X,<5>) space, namely
(16-75)
Fit a biquintic polynomial to the (ueff, v *f) values at time step N on
the 36-grid-point square centered at the cell nearest (A.,, <)>,). Note that the
polynomial must be in terms of (\,) rather than (x,y). With this polynomial
and a linear interpolation in time between NAt and (N-l)At, estimate the
values of (ueff, veff) at- the point (A.-^,^) at time NAt-AT , i.e.,
U
- AT), Vef fO^pj, NAt-AT).
(16-76)
C"V
eff ..
Now compute the new point
\9 = X, + [-u_ff (A,,,,j,NAt-AT)AT]
1 L eft 1 1
eff (^1»<()1»J»NAt-AT)AT].
Using the biquintic space approximation and the linear time interpolation
again to estimate u ~ and v -f at (A^, (fp. NAt-2At), compute finally
(A*=) \RT. = LAMBT(I,J,j,N) = \? + [-uoff(\?,<)>?, j,NAt-2AT)AT]
BTj 2 eft Z 2 (16-77)
(
-------
The diffusivity fields K*. and K*. should also be recorded on the
b-matrix tape for all I,J,j and all times steps N=1,...NMAX.
The b-matrix compiler is illustrated in Figure BMC-1, and each of the
variables in the MIF is defined in Table BMC-1.
292
-------
INPUT
(MIF)
BMC
OUTPUT
("b-matrix" tape)
Ml!
CD ^ W
IK-J=
.2 «> o
Hi-
•g at2
«£.S 2,
« S E S
5s^l
•*• x: "® m
5 w 5 O
o — w O
S So «
^r^ °
3 O CO O
= ^
C/> »- O XI
m
£
3
O)
293
-------
Table BMC-1. Definitions of parameters in the
Model Input File (MIF).
Parameter
Definition
JTO
o
Tl
w,
w
Thickness of model layer n, n=l,2,3.
Average east-west wind component in model layer n=l,2,3.
Average north-south wind component in model layer n=l,2,3.
Horizontal eddy diffusivity in model layer n=l,2,3.
Fraction of the top surface of Layer 0 penetrated by
terrain in a given grid cell.
Fraction of the top surface of Layer 1 penetrated by
terrain in a given grid cell.
Fraction of sky covered by cumulus clouds in a given cell.
Plume entrainment velocity in Layer 0.
Fraction of Layer 0 occupied by line and point source plumes.
Deposition velocity of pollutant species X.
Fraction of the top surface of Layer 0 covered by ascending
fluid.
Mean speed of upward moving fluid on top surface of Layer 0.
Mean speed of descending fluid on top surface of Layer 0.
Composite vertical turbulence parameter on top surface of
Layer 1.
294
-------
Table BMC-1. (continued)
Parameter
Definition
'Ira
WD1
"2
"c
fe/Ae
x
X
1
X
2
x
oo
Composite vertical turbulence parameter on top surface
of Layer 1.
Mean vertical air velocity on top surface of Layer 1
induced by flow divergence in Layer 1.
Mean vertical air velocity on top surface of Layer 2.
Cumulus updraft velocity scale.
Turbulent entrainment velocity at mixed layer top.
Fraction of cumulus cloud volume flux drawn from Layer 0.
Local time rate of change of the elevation of the top
of Layer 2.
Volume flux across top of Layer 3.
Emission rate of surface sources of pollutant x.
Emission rate of sources of pollutant x in Layer 1.
Emission rate of sources of pollutant x in Layer 2.
Concentration of species x in air above the top of model
Layer 3.
295
-------
REFERENCES
Artoz, M. A. and J. C. Andre, (1980): "Similarity Studies of Entrainment in
Convective Mixed Layer", Boundary Layer Meteor., Vol. 19, pp 51-66.
Brost, R. and J. C. Wyngaard, (1978): "A Model Study of the Stably Stratified
Planetary Boundary Layer", J. Atmos. Sci., Vol. 35, pp 1427-1440.
Bullock, 0. R. (1983): "Spatial and Temporal Interpolation of NEROS Radiosonde
Winds", M.S. Thesis, Dept. of Marine, Earth and Atmospheric Sciences,
N.C. State University, 105 pages.
Demerjian, K. L., (1983): personal communication.
Demerjian, K. L. and K. L Schere (1979): "Applications of a Photochemical Box
Model for Ozone Air Quality in Houston, Texas". Proceedings, Ozone/
Oxidants: Interactions with the Total Environment II, Houston, TX, 14-17
Oct. 1979, APCA, Pittsburgh, PA, pp 329-352.
Godowitch, J. M. and J. K. S. Ching (1980): "Formation and Growth of the
Nocturnal Inversion Layer at an Urban and Rural Location", Proc. Second
AMS Joint Conference on Applications of Air Pollution Meteorology, New
Orleans, LA, pp 165-172.
Golder, D., (1972): "Relations Among Stability Parameters in the Surface
Layer", Boundary Layer Meteor.. Vol. 3, pp 47-58.
Haltiner, G. J., (1971): Numerical Weather Prediction, John Wiley and Sons,
New York, NY. 317 pages.
Hill, A. C. and E. M. Chamberlain, (1976): "The Removal of Water Soluble Gases
from the Atmosphere by Vegetation", Proc. Symp. on Atmosphere - Surface
Exchange of Particles and Gases, ERDA Symp. Series, NTIS, pp 153-170.
Holtslag, A. A. M. and A. P. Van Ulden, (1983): "A Simple Scheme for Daytime
Estimates of the Surface Fluxes from Routine Weather Data", J. Climate
and Applied Meteor.. Vol. 22, pp 517-529.
Jones, F. L., R. W. Miksad and A. R. Laird, (1981): "A Simple Method for
Estimating the Influence of Cloud Cover on the N0? Photolysis Rate
Constant", J. Air Poll. Cont. Assoc., Vol. 31, pp 42-45.
Lamb, R. G., (1983a): "Air Pollution Models as Descriptors of Cause-Effect
Relationships", Atmos. Envir., (in press).
Lamb, R. G., (1983b): "Theoretical Issues in Long Range Transport Modeling",
Preprint Volume, AMS Sixth Symposium on Turlulence and Diffusion, Boston,
Mass., pp 241-244.
Lamb, R. G. (1983c): "Causality and Atmospheric Phenomena" (in preparation).
296
-------
Lamb, R. G. (1983d): "A Regional Scale (1000 km) Model of Photochemical Air
Pollution. Part 1: Theoretical Formulation", EPA report EPA-600/3-83-035.
226 + xi pages. NTIS-PB83-207688.
McMahon, R. A. and P. J. Denison, (1979): "Empirical Atmospheric Deposition
Parameters - A Survey", Atmos. Environ., Vol. 13, pp 571-585.
Melgarejo, J. W. and J. W. Deardorff, (1974): "Stability Functions for the
Boundary-Layer Resistance Laws Based Upon Observed Boundary-Layer
Heights", J. Atmos. Sci., Vol. 31, pp 1324-1333.
Nieuwstadt, F. T. M. and H. Tennekes, (1981): "A Rate Equation for the
Nocturnal Boundary-Layer Height", J. Atmos. Sci.. Vol. 38, pp 1418-1428.
Sheih, C. M., M. L. Wesely and B. B. Hicks, (1979): "Estimated Dry Deposition
Velocities of Sulfur Over the Eastern United States and Surrounding
Regions", Atmos. Envir., Vol. 13, pp 1361-1368.
Wesley, M. L., (1981): "Turbulent transport of ozone to surfaces common in
the eastern half of the United States" submitted to
Advances in Environ. Sci. and Tech., Vol. 12.
Zeman, 0., (1979): "Parameterization of the Dynamics of Stable Boundary
Layers and Nocturnal Jets", J. Atmos. Sci.. Vol. 36, pp 792-804.
Zilitinkevich, S. S., (1972): "On the Determination of the Height of the
Ekman Boundary-Layer", Boundary Layer Meteor., Vol. 3, pp 141-145.
297
-------
TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing)
1. REPORT NO.
4. TITLE AND SUBTITLE
A REGIONAL SCALE (1000 KM)
AIR POLLUTION Part 2. Inp
2.
MODEL OF PHOTOCHEMICAL
Lit Processor Network Design
7. AUTHORIS)
Robert G. Lamb
9. PERFORMING ORGANIZATION NAME AND ADDRESS
Same as Block 12
12. SPONSORING AGENCY NAME AND ADDRESS
Environmental Sciences Research Laboratory— RTP, NC
Office of Research and Development
Environmental Protection Agency
Research Triangle Park, North Carolina 27711
3. RECIPIENT'S ACCESSION-NO.
5. REPORT DATE
6. PERFORMING ORGANIZATION
8. PERFORMING ORGANIZATION
CODE
REPORT
10. PROGRAM ELEMENT NO.
CDVIA1A/02-1335 (FY-84)
11. CONTRACT/GRANT NO.
13. TYPE OF REPORT AND PERIOD COVEF
In-house
14. SPONSORING AGENCY CODE
EPA/600/09
16. SUPPLEMENTARY NOTES
16. ABSTRACT
Detailed specifications are given for a network of data processors and submode
that can generate the parameter fields required by the regional oxidant model for-
mulated in Part 1 of this report. Operations performed by the processor network
include simulation of the motion and depth of the nighttime radiation inversion la>
simulation of the depth of the convective mixed and cloud layers; estimation of the
synoptic-scale vertical motion fields; generation of ensembles of layer-averaged
horizontal winds; calculation of vertical turbulence fluxes, pollutant deposition
velocities, parameters for a subgrid-scale concentration fluctuation parameter!'zati
scheme; and many other functions. This network of processors and submodels, in
combination with the core model developed in Part 1, represent the EPA's first-
generation regional oxidant model.
17. KEY WORDS AND DOCUMENT ANALYSIS
a. DESCRIPTORS
13. DISTRIBUTION STATEMENT
RELEASE TO PUBLIC
b.lDENTIFIERS/OPEN ENDED TERMS
19. SECURITY CLASS (This Report)
UNCLASSIFIED
20. SECURITY CLASS (This page}
UNCLASSIFIED
c. COSATI Field/Group
21. NO. OF PAGES
22. PRICE
EPA Form 2220-1 (9-73)
------- |