600383035B
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
-------
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
i
After the model described in Part 1 of this report was formulated, a
draft of an instruction manual was rather hastily prepared to guide computer
Drogrammers 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
c\.
abjective is to provide in conjunction with Part 1 -of detailed description of
--hat 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.
-------
CONTENTS
Deface i i i
•stract iy
gures vi i i
Dies xi
iknowledgements xiii
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 monitoring data 34
-------
CONTENTS (continued)
Stage 1C: Estimating initial conditions from
monitoring data — 34
Stage UBC: Upper boundary conditions 36
4. Processor P3 39
Genera] 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 friction force 61
The Coriol is force 65
The momentum equations 65
The fluid depth equation 67
Simplified model equations 76
Solution of the u, v, and zwe 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
Appendix B to Section 7 110
Calculation of the BAR variables 120
Calculation of the PRIME variables 122
Calculation of the dependent variables uX and h' 124
7. Processor 8 131
Introduction 131
Derivation of basic equations 131
Stage 2Q 147
Stage PATH 161
Stage WEWC 165
Stage W2 172
VI
-------
CONTENTS (continued)
Stage WV 182
. Mode 0: Ap!=0 : 182
Mode 1: Ap^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
Development 254
Stage K 254
Stage WWO 261
Stage WW1 264
12. Processor P15 271
Development 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
-------
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 ,;
f
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 zys and zt ................................ 62 3
s»
7-4 Projection of the horizontal rectangle (Ax, Ay) onto the I,.
surface H.,, centered at point Q
»5
7-2 Illustration of the 54 5 rain 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 T and S are taken to derive a biquintic expansion ":
of F(x',tn) about the point (IST,JST) (see Eq. 7-49) 104
7B-1 Flow chart of FIOMOO operations 119
7B-2 Grid network on which 5, jj , p. , S , S , and S. are computed.
Different spatial derivative operators Ax and Sy 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 model domain 129
vm
-------
FIGURES (Continued)
Number Page
7-6 Schematic illustration of Processor P7 and its input and
output interfaces with the processor network 130
3-1 Illustration of surfaces H2 and HS during situations in
which convective clouds cover only a portion of the
modeling region 132
>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 < tj. .. .7: 141
->-3 (a) Idealized profiles of mixing ratio q and potential
temperature 9 in dry, convective conditions.
(b) Second derivatives of the profiles illustrated in
panel a 150
y-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
3-6 Schematic illustration of Processor P8 and its input and
output interfaces with the processor network 196
3-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 aT1
and OTO It 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
IX
-------
FIGURES (continued)
Number
- •
12-2 Illustration of Processor P12 and its input and output
interfaces with the processor network
15-1 Illustration of Processor P15 and its input and output
interfaces with the processor network
• I
BMC-1 Schematic illustration of the B-Matrix Compiler (BMC). -i
The input interface is the Model Input File (MIF). The -|
output of the BMC is the "b-matrix tape" which is read g
by the model code CORE (see Figure 1-1) ..................... . :1
-------
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 K 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
xn
-------
ACKNOWLEDGEMENTS
I am indebted to Dr. Don Pelles 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.
xm
-------
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
-------
ar^^tl
-------
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
-wo vectors P and Q, each of length N, where N is the total number of chemical
jpecies simulated. The n-th element of P is the net rate of production of
:,pecies n due to source emissions and chemical reactions among all other
--.pecies, and the n-th element of £ is the net rate of destruction of species n
;:ue to its chemical interaction with all other species. Thus, any chemical
nineties mechanism can be incorporated into the model as long as it is
.xpressed in a form that is compatible with the vector interfaces that link
:ORE 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-l.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
-------
-atnematics contained within the module CORE were given in Section 9 of Part
We will not elaborate further on this part of the model system in this
•'sort other than to describe in more detail the numerical scheme that we
;oiy to the transport and diffusion portion of the model equations. These
:tails 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
: regional model will be applied. Each point in the figure represents the
;'er of a grid cell in which values of all pollutant species concentrations
. computed by the model.
:MARY OF PROCESSOR FUNCTIONS
The functions performed by each processor in the model network,
•ustrated earlier in Figure 1-1 are summarized below.
'rocessor 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.
-------
"^•%9
CM
I
01
3
O5
-------
Processor P3
Prepares surface meteoro.logical data for use in other processors.
Processor P4
Estimates surface roughness, Obukhov length, surface heat flux, and
friction velocity.
'rocessor 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.
!
I
8
-------
£'JMMARY OF MODEL- EQUATIONS
The equations upon which the regional -scale mode is based were derived in
:.rt 1 of this study. The basic forms of the governing equations in each of
- a model's four layers are summarized below.
/er 1
nere
+ M 3i + u 3i
M\ 9A ^ 30
* V4 C + Fo,l - Fl',lJ
a cos<(»
" = a
a = earth radius at MSL (in meters)
\ = longitude
0 = latitude
A = a2 (Ad>AX) cos(})
A((),AA = latitude, longitude grid cell dimensions
= constants (A0 = J° A\ = V°)
0
-------
X+AA./2 <|>+A/2
- MAX [ZT(\',d>'),
subgrid scale fluxes of c (see Part 1, Section 7) x
'c'>, '
-, = 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'CTT1) [(1 * 2)WLn
p = deposition velocity of species c
a-n(A,<|>,t) = fraction of surface H, penetrated by terrain
w zn wm - w
+ ^ + erf C-^-
W W
j, - \ wm " ^1» neutral and unstable conditions;
wR1 - | ui i
Hi> ( = given inversion layer depth growth rate), stable
10
-------
wn, = mean vertical velocity on H, (terrain induced component
U1 excluded) *"
cr = rms vertical turbulence on H, ( = 0 in stable cases)
w = threshold of cumulus "root" updraft velocity on
o , w (see Layer 2 equations)
c c
•) I _+2
erf(x) = — 1 e l dt
Fo,l =
ov = fraction of surface H in given grid cell penetrated by
terrain
Ozone:
Fo(03) = <03>1 W-X- " V+
- (1 - a)
0, if
3V ^ " NOV , otherwise
PO + v
Nitric Oxide:
FQ(NO) = 1
- NO
- a)
- (1 - a)
- 0) , if
, otherwise
-------
Nitrogen Dioxide:
FQ(N02) = lW_A_ - N02
- (1 - a)
- o)
, if SNQ > v£03
2 , otherwise
all other species, x-
F0(X) = iW_\_ - xw+\+(l - 4)(1 - a) - (1 - a)Kx * ? )(1 + p
a =
tjjcr
c
ft
w2 - w
-i - ffl/\n
1 - ac f0/AO
See Layer 0 equations for specification of O,, NO, N02, x» w+, X+, etc.
Layer 2
at
2
2
2
:'>,
l,2 - 2,
= 2
V2(X,(p,t) = a2cos4> {22(Xl,0',t) - MAX CzT(X',(i)!)
A-M/2 (JJ-A0/2
',f ,t)]}dfdA'
12
-------
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.
P o_(w2 - w ) _ aL f
2,2 = 1-a 2 1* ' A6 [~ (1 ' ac)3 ' CTcCc + 2]
Fl,2 = anP2 + (1 " CT1^1W1 + (i " 2)Wlm + (2 "
Wl'Wlm' see Layer 1; Dl' see Part 1> Eq- 4~
a = fractional area coverage of cumulus clouds
w = mean upward velocity in cumulus clouds
fe
•5 - turbulent entrainment rate of inversion air into mixed layer
« - mean vertical velocity (mean of both terrain induced component wT2
and divergent component w)
-~ = z0 = local time rate of change of mixed layer top
Ot L.
= fraction of cloud air from surface layer
0
-------
w2 - w
vc = - -f
c, c1, 4, w+, A.+ = See Layer 0 equations
Layer 3
9A (p 30
33 . ,, 83 . A r
_
2,3 3,3
3
V3(A,<|»,t) = a2cos<|»
Z2(\,,t) = mixed layer top
f.t) = model top
subgrid scale flux of c (see Part 1, Section 8)
3 = all chemical (wet and dry) rainout and washout processes
3, 3 - layer averaged winds
F W2 " Wc - 322
2,3 = ac C 1 - a + fe/^](cc - 3) * 3 ^
14
-------
3 , if dfij/dt < 0 ;
_ .
3) dH3/dt + 3 — , otherwise
CM = concentration of species c above z3
dhL/dt = given volume flux through model top surface
c = (1 - 4<)2
where c1 and c are defined with the Layer 0 variables.
er 0
<03>0 = (1 - C)
0 w.\,<03>1
3
( °. if SNO > vhj
°3 'I
3V " W^ ' otnerwl'se
0 = (1 - C) NO
N° =
15
-------
+ v(NO - 03) , if SNO >vC03
0 =
NO,
NQ2v + 03v
, otherwise
, if
'NO'" 3
2 , otherwise
Species x other than 03> NO, and N02:
o = (3
X =
w+A+ * (1 - OP,
Parameters:
\+ = %[1 - erf (-r
0 N-J
16
-------
O '9
W Z^
0 - 0
w. = -*0 - — exp
wo
[1 - erf (
wo
v = u* (plume entrainment velocity)
£ = plume volume fraction
-------
SECTION 2
PROCESSOR PI I
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 Im "vectors"
[SPD,DIR,T,TD,p]in), 1=1,...Im (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 J
(millibars) at a given station m. The number of observation levels I in each I
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. f
\
Step 1
Convert the speed and direction (SPD.OIR)^ into cartesian components
(u,v) at each level i and at each station m by
18
-------
/here
9. = (DIR. )-(at/360). (1-3)
i ni ini
hese values must be processed further (in step 10) before recording in the
'IF.
•:tep 2
Convert the temperature and dew point depression TD into the water vapor
.nixing ratio q (mass of water per total mass of air) at each level i and
station m by
qim = -e (1"4)
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)
D- = T. - TD. + 273 (dew point temperature at level i, station m)
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 . i and station m by
where e is found from (1-5) by replacing the dew point I-, with the
temperature T. .
Step 4
Convert the pressure p^ at the levels i=l,2,...I of the upper air
observation at station m to altitudes z. (AGL) by first converting the
temperature measurements T. to virtual temperatures
). d-7a) j
Then |
where
p. = pressure at level j (j=0 represents ground level) at station m;
z. = elevation (AGL) of observation level j at station m;
g = 9.8 m sec" (gravity);
-2 _1
Rd = 287 m2 sec K (gas constant for dry air)
T . = virtual temperature (°K) at level j, station m.
20
-------
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_-,m = sea-level pressure at station m (in millibars)
and we calculate the geopotential height Qm (m2/sec2) of the 1000 mb surface:
*om ' 92m * RdTvoJn OTJ • d'9
Using the *om value from the previous observation time, say, t-At, compute
Record 4>om and $Qm 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
-------
I
where
z = zm + nAz , n=0,...
Az=50 meters, and
5000-z
m
SO
(1'10e)
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 T f qm, TDm, and RHm values for each of
the z given by (l-10d) in the PIF.
Step 7
The pressure-height pairs1 (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.,
=>
where
z=nAz, n=0,l,...100;
(1-lla)
(1-llb)
and Az=50 meters. In performing the interpolation (1-lla), an exponential
formula should be employed based on the hydrostatic condition
Ap_ _ _
(1-12)
22
-------
"ecord the pm values at each z given by (1-llb) in the PIF.
;:ap 8
Using the P_(z) profile obtained in Step 7, determine the geopotential
-2
i2sec ) of the 850, 700, and 500 millibar surfaces at each station m. For
/700)m and *(500)m*
Step 9
Using the virtual temperature profile Tvm(z) from step 6, Eq. l-10a
above, and the pressure profile p (z) from Step 7, Eq. 1-11, compute the j
potential temperature profile 6m(z) for each level z defined by (l-10d):
V« - Tm WCgs,)-
'I
I I
where p is in~millibars, as above, and Tm is in degrees Kelvin. Now compute
the air density at each of the same levels z using
23 !
-------
n -2 „ -1
where pm is in millibars and Rd=287 m2sec °K , Tvm is in °K, and p=kg m
From 9 and p compute the profile a (z) of static stability from
- [8m ( z+Az ) -6m ( z- Az ) ] • 10"
CTsm(z) =
_3 _2
With p in millibars and p in kg m , a has units of m4sec2kg . Record the
profiles 9 (z),p (z) and cr (z) in the PIF for each rawin station m.
Step 10
Interpolate the wind components uim and v^ of Eq. (1-2) above to obtain
the profile
vim
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
i.nd 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
Description
Source
Step
Output
Variable
Description
SPDim
wind speed (ms ) at RAW
level i at rawin
station m.
Wind direction RAW
(compass degrees) at
level i at rawin
station m
u. (t.) east-west wind
component at level
i, station m, hour t^.
v. (tk) north-south wind
component at level
i, station m, hour t^.
pressure (mb) at RAW
level i, station m.
dew point depression RAW
(°C) at level i,
station m.
temperature (°C) at RAW
level i, station m.
q- (tk) mixing ratio
(dimensionless) at
level i, station m,
hour t^.
Tn. (tu) dew point temperature
mm K at level i, station m,
hour t,..
pim
-------
Table 1-1 Summary of input and output variable
of each step of Processor PI. (Continued)
Input
Variable Description Source Step
zm
in
p
worn
V Ulll
vim
V 1 III
Him
i in
RHim
i HI
2im
TDim
p.
zim
1 III
"™(z)
Pm(z)
m
T (z)
vm
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
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
*om(V
Ulll K.
*om(tk}
Ulll N
Tvm(z'V
VIII ls\
qm(z,tk)
III N.
RHm(z,tk)
III N
TDm(z'V
Pm(z,tk)
III K
!(850)m
!(700)m
*(500)m
e.(z)
in
Pm(z,tk)
m K
Description
geopotential (m2s~2)
of 1000 mb surface
at station m, hour tk<
time rate of change
(m2sr3 ) of geo-
potential of 1000 mb
surface at station m.
Temperature, mixing
ratio, relative humidity,
and dew point profiles
at station m resolved
to Az=50m as prescribed
by Eqs.
(l-10d,e) at hour tR.
pressure (mb) at
elevation z(msl) above
station m, resolved to
Az=50m as prescribed
by Eqs. (l-10d,e) at
hour t..
geopotential (m2sec~2 )
of the 850, 700, and
500mb surfaces at
station m.
potential temperature
(°K) at elevation z
over station m, resolved
to Az=50m as prescribed
by Eqs. (l-10d,e).
o
air density (kgm~ ) at
elevation z (resolved
as above) over station m,
hour t. .
sm
27
:. ) static stability
K (m4s2kg-2) at elevation
z at station m, hour t..
-------
Table 1-1. Summary of input and output variables
of each step of Processor PI. (Concluded).
Input
Variable
Description
Source
Step
Output
Variable
Description
uim
im
see above
see above
Step 1 10
Step 1
u(x ,z,t.> east-west and nor
~m K south wind
) at elevation 2 (a*
given by Eqs. 1-K
at location x of
station m, at hour
V
^m''
[also
28
-------
YYYYYYYYYYYYYYY
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
initia!7 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"1.)
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
"I?-
scheme for use in treating the pollution chemistry. Proceeding along lines ^|
similar to those described there we could formulate an approximate way of f
extracting volume-averaged concentration estimates from point measurements. $
f
However, we will not attempt this here because the improvement in accuracy |
gained in the initial concentration fields would probably not be significant ;|L
enough to warrant the development and implementation efforts. Perhaps future ||
modeling studies can investigate this problem in detail if the need is great. "I
1
;.K.
In the remainder of this section we outline a procedure for obtaining J
-Tt
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 Os, 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-(ppm) .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 %. (*B, 4B, t)(ppm) at each grid point (Ag, <|>B) on the
boundary. Here (A,) denotes the longitude and latitude of a cell center
on the boundary of the regional model domain.
(3) Using the functions x ^ obtained in step 2, estimate layer averaged
concentrations along the boundaries as follows:
= B Xj (Xo, +B, t), 11=1,2,3 (2-1)
" ° i=l,...IMAX
where the B are empirical constants to be estimated from the NEROS field
n
experiment data.
(4) For each hour tm, each boundary point (\B, <)>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 monitoring 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 1 and 2 can be
equated with minimum error.
I) 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-j(X, 4, 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:
Xi(A, 4>, V (2-2)
2 = X1(A, , t0) (2-3)
3 = SjXfC*. + , V (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 x^
for., all cells and all hours. Record these in the ICBC portion of
MIF(units = ppm for all species), i.e.
t} = XCLEAN, i
all (X,
-------
Table 2-2 Summary of input and output variables of each
stage of Processor P2.
Input
Variable
Description
Source
Stage
Output
Variable
Description
X,'i,(t ) concentration (ppm)
* m of i-th pollutant
at hour t measured
at surface monitoring
station k=l,...K.
RAW
LBC
average concen-
tration (ppm) of
i-th pollutant
in layer n=l,2,3
at boundary cell
(\B, <)>g) at hour
m'
Xik(tm) see Stage LBC input.
RAW
ic
, n
average concentration
(ppm) of i-th pollut-
ant in layer n=l,2,3
in grid cell (\,$) at
the initial instant
v
see
RAW
UBC
, 0>,
concentration (ppm)
of i-th pollutant
at top surface of
model over grid
cell (\,4>) at hour
m
37
-------
CM a.
o- i-
1YY
c
A
E
m
c
A
o
V
3
Q.
*j
3
O
•o
c
CO
3 .
0. j*
.£ 5
51
_ Q)
= 5
« o
^
o
CO w
8 2"
0)
Is
s °
2.S
15
CN
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,T,TD,P]n
where n denotes the surface station, whose location is x ; and the other
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:
un = 0(xn) = - SPDn-sin 9n (3
vn = v(xn) = - SPDn-cos 6n (3
39
-------
where
9n = DIRn*(27t/360)-
The components (&n,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
q = — (=mixing ratio) (3-2a)
pn n
where
111
R = 0.461 joule g"1 °K~1
L - 2500 joule g"
eQ = 40mb (saturation vapor pressure at temperature T )
TQ = 302 °K
Pn = station pressure (mb)
and
TDn = Tn " TDn * 273 (=dew point temperature) (3-2b)
p -e
RHn = qn On62||| (relative humidity) (3-2c)
sn
where e is obtained from (3-2a') by replacing the dew point temperature T-
in the formula with the dry bulb temperature Tn (expressed in °K);
TV = T (1 + 0.61qn) (=virtual temperature) (3-2d)
40
-------
where T (and hence T ) is in °K.
The mixing ratio q (dimensionless) and the relative •humidity RH
(dimensionle'ss) 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 <°c> = V°
Tvn<°c> = Tvn<°K>
where TQn (°K) is from Eq. 3-2b and Tyn (°K) is from (3-2d).
Step 3
Let 2n be the elevation (meters, MSL) of surface station n. Compute the
geopotential of the 1000 mb surface
*on ' 9*n * RdTvn ln I5B5 <3'3>
I
where Tyn is from (3-2d) above (in °K) and Rd = 287 m2s" °K~ . Using the
value of * from the previous time interval, estimate the time rate of change
of 4> by
on
*on = C*on(t) ' Oon^^"1 (3"4)
Compute the sea-level pressure p,.-, at each station site x ,n=l,2,...N as
follows.
p2- ] (3-5)
vn
41
-------
where p 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
6n = Tn (M2}0.286 (-p0tential temperature)
^
(3-6)
10~
-2
vn
(=density, kg m
r3)
(3-7)
where T is from Eq. 3-2d and is in degrees Kelvin (°K), p is the station
pressure in millibars and
Rd = 287 m2sec
"2
"
Record 0n(°K) and pR (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
SPD
n
DIR
n
Tn
Description
surface wind speed
(m/s) at station n
surface wind direction
(compass degrees) at
station n.
surface temperature
fOr^ = +• <-4-^-f -inn n
Source
RAW
RAW
RAW
Output
Step Variable
1 u (t )
W
2 «n
-------
YYYYYYYYYYY
(A
03
O
CD
k»
CD
b
Q.
^*
3
O
TJ
C
CD
**
3
Q.
w
^rf
•o
C
CD
CO
Q.
w
O
(A
(0
0)
O
- o
o I
s ®
.2 c
i: co
(O (0
18
.2 o.
(D O
O>
iZ
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 (i, ±) 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 1
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 (un>v_) at
each station and at each hour t , estimate the wind speed u in the given
grid cell at hour t by performing a weighted average of the observations
(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
|u| m < _x
4, otherwise
Next, using the fraction cr^y (x,tm) of local sky coverage by all cloud types
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 < 0.5
\ nighttime hours only
.-2, a < 0.5
47
-------
«- -.«. ^ of
greater than 30%, i.e. ,
of C to usage (
> 0.30, C (x) > 1
6~ ~ (4-3a)
- --
if T(x,7)>0.6, Ce(x) = ' daylj9ht hours
L 0, night.
m J
L = [ (v * v3) • *o(bl"' "
where
«1 = 0.004349,
a2 = 0.003724,
4 = 0.5034,
b2 = 0-2310,
b = 0.0325.
s
-SignCCJ ,
48
-------
where
1, if Ce > 0;
0,.ifCe = 0;
-l, if Ce < 0,
Step 2
Estimate the effective surface roughness ZQ in each cell of the NEROS
grid using the following expression
10 10
ZQ(X) = I T(x,n)zo(n)/I T(x,n) (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, Mixed Forest Land (including Forested Wetland)
2. Agricultural Land 7. Water
3. Rangeland 8. Land Falling Outside the Study Area
4. Deciduous Forest Land 9. Non-Forested Wetland
5. Coniferous Forest Land 10. 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):
ZQ (1) = 0.7 i**1*
(2) = 0.2
(3) = 0.1
(4) = 1.0
(5) = 1.0
49
-------
(6) =0.7 (*2)
(7) = 0.05
(8) = 0.00^*3)
(9) = 0.3
(10) = 0.15(*4)
Notes:
*1. Sheih, et al. (1979) recommended a value of Z0=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, tm) using
k u
u* =
G(z,L,zQ)
(4-6)
where k = 0.4 is the von Karman constant; z is the elevation of the surface
wind observations above ground—use
z=10m;
(4-7)
and G is the similarity function
G = In
o
if 1/1 = 0
(4-8a)
G = In
MIN(5.2f,5.2) , if z < 6L and L > 0
-1 -1 (M)(40+1)
G = 2(tan i | - tan i 4n) + In C/ui\/t -1\3* if L < °-
(4-8b)
(4-8c)
The Obukhov length L is from (4-4), ZQ is the local surface roughness from
(4-5), and
4=
Z+Z
40 = a-15-j)
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)
_p
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
Variable
Q (t )
Description Source
east-west surface _n P3
Step
1
Output
Variable
Ux,t )
Description
Obukhov lengl
:h (meters).
wind component (ms )
at surface weather
station n, hour t .
north-south surface P3
wind component at
station n, hour t .
location of surface PIF
weather station n
land use fraction PIF
of category j in
NEROS cell centered
at x.
aCT(x,t ) fractional sky PIF
coverage, total of ,
all cloud types over
cell centered at x
at hour t . ~
m
in NEROS cell centered at
x, at hour t .
~ m
iu (x,t )| w.ind speed (m/sec) in NEROS
~ cell centered at x, at hour
t (for use in thTs Proces-
sor only).
see above
PIF
effective surface
roughness (m) of
NEROS cell centered
at x.
u (x,t )|wind speed
m
L(x,tm) Obukhov length (m)
***
surface roughness
On)
Step 1 3
Step 1
Step 2
friction velocity
(m/s) in NEROS
cell centered at
x, at hour t
~
.
m
Tfr,(x,t ) surface virtual
w n ^* m . . »r* — •
vn
temperature (°C)
at surface weather
station n, hour t
m
P3
surface,heat flux
(m °Ks"1) in NEROS
cell at x, at hour
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*(x,t ) friction velocity Step 3
~ m (m/s) in cell at
x at hour tm.
~ m
L(x,t ) Obukhov length Step 1
~ ra (m) in cell at x
at hour t.
54
-------
YYYY
TJ
C
CO
3
a.
c
l
•*- s.
o *-
I!
S2 «
» 8
3
® a
.c <—
U 3
(/) O
O)
iZ
v V V V V V
55
-------
. • 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 Hys 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 sealevel. 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 1 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-^Pg-
Figure 7-1. Illustration of the variables used in the flow model.
The x-component of the horizontal acceleration of a (2-D) fluid parcel in
the cold layer is
d.u
(7-2)
where m is the mass of the parcel in question and Fxp, FXC> 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
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/8z = -pQg.
(7-3)
Figure 7-2. Air parcel considered in the force balance analysis.
58
-------
And since pQ is assumed to be constant, at any level 2 within the parcel the
pressure is
.P = Pvs + P09(zvs -2) (7-4)
where p is the pressure at the top of the parcel, that is on H .
V o V 5
The total force on the left face of the parcel due to the hydrostatic
pressure is
2vs
fXL = Ay f pdz'
2t
Substituting (7-4) into this expression and making use of the definition
h(x,y,t) = 2vs(x,y,t) : 2t(x,y) (7-5)
we get
Keeping in mind that the pressure force is exerted inward on all faces, the
force on the right face of the parcel is
3p
fXR = " (pvs + ~~^ A2vs)(h+Ah)Ay - o0gAy(h+Ah)z/2 (?.7)
82
The force exerted by the ground on the parcel is directed normal to the
lower face which is inclined at an angle 6T 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 et
= - (P0gh + PVS)(A*2 + A2t2)35 Ay sin et (7-8)
= - (Pogh + Pvs)Ay Azt
In a similar manner we find that the force on the top face of the parcel is
fxvs
On combining (7-6) - (7-9) and noting that 3p /3z = -p g we obtain the total
x-component of the pressure force on the parcel: f
s
2
V
Fxp = fXL + fXR * fxt + 'fxvs |
= Ay(A2vs - Azt - Ah)pvs + pmghAyAzys (7-10) |
- P0ghAy(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
^ = po'pm . (7-12)
Now the mass HI of the fluid parcel is simply
m = poghAxAy,
60
-------
hence
i F = - £& g vs . (7-l3a)
m xp p b 3x
By analogy the y-component of the pressure Induced acceleration is
3z
1 r AD „ VS
F-- 9 -—
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 face
of the parcel is
fvl = Ayht (7-14)
AL XX
61
-------
where T is the net flux of x-momentum in the x-direction across the left
^A
face of the parcel. The force on the right-hand face is
fXR =
c *!L2
xx ax
Ax)(h
8z
Ax)Ay
,»&*
(7-15)
Figure 7-3. Friction forces on a fluid parcel of horizontal
dimensions (Ax,Ay) bounded by zvs and z^.
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
62
-------
faces of the parcel are approximately AxAy, we get
Fxf = -
3l
- AX
32
8Z\
W
55s 5
atvx 9zv
_*X _X
(7-16)
In the limit as Ax and Ay become very small, the mass m of the parcel is
m = p hAxAy
and hence
F - - -
m xf pQ
(7-17)
8y h v zx
(7-18)
T +T
xx h 3x yx h 8y
For the lateral stress components i and T we will adopt the gradient
yv/\ Jf ^
transfer forms
P0Txx
P0V
§U
8x
(7-19)
(7-20)
where K and K are elements of the eddy vicosity tensor, to be defined
A^\ — y
later. We will treat the stress T^ on the upper surface of the fluid as a
fm1\
prescribed variable. The surface stress T will be given later.
63
-------
By analogy with (7-18) we have
i yf ~ p_ •• 3x ay ' h vtzy tzy>
0 (7-21)
* Txy h 3x + Tyy h 8y •*
with
-I T = - K f§H + QL]
p xy xy 3y ax-*
p7 V = "V ay'
The stress t^J is considered to be a prescribed variable; T* and T* are given
AY zy zx
by
Tzx = " P0U*COS 9 (7-25a)
Tzy s ~ P0u*s1n 6 (7'25b)
where u* is the friction velocity which we shall express in the form,
following Melgarejo and Deardorff (1974),
u* = CQ(u2 + v2)%; (7-26)
where C is the surface drag coefficient and 9 is the wind direction,
8 = tan"1 J. (7-27)
From (7-25) - (7-27) we obtain
p* Tzx = " CD (u2 * v2)35u (7"28a)
64
-------
r Tzv = - CD (lj2 + v2)3V (7"28b)
^ y
According to Melgarejo and Deardorff,
Cn = k2 [(In(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
1 Fxc = +fv (7-30)
I Fyc = -fu (7-31)
where f is the Coriolis parameter (=2Q sin 4>, where <)> = latitude and fl = earth
angular speed of rotation = 211/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
-------
sz 2 2
9u * ..3" . ,,3u _ Jte vs - .x 3fcu » 3 u
at + "37 + vay - - IT + f v * Kxx 5 + Kyx
8Kxx . Kxx 8zvSl au . r8Kyx . ^x 8zvs, 3u
TT * IT -aT] Sx * c~ay + "T "ay-3 ay
VS
K 3z T
_J2i vs §v zx
h ay ax " hpQ-"
For convenience we will express the shear stress at z in the form
p Tzx = ~aw (UM " u) (7"33)
where a^f is a vertical exchange velocity scale at z = z „ and Uu is the
w vs 1*1
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
§H + r,, xx I^x 9zvs1 3u . ru 8Kyx x ySl 3u
3t LU 3x h 3x J 3x LV 3y h 3y J 3y
(7-34)
Jt* azvs av, . i 8psi
h 3y 3xJ p ax
66
-------
Likewise we have
§v + [u . xy_ . ys, ay + [v . . vs 8v
3t LU ax h ax J ax LV ay. h ay J ay
? 3* VS
V G
r-q (£E) _JL» - fu + r a V.. + -^ — + K
Prt 9y h w M 3x ay xy
h 3x 3yJ p 3y
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 tQ be (XQ, yQ, ZQ). If the surface is moving
with a velocity vys = (uys, vvs> wys) at the point Q, then at the later time
t^ + At Q will have the coordinates
o
(xr yr zx) = (x0 + uvsAt, yo -H vvsAt, ZQ -H wysAt). (7-37)
67
-------
Since by definition Q lies on the surface HVS at all times we must have
H( = 2(x' " z = ° (7'38)
vso' o
' Hvs(xl* VAt) * 2vs(xl>yi>VAt) 'zl = ° (7"39)
If At is sufficiently small we can write
S 2vsCxo'VV + T (xl ' V
(yry0) + At
ay -1 ° at
where all derivatives are evaluated at (XQ, y , tfl). Substituting (7-40) into
(7-39), subtracting (7-38) from (7-39), and making use of (7-37) we get
vs VS VS
-IT = IT "vs + -af V * -5T - wvs
and upon taking the limit At*0 we obtain
The unit, outward normal vector to the surface H at Q is
-—l2Hu« (7-43)
where V is applied at Q. Let As denote the area of the projection on H of a
horizontal rectangle (Ax,Ay) centered at the coordinates (*0,y0) °f Q (see
Figure 7-4). It is evident from the figure that for sufficiently small Ax and
Ay, the area As is
68
-------
As *
ayj
(7-44)
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
az
—
3x
vs
«i
(7-45)
then (7-44) can be approximated by
3z 2 3z 2 ,
As ~ AxAy [(-a~) + C-sr5) + 1] . if 7-45 holds.
(7-46)
Let vf be the velocity of the fluid at the point Q on Hvg. Then the
normal downward component of the fluid velocity relative to HVS at point Q is
= ' (Jtf - vus) • n
fv «VH - v-'VH ]
iu~vs ~ vs ~f ~ vsj
(7-47)
69
-------
Making use of (7-42) we can write (7-47) in the equivalent form
"HOT? . .. <™>
where
J* = _§. + w .n c,
Ht ~ At «!*•? ~ *
ui» a u »»| <«r
The total downward volume flux of fluid crossing As is vAs. Using
(7-46), (7-48) and
32... 2 8z 2
we have
dHvs
vAs = AxAy - (7-51)
Thus, the net downward flux or fluid per unit horizontal area is
n «^S (7-52)
nvs ^JT c/ b^J
or
where (u, v, w) is the fluid velocity on H
-~T-~ - — V5
An expression for n, 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., radiative cooling on H and heat transfer to
the terrain surface Ht, we can write the first law in the form
70
-------
d6 _ 86 . ,86 . 86 , 86 _
dt - 8t * "Bt * V8y + "§! *
where 6 is the potential temperature defined as
(7-54)
e = T(-
P
P R/CD
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)
dH,
vs
(7-56)
vs
WH.
~ ~ t
— vs
= 0
t vs
where T~J denotes averaging over the terrain surface H., T~7 denotes averaging
over the virtual surface H ,
and
n(x,y,t) = zvs(x,y,t) - zt(x,y)
5 H f f <*2da
(7-57)
(7-58)
and A is the area of a grid cell in the model
Look at the first term in brackets in (7-56). In analogy with (7-52) we
have
dH,
(7-59)
71
-------
where nt is the downward fluid volume flux at the ground. Thus,
t
di-L
= -^ (7'60)
where w'9' is the kinematic heat flux at the ground.
Consider now the second term in brackets in (7-56). Since 8zt/3t = 0
__t
and since Ht = 0 (there is no net flux of air through the terrain), we have
t
t
= o. (7-61)
The third term in brackets in (7-56) can be written with the aid of
(7-52) as
vs
vs
and the fourth term becomes, using (7-52) and (7-53),
vs 32
(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
- d , _^ vs
~ <0> + i r-w'61 - 6n + <6>n ] = 0 (7-64)
at h o vs vs
where
d,
(7-65)
72
-------
Following Zeman (1979) we shall assume that within the cold fluid, i.e.,
-r.e nighttime stable boundary layer, the potential temperature has a linear
ariation with height, namely
e = eh - (eh - e0)(i - ^) (7-66)
/here 8. is the value of 9 at the elevation h above ground and 9 is the
potential temperature at the surface. At elevation h the turbulent heat and
Tiomentum 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
due to the motion of the surface itself. In this case we have
vs
if%s>0
Also, on integrating (7-66) in the manner of (7-58) we get
<9> = h (9h + 9Q) (7-68)
and hence
d,, du du
n ^rt^ _i_nA.inrt r-
dt ^ dt wh ^ dt o*
Substituting (7-67) and (7-68) into (7-64) we obtain
"« Veodt 9h-8o
Making use of (7-69) in this equation we obtain
Hvs = B(he-h) (7-71)
where
d
+ V (7'72a)
73
-------
2wTern
h = - - £- . (7-72b)
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
= B'(h-h) (7-74)
where
• /•* o^
pi - _ 4/3 _ 0
~ ® 3t
ne - c4 resin cosa (?.75b)
where f is the Coriolis parameter, G is the geos trophic wind speed, a is the
angle between the geos trophic 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-71j - (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
jSrgr j s CUL5 fG2sina cosa^ (?.76)
To determine whether this is a plausible relationship, we note first that
for large values of h/L, where L is the Obukhov length
L = - — y* - (7-77)
gk vT^e
Brost and Wyngaard (1978) found
G2sinacosa = u£ ()2 (^) (7-78)
where k = 0.4 in the von Kartnan constant. Substituting (7-77) and (7-78) into
(7-76) we find that (7-76) implies
h = 0.8 (O (7-79)
But this is just the formula that Zilitinkevich (1972) derived for h from
u I i,
similarity theory [in particular h = c5(-4=) ]; and thus we conclude that
(7-76) is a plausible expression, at least for large h/L.
Having an expression now for n. , 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
t VS
with respect to z, we get
75
-------
"C+• "t - "„• '
Since there is no flux F of fluid through the ground, we see from (7-61)
that
az. azt BZ+
5T + "5F * »5T - wt ' ° (7'82)
Solving (7-82) for w. and using the result in (7-81) we obtain
»«• 4r+ v§r - hC§J+1] (7'83)
where h is given by (7-57). Substituting (7-83) into (7-53) we obtain finally
§£ + ,lll * »?Jl j. Kr9." j. 9Yi = « (7-84a)
at
where
h = zvs - zt (7-84b)
This is the model equation that we shall use to obtain zys.
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
o^s~ 0. (7-85)
Let us also assume in this "first generation" model that the lateral eddy
76
-------
.iscosity tensor K is also zero. With these assumptions (7-34) and (7-35)
-educe to
2 /, rt ! 3z
at ' "ax -ay H H CD (u +v >* = '* pn ~*T
1 3psl
pj IT * fv (7-86)
9v * ,^v . W8v . v r2 ,2 . ,,2^ _ „ Ap 3z
at + "a* * Va7 + H CD (u + v } "g'S
"aPsi
: -*~- - f" (7-87)
where p , is the sea-level pressure and h and z are given by (7-84), which
91 V o
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 , , Zj and nvs- Earlier we described how n
is related to the surface heat flux w'9' and the temperature distribution 6
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) n - Eq- 7-71,72 with 6. and 6 from surface and upper air mete-
orological data (Processors PI and P3); w'6' = Q = kinematic
heat flux at the ground (from Processor P4);
(2) Zy(x,y) = 30 x 45 min (lat-lon) smoothed terrain heights (meters,
MSL, from P7)
(3) Cn = k{[ln(-£) -b]2 + a2}"35 (7-89a)
2o
where k = 0.4 is the von Karman constant; h is the local
solution of Eq. (7-88); ZQ is the surface roughness (m), from
P4;
f5, if L>0,
/_-!, otherwise;
a
(7-89b)
•T
(.4,
If L>0,
(7-89c)
'otherwise;
and L is the Obukhov length (from P4). Note that L, ZQ, 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)
where 6k and 9 are from Processors PI and P3;
ho
(5) pQ = (7-90b)
where pgl is the sea-level pressure (from P3); Ty is the virtual
temperature (from P3); and R is the gas constant.
78
-------
_
(6) g and f are gravity (9.8 m sec ) and the Coriolis parameter (=
2Qsin<|0, respectively. H
local latitude in radians.
2Qsin<|0, respectively. Here Q = 2n/(24-60'60) sec" and is the
.ater we outline the stages of calculations that are needed to compute the
parameters listed above and the additional quantities, not used in the flow
i
.simulation, that P7 must provide to the model and to other processors in the j
network. \
Solution of the u, v, and z equations
Each of the Eqs. (7-86) - (7-88), which govern u, v and zvs>
respectively, has the form
i/|y + br = 5 u~si; t
!
ii
where U and V are functions of r. 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 F at the beginning of the interval, t say. In this case
the solution of (7-91) can be expressed in the closed form
F(x,t0+At) =Jp(x,t0+At|x',t0)F(x;t0)dx'
(7-92)
VAt
p(x,t+At|x;t')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,tjxlyit0) = 6[x-(xlfx)]5[y-Cyl+y)]e~b(t~to) (7'93)
where 6[x] is the delta function and
t
x = x(t|x',y',t0) = U[x' + x(t'), y' + y(t'),t']dt' (7-94a)
= y(tjx',y',t0) = Vfx1 +x(t'),y' + y(t'),t' ]df . (7-94b)
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
if] TT At <
-------
F-
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
t
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.
Stage 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 zt(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 (I,J) (column I, row J) of the NEROS region. Then
zt(i,j) (7-97)
) L
where the summation is over the 54, 5 min x 5 min cells surrounding grid point
(I,J) (see figure 7-5). Note that the 54 cell smoothing area used in the
definition of zt(I,J) overlaps the smoothing areas associated with the 8 NEROS
grid points nearest the point (I,J).
82
-------
The input and output requirements of Stage ZT are summarized in Table 7-1
later in this section.
®..
®_
1
(I.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 (7-98)
1 n * at time tn;
.0, otherwise.
Under this definition the magnitude of
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
5Pl(tn) =
+1, if Ap
-1, if Ap^t
0, otherwise
) = 0 and Q(tn) < 0;
= 1 and TDEF(tn) <0;
H(tn)-nQCi,J,tn>
(7-99)
is the surface kinematic heat flux (°Kmsec ) summed over all grid points (i,j)
of the model domain, and TOEF is the mean "temperature deficit" of the cold layer
which we approximate by
2
max
m
(7-100)
In this expression hmav is the depth of the cold layer at the time the surface
UlClA '
heat flux reverses, i.e.,
h(i,j,tm), if Q(i,j,tra) > 0
< 0;
(7-101)
0, otherwise
The value assigned to h before the surface heat flux reversal, in this case
max
0, is immaterial since the value is never used. In (7*100), the variable q is
given by
•{
Q(U,tm), if Q(i,j,tm) > 0;
0, otherwise
(7-102)
84
-------
The constant a that appears in (7-100) is a temperature gradient that we
Define in Stage ETA.
(2) With 6p}(tn) determined by (7-99) a functional form for Ap1(tfl) that
s consistent with (7-98) is the following:
Apl(tn) = Apl(tn-l) + 6pl(V (7-103a)
v/ith initial value
Ap1(tQ) = 0 (7-103b)
and
t0 = 1200 EST. (7-lOSc)
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,(tn),
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 .FLOMQD, 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/pQ, 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 8h from the observed
ground-level temperatures 6 assuming a constant temperature lapse f
1
Q/3
|f = or (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
6h = eo * ffh (7-106)
with cr given by (7-105); and hence from (7-90a)
o o
We will estimate the surface temperature 6 at each grid point (I,J) using the
virtual temperature observations T , n=l,...N made at the N surface weatl-.sr
stations. We assume here that the T (°C) are available from Processor P3 for
each hour. In this case 9 (I,J,t ) is computed as follows:
86
-------
(7-108)
n=l
where
rn - <7-109)
is the distance between grid cell (I,J) and the site (x ,y ) of surface
station n. Equation (7-108) should be evaluated at each time step tm that
The inversion layer growth rate parameter nvs 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 rj 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
nvs(I,J,tm) = max {0, B[he - h(I,J,tffl)]} (7-110)
where
87
-------
and
2Q(I,J,tJ
In these expressions we can use the approximations
t
(7
8h
3t z At
where 8Q 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
hg. 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 or acquired from wind observations
would probably be erroneous.
In summary. Stage ETA provides values of Ap/pQ (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 6Q 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,(t ) = 1-
We assume that the sea-level pressure data p , measured at each of the
5 i y 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:
n=l "
where
rn = [(IAX - xn)2 + (JAy - yj2]* (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
n=l
where r is given by (7-115) and pn is the density at surface station n (from
P3).
89
-------
At each grid point we require values of the functions
. 3p .
Px .5! -jfl (7-117.)
"y ' p - (
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'tm) =
- Psl(I-l,J,tm)] - l/12[psl(H-2,J,tm) (7-118)
.) « A[p0(I,J,tm)]"1{2/3[Psl(I,J-H,tm)
- Psl(I,J-l,tra)] - l/12[psl(I,J+2,tm) (7-119)
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
PQ 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 j.
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
-------
h(I,J,t ) j 2 -Sj
=k{[ln ( )-b]Z +a2} ' (7-120)
iere k=0.4; z is the surface roughness (m) of cell (I,J); h is the depth (m)
f 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 CD(I,J,tm) 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)
2_2
^-j- = ~PX + fs (sine) (7-121)
9
= -p - fS (cos6) (7-122)
91
-------
where s = (u2 + v2)^ is the flow speed and 6 is its direction. Cross-
multiplying (7-121) by (7-122) we get
P sine -. P cose = fs (7-123)
x y
and on squaring (7-121) and (7-122), adding the results and jnaking use of
(7-123) we obtain
s4C 4
-^-8- + f2s2 - (Pv2 + P 2) = 0. (7-124)
h^ 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 e.
Next, we define the time t«g:
t«Q = 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^Q, Stage IBC computes initial values as follows:
h(i,j,tK£)) = h0 » 30m, (7-126)
u(i,j,tKO) = s(i,j,tKO)cose(i,j,tKO) (7-127)
v(i,j,tKO) = s(i,j,tKO)sine(i,j,tKO) (7-128)
92
-------
here s(i,j,tKO) is the solution of (7-124) at grid point (i,j) based on
alues of PX, P , CQ and h (from 7-126) at that point; and 6(i,j,tKO) is the
:orresponding solution of Eq. (7-123).
if t™ > ti/n and Ap-,(t,) = 1, Stage IBC computes
m _ t\u j. in
boundary values of u, v, and h as follows:
flBC^'J*V = %s(1,J.tB) (7-129)
flBC<1'J'V = S£ tu(i>J,tm) - ud.j.t^)] (7-130)
*Bc(1,»J'V = & Cv(i>J''V ' v(1>J>Vi)]- (7-131>
In these equations (i,j) are boundary points only. Also, nvs is given by
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 Ap1 = 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.
h1(i,j,tm)=
h(i,j,t ), if Ap,(t ) = 1;
m i m (7-132)
.max{100, min{500, L(i,j,tm) }},
if Ap^tJ - 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 zys is computed as follows:
z-rd.j) + h^(i,j,tm), if Ap, (tm) = 1;
Since Zy is given in meters above sea-level, the values of z^ 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:
z-,(i,j,t )g ...
•Pl^'J'V = Psl(1.J,tn)exp[-i —] (7-134)
RdT(1,J,V
where g = 9.8m sec" , Rd = 287m sec"2 °K , Z-^ = ZT + hp and
T(1,j,tB) = %[2eo(i,j,tB) + 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. (7rl34) with z1 replaced by zys (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
h0(i,j,tm) = h^i.j.y/lO. (7-137)
The inputs and outputs of Stage H1HO are summarized in Table 7-1.
Stage SIG
This stage estimates the fractions a™ and QJQ of surfaces 2i(x>y.tm) ancl
[Zy(x,y) + h0(x,y,tm)] that are penetrated by terrain.
Referring to Fig. 7-5 which illustrates the process of computing Zj in
Stage ZT, we define
95
-------
ov,(I,J,t ) = 4 £*0',j,t ) (7-138)
11 m 54 1,jeD(IfJ)n ,
where X is defined as
f l,if zf(i,j) > z,(I,J,tm); (7-139)
A / ^ * ^ \ — ^ ™ "*
&\ • »J j » y " '
m (_ 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
where a^ 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 97 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
surface kinematic ,,
heat flux (°K m s"i;
at grid cell (i,j)
at hour t
P4 DELRO Ap,(t )
L
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
-------
Table 7-1. (Continued)
Input
Variable
Description
Output
Source Stage . Variable
Description
h(I,J,tm)
depth (meter) of
cold surface layer
in cell (I,J) at
hour m.
Stage PCD Cn(I,J,t )
FLOMOD (cont.) u m
drag coefficient
(dimensionless)
in cell (I,J) at
hour m.
z0(I,J)
L(I,J,tB)
surface roughness
(m) in cell (I,J)
Obukhov length (m)
in cell (I,J) at
hour m.
P4
P4
p , (i,j,t ) sea-level pressure ii
cell (i,j) at hour t
•yUlJ.«.BJ
horizontal pressure Stage IBC
force term (see 7-118) PCD
horizontal pressure Stage
force term (see 7-119). PCD
h(i,j,tKQ)
initial depth (m)
of the cold layer in
cell (i,j)
initial east-west
flow speed (m sec
component in cell
_l
CD(i,j,tm) drag coefficient
.-.-O'j.t ) cold layer growth
vs m rate (m/sec)
Stage
PCD
Stage
ETA
v(i,j,tKQ)
m
lf V
W'-J>V
initial north-south
flow speed component
in cell (i,j)
time derivative (m/sec
of cold layer depth
at boundary cell
(i,j) at time tffl
time derivative of
the east-west flow
component on the
boundary at time t
m
time derivative of
the north-south flow
component at boundary
point «(i,j) (units:
m sec )
Zj(isj) elevation (m, MSL)
of smoothed terrain
in grid cell (i,j)
Stage H1HO
ZT
h-i(i»j»t ) thickness (m) of
1 m Layer 1 in cell (i,j)
at time t.
m
98
-------
Table 7-1. (Continued)
Input
Variable
Description
Source Stage
Output
Variable
Description
hcu.y
depth (m) of cold
surface layer
sl
m
L(i,j,tm) Obukhov length (m)
Pci(1»J»tm) sea-level pressure
(mb) in cell (i,j)
at hour t
m
surface temperature
(°K) in cell (i,j),
hour t
m
Stage H1HO zu,(i,j,tm)
FLOMOD (cont.) V5 m
P4
Stage
PCD
Stage
ETA
elevation (m, MSL)
of virtual surface
in cell (i,j), hour
V
elevation of surface
z, in pressure (mb)
coordinates.
elevation of the
virtual surface z
in pressure (mb)
coordinates
indicator of
presence of cold
surface layer
Stage
DELRO
depth (m) of Layer
0 in cell (i,j) at
hour t .
m
zt(i,j) smoothed terrain Stage SIG
elevations (m, MSL) • ZT
zt(i,j) 5 min x 5 min averaged RAW
terrain elevation
(m, MSL) of cell (i,j)
z-,(i,j,t_) elevation (m, MSL) Stage
1 ra of top surface of H1HO
Layer 1.
fraction (0<0y-. . cell averaged east-
west flow speed
component (m/sec)
at time t in the
cold laye?1
99
-------
Table 7.1. (Concluded)
Input
Variable
Description
Output
Source Stage Variable .
Description
CD(i,j,t_) drag coefficient
Ap/P00'.j»t )buoyancy parameter
Stage
PCD
Stage
ETA
vL
except nortn-so
component
nwe(1»j»t ) cold layer growth
vs m rate
Ap1(t ) indicator of presence
of cold surface layer
h(i,j,tKn) i
u(i,j,t^)| initial values of
v(i,j,t.,n) i h, u, and v
Stage
ETA
Stage
DELRO
Stage
IBC
,.
m
time rate of ,
change (m sec )
of cold layer
depth at boundary
point (ijj1) at
time t
m
"2
acceleration (m sec
of east-west flow
component at boundary
point (iij1) at time
)
same as UBC except
north-south flow
component.
Stage
IBC
Stage
IBC
Stage
IBC
100
-------
Appendix A to Section 7.
Here we describe a numerical procedure for deriving solutions of tf
differential equations of the form it
'5;
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
r(x,t.) = p(x,t, | xl
~' ly rv~' 1 ' ~»
*0
where
t
b(t")dt"] (7-A3)
t1
t
•
x = x(t | xjyjt1) = | UCx'+xCfhy'+yCt'^.t'^dt" (7-A4)
t1
y = yCtlx'.yU1) = IvCx'+xCfhy'+yCfht'^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
x = (lAx,JAy) , I,J=1. ..
and only at the discrete time intervals t =nAt, n=l ..... Furthermore, in the
• n
situations of concern to us the spatial variations in S are of a scale much
larger than UMAXAt or V^At and the temporal variations are generally
slower than the time step At. Under these conditions we can express (7-A2) in
the approximate form
p(lAx,JAy,tn+1l xltn)F(xJtn)dx' (7-A7)
J
where
rjj = r(!Ax,JAy,tn) ' (7-A8)
and
F^itn) = r(xjtn) + S(x;tn)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',tn) 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
x1 ss x* = lAx-x
y1 - y* = JAy-y \
102
(7-A10)
-------
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 x1 = (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* = lAx - U(lAx,JAy,tn+1)At 1 (?_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
" n
define the new coordinates
n - (7-A14)
Note that the origin of the (n,£) system is the cell center (1ST,JST). We can
103
-------
now express F(x' ,tn) by the biquintic (5th order) Lagrange polynomial
6
TT
6 '
(7-A16)
O G G G
O G O O
o o o o/ o o
• • •
0
o •
G •
G G
G G G
G G G G G G
/ Back trajectory starting
from (I,J) at time t +,,
going back in time
totn.
GRID POINT (I,J)
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
(7-A17)
:Jj = F(iAx/2+IST,jAy/2+JST,tn)
(7-A18)
104
-------
and
6 ,
Fl ( ) = product over all k except k=i,
k=l .
'
J] ( ) = product over all 1 except l=j.
1=1
Note that (a. ,b.)> i,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, =11 (n-ak) (7-A19)
1 k=l K
'a-bn) (7-A20)
' I
A- '
L, -I! (arak) C
1 k=l 1 k
M, =n'(brbi)- (7-A22)
J 1=1 J '
The last two parameters can be evaluated directly using (7-A17). We get
o
L2 = M2 = 768 '
L3 = M3 = - 384 I (7-A23)
L4=M4= 384 r
L5 = M5 = - 768
= 3840
_J
105
-------
Now (7-A16) can be written
6 6 n X4Y,
F(x',tn)=2; I F.. ^ (7-A24)
1=1 j=l J i j
Now look at the expansion of X-:
Xj = n,5 - 5n.4 - lOn.3 + 50n2 -t-9n-45 (7-A25)
X2 = n5 - 3q4 - 26n.3 + 78n2 + 25n - 75 (7-A26)
•
Xg = n5 + 5n4 - lOn.3 - 50q2 + 9n + 45. (7-A27)
Thus, we can express the X- in the alternate form
5
X1 =2 a^n* (7-A28)
where the coefficients a-a are the known values given in (7-A25 - A27).
Similarly
Yj =2- ajp|p. (7-A29)
Substituting (7-A28) and (7-A29) into (7-A24) and we obtain
5 5
106
-------
6 6
n -o - n a.
where
A =2- 2, F..-12JE . (7-A31)
«P i=l j=l !J L.H.
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,4)« We get
oo
-08
CI.*iVl'WJ'y''Vl'V5 At (7-A32)
F(n,4,tn)dnd4
where
(7-A33)
= _JL exp [- (Jtnr'V)1 3 (7-A34)
2cz
(see Eq. 7-A3); and b* = b(x*fy*,tn).
Let
ex = x* - 1ST (7-A35)
Ey = y* - JST. (7-A36)
Then, representing x1 in (7-A33) in terms of n we obtain (from 7-A14 and
7-A10)
<*' exp [- -]
107
-------
where
ex = ey/(Ax/2)
/ (7-A38)
a = a/(Ax/2)
(7-A39)
An expression similar to (7-A37) describes * .
Substituting (7-A30), (7-A37) and the analogous expression for , 1nto
(7-A32) we obtain *
or=0 ^=0 "P ore (7-A40)
where
(7-A41,
v = i f.R r -" -v
" 4 exp [- Ji_ T dj
J 2&2 J d^ (7-A42)
-08
These last two integrals can be evaluated analytically
and we get
ex
(7-A43)
X
108
-------
Similar expressions give v , 8=1,...6, except e 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" is
derived from (7-A31) using known values of F at time level n; and where the X
and v parameters are given by (7-A43).
109 i
-------
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
it *•$!+ *l *ii c§ w***-
(7-Bl)
3v . 3v , 3v . v
—~- + UK— + Vs— + T-
8t ^x 9y n
(7-B2)
110
-------
j
if
where
D_ D
px - p§! ' Py
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 (u1, v1, h') 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
* * = * ' * -' * "'
Now let u and v be solutions of
= -Pv + fv (7-B7a)
n u x
111
-------
= -P - fu. (7-B7b)
and assume-
. 32u] (7-B8)
ax
8u _ 8v ~ n ~ J
it " at " u " l
Subtract (7-B7a) from (7-B6) and use (7-B8):
From (7-B2) and (7-B5b) we have
^zt - »2u
I. - p fu . fu, . K[|_v
8V,
II I
2 J
112
(7-B9)
-------
Using (7-B7b) and (7-B8) we can reduce (7-B10) to
Now combine (7-B5c) and (7-B3):
8h 8h' . -ah . -ah1 . nl3h . Itl8h'
"ix + U3*+ U9^
ay 8y ay
I £i>'
(7-B12)
ay
where
(7-B13)
and A is the model domain. We define
iH = n. (7-B14)
113
-------
Subtracting this equation from (7-B12) we have
(7-B15)
where we have assumed
-- 0. (7-B16)
Now we write the basic equations (7-B9), (7-B11) and (7-B15) in the following forms:
«i) » + Su + Si + K(§3jSf +ra) (7-B17)
= + Sv + s; + K(+ ) (7-B18)
Sfc * Kh(|^' * ) (7-B19)
where
pu = 8u/3x (7-B20)
- Pv s 9v/3y (7-B21)
p = 8u/8x + 8v/8y (7-B22)
114
-------
]| 1f u. t Q;
u n h
. . (7-B23)
C2 _ ,
,-fi(u2+v2) , otherwise
C2
TTT ErC^+v2)*5 - * (u2-^2)*5] , if v1 * 0;
v n h (7-B24)
6' =
HV
C2
-!i(u2+v2)^ , otherwise.
h
9u' , 9v'
" (7-B26)
(7"B29)
(7-830)
(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=IAX
y=JAy
A („ \ - 9.3
Ay(gI,J)
^y
<=IAX
y=JAy
We define several different AX and A operators, each of different orders of
accuracy.
(7-B33)
(7-B34)
(7-B35)
glfj.2)3/(12«y) (7-B36)
9IjJ.1)/(26y) (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
-------
,J * 169gI+l,J + 105gI+2,J ' 160gI+3,J
)J + gI+6jJ]/(1205x) (7-B38)
= C'171gI,J + 16"l,J+l + 105gI,J+2 - 160gI,J+3
+ 65gI,J+4 - 9gI>J+5 + gM+6]/(120Sy) (7-B39)
[171glfj
(7-B40)
glfj.6]/(1206y) (7-B41)
The operators defined above are used to compute the variables p , jj , ... S1
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, pv, jj^, S , Sy and S^
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 zt, cooling n'» fitc. 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, v', 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 = h' =
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'(l,J) = v'(2,J) = v'(3,J) = 0 } if [u(4,J) + u'(4,J)]> 0
h'U.J) = h'(2,J) = h'(3,J) = 0j
At points of outflow, we use the following extrapolation, illustrated for the
case of point (4,J):
118
-------
READ FIXED
FIELDS
zt.f.
INITIALIZE
u'=v'=h'=0
READ
u. v, h. CD. g'. n. n'
*
r
COMPUTE
BAR_ VARIABLES
/Jv. 3u. P_h
Su. Sy. Sh
'
r
COMPUTE
PRIME VARIABLES
/J'u» P'v> P'h*
S'u. S'v. S'h.
\
r
CALL SOLVER
COMPUTE NEW
u'. V. h'
N
to
c
o
*•*
m
I
o
Q
O
O
CO
u
CO
3
O>
119
-------
u'Cl.J) = u'(2,J) = u'(3,J) = 4U H
v'(l,J) * v'(2,J).= v'(3,J) =4V
h'Cl.J) = h'(2,J) = h'(3,J) = 4n
if [u(4,J) + u'(4,J)]< 0
where
= - 2u'(5,J) + (l/2)u'(6,J) + (5/2)11'(4.J)
etc.
Calculation of the BAR variables.
Figure 7B-2 shows the grid on which the BAR variables p , S , etc. are to
be evaluated and the derivative operators that are to be employed at each
point.
From (7-B20) - (7-B22)
PhU,J) » Pud,J)
(7-B42)
(7-B43)
(7-B44)
and from (7-B26) - (7-B28)
120
-------
(7-B45)
SV(I,J) = G
(7-B46)
(7-B47)
in evaluating the six expressions above, the operators A and A should be
selected as follows (see Fig. 7B-2):
, if 1=66;
, if 1=2 or 1=65;
, if 1=3 or 1=64;
, if 1=4-63.
, if 0=1;
, if 0=48;
, if J=2 or 0=47
, if 0=3 or 0=46
, if 0=4-45
here 2A , 2AV, etc. are defined by (7-B32) - (7-B41),
x y
(7-B48)
(7-B49)
121
-------
H»'
17 1
rtfi
AC
j
44,
tf
;
5*
4,
«r
*
*x
2'
f
»
'
k '
1 J
i i
4/
,
• i
i 1
1 1
>
_/
e'
^x
»
» 4
' '
k
>x-:-x^^X:x".vX::::;::::X:::::X:::::::X:X:
;:;:g;:;i:;:;:x:;:i:;x::;:;x;:;:;:ix;x;:;:;:;:;:;:;:;
X:X:X:::::X:::::X:::::::::X:X:::::X::;:::::;:::::::::::
; * «.
•
K
« -r *
"^•'•"•'•'•^'••-'•'•'•'•^.•r'-v.xvX'X^'X'X'X'X*
'
r
« » eA.
4Ay
0"U r+^4 yy.*^*^*r1*?*^*.*?^*
I
;.;*X«|iX'
(11
• ::*X::
ill
«
1
•'.•'.•'•'. :i
::::::::::*
I'X'I-"*
•X X-ll
:::;|:::::i
.;.;. w *
-»»••
et
r
.
r •»«•»•
***•••
1
i
. . I •
I
I
'
*
(
I
I
I
\A
60 61 62 63 64 65 66
Figure 7B-2. Grid network on which p , p , p., S , Sy, and S.
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
-------
Pid.J) = PV(I>J> = Ph(I*J) = su(I»J) = sv(I>J) = sh(I>J) = °' (7"B50)
if I = 1,2,3,64,65, or 66; or 0=1,2,3,46,47, or 48.
all other grid points, namely 1=4-63, J=4-45 compute the PRIME variables as
llows:
where
i j * v! j> " -
'I,J I,J ' ' hj
if u' * 0
(7-B51)
(7-B52)
C2 v v
VI,J hI,J ' ' hI}J
where
(7-B53)
(7-BS4)
(7-B55)
123
-------
In equations (7-B53,. . ,7-B56) the derivative operators A and A should be
X }r _
selected as follows when they operate on PRIME variables:
EAX if 1=63;
\ 2Ax if I=5 or I=62; (7'B57)
4AX if 1=6 or 1=61;
( _ 6AX if 1=7-60.
if J=45;
if 0=5 or J=44; (7-B58)
if J=6 or J=43;
if J=7-42.
In Eqs. (7-B53,..,7-656) use gAx and 6\ in all 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,£6x5,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.u1 and v' (see 7-A12, A13); 4gxg 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 h' 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,e6x6,K)
where
Step 3.
u'(I,J,N-H) = C*(I.J)*exp(-B*(IfJ))
2. v'(I,J,N+l) in 60 x 42 model domain.
Same steps as in u1 calculation except replace Pu> B^, u'(N), Su and
su by P«» Pi» V'(N)> "SM and Si> respectively.
125
-------
3. h'(I,J,N+l) in 60 x 42 model domain.
Same steps as u', except replace pu, p^, u'(N), $u, and S^ by ph,
Pi, h'(N), S. 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) -i- 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-H) 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-H) ? u'(2,J,N-H) = 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) + (V2)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-H) = u'(65,J,N+l) = u'(66,JfN+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 h1. )
6. u'CN+l), v'(N-KL) 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) H- v'(I,4,N)) < 0, then
-2u'(I,5,N)
Similar expressions hold for v1 and h1.
7. u'CN+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 h'.
127
-------
If (v(I,45,N) + v'(I,45,N)) > 0, then
= u'(I,47,Nn) = u'(I,46,
(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,J,N+l) = l/2[u'(4,l,N+l) + u'CM.N-KL)]
1=1,2,3 and J=l,2,3.
and similarly for v1 and h1. We use a similar method to compute the dependent
variables in the northwest, northeast, and southeast corner zones.
128
-------
}
?J .
t WE
j BOUN
6» •
i ZO
•i -1
• i
ST
DARY
•
UE
•
41— JL—
3 S&igSS*
tt-im
s SOUTHWEST*
&? CORNERS*
2 <&&ffMSSSS£&i&t'.
••::'-':X' 7nwc ::'-v">
&£•& tUIHC v.:-.-.-.
, . 4 (9GRI.pPpl.NTS>
l«1 2
3
<
1
4
h' •.
*, !•!•!•'
•^»-
•:«•:-:•:•:•:•:•:»:•:•:•:•:•:•:•:»:•:•:•:••• •:•:•:•:•: •
S.MODEL DOMAINS x::;:.;::
:»:•:•.•: :•: :•:•:•:•:•:•: :•:•:•:•:•:•:-.•.•:•••:• •:• •:•'•••:
• • •
SOUTH BOUNDARY ZONE
• • •
567
Figure 7B-3. Illustration of the southwest corner zone and the
west and south boundary zones of the model domain.
129
-------
5YYYYYYYYYYYY
•-
si
O. 0)
u. "
O g
8 S.
o ~
O <»
7 O
2 «
11
If
e o
£ -o
o c
(A (0
CO
I
r-
S
I-
D
Q.
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 Ha. 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 bases, 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 z3 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 Hs 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
I-CT
.
V
(8-2)
Here w is the entrainment velocity (represented in Part 1 by f8/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); vj,. is the
horizontal wind velocity at elevation z2; and ac is the fractional area
covered by convective (cumul.us) 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 crc(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 31nVi A
at <^j + j "lit1 + SJH-SH^J * V I Fj-i. j - Fj,j3 = ° ^
J
where . is the cell averaged mixing ratio in layer j, A is the horizontal
w
area of a cell, V. is the volume of a cell lying in layer j and F. v is the
J J > K
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.,,. Our interest
is in 3 during periods when cumulus clouds are present (i.e., a / 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
we)(qc-3) + 3Z2 - (q* - 3) H3 (8-4)
-sZ3l = 0
133
-------
where z2 = 3z2/9t, z3 = 3z3/3t and
H3 = 23 + V3'2z3• - wa. (8-5)
Also, q represents the mixing ratio in cumulus clouds and a is the mixing
V* ^^
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
8 1 W2 -We
C (8-6)
•
and after multiplying by (z3-z2) and rearranging terms we obtain
Q-O" *M ' Q O
c c fln« _ c c
l-"c C " °C"° "' l~°e (8-7)
•
'R^ZS * ^(Q^ ~ s)
where
M3(x,t) = q(x,z,t)dz . (8-8)
z2(x,t)
and
(8-9)
At each of the sites x^ of rawin stations we have available measurements
~m
of q(z), V2u(2) and 8(z) at each observation time. We also know o 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 w and w , 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
-ith (8-2) to obtain estimates of z2, w and w throughout the modeling
lomain.
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 tx
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
l£/- wc - we(l-ac)3 - weqcac]dt = M3(ti) - M3(t0)
to
(8-10)
tl q a.
+ [w2(3 * p^T ) • Vs + w3(qOB-3)]dt
where all variables are evaluated at a given rawin site x . We assume that
~m
only w and w are unknowns in this equation. For all other parameters we can
C c
use the following expressions.
First, we adopt the approximation
qc ~ q(z2). (8-11)
Then
qc(t) * q(z2,t0) * qc(t-to) (8-12)
135
-------
where
qc = Cq(z8,ti) - q(z2,to)](t1-t0)"1 (8-13)
and
q(z2,t0) = q(xm,z=z2(to),t0).
Next, we approximate z2 and Zs by
= F(q(2B,zftn)> e(xm,z,tn)) iFO.l (8-14)
where F is a function that is defined later; and
z2(t ) + 100 meters, if a (x _,t )=0;
C m n n=0,l (8-15)
from satellite data, otherwise
Z3(t) = Cz3(ti) - za(t0)] (trto)"1. (8-16)
Since (8-4) bears the implicit assumption that HS > 0, which we adopted
in anticipation of application to situations where convective clouds are
present and growing, we must ensure that the value of £3 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 tj, but not both.
In this instance we can adjust the estimate of z2 at the hour that a = 0
to achieve is > 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 q^ by
q. - q(z3), . (8-18)
and from this we assume in analogy with the form we adopted for q , i.e., Eq.
8-12,
q.Ct) = q(z3,t0) + qjt-to) (8-19a)
where
q. = CqUa.ti) - q(z3,t0)] (trto)"1. (8-19b)
From (8-8) we get
M3(xffl,t1) = q(xm,z,t1)dz (8-20)
and
q(xffl,z,t0)dz (8-21)
Similarly,
3 = 3 + 3(t-t0) (8-22)
where
3 = [3 - 3]/(ti-t0) (8-23)
and
Ma(tn)
The vertical air speeds w2 and w3 that enter into (8-10) will be
approximated by
137
-------
w2(xffl,t) = wCx^CtO.t) , to < t < tu (8-25)
wsQ^.t) = w(xm,z3(t),t) , t0 < 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 z2(t) and z3(t) at which
w(t) is evaluated in (8-25) and (8-26) are taken to be
*2(t) = Z2(xni,t0) + Z2(t-t0) t0 < t < tj (8-27)
where
22 = [Z2(Xffl,t1) -Z2(Xfn,t0)]/(t1-to);
2a(t) = Z3(xra>t0)^ za(t-to) (8-28)
with z3 given by (8-16).
Now that we have formulated approximate expressions for each of the
parameters in (8-10) except w and wfi, we can define a constant wc which
satisfies
wc ) j~ dt = ] [3 + (Jc(qc-s)]wedt
(8-29)
to c 'to
A(xffl; to.
138
-------
where WG = w^; t0,ti) and
f'1
= Ms(ti) - M3(t0) + [w2(3 (8-30)
to
00,23 +
The parameter w is the effective cumulus cloud updraft speed w character! s-
c 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 >e we
assume
Wc0<>t) = Wr0$m; to.ti) t0 < t < ti, and
u ' " (8-31)
X->L, < 6
I ~ ~m I
where w satisfies (8-29). Note that (8-29) relates wr to the unknown
U C
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 cr
- "
c
and integrating it from to to ti along the trajectory ending at xm at time
we get
139
-------
tl tl
w C a.
-nGi.W-c. .
to to
(8-33)
+ 0) wedt
to
where 0 indicates that all parameters in the integrand are evaluated at
space-time points along the trajectory between x^ and x (see Figure 8-2).
Since xm is the site of a rawin station and ti is an observation time, we can
get Z2(xm,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(xi,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
= Az2(xm;t1,t0) ±
oz2(xin;t1,t0)
(8-34)
where Az2 and 6z2 are known, partly from (8-14). We can now write (8-33) in
the form
ti ti ti
f fas
AZ2 ± dz2 =6 -r^-dt - wr(t> TT§- dt +(h w dt. (8-35)
T A c J c T
t0 t0 ''to
The cloud cover fraction or(x,t) is known everywhere and we assume that w2 is
also available at all points and hours. Thus, (8-35) contains only 2
unknowns: w and w .
c e
140
-------
Figure 8-2. Illustration of the air parcel trajectory
that arrives at point x at hour ti. Point
xi denotes the parcel's location at time
Following a number of previous investigators (see, for example, the
*aview 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 =
z=z2
A6
(8-37)
where A6 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 v = d8/dz evaluated at z=Z2+e and h=z2-zt, where z^ is the local terrain
elevation. Eq. (8-38) suggests that within the vicinity of rawin station x .
we can approximate w by
WpQS.t) = G(xm;t0)t1)Q(x,t) , t0 < t < ti; (8-39)
x-x_ < 6
~ ~m
where G is an unknown function that we suspect from (8-38) is of the order
(8-40)
where T is some instant in the interval to to ti.
Substituting (8-39) into (8-35) we obtain
l j !
I * r o- r
Az2 ± 6z2 =014- dt " wr CpT~ dt + G(pQdt, (8-41)
7 c J c j
to 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
-------
AZ2 ± 6z2 =
dt -
[G J
[3 + ac(qc-3)]
(8-42)
ti
+ A] + G(^ Qdt
i
to
where
dt
ti
to
(8-43)
dt
We repeat that in our notation^ C^t denotes integration of £ evaluated at the
jr notation\ £dt
£tl * o
>(u Cdt denotes
J>tn
fixed point *_ whileQ Cdt denotes integration of C(x,t) along the space-time
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 - 6Z2 < AZ2 < AZ2 + 6z2 (8-44)
there corresponds a G given by
G = [AZ2 - ^) ^_ dt
tl
Qdt
[3 + ac(qc -
(8-45)
-i
1
vl
143
-------
Let us assume that the set of G associated through (8-45) with the set of
defined by (8-44) lies on the interval
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, ac (8-47)
Suppose that through a process like that described above we have managed
to extract from (8-45) an interval of G that contains some positive values.
Our~interest is only in the positive values, i.e., the G in the interval
144
-------
To each G in this interval there corresponds a w through relationship
(8-29), namely
(8-49)
"c = CWF§ dt]" {A + G [3 + CTc (qc
C
Q(xm,t)dt}
Let the interval defined by (8-48) and (8-49) be represented by
wc - 6wc < wc < wc+ 6wc . (8-50)
Like G, wc 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
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, za. (8-51)
At the conclusion of this operation we obtain a positive set of WG 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 v^Cx.t), w2(x,t),
and the cloud cover distribution cr , we can solve (8-2) for 22 at each grid
point and each hour t in the interval t0 < t < ti. By repeating the entire
process for the next observation interval ^ < t < t2, we can obtain 22 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 22 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 22, w2, Q and the
other governing fields.
Second, our interest in regional modeling is with historical situations
where available observations exist from which 22 can be inferred, at least
approximately, not only at the initial moment to of the simulation period but
also at discrete intervals throughout it. In general, the 22 predictions of a
model initialised 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 22 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
(8-52)
v v ^
j
From this equation one finds that the Fourier amplitude A of spatial variation 1
in the scalar c of wave number k = 2n/\ decays in time at the rate ;
|£ = -Kk2A. (8-53)
i
Sources of c can generate small scale variations but turbulence always :;
destroys them. Tj
_<
These observations suggest that if second and higher order derivatives of ]
a conservative, scalar quantity remain large for an extended period of time at j
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 |
i
(8-53), must be very small. Otherwise, the small-scale variations that cause .1
V
large values of these derivatives would be eradicated. •
Thus, our basic premise is that in cloud free conditions we can estimate
the ejevation at which the vertical diffusivity becomes vanishingly small by
examining the vertical profiles of the second, and psrhaps 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 0 in dry, convective conditions illustrated in Figure
8-3a. In Figure 8-3b we show the corresponding profiles of d2q/dz2 and
d26/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 Zi, z2, ... Zj that are not necessarily equally spaced. Let us denote
the £ values at these points by
C1 =£(xm, z.,t0) , i=l,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 £ as well as other properties of
interest at the desired elevation z" by expanding £ in a polynomial about z".
Thus, let
£(z) = a0 + axn + a2n2 + a3n,3 (8-56)
where
n = z-z". (8-57)
149
-------
\q
\
\
\
(a)
Figure 8-3. (a) Idealized profiles of mixing ratio
q and potential temperature 9 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
-------
T-epresentation is obtained by using the four successive measurements of £,
;hown 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 zlf z2, z3 and z4, as indicated in the Figure;
ind let the associated £ measurements be designated £x ... £4. In this
notation we have from (8-56) ijlj'
if"
ill
.(!••.«
Figure 8-4. Illustration of the 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".
£i = a0 + ajHi + a2qf + a3nf
• •
(8-58)
• *
£4 = a0
where
0,. = z.-z". (8-59)
Solving the system of equations (8-58) for the a's is straightforward and we
obtai n
151
I
H
t
i-('
iip
«!• (•
-------
a _ i
2 "
(8-60)
_
31 '
(nf-ni)
- (n2-ns)(ni-ni)
A22 =
= (Ci-CaXni-ni) - (C2-Ca)(nf-ni)
B2 = (Ca-C»)(ni-nl) -
(8-61)
(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 £:
= a0
32
(8-70)
(8-71)
z=z'
152
-------
= 2a2
z=z
(8-72)
923
= 6a<
z=z
z"+6z
z"-6z
(8-73)
= 2(6z)a0 + f(6z)'
(8-74)
Now that we have formulated a method of estimating the derivatives of
^arameters measured at discrete points, we can proceed to estimate T.-L.
Let
v
ft
ft
7
zz
(8-76)
The steps for Stage ZQ 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 qm(z,to) and 8m(z,to), respectively,
compute
qZ2(xm,zk,t0)
ezz(xm,zk,t0)
k = 1, ... 60
at rawin station m (=1 on the first pass through this stage) at hour to, where
153
-------
ZL. = z+(O + 25m + kAz (8-77)
x t ~m
and Az = 50 m.
• (2) Next form the products
Pk = ^22(l^'Qzz(z^ ' k = lj ••' 60 (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).
,10) = F(q(xra,t0), e(xffl,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 22 °f the mixed layer top
is likely to lie.
154
-------
We anticipate that a measure of 6z2 is the width of the interval
centered at z=zk* in which P., 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 6z2. At this time, no
information is available on which to base a quantitative rule for
estimating the range of z2 values. Therefore, we will assume for
now that
6z2 = 100m (interim assumption) (8-82)
and after we have gained experience with actual soundings, we will
attempt to formulate an empirical rule for estimating 6z2.
(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 z2. The elevation of the LCL is found
first in pressure coordinates as follows.
Let qn and 6n 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
Mn—- (8-83a)
(cf Eq. 3-2a), and its temperature will be
155
-------
(8-83b)
(cf £q. 3-6). The elevation p,c, 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 = eoexp[j^ - ^)] (8-83c)
o
where eQ = 40mb, TQ = 302°K, L = 2500 joule g"1, and R=0.461 joule
g °K . Hence, p is the solution of the equation
P°,n
/ i -i
^3U2 - T" (8'83d)
0.622 + qn
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 which the left side of (8-83d) first
exceeds the right side.
Using the pressure-height functions p (z,t ) available from
processor PI, convert p,£L into elevation z.£, (m MSL) as follows
zLCL(xn,t0) = z* (8-84)
where z* is the elevation for which
p(xn,z*,t0) = PLCL(~n)to^ (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
1 v xjrWz.to)
P(xn,z,t0) =! - ; - - (8-86)
"
At the end of this step we have values of ZLQL(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 ZLCL given by (8-84). As we noted
earlier, the decision rests solely on whether ovQSm.to) ls greater
than or equal to zero. In particular
, to) if CT(X,t0) = 0;
(8-87)
zLCL(xm,t0) (Eq. 8-84), otherwise.
where
zK.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 z2 and/or z,/»i . The first is
ffc(~m'to) = ° 2nd z|_CL(Vto)
-------
and the second is
ac(Vto) * ° ^ zLCL^m)to) >[z2^m>to) + 6z2J . (8-90)
The first condition indicates that clouds are not forming even
though our estimates of z,p, 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 z3. As we noted earlier (see Eq. 8-15)
pz2(xm,t0) + 100 , if ac(xm,t0) = 0;
Z3(xm,t0) = J (8-91)
-zTCU(xm,t0) , otherwise
where Zj^,, 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(xffl,t0) = q(xfn,z2,t0) = a0jq(z2,t0) (8-92)
where aA is the coefficient a,, of the expansion of q about the
o,q o
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)
qJXm.to) = qfem'Zs.to) = ao,q(l3'to) (8"93)
where a is obtained as in step (8-8) except the expansion is
°»M
about the elevation z"=Z3.
(10) From (8-8) and (8-74) we have
M3(xm,to) = Z*[Az ao>q(zk;t0) + ^ (Az)3a2)q(zk,to)] (8-94)
where k2 is the altitude interval given by (8-77) that is nearest
Z2(xfn,t0), i.e.,
|22(xnj)t0) - (zt(xm) + 25 + k2Az)| = minimum (8-95)
and similarly k3 is the integer that minimizes
UsC^.to) - Ut(xm) + 25 + k3Az)| = minimum; (8-96)
and Az= 50m as in (8-77). The coefficients aQ and a£ in (8-94)
are derived from the q sounding data using formulas (8-60) - (8-69).
(11) From (8-9)
(8-97)
(12) Repeat steps 1-8 at rawin station m for the next observation hour tj.
to obtain z2(xm,t1), Z3(xm,t1), q^.ti), q^C^.tj),
(xm,ti) and 3. (In general, ti = t0 + Atm, where Atm is
159
-------
the interval between observations at station m; so we should
actually write t. rather than ti 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
lA.t) = qc(xm,to) + qc(t-t0) 1
-to) >t0 < t < tx (8-98)
3(t-t0) J
where
^c = ^c^'1^ ' Ictem'^Kti-tor1 (8-99)
with similar expressions for 6^ 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 ' '
t02 t12
Station 1 ' '
t t t
om im . 2111
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 t0 < t < tx 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 .t^
during the interval t0 < t < tj. of the air parcel that arrives at station x^
at hour ti. By definition
x(ti;xm,ti) = xm (8-100)
** ^ * * /N'fiI /wifl
and by previous declaration (see Figure 8-2)
(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-z+ at XM at time ti is affected by the
* t ~m * J
surface heat flux in a plume-shaped area whose width a increases from zero at
x at t=tt to a value of the order
~m *•
a (x^) ~ (t!-t0)Av (8-102)
at the upstream starting point of the trajectory that ends at x^ at ti. 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 x! and x. we can assume that the difference in the
^MH *iM
integrated heat flux along the trajectory extracted from the winds v2n at the
level z2 and that along the "correct" trajectory is absorbed in the function
G.
162
-------
To compute the trajectory xi 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 v,, at any given
space-time point (x,z,t) in the NEROS domain. Thus, when we write w(x,z,t) or
vH(x,z,t) it will be understood that these values are available from the
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:
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/tjx^.ti), 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 < tx and
express the results in the following functional forms:
<-!
4 F§C
"to ^ J=°
t-
f
to
where
CTcm(t) = ^(x^jx..^).!) (8-106)
QBtt> = Q(x(t;xfl|ftl),t) (8-107)
where
,t1)]. (8-108)
(3) Now compute the following path integrals:
, - -c " jt0 l-^Cti-jAt) (8-109)
to
Qdt = XQin(tl) = AtCtrJAt) . (8-111)
J = (t1-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 TW, Ia and IQ (Eqs. 8-109 - 8-111) for
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;t!,to), which will
provide the entrainment velocity field w (x,t) through (8-39); and equation
(8-49) for wc(xm;ti,t0) which will provide the cumulus updraft speed wc(x,t).
(1) Using the estimates (8-98) of the temporal variations in qc, and
3 at each upper air site x^ in the period t0 < t < tx, and the
surface heat flux QCx^t) in this period, compute the following
integral:
J
= AtE[3 + o-^.t,
(2) Compute the parameter (see 8-43):
qc(Xn,.'ti-JAt)a (x tj-jAt)
-
where j(t!) 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 ti.
""ID "Til . "Til
As we .noted earlier, we have an estimate of z2 at xm at ti but no measurement
* ~m i
of z2 at x' at t0.
*wm
III
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 ±6z2(x|n,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
* ~m
will not represent the lifting condensation level, which is the altitude of z2
wherever ac(x,t0) $ 0. Thus, we assign Z2(xm,t0) values according to the
following rule:
' 1f ac(*m'to) = 0; (8-115a)
, otherwise (8-115b)
166
-------
vhere
ZLCL is glven by (8~84)> Z2(x-j,t0) is given by (8-87) and the summations in
(8-115a) are over all M rawin stations while in (8-115b) they are over all N
surface weather stations.
We apply the same |r| interpolation formula to estimate 6z2 at x' At
all times t = tls 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(x||J,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 ti < t <_ t2, we assume 6z2(t!)=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 t0 is the initial instant of the regional model simulation
period, determine Z2(xm,t0),m=l,. . .M from (8-87) (Stage ZQ) and
z2(xi,to) from (8-115). Then estimate 6z2(xm,t0),m=l,.. .M (see
step 4, Stage ZQ) and subsequently 6z2(xjJ1,t0) using these
values and (8-115).
b. If t0 is not the initial instant of the model simulation, then
is available at all grid points x from Stage Z2,
described below; and 6z2(x|j1,to) = 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
and (x_,ti). We find after some thought that
6z2(x;,to>] (8-117)
(8-118)
(5) We now solve (8-45) for the minimum and maximum values of G
(see 8-46) in the vicinity of x during the interval t0 < t <
^ Ifl.(ti) or q^x^) = 0 and IQm ^ 0;
.t^rtn - I^CtO* ^ClQ.Cti) - (8-119)
c
I(xm;t1,to)/qc(xin;t1,t0)3"1, if
W*^ aafi W an- iQm ^ o;
1f Vtl} - °'
W1'th "*<*M rePlaced ^ ^(x.t,)- (8-120)
In Eq. (8-119), AZ2( )m1n is from (8-117), AZ2( )max is from
(8-118), I^tti) is from (8-109), T^Cti) is from (8-111),
Jcm(t1) is from (8-110), I( ) is from (8-113), qc( ) is from
(8-114), and A is given by (8-30).
168
-------
(6) Next we must check whether the upper bound G__v on G is
lUoX
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, qc, etc. in (8-120)
yield a value for G_,v that is positive, and preferably of an
ilicLX
order of magnitude consistent with (8-40). Tne 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_,t1)_,v > 0, go to step 7; otherwise return to
/VHI IHQ.X
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 < G(xm,t1) < G(xm,tl)max (8-121)
where
Go = -axCO.GCx)]. (8-122)
We now compute limiting values for the cloud updraft velocity
parameter w .
°ml • 1f ^
1f '^
169
-------
where
- M3(.xm,t0)
ti-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,
qc, 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 .m.
am
"ctem'tiW = EP- (8-123) with GO replaced by
c ~m max (8-125)
(8) If wc^'tl^max - ° and Iom^tl^ — qc^m;tl
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 rC^^i^ax- 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 w ( ) . < w < w ( ) has been found
C ml n — C "" C mOlX
that lies at least partially on the positive real axis, proceed
to step 9.
170
-------
(9) At this point a range of values of G(x m,ti) has been
established (in step 7, Eqs; 8-121 and 8-122); and a range of
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
ti:
G(x_,ti) = value in the range (8-121) that
~m satisfies (8-40) closest; (8-127)
w (x.ti) = solution of (8-49) with G
c ~n . given by (8-127) and all (8-128)
other parameters with values
determined in step 8.
We assume that these values are constant at x during the
~in
entire interval t0 < t < tlf specifically, we assume
G(x,t) = G(x,t1) 1
~" "* to < t < tl (8-129)
wc^m,t) = w^x^.tj J
These values should be recorded for each hour in the interval
t0 to ti for site x. Note that tlt which is the hour of the i
~m j
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 WG.
171 ] |
i :
i 'I-
-------
(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 to to to + T.
(11) Repeat steps 1-10 for each rawin station m=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:
M
I
-•
G(x,t0) = - (8-130)
I |x-xj G(x.to)
' -- ffl ~flr u
where x ranges over all grid points.
(2) Convert the G field into the wg field as follows (see 8-39):
wp(x,to) = 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 wfi(x,t) (designated f6/A9 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 12 field at each grid point using the
interpolation scheme (8-115) employed earlier in Stage PATH.
Specifically,
f Eq. 8-115a, if or(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 WG 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 WG and cloud depth
using all wc, 22 and Z-J-QJ 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 WG values obtained in Stage WEWC
for various values of the cloud thickness h3 = ZJQJ - z2, say values
of hs at intervals of 200m; and call the resulting function W(h3).
Then we assume
= (l-a(x)MzTCU(x,t0)-z2(x,to))
(8-133)
173
-------
where
W(hs) = w for clouds of depth h3 (an empirical (8-134)
function);
a(x) = distance weighting function
(8-135)
= exp[-|x-xJ/SIG]
x, is the rawin site closest to x where w,. ? 0; SIG is a distance
~ni "* c
constant, e. g. , 50km; and
"cLOCAL(*>to) = "A'^' <8'136)
Formula (8-133) is an heuristic expression that assigns w a value
that is determined partly by any wc 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 v^ required in Eq. (8-2) at time
t0, we will use
= [Z2(x,t0)-zt(x)]"1 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 v-Xx.z'.t) at elevations z. < z'< z.
~n ~ t •— —
(6) The vertical air speed at each grid point x at the initial moment t0
is obtained from the function routine WV as before, namely
to) (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
AtS(x;to)]dx'
where
S(x',to) = [W2(xjto)-o (xjto)w (xjto)]-
L c (8-140)
(l-ac(x',t0)) + we(x',t0)
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); WG is
from (8-133); WQ is from (8-131); and
€
p(x,t|x',t') = known function of v2H- (see Chapter 9, Part 1)
(8-141)
175
-------
* .*
(8) At this point we have computed z2(x,t0+At). From this field we*
construct z3 in the manner of (8-91), namely
Z3(x,t0+At) =
Z2(x,t0+At) + 100, if ac(x,t0*At) = 0;
zTn,(x,to+At), otherwise
(8-142)
-)
where Zypy is the known elevation of cumulus tops. This calculation j
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 f'
each grid point at time t0 + At as follows:
H3(x,t0+At) = Z3(x,t0+At)
• -W(x,z3(x,t0+At),t0+At)
(8-143) 3
where
Z3(x,t0+At) =
At)-z3(x,t0)
At
and
= v(x,z3(x,t0+At),t0-HAt)
calculated from Vs.. and z3 in the
manner of (8-A2). "See the Appendix
to this section.
(8-144)
(8-145) f
(8-146)
.-ft
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) = H3(x,t0+At) (8-147)
(10) We must also output the local time derivative z2 at each grid point
and each hour:
- Z2(x,t0)
Z2(x,t0+At) = . (8-148)
assume
Z2(x,t0) - 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) - h^x.to+At) (8-151)
\
where z2 is from (8-139) and ZT and hj 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,to+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 Y(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 .¥' of this parameter based strictly on heuristic
notions. In Processor P12 we test whether the V values derived
here satisfy criterion (5-47) of Part 1. If they meet this
condition they become the final estimate of the parameter V.
Otherwise ¥(x,t) is given the largest value that satisfies (5-47).
As a rough, heuristic approximation we shall assume
f
[-
2h3(x,t0+At)
r (x,t0+At) = 0.5 exp -[22(x,to+At)-ZT(x) 3 + "
Under this assumption, V •* 0.1 as the depth of cumulus clouds
becomes large compared to the depth (z2-z-j.) of the subcloud layer.
For shallow cumulus, V •» 0.6. The V 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)
<«3 = f(YH-v)dp (8-155)
179
-------
<6(x,y,t)>2 =
P2(x,y,t)
, if
Pvs(x,y,t)
P2(x,y,t)
> if
(8-156)
PjCx.y.t)
<5(x,y,t)>1 - p .* j (VH-v)dp , defined only when Apx=0 (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 Pj^o. P^^i or PVS* and pl0 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 V.,*v. We leave the
method of computing the divergences (£H*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 pn that bound
each layer by (see Eqs. 11-Ib and 11-40)
180
-------
= ^- [Qn(x,y,t) - tt^Cx.y.t)] (8-158)
^
where
Apn = iPn0^'^ " Pn-i^^'^l (8"159)
and uin 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
u>n(x,y,t) - wn(x,y,t) (8-160)
where u>n denotes the local value of tu on pressure surface p at the cell
center (x,y).
By definition
» 5 If + W
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
Two of the fields that we need in this processor, P8, are w. and w*, 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 u3 in the form
+ u)Q (8-163)
where WQ is the value of tu at ground level.
To estimate u>0 we note first that the vertical velocity w at ground
level is
wrt & v •VMZ.,. (8-164)
0 ~o ~H T v '
where z-p 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*PQ- We also expect the last term to exceed the magnitude of
the local time rate of change of pQ. Thus, to good approximation we can
assume
(8"165)
182
-------
Combining (8-162), (8-163), and (8-165) we have
wi * TTT; C~af "*" Xq'SwPq ~ Ap.<6>, - AF
o p^cj ot ~o ~n J o o
••• P.gv_'VuzT], mode 0. (8-166)
where
AP-L = jp^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 u> estimates.
To estimate the densities p, and p , we use
Pn(x,t)
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^Zj in Eq. (8-166), it will be necessary to
fit a polynominal of at least 3-rd order to the pressure anditerrain data (see
Appendix to this chapter).
The expression for Wg is obtained in the same manner that (8-166) was
derived. We get i
183
-------
W2 =
(8-168)
'* pogWT]' mode °
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 I (see Eq. 3.1b of Part 1). We denote this
velocity component by WQ, where
wl = WT1 + "Dl = pTg caT
(8-169)
Since by definition WQ.^ is the divergence induced vertical motion on p,,
the terrain induced component is just the right side of (8-169) with <6>,=0.
Thus ,
^PT
wAi = ' 7TZ <6>i « mode °- (8-170)
and Ap, is given by (8-166a). In order to emphasize that this value applies
in mode 0 only, we designate it W.U at the output of Pll. Processor P12 will
take this field as input and generate the final, complete function WD,
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
vith terrain features that extend up through the inversion layer but which
otherwise is the top of the inversion layer itself. In this situation the
3xpression for u>3 becomes
u)3 = Ap3<6>3 + Ap2<6>2 + u)vg. (8-171)
In order to obtain an estimate of w 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
3
where (%»v0twQ) is tne fluid velocity in the cold layer; and nvs 1S tne
cooling rate, a known function of space and time. Just above zvs> i.e., at
the base of Layer 2, the downward volume flux is
dt
Continuity of the flux across Hys requires that (8-172) and (8-173) be equal,
which is true if
185
-------
Ji *y 5^7 5^7
aT* + "vs-aT * vvs^ ' -W
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
"Vs ' 'Pvs^vs (8'175)
where wys is given by (8-174).
Combining (8-162), (8-171), (8-174), and (8-175) we get
W = C *
3 = pTg C§r * W3 ' Ap3<6>3 *
P7 C-5T ^vs'Vvs - nv$L n°del (8-176)
and similarly
«, = -i- c^
2 p-g L at
(8-177)
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 s-
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
q (z,t)
Description
mixing ratio at
Source
PI
Stage
ZQ
Output
Variable
z2(x ,t)
Description
elevation (m MSL)
Qn(t)
en(t)
oc(x,t)
elevation z (m MSL)
at hour t over rawin
station m.
potential temperature
(°K) at elevation z
(m MSL) at hour t over
rawin station m.
PI
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.
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. ~
of top surface of
model layer 2 over
rawin site x at time
t of rawin soundings.
expected range of
values, centered at z2
at which actual mixed
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 mixing
ratio entering cumulus
clouds over rawin
station m at time t
(available at 30 mi nut
intervals).
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
Output
Source Stage Variable
Description
zTCU(x,t)
elevation (m MSL)
of highest cumulus
cloud tops in grid
cell centered at x
at time t. ~
RAW
ZQ, 3 average mixing ratio
(Cont.) ~ in Layer 3 over rawin
site m at time t
(available at 30 min.
intervals).
z, p,(x t) elevation (m MSL) of
LLL ~n' the lifting condensa-
tion level over sur-
face station n at
time t (required
only at the hours
t of rawin observa-
tions)
vH(x,z,t)
w(x,z,t)
Q(x,t)
a (x,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 j. P4
heat flux (m°K sec" )
in cell at x at time
*w
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-111) units =
m °K.
189
-------
Table 8-1. (continued)
Input Output
Variable Description Source Stage Variable Description
)>3 average mixing ratio Stage WEWC G(x ,t) entrainment veloc
in Layer 3 over ZQ scale factor (uni
rawin station m at °K~1) (see 8-39)
time t (available at rawin site x at
30 min. intervals). hourly intervals
q (x ,t) mixing ratio Stage wc^m'^ 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).
Q(x,t) surface kinematic P4
heat flux (m °K
sec- ) in grid
cell at x at time t
(at hourTy intervals).
o\.(x,t) fractional coverage RAW
~ of cumulus clouds
cumulus cover path Stage
integral PATH
) elevation of top Stage
surface of Layer 2 ZQ
at rawin site x^
at times t of
rawind soundings.
,t) lifting condensation Stage
level at surface ZQ
weather station n
at times t of rawin
soundings.
w path integral Stage
leading to rawind PATH
site >< at time t
of rawrn sounding
^nn/t) surface heat flux Stage
^m path integral PATH
leading to site
>c at hour t of
rawin sounding.
190
-------
Table 8-1. (continued)
Input
Variable
M'V
Description
vertically inte-
grated water mixing
ratio in Layer 3
over rawin site x
(at sounding hours t
only).
Source
Stage
ZQ
Output
Stage Variable
WEWC
(Cont. )
Description
vertical air speed Stage
(m/sec) at (x,z,t) WV
water mixing ratio Stage
above Layer 3 at ZQ
rawin site x_ (at
t = 30 nrin 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
KSur t.
TCU
,t) hourly values of Stage
the entrainment WEWC
velocity scale factor
at rawin site x .
~m
,t) hourly values of Stage
the cumulus updraft WEWC
speed (m/sec) at
site x .
(x,t) 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 t1!1
W2
w (x,t) = entrainment velocity
fe/A8(x,t) at hour t at grid cell
centered at x (m/sec).
wr(x,t) cumulus updraft speed
" ~ (m/sec) at hour tin
grid cell centered at
x.
wz(x,t) 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 H£ (top of
Layer 2) in grid cell
centered at x at hour
t.
191
-------
Table 8-1. (continued)
Input
Variable
zt(x)
Description
mean terrain
elevation (m MSL)
in grid cell
centered at x.
/w
Source Stage
P7 W2
(Cont. )
Output
Variable
z3(x,t)
Description
elevation (m MSL)
surface HS (top o
Layer 3) in grid <
centered at x at t
t.
w(X,Z,t)
Pm(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
elevation z (m MSL)
at hour t over rawin
station m.
PI
depth (m) of Layer 1 P7
in grid cell centered
at x at time t.
Ha(x,t) volume flux (m/sec
through top surfac
Layer 3 in grid ce
at x at hour t.
*%*
Z2(x,t) local time derivat
of elevation 7.2
(m/sec) at grid
cell centered at x
at hour t.
ns(x,t) thickness (m) of
Layer 3 at hour t
in grid cell center
at x.
thickness (m) of La
2 at hour t in gric
cell centered at x.
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 center
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
v«>
Description
terrain elevation
Source
P7
Stage
WV
Output
Variable
vH(x,z,t)
Description
horizontal wind
O.Cz.t)
On(t)
Ap-^t)
centered at x.
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.
w(x,z,t)
Wrj-,(x,t)
<5(x,t)>3
~
2
<6(x,t)>1
[VH=(U,V)] at (arbitrary
eTevation z (m MSL),
time t at site x.
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 surface
of Layer 1 (Defined
for daytime hours only. )
average horizontal wind
divergence (sec-1) in
Layer 3 in grid cell
centered at (x) at
hour t. ~
same as <6>, except
applies to Layer 2.
Same as <6>3 except
applies to Layer 1
(values of this
quantity are computed
only for daytime hours.)
same as <5> except
applies to Layer 2.
same as <5> except
applies to Layer 1
(values of this
quantity are computed
only for daytime hours. )
p,(x,t) same as p, except Stage
L ~ elevation^of top W2
of Layer 2.
193
-------
Table 8-1 (Concluded)
Input
Variable
Description
Source Stage
Output
Variable
Description
"vsv~'
same as p- except
top of Layer 1.
same as p3 except
top of raaiation
inversion
elevation (m MSL)
of virtual surface
in cell at (x) at
hour t. Note: z
ZT when
= 0.
P7
P7
P7
vs
ntf,.G<,t)
'vs
growth rate (m/sec)
of the radiation
inversion layer
depth.
P7
194
-------
Appendix to Section 8
Some of the equations in this Processor, such as (8-166), contain terms
of the form
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
(rYHP>i,j = ui,A(pi,j> * VI,A(PI.J} (8"A2)
where
VpI,J) = I (PI+I,J " PI-1,J} " 12 (PI+2,J ' PI-2,J}
and
Ay (PI,J) = I (PI,J+1 ' PI,J-1} ' IZ (pI,J+2 ' Pl.J-2^
Here PT , is the value of p at point (I,J), i.e., column I, row 0, with rows
l,vJ
parallel to the x axis and columns parallel to y.
195
-------
YYYYYYYYYYY YYY
N
)
i
9
1
3
^ 1
lyl
"rV
* A" ,
«• S •
3 *! .
T S '
K
j*
M
I
y
5
V
^ i<
i .5
9 15
3
Q.
O
•o
CD
«•«
3
Q.
C
V)
•**
•o
C
CD
00
I
o >
0) C
O w
o o
«• w
Q. co
(D
o a
•=•§
o w
t£ O
(D u
E 42
CD w
CO
ob
£
O)
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
af (YP) + V * af ^ ' I + 5 ' ? • <9-4>
Thus, the equation governing the mixing ratio Y has the same form as the
equation governing the concentration c (i.e., Eq. 2-1 of Part 1) except that
the inhomogeneous terms S, R and W are divided by the local air density.
Notice that Eq. 9-4 indicates that material is well-mixed (i.e., 9Y/3t = 0)
vertically when the mixing ratio is constant in z. However, the concentration
form of the continuity equation yields the contradictory result that the
well-mixed state is achieved only when the concentration c is independent of
height. The latter is in error because it fails to take into account the
effects of gravity on the mass distribution of material in the atmosphere.
The mixing ratio form of the equation handles this effect implicitly through
the link with the air density p. In short, Eq. (9-4) is the proper form of
the continuity equation for use in applications to atmospheric layers deep
enough (> 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
reacts, k is the rate constant (units of mole m3sec ) of this second-order
reaction, and k-j 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 use 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 ejevation (m MSL) of that station. Now define
m om
p|"(2)dz (9"8)
zm
"on, + hlm
hom
.
-------
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 wit.h the intervals at which T (z) and Pm(z) are available (from PI).
Similarly, we define
(m om 1m
Ttfn,U)dz (9-12)
hom
(9-13)
3n, '
• + hom +hlm +h2m
The temperature and density profiles T and pm are available from PI
each hour but they give values only at the measurement station sites >r.
Therefore, the layer averaged density and temperature values derived from
(9-8) - (9-14) must be interpolated to each grid cell location. Use the r
weighting scheme as follows:
^ >1>n»
= —S^j n = 0,1.2,3 (9-15)
nil II ^
Z rm
m=lm
=-=± n = 1,2,3 (9-16)
M .1
rn^l^
201
-------
where
rm = K^-V2 + (Jty-yB)2]. (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 n 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
k* = 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 n 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
•here ka(s) is the "clear sky" photolytic rate constant for reaction a and E
cc) is an empirical function of the cloud cover cc = cc(h), which is the
-raction (0 < cc £ 1) of the sky covered by clouds of height h. We assume
further that the clear sky constants kffl are contained in the chemical kinetics
subroutine of the regional model code in the functional forms k (s 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 (1,0) at hour t (9-21)
(dimensionless).
The zenith angle can be obtained from standard astronomical formulas given
the latitude (JA) and longitude (IA\) 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 tm. We will not elaborate on the computation of and E here.
m &
203
-------
In summary, the solar zenith angle <]> is used in each grid cell and hour
to determine the clear sky rate constant kff for each photolytic reaction a.
These values are multiplied in turn by the cloud cover factor E for that cell
and hour to arrive at the corrected rate constant
The photolytic constants are not modified by the density correction terms
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
vm
(z,t.)
virtual temperature
(°c) at elevation z
(m MSL) at hour tk
over rawin station
PI
DEN
0
m
same as T except
density (KSJnr3)
elevation (m MSL)
of rawin station m
h0(I>J»t.J 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 m"3)
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, cell
(I,J), hour tk (for
temperature correction
of rate constants in
model)
h-,(I
h2(I
h,(I
,J,t.)
N
,J.tk)
,J,tk)
fv
same as
Layer 1
same as
Layer 2
same as
Layer 3
hft except
0
h except
h except
P7
P8
P8
2 same as
Layer 2
o same as
* J Layer 3
, except
JL
1 except
cc(x ,h,t ) fractional sky
coverage of clouds
of height h (low,
middle, and high)
at surface weather
station n at hour
RAW ZEN »sa,J,y
SCI.J.V
solar zenith angle
(degrees) at grid
cell (I,J) at hour
t (used to obtain
pnotolytic rate
constant values).
cloud cover cor-
rection factor
(dimensionless)
for photolytic rate
constants in cell
(I,J) at hour t
205
-------
D
O) O.
Q. H
YY.YYYYYY
A
OJ
A
cj.
V
n
A
«x
V
A
V
CM
A
V
en
A
V
Q. ^
S- O
O 2
(O -M
10 0)
+J
V> JC
3 4->
I/I
U 9)
•p- o
-»-> 10
§<4-
QJ
f -M
U C
«/) -i-
I
cn
01
206
-------
SECTION 9
PROCESSOR P10
DEVELOPMENT
Processor P10 transforms the emissions inventory into the source strength
«v
functions S, S^ and S«. 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 Ah. (t ) of the k-th source at hour t as follows.
l\ HI Hi
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),tra] < 0;
Ft 1/3 (10-1)
3 (^|) , otherwise.
Here L[I,J,t ] is the Obukhov length (meters) in cell (I,J) at hour t , and
VTak 2
Fk = g[ K alC 3 («Dkwk. (10-2)
* Tk H- 273 * K
In this expression TV is the exhaust temperature (°C) of the stack, D^ is the
-1 -2
stack diameter (m), wk its exhaust velocity (m sec ), g = 9.8 m sec is
gravity, and
N _2
VnkTvn
-------
is the estimated surface air temperature (°C) at the location of source k. In
(10-3)
(*n»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)
u = [(1)2 +
(^(KlO.JdO.t^)2]*
where 1 and 1 are the Layer 1 averaged wind components (m sec"1).
Finally, the parameter s in (10-1) is approximated by
.2
s =
2 ~2 ~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
m
where z . 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 and 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 OELH.
210
-------
Stage S
Having estimated the effective discharge heights of the major point
*v
sources in Stage DELH, we can now compute the source strength functions S, S^
and S2.
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 ow 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,tm) 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
m
pollutant have the following forms.
LAYER 2
H1
_^^f^fm^:^-:-'-
^...............-rrr^'^v^vvK--- -
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 (1 u
-------
K
k=l .
1 = C(aTO-an)E + 2 E,
k=l
k=l
where CTTO, E, HI§ and h2 are functions of (I,J,tm); and ER and I
depends on t . The latter is defined by
.0, otherwise.
< «k < (Zj -
where
^ = 2T(I(k),J(k)).
z =
zT(I(k),J(k))
h1(I(k),J(k),tm)
h2(I(k),J(k),tm)
and h is the thickness of Layer n. Also
(10-
212
-------
A2 = A (10-11)
Al ~ (1 " aTl)A (10-12)
A0=(l-aTO)A • ' (10-13)
A = a2(cos
S(I,J,t) T
(ppm m/sec)
'V W.J.V'o
where is the average air density (units = 10 moles m ) in
Layer n, time tm 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 result*
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
Nmj(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
Cd.J.t.) = C(Nmp +Nmj)(h0) + Lh0]/A .(10-18)
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
LCI.J.V
Description
Obukhov length
(m) in cell (I,J)
at hour t
Source
P4
Stage
DELH
Output
Variable
• w
Description
effective source height
(m) of k-th major point
source at hour t .
.
m
virtual temperature
(°C) at surface
weather station n
at time t.
.
, Layer 1 averaged
east-west wind
(m/sec) in cell
(I,J), hour tffl.
ra
same as j except
" north-south
component
temperature (°C)
of exhaust gas
of major point
source k at hour
V
diameter (m) of
stack k.
P3
Pll
Pll
.RAW
RAW
wt^tm^ exhaust velocity
m (m/sec) of stack RAW
k.
[I(k),J(k)] grid cell coordinates RAW
(row J, column I) of
major point source k.
z k stack height (m) of RAW
source k
E,/t ;°0 emission rate (moles/ RAW
K m hr) of species a at
time t from major
point source k
S(I,J,tm;cO
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
Description
Output
Source Stage . Variable
Description
E(I»Jft_;ot) emission rate
m (moles/hr) of
species a from all
except major point
sources in grid cell
(I,J,) at hour t,
ovn(I,J,t ) fraction (non-
m
JTO
m
dimensional) of
model surface H
penetrated by
terrain.
crT1(I,J,t ) same as ovn except
IJ- m applies tiusurface
zt(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)
Ml.J.tJ thickness (m) of
Layer 1 at time
tm, cell (I.J).
m
h2(I,J,t ) same as h, except
applies to Layer 2.
. (t )
effective source
height of major
point source k
2
mean air density
(moles m~3) in
Layer 2 at hour
tm, cell (I,J)
same as 2
except applies to
Layer 1.
RAW
S
(Cont.)
P7
P7
P7
P7
P7
P8
DELH
P9
P9
S2(I,J,tm;a)
emission rate (p
sec'1) of specie
at hour t from
sources in Layer
cell (I,J).
emission rate (pi
seer1) of specie
at hour t from .
sources in Layer
cell (I,J).
total number of i
point sources in
0, cell (I,J) at
V
216
-------
Table 10-1. (Completed)
Input
Variable
Description
Source Stage
Output
Variable
Description
0
same as 2
except applies to
Layer 0.
P9
L(I,J) total length (m)
of major line
sources in cell
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 tm
h (I,J,t ) depth (m) of Layer
0 m 0 in cell (I,J)
at time t,.
m
RAW
ZTA
.fraction (0<£<1) of
Layer 0 in cell
(I,J.) filled by line
and point source plumes
at hour t
.
m
RAW
Stage S
P7
217
-------
YYYY
•i
03
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 b
continuous functions fit to discrete meteorological data, or by flow field
derived from mesoscale meteorological models. In this study we use neithe
of these approaches. Following the consideration presented in Chapters 6 ani
7 of Part 1 of this report and in Lamb (1983a), we use wind observations an<
physical principles jointly to define a set or family of flow fields eact
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 pff and pR (pa > pft, i.e., surface a has the lower elevation)
Pp(x',y')
<())(x>y)> = _£. \ l |(x;yl,p)dpdx1dyl (11-la)
P Vap
A Pa(xly')
where A denotes the rectangular domain x - < x'< x + *, y - ^ < y1 < y +
of a model grid cell; and $ is an arbitrary scalar or vector function. The
averaging volume Vtfg is given by
Pp(xiy')
Jy1 ~ (11-lb)
A pa(xly')
Note that since p^ > p0, both V _ and the integrals in (11-la) are negative
U p ™"""""" up --
(assuming i)»0).
Suppose that observations of v are available at N arbitrary sites on
surface pfl, i.e.,
222
-------
n=l,2,...N.
where (x ,y ) denotes the location of observation n. (Throughout this
section, we shall use the caret (~) 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
vm(p,t) = v(xra,ym,p,t), m=l,2...M. (11-3)
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
^(Xjy.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 g^ =
u *" —U
(gua'gva) where
n fy v n
gvcr(x,y,p,
223
-------
VERTICAL SOUNDINGS
SURFACE OBSERVATIONS1
AT N SITES
"9ure
A.o.^' th° H°"-*
Available on the Bottom SurS o? Layer * "" 9re a'S°
Assuming that j and v vary by only «„ fract1ons of themselves ^
ranges over the area A of any grid „„, we get fron (u.1) ^ ''
0 *
1 ga(x,y,p,
* Pg J "
t)dp
(11-6)
or
p * ^(x.y.
5»/y^'x»Jr » ''Jo
P
where
= v(x.y,paft),
(11-7)
(11-8)
224
-------
i
and 0 is defined as in (11-1).
-a P
] . At the M sites of the measurements of the vertical wind profile, we can
i • .
! evaluate D. We have
: ~ P
.y t)>
8 at each point (x,y) in space by
interpolation of the "observed" values 0. That is, we assume
sown p
p *
m
m=l m ~
_
where W (r)is a weighting function, e.g., r , to be selected later, and r is
HI *** *^* *"**
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-
proximate the layer average winds at all those points (x ,y ), n#n, 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
8 is given by (11-10), and the vector product on the right side of (11-11)
is defined by
E (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
z k-Vxv =
ax ay ap
U V W
(11-12)
and the horizontal divergence 6 is
A _ 3u . 3v .
6 = +
(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
« f** **
(11-12)
t
ay
(11-14)
226
-------
From Part I, page 21, we have
' v - ' + 1n
. —aT* —~*
A /.. a
where A is the horizontal area on which < > is defined; V _ is defined by
(11-lb);
and
CxiylPtxly1 ))dx'4y- (11"16)
A
with a similar definition of the surface average
Assumption: On pa, v(x,y,pa) - and on p
v(x,y,pft) - . (11-17)
~ p ~
In this case we have
- \ \ 9pcr dx'dy1,
" A J J ax1"
227
-------
Note that
_3_ I 1 I dpdx'dy1
' J * -» /
A Pa(
IT al 11 CPB(xly')-p (xjy1)] dx'dy' (11-19)
where the domain A is a function of x inasmuch as A = x-Ax/2
-------
Using Helmholz theorem that any vector can be written as the sum of a
rotational and a divergent component we can write
v = v... + v
~ ~¥ ~X
from which we have
(11-23)
*• . (11-24)
Since by definition vju is nondivergent, we can express it in terms of the
**T
stream function V:
= k x W (11-25)
and since v is irrotational we can express it in terms of the velocity
A
potential x:
<> = VX. (11-26)
From (11-25) we get
< V =
i i *
001
ajF &¥ 8V
3x 3y 3p
(11-27)
and from (11-26)
=
~ ax ~ ay &'
(11-28)
Combining~(ll-24), (11-27), and (11-28) we get
(ll-29a)
229
-------
<»» • g * fjj * ; . (11'29b)
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 V and x that are consistent both with physical
laws and with our observations rio,n=l. ..N (see Equation 11-11).
•~~— - *" n n UD
(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
f~ + (v-VH)v + ur~ + f k x v = - V4> + Kv2y + FH
(11-31)
•Vu-v + 9u)/3p = 0 (11-32)
"^n **•*
230
-------
where u> is the "vertical velocity" in pressure coordinates, i.e.,
io =
is the geopotential height, f the Coriolis parameter, K the "eddy
viscosity", F,. 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
!! * $ * 4i! .„.-» + KtSjj . , . FHX (n.34)
3t 3x 3y 3p 3y «, 2
9x oy
5 = - au)/ap (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
(11-37)
r3u> 3v 3iu 3un _ „ r32t . 32g-, . 9FHy - 9FHx
C5x 3p " 3y §p] " K [2 + 2] 3x "3y~
231
-------
The equation governing <6> is found by averaging (11-36) in the manner
of (11-1):
pp
|dpdx'dy'
P«
<6> = r II K^iy') ' «"B(x',y')]dx'dy' (11-39)
otp J •* p
where u»a denotes u(x;y;Pcf). Using (11-16) we can write (11-39) in the final
form
= [3°- ;#]. (11-40)
Later when we apply (11-40) we will describe approximations of 3", 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
3v 8u ~
§p' 9p ~ °' (under assumption 11-17) (11-41)
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
will also omit the "beta" term 3f/3y in (11-37). With these approximations we
can write (11-37) in the simplier form
232
-------
(11-42)
where we have used the approximation (see Haltiner 1971, page 151)
FH = - gI (11-43)
and T 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:
• K <«> + «> a! ln §1
(11-45)
. 8p 3p
V
+ ^ ln(Vap)
Vap 3y 8y
233
-------
are quite complex and we will approximate them simply by
9x2 > »? (11-48)
Combining all the above and assuming <£6> = <£><5>> we get
a
+ 'YH<£> - <£> Cg|lnVaB + .VHlnV * 2<6>]
;•>
(11-49)
where
« «'' (11-52)
and
In deriving (11-49) we aSsumed that j2 = 0 on pp ,„ anticipation of the
eventua! use of this equation In applications where the upper surface p Is
above the friction ,ayer. We win use Haltlner's approximation (page 152) for
TZ, namely
234
-------
tz = P0CD (11-54)
where pQ is air density in. layer a,p, in particular the mixed layer; CQ 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
2
at ^ + J <£> + J <£> + <£>G - K[-2« + i-» ] = H (11-55)
y ax ay
where
H E - f<6>
Va6
3
t • ^ ' - • (11-58)
We are now ready to collect the equations that we will use to generate
the V and x fields. First, let the superscript denote the time variable,
i.e.,
235
-------
= , i.e.,
(11-60)
00
J = Af(V(k)elJS'~
-------
Note from (ll-61a) that
JK/V3 =
ax <•
••
-00
y
(11-64)
£, J»
9 J ^'X
k5r(Js)e dk dk ,
-09
(11-65)
and in general
an
axr
n J
k) 4
(11-66)
Now, replacing each term in (11-60) by its Fourier transform and making use of
(11-63) and (11-66) we get
+ k2)|J(2n)2
(11-67)
00
J
-OB
237
-------
where
GJ <-» GJ (ll-68a)
HJ «-»• ffJ (ll-68b)
-» J (ll-69a)
-» J (ll-69b)
We want equations for V and x- Thus, we obtain from (ll-30a) and (11-66)
(11-70)
where
^(x) ^ AJ(k). (11-71)
Using (11-70) in (11-67) we obtain one of the equations governing 4*. namely
y k2 + -2
X
oo
f
At
9
(11-72)
238
-------
The transforms U and V of the velocity components and can now be
expressed in the form (see 11-29)
= -ik AJ + ikxBJ + GJ6(jk) (ll-73a)
= ik A° + 1k 8° + vJ6(k) (ll-73b)
x y ~
where
XJ(x) *-» SJ(k) (ll-73c)
and S(k) is the delta function of the wave number vector k. Then from
(ll-30b), (11-40), and (11-66) we find that BJ must satisfy
(11-74)
where
A a P
<5> = ^— — [ui (x.JAt) - w (x.JAt)] (11-75)
"
and u> is considered to be a known variable.
Finally, the observations nQ given by (11-11) must be satisfied.
** ~n UD
Hence, from (11-29), (ll-61a), (11-66), and (11-11) we have
i k*x
9 f f C-k AJ^k) + k¥BJ(k)]e ~ ~ndk = nR - u(JAt)
(2n)2jl y ~ x ~ ~ ~n OB (u.76)
-09
n=l,2...N;
239
-------
•X
~n dk =
(271Y" '
" • (11-77)
'-00
n=l,2...N
where x = (x^,y ) are the observation sites described above, and u and v are
~n n n
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 8, 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 hyperplane 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
•s.
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
fluctuations in the flow in proportion to |jk| , 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(ik) ~ | jk | . 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 eQ;
p(v(x,JAt),X£D) =-j ~ ~ (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 Q 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:
" (ll-78b)
where E is the kinetic energy of the flow field v integrated over the entire
model domain 0. 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 Q 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
I
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 Stoke' s theorem, we can express (11-79) in terms of
the transform of the stream function:
ik L ik L k2 k2
A(k)[(e x x-l)(e y y-D( Vky )3dk dkv = c(t) (11-80)
KxKy x y
Thus, the initial version of processor Pll will be based on the equation set
Eq. (11-81)
Eq. (11-74)
Eq. (11-76)
Eq. (11-77)
(11-81)
244
-------
Each solution set [A (Ik), B^jk)] derived from (11-81) yields layer averaged
winds [, ] through the following equations, whose origin
is (11-29):
i k*x
= —!-»\\ C-kvAJ(jk) + k RJ(J$)]e ~ ~dk + u(JAt)
-09
00
i i ik-x
= - 2UCkxA(JS) + kyB(JS)]e ~ ~dk + v(JAt) (ii-82b)
-00
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 W3. (We should add that even
though the members of VT are not explicitly related to those of W , KAJ, 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), xeD, 0 < t< T, by "driving" the dispersion model with M wind fields
[, ] , xeD, J=0,l,. . . JMAX, u=l,...M. Each of the M flow
fields is created by selecting a and a at random from
the ensembles W 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
(VI_JV|_) f°r 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:
m=l,2,...M
mode 0 and
(11-83)
mode 1
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).
~in in
(2) Since the bottom surface of Layer 2 is pvs in mode 1 and p1 in mode 0,
the corresponding expression for the Layer 2 observed winds is
247
-------
- 1-
(11-84)
mode 0 and
x.t) . mode !;
(11-85)
(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 no ,:
"* UD , j.
(11-86)
mode ° only
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
-------
<"(&,»t)>OBt2
'm> ' I BF1....H
(11-87)
u(.xm,p.lc
~m VS-"r •' • ' BF1....H
mode 0 only
^•n't).^
In mode 0, p 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
g-ii 2 rnm (11-88)
where
W^rnm^ = ^xn"xm^ + ^rf^^"1' (11-89)
A similar operation yields 2- We can now define pseudo Layer 2 winds
over each surface station as follows
mode 0 only. (11-90)
OB)2 = 2v(xn,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, -J
1 VL I mode 1
1 = VL J
where yi and wi 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
P3(x,t)
p«(x,t)
p (x,t)
Description
top surface of Layer 3 in pressure
(mb) coordinates
same as above, except top of Layer 2
pressure (mb) at the virtual
Source
P8
P8
P7
surface
t) top surface of Layer 1 in pressure P7
(mb) coordinates
u(xm,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)
vCx^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.
<5(x,t)>2 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
<5(x,t)>,~ Same as above except applies to Layer 1 P8
v. Average east-west wind component in P7
Layer 1 during mode 0 (generally night-
time) hours (m/sec)
V| 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) ~
2 same as above, except north-south wind component
>$,t)>2 Vertically averaged wind (m/sec) in Layer 2 at grid
cell centered at x at hour t (east-west component)
*.,t)>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
-------
YYYYYY
o -n
•Si
0) U.
2-S g
O l
o)
0)
CQ u (/)
^ CI c
•* ss
•o o> ~
C CD w
O 5 03
ill
o o c
CM —
0.
-------
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 W, 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 ensemble averaged concentration. This requires
the model to be run a number of times, each time using a different flow field
[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.
ENVELOPE01 ENSEMBLE
AVERAGE PLUME (I)
CENTERLINE OF ENSEMBLE PLUME WIDTH (o),
AVERAGE PLUME CONTROLLED BY K%
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 K_ 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 = WA, where
w* = (HQg/0)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
where w is the average upward air speed within clouds and
H = H2 + H3 (12-5)
258
-------
where FL is the depth of the convective layer below cloud base and H3 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
V.J.V
Kjd.J.V
K30 <0 and <0 and
CTC = 0
-------
The relationships between H£ and H, 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 w (i,j,t ) required in the
expressions for K is an input parameter to processor P12; and w* should be
computed by
= (H2Qg/90)1/3 (12-7)
where H^ is given by (12-6a); Q is the surface heat flux, an input parameter;
-2
g = 9.8m sec ; and eQ is the surface temperature interpolated from the
surface virtual temperature measurements T as follows (compare with Eq.
7-108)
n=l
where
rn = C(iAx-xn)2 + (jAy-y^2]15 (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 \+, 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
•5. -
= h0/L
m
Of.
mj
2/3
"%
fm =
(l -164) , if 4 < 0;
i + 54 , if 4 > o
0.4 [0.4 + 0.6exp(44)], if | < 0;
o.4 [1.0 * 3.4| - 0.254], if 4 > 0.
(12-10)
(12-11)
(12-12)
and u* is the friction velocity. In Eqs. (12-10) - (12-12), the variables
o , u^, and 4 are all functions of (i,j,t ).
Next we approximate the time derivative of h by
h0(i,j,tffl) = [h0(i,j,tm) - h^lj.v
where At = 3600 sec is the time interval at which fv and other variables in
o
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, ai
8-13 of Part 1). .
h00',j,tm)
)] (12-1<
where erf is defind by
2
erf(x) =-£ fe"X dX. (12-16
V* J
awn fin h" hn
w+(U,V =-r2 «p(-A-) - f [i-erf(-a-)] (12-17
V2H 2awo Z V2awo
where aWQ, given by (12-10), and hQ, given by (12-14), are functions of
w_(i,j,tra) = n0(i,j,tm) + w+(i,j,tm) (12-18
v(i,j,tm) = u*(i,j,tm) (12-19;
The parameter v is the entrainment velocity of ambient air into plumes fror
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 4* defined in Part 1, Eq. 5-44. Interim values ¥' ol
this parameter were developed in Processor P8, but we must test whether the>
satisfy criterion (5-47) of Part 1, namely
262
-------
w A (l-0Tn)
v|»'< * + IU (12-20)
c c
where w+ and A+ are the parameters that we just computed in (12-17) and 12-15)
above, and where
V,. = - 4-r - ffl/A8 (12-21)
c l-ac 9
(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 A+) is an input to this
processor. Consult Table 12-2 for definitions.
To obtain the final values of ¥ 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,t) defined by (12-21) using the input fields w2(i,j,tj,
m ni
w (i,j,t ), etc. Next, compute
V'Ti i t ) = , . . . . — (12-22)
¥(1'J'V a^i.j.t^d'.j.tj ^Z2Z)
where w+ and X+ are from (12-17) and (12-15) and OVQ and a are inputs to this
processor. Finally, we have
Y(i,j,tffl) = «1n{f'(IJ.t^.rd.J.t^)} (12-23)
where V is the interim estimate of V that is input from processor P8. The
values of V 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 WQ, (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.
wn;i(i»J»'tlJ = vertical velocity at level z., during neutral and
unstable conditions (input from P8)
0.6w*(i,j,y if Ui,j,tm) < 0;
0 , if UU,tm) > 0
(12-24)
,-1
zjdj.v = cz^uiv - z1(i,j,tra.1)](At)-
where (12-25)
where At = 3600 sec is the time interval at which h, and hQ values are
available.
Let wm be the solution of the following transcendental equation:
i
_ 2 — —
ff , (w - *„,) w.:, wm- w^ ^w,.
J-i exp [- m D1 ] * -Si Ci-erf(-i«2i)] = ::£_£
2aJ Z /2awl 1-aTl
wl W (12-26)
264
-------
where erf(x) is defined by (12-16), ¥ is from (12-23), etc. (Note that awl,
*D1' ^' ac' *cf CTT1 and hence wm are a11 functions of (i,j) and tm). Now let
At each grid point (i,j) and each hour tm compute
(12-27)
+ C1 +
°T1 L
wl
i.j.y , if Ui,j,tm) > 0;
(12-28)
"m ^*
erf(-51) , if L(i,j,tm) < 0
'wl
nvs(i,J,tm) , if Ui,j,tm)> o.
(12-29)
, if U1,j,tn) < 0;
i,j,tm) , if
i,.i.t_) > 0.
(12-30)
265
-------
To evaluate the function erf in (12-28 and 12-29) in the limit as awl -»
0, the following procedure should be used.
if awl = 0, erf (——) = SIGN(A)-1 (12-31)
awl
This condition should be encountered only rarely since a - = 0 when L > 0 (see
Eq. 12-24) and in this case W, and W... 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
Description
Source
Stage
Output
Variable
Description
hxd.j.y
thickness (m) of
Layer 1 in cell
(i,j) at hour t
m
P7
K-iO'.Jit ) horizontal eddy
x m diffusivity (nT/sec)
in Layer 1, cell (i,j
hour t_
m
ho(i,j,t ) thickness (m) of P8
m Layer 2
h,(i,j,tm) thickness (m) of P8
J m Layer 3
Ko(i,j,tm)
^ m
»_
m
same as K, except
Layer 2 x
same as K, except
Layer 3 *
fraction (0
-------
Table 12-2 (Continued)
Input
Variable
Description.
Source Stage
Output
Vari ab 1 e
Description
uu.V
Obukhov length (m)
in cell (i,j) at
hour t
P4
WWO
(cont.
average speed (m/
of fluid moving u
relative to surfa
thickness (m)
of Layer 0 in
cell (i,j) at
hour t
m
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
vertical velocity
parameter (m/sec)
on surface HI
vertical velocity
parameter (m/sec)
on surface H,
w*(i,j,t ) convective velocity Stage K
m
scale (m/sec)
h-,(i,j,t ) thickness (m) of P7
m Layer 1 at (i,j,tm)
thickness (m) of
Layer 0
P7
vertical velocity
(m/sec) on surface
H, due to horizontal
divergence in the
flow field.
LCU.V
«Ki,J,t )
(i.j.t )
Obukhov length
scale (m) in cell
t.
m
cumulus flux
partition function
(dimensionless)
fraction (0<0 <1)
of sky coverea by
cumulus clouds in
cell (i,j) at hour
P4
Stage
WWO
RAW
w (i,j,t ) cumulus updraft
c m velocity scale
^m/sec)
dT1(i,j,t) fraction (0
-------
YYYYYYYYYIY
3
Q.
3
O
•a
c
CO
a
c
l
5! 8
n- 8
te w
O o
(A L.
(o a.
® a>
c'1
Is
2 S
^s *•
11
2
O)
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 p 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, N02, and
Oq 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) = Z T(x,n)r (cr,L)/I T(x,n) (15-1)
n=l n=l
where rn(or,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 rR 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 ~~**
103
103
10
103
VMM*
0
0
800
103
)
> use 0.0 if
[ L > 0 and
RH(x)> 0.9
)
-------
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 '
104 3
2.10d
0.0
100
100
400
200
200
300
400
300
1<>4 3
2.10J
0.0
250
200
400
400
300
1.5-103
1.5-103
1.5-103
"4 3
2.10^
0.0
300
350
Resistances for other pollutant species are not available in the current
literature. Deposition velocities for several pollutants other than S02 and
03 onto "vegetation11 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 Ar2(a) 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
00
°3
SO,
CO
NO
PAN
NO,
Ar2 C=r2(03)-r2(a)](sec/m)
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)
a 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, ML 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;a) = 0.4u*[ln(2/zo)+ 2.6 + 0.4u*r(x,t;a) - Y^"1 (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 AZQ - 10m (15-4)
275
-------
see Part I, Eq. 5-6b); and the function ?c is given by (Sheih et al., 1979)
{0.598 + 0.3901n (-2/L) - 0.090[ln(-z/L)]2}, 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.
.nput
'arametier
Description
Source Step
Output
Parameter
Description
3Hn(t)
L(x,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. -
£(x.,t;a) effective deposition
resistance (sec/m) of
cell centered at x
to deposition of species
a=0,, NO, N09, etc. at
hour t. *
z (x) 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.
s*
r(x,t;a) deposition resistance Step 1
(see Step 1 output).
p(x,t;a) effective deposition
velocity (m/sec) of
species a=02, NO, N02,
etc., in ceil centered
at x at hour t. Reference
level z=10m)
277
-------
Y
XI
-~ **
*
in
g
to «
(0 3 O
111
to
3
OJ
iZ
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 T-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 U« which is defined as follows:
_
0
0, otherwise
where ?Mn is the local emission rate of NO in Layer 0 and
Spec1es
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
H (CTT(faTl)Px "
ej 5 w.\. - w+\+(l-|)(l-a)Qx (16-4)
Ozone (03)
bn = *i + U-<5Vo)J> if uo -i-
0 0 2 (15'5)
3+ h + Q-a )fe 3- (1"q)" tQQ >T if U =0
nl u CTTO;teO 3 °
"12 '
b13 = 0 (16-7)
^ (l-aTO)(l-n)uSNO
9103 h1OQ + u)
280
-------
Nitric oxide (NO)
"ll • th* l * a-aTO)(e0- a. 0£Q )] (16-9)
l TO0
J. rNO
b!2 = H
b!3 =
QINO - SNO 'Vr 1f uo =
» if
— AHV
where
. £NO .
-S uS
Nitrogen dioxide (N02)
NO (16-13)
° (16-14)
'u ^' (16-15)
b 1 N02
bl1 - Hfew
•*•
b!2 - K^Tl^lm (16-15)
b13 * 0 (16-16)
TgiNO + SNO <03>i ' 1f Ut
= 1N°2 N0£ 3 1 o
1 ' n** U = 0
91N02 ' Uo °'
281
-------
where
NO
glN00 2 Sl * h.. o" + u) u5NO, . (16-18)
L 1 NQ2 2.
N02 (l-CTTO)(l-a)u
5 Sl + h,(f3Mn + u) (5NO
g1NO =
2
species y (excluding 03, NO and N02)
b22 = H
r (W2-W
c
b!2 = ^-"Tl^lm
b13 = 0 ' (16-23)
V i II ^w
«1X = Sl * h^u^p) uSX (16'24)
all species
b9i = \ (1'°n)(wi'l«'i»'Wm) (16-25)
^i n2 ii i im ui
b = - -La-a ) ^ (16"27)
D23 h,U CTCJ A6
Z (16-28)
282
-------
Ozone (03)
. if
b.3I =<
= 0.
(16-29)
b33 = H
> 1f Uo = 1:
n** = i r "0
30~ h
+ w:3L if u = o.
(16-30)
(16-31)
(16-32)
Nitric oxide (NO)
NO
9n =
if u =
93NO • 1f Uo =
where
* =
3NO
NO
(16-33)
(16-34)
(16-35)
(16-36)
(16-37)
283
-------
n** - H II r^
93NO h3H3U3C»
(16-38)
(16-39)
Nitrogen dioxide (N02)
NO,
(16-40)
b32 = E
(16-41)
(16-42)
C0o>i, if U
2 JHU2 31' 0
• 1f
(16-43)
where
2 "3
(16-44)
1 . N°9
** - I rij ii /» 2
(16-45)
9
3N0
J3
2 3N
(16-46)
284
-------
species x (excluding 03> NO and N02)
(I6.48)
(16-49)
where
l, if H, > 0;
(16-51)
0, if H, < 0.
™~
M= a [-f-r + fe/A83 (16-52)
i CTC
H3 = dH3/dt (16_53)
PREPARATION OF TERMS IN THE f-EQUATION
In designing the algorithm that treats the advection and diffusion
processes, we assumed that these phenomena could be expressed in the form (see
Eq. 9-5 of Part 1)
K =° (1654)
285
-------
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
_
sf
3
81 nV
8 .
v>
31 nV
(16-55)
where u^ = (a cos 4») and |j. = a are metric factors and a is the radius of
the earth. The fluctuation quantities u1, v1 and c' in (16-55) represent
deviations of the value of the given parameter at any point in space from the
local cell averaged value, denoted by < >.. 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.L
(16-56a)
a.
V
(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
3K. 3.
t
31nV
31nV. 3K. 3.
OA HA O/K J a/K
286
-------
2 2
a .
' — = ° (16'57)
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,4>) 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 ra s ; and
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, u,. ls tne
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 (Dj)"1 (16-58)
where *, is the latitude in radians of row J, but u. is a constant, namely
»l = 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 (\,) space for each layer j=l,2,3 and grid
point (I,J) as follows:
, (16-61)
,.., (16.-62)
KJjd.J) = &ix(I,J)]2KjU.J), (16-63)
KCia>J) = MJKXI.J). (16-64)
y j T j
The K*. and K*. fields should be recorded directly in the output file of the
xj yj
BMC for input into the model.
Step 2
Take the natural log of each cell volume V. in each layer at each grid
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
^ InV-d.J) - DXLVd.J.j) =
(16-65)
where
6\ = \ (35^) = model grid interval in A.
On the western boundary where 1=1, use
OXLV(I,J,j) = [ln(V(2,J))-ln(V.(l,J))]/6\ . (16-66)
and on the eastern boundary
DXLV(IMAX,J,j) = Cln(V..(IMAX,J))-ln(Vj(IMAX-l,J))]/6\ (16-67)
In the same manner define
4 InVjd.J) -.DYLV(I,J,j) = [InCVjCI.J+l)) (16.68)
where
°~ -id 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).
Form approximations of the derivatives of K*. and K*. as follows:
xj yj
9Kxi
^ - DXKXd.JJ) =[K*.(H-1,J)-K*.(I-1,J)]/(25X)
(16-70)
289
-------
3KIi
gj*l * DYKY(I,J,j) * EKJJ.(I,J+l)-K*j(I,J-l)]/(26)
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.
ueffd,J,j) = * - Kjj(I,J)*DXLV(I,J,j)
- DXKXd.J.J) (16-71)
veffd,J,j) = * - K*.jd.J)*DYLVd,J,j) (16-72)
- DYKY(I,J,j)
Values of ugff and vgff 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 d.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=tQ+NAt. Consider a given grid cell (I,J,j) at a given time step N.
Compute
290
-------
ueff (M.J.^AT (16-73)
d,J,j,N)AT (16-74)
where AT = At/n, say At = At/3 = 10 minutes. These increments define a point
(A.O in (A,) space, namely
A, =
(16-75)
Fit a biquintic polynomial to the (u «, veff) 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 (A,<|>) rather than (x,y). With this polynomial
and a linear interpolation in .time between NAt and (N-l)At, estimate the
values of ("eff, veff^ at tne P°int (^.p^) at "time NAt-Ai , i.e.,
ueff C^.^.J.NAt-Ai), Veff(\1,(01,j,NAt-AT).
Now compute the new point
\? = X* + [-u ff (\1,1,j,NAt-AT)AT].
Using the biquintic space approximation and the linear time interpolation
again to estimate u -- and v .- at (X-. 4>o' NAt-2Ai), compute finally
=) XR . = LAMBT(I,J,j,N) = \9 + [~u ,f(\?,4>0,
BTj Z eff Z *. (16-77)
= PHIBT(I,J,j,N) = «|>2 + [-veff(\2,4»2,j,NAt-2AT)Ai].
These are the coordinates of the "back track" point associated with cell
(I,J,J) at time NAt when the subinterval AT = At/3. This is the point
indicated by (x*, y*) in Figure 7-A1 of P7. Values of LAMBT and PHIBT should
be recorded on the b-matrix tape, i.e., the output file of the BMC, for all
I.J.j and all time steps N=1,...NMAX.
291 '
-------
The diffusivlty fields K*. and K*. should also be recorded on the
xj yj
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)
~± -a
f) to CO
•S £• ®
«2 t- '-
CQ ^
-------
Table BMC-1. Definitions of parameters in the
Model Input File (MIF).
Parameter
Definition
n
JTO
JT1
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
'1m
W
D1
"c
feMe
*
sx
52
08
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.
Colder, 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 pagesT
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 N09 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 Turlulerce 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 Inurucnons on the reverse before completing)
1. REPORT NO. 2.
4. TITLE AND SUBTITLE
A REGIONAL SCALE (1000 KM) MODEL OF PHOTOCHEMICAL
AIR POLLUTION Part 2. Input Processor Network Design
7. AUTHOR(S)
Robert G. Lamb
9. PERFORMING ORGANIZATION NAME AND ADDRESS
Same as Block 12
12. SPONSORING AGENCY NAME AND ADDRESS
Environmental Sciences Research Laboratory—RTF, 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 CODE
8. PERFORMING ORGANIZATION REPORT IN
10. PROGRAM ELEMENT NO.
-CDWA1A/02-1335 (FY-84)
11. CONTRACT/GRANT NO.
•
13. TYPE OF REPORT AND PERIOD COVERE
In-house
14. SPONSORING AGENCY CODE
EPA/600/09
15. SUPPLEMENTARY NOTES
i
16. ABSTRACT
Detailed specifications are given for a network of data processors and submodel
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 lay«
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 parameterizatic
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
3. DESCRIPTORS
13. DISTRIBUTION STATEMENT
RELEASE TO PUBLIC
b. IDENTIFIERS/OPEN ENDED TERMS
19. SECURITY CLASS (Tilts Report)
UNCLASSIFIED
20. SECURITY CLASS (This page)
UNCLASSIFIED
c. COSATI Field/Group
21. NO. OF PAGES
22. PRICE
EPA Form 2220-1 (9-73)
------- |