EPA-R4-73-025a
May 1973 Environmental Monitoring Series
Tests of an Urban
Meteorological-Pollutant Model
Using CO Validation Data
in the Los Angeles Metropolitan Area,
Volume I
Office of Research and Monitoring
U.S. Environmental Protection Agency
Washington. D.C. 20460
-------
EPA-R4-73-025a
Tests of an Urban
Meteorological-Pollutant Model
Using CO Validation Data
in the Los Angeles
Metropolitan Area,
Volume I
by
Joseph P. Pandolfo and Clifford A. Jacobs
The Center for the Environment and Man, Inc.
275 Windsor St.
Hartford, Connecticut 06120
Contract No. 68-02-0223
Program Element No. A11009
EPA Project Officer: Kenneth L. Calder
Meteorology Laboratory
National Environmental Research Center
Research Triangle Park, North Carolina 27711
Prepared for
OFFICE OF RESEARCH AND MONITORING
U.S. ENVIRONMENTAL PROTECTION AGENCY
WASHINGTON, D.C. 20460
May 1973
-------
This report has been reviewed by the Environmental Protection Agency
and approved for publication. Approval does not signify that the contents
necessarily reflect the views and policies of the Agency, nor does mention
of trade names or commercial products constitute endorsement or recommenda-
tion for use.
ii
-------
ABSTRACT
The urban boundary-layer model, described In a previous report (Pandolfo,
et al., 1971), was modified and used in forty test runs. Many of the runs
varied the meteorological input about a standard (observed) set. It has,
therefore, been demonstrated that an economical, objective, physically con-
sistent, and precisely specified (though with some arbitrary elements) pro-
cedure has been achieved for obtaining and predicting the three-dimensional
meteorological fields needed. In several of the runs, the input topography,
land-water distribution, and other physical characteristics of the underlying
surface was varied. The results demonstrate that ready generalization to
other regions can be expected.
The modelled region was simulated with relatively coarse (8-mile) grid
spacing. This is in contrast to other models which deal with pollutants
only, and which are based on two-mile grid spacing (Sklarew, et al., 1971;
Roth et al., 1971). Nonetheless, the temporal and spatial variations of air
temperature, humidity, and wind, are simulated with an encouraging degree of
realism. Temporal and spatial variations of CO are also simulated fairly
realistically, with somewhat less accuracy than in the model described by
Roth, et al. (1971), and with accuracy equivalent to that shown by Sklarew,
et al. (1971). It is reasonable to expect improved simulation accuracy with
the finer horizontal resolution used in these other models, and the perform-
ance of such simulations with this model is strongly urged.
iii
-------
TABLE OF CONTENTS
Section Title Page
1.0 INTRODUCTION 1
2.0 THE PLANETARY BOUNDARY LAYER MODEL 3
2.1 General Description 3
2.1.1 Designation of Land and Water Grid Points 3
2.1.2 The Terrain-Related Coordinate System 3
2.1.3 Symbols 5
2.1.4 The General Conservation Equation . 9
2.2 The Conservation Equation and Auxiliary Equations 9
2.2.1 The Atmospheric Layer 9
2.2.2 The Water Layer 11
2.2.3 Dependent Variables in the Soil 14
2.2.4 Boundary and Interface Conditions 14
2.2.4.1 Atmosphere , 14
2.2.4.2 Water 14
2.2.4.3 The air-water-soil interface 15
2.3 Vertical Diffusion Coefficients 20
2.4 Solar and Infrared Radiation 24
2.5 Computer Considerations 25
3.0 TEST SIMULATIONS WITH THE BOUNDARY LAYER MODEL 27
3.1 The Input Data 27
3.1.1 The Modelled Region 27
3.1.2 Physical Parameters of the Interface 33
3.1.3 Initial Fields 33
3.1.4 Boundary Conditions 41
3.2 The Experiments 43
3.3 The Results 43
3.3.1 Results of the Base Run (Run 37) 43
3.3.1.1 Spatial and temporal variations of the temperature 43
3.3.1.2 Spatial variations of wind and humidity 50
-------
TABLE OF CONTENTS (Continued)
Section Title Page
3.3.1.3 The eddy diffusivity 52
3.3.1.4 The vertical velocity 57
3.3.1.5 The pollutant concentrations 62
3.3.2 Sensitivity of the Results of Changes in the Terrain
Slope 64
3.3.2.1 A simulation with a level land surface (Run 38) 64
3.3.2.2 A simulation with doubled terrain slopes (Run 30) 72
3.3.3 Sensitivity of the Results to Other Modelling Changes 75
3.3.3.1 A simulation with the free atmospheric K . = 101*
cm2/sec (Run 39) mln 75
3.3.3.2 A simulation with doubled terrain slopes and
K . = Wk cm2/sec (Run 26) 75
mln
3.3.3.3 A simulation run with only land grid cells (Run 33) 78
3.3.3.4 A simulation run with observed winds replacing
model simulated winds (Run 31) 78
3.3.3.5 A simulation run with increased soil moisture
parameter (Run 19) : . 81
3.3.4 Comparison with the Results of Other Models 81
4.0 SUMMARY AND CONCLUSIONS 87
5.0 REFERENCES 91
6.0 ACKNOWLEDGEMENT S 93
APPENDICES
A Specifications for the Computer Program 95
B A Terrain-Based Coordinate System 135
C Atmospheric Radiation Computations for an Urban
Meteorological Pollutant Model 139
D The Non-Hydrostatic Model 159
E Advection Modifications 161
vi
-------
LIST OF TABLES
Number Description Page
1 Extinction coefficients in the ocean. 18
2 Heights (depths) of vertical grid levels in meters. 30
3 Physical parameters of the interface. 33
4 Initial vertical profiles of meteorological variables
used in the base run at grid point 3,2. 36
5 Horizontal gradients of meteorological variables input
to the base run. 37
6 Vertical profiles of wind derived from the San Nicolas/
Vandenberg RAWINS and surface observations (in cm/sec). 38
7 Horizontal gradients of wind as derived from the San
Nicolas/Vandenberg RAWINS and surface observations. 40
8a Test runs with K . = 103 cm2/sec , and variable
, min ,.
slopes. 44
8b Test runs with K' = 101* cm2/sec , and other model-
ling variations. mln 44
A-l Definitions of generalized coefficients. 122
A-2 Identification of quantities. 128
C-l Infrared diffuse transmission functions. 147
C-2 Transmission equations for clouds. 153
LIST OF ILLUSTRATIONS
Figure Caption Page
la A schematic diagram of a terrain-relative coordinate
system superimposed on a lake and surrounding land. 6
Ib A schematic diagram of the region shown in Fig. la as
it might be approximated through the finite-difference
techniques used by the model. 6
2 The locations of the horizontal grid used by Roth et
al. (1971) and of the coarser, smaller area grid used
in this report are shown. 28
vii
-------
LIST OF ILLUSTRATIONS (Continued)
Figure Caption
3 Modelled distribution of land and water grid cells in
the region. 29
4 Contours (labeled in feet) of smoothed elevation used
to calculate terrain slopes input to the model. 31
5 Values of terrain slope calculated from smoothed ele-
vation, where the upper values are the eastward compon-
ent of the slope (9£/9x); and the lower values (italic)
are the northward components of the slope (9£/9y). 32
6 Location of surface weather observing stations in
modelled region. 35
7 Observed and grid-point initial values of CO (ppm)
over the modelled region. 42
8-11 Regional distributions of observed surface tempera-
tures (bold face) and simulated temperatures at the
four-meter grid level (italic) for the base run (Run 37).
Both fields are shown at six-hourly intervals. 46
12-15 Time series of observed surface temperature at the
four weather stations in the modelled region; and of
simulated air temperatures at the 4m grid level and at
horizontal grid points in the vicinity of the station. 47
16 12P/29 September 1969. Vertical sections along the
.SW-NE diagonal of the modelled region showing simulated
temperatures (°C) and simulated winds. Note that winds
are plotted as on a conventional weather map. 48
17 18P/29 September 1969. Same as Figure 16. 48
18 OOP/30 September 1969. Same as Figure 16. 49
19 06P/30 September 1969. Same as Figure 16. 49
20-23 Maps of observed surface winds and specific humidity
(g/kg), and of simulated wind and humidity at the 4m
grid level, at six-hourly intervals. Observed humidity
values are in bold face, and the simulated humidity
values are in italic. 51
24 06P/29 September. Vertical sections along the SW-NE
diagonal of the modelled region showing simulated
mid-layer values of diffusivity (cm2/sec). 53
viii
-------
LIST OF ILLUSTRATIONS (Continued)
Figure Caption Page
25 08P/29 September. Same as Figure 24. 53
26 10P/29 September. Same as Figure 24. 54
27 12P/29 September. Same as Figure 24. 54
28 18P/29 September. Same as Figure 24. 55
29 OOP/30 September. Same as Figure 24. 55
30 06P/30 September. Same as Figure 24. 56
31 12P/29 September. Vertical sections along the SW-NE
diagonal of the modelled region showing simulated values
of the vertical velocity component calculated from the
horizontal wind divergence in terrain-parallel surfaces
( WD , cm/sec). 59
32 18P/29 September. Same as Figure 31. 59
33 OOP/30 September. Same as Figure 31. 60
34 06P/30 September. Same as Figure 31. 60
35 12P/29 September. Vertical sections along the SW-NE
diagonal of the modelled region showing simulated values
of the vertical velocity ( w , cm/sec). 61
36 06P/30 September. Same as Figure 35. 61
37-41 Maps of observed surface layer CO concentration (ppm)
and simulated concentrations at the 4m grid level.
Maps are shown at two-hourly intervals. 63
42-46 Time series of observed surface layer CO concentra-
tions (ppm) at monitoring stations; and of simulated
concentrations at the 4m grid level and at horizontal
grid points in the vicinity of the stations. 65
47-51 Time series of observed surface layer CO concentra-
tions (ppm) at monitoring stations; and of simulated
concentrations at the 4m grid level and at horizontal
grid points in the vicinity of the station. 66
52 06P/29 September. Vertical sections along the SW-NE
diagonal of the modelled region showing simulated CO
concentrations (ppm). 67
ix
-------
LIST OF ILLUSTRATIONS (Continued)
Figure Caption Page
53 12P/29 September. Same as Figure 52. 67
54 18P/29 September. Same as Figure 52. 68
55 OOP/30 September. Same as Figure 52. 68
56 06P/30 September. Same as Figure 52. 69
57-60 The time series of Figures 42 through 46 are repeated.
Also shown are simulated time series for Run 38, in
which the terrain slopes were set at zero everywhere
in the modelled region. 70
61-66 The time series of Figures 47 through 51 are repeated.
Also shown are simulated time series for Run 38, in
which the terrain slopes were set at zero everywhere
in the modelled region. 71
67-71 The time series of Figures 42 through 46 are repeated.
Also shown are simulated time series for Run 30, in
which the terrain slopes were doubled everywhere in
the modelled region. 73
72-76 The time series of Figures 47 through 51 are repeated.
Also shown are simulated time series for Run 30, in
which the terrain slopes were doubled everywhere in
the modelled region. 74
77-81 The maps of Figures 37 through 41 are repeated. Also
shown are changes in simulated CO concentrations in-
duced in Run 39, in which the free atmospheric minimum
diffusivity was increased to 101* cm2/sec . 76
82-86 The maps of Figures 37 through 41 are repeated. Also
shown are changes in simulated CO concentrations in-
duced in Run 26, in which the terrain slopes were
doubled and the minimum diffusivity was changed simul-
taneously. 77.
87-91 The maps of Figures 37 through 41 are repeated. Also
shown are changes in simulated CO concentrations in-
duced in Run 33, in which all horizontal grid cells
were designated as containing a land-air interface. 79
92-96 The maps of Figures 37 through 41 are repeated. Also
shown are changes in simulated CO concentrations in-
duced in Run 31, in which observed winds were used
throughout the simulation, rather than model-simulated
winds. 80
-------
LIST OF ILLUSTRATIONS (Continued)
Figure Caption Page
97-101 The maps of Figures 37 through 41 are repeated. Also
shown are changes in simulated CO concentrations In-
duced in Run 19, in which the surface moisture parame-
ter was increased to 0.2 in all land grid cells of the
model. 82
102-106 The time series of Figures 42 through 46 are repeated.
Also shown are simulated time series obtained by Roth
et al. (1971). 84
107-111 The time series of Figures 47 through 51 are repeated.
Also shown are simulated time series obtained by Roth
et al. (1971). 85
A Schematic diagram of vertical grid mesh. 96
B Schematic diagram illustrating the relationship between
the gradient in a surface parallel to the terrain
(z constant), and the gradient in a horizontal surface
(z constant). 136
C Vertical grid for atmospheric radiation computation. 141
El Advective transfers for conventional Eulerian scheme
with three wind trajectory azimuths. 162
E2 Concentration contour map for constant emission from
a point source using upwind differences on Eulerian
grid. Wind from ~275° (~W). 163
E3 Concentration contour map for constant emission from
a point source using upwind differences on Eulerian
grid. Wind from ~284° (~WNWbyW). 163
E4 Concentration contour map for constant emission from
a point source using upwind differences on Eulerian
grid. Wind from ~293° (~WNW). 164
E5 Concentration contour map for constant emission from
a point source using upwind differences on Eulerian
grid. Wind from ~304° (~WNW by N). 164
E6 Concentration contour map for constant emission from
a point source using upwind differences on Eulerian
grid. Wind from ~315° (~NW). 164
E7 Advective transfers for modified Eulerian scheme
allowing transfers to diagonal neighbors. Three wind
trajectory azimuths shown. 166
xi
-------
LIST OF ILLUSTRATIONS (Continued)
Figure Caption Page
E8 Typical allowable advective transfers for 16 point
(two-cell) modified Eulerian scheme. Five-wind tra-
jectory shown. 167
E9 Three typical trajectories in 16-point, nearest-grid-
point Eulerian advection scheme. 168
E10 Concentration contour map for constant emission from
a point source using Lagrangian advection scheme with
nearest grid-point interpolation (fully conservative).
Wind from 360° (N). 173
Ell Concentration contour map for constant emission from
a point source using Lagrangian advection scheme with
nearest grid-point interpolation (fully conservative).
Wind from 337P (WNW). 173
E12 Concentration contour map for constant emission from
a point source using Lagrangian advection scheme with
nearest grid-point interpolation (fully conservative).
Wind from 315° (NW). 173
E13 Excerpt from computer concentration map showing
south-central Connecticut. 175
xii
-------
1.0 INTRODUCTION
The urban boundary-layer model extensively described in a previous re-
port (Pandolfo, et al., 1971) has been tested with data from Los Angeles,
California. The tests were intended to estimate the degree to which observed
spatial and temporal variations of meteorological conditions and concentrations
of a stable air pollutant could be realistically simulated by the model.
The reader interested in more background discussion of the model is re-
ferred to that report. In Section 2.0 of the present report, we do, however,
give the set of basic equations again, since these have undergone two major
modifications and many minor modifications and corrections. In Section 3.0,
the results of several test experiments will be presented. These results are
summarized in Section 4.0, and possible steps to improve the performance of
the model are suggested. The Appendices give detailed computing formulas used
in the computer program, and a discussion and derivation of the modified co-
ordinate system used in the present version of the model. Appendix C, in
particular, gives the latest version of the formulas used for radiation cal-
culations in the model computer program. Much of the revision and further
development of these formulas carried out since our previous report were not
included in the scope of this contract (see Scope of Work), but were carried
out under separate sponsorship by one of our associates, Or. Marshall Atwater.
Appendix C is included here to provide the reader with complete documentation
of the model in convenient form.
Scope of Work
1) The Contractor shall continue the analysis of the simulation
results obtained under Contract No. CPA 70-62 and shall modify
the previous formulas and computer programs to accept, as a
model input an average terrain slope vector for each grid
square. This vector shall be used to calculate topograph-
ically induced vertical velocities and.to modify the calculated
values of incident solar radiation at the lower atmospheric
interface.
-------
2) The simulation model, as modified under 1) above, shall be
applied to approximately six days of meteorological and air
quality data selected by the Project Officer for the Los
Angeles Basin of California. Evaluation of the model shall
be based upon the ability to predict the following conditions:
a) The spatial and temporal development of the meso-scale
meteorological fields that are considered in the model.
b) The carbon monoxide (CO) concentration field when the
model equations for the pollutant field are integrated
simultaneously with those predicting the development of
the fields of the relevant meteorological variables.
c) The CO concentration when the observed meteorological
fields are treated as specified inputs to the pollutant
prediction model.
-------
2.0 THE PLANETARY BOUNDARY LAYER MODEL
2.1 General Description
2.1.1 Designation of Land and Water Grid Points
The set of vertical grid points associated with a grid cell in horizon-
tal space is referred to as a station. Four types of stations can be accom-
modated by the model. The stations may be characterized by their interfaces
and restrictions placed on the subsurface velocity components as follows:
o coastal corner stations have air-water interface and water velocity
components identically zero;
o coastal stations also have an air-water interface and the one water
velocity component normal to the coastal wall of the grid cell equal
to zero;
o water stations have an air-water interface but have no restrictions
on.the water velocity components directly associated with the sta-
tion designations;* and
o land stations have a land-air interface and zero subsurface velocity
components.
2.1.2 The Terrain-Related Coordinate System
The coordinate system used in the model is based upon a right-handed
Cartesian coordinate system with the positive x , y , and z directions
taken as east, north, and up, respectively; and with the plane 2=0 taken
to coincide with some prescribed absolute reference height (e.g., the height
of the water surface). The effects of topography are incorporated into the
equations by redefining the vertical coordinate to be
z = z - 5 (x,y) (1)
* Subsurface velocities at water stations might be restricted by the bot-
tom configuration of stations immediately surrounding the water station
(see Section 2.2.4.2).
-------
where €(x,y) is the elevation of the land or water surface above the
reference level, so that the surface z = 0 coincides with the interface
(i.e. air-land or air-water) (Gerrity, 1967).
The choice of the coordinate system leads to mean velocity advection
terms in the prediction equations that may be written in two alternative
forms
V'VX + w|j = V-VQX + w || . . (2)
where X is the dependent variable and the components of v X are defined
in Appendix B. The vertical velocity relative to an absolute reference
surface, W , is given by the vertical velocity relative to the terrain, w ,
plus the terrain induced divergence (see Section 2.2.1 and Appendix B).
In the present version of the model, the absolute reference level
z •> 0 is taken to coincide with the water surface. Thus, £(x,y) «= 0
for stations which contain an air-water interface, viz. coastal, coastal
corner, and water stations.
A schematic diagram of this terrain-relative coordinate system super-
imposed on a lake and surrounding land is presented in Figure la. The
method of solution for the equations used in the model, i.e. finite-differ-
ence techniques, do not resolve the degree of detail shown in Figure la
but treat the lake and surrounding land as composed of discrete "blocks"
centered about a grid point at which values of the dependent variables are
determined. The values determined at each grid point are therefore consid-
ered as representative of the. dynamic or thermodynamic state of the entire
block. Figure Ib is a schematic representation of the region shown in
Figure la as "seen" by the model.
-------
2.1.3 Symbols
The symbols used in this section are defined below with their units
(n.d. is nondimensional). Occasionally a symbol, used once, will be de-
fined when used.
a Surface albedo (n.d.)
A. Internal source (sink) for the dependent variable. Units are
determined by the i(th) variable.
B(x,y) Bottom depth for water stations (cm)(see Figure la,b).
c Specific heat of air (cal/g/°K).
Si
c Specific heat of subsurface (cal/g/°K).
5
c Specific heat of water (cal/g/°K).
D(x,y) Depth at bottom of boundary layer (cm)(see Figure la,b).
t . f.
E Surface evaporation (g H.O/cmVsec).
f Coriolis parameter (sec** ).
F Flux of pollutant at surface (g poll/g air/sec).
F Flux of pollutant at elevated level (g poll/g air/sec).
F+(z) Downward infrared flux at height z (£y/sec).
g Gravitational acceleration (cm/sec2).
h Priestley constant (n.d.).
H(x,y) Height at top of boundary layer (cm).
I(z) Incident solar radiation at height z (£y/sec).
k von Kartnan constant (n.d.).
k1 Kitaigorodsky constant (n.d.).
K Vertical diffusion coefficient for conductivity-diffusivity
(cm2/sec).
K. Generalized vertical diffusion coefficient (cm2/sec).
K Vertical diffusion coefficient for momentum (cm2/sec).
m
L Latent heat (cal/g).
-------
Figure la. A schematic diagram of a terrain-relative coordinate
system superimposed on a lake and surrounding land.
Figure Ib. A schematic diagram of the region shown in Fig. la
as it might be approximated through the finite-difference
techniques used by the model.
-------
M Halstead moisture parameter, 0 <. M <. 1 (n.d.).
P. Concentration of pollutant 1 (yg/m3).
P- Concentration of pollutant 2 (yg/m3).
P Precipitation rate (cm/sec).
q Specific humidity (g/g).
q Specific humidity at saturation (g/g).
8
R Incident solar flux at the interface Uy/sec).
8
R Artificial heat source due to combustion (£y/sec).
Ji
Ri Richardson number (n.d.).
s Salinity (g/g).
S. Solar radiation source term (°K/sec).
S_ Infrared radiation source term (°K/sec).
S Pollutant source term (g poll/g air/sec).
P
S Vertical gradient of the scalar value of the average orbital
W velocity in a turbulent wave (sec""-1).
t Time (sec).
T Temperature (°K).
T Temperature of rain (°K).
u Eastward component of velocity (cm/sec).
u Eastward component of geostrophic velocity (cm/sec).
O
v Northward component of velocity (cm/sec).
v Northward component of geostrophic velocity (cm/sec).
O
.Q 5 Wind speed at 19.5 m (cm/sec).
V Three-dimensional velocity vector (cm/sec).
» •
\?2 Two-dimensional (horizontal) velocity vector (cm/sec).
w Upward component of velocity relative to the terrain (cm/sec).
W Upward component of velocity relative to an absolute reference
surface (e.g. mean sea level).
-------
X Generalized dependent variable.
z Height relative to the terrain (cm).
2 Height relative to an absolute reference surface (cm).
Z Solar zenith angle (deg).
9/8x Partial derivative in the x-direction taken along surfaces
parallel to the terrain.
3/3y0 Partial derivative in the y-direction taken along surfaces
parallel to the terrain.
VX. Generalized three-dimensional gradient for dependent variable X.
relative to the absolute reference surface (see Section 3.0).
V0X Generalized three-dimensional gradient for dependent variable X^^
relative to the terrain parallel surface (see Section 3.0).
a Monin-Obukhov constant (n.d.).
T Dry adiabatic lapse rate (°K/cm).
6 Steepness of characteristic turbulent wave (height to length
ratio)(n.d.).
At Time increment (sec).
n Departure of the water surface elevation from a mean level (cm).
X Wavelength of turbulent wave (cm), the roughness height is 6 • 0 .
€ Elevation of the terrain above an absolute reference surface
(e.g. mean sea level)(cm).
p Density (g/cm3).
o Stefan-Boltzman constant (£y/°KI*/min).
o_ Density parameter of sea water at atmospheric pressure (g/cm3).
T Transmissivity of radiation (n.d.).
The following subscripts are used to further define variables.
a Atmosphere.
g Geostrophic.
I Interface.
-------
0 Relative to the terrain.
s Subsurface (water or soil).
w Water.
2.1.4 The General Conservation Equation
The basic equation for the dependent variables takes the form
aXi
at ' * *
-------
where all symbols have been defined in Section 2.1.1.
SM 1 3T I
& VRi.z) ll
L J
+ Sl + S2 + W
(6)
The humidity equation is written
I
] •
The equations for the prediction of pollutants are written
3Pi •»• 3
TT" B " V*V0Pi + 82.
i - 1,2 .
The formula for the vertical diffusion coefficient and the Richardson num-
bers are given in Section 2.3.
Geostrophic winds at any height in the layer are computed from the upper
boundary value and from the horizontal temperature gradients by use of the
integrated form of the thermal wind equation (see, e.g., Hess, 1959, eq.
12.11)
H
V" V ' T(H) r f /
J
z
and
H
T(z) g-T(z) f 1 3T
M - V^ J F a«
V«> - v (H) • ±£(. - *1L^ I ± |i dz . (10)
o o
The horizontal temperature gradients"in eqs. (9) and (10) are relative to
the absolute reference surface and are related to the horizontal tempera-
10
-------
ture gradient relative to the terrain through the equations (see the Appen-
dix) :
3T 3T_ 3T 31
and
_
3x 3xQ 32 3x '
3T 3T_ 3T i£
3y D 3y0 ~ 32 3y '
Vertical velocities at any height relative to the terrain are computed
from the horizontal velocity gradients by use of the continuity equation
for an incompressible fluid
W(z)
/f3u . 3v 1
(^+W
.The vertical velocity relative to the absolute reference surface is com-
puted from the equation
H
The vertical velocity is, therefore, assumed to be 2ero at the interface.
It is not necessary to specify w at either the upper or lower boundaries
(Gerrity, 1967).
2.2.2 The Water Layer*
The time dependent equations for the u and v components of the
current are written (Pandolfo, 1969) as
* Recall that in the water layer the absolute reference level z * .0 is
taken to coincide with the water surface and therefore z = z and
11
-------
and
3t ° 0 3z
iX = _ v-7 v + —
3t 0 dz
r
K (Ri
I m
r
I V fP-f
1 •>• \A\1
I m
J3u
>Z;J3Z
.]3v
>Z;j3z
+ f[v - v (z)]
O
- f[u - ug(z)]
(15)
(16)
The temperature in the water layer is computed from Pandolfo (1969)
3t
The salinity is computed from
-L
3z
9T
3z
+ S,
(17)
i£ „ _^.v s + A.
at v vos + 3z
r
ke(Ri,z)
Vi . *
3s
dz
(18)
with the diffusion coefficient and Richardson numbers defined in Section 2.3.
The geostrophic currents at any depth are computed from the analogue to the
thermal wind equation for the atmosphere, viz.,
and
ug(z=0)
= -v (z=0)
g
z
g /" I IP.
r y p 3y
5.
g / I IP
f y P 3x
z»0
dz
(19)
dz
(20)
The horizontal density gradients are computed from the horizontal temperature
and salinity gradients based on the relationships presented by Cox et al.
(1970).
In the RIGID LID version of the model, the terms u (z=0) and v (z=>0)
8 g
in eqs. (19) and (20) are prescribed as a function of time. The FREE SUR-
*
FACE version of the model predicts these values through the equations,
12
-------
and
V2-°> • - £ * (21)
Vz=0> '• £ U • <22)
Vertical velocities at any depth in the RIGID LID version are com-
puted from the continuity equation
•-/
WO) - v(>) - - / [^ + P da (23)
where W(2=I) « 0 (the vertical velocity at the air-water interface) .
Therefore, W , at the lowest grid level in the water, is not prescribed
but predicted through eq. (23) and is physically analogous to an open bot-
tom boundary in the water layer capable of absorbing gravity waves.
In the FREE SURFACE version, the vertical velocities at any depth
are computed by integration of the continuity equation which yields the
result
W(.) - w(z) - -+d2 (24)
B
where use has been made of the fact that W(z=B) » 0 . The vertical velo-
city at the interface computed from eq. (24) is used to predict the water
surface elevation through the equation
- W(a-I) . (25)
Equations (24) and (25) lead to a free water surface elevation and to the
v
introduction of undamped gravity waves at the water-air interface.
13
-------
2.2.3 Dependent Variables in the Soil
Temperature is the only time dependent variable in the soil and is
governed by the equation
f
where K (z) is the thermal diffusivity, which is specified as a function
of height. No other variables are simulated in the soil.
2.2.4 Boundary and Interface Conditions
2.2.4.1 Atmosphere
The dependent variables in the atmospheric layer are specified at the
upper boundary. The velocity field must be nondivergent but can be pre-
scribed as a function of time at the upper boundary. At lateral inflow
boundaries, the gradients of the dependent variables are prescribed rather
than calculated in the program.
2. 2. A. 2 Water
The temperature, salinity, and two pollutants are fixed at the bottom
water boundary. When the model includes land underlying the water, the
vertical temperature profile is specified, and the salinity and pollutants
are taken identically zero. The horizontal velocities are subject to kine-
matic boundary conditions with slip allowed at lateral boundaries. These
conditions are expressed as follows for each of the four types of water
stations:
u (z) = 0
w
v (z) = 0
5 D(x,y) <_ z <. 0
. . „/ x a coastal station
u (z) = 0 ; D(x,y) <. z <_ B(x,y)
w
v (z) «= 0 ; D(x,y) <. z <. 0
w
Coastal corner station (27a)
which is oriented parallel (27b)
to the latitude lines
(east-west coast)
14
-------
uw(z) = 0 ; D(x,y) <_ z <_ 0
vw(z) - 0 ; D(x,y) <. z <. B(x,y)
uw(z) = 0 ; D(x,y) <. z <. Bjj
vw(z) = 0 ; D(x,y) <. z <_ B^
For a coastal station
which is oriented perpen-
dicular to the latitude
lines (north-south coast)
For a water station
(27c)
(27d)
where BA is the bottom depth of the shallower of two stations immediately
north or south of the water station, and B' is the bottom depth of the
shallower of two stations immediately east or west of the water station.
Thus, B* and B' depend on the system of grid points to be chosen.
These definitions imply an additional constraint on a water station: that
a water station must always be surrounded by stations that have air-water
interfaces (i.e. coastal corner, coastal, or another water station). By
accounting for the bottom topography of the Immediately adjacent east-west
and north-south stations, the model is able to apply the proper kinematic
and slip boundary conditions for a variable-depth water basin.
Both versions of the model can accommodate open lateral water bound-
aries (therefore, the model can handle inflows and outflows to a lake basin).
The use of this type of boundary requires the gradients of the dependent
variables used in the prediction equations at an Inflow boundary station
to be specified rather than the calculated gradients. An additional re-
quirement is the specification of the water surface elevation outside all
open lateral water boundaries. This condition allows the water velocities
at open boundaries to be determined by the model.
2.2.4.3 The air-water-soil interface
Values of the dependent variables are obtained at the interface and
are highly dependent on the nature of the subsurface (water or soil). In
the case of water, several equation parameters are computed, while for soil,
15
-------
they are specified. These parameters Include albedo, a moisture parameter,
specific heat of the subsurface material, density of the subsurface material,
and the roughness height.
It is assumed that there are no discontinuities in velocity or temper-
ature. This leads to the interface velocity condition over water given by
#!>. • (Vl)w • ' (28)
For a soil subsurface, a no-slip condition is used, i.e. V, = 0 , at z •»
ZQ ; i.e., the velocity vanishes at some height above the interface. For
a water surface, the additional condition imposed is that the horizontal
shearing stress is continuous,
p K . K
pa m 3 •
The interface specific humidity qx is computed from
q_ = M(q) + (1 - M) q , ! (30)
x s «*•
where q. is the humidity at the first atmospheric grid level. The parame
ter is 1 for water and 0 for dry soil. For moist soil cases, M takes
on intermediate fractional values. Halstead et al. (1957) discussed the
empirical evaluation of M when it falls into this intermediate range.
The evaporation, E , is computed from
E - p Ke I1 • (31>
Ha^a 3z
The interface salinity, when required, is computed from
16
-------
E - P
<32)
The interface temperature is computed by an iterative solution to the
heat balance equation (see e.g., Pandolfo, 1969) by
o = R + F+(O) -[L-E](TT) - oil! - P c K6a
8 J- A c* cl a
(TT) + c p P (T - T_)
v Ix wkw r • r I
(33)
where the subscripts on the temperature gradients indicate derivatives taken
on the atmospheric side of the interface "(+)( or the water/soil side of the
interface (-). The assumption that the temperature is continuous at the
interface is applied in the calculation of the Individual terms in the ba-
lance equation.
The solar radiation term is computed from
Ra - (1 - a) 1(0) - TO 1(0) , (34)
8 S
with T the transmission of solar energy through the water interface (zero
O
for a soil interface) and 1(0) the solar energy at the surface. The solar
energy reaching the surface is computed by methods described in a forthcom-
ing report (see Section 2.4). The albedo is specified when the subsurface
is soil and is computed from
a - - .0139 + .0467 tanZ ; (35)
a >_ .03 when the subsurface is water (Sverdrup et al., 1942).
The amount of solar radiation absorbed by the water is used to complete
the source term for solar radiation in the water layer. The source term is
17
-------
computed from Hess (1959) as
1 3Ia
p c 3z
Kw w
(36)
According to Beer's law, the intensity of solar radiation at any level
is given by
I(z)
1(0)
8e°Z
, (z < 0)
(37)
where all radiation and absorption terms have been Integrated over the solar
spectrum. Because of selective extinction in the water, the value of K
is computed as an empirical function of depth from observations.
The absorption coefficients are computed or linearly interpolated from
data for the oceans, presented by Sverdrup et al. (1942) and Kondrat'yev
(1969), and are shown in Table 1. For Lake Ontario, an empirical relation
I
TABLE 1
Extinction Coefficients in the Ocean
Depth
< 2 m
2 m
5 m
10 m
20 m
50 m
100 m
based on data presented by
K «
w
. Kw
K - 1.034 - .5698 im
.6370
.3510
.2350
.1650
.1160
.0977
Beeton (1962) is given by
0.524 - 0.0l4z + .0034z
(38)
18
-------
where z is the depth in meters. This formula may turn out to be more gen-
erally representative of inland water than are the oceanic tabular values.
The infrared radiation flux received at the interface is calculated by
methods described in a forthcoming report (see Section 2. A).
The seventh term includes the effect of rain on the heat balance tem-
perature. The temperature of the rain is assumed to equal the height aver-
age of wet bulb temperature in the layer of air near the surface. While
the thickness of the layer appears to be realistic for the tropical oceans,
the validity of this thickness has not been tested in other regions.
The final term is an artificial heat ..source (heat added per unit area
per unit time) due to the activities of man by combustion of fossil fuels.
When the subsurface is soil, the albedo, the density, the specific
heat, the thermal diffusivity, and the moisture parameters are specified.
They are computed for water.
The value of pollutant concentration at the interface is calculated
from a flux value using the formula
(39)
where i •» 1,2 and F is the surface flux of pollutant as given in the
inventory. Note that K (z-r.i**) and Pi^2T+l>t:^ are calculated bY the
model.
The flux of pollutant into the atmosphere is not restricted to a
source at the ground and can occur at various elevations. Therefore, the
model has been developed to include elevated sources. The emission strengths
from the source inventory are divided and assigned, according to the nature
of the source, to points on the vertical grid. The strength at the inter-
19
-------
face is set by eq. (39). At the higher elevations, the pollutant is assumed
to be instantaneously mixed in the layer that is bounded by the midpoints
of the layers adjacent to the grid level. This results in a source term
S - —• . (40)
p pAz
2.3 Vertical Diffusion Coefficients
The vertical diffusion coefficients are internally computed in the
model as a function of the vertical gradients of wind, humidity, and tem-
perature. The method has been previously described by Pandolfo (1969) and
Pandolfo and Brown (1970).
When the subsurface is water, the formulas for the Richardson number
and vertical diffusion coefficients are modified to include mixing length
variations due to the presence of a "mean turbulent wave" at the interface.
The formulations used here are based on work by Kitaigorodsky (1961). The
parameters needed include 6 , the wave steepness (height to length ratio),
and X , the wavelength. The wavelength is found using the wind speed at
19.5 m (Pierson, 1964) by
X .- 28.03 x 10"** (V19 5)2 . (41)
Pierson also suggested that a constant steepness is appropriate for many
applications in "fully developed" seas. The value 0.055 is used in the
model when "fully developed" waves are indicated and is zero when the sur-
face is soil. .
Note that the modelling formulas shown are for "fully developed" seas
and are applied when over water to bring the parameterized wave state into
20
-------
immediate and complete adjustment with the low level wind. The method needs
revision for the simulation of duration-limited or fetch-limited cases. For
example, the wavelength could be reduced by a given fraction.
The relevant wave property used in the calculation of the Richardson
numbers and exchange coefficients is the vertical gradient of the scalar
value of the average orbital velocity in the turbulent wave given by
w
8
_x_
2ir
-z2w/X
(42)
In essence, this velocity shear, S , is- added to the shear of the mean
velocity when calculating Richardson numbers and vertical diffusion coeffi-
cients.
Richardson numbers are not computed at the first grid levels adjacent
to the Interface in either the atmospheric or water layers for computational
convenience. This does not lead to unrealistic results for a shallow enough
grid layer. Therefore, the Richardson numbers are assumed to be zero at
these grid levels in the following computations for the water case.
Whenever the subsurface is land, the term S is zero.
w
The Richardson numbers in the atmosphere are computed as a function
of height from
Ri
-2
(A3)
where Ri <. l/|a| . The Richardson numbers in water are computed from
Ri
8
-2
(44)
21
-------
The density gradient (3p /3z) is computed (e.g., Sverdrup et al. , 1942)
from
_3 SOT -3 f 3O«r af StJ-p a
10 ~- - 10 hrr!7+^- , (45)
dz I 3T 3z 3s 3
where the values of 3a_/3T and 3aT/3s are obtained by differentiation •
of an empirically obtained formula for a_ given by Cox et al. (1970).
The exchange coefficient formulas for the atmospheric layer are con-
sistent with the profile formulas discussed by Pandolfo (1966) , except that
the shear S is added to the shear |3V./3z| .
V» £*
Thus, for stable conditions, the eddy viscosity K is given by
Ri >_ 0 ,
where k = 0.4 ; and the eddy conductivity-dif fusivity is assumed to be
K - ' K . (47)
em
For lapse forced convection conditions, the formulas are
' (48)
(49)
and
Ke . = Km ll - aRi ,
where -.048 < Ri < 0 . For lapse free convection conditions, i.e. for
Ri <. -0.048 , the formulas are
22
-------
and
where
If B If I—I
m e IJ
«*/3 -2/3
C = 3k' h ' , (52)
and h Is the Priestley constant. See Pandolfo (1966) for more detail on
the formulas and constants used in lapse conditions.
Computed values for the eddy exchange coefficients for extreme stabil-
ity are meaningless whenever they fall below estimates for minimum values
in the atmospheric boundary layer (Pandolfo, 1969). The minimum values
used for the atmosphere are:
101* £ K <. 107 2 >. 100 m ,
Cfffl • .
A < K < 107 z < 100 m .
~~ e,m ~
where A •» A2 - z 2 for a land subsurface or A » (H)2 for water sub-
surface. Here ZQ is the roughness height, and H..- is the 1/3 highest
wave (see eq. 67). :
The exchange coefficient formulas used for the oceanic layer are the
modified Rossby-Montgomery formulas (after Kitaigorodsky, 1961) for stable
stratification: .
and
( aVo "\( m 1-3/2
(54)
where Rl ^. 0 . For these formulas, Kitaigorodsky assigned an estimated
value of k12= 0.02 . For lapse forced convection conditions (-.048 <
23
-------
Ri < 0), the formulas are given by eqs. (48) and (49) with the Kitaigorodskyr
constant (k..), replacing the von Karman constant (k).
For unstable stratification in the oceanic layer, the exchange coeffi-
cients are given by
and
fa) i . -1/6
v • .Kelt] W • <56)
where Ri <. -0.048 , and
_ _ *t/ 3 ~ fc/ 3 ^ _ _ k
C • 3k h ' (57)
Limiting values of the water exchange coefficient for stable conditions,
based on values that correspond to molecular viscosity, and conductivity-
diffusivity, are
K i. 0.14 cm2/sec
m
K £. 0.0014 cm2/sec ,
respectively.
2.4 Solar and Infrared Radiation
A discussion of the methods for obtaining the solar and infrared radi-
ation source terms in the model will be presented in Appendix C, written by
Dr. Atwater, under separate sponsorship.
24
-------
2.5 Computer Considerations
The model program is written to take any arbitrary number of horizontal
grid points in any rectangular array without further reprogramming (see Vol.
II for further detail). The practical restriction on the horizontal resolu-
tion then becomes the duration (and cost) of a single simulation run. With
the grid here used (8-mile horizontal grid mesh) we ran 40 test runs with the
hydrostatic model and several test runs with the non-hydrostatic model. It
was realized that the 8-mile grid resolution would fail to capture smaller
horizontal scales which contained important variability in the emission
source patterns, but it was decided that at this state of development it was
preferable to gain experience with many variations of the meteorological in-
put. An alternative course which was rejected would have allowed, within the
same computer budget, something like three test runs with the 2-mile horizon-
tal grid mesh more commonly used in pollutant models. Still another alterna-
tive is possible because of the program option which allows the meteorologi-
cal variables to be input. This would run the full model (meteorological
plus conservative pollutants) on a relatively coarse horizontal grid (say,
with 8-mile mesh spacing), and subsequently run a fine resolution (say, with
2-mile grid mesh) simulation, with the meteorological solutions of the first
run input and interpolated to the fine grid.
We have estimated the time and costs required to run the full model
(meteorological plus conservative pollutants) over a 24-hour simulation
period on the 2x2 mile horizontal grid. This would cover the Los Angeles
region with a grid network consisting of 25*25 (horizontal) *15 (vertical)
grid points. Our experience with the simulations run thus far yields a
running time of about 2 seconds (IBM 360-65) per mesh point for simulation
over a 24-hour period. For this computer system, therefore, a 24-hour simu-
25
-------
lation would take approximately five hours at a cost of about $1,000*.
Relative time on other systems available in the Hartford area are 1/3
(IBM 370/155; UNIVAC 1108) and 1/33 (CDC 7600). Relative costs* on these
systems are 0.58 (IBM 370/155), 0.85 (UNIVAC 1108): and .41 (CDC 7600).
* These costs are based on off-shift rates on systems available to CEM
locally in the Hartford area.
26
-------
3.0 TEST SIMULATIONS WITH THE BOUNDARY LAYER MODEL
3.1 The Input Data
3.1.1 The Modelled Region
The decision was made early in the course of the work to begin with
experiments performed on a relatively coarse horizontal grid for reasons
of economy. Also, as implied by eqs. (9) and (10)(Section 2.2.1), the hy-
drostatic approximation is applied in the present version of the modelj
while this approximation should be valid for wind systems resolved on hori-
zontal grids with mesh spacing of the order of 10 km, its validity becomes
questionable with much finer resolution. The computer program for a non-
hydrostatic version of the model has been checked out recently, but was not
available at the time of the work reported here. The modelled region was
chosen to cover a smaller area than that used by Roth et al. (1971), with
a horizontal grid cell dimension of 8 miles, rather than their 2 miles.
The two grid arrays, and their location in the Los Angeles region, are
shown in Figure 2.
The grid chosen resulted in the distribution of land and water grid
cells shown in Figure 3. The distinction between different types of water
grid cells is used in the application of barrier boundary conditions in
the prediction of subsurface flow. Of greater relevance in this project
is the basic distinction between land and water grid cells, and its effect
on meteorological conditions in the region. This will be seen later during
the discussion of the results of one experiment in which all grid cells in
the modelled region were treated as if they were land cells. It should be
noted that the greatest ambiguity in the definition of a grid cell type
arose in the case of grid cells 3,1 and 5,4.
27
-------
25
24
23
22
21
20
19
18
17
16
15
14
13
12
11
10
9
8
7
6
5
4
3
2
1
*
;=*=
SrS
'AH
'-f~-a
f^fj
PBjy
•i
*'•»
z
AJIfo
*-s^.-
-CA
^ -—
^
*AC
0 Vj
!ITA
KOU*
,
rs
. M
\
J>r(
1 * N
[ San Fe
LLLey
1
ION
TAI
~*"
anta
onlc
\
»
°C
CA
^S
i-t--^
a
%
\
\
J
\
Sty
•nan
-^-^-^
r
x
*
Los
• Tnt
Air
•t
'•••
^k\
ifcT"
2
tiffH^h
ft
Ho
lyw
Angel e
ernatloi
•port
5*
Nii^v;
- T\1
tjjp;
c,\XV>t
*•„..
"^^
)
>od
iai-
?,
W
if?
.wij)
3
-^
• c
^x
Pa
ownt
os A
Lon
!^
3 \^
sade
own
ngel
-Dow
g Be
t^^V
na«
es
ncjr
ach
^
1 — ^
s
s
AH £
4
•^
i
^
S
\
yl'.
Ni
•S_
--.
*ABR
^~
.^^
ffiz,
"X.-.
PU
^
n
-------
L
\
9\
\
9«
12.87 km
8 mi
L
L
L
L
V
L
L
L
L
9o
L
L
L
L
\
L
L
L
L
L
\
i Key
L Land cell.
Qn Coastal cell containing barrier walls
O to both N-S and E-W subsurface flow.
Q Coastal cell containing barrier walls
1 to N-S subsurface flow.
Q Coastal cell containing barrier walls
^ to E-W subsurface flow.
Figure 3. Modelled distribution of land and water grid cells
in the region.
29
-------
The model requires smoothed terrain heights and slopes to be input.
Averaged elevations over 64 square-mile squares, centered on and staggered
with respect to the basic grid array, were calculated from the 4 square-
mile average elevations given by Roth et al. (1971, Fig. C-9). Contours
of the smoothed elevations thus obtained are shown in Figure 4. Figure 5
shows numerical values of terrain slope for the land grid cells. In one
of the experiments to be described later, the land surface was treated as
if it were a level surface with zero elevation everywhere. In another, the
terrain slopes of Figure 5 were doubled.
The vertical layer modelled extends from 200 m below the air-soil(water)
interface to 1500 m above this interface. The grid levels used are listed
in Table 2.
TABLE 2
Heights(depths) of vertical
grid levels in meters
1500.00
900.00
650.00
450.00
300.00
150.00
75.00
25.00
4.00
0.00
- 0.05
- 0.10
- 0.15
- 0.20
- 1.00
- 10.00
- 50.00
-100.00
-200.00
30
-------
\ ',
\ \
\ I
WOO
-—''.
/
'& ,'
\
i
\ M
\
\00
"
( ]
\
Figure 4. Contours (labeled in feet) of smoothed elevation
used to calculate terrain slopes input to the
model.
31
-------
-4.18xlO~3
-2.58*10~3
-1.78xlO"2
2.3SxlO~2
•
•
•
4.97x10 3
0.7?*10~2
-2.77xlO~3
l.OOSxlO'2
2.985xlO~3
1.705*10~2
1.61xlO~3
2.80*20~3
•
1.04xlO~3
3.775*20~2
1.825xlO"3
0.98xlO~2
-2.605xlO~4
3.08*10~3
-0.57xlO~3
1.655*10~3
9
1.215xlO~2
1.20*10~2
0
0.715*10~2
4,995xlO~3
4.095xlO~3
4.735xlO~4
2.375*10~3
0.52xio"3
4.26xlO~4
-2.745xlO~3
0.575X10'1
0.63xlO~2
3.:7xJ^~5
0.785xlO~2
0.78*20~2
3.41xlO~3
3.41*20~2
1.16xlO~3
2.26*20~3
Figure 5. Values of terrain slope calculated from smoothed
elevation, where the upper values are the eastward
component of the slope (8£/3x); and the lower
values (italic) are the northward components of
the slope O£/3y).
32
-------
3.1.2 Physical Parameters of the Interface
A brief discussion of the considerations involved in prescribing these
parameters is given in Pandolfo et al. (1971, Section 3.2). In the Los
Angeles experiments, the values used in the base run (Run 37) are shown in
Table 3.
TABLE 3
Physical parameters of the interface
Parameter
Albedo (n.d.)
Soil heat conductivity (cal/cm sec °K)
Soil specific heat (cal/gm °K)
Soil density (gin/cm^)
Roughness height (cm)
Moisture Parameter (n.d.)
Artificial heat source (mJly/sec)
Land
.20
.005
.26
2.00
100.00
0.1**, 0.0
0.0
Water
*
*
1.0
1.0
*
1.0
0.0
* These values are internally computed for water grid squares.
** This value was used for grid squares 2,1 and 2,2, which cover
the seaward slopes of the Santa Monica Mountains.
In one of the experiments to be described later, the moisture parame-
ter was changed to 0.2 at all land stations.
3.1.3 Initial fields
Typically, the meteorological observations available for simulations
of this kind consist of a single vertical sounding, or at best, a small
number of such soundings, in a much denser net of surface observations.
The initial fields of meteorological variables are, therefore, constructed
in the computer program from input vertical profiles at an arbitrarily spe-
cified horizontal grid point (chosen to coincide as closely as possible
with the location of the observed sounding); and horizontal gradients (in
33
-------
terrain parallel surfaces) of these variables at selected heights above
the terrain. The construction is by linear extrapolation in the horizontal.
The program is written so that it is also possible to input meteorological
variables in this manner at other times during the period of simulation,
and to circumvent the model prediction of wind, or temperature, or humidity,
or any combination of these three variables.
In the Los Angeles observations, temperature and humidity profiles
were measured at Los Angeles International Airport twice daily at 14GMT
and 18GMT. Therefore, 14GMT (06 PST), 29 September 1969, was chosen as
the initial time for the experiments. The nearest wind soundings were
measured twice daily at OOGMT and 12GMT at San Nicolas Island and Vandenberg
Air Force Base, well outside the modelled area. Hourly surface weather
observations are available at four locations within the modelled region
(viz., USAS, Los Alamitos (NTB); Daugherty Field, Long Beach (LONB); Bur-
bank (BURK); and Los Angeles International Airport (LAX)). The location
of these stations within the modelled region is shown in Figure 6.
The base run (Run 37) used vertical profiles of temperature and humi-
dity derived from the LAX sounding at 14GMT, 29 September 1969. These are
shown in Table 4, and were input at grid point 3,2. Horizontal gradients
of temperature and humidity at the surface were obtained by fitting a plane
to the surface observations at the four stations shown in Figure 6, and
observations at two additional stations outside the region (viz., Van Nuys,
California, and Ontario, California) and are shown in Table 5.
The large-scale vorticity of horizontal wind near the surface was ob-
tained from measured winds in the same way, i.e. by plane fit to the obser-
vations, and is also shown in Table 5. The divergence of the initial wind
34
-------
BURK
LONB
Figure 6. Location of surface weather observing stations
in modelled region.
35
-------
TABLE 4
Initial vertical profiles of meteorological variables
used in the base run at grid point 3,2
z(m)
1500.00
900.00
650.00
450.00
300.00
150.00
75.00
25.00
4.00
0.00
- 0.05
- 0.10
- 0.15
- 0.20
- 1.00
- 10.00
- 50.00
-100.00
-200.00
u(cm/sec)
490.00
286.74
402.94
482.55
680.09
740.18
582.33
254.97
94.36
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
v (cm/sec)
-150.00
-171.83
-185.28
-171.29
- 58.34
42.14
90.18
79.97
34.46
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
T(°K)
295.51
300.71
302.67
302.55
298'. 78
296.82
295.29
292.91
292.35
291.11
303.16
303.16
303.16
303.16
289.98
286.86
284.66
282.66
281.56
q(g/kg)
3.650
4.101
4.486
5.974
9.442
10.881
11.071
11.085
11.087
11.087
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
36
-------
TABLE 5
Horizontal gradients of meteorological variables input to the base run.
( ) duf -!\ 3u/ -1\
1500
900
650
450
300
75
2
0
-200
0
0
X
0
0
0
0
0
0
.48867x10"'*
.15534X10"1*
X
.37217x10"'*
-.38835xlO~5
-.35599X10"1*
-.27710x10"'*
0
0
£<-"> £
.22654xlO~5
.6l489xlO~5
X
.14240X10"1*
-.16829X10"1*
-.18770X10"1*
.18854X10"1*
-(sec'
0
0
X
0
0
0
0
0
0
l) -jgrcm"1)
0
X
0
X
X
X
,12618xlO~6
0
0
•=— (°cm )
ay
0
X
0
X
X
X
.44229xlO~6
If*
-.71197xlO~7
-.19094xlO~6
X
-.20388xlO~6
-.19741xlO~6
-,17799xlO~6
-.32151xlO~6
0
0
3y
-.16181xlO~6
-.48867xlO~6
X
-.50485xlO~6
-.49515xlO~6
-.50809xlO~6
-.70550xlO~6
0
0
* Units in g(kg/cm) .
-------
field was set everywhere at 0 In the base run. The vertical profile of
wind for the base run was taken from the output of an auxiliary run which
had its initial profile derived from the 12GMT average of the San Nicolas
and Vandenberg RAWINS above the surface layer and which was run through a
24-hour simulation period. The initial vertical wind profile for this run,
derived from the RAWINS, is shown in Table 6.
TABLE 6
Vertical profiles of wind derived from the San Nicolas-Vandenberg
RAWINS and surface observations (in cm/sec)
z(m)
1500.0
900.0
650.0
450.0
300.0
150.0
75.0
25.0
4.0
12GMT/29 Sep.
u v
490.0 -150.0
415.0 -325.0
255.0 -230.0
110.0 -135.0
230.0 0.0
330.0 110.0
375.0 155.0
409.0 189.0
420.0 210.0
OOGMT/30 Sep.
u v
20.0 - 50.0
145.0 - 65.0
250.0 -160.0
310.0 -150.0
320.0 0.0
335.0 130.0
345.0 205.0
348.0 287.0
350.0 305.0
12GMT/30 Sep.
u v
350.0 -205.0
310.0 - 22.0
210.0 -105.0
135.0 - 15.0
85.0 - 10.0
35.0 - 10.0
10.0 - 5.0
90.0 55.0
165.0 95.0
The values of initial horizontal gradients at heights above the sur-
face layer must be guessed, because only one vertical sounding was available
within the modelled region. For the base run, the vorticity of the wind at
the upper boundary (1500 m above the terrain)was derived from the observed
average profile by assuming that the winds aloft were homogeneous in planes
parallel to sea level. The humidity at the upper boundary was assumed to
be homogeneous in surfaces parallel to the topography. The temperature was
also assumed, to be homogeneous in surfaces parallel to the topography, but
38
-------
only at grid heights of 300 m above the terrain. (This level defines the
top of the inversion in the LAX temperature profile.) Values at levels
intermediate between the surface and the level aloft at which horizontal
gradients were specified were obtained by linear interpolation in the ver-
tical. Many runs were made in which these guesses were varied; these will
not be described here, but suffice it to say that those that departed most
from this set of specified values gave grossly unrealistic results. These
values of gradients aloft are also shown in Table 5.
For one of the experiments, the model was not allowed to predict winds.
Instead, wind fields were derived from the RAWIN observations at OOGMT and
12GMT, 30 September and were input to the model. Wind fields at intermedi-
ate times are obtained in the computer program by linear interpolation in
time. However, in this run, the divergence derived from the observations
was kept in the model program. Vertical profiles of wind are also shown
in Table 6. Divergence and components of the vorticity as input to the
model for this run are shown in Table 7.
The observations of CO concentration available for the construction
of initial fields are even less extensive, in some ways, than the meteoro-
logical observations. There were no vertical profiles available at initial
time, and no measurements made over the water area. For this variable,
however, the computer program was written to allow the specification of a
separate vertical profile for each horizontal grid point. It was assumed
that, initially, the surface concentration at a given horizontal point was
valid at all vertical levels except at the upper boundary. The initial
surface concentrations at CO observing stations were obtained by averaging
hourly average concentrations centered at 0530 PST and 0630 PST, 29 Septem-
ber as given by Roth et al. (1971, Fig. 4). These were transferred to
39
-------
TABLE 7
Horizontal gradients of wind as derived from the San Nicolas-Vandenberg RAWINS and surface observations.*
12GMT/29 September
OOGMT/30 September
z(m)
— —— . TTT!
1500
900
450
300
75
2
0
-200
9u
9x
.19094X10"1*
.58253xlO~5
.14240x10"'*
.lossexio"4
-74434X10"5
-87084xlO~"lt
0
0
Su
9v
.48867x10~lf
. 15534x10^
.37217xlO~'+
-.38835x10"^
-.35500X10"1*
-27710X10"4
0
0
9v
9x
.22654xlO~5
.6l489xlO~5
.14240x10-^
-.16829X10"1*
-.18770xlO~'t
.13854xlO~I+
0
0
9v
9v
.41780xlO~5
.14563X10"1*
-.37864x10^
-.44984x10^
-50162xiO~'t
-ssoigxio'4
0
0
9u
9x
0
-.906xlO"5
-.971xlO~5
-.110x10"^
-.518xlO~5
.145X10"1*
0
0
9u
9v
0
-.227xlO~U
-.259x10^
.168x10^
-.104xlO~I+
-. 262x10^
0
0
9v
9x
0
.324xlO~5
.906xlO~5
.647xlO~6
-.168x10^
-.244X10'1*
0
0
9v
3v
0
.583xlO"5
-.518xlO~5
-.SSOxlO'4
0.570xlO~5
,248xlO~'t
0
o
12GMT/30 September
z(m)
—g-
1500
900
450
300
75
2
0
-200
3u
3x
0
.647xlO~5
.129X10"1*
.117xlO~4
.117x10^
-.217x10"^
0
0
3u
9v
0
.712xlO~5
.SllxlO'4
.272X10"1*
^esxio'1*
.493xlO~5
0
0
9v
9x
0
-.435xlO~5
-.149x10^
-.117xlO~lt
-.647xlO~5
-.517x10^
0
0
9v
9v
0
0.324xlO~5
-.350x10^
-.201X10"1*
-.777xlO~5
-.855xlO~8
0
0
* All units in sec
-1
-------
the model grid in the manner shown in Figure 7. Over the water, the value
of 3 ppm was estimated as being most consistent with the lateral boundary
conditions used by Roth et al. (1971). The horizontal gradient of this
variable was assumed to be zero at all inflow lateral boundary points
throughout the simulation. Again, there were several runs made varying
these input quantities, which are not described further in this report, and
which resulted in simulations that were less accurate than those shown here.
3.1.4 Boundary Conditions
The values of horizontal gradients of meteorological variables (see
Table 5) were used to estimate horizontal advection of these quantities at
inflow lateral boundaries throughout the simulation period. At the upper
and lower boundaries, values of these variables were prescribed as follows.
At the upper boundary (1500 m) the winds at grid point 3,2 were linearly
interpolated and extrapolated in time using the values shown in Table 6.
Nondivergent upper boundary wind fields were then linearly extrapolated
from this grid point value using vorticity components shown in Table 7. All
other meteorological variables are held constant in time at the values shown
in Table 5 at the upper and lower boundaries.
CO concentration was input as 3 ppm at the upper boundary. The hori-
zontal gradient was assumed to be zero at inflow lateral boundaries. At
the air-land interface, the vertical flux of CO was specified as in Roth,
et al. (Appendix A, and p. 22). Their fluxes, defined on the fine-grid
shown in Figure 2, were averaged to the coarser grid shown in the same
figure. At the air-water interface', the vertical flux of CO was assumed
to be zero.
Other boundary and interface conditions required by the model are
implicitly specified by the information given in Section 2.0 and Table 3.
41
-------
{RESD
11
,
11
WEST
v*'
\
\
\
3 \
!<
\
3
BURK
*14
t
14
7
LEUX
* .
10.5 11
10
~^J
3
12
CAP
*
5 '6
*VER
10.5
11
LONB
*8
8
3
11
ELM
*
11*
WHTR
* 11
11
1?
.
Key
.10
AZU*
10
11"
11
if
11
X.
% Observed
• Grid Point
Figure 7. Observed and grid point initial values of CO
(ppm) over the modelled region.
42
-------
3.2 The Experiments
There were a total of 40 test simulations carried out during the course
of the project. With minor input changes, this is the practical equivalent
of two or three simulations on the 2x2 mile mesh used in the other pollution
models referred to previously. The present version of the computer program
has not been designed for optimal efficiency, and the simulation of a 24-hour
test period takes, in general, about 20 minutes of CPU time on the University
of Connecticut shared-time 360/65 computer system. Computer costs contri-
buted to this project by the University of Connecticut exceeded the computer
budget requested from the sponsor by a significant amount.
A major problem in the preparation of this report was, therefore, the
selection and presentation of interesting portions of the test output. Those
runs selected for discussion in this report are listed in Table 8.
3.3 The Results
3.3.1 Results of the Base Run (Run 37)
3.3.1.1 Spatial and temporal variations of the temperature
The model simulation of the large spatial and temporal variations of
temperature typical of the modelled region and period, though not completely
accurate in detail, were quite encouraging. Realistic simulation of these
variations is important because of the indirect influence of the temperature
stratification on the evolution of pollution concentrations. It is also im-
portant as an indicator that the basic mechanisms affecting the distribution
and evolution of all the dependent variables (viz., the horizontal and ver-
tical advective transports, and the vertical mixing), in the presence of a
strong surface source, are simulated with some approximation to reality.
43
-------
TABLE 8
a) Test runs with K . = 103 cm2/sec , and variable terrain slopes,
min
Run Meteorological
Number Input
Properties of the
Modelled Region
Interface
Parameters
37
(base run)
38
30
As in Tables 4,5
As in base run.
As in base run
Smoothed terrain
slopes, both land
and water contained
in region.
Modelled region
contains level land
surface and water.
Modelled region
contains doubled
terrain slopes.
As in Table 3
As in base run.
As in base run.
b) Test runs with K
min
cm2/sec,and other modelling variations,
Run
NumbfciT
39
26
33
31
19
Meteorological
Input
As in base run.
As in base run.
As in base run
As in base run
but with observed
winds input from
Tables 6,7.
As in base run.
Properties of the
Modeled Reg.ion
As in base run.
As in Run 30.
As in Run 26, with
only land.
As in Run 26.
As in Run 26.
Interface
Parameters
As in base run.
As in base run.
As in base run.
As in Run 26.
Land moisture
parameter . 2 .
44
-------
Figures 8 through 11 compare the regional distributions of simulated
4 m air temperatures, and observed surface temperatures at the four stations
within the modelled region, and at six-hourly intervals after initial time.
Figure 8 compares the distributions at noon (PST), which is the time of maxi-
mum regional variation in this set of maps. The simulated range of temporal
and spatial variation of surface temperature is generally realistic, with
the largest discrepancy between simulated and observed temperatures evident
in the southern part of the modelled region. Grid cell 5,4, as has already
been mentioned, was treated as a land grid cell. This could account for the
over-prediction of surface temperatures at later times (especially at 18PST)
evident in the vicinity of stations LONB and NTB. In any case, use of a
finer horizontal grid mesh should be the first measure attempted in future
experiments to try to improve the simulation accuracy.
Time series of observed surface temperatures at the four observing sta-
tions within the modelled region are shown in Figures 12 through 15. Also
shown are time series of simulated 4 m air temperatures at the closest ad-
jacent grid points to each station. Again, the simulated time and space
variations are encouragingly realistic, with the most noticeable error evi-
dent in the overprediction of evening temperature in the vicinity of the
two southernmost stations, and the delayed time of maximum simulated tem-
peratures at LONB and LAX.
Simulated temperature distributions in the n-z plane, where n is
directed along the map diagonal extending from the southwest corner to the
northeast corner of the modelled region^are shown at six-hourly intervals
in Figures 16 through 19. The large-scale topographic variation is also
shown in these figures. The land-sea horizontal temperature contrasts of
the daytime hours are seen to extend to heights of about 400 m above the
45
-------
Figure 8. 12P/29 September
99
V
\LAX
\©
y
A
55 /
\.
BURK
96
65
75
75
/
-\/
65
96
94
83
LONB
•
64
93
99
89
90
\
S^X
93
99
92
92
89
Figure 10. OOP/30 September
72
v
\LAx
•\67
63 \
63 f
\
63
BURK
71
68
64
62
-^-J
63
69
*
59
55
LONB
64
57
70
68
NTB57
64
\,
69
71
m
69
69
57
Figure 9. 18P/29 September
•
92
^77
Aw
56 \
i
• >
o 5 {
\
•
53
BURK
73
75
71
s
^J
55
92
89
82
LONB
75
55
91
92
65
NTB65
©"
76
77V
•
32
S'S
89
84
V
Figure 11. 06P/30 September
56
\55
\LAX
.\6°3
52 \
52 /
\
61
BURK
©65
57
55
52
61
/
~w
61
64
65
62
LONB
52
74
Key
60
65
63
58
52\^
61
55
63
62
62
© Weather Station
98 Observed temperature at station
36 Temperature predicted by base run (37)
Figures 8-11. Regional distributions of observed surface temperatures (bold
face) and simulated temperatures at the four-meter grid level (italic)
for the base run ( un 37). Both fields are shown at six-hourly
intervals.
46
-------
-co
06B79
14R79 22B79
Time
06FH9 14FV29 22P/29 06P/30
Time
Figures 12 - 15. Time series of observed surface temperature at
the four weather stations in the modelled region; and
of simulated air temperatures at the 4 m grid level
and at horizontal grid points in the vicinity of the
station.
47
-------
.p-
oo
1800-
1400-
1600-
10 20 30 40 50 60 70 SO 90
D|km|
Figure 16. 12P/29 September 1969
0 'Ifl 10 JO 40 SO fO 70 80 »0
D{km|
Figure 17. 18P/29 September 1969
Vertical sections along the SW-NE diagonal of the modelled
region showing simulated temperatures (°C) and simulated
winds.. Note that winds are plotted as on a conventional
weather map.
-------
VO
1600-
1400-
4 .X^ I I I • I.I I I •
^ 10 20 39 40 50 fO 71 II 10
0 10 20 30 40 50 CO 70 SO SO
Figure 18. OOP/30 September 1969
Figure 19. 06P/30 September 1969
Vertical sections along the SW-NF, diagonal of the modelled
region showing simulated temperatures (°C) and simulated
winds. Note that winds are plotted as on a conventional
weather map.
-------
terrain, and to be replaced, as night wears on, by a more horizontally uni-
form pattern (in surfaces parallel to the terrain), with intense vertical
stratification. An interesting detail is the development of a pronounced
temperature maximum, extending through the lowest 400 m (above the terrain),
at the base of the steepest slope in the section, and a corresponding low-
altitude temperature minimum on this slope. This reflects a not unreason-
able diurnal variation in the location of the zone of maximum horizontal
temperature gradient, which is at the coast in the daytime, and at the base
of the San Gabriel mountains at night.
3.3.1.2 Spatial variations of wind and humidity
The model used for these runs calculated the vertical variation of
horizontal pressure gradients (geostrophic wind) from the imposed large-
scale, 'time constant horizontal temperature gradients listed in Table 5.
The use of more extreme (small scale) temperature gradients generated by
more local land-water and terrain differences appears to be unsuitable in
this portion of the calculations, as long as the hydrostatic approximation
is used (see Pielke, 1972, for a recent discussion of the effect of this
approximation in meso-scale numerical models). We have formulated and coded
a non-hydrostatic version of the model, but have not yet tested it. There-
fore, one would expect the simulated wind fields to exhibit smaller temporal
and spatial variations than those observed.
This is borne out by inspection of Figures 20 through 23. Although
the model captures the general variations of wind speed seen in the observed
surface winds (viz., in general,light speeds with a tendency toward higher
speed at noon, and lower speed through the evening and night hours, and
the general westerly component observed at noon), the model-generated winds
are much more uniform in direction than the observed winds.
50
-------
Figure 20. 12P/29 September
Figure 21. 18P/29 September
6.5
\?7.8
14. 8\
13. 3/
\
14.6
BURK
©9.1
/j >y\
o , o
/17.9
15.6
15.9
./
~^f
*I2.9
^9.4
/14.9
15.3
9.0
•16.8
12.7
/iz.s
/15.
15.9
16.3
NTEj
9.9
I2.2
13.5
xfs.
1 5. 6
/14.0
BURK
11.7
14.5
11.7
14.1
'14.9
LONB
12.8
/14.
14.7
15.6
NTE
7.4
/IS.O
'16. 0
IS. 2
Modelled • Observed
• 0-2 mph ®
•^ • 3-7 mph @—«:
' • 8-32 mph @ 1
Figures 20 - 23. Maps of observed surface winds and specific humidity (g/kg),
and of simulated wind and humidity at the 4m grid level, at six-hourly
intervals. Observed humidity values are in bold face, and the simulated
humidity values are in italic.
51
-------
It is also interesting to note in the cross sections (Figures 16
through 19) that there is a general tendency for the simulated wind speed
to increase significantly above the 4 m grid level. This suggests that the
transport of pollutant may be determined to a greater extent by the flow in
the layer between ten meters and a few hundred meters elevation, which is
not described by conventional observing networks, than by the winds at nor-
mal anemometer levels.
Comparison of observed and simulated fields of specific humidity shows
that the model (notwithstanding the coarse horizontal grid and subsequent
smoothing of terrain features) is able to generally simulate the large
contrasts in humidity evident between coastal portions of the region and
inland portions. This also suggests that the surface moisture parameter .
used over inland regions (see Table 3) is not drastically unrealistic.
3.3.1.3 The eddy diffusivity
The eddy dif fusivity-conductivi.ty formulas given in Section 2.0 may be
difficult for the reader to interpret in terms of the values and patterns
of diffusivity implied. Of the many ways to present the three-dimensional
fields of diffusivity developed by the model, we have again chosen to make
use of sequential (in time), n-z cross sections directed along the diag-
onal of the modelled region that extends from the southwest corner to the
northeast corner. These are shown at two-hourly intervals from initial
time (06PST, 29 September 1969; Figures 24 through 27) to noon of the same
day; and at six-hourly intervals thereafter to 24 hours after initial time
(Figures 28 through 30).
The somewhat arbitrary value for the minimum free atmospheric value
imposed, viz. 10 cm2/sec , dominates throughout most of the vertical slice
shown, at most of these times. One of the experiments to be described later
52
-------
(Jl
1600-
WOO-
40 SO
D{km|
Figure 24. 06P/29 September
90
1600-
1400-
1200-
1000-
800-
co
E
600-
400-
200-
5.1
10s
103
10'
10'
103
2xioL^
io»
^c-T^IO^vS
x
10
I
20
i
30
40 SO
D|km|
60
70
SO
90
Figure 25. 08P/29 September
Vertical sections along the SW-NE diagonal of the modelled
region showing simulated mid—layer values of diffusivity
(cm2/sec). The sections are shown at two-hourly intervals.
-------
1/1
-ts
o-
90
Figure 26. 10P/29 September
Figure 27. 12P/29 September
Vertical sections along the SW-NE diagonal of the modelled
region showing simulated mid-layer values of diffusivity
(cm2/sec). The sections are shown at two-hourly intervals.
-------
Ln
Oi
1600-
1400-
0-
90
Figure 28. 18P/29 September
Figure 29. OOP/30 September
Vertical sections along tbe SW-NE diagonal of the modelled
region showing simulated mid-layer values of diffusivity
(cm2/sec). The sections are shown at six-hourly intervals.
-------
1600-
1400-
90
D{km|
Figure 30. 06P/30 September. Vertical sections
along the SW-NE diagonal of the modelled re-
gion showing simulated mid-layer values of
diffusivity (cm /sec).
56
-------
shows that the accuracy of those simulated ground-layer concentrations, for
which validation data are available, is not very sensitive to an order of
magnitude increase in this value. Otherwise the simulated behavior of this
parameter, while not verifiable in detail, and (presumably) not accurate in
detail, is not unreasonable in general.
The major expected features are the development during the morning
hours of a layer of more intense mixing below the inversion base over land
which extends by noon to a height of 300 to 400 m above the surface; its
disappearance by late afternoon, and throughout the nighttime hours; and its
absence over water. A questionable detail is the development of an elevated
layer (based at 600 m and above) of intense mixing, separated from the sur-
face layer by an intermediate minimum of less intense diffusion over the
northeastern part of the cross section. It may be determined in the future
that this is not a realistic feature. For the present, as will be seen in
Sections 3.3.1.5 and 3.3.3.2, there is some indirect evidence in the accuracy
of simulated pollutant variations at AZU and BURK, that this is a realistic
feature of the simulated meteorology.
3.3.1.4 The vertical velocity
The simulated vertical velocity, like the simulated diffusivity, can-
not be directly compared to observations. Thus, we can only infer the
validity of the simulated patterns from the accuracy of simulation of other
variables, which are sensitive to errors in the vertical velocity to an
undetermined extent.
The formulation used in this model (see Section 2.0) separates the
calculated vertical velocity into two terms, one evaluated from the diver-
gence of the horizontal wind in terrain parallel surfaces; and one evalu-
ated from the scalar product of the local horizontal wind and the terrain
57
-------
slope. One would expect the simulated vertical velocity, therefore, to re-
flect the terrain scales defined by the horizontal grid; and to contain more
spatial variability when finer grids are used. Confirmation of this expec-
tation must await future experiments. Nonetheless, some readers may be in-
terested in the general characteristics of the simulated vertical velocity
fields in the coarse grid used in these experiments and we will briefly de-
scribe them here.
The divergence-induced component of the vertical velocity is shown in
the previously described cross section, at six-hourly intervals, beginning
at noon PST, in Figures 31 through 34. A general low-level pattern of sub-
sidence over the water, upward motion on the coast, subsidence at the foot
of the large slopes inland, and a more variable cellular structure on the
steeper slopes is evident. Vertical velocity contributions from this term
are of the order of one centimeter per second at a few hundred meter
elevation, and appear to be more intense in the daytime.
The terrain-induced component of the vertical velocity, when added
to the patterns above, produced total vertical velocities which vary spa-
tially in the same general manner, but which exhibit much less temporal
variability. Two of these fields are shown in Figures 35 and 36, at noon
PST on 29 September, and 06 PST on the 30th. This component, as would be
expected, has no influence over the water, has at least equal influence
over the coastal plain, and is definitely dominant over the steep slopes,
producing total vertical velocities of order 10 cm/sec in the northeastern
portion of the cross section.
The effects of the vertical velocity on the simulated thermal struc-
ture, and hence on the simulated diffusivity and subsequently the carbon
monoxide concentrations, should therefore be most pronounced in this part
58
-------
Ui
VO
1600-
1400-
1200-
1000-
800-
600-
400-
200-
30 40 SO 60 70 80 90
D|km|
Figure 31. 12P/29 September
1600-
«00-
1200-
1000-
800-
co
E
600-
400-
200-
5,1
3.3
2.4
1.5
1-1.4
•1.2 /
\l
t ' ' I'?
-0.6 \ I ' / I /
10 20 30 40 SO 60 70 80 90
Figure 32 . 18P/29 September
Vertical sections along the SW-NE diagonal of the modelled
region showing simulated values of the vertical velocity
component calculated from the horizontal wind divergence
in terrain-parallel surfaces (wp, cm/sec). The sections
are shovm at six-hourly intervals.
-------
1600-
1400-
i i i
10 20 30 40 50
D{km}
60 70 SO 90
Figure 33. OOP/30 September
1600-
1400-
10 20 30 40 SO 60 70 SO
SO
D|km|
Figure 34. 06P/30 September
Vertical sections along the SW-NE diagonal of the modelled
region showing simulated values of the vertical velocity
component calculated from the horizontal wind divergence
in terrain-parallel surfaces (WD, cm/sec). The sections
are shown at six-hourly intervals.
-------
1600-
1400-
10 20 30 40 SO 60 70 SO
0-
90
Figure 35. 12P/29 September
1600
1400-
1200-
1000-
800-
ro
E
400-
200-
600-
60 70 80 90
Djkrn}
Figure 36. 06P/30 September
Vertical sections along the SW—NE diagonal of the modelled
region showing simulated values of the vertical velocity
(w, cm/sec).
-------
of the modelled region. It has already been pointed out that the detailed
validity of the simulated diffusivity and vertical velocity fields must re-
main open to question until further validation has been carried out. It is
reassuring, however, that those portions of the grid which contain the most
intense and most complicated fields of these variables (i.e., the northeastern
portions) are also the portions in which the most accurate simulations of
ground-level concentrations of CO were produced. This will be shown in the
following section.
3.3.1.5 The pollutant concentrations
The CO concentration is the key variable simulated by the model. Its
spatial and temporal variations should be simulated with sufficient accuracy
to justify further development of the model in this application, and/or the
use of the temperature, wind, and diffusivity fields produced by this model
as input to other, more complex, models of pollutant transport, diffusion,
and chemistry (e.g., that of Eschenroeder and Martinez, 1971).
Maps showing observed ground-layer concentrations and simulated concen-
trations at the 4 m grid level, and at 2 hourly intervals from 08PST to 16PST,
are shown as Figures 37 through 41 inclusive. Most discouraging are the
large underestimates of CO concentration in the coastal and downtown dis-
tricts at 08 and 10PST, the. period of morning peak concentrations. Most
encouraging are the realistic simulations at Burbank, Azuza, and El Monte
throughout most of the period, in the region of most rugged topography, and
the region farthest from the upwind lateral boundaries of the model. Also
one might be encouraged by the fairly accurately simulated spatial variations
subsequent to 10PST, when it is remembered that these were obtained with a
relatively coarse horizontal grid.
The results are also presented as time series of observed CO concentra-
62
-------
Figure 37. 08P/29 September
Figure 38. 10P/29 September
•RESD
13.5
9.0
WEST
. #»
\ 4.0
\
3*.3\
'••<
3.0
BURK
#17.5
9.7
6.*2
LENX
* .
4.8
4.*7
^J
3.3
10.8
^7.3
17.5 '
#16
VER
7-7
LONB
#12
e'.7
3.4
77.2
ELM
*
1
s'.6
WHTR
70.7
. e'.2
\j
77.2
AZU#*
lOi
9.*9
' 77.0
70. 2
7.2
^.
•RESD
I
8.4
WEST
. *»
\ 4.0
\
\
A
3*0 ^
3.0
BURK
#10
«.8
•
4.8
LENX
5.5 3*. 2
3,7
/
"^J
3.5
8.8
CAP
# e-t
#12
VER
4to
LONB
#7.5
3.8
^-v— -v.
3.5
9.4
ELM
#
6,9
WHTR
» *"
5.J
4.8
\
9.8
AZU#*
9
7.2
.:.
8.2
/•
Figure 39. 12P/29 September
Figure 40. 14P/ 29 Setpember
rRESD
(.5
10.1
WEST
\«
,1
\
.\
2.9 \
A
':'<
\
3.0
BURK
#10.5
70.0
*
4.8
LENX
* .
3.2
3.7
^J
3.6
8*7
CAP
# 5.9
LS '
#5
VER
3.6
LONB
#6
3.5
-— i -^
3.6
8.7
ELM
*
•
5.3
WHTR
. #2
4.6
4.7
\
'\
8.0
AZU#*
(5
•
ff.S
5.3
5.0
*
4.0
•RESD
3
3.5
WEST
\ '*'
\'
3.0 \
\
2.9 /
2*9
BURK
#7
4.7
•
S.I
LtNX
»»
*
3.7
^J
3.*9
8.2
CAP
# 4.7
" '«
VER
3.9
LONS
3.6
3.9
5.2
ELM
#
3.5
WHTR
. #1
4.4
4.*7
3.3\
3.8
AZU#"
u
3.3
4.6
4.5
3*8
Figure 41. 16P/29 September
•RESD
5
3*. 3
WEST
. *«
\4
\
A
"\
\
2.9 y
V
2.9
BURK
3.3
•
5.6
LENX
# .
4.0
3.3
y
^J
4! 2
3.2
CAP
# 3.8
#3.5
VER
,
4.4
LONB
#5
-v— --_
4*2
3.7
ELM
f
•
3.3
WHTR
.*!
4.6
•
4.4
X
3.0
AZU#'
•
3.2
..
3.9
4.2
4,*0
Kejr
# CO Monitoring Station
9.8 Observed CO Concentrations at Station
9.8 CO Concentrations Predicted by Base Run 37
Figures 37 - 41. Maps of observed
surface layer CO concentration (ppm)
and sumulated concentrations at the
4m grid level. Mpas are shown at
two-hourly intervals.
63
-------
tions at ten monitoring stations at two-hour intervals. Time series of
simulated concentrations at the nearest grid point are also shown in this
set of figures. The same general features are evident in these time series
(Figures 42 through 51). In the time series, the simulated ground-level CO
concentrations are also shown for twelve hours beyond the time of the last
time of observation.
Simulated pollutant concentrations in the southwest-northeast vertical
slice used previously, and at six hourly intervals, are shown in Figures 52
through 56 inclusive. The unrealistic initial vertical distributions (viz.,
CO constant in the vertical to one grid level below the upper boundary), that
were used for want of better information, are shown in these figures to be
rapidly modified by the model. The simulated patterns are characterized by
isopleths of CO concentration that are oriented, in general, parallel to the
terrain. There is no reason to believe that these simulated patterns are
less realistic than the arbitrarily specified initial pattern.
3.3.2 Sensitivity of the Results to Changes in the Terrain Slope
3.3.2.1 A simulation with a level land surface (Run 38)
There are two experiments listed in Table 8 in which the input terrain
slopes were varied. Even a cursory analysis of the modelling equations shows
several direct and indirect linkages between this input parameter and the simu-
lated ground-level CO concentrations. It was the purpose of this experiment
to illustrate how these complex, and perhaps counteracting, influences affect
the pollutant concentration in this one case.
• Figures 57 through 66 repeat the time series of observed 4 m level CO
concentrations, and the base simulation (Run 37), already seen in Figures
42 through 51. However, in the new figures, time series of simulated CO con-
centrations for Run 38 are also added. Inspection of these figures shows,
64
-------
T»
8 ,
1"
8 •
OSP/2J UP, 29 IIP. 23 HP/3J
Time
HfVJt «P/»
_1 I
T..
I I
HP/29 UP/21
JlP/n IlP/2t
Time
Figures 42-46. Time series of observed
surface layer CO concentrations (ppm)
at monitoring stations; and of simulated
concentrations at the 4m grid level and
at horizontal grid points in the vicinity
of the station.
«p/» up/a
iimi
Time
linn OB» *rr»
-------
48
Run J7,
Grid Point 2,4
nan IORM IIP/»
Time
mr>'»
Figure 47 - 51. Time series of observed
surface layer CO concentrations (ppm)
at monitoring stations; and of simulated
concentrations at the 4m grid level and
at horizontal grid points in the vicinity
of the station.
-------
0-
1600-
KOO-
1200-
1000-
m
600-
400-
200-
0-
5.1
4.2
3.3
, 3.2
. .3.4
. 3.5
. .2.9
I. . 3.0
10
20 30 40 50
Dlkm|
60 70 SO 90
Figure 52. 06P/29 September
Figure 53. 12P/29 September
Vertical sections along the SW-NE diagonal of the modelled
region showing simulated CO concentrations (ppm), at six-
hourly intervals.
-------
00
1600-
1400-
10 20 30 40 SO fib 70 80 90
D|km}
Figure 54. 18P/29 September
1600-
1400-
10 20 30 40 SO 60 70 81
0-
90
Figure 55. OOP/30 September
Vertical sections along the SW-NE diagonal of the modelled
region showing simulated CO concentrations (ppm), at six-
hourly intervals.
-------
1GOO-
1400-
1200-
1000-
800-
1
<0
E
coo-
400-
200-
0-
5.1
T3.0
• -3.0
• 3.0
. 3.0
_ 3.0
..2.9
2.9
• -3.0
• 3.0
••3.1
-•3.0
..3.0
2ft
1.5
3.0
i
20.
30
40 50
Dfkm|
70 SO $0
Figure 56. 06P/30 September. Vertical sec-
tions along the SW-NE diagonal of the
modelled region showing simulated CO con-
centrations (ppm).
69
-------
1"
I ,
I .
1
S <
mm umt arm
Time
ll
tam
8 ,
1IFV2I
Time
2JP/H OlMO
Figures 57 - 60. The time series of
Figures 42 through 46 are repeated.
Also shown are simulated time series
for Run 38, in which the terrain
slopes were set at zero everywhere
in the modelled region.
ami m*n am
lima
asm ' •ma tarn
-------
I «
!
8 ,
I
Q.
8 «
_ _ _ _ _
Ufia 5*78 ISvs 55H35 53B35 uma
Time
*j 'Old Pout J.2
I ft--* »;„ }l
66
Run )>.
G<-U Point 4.3
ft--* Run U
g I _ I _ I - 1 - 1
MOD w/n wn IIPJI nfWJ nf>Ju
Tune
Uftl) HP/I! UKrH llfira ami U^M OIF>M
Figures 61 - 66. The time series of
Figures 47 through 51 are repeated.
Also shown are simulated time series
for Run 38, in which the terrain
slopes were set at zero everywhere
in the modelled region.
-------
first, that the terrain effects are minor at most locations during the morn-
ing and early afternoon hours. The presence of terrain slopes appears to
decrease the afternoon and evening concentrations significantly only at CAP,
VER, and AZU. It is noteworthy, however, that even at the other locations
where the effect is much smaller, it is, in general, of the same quality
(viz., the presence of the terrain slopes tends to decrease afternoon and
evening concentrations of pollutant in this case).
3.3.2.2 A simulation with doubled terrain slopes (Run 30)
It is reasonable to expect, with an increase in the horizontal resolu-
tion of the modelled grid, that larger terrain slopes (less smooth topography)
will be input. One might expect more extreme meteorological behavior with
such input: and, of course, one might also be interested in the consequences
as reflected in the simulated pollutant concentrations.
Run 30 was,therefore, run with doubled terrain slopes. The simulated
meteorology was well behaved, and the resulting surface temperature predictions
(not shown here) were slightly, but noticeably, different in the northeast
portions of the region. Figures 67 through 76 again repeat the time series
of CO concentrations shown in Figures 42 through 51, but also show simulated
time series for Run 30.
Again, the terrain effect is small at most places and most times. And
again, the effect of more rugged terrain is to decrease the evening concentra-
tions, where it is more pronounced (viz., CAP, VER, ELM).
However, it should be noted that at AZU, Run 30, with the doubled ter-
rain slopes, produces higher afternoon and evening concentrations. These
higher afternoon concentrations, as at BURK, are more consistent with the
observations. This reinforces the expectation that increasing the horizontal
grid resolution, with a consequent increase in the resolution of the terrain
72
-------
U)
*»
Time
EFya «mo «#»
tsera ami imn up/a wva OKJ> ttmt
Figures 67 - 71. The time series of
Figures 42 through 46 are repeated.
Also shown are simulated time series
for Run 30, in which the terrain
slopes were doubled everywhere in
the modelled region.
tutu itf>/it MF*» uf»:i
Time
-------
o
U I
I
a
8
U •
turn UBB jip/n i«p/»
Time
I"
8 i
«m» Kfws up/a ima tumt otwo
Time
Figures 72 - 76. The time series of
Figures 47 through 51 are repeated.
Also shown are simulated time series
for Run 30, in which the terrain
slopes were doubled everywhere in
the modelled region.
itmo torn vsta \mn ami am urns arm
Time
-------
slopes, will lead to more accurate simulations.
3.3.3 Sensitivity of the Results to Other Modelling Changes
3.3.3.1 A simulation with the free atmospheric K . = 101* cm2/sec
- - min -
Figures 77 through 81 repeat the distributions of observed low-level
CO concentrations at two hourly intervals, and the simulated patterns of CO
concentration at the 4 m grid level obtained in the base run (37) . Also
shown in these figures are the changes in simulated values obtained in each
grid cell by increasing the free atmospheric value of K to 10Vm2/sec .
min
It is evident from these maps that this modelling change had, in general,
an insignificant effect on the accuracy of the simulated concentrations, dur-
ing the period for which validation measurements were available.
3.3.3.2 A simulation with doubled terrain slopes and K . ° IP1* cm2/sec
(Ruil 26) = - min
Figures 82 through 86 also repeat the maps of observed and base run simu-
lated CO concentrations first shown in Figures 37 through 41. Also shown on
these maps are the changes in simulated values obtained with the simultaneous
change of two modelling parameters which were changed individually in Runs 39
and 30 described above. The results are as would be expected from the results
of the individual runs. However, it is more easily seen on these maps that
the improvement in simulation accuracy obtained by using less smooth terrain
slopes is confined to the northern boundary of the modelled region as might
be expected.
Less easily explained is the fact that this improvement is apparent only
in the early afternoon.
75
-------
Figure 77. 08P/29 September
Figure 78. 10P/29 September
•RESD
13.S
9.0
(+0)
WEST
. *ir
\ 4.0
\(+0)
\
A
3.3 \
(+0) \
3.0 S
(0) \
X
•
3.0
(0)
BURK
#17.5
9.7
(+0)
5.2
(+0)
LENX
* .
1.8
(-0)
I*
4. J
(-0) x
^J
•
3.3
(-H))
JO.C
(+0)
CAP
* 7-3, „
I7.S * <-°
#16
VER
•
7.7
(-0)
LONB
#12
6./(-0)
-~v— ^
3.4
(-K))
(-K))
JJ.2
ELM
#
s
6.6
(-0)
WHTR
. #12
JO.l
(-0)
•
8.2
(-0)
N
4^v
(•fo) \
11.2
AZU#"(+0)
10.5
•
9.9
(-0)
>
11.0
(-0)
•
JO. 2
(-0)
7.2
v (-0)
:RESD
7
' •
(40)
WEST
\4.0
\40)
3.*2 \
(+0) \
• A
3.0 S
(0) V
3.0
(0)
BURK
#10
5.*6
(40)
•
4.8
(40)
LENX
5-5 3*2
(+0)
3.7
^J
*
3.5
(40)
6* 5
(40)
CAP
#6.3
•(-0)
#12
VER
'•
4.0
(+0)
LONB
#15
3.5(40)
3? 5
(40)
(+0)
•
ELM
#
6*9
(-0)
WHTR
. #12
5.5
(+0)
4.6
(40)
X
(40) \
9.8
AZU*'(-K»
9
7.2
6.9
(+0)
6.2
(+0)
•
4.2
v (+0)
Figure 79. 12P/29 September
•RESb
4.5
10.1
(-0)
WEST
. #«
\(+0)
\
A
2.9\
(40) \
2*9 y
(40 A
>.
3*0
(40)
BURK
#10.5
io.o
(-0)
4.6
(+0)
LENX
# .
3.2
(40)
3.7
(40) y
^vj
3! 5
(40)
6*7
(41)
CAP
# 5.9 *
15. '(40)
#5
VER
.
3.5
(40)
LONB
#E
3.5 (40)
— Ir^^,
3*6
(40)
8.7
(+1)
ELM
5.3
(+0)
WHTR
. *2
4.6
(40)
4.1
(40)
\
(40) \
6.0
AZU*'(40)
6.5
6.5
(40)
.
5.3
(+0)
5.0
(40)
4.0
(40)
Figure 80. 14P/September
hRESD
3
3.5
(-0)
WEST
\4*.2
\(+0)
\
A
3.0 \
(-to) \
2.9 /
2.9
(40)
BURK
#7
4.7
(-0)
5.7
(•HO-
LE NX
# .
3.5
(40)
3.7
(40) /
"^J
3.9
(40)
6.2
(40)
CAP
#3
VER
•
3.9
(40)
LONB
#5
3.6(40)
~T/
3.9
(40)
(43)
5.2
ELM
#
3.5
(40)
WHTR
. #1
4.4
(40)
4.7
(40)
. .
3.5
AZU#'(40)
3.3
(•K»
•
4.6
(40)
4.5
'(40)
3.6
v (40)
Figure 81. 16P/29 September
iRESO
3.3
(0)
WEST
#4!
\ 4.4
\
A
3.A
\
• A
2.9 /
(+0) \
V
2*9
(+0)
BUIRK
3.3
(0)
t
5.6
(40)
LENX
* .
4.0
(40)
•
3.3
(40) y
•vj
4*2
(40)
3.2
(•H»
CAP3 .
* .CO)
#15
VER
.
4.4
(+0)
LONB
#5
3.9 (40)
— v-'^—
4.2
(40)
(40)
3.7
ELM
#
2
^
3.3
(40)
WHTR
.*'
4.6
(+0)
4.4
(+0)
\
(40) X,
3.0
AZU#* (40)
^
3.2
(40)
.
3.9
(40)
4.2
(40)
4.0
^ (+0)
# CO Monitoring Station
18 Observed CO Concentrations at Station
9.5 CO Concentrations Predicted by Base Run 37
(0) Change In CO Concentration for Run 35
Figures 77 - 81. The maps of
Figures 37 through 41 are repeated.
Also shown are changes in stimulated
CO concentrations induced in Run 39,
in which the free atmospheric minimum
diffusivity was increased to 10l*cm2/sec.
76
-------
Figure 82. 08P/29 September
•RESD
IIS
9*0
<-K»
WEST
#li
\4.0
\(-K))
\
A
3.3 \
(0) \
*
3.0 /
(o) \
\
•
3.0
(0)
BURK
*ns
•
a. ?
OK))
•
6.2
(+0)
LENX
* .
4.9
(+0)
•
4.1
<-K»
^
•
&J
10, e
Oo)
CAP
* 7.3
I7.S * (o)
*16
VER
•
7.7
(0)
LONB
*12
6. 1 (-K))
-~V -s_
fe!
.00)
11.2
ELM
*
1
•
8.6
(-0)
WHTR
. *12
10.1
(-0)
8.*2
(+0)
\
11.0
(-0)
10.2
(0)
&!
^.
Figure 83. 10P/29 September
•RESD
I
8*4
01)
.IS
\4.0
\OO)
A
OO) \
i.o /
(0) C
3*0
(0>
BURK
*10
8.*8
(•TO)
•
4.8
OO)
LENX
*
U 3.*2
(-K))
3*1
"V
8*8
OO).
CAP
* e-.*(-o:
VER
4.0
Oo)
Lf}|
3.8*OO)
3*5
(-K))
(-K))
9.*4
ELM
*
•
6.9
(-2)
WHTR
£.5
OO)
4*8
OO)
V
^^
9.8
AZU*' (-K>)
•
7.2
(+D
6.0
(•to)
6*2
(-to)
•
Figure 84. 12P/29 September
Figure 85. 14P/29 September
•RESD
(.S
10.1
Oo)
WEST
V3**
\*0)
2.*9\
(+0)\
£<
\
s'.o
(0)
BURK
*10.S
10.0
(+0)
4.8
(+0)
LENX
* .
3.2
(-0)
3*.l
OO) y
"^J
3*6
(•w)
a'?
CAP
« {>.9(-o)
«•$
VER
•
3.6
OO)
LONB
*6
3. 5* OO)
J.*6
8*1
ELM
8.3
(-0)
WHTR
4.*6
(-0)
4.1
(0)
V
8.0
AZU** OO)
6.S
0.5
(-0)
5*3
(-0)
5.0
(-0)
4\0
•RESD
•
3.5
(-0)
WEST
3.*0\
\
2.9 y
O) \
^
2.9
OO)
BURK
O3)
5.3
(0)
LENX
*
» 3.'5
(-0)
3.1
^vj
3.9
OD
6.2
01)
CAP
£ ^(+1)
VER
3.9
OO)
LONB
*s
3.SOO)
3.9
(-3)
£.2
ELM
*
3.5
(+0)
WHTR
. *l
4.4
OO)
4.1
(-0)
3.3\
3.8
AZU** OS)
IS
3.3
(+1)
4.6
(-0)
4.5
(-0)
3.8
v (-o)
Figure 86. 16P/29 September
rRESD
S
3.3
(0)
WEST
. **!
\(+o)
\
A
3.1 \
OD \
'
(i/
V
2*9
(-K))
BURK
*s
3.3
(0)
5.6
(+0)
LENX
* .
4.0
Oo)
3*. 3
00) y
~vj
4*2
01)
3.2
02)
CAP
? 3« 8 (+3
*15
VER
(-0)
LONB
#5
3.9* (+0)
—V-^-s.
4! 2
(-K))
(+4)
3.1
ELM
*
2
L ^
3.3
(+0)
WHTR
.*>
4.6
(-W)
(+0)
\
3.0
AZU** 03)
^
3.2
(•w)
.
3.9
OD
4.*2
(+0)
4.0
[^ (-0)
* CO Monitoring Station
9.S Observed CO Concentrations at Station
9.8 CO Concentrations Predicted by Base Run 37
(0) Change In CO Concentration for Run 26
Figures 82 - 86. The mpas of Fig-
ures 37 through 41 are repeated.
Also shown are changes in simulated
CO concentrations induced in Run 26,
in which the terrain slopes were
doubled and the minimum diffusivity
were changed simultaneously.
77
-------
3.3.3.3 A simulation run with only land grid cells (Run 33)
Figures 87 through 91 also repeat the maps used above. Also shown in
these figures are the changes produced in the simulated CO concentrations of
Run 33, when all grid cells are treated as land.
It is seen that the presence of a water surface within the region leads
/
to, in those cases of significant change, higher values of simulated pollu-
tant concentration. This is not surprising since the maps are exclusively
for daylight hours, and one would expect the presence of water to lead to
increased stability in the lower atmospheric layers during the day for the
time of year simulated.
Significant changes are both more widespread and more frequent than those
obtained by varying the input terrain slope. However, the maximum changes
produced by the elimination of water grid cells are smaller in magnitude than
those seen in Figures 82 through 86.
3.3.3.4 A simulation run with observed winds replacing model simulated
winds (Run 31)
The changes in simulated CO concentration obtained by using linearly in-
terpolated ground-level winds, instead of allowing the model to predict winds,
are shown in Figures 92 through 96.
Only the first two maps, for 08P and 10P/September, show a definite,
though minor, tendency toward more accurately simulated pollutant concentra-
tions. The later maps in the series (12P through 16P) show, in many places,
a decrease in simulation accuracy. This suggests that model-simulated winds
may be at least as suitable as winds obtained from conventional observing
networks for the prediction of pollutant transport and diffusion.
78
-------
Figure 87. 08P/29 September
Figure 88. 10P/29 September
•RESD
I3.S
9.0
(-0)
WEST
. *n
\ 4.0
\it-o)
\
\
3.3\
(-0)\
3.0 /
(40 >.
>
z.o
(0)
BURK
*17.S
9.7
(-0)
5.5
(+0)
LENX
*
4.8
(+0)
0
4.7
(+0) ^
~v
3^3
(-0)
70.8
(-0)
CAP . ,
# '•"
175 *<•«»
#16
VER
7.7
(+0)
LONB
#12
8.7
--J^k.
3.4
(-0)
77.2
<-°>ELM
#
I
8*8
WHTR
70.7
(+0)
*
8.2
(+0)
\
(-o)X
77.2
AZU#* (-0)
105
9."0
(+0)
77.0
(+0)
70.2
(+0)
7.2
^ (40)
•RESD
7
8.4
(-0)
WEST
*s
\-0)
\
A
3.2 \
(-0) \
1
3.0 /
(+o) \^
X
3.0
(-0)
BURK
#10
•
8.8
(-0)
4.8
(+0)
LENX
# .
U 3.2
(+1)
3*7
-^Cj
3.5
(-0)
*
8.8
(-0)
CAP
#8.3
•(-0)
VER
.
4.0
(+1)
LONB
#7.5
3.*8
3.5
(-1)
9.4
(+0) ELM
*
6,*9
(-2)
WHTR
. *«
5.5
(+0)
4.8
(-0)
\
(-H3) X
9.8
a
7.2
(-D
.
8.9
(-0)
•
4.2
v (40)
Figure 89. 12P/29 September
Figure 90. 14P/29 September
:RESD
15
70.7
(+0)
WEST
*l!
1 •
\4.7
\(-0)
\
A
2.fl\
(•»0)\
1
2*9 y"
(+0)\
\
3'.0
(+0)
BURK
*I05
10.0
(+0)
;.s
(•w)
LENX
* .
3.2
(+0)
|3.*7
(+0) y
^J
3."e
(-D
8*7
(-0)
CAP
* 5.9
" '^s0
VER
3.S
(-0)
LONB
*6
•
3.5
— v— <*£
3.V.
(-1)
•
8.7
(-0)
Ea
8.3
(-2)
WHTR
. **
4.6
(-0)
•
4.7
(-0)
r\
(0)\
a.o
AZU*' W)
6.S
8.5
(-D
5.3
(+0)
5.0
(-0)
4.0
(+0)
•RESD
3
3.S
(-0)
WEST
\-0)
3.0\
(-0) \
2.9 /
(-0) \
%
2*. 9
(-0)
BURK
#7
5*7
(-0)
LENX
s,-s
(-0)
3.7
(-0) s
"^-J
3.9
(-1)
a'.z
(-0)
CAP
I* '••{-«
#3
VER
3.9
(-1)
LONB
#5
3.*9
(-1)
5.'2
(-0)
ELM
#
3.5
(-0)
WHTR
. #1
4.4
(-1)
4.7
(-1)
\
(0) X
3.8
IS
3.3
(-1)
4.8
(-1)
4.5
(-1)
3. 'a
Figure 91. 16P/29 September
•RESD
S
3.'3
(-0)
WEST
\&T
(-0) \
2*9 y
(-0)\
^
2*9
(-0)
BU^RK
3/3
(-0)
5*8
(+1)
LENX
(:D
3.3
(-0),
"V
4.2
(-2)
3*2
(-1)
CAP , .
# 3.8
t ' (-2)
#15
VER
4! 4
(-1)
LONB
*s
•
4.2
(-2)
•
3.7
(1JM
f
3? 3
(-0)
WHTR
.*»
4.8
(-1)
4*4
(-1)
(:1)X,
3.0
AZU#*(-0)
3*2
(-0)
c-i)
4*2
(-1)
4*0
(-1)
# CO Monitoring Station
U Observed CO Concentrations at Station
*.8 CO Coneantratlona Predicted by Base Run 37
(0) Change In CO Concentration for Run 33
Figures 87 -• 91. The maps of Fig-
ures 37 through 41 are repeated.
Also shown are changes in simulated
CO concentrations induced in Run 33,
in which all horizontal grid cells
were designated as containing a
land-air interface.
79
-------
Figure 92. 08P/29 September
rRESO
13.5
9.0
WEST
. *ri
S. 4.0
\
\
3:3\
(-0) \
3.0 y
(•K» \
N
3*0
<-K»
BURK
*17.5
9.7
(+3)
8*2
(+2)
ILENX
*
4.*8
(+2)
•
4.7
(+2) ^
^V
3. *3
(•K»
30.6
(+2)
CAP. ,
* 7.3
17.5 ' (+D
*16
VER
i*.7
(+3)
LONB
#12
~v-\.
3*4
(•to)
37.2
(+DELM
#
I
a. e
(+D
WHTR
jt 15
70*. 7
(+2)
•
8.2
(+2)
\
77.2
AZU** (-0)
10.S
9*9
(+1)
77*. 0
(+1)
•
30.2
(+1)
7.2
(+3)
\.
Figure 93. 10P/29 September
'RE5D
I
8.4
(«)
WEST
. *<
N. 4.0
V
\
3! A
\
<-K>)^
3*0
(-0)
BURK
*10
8.6
(+2)
(+1)
LENX
•»
IS 3*2
(+1)
3.*7
(+1) ,
^J
3.*5
(-0)
8.8
(+3)
CAP 6 3
*12
VER
4.*0
(+2)
LONB
*7.5
•
3*5
(+0)
9.4
<+3> ELM
*
5.9
(-1-1)
WHTR
5.*5
(+4)
(+3)
\
(+1) \
9.6
AZU** (+2)
9
7.2
(+2)
6?9
(•«)
(+3)
4.2
J-t-3)
Figure 94. 12P/29 September
•ftESo
(.5
(io)
WEST
. **!
\(-0)
2r9\
\
2:9 y
(-0)\
X
3:0
(-0)
BURK
*10.5
il.o
(-0)
LENX
* .
3.2
(+0)
13*7
(HO)
^J
3.*5
(-1)
8.*7
(+2)
CAP
* 5.9
15 ' (-«)
#5
VER
3.*6
(-1-1)
LONB
*6
•
3.*6
(-0)
8.*7
(+3)
*
ELM
S.*3
(+2)
WHTR
4*6
(+3)
4.*7
(+1)
V
8.0
AZU** (+3)
6.5
6*5
(+3)
•
5.3
5.0
(+3)
4*0
Figure 95. 14P/29 September
•RESO
3
•
3.5
(+7)
WES1
V ' **
\4.2
\(-0)
3.0\
(+0) \
1
i.9 //
(-o)V
V
•
2.9
C-0)
BURK
#7
•
4.7
(-W)
5.*7
(-0)-
^ENX
" 3.*5
(-0)
•
3.7
(-0) j
r^J
•
3.9
(-1)
8*2
(+1)
CAP
* 4.7
IS *(+!)
*3
VER
•
3.9
(-0)
LON.B
*&$.
•
3.0
(-1)
5.2
(+2)
ELM
•jt
3.*5
(«)
WHTR
. *1
4.4
(+1)
•
4.7
(+1)
3.3S.
(-0) \
3.8
AZU* (+1)
IS
3.3
(+1)
•
4.5
(+2)
•
4.5
(+2)
3.8
(•••I)
Figure 96. 16P/29 September
rRESD
$
3*3
(+7)
WEST
\-;C
\
\
3.*7\
(-0)\
2.*sy
(-0) \
X
2.*9
(-0)
BURK
*5
3.*3
(+7)
5*6
(-1)
LENX
4 *0
(-1)
3.3
(-0) ;
^J
•
4.2
(-2)
3.2
(+5)
CAP T ,
Jt 0.0
{ • (-0)
*15
VER
4/4
(-1)
LONB
*s
3.9
•
4.2
(-1)
3.7
(+2)ELM
*
2
3*3
(+1)
WHTR
*1
(-0)
4.4
3.X.
.3.0
AZU** (+D
S
3?2
3T9
4.2
(+1)
4*0
* CO Monitoring Station
9.S Obierved CO Concentrations at Station
9.8 CO Concentrations Predicted by Base Run 37
(0) Change in CO Concentration for Run 31
Figures 92 - 96. The maps of Fig-
ures 37 through 41 are repeated.
Also shown are changes in simulated
CO concentrations' induced in Run 31,
in which observed winds were used
throughout the simulation, rather
than model-simulated winds.
80
-------
3.3.3.5 A simulation run with increased soil moisture parameter
(Run 19)
The changes in simulated CO concentration obtained by using an increased
surface moisture index over the land areas are shown in Figures 97 through
101.
The most widespread and largest changes appear in the earliest maps.
These changes decrease the simulation accuracy significantly.
3.3.4 Comparison with the Results of Other Models
A previously described model (Roth et al., 1971), applied to the same
observational period, differs from the model used here in three major re-
spects. First, the three-dimensional winds and diffusivities are input to
Roth's model, and are obtained from analysis of data observed throughout the
simulation period. Second, a finer and more extensive horizontal grid was
used by Roth et al. (see Figure 2). Third, the source distributions for CO
used by Roth, besides being defined on the finer grid (and presumably cap-
turing more of the detail in vehicular traffic patterns), also include the
regional airports in the source patterns. While these make a minor contri-
bution to the region's total CO inventory, they may be locally and sporadic-
ally important in producing major variations.
We omitted the airport sources in the simulations described here on
the assumption that their contribution would be minor as long as the coarse
horizontal grid was used. This was borne out by a test simulation (not de-
scribed elsewhere in this report) in which the airport sources were included
with their contributions spread over the individual 64 square-mile grid cell
in which each is located. The simulated concentrations were essentially the
same as those shown for Run 37.
81
-------
Figure 97. 08P/29 September
Figure 98. 10P/29 September
•RE5D
US
9.0
(-0)
WEST
\ 4.0
\(+D
\
A
3.3\
(+0)\
1
3*0 /
(+0)\
X
3.0
(+0)
BURK
*17.S
9*7
(-1)
6.2
(+1)
LENX
* .
4.0
(-0)
4.1
(-0)
3*3
(+0)
10.6
(-3)
CAP
# 7.3
17.S V+l)
#16
VER
.
7.7
(-2)
LONB
*12
-Y^i,
3.4
(•to)
11.2
ELM
*
s
8.6
(-1)
WHTR
10.1
(-4)
8.2
(-3)
\
(-D X|
11.2
AZU* (-*>
105
9.9
(-3)
* •
11.0
(-5)
10.2
(-5)
7.*2
,(-3)
rRESD
8.4
(-2)
WEST
. *»
3.2\
(+0) \
3.0 y
(+0) \
3.0
(-0)
BURK
*IO
8.6
(-2)
4.0
(+1)
LENX
«»
(-to)
3*1
(+0) /
"^J
3*5
(-0)
0.0
(-2)
C*AP ».*
•(•to)
#12
VER
4*0
(-0)
LONB
*IS
3,*0
3.5
(-0)
9.4
<-2)ELM
*
6.9
(-1)
WHTR
. *U
5.5
(-1)
4.'0
(-0)
\
(-0) X
9.0
AZU*'("3)
9
7*2
(-2)
•
6.9
(-2)
6.2
(-2)
4.*2
v (-1>
Figure 99. 12P/29 September
•ftESB
I.S
10.1
(+3)
WEST
\%)L
\
A
2.9\
(+0) \
1
2*.a y
(+0)\
\
3*0
(•W)
BURK
. *lfl.S
10.0
(-3)
4.0
(+0)
LENX
* .
3.2
(•to)
3.1
(+0) y
^J
3*0
(-0)
8.*7
(-2)
CAP
# 5.9
is • (-to:
*s
- VER
.
3.6
(40)
LONB
*6
3.5
3*0
(-0)
0.*1
(-1)
*
. ELM
6.3
(+0)
WHTR
• *2
4.6
(-0)
4.1
(+0)
3.2\.
(-o) N
0.0
AZU# (-2)
6.S
6.5
(-0)
.
5.3
(-1)
5.0
(-0)
4.0
^ (-0)
Figure
100. 14P/29
September
•RESO
3.5
WES1
v. '**
\4.2
Vo)
s!o\
2.9 x/
2.9
(•to)
BURK
*1
4.1
(+0)
•
5.1
(•to)-
LENX
*5 3*5
(•to)
3.1
S^
3.9
(-0)
8.2
(-2)
CAP
#4.1
" V?1
VER
•
3.9
(-0)
LONB
3.6
•
3.9
(-0)
5.2
(-2) ELM
*
.
3.5
(«)
WHTR
. *1
4.4
(+0)
4.1
3.6
AZU-jf(-O)
IS
•
3.3
4.*0
(-0)
4.5
(•to)
3.0
„ (-0)
Figure 101. 16P/29
September
»RESO
3'.3
0-5)
WEST
N"*45
A
\
(•K>) ^
V
2*9
(•to)
BURK
3*3
(+5)
5.6
(+0)
LENX
* .
4.0
(+0)
3.3
(+0)
^J
4.2
(-0)
3.2
(+2)
CAP
* 3.8
*i+V
• VER
4*4
(+0)
fsa
4.2
(-0)
3.1
<-°>ELM
f
•
3.3
(+3)
WHTR
.*»
4.0
(+0)
4.4
X
(•to) X
3,0
AZU$ (+0)
3.2
(+3)
3.9
(•to)
(+i)
4.0
v (+0)
# CO Monitoring Station
U Observed CO Concentrations at Station
9.0 CO Concentrations Predicted by Base Run 37
(0) Change In CO Concentration for Run 19
Figures 97 - 101. The maps of Fig-
ures 37 through 41 are repeated.
Also shown are changes in simulated
CO concentrations .induced in Run 19,
in which the surface moisture parame-
ter was increased to 0.2 in all land
grid cells of the model.
82
-------
Figures 102 through 111 repeat the time series originally shown in
Figures 42 through 51, but also show time series simulated by the model of
Roth et al. (1971). In general, the simulations of Roth et al. are more
accurate, particularly for the monitoring stations WEST, CAP, VER, and WHTR.
Our model is clearly more accurate only at AZU, particularly when the doubled
terrain slopes are used (also see Figure 69). The simulated time series are
remarkably similar at RESD, LENX, and LOME. This could be optimistically
interpreted in terms of the general consistency of the modelling concepts
used in these models, considering the large number of differences in the
detailed modelling assumptions.
Another previously described model (Sklarew, et al., 1971) has been
applied to a different observational period. It also differs from the model
used here in three major respects. First, much less complex three-dimensional
wind and diffusivity fields are input to the model. Second, a finer and more
extensive grid, on the same scale used by Roth, et al., is used to define
the automotive source fields, and for the prediction. Third, a highly sophis-
ticated computational technique (PIC) is used for the calculation of pollu-
tant transport in order to eliminate numerical (calculational) diffusion.
(The interested reader will find an extended discussion of this suprious
diffusion as it appears in our models in Section 2.5 of our previous report
[Pandolfo, et al., 1971]). We have prepared, under the present contract,
a version of the PIC formulas for inclusion with our models. We believe
that other modifications should, at this stage, take precedence in attempts
to increase the realism of the simulation.
Our basic reason for this belief is that the simulation outlined by
Sklarew, et al. (1971), and shown in Figure 31 of their report, is not
obviously more accurate than that shown here.
83
-------
1
•era wni up/a upm «n>» tutm
TlKM
f •
!
8 i
Mmt mm
8 ,
Figures 102 - 106. The time series
of Figures 42 through 46 are re-
peated. Also shown are simulated
time series obtained by Roth et al
(1971).
Kfl-a own ma
Tkm
-------
1
00
tun wm> tmn arm arm asm
tan
Tkw
Figures 107 - 111. The time
series of Figures 47 through
51 are repeated. Also shown
are simulated time series
obtained by Roth et al.
(1971).
ami
Itow
nrr»
-------
4.0 SUMMARY AND CONCLUSIONS
The urban boundary-layer model described in a previous report (Pandolfo,
et al., 1971) was modified and used in test runs to estimate the degree to
which observed spatial and temporal variations of meteorological conditions
and concentrations of a stable air pollutant within a complex urban region
can be realistically simulated.
The major modification relevant to this project was the inclusion of an
average terrain slope vector for each grid square in the input to the model.
This slope is used to calculate topographically induced vertical velocity
and its effects on the temperature, humidity, and pollutant within the region.
The modelling equations, as modified, are presented. A series of simulations
was conducted on a data base obtained from the results of observation pro-
jects at Los Angeles, California, on September 29-30, 1969. Forty full test
simulations were run over the span of a few months. The observed meteoro-
logical parameters form a sparse set of data in the four-dimentional (space-
time) mesh needed to define the relevant meteorology for this application.
Many of the runs varied the meteorological input about the standard (observed)
set. It.has, therefore, been demonstrated that an economical, objective,
physically consistent, and precisely specified (though with some artibrary
elements) procedure has been achieved for obtaining and predicting the three-
dimensional meteorological fields needed. This procedure provides a neces-
sary supplement to the pollutant transport-diffusion portion of the model
developed here; or to other pollutant models developed elsewhere. Hereto-
fore, these other models (e.g., Sklarew, et al., 1971; Eschenroeder and
Martinez, 1971; and Roth, et al., 1971) have relied on incomplete meteoro-
logical fields subjectively and laboriously obtained, that are difficult to
generalize to other weather situations. In several of the runs the input
87
-------
topography, land-water distribution, and other physical characteristics of
the underlying surface was varied. The results demonstrate that ready gen-
eralization to other regions can be expected.
For purposes of economy in these initial tests, and for reasons of
physical consistency, a major simplification was made, and the modelled
region is simulated on a relatively coarse 5x5 horizontal grid with 8-mile
grid sapcing. This is in contrast to other models which focus on the pre-
diction of pollutants only, and which are based on finer horizontal grids
with 2-mile grid spacing (Sklarew, et al., 1971; Roth, et al., 1971). None-
theless, the temporal and spatial variations of those meteorological condi-
tions for which observations are available, via., air temperature, humidity,
and wind, are simulated with an encouraging degree of realism. Temporal
and spatial variations of CO are also simulated fairly realistically, with
somewhat less accuracy than in the model described by Roth, et al. (1971),
and with accuracy equivalent to that shown by Sklarew, et al. (1971). It
is reasonable to expect improved simulation accuracy with the finer hori-
zontal resolution used in these other models, and the performance of such
simulations with this model is strongly urged.
Slight differences of simulation accuracy of surface-layer CO concen-
trations are seen when the model is modified to neglect the presence of the
ocean, to neglect the terrain slope, and to increase the soil moisture parame-
ter. On balance, these changes appear to decrease the simulation accuracy.
More significant increases of simulation accuracy at two of the sta-
tions in the vicinity of the San Gabriel Mountains in the early afternoon,
are obtained when the smoothed terrain slopes are doubled in magnitude.
This reinforces the expectation of improved simulation accuracy with a
finer grid, and, consequently, less smoothing of the topographic features.
88
-------
Noticeable changes in simulation accuracy are seen when observed winds
are used to circumbent the model prediction of wind. The increases in accur-
acy at some stations and times appear to be balanced by decreases in accuracy
at other stations and times.
Insignificant changes in simulation accuracy are seen when the minimum
free atmospheric value of diffusivity is increased by an order of magnitude.
A major deficiency of the simulated CO variations is the underestimation
of morning peak concentrations in the earliest stages of the simulation at
many horizontal locations. A similar deficiency was noted in the model
simulations of Sklarew, et al. (1971) and was eliminated by modification of
the modelled vertical distribution of vertical eddy diffusivity in that study.
The model of Roth, et al. (1971) does not exhibit this deficiency, but does
overestimate morning peak concentrations at other stations.
Before other measures to eliminate this deficiency are investigated, a
test simulation should be performed with the refined grid (2x2 mile) as used
in other models, and with the addition of airports to the source inventory.
This is expected to yield the greatest improvement in accuracy (through
the more detailed definition of the emission sources and topography) at no
development cost, with only increased, and easily bearable, computer costs
added.
Other recommended changes, in order of decreasing priority, are:
1) inclusion of additional reactive pollutants along with their
associated photochemical formulations (e.g., as outlined in
Roth, et al., 1971);
2) generalization to a non-hydrostatic formulation (see Appendix D);
3) the inclusion of a (humerically) more accurate differencing
scheme for the horizontal transport with explicit treatment
of the horizontal diffusion (see Appendix E); and
89
-------
4) optimization of the computer program to eliminate details that
are superficially unnecessary in this application, e.g., the
requirement to carry as many grid levels in the soil as are
needed in the water portions of the model.
Implementation of item (1), according to rough estimates based on Roth,
et al. (1971) and Sklarew, et al. (1971), would increase the running time by
an order of magnitude. However, such time, and the cost implied, is feasible
for the full range of applications in air quality simulation on more advanced
computers for which the model is intended than the IBM 360/65 used in this
project. Relative to this increase, the increases (or savings) in computa-
tional costs resulting from the implementations of items (2), (3), and (4)
become trivial. The development of the necessary formulations for (2) and
(3) have been carried out (Appendices D and E). In this context, it is our
opinion that the kinds of programming changes that might be carried out under
(4), though superficially appealing on aesthetic grounds, and easily imple-
mented, should take on the lowest priority.
In summary, we present here what is to our knowledge the only documented
and validated three-dimensional model. We know of no other model that has
achieved multi-day predictions or simulations of both meteorological vari-
ables and conservative pollutants on the spatial scales considered here. As
is inherent in any complex model, there are apparent variations in the degree
of accuracy of treatment of the many processes and parameters included in
both the input and formulation. It is, however, unreasonable to expect agree-
ment among scientists as to exactly which components are treated too crudely,
and which are treated too accurately, until much more experience has been
gained with the model. It has also been demonstrated that many useful re-
sults can be provided in this application while this experience is being
gained.
90
-------
5.0 REFERENCES
Beeton, A.M., 1962: Variations in radiant energy and related ocean temper-
atures. Proc. Fifth Conf. on Great Lakes Res., 1962, pp. 68-78, Univ.
of Michigan, Great Lakes Res. Div., Ann Arbor, Michigan.
Cox, R.A., M.J. McCartney, and F. Culkin, 1970: The specific gravity/
salinity/temperature relationship in natural sea water. Deep-Sea Res.,
17, 679-689.
Eschenroeder, A.Q. and J.R. Martinez, 1971: Concepts of 'Photochemical Smog
Models. Tech. Memo. 1516, General Research Corp., Santa Barbara,
California, 123 pp..
Gerrity, J.P., Jr., 1967: A physical numerical model of the prediction of
synoptic-scale low cloudiness. Mo. Wea. Rev., 95, 261-282.
Goody, R.M., 1964: Atmospheric Radiation, Oxford at the Clarendon Press,
London, 436 pp.
Halstead, M.H., R. Richman, W. Covey, and J. Merryman, 1957: A preliminary
report on the design of a computer for micrometeorology. J. Meteor.,
14, 308-325.
Hess, S.L., 1959: Introduction to Theoretical Meteorology, H. Hold & Co.,
New York, 362 pp.
Kitaigorodsky, S.A., 1961: On the possibility of theoretical calculation
of vertical temperature profile in upper layer of the sea. Bull.
(Izvestia) Acad. Sci., USSE, Geophys., Ser. 3, 313-314.
Kondrat'yev» K., 1969: Radiation in the atmosphere. Int'l. Geophys. Ser.
12, Academic Press, New York.
Pandolfo, J.P., 1966: Wind and temperature profiles for constant-flux
boundary layers in lapse conditions with a variable eddy conductivity
to eddy viscosity ratio. J. Atmos. Sci., 23, 495-502.
, 1969: Motions with inertial and diurnal period in a numerical model
of the navifacial boundary layer. J. Mar. Res., 27, 301-317.
, and P.S. Brown, Jr., 1970: Simulation Experiments with a Numerical
Sea-Air Planetary Boundary Layer Model and Its Extension to Three
Space Dimensions. Vol. I, TRC Rpt. No. 7499-386a, The Travelers Re-
search Corporation, Hartford, Connecticut, 95 pp.
, M.A. Atwater, and G.E. Anderson, 1971: Prediction by Numerical
Models of Transport and Diffusion in an Urban Boundary Layer, Final
Rpt., Vol. I, The Center for the Environment and Man, Inc., Hartford,
Connecticut, 139 pp.
Pielke, R., 1972: Comparison of a Hydrostatic and an Anelastic Dry Shallow
Primitive Equation Model. Tech. Memo. ERL OD-13, U.S. Dept. of Comm.,
NOAA, ERL, Boulder, Colorado, 47 pp.
91
-------
Pierson, W.J., 1964: The interpretation of wave spectra in terms of the
wind profile instead of the wind measured at a constant height. J.
Geophys. Res., 69, 5191-5204.
Roth, P.M., S.D. Reynolds; P.J.W. Roberts, J.H. Seinfeld, 1971: Develop-
ment of a Simulation Model for Estimating Ground Level Concentrations
of Photochemical Pollutants. Final Rpt., 71 SAI-21, Systems Applica-
tions, Inc., Beverly Hills, California, 51 pp., Appendices A-F.
Sklarew, R.C., A.J. Fabrick, and J.E. Prager, 1971: A Particle-in-Cell
Method for Numerical Solution of the Atmopsheric Diffusion Equation,
and Applications to Air Pollution Problems. Final Rpt., 3SR-844,
Systems, Science, and Software, LaJolla, California, 163 pp.
Sverdrup, H.U., M.W. Johnson, and R.H. Fleming, 1942: The Oceans, Prentice-
Hall, Inc., New York, 1087 pp.
92
-------
6.0 ACKNOWLEDGEMENTS
The authors wish to acknowledge the thorough review of the equations
and formulations by Mr. Philip S. Brown, Jr. The collection and preparation
of the input data owes much to the guidance and labor of Mr. Gerald Anderson.
The computer programs were prepared and run by Mr. Joseph Sekorski, with
consultation and guidance by Dr. Marshall A. Atwater. Ms. Margaret G. Atticks
carried out the formidable task of typing and figure preparation entailed in
this report.
The University of Connecticut Computer Center contributed significantly
more, in computer time, than the computer project budget.
93
-------
Appendix A
Specifications for the Computer Program
A.I Introduction
The purpose of this section is to describe the set of finite-difference
equations which are solved in the computer program. These equations are the
analogues of simplified forms of the prognostic differential equations (see
Section 2.0) for various scalar properties in the planetary boundary layer
of the atmosphere, ocean, and soil system. The scalar properties include
wind and current components, temperature, specific humidity, salinity, and
pollutants. The units of specific humidity and salinity in this section
are g/kg rather than g/g as in Section 2.0. Much of these specifications
is based on the specifications and modifications in the Appendices to
Pandolfo et al. (1971).
The upper and lower boundary values required for the numerical solu-
tions of the equations are to be input or computed from input. The equa-
tions are to be solved on a specified, arbitrarily spaced, vertical mesh,
which is also input. Figure A is a schematic diagram of a vertical grid
mesh. This figure should provide the reader with a better understanding
of the indexing necessary for the presentation of the finite-difference
equations. The two vertical grid points associated with the interface
represent the air and land (water) sides of the interface.
Volume II of this report contains the FORTRAN program listings for
all the subroutines and flow diagrams for the MAIN program of both versions
(RIGID LID and FREE SURFACE) of the model.
95
-------
SN-1
6N-2
-1+1
'1-1
'1-1
'1-2
'1-2
Az
1+1
> Air Layer
Az
l-2
Az
Interface zi_ = zj. «
» Land or Water Layer
Figure A. Schematic diagram of vertical grid mesh.
96
-------
A. 2 Input
The input is specified either for the entire grid, for one location in
the grid, or for each station of the grid. Input specifications for both
versions of the model are contained in Volume II.
A.3 List of Additional Symbols
A.3.1 General Parameters
I Number of vertical grid levels from the bottom of the bound-
ary layer to the air side of the interface.
I Number of vertical grid levels from the bottom of the bound-
ary layer to the water or soil side of the interface.
J Number of stations in the x-direction.
max
K Number of stations in the y-direction.
max
N Number of vertical grid levels within the boundary layer.
z Heights of input gradients (cm).
O
ZW Height of vertical grid level used in computation of wave
dimensions (cm). This height is 1950 cm or the lowest grid
level above 1950 cm.
A Horizontal grid length.
e Infrared emissivity of surface.
S
A.3.2 General Grid Location
J0»fe0 The horizontal grid indices corresponding to the x,y loca-
tions (respectively) of the input initial vertical profiles
of the dependent variables.
PX(Z ,0) Initial values of Pollutant 1 (g/100 g); zj^ < z s. ZN .
P2(z.i»0) Initial values of Pollutant 2 (g/100 g); Zj. < z <. ZN .
q(z ,0) Initial values of specific humidity (g/kg); ZL. < z <. SL. .
s(z. ,0) Initial values of salinity (g/kg); z, £. z <. z_ .
T(z ,0) Initial values for temperature (°K); z <. z <_ z .
97
-------
u(z ,0) Initial values of eastward velocity components (cm/sec),
zx < Zj < ZN .
v(z ,0) Initial values of northward velocity components (cm/sec),
z, < 'z . <. z
1 - j n
A. 3.3 Input at Each Station
h- . ,h_ . Fractional strength of emission for pollutants 1 and 2 at
1 * the i(th) height level.
z ) Soil thermal conductivity for each soil layer (cm2/sec) ,
J z < 0 .
n Number of geos trophic wind and current input values.
O
S1 (t) Source strength for pollutant 1 (yg/m2) .
S-Ct) Source strength for pollutant 2 (yg/m2) .
U (zj ,t) Eastward components of the geos trophic current at the inter-
. face (cm/sec).
V (zj ,t) Northward components of the geostrophic current at the inter-
* face (cm/sec) .
U (zvr>t) Eastward components of the geostrophic wind at the upper
° boundary (cm/sec) .
V (zM,t) Northward components of the geostrophic wind at the upper
° boundary (cm/sec) .
p Density of soil (g/cm3).
s
A. 3. A Computations of Program Constants
2
c - 3(R )/(1 + R )
c c
c - .239
a
cw - .935
C - 3k"/3 IT2/3
f = 14.585 sin* x 10~5
g = 980
h = k2(3/c)3/2
k - .4 •
98
-------
kl
nT
r
TT
P,
W
.14
3/2
1/2
- l/7a = .048
-3
10/3
10
.98 x 10"
3.1416
i-«f
= 1.354 x 10
-12
A.4 Parameters Constant Over Entire Grid
The height grid parameters are computed according to the notation
ZJ + 2J+1.
» J
, J
.....
(58)
..... N
where z. i z <. z.. , and the absolute value is used because depths that
belong to the interface are taken to be negative. The value of z_ will
always be zero since Zj =• Zj = 0 . The distance between vertical grid
points is given by
(59)
z. - z.^ , j - I+l N
The input horizontal gradients
99
-------
31
3s , , 3s
te (V ; 37
3u , . 3u , x
^(«j) ; ^(.j) Zl < Zj
3v > » 3v
zi-zj ~
zi-zj -ZN
are to be computed through linear interpolation in the vertical at levels
for which they are not specified. Linear interpolation is performed if
more than one gradient of a dependent variable is input within a layer.
If only one gradient changes at the interface, the input value is given
at Zj_ and Zj for both the subsurface and the air. These gradients
are independent of horizontal location and will be used at the horizontal
grid boundary, whenever necessary.
A.5 Initial Profiles
The initial profile of the variables u,'vtT,q,s,P., and
P. are computed from the input vertical profile at grid location SQ^-Q •
The initial profile (t » 0) of X at horizontal grid point (j,k) and
height n is computed from
100
-------
'»' <> • • • • (60)
where the gradients are the input gradients. For u , v , and I , the
added restriction
(Xi)I~ - (X )^" (61)
is used. At each land station, a subsurface temperature profile is input
for each station and the subsurface profiles of u , v , s , P, and P^
are set to zero.
A. 6 Eddy Exchange Coefficient
The eddy exchange coefficients are computed at each time step for each
station. The subscript m is used to indicate eddy viscosity coefficients.
Eddy conductivity-diffusivity is subscripted with e . Special formulas
for the exchange coefficients at the z_ . and z - levels are given in
Section A. 10.
If the subsurface is water, a relevant wave property used to determine
the mixing is the vertical gradient of the scalar value of the average orbi-
tal velocity in a turbulent wave. For a given non-zero \ and 6 , this
wave property, S , is given by
—Z4 *ZlT/A /£1\
e J , (62)
where z- <. z. <. z . If no characteristic wavelength is specified (viz.,
X • 0 for a water subsurface), 6 is set equal to .055 and S is com-
puted from the average wave height, H , based on the formulas,
101
-------
W(ZW, t)
w(zw_lt t)
W1950(t)
tu2(ZW, t) + v2(ZW, t)]1
[u2(ZW_r t)
W(ZW,
Ll» t)]1
(63a)
(63b)
(63c)
where ZW - equals the vertical grid level immediately below ZW , and
An 1950/ZW
ZW
(64a)
ZW
-1
ZW
1950
ZW
(64b)
ZW
-1
From the wind speed at 1950 cm, the average wave height can be computed
H(t)
1.54165
[W1950(t)] .
(65)
Here S ' is computed:
•*«;
0
H
[««2lT|
6 J • J
•
»
H -
0
0
(66)
where z.. <. z < ZN . To set a lower limit of the atmospheric exchange
coefficient (see Section A.6.1.4), and for printout use, the average height
of the 1/3 highest wave is computed by
H1/3(t)
.8727 x 1.5989 H(t) .
(67)
If the subsurface is land, S (z.) = 0 for z.. <. z. <_ z^. .
102
-------
Coefficients used in the eddy exchange formulas are computed from one
^
of the following sets of equations according to the value of X ;
«• for X ^ 0 ,
(68)
f *,)2
V± \ _ k2U 4. ££. . £ < 7 < 7
V in 2J * i - j - 1-1 '
for X = 0 ,
VA \
Zj,t)
A. 6.1 The Atmospheric Layer (z_. 0 <. z . <. sL)
" - J N
N
(69)
2
Note that the computed quantities S , L , f , Ri , K ,K are under-
stood to be dimensioned variables with the coordinates (z.,t):
(70)
- T(zA
(71)
Az.
h _ r
-------
If T(S + S ) - 0 , check L ; if L < 0 , go to A.6.1.3; if L = 0 , set
w
Ri = 0 ; and if L > 0 , set Ri - |l/a|. If Ri >. 0 , proceed as in
A. 6. 1.1; if -.048 < Ri < 0 , proceed as in A.6.1.2; if Ri <. -.048 .pro-
ceed as in A. 6.1.3.
A. 6. 1.1 If Ri L 0
K = K (74b)
e m
and proceed to A.6.1.4.
A.6.1.2 If -.048 < Ri < 0
(75a)
K = - - (75b)
(1 - aRi)
and proceed to A.6.1.4.
A.6.1.3 If -.048 >. Ri
<76a>
-1/2 ,A
' l"l •• *. lQk
m,e j
rA ^K<107 ;
104
-------
where A = A2 - zQ2 for a land subsurface, or A = ^i/^2 ^or a water
surface.
A. 6. 2 The Oceanic Layer (z <. z < z .)
Note that the computed quantities S , L , Ri , s , f , 3a_/3T ,
3a_/3s , K , and K are understood to have the coordinates (z ,t).
i e m J
s(z ) + s(z )
8 = - m - L (78)
T(z ) + T(z
f ^ - - L. . • (79)
Compute
3 a
•—- - .05819402 - 2(.00811465413) f
oi
+ 3(.476600414 x 10~**) f2 + 2(.389187483 x lO"1*)!? (80)
•
- .25310441 x 10"2 s + .28797153 x 10~5 s2 ;
= .79701864 - .25310441 x 10 T
+ 2(.131710842 x 10~3) s + .389187483 x lo"1* f2 (81)
+ 2(.287971530 x 10"5) si - 3(.611831499 x 10"7) s2 ;
L - -IO-^-H—ei- M + --L i—Jz±- H> ; (82)
AZJ J ds
and
-------
g«L
(83)
If (S + S )2 = 0 , check L ; if L < 0 , go to A.6.2.3; if L - 0 , set
16
Ri = 0 ; and if L > 0 , set Ri = Ri (Ri = 10). If Ri >_ 0 ,
' max max '
proceed as in A. 6. 2.1; if -.048 < Ri < 0 , proceed as in A. 6. 2. 2; if
Ri < -.048 , proceed as in A.6.2.3.
A. 6. 2.1 If Ri LO
(84a)
(84b)
and proceed to Section A. 7.
A.6.2.2 If -.048 < Ri < 0
- aRi]
"2
(85a)
and proceed to Section A.7.
A<6.2.3 If RI <.-.Q48
- o
-1/6
(S
; (s + sw)2 = o .
106
-------
A.6.3 The Soil Layer
When the subsurface is soil, a soil conductivity coefficient is pre-
scribed, K^ , and used in making soil temperature predictions.
A.7 Infrared and Solar Radiation Computations
Program specifications for the computation of infrared and solar radi-
ation will be given in a forthcoming report to be written by Dr. Marshall A.
Atwater at the Center for the Environment and Man, Inc. under IFYGL sponsor-
ship.
j
A. 8 Solar Absorption in the Ocean
Solar radiation absorbed by the ocean, IQ , is dependent on the inci-
dent solar radiation at the interface, iCzj.) , and the surface albedo.
Therefore, the amount of energy available for heating the ocean is given
t>y
IQ = I(ZI+)(1 - a) (87)
where a is the albedo computed from
a .0139 + .0467 tanZ (88)
with a minimum value of .03. This equation is based on results of Sverdrup,
et al. (1942).
The transmission at any level z is given by
I - In exp(- K\z\ secZ) , (89)
where < is obtained by linear interpolation using Table 1 (Section 2.2.4.2),
107
-------
The radiation absorbed in each water grid layer is calculated from the
difference in radiative flux between the top and the bottom of the water
layer. In the uppermost layer,
' xo - xi-i (90)
and the layer B+l (the bottom layer)
In the remaining layers , AI is computed similarly from
The heating rate at level i (i - B+l, . .., I_) is computed from
fdT] Pw
"i. - s,^' (93)
where heating is applied to the midpoint of the layer.
A. 9 Geostrophic Wind (Current) at Each Level at Each Time Step
In the following two sections (A. 9.1 and A.9.2), U (z -,t) » v
U (ZT. ,t) , and V (zj ,t) are computed by linear interpolation from the
o o
values of Ug(zN,t1) , V («Ntt±) , U (a^.t^) , and Vg(zI_,tjL) at input
times t. for the RIGID LID version of the model. The values of U (zi_,t
and V (zi_»t) are computed from the slope of the water surface in the
O
FREE SURFACE version for a land subsurface, U (ZT tt) and V (ZT ,t)
g . L- g A~
are zero.
108
-------
A.9.1 The Atmospheric Layer (j = 1+1, ... N-l)
and
U (z ,t)*T(z ,t) g*T(z ,t)
"g^j*1"' T(Z ) f
2T AZ 1 ST 1
.... fT i i - ,, f-r ~\\
•7, kz.; i 2/ v ^v —ll
V (z ,t)-T(z ,t) g-T(z ,t)
wrTt-^cB J -J
viz. , t; = - "" -•
8 ^ T^2N^ f
3T , ^ . Azi-l 3T , J
« (9 1 4- — — ^— — 1-7 11
N , r Az.
r 1 i
^ 2 T2/ ^
i=j+l VT ^z^^;
N , f AzJ
r 1 i
A 2 2/
where
3T
3T
(94a)
(94b)
(95a)
3T
3T
3x
(95b)
Since vertical grid points are most often unequally spaced and the vertical
temperature gradient is approximated by the following weighted center differ-
ence ,
Az,
(95c)
II
9z
^fedr^r'^l + (;
M+i
Az,
The horizontal gradients of the temperature relative to terrain are space cen-
tered gradients and are calculated from formulas presented in Section A.13.2.
A.9.2 The Water Subsurface (j = 1, ... 1-1)
10
r
, (96a)
109
-------
3T
«Jt.) . (96b>
The horizontal gradients of water temperature and salinity indicated in
these formulas are vertically averaged; e.g.,
3T . /-3T 3T
W /* N 1 W , v . W, N /m\
(97)
Using the results of eqs. (96a) and (96b), compute
.t) - v(Zl.,o-f
J 6 ^ KW
1-1 Az. '
I — |fi-(z.,t) , (98a)
j j P 3y i
i-j Kw '
(98b)
The values of U (ZT ) and V (ZT ) in the FREE SURFACE version of the
g i- g i_
model are calculated from
f S
For a soil subsurface, U and V are, of course, zero,
g 8
A.10 Interface Values
The boundary and interface values are computed for each station at
each time step.
A. 10.1 Humidity . ' " '
The saturated humidity at the interface is computed from
(100)
where
110
-------
7'5
EI+,t) - 35.66
The interface humidity is computed from
q(zI+,t) - M«qg(zI+,t) + (1 - M) q(zI+1,t) (102)
.where M is input for a soil subsurface and equals 1 for a water subsurface.
A. 10. 2 Winds, Currents, and Salinity
Whenever & j 0 , always a water subsurface, the computations proceed
for the interface value for u , v , and the salinity as follows :
AZT+1 AZT-1
U - - TFT? - \ - TT-7? - T^ - (103a)
PAKm('l+l) + pwKm('l-l)
AZI+1 AZI-1
AZT4-1 AZT-1
v(Zl ) - - -SJ» - - - _—- - p - (io3b)
pAKm(8I+l) , pwKm(zI-l)
Azl+l A2I-1
- 10
S(ZT )(E - PJ
— =• (105)
1 - 3(2 * 10
111
-------
The values of K and ,K at the grid level adjacent to the interface
are given by
(107b>
where Kfl is computed from eq. (68) if X ^ 0 , or from eq. (69) if X « 0
(these coefficients are a function of time for X = 0 ). Then proceed with
Section A. 10.3.
When 6 » 0 and the subsurface is water, the interface values of u ,
v , and s are computed as follows. Compute the value
B
(108a)
Then compute
B2 • in
(108b)
(109b)
If X < 0 , compute
S(z.,)
u*
and
u*
(110)
(111)
and proceed with computations from the computation of u(z_) below.
112
-------
If X > 0 , compute
* ""* 1
(112)
and
S(zx_)
(113)
The computation proceeds with the calculation of u(z ) and other variables
as follows:
(lUb)
-1
- q(2lf)] B/ ; K
u*
(115)
- 10
(116)
F(ZIJ
u*
(117)
s* -
(118)
and
S*'B
(119)
When the subsurface is land, 6=0 and X < 0 (i.e., X = -z«)-
113
-------
The computations follow the previous section for the water subsurface,
except that
S(2l_) - S(zI+) = 0 . (120)
All values of u and v at heights at and below the interface are zero.
This leads to the result
' • <121>
Similarly,
s(Zl_) « 0 . (122)
The remaining equations are unchanged.
A.10.3 Temperature
The interface temperature T is computed by using the heat balance
O • •
equation
I*
(123)
t+At
where R is the atmospheric radiation incident on the interface. The term
RU is defined as the net solar energy absorbed at the interface minus the
amount transmitted (if any) to lower layers (viz., water subsurface).
The incoming radiation, both solar and infrared, is computed by methods
which will be described by Dr. Atwater in a forthcoming report under IFYGL
sponsorship.
The heat of vaporization, L , is calculated from
L - 597.3 - 0.57 (T - 273.16) - 753.0 - 0.57 T C124)
8 8
114
-------
for T > 273.16 . If T < 273.16 , then L = 677.0 .
8 ~ 8
The value E is computed in Section A. 10. 2. The air heat flux is
calculated from
j - T(z
A ' HA ' -' +r • (125)
and the surface flux is calculated from
S - H± - Pj^VW—^ I ' (126)
where the i subscripts in eq. (126) are taken as w or s according to
the subsurface.
Compute the heat loss or gain to the interface from rain, P_ ,
PT ' "wcw Pr(Tr - V • (127>
where P is the precipitation ratio, and T is the rainfall, temperature
assumed .to equal the average wet bulb temperature near the sea surface.
All the parameters in eq. (123) have been previously defined in terms
of (T»)t .Since eq. (123) is nonlinear in T , an iteration procedure
is required at each time step. In the iterative procedure, R-. , R- , R ,
n a x
and L will remain constant. The values of A , S , P_ ., and E will be
recomputed using the value of T from the previous iteration. Calculation
O
of the right-hand side of eq. (123) using these quantities yields a provi-
sional value for [T (n)]f..At. for each iteration, n = 1, 2
For the initial iteration, the balance temperature used for the right-
hand side of eq. (123) is obtained by
= (Vt ' (128)
115
-------
In the special case when the initial value of T (1) leads to a nega-
O
tive value for the expression in brackets in eq. (123), a new approximation,
T (2) , is defined by
O
0.001 , an iterative procedure is used that assumes the
temperature at the preceding time step or iteration to be near the final
value. Thus, when the temperature used on the right of eq. (123) is too
low, the computed temperature will be too high and vice versa.
This procedure compares the signs of A and A , . Until the first
change in sign, the temperature used on the right-hand side in eq. (123)
for the following iteration is computed from
T (n+2) - T (n+1) + A . (132a)
g g n
116
-------
After a change in sign has occurred, it is assumed that the new balance
temperature is within the range T (n) to T (n+1) . Signs of A and
B o
A , are compared. When the signs are the same, the next balance tempera-
ture approximation is computed from
A
Tg(n+2) = T (n+1) + y- . (132b)
When the signs are opposite, the balance temperature approximation is computed
from
A
T (n+2) = T (n) + -^ . (132c)
g g 2
The iteration is then incremented and eq. (123) is then solved using the
latest approximation. A maximum of 25 iterations is used, although generally
less than 10 are required.
In the iterations, the following equations may be used:
/T (n+1) - T(z_ ,h
n)[ T8(n) - T(.^") J ' <133a)
1(z ) - T (n+l
- T n+
) -£(.) J -
q(z) * q(n+1)
where q (n+1) is calculated from the saturated specific humidity for
s
TB(n+l) .
O
At each time step after all iterations have been completed, a time
sum of each of the components of the heat balance is computed from
117
-------
IX At , (134)
where X is Rp , R , A , S , L*E , P , R , and aT1* .
tl a XX g
A. 10. 4 Air Pollutants
If pollutants are present, the value of the pollutant concentration is
computed as follows.
Compute the source emission S.. (t) as a function of time by linear
interpolation from the input values of S- (t) . Then the interface value
for the first pollutant is computed from
S (t) Az h
P1(zLf,t) = -^ i±L_iL + P1(zI+1,t) . (135a)
paKe(2I+l)
Similarly for the second pollutant, if present, compute
S9(t) Az . h,
P2(Zl+,t) - — L * + P2(zl+l't) * (135b)
paKe(2I+l)
A.10.5 Computation of Interface Fluxes for Printout
If 6^0, the computation proceeds as in Section A. 10.5.1. If
6 - 0 , the computation proceeds as In Section A. 10.5.2.
A.io.5.1 if-ay o
(136>
T(z ) -
, ) - T(z .)
(138)
118
-------
where S(z ) is computed from eq. (109a) and an analogous equation with
E and F computed in Section A.10.2 for 6 J 0 .
A.10.5.2 If 6 • 0
> - pa-u*2 , ' (139)
HA(zI+1,t)
(141)
where the i subscripts in eq. (141) are taken as w or s according to
the subsurface. E and F are computed in Section A. 10.2 for 6 = 0, and
u* is computed from eq. (110) or (112) for X < 0 and X > 0 , respec-
tively.
A. 11 Upper and Lower Boundary Values
A. 11.1 Winds and Currents
The winds and currents at the upper and lower levels are computed from
the geostrophic winds calculated in Section A. 9 or are prescribed through
linear interpolation of the input values of U (zN,t ) and V (z ,t.) .
The upper level winds are given by
U(zN,t) - Vg1:)
V(zN,t) = Ug(zN»t> » <142b)
where linear interpolation with respect to time is used to obtain values of
119
-------
U (Zjj.t) and V (z^t) at times for which they are not specified.
In the RIGID LID version of the model, at stations with a water sub-
surface, the currents are computed at the lowest grid level, z, , from
u(zlft+l) - u(Zl,t) + At.f[v(zltt) - V (zrt)] (143a)
v(zltt+l) = v(zrt) + Aff[U (Z;L,t) - u(zrt)] (143b)
where V (z..,t) and U (z, ,t) are obtained from the equations in Section
A. 9. The boundary conditions impose restrictions on the solutions of eqs.
(143a) and (143b); viz., for a coastal corner station,
u(Zl,t) - -0 (144)
v(Z;L,t) - 0 ; (145)
for a coastal station with a coast oriented in an east-west direction, eq.
(144) , and for a coastal station with a coast oriented in a north-south
direction, eq. (145).
The FREE SURFACE version of the model sets the current velocity compon-
ents at water stations equal to zero at and below the bottom boundaries, i.e.,
where J « 1 , . . . B .
A. 11.2 Temperature, Humidity, Salinity, and Pollutants
The upper and lower boundary values (at z = z., , ZN ) are held con-
stant for T , q , s , P. , and P. .
120
-------
A.11.3 Printout of Boundary Fluxes
T(zN,t) - pa.Ko(zN,t).S(£N,t) (147)
T(Z ,t) = p -K (z ,t)«S(z.,t) (148)
x w . m l i
,.) - T(zM
T(z ) - T(z.
(150)
E(zN,t) = -Pa-Ke(zN.t) (151)
(152)
where the i subscripts in eqs. (150) and (152) are taken as w or s
according to the subsurface.
A.12 Solution of the Finite-Difference Equation
The following sets of equations are solved simultaneously, one station
-at a time, for each station in the grid through an implicit three time level
differencing scheme which makes use of Gaussian elimination techniques for
tri-diagonal and block tri-diagonal matrices. Therefore, there are no sub-
scripts indicating horizontal grid locations; the subscript j denotes
grid location in the vertical and superscript t designates the time level,
viz., t -. 1 » past , t » present , t + 1 - future . The matrix of coef-
ficients contains two less rows than vertical grid levels in a layer (air,
land, or water); this is because the values of the dependent variables are
121
-------
TABLE A-l
Definitions of Generalized Coefficients
Quantity
Ke(Zj)
+ v^l -f«
2r
TA(2H-1>
l
Ke(2J>
s(z2)
2sihi*
1+ - J) £.6 A - (J - I_) £.6
6 SPl(z ) - 0 t > 6
P±(z2)
t Heating or cooling due to solar radiation.
+ Heating or cooling due to infrared radiation.
-------
known at the interface and the upper and lower boundaries of the boundary
layer. The relationship between the X 's and the X 's is clarified in
Table A-l.
A.12.1 T , q , s , P.. , and P _
The finite difference equations for T , q , s , P. , and P_ take
the form
d
1
+ C1X2
d2 = a2*l + b2X2 + C2X3
(153)
d * a
n
where the x^ 's denote the unknown values of T f q , s , P. , or P- ,
and the coefficients are known quantities.
Generalized coefficient formulas are
a - 0 (ISA)
J v ""j+l
where j •» 2, ..., n ;
- 0.5v
123
-------
where J « 1, . .. , n-1 ;
„ C5
where j = 1 ..... n ;
d . C6-X +
1 2
W
c - 0 ;
n
2 At t
Z X*
- ,. . .
, (Az- + Az., ) 1
2At
n+l (Az
n+1 n
n+1
(157)
(158)
5— + At- A (z.)
2 a 2
(159)
Ka(2n+l) _t
Az ) Az ' Xn+2
A (z ^.,
a n+17
- ^ w .. [ -. - —7 — • X^. ;
n+1 [ Az . , + Az n+2l
v '
(16Q)
fi.x
t-l
(161)
where j = 2, ..., n-1 . In the initial time step, the constants C5 = 1 ,
C6 » 1 , and 6 - 0 . In all subsequent time steps, these constants are
C5 - 1.5 , C6 - 2 , and 6 = 1 .
Table A-l gives specific symbols in our previous notation for the unde-
fined quantities in our generalized equations above. All horizontal spatial
derivatives of dependent variables in Table A-l are one-sided differences
(see Section A.13.2).
Once these coefficient values are computed, the solutions to the finite-
difference equations can be computed recursively from the formulas:
124
-------
c* - - r ; (162)
1
d* » + ; (163)
c* = cit 3 + b , (164)
where j * 2, ...» n-1 ;
d* - J ^ +^\>l » (165)
where j = 2, ..., n-1 ;
d - a d* .
n n n-l /-\e.e.\
a c* . + b ; (166)
n n-1 n
where J « n-1, n-2, ..., 1 .
A.12.2 u and v
'The finite-difference equations for u and v take the form
Au.+Bu +Cu - f (168)
where the vectors are defined:
u.
(169)
125
-------
- u
- v
6*u
t-1
j+l
t-1
(170)
where J = 2, ..., n-1 for eqs. (169) and (170);
. t-1
0 *U« IXV*.. / U- t-
C6'u,: - AffV (2) f— + . . . .—
2 g 2 Az Az2+Az
fi.vt-l K(j z x
C6-v! + AffU (2) 1— + -~
2 g 2 Azr
Az-
(171)
fi*u Kfz 1
t . -TT / ..•>. n+1 n+1.
•w.j.1 - AffV_(n+l) = +
g
C6-vV, + AffU (n+1) - — 0
n+1 g 2
Azn+l
K(2n+l>
Az
n+l
Azn+l Azn+Azn+l
- w
At
n+1 Az + Az
n+1
u
n+2
n+2
(172)
and the matrix elements
a.
0 ;
(173)
a.
2At
rK(2J
Az + Az
+ 0.5
(174)
where j « 2, ..., n ;
BR.
C5 +
2At
AZ
Az
j+l
(175)
126
-------
BG - - f'At , (176)
where j = 1, ...n for eqs. (175) and (176);
9Af. K(£ )
c. - - -r-^TI IT- ; (177)
1 Az., + Az2 Az2
CJ = - A2 ^, |-7r1Z^--0.5wJ.J , (178)
where j • 2, ..., n-1 ; and
c - 0 . (179)
n
Terms A and C. are diagonal matrices where
a. 0
(180a)
0 v'
(180b)
and B. has the form
IBR./R2. -BG./R2.
33 J J | , (180c)
-BG./R.2 BR./R2
with
Rl - BR2 + BG.2
and
R2 - BR.2 - BG2 .
j
J
127
-------
The constants* C5 , C6 , and 6 are defined as before. Quantities
for the upper and lower atmospheric and oceanic sublayers are given in
Table A-2 below.
TABLE A-2
Identification of Quantities
Quantity
K
n
£. z <.
K
m
K
m
Recursion formulas for the solutions of the finite-difference equa-
tions are: . .
E
Cl '
(182)
(183)
E,
,-1
(18A)
where j •» 2, ..., n-1 ;
(185)
u
(186)
(187)
128
-------
\
v where J « n-1, n-2, ..., 1 . The solutions to the water velocities are
subject to the boundary conditions stated in Section 2.2.A.2 of the text.
A. 13 Horizontal Advection
The horizontal array of grid points is to consist of a rectangular array
of equally spaced grid points with no limitation within the formulation of
the FORTRAN Program on the number of stations in the grid. The horizontal
gradients for each of the unknown variables, u,v,T,q,s,P,, and
?„ are computed according to the scheme outlined in the model formulation
of Gerrity (1967).
In the FREE SURFACE version of the model, the horizontal gradients of
the dependent variables are calculated every time step, but in the'RIGID LID
version, the number of basic time steps to elapse before recomputation of
the horizontal gradients is calculated In the following way. The entire
three-dimensional array is to be scanned for the maximum absolute velocity
component value, c . The number of basic time steps to elapse is then
N - : . (188)
U .01) |c |-At
A.13.1 Horizontal Advection Gradients
The horizontal gradients of a variable X. in the x (£ • 1) and
y (ft • 2) directions are computed for station j,k at each vertical grid
level n using the form (6 X /Ax )" , where ( 1 * j ,«c * 1 j ,K
129
-------
for u" > 0
(189a)
for v
(189b)
Whenever one of the grid points to be used in the computations above (e.g.,
(Xj). n) falls outside of the horizontal array, the horizontal gradients
J »u
cannot be computed and are, therefore, prescribed through input to the pro-
gram (see Section A.4). These one-sided horizontal gradients are used in
the advection terms, particularly in Section A.12.
A. 13.2 Space-Centered Gradients
Space-centered gradients are computed from (6JlXi/2Ax£) ^ where
and
k-1 C190b)
for all internal points. At boundary points where u and/or V are di-
rected into the grid, the respective space-centered gradients are not cal-
culated but are prescribed through the model input gradients. If, on the
other hand, the predicted u and/or v components at boundary points are
directed out of the grid, the space-centered gradients are taken equal to
the one-sided horizontal advection gradients which are calculated by the
formulas in Section A.13.1.
130
-------
In the FREE SURFACE version, space-centered gradients are used to cal-
culate the terms 9n/3x and 3n/3y appearing in the equations used to com
pute the geostrophic currents, viz., eqs. (99a) and (99b). The finite-
difference forms of these derivatives are (6£n/2Axn) . where
J »*
-l. k
and
The value of n at the land stations is taken as zero. If the above differ-
ences are taken at open water boundaries, n is specified outside the grid.
A. 13. 3 The Terrain Relative Vertical Velocity
The vertical velocity relative to the terrain is prescribed or computed
at each vertical grid point. In the RIGID LID version, computations are
performed as follows :
« water layer
where B <. j <. 1-1 , and wCO » 0 ..
u
o air layer
»
-------
The FREE SURFACE version uses the equations as follows:
o water layer
where B i. j <_ I , and w(z_) = 0 .
D
air layer
k=1-l jf 3u 3v
k=I+
where I <. j <_ N , and at the interface
w(zj;_) = w(zLf) - 0 . (195)
The water surface elevation at time t , T) , is given by
n* - n^1 + w(Z]:_)At . (196)
The vertical velocities and water surface elevations are computed when-
ever the horizontal advection gradients are calculated, i.e. every time step
in the FREE SURFACE version and every N time step in the RIGID LID version.
A. 14 Output
The output will be written on a tape or disk for later use. These re-
cords will show the height and time variations of the following parameters:
K eddy diffusivity (cm2/sec). •
P. pollutant 1 (yg/m3) if present.
P- pollutant 2 (pg/m3) if present.
q humidity (g/kg), z >. zI+ .
132
-------
Rl Richardson number.
s salinity (g/kg), z <. zj_ .
T temperature (°K).
u eastward wind/current component (cm/sec).
v northward wind/current component (cm/sec).
V wind/current speed (cm/sec).
w vertical velocity relative to the terrain (cm/sec).
(3T/at) infrared and solar heating rates (°K/day).
The following parameters will also be printed out. They will all be
in one 600-word record for each station:
a albedo (n.d.).
h height of wave (cm).
M moisture parameter (n.d.).
R artificial heat source (£y/sec).
6 steepness of characteristic wave (n.d.).
X wavelength of turbulent wave (cm) if water, if negative, indi-
cates land and is -ZQ , the roughness height (cm).
p soil/water density (cal/g/8K).
I (N) upward solar flux at z (mfcy/sec).
R,(N) downward infrared flux at z = ZN (mS,y/sec).
R (N) upward infrared flux at z = ZN (m£y/sec).
R (N) downward solar flux at z = ZN (mjly/sec).
R,(0) downward infrared flux at surface (mX,y/sec).
R (0) upward infrared flux at surface (may/see).
R (0) downward solar flux at surface (may/see).
E(0,t) eddy water-vapor flux into the atmosphere at the interface
g/cm2/sec).
133
-------
F(0,t) eddy salinity flux into the ocean at the interface
(g/cm2/sec).
H.(0,t) eddy heat flux into the atmosphere at the interface
(cal/cm2/sec).
H (0,t) eddy heat flux into the ocean at the interface (cal/cm2/sec).
s
i(0,t) stress at the interface (dynes/cm2).
E(zN,t) eddy flux of latent heat at the upper boundary.
F(z..,t) eddy flux of salinity at the lower boundary.
H (zN,t) eddy flux of heat at the upper boundary.
AW
H (z. ,t) eddy flux of heat at the lower boundary.
S X
t(zN,t) stress at the upper boundary.
T(z.,t) stress at the lower boundary.
The following time sums will be printed and stored on tape (units B fcy),
A(t) heat flux into the atmosphere.
b(t) outgoing infrared flux at surface.
[L-E](t) latent heat flux.
P(t) heat flux due to precipitation.
R (t) incoming infrared flux at surface.
a
R8 (t) solar flux used for heating interface.
R (t) artificial heat source.
A
S(t) heat flux into, the subsurface.
Additional printed variables for the FREE SURFACE version include:
n(x,y,t) departure of the water surface from a mean state (cm).
r\
T~~(x>y»t) slope of the water surface in the x-direction (n.d.).
3X
T^-(x.y.t) slope of the water surface in the y-direction (n.d.).
3y
-------
Appendix B
A Terrain-Based Coordinate System
The choice of a terrain relative vertical coordinate axis allows for
a convenient incorporation of topographic effects in the solutions of the
equations of motion. An additional property of this coordinate system is
its compatability with surface data observation networks, since station
observations are obtained in constant z surfaces, while horizontal mesh
distances are given in constant a surfaces (see Page 8 of the text for
a definition of z and 3 ).
The horizontal derivative of the dependent variable X , most readily
obtained from observations, may be defined (see Fig. B) as
||- = Lim |-^rHH . d97)
9X0 Ax ->- 0
where z is fixed. The relevant horizontal derivative in the coordinate
system relative to an absolute reference surface may be defined as
f - Lim |-^r^l . d98)
* Ax -»• 0
where z is fixed. Equations (197) and (198) together may be written
,y _ v v _ v .
^v IA1 An Al A1 A r I
|X . Lim M—2. - -I—*- M. . (199)
3x Ax Ax Az
Taking limits sequentially and using the definition given in eq. (197),
« - .22L _ 1*. li , (200)
3x 3x 3z 3x ' V '
135
-------
Z.J3
Ax
As = -AC
Figure B. Schematic diagram illustrating the relationship
between the gradient in a surface parallel to the
terrain (z constant), and the gradient in a horizontal
surface (2 constant).
136
-------
and similarly,
3X 9X 3C.
3yQ " 3z 3y
(201)
These horizontal derivatives could have been derived through an appli-
cation of the chain rule, viz.,
X <= X(x,y,z(x,y,s))
and
z » a -
(202)
(203)
Therefore,
3X 3X
3x 3x
ft /* ^\»» ** T* «»•*+•
3X 3z
k ,. 3z 3x
rm S*nv* ** f* f\v* +
(204)
3X
3x
_ Mil
3z 3x '
137
-------
Appendix C
Atmospheric Radiation Computations for an Urban
Meteorological Pollutant Model
M. Atwater
C.I INTRODUCTION
The purpose of these specifications is to define, in detail, the radi-
ative computations used in a generalized three-dimensional model of the
atmosphere and land or water subsurface. Section C.2 will list the symbols,
general definitions and assumptions used in the radiation model. Symbols
may vary slightly from the symbols used in other specifications given pre-
viously. Section C.3 will describe the infrared computation and Section C.4
will describe the solar radiation computations. Section C.5 will include
the effects of terrain on both solar and infrared computations. However,
it should be noted that the terrain effects have not been included in the
results of the final report. Their effects will be tested in the. future.
The vertical profiles of pressure, temperature, and humidity to level
N (level H in previous specifications) are obtained from the values at
each time step for each horizontal grid point. Input values of pressure,
temperature, and humidity are specified for all levels above level N and
are constant over the entire grid. Aerosols are also specified above
level N . Aerosols can be specified to be the first pollutant below
level N and would be computed within the model. This was not done dur-
ing the present contract period, clouds are input at any level, above or
below level N , and are linearly interpolated in time over input values.
139
-------
C.2 GENERAL DEFINITIONS AND ASSUMPTIONS
C.2.1 Vertical Grid and Symbols
The radiation computations described here use a vertical grid shown
in Figure C at each horizontal grid point. There are n levels from the
o
surface to the top of the atmosphere, and there are N (1 <_N < n ) levels
• J S
where the net infrared and solar fluxes are computed. If N >_ 3 , radiative
heating rates are computed for all levels from level 2 to N - 1 .
The symbols will be separated into two types. First 'are symbols that
have well defined meaning in relation to the radiation model. The second
type are symbols used only once or twice, and often without any specific
physical meaning. They are mainly used to make the equations more easily
read or for ease of programming.
C.2.1.1 Radiation and Model Variables
a surface albedo (n.d.).
B aerosol absorption coefficient (cm ).
cl
B aerosol scattering coefficient (cm ).
8
c specific heat of air (.24 cal/gm/°K).
c incoming solar flux on horizontal surface at top of atmosphere
v
(n.d.).
F (z) net flux at level n (Ay/sec).
g gravitational acceleration (- 980 cm/sec^).
h local hour angle (deg).
I,* solar constant (= 1.95 Ay/min).
M solar optical path length (n.d.).
n number of cloud layers.
p pressure (mb).
P concentration of aerosols (ygm/m^).
140
-------
Level
Layer
Q)
0)
0)
0)
N
ns-1
N - 1
Land or Water Surface
Figure C. Vertical grid for atmospheric radiation computation.
141
-------
q specific humidity (gm/kgm).
R gas constant for dry air (2.87 x 106 erg/ g/°K).
R
-------
e emissivity due to water vapor (n.d.).
Cp fractional cloud amount for cloud layer SL (n.d.).
H azimuth of terrain slope (deg),
K infrared absorption coefficient for aerosol (cm ).
—12 u
a Stefan-Boltzmann constant (1.354 x 10 i,y/°KH/sec).
T transmissivity due to carbon dioxide (n.d.).
T transmission between levels n and m (n.d.).
n,m
r(n,m) value in transmissivity matrix (n.d.).
<|> latitude (deg).
iji transmissivity of solar flux through a cloud of given type
with a given path length (n.d.).
3H/8x, .. fc . _ , , .
su/1 f gradients of terrain H (n.d.).
on/oX ' . •
9T/3t radiative heating rate for layer (8K/sec).
C.2.1.2 Temporary Symbols
The following symbols have little, if any, significant meaning other
than to make the equations more readable and easier to program. Therefore,
only the variables will be given and the equations that define the variables.
Variable Equation Defined
Am 11(+4) * * Where eq. 11(+4) is
the fourth equation
E 15 after equation (11).
This procedure in
F 11(+10) identifying equations
follows throughout
G 11(+8) Appendix C.
R* Id, 4(+l)
R t 6d, 8
R t 14
s
T 1K+6)
T2 11(+7), 12(+3)
143
-------
Variable Equation Defined
Ta 13
u 7 (+2), 11 (+1)
cl
1C-3)
1C-3)
V 10 (+3)
W 11 (+9)
Y 11 (+5)
e 5, 3, 6c, 8, Ic
a
e Ib, 3(+l), 7(+l), 6b
c
C.2.2 General Assumptions . ,
The radiation model described here is designed for use within the bound-
ary layer of the atmosphere. Its applicability is probably valid within the
troposphere.
The radiative fluxes, both upward and downward solar and infrared, are
computed at each time step. Both solar and infrared use empirical equations
to cover their respective parts of the .spectrum. In addition to standard
absorbers and scatterers, clouds and pollutant or natural aerosols are in-
cluded in the model. Using appropriate radiation coefficients, a pollutant
gas can be added to, or substituted for, the aerosols. The influence of
clouds is estimated using the following assumptions.
1) The amount of coverage of two or more layers is equal to the
product of the individual coverage. This is equivalent to assum-
ing a horizontally uniform distribution of clouds.
2) The thickness of the clouds below level N is assumed to be one
layer.
3) The clouds are black for atmospheric radiation.
144
-------
The method used here follows that outlined by Manabe and Strlckler (1964).
C.2.3 Optical Path Lengths
Increments of the optical path length for each of the absorbers and
scatterers are computed for n - 1 layers of the atmosphere.
8
The path length increments for water vapor are computed for
Au
2g
The path length increment for carbon dioxide is
Auc± = .4148239
which assumes 320 ppm for the atmospheric carbon dioxide.
The path length increment for aerosols is computed from
)±
where P is the aerosol concentration. Above level N , the layer thick-
ness needs to be computed. This is done by using the hydrostatic equation,
which results in
zi+i - zi = t ?T (pi " pi+i} •
In all cases, the path length increments are computed from the surface
i = 1 to the top layer of the atmosphere, i.e. i = n - 1 . The layer
s
thickness is only needed for i >. N .
145
-------
C.3 INFRARED RADIATION
C.3.1 Transmissivity Matrix
A transmissivity matrix for dimension N by n - 1 is computed at
8
certain specified intervals, currently six hours, to reduce the computer
time. The indices used with the transmissivity matrix will be m with a
range of 1 to N and n with a range of 1 to n - 1 .
s
For each level from 1 to N , downward and upward infrared radiation
will be computed. As the transmission is a function of the total path length
from the level m to a given layer n , the path length is computed first,
and then the transmission is computed. The total path length for water
vapor u is computed from
n-1
u = ) AVL. , n >. m
n , i
i=m
m-1
un " I A"*! » n < m
i=m
Identical equations except for the replacement of Au with Au are
used for the total path lengths for carbon dioxide. From the total path
length, the transmission (or emissivity) is computed. The transmissivity
(or emissivity is computed from equations given in Table C-l for water
vapor and carbon dioxide.
Then each value of the transmissivity matrix is computed from
T(n,m) = .815 - e + T
c •.
The above procedure is repeated for all valid values of m and n
except m = n . The matrix i(n,m) is used for later computations.
146
-------
TABLE C-l
Infrared Diffuse Transmission Functions
Water Vapor
Functions derived from data of Kuhn (1963) for
path length (u in precipitation - cm).
e •- 0.11288 log (1 + 12.63 u) *
e - 0.104 log (u) + .440
e = 0.121 log (u) + .491
e = 0.146 log (u) -»• .527
e -, 0.161 log (u) + .542
e = 0.136 log (u) + .542
Carbon Dioxide
Presented by Shekhter and revised by
Kondrat'yev (1969) (u in atm/cm),
e*p(-.3919 u >
Log u
<_ - 4
£ - 3
1 -'1.5
<_ - 1
£°
> 0
* Presented by Zdunkowski arid Johnson (1965).
C.3.2 Downward Infrared Radiative Flux
Computations for the downward infrared radiative flux are performed for
m levels starting at the surface. The initial computation at each level m
is to determine the first cloud level above level n at which a non-zero
cloud amont occurs, and this level is saved.
The following values are initialized as:
i(m,m) = 1
e «= 1
C(la)
C(lb)
147
-------
e^ - 1 C (lc)
R4- - 0 C (Id)
R 4- - 0 C (le)
c
Then for a given level m , all layers above m to the top of the
atmosphere (n = n -1) will be used sequentially to compute the downward radi-
s
ation at the level m . Therefore, the following computations are repeated
for each level n until n = n starting with n = m+1 . It should also
5 •
be noted in the folloiwng equations where the same term appears on both
sides of an equal (=) sign that "=" does not mean equal but means "is replaced
by" as is standard in FORTRAN programming.
The optical path length for aerosols is computed from
u - * Au . C(2)
If a cloud base exists at level n-1 , then
R 4- - R + + aT1* e, e (e ) C (3)
c c n A c a n_i
where e. is the amount of cloud for layer I . Then
and the £ is then incremented by 1 for the next cloud level, if any.
Then the transmission is computed using values of the transmissivity
index or
T = e"-* T(n,m) . . . C (A)
n , m
Then, R+ is incremented by
T + T
R4- = Ri + a -S - e e
c a
148
-------
e = T , - T C (5)
a n—j.,m n,m
All computations from C (2) to C (5) are repeated until n = n . Then
s
the downward flux at level m is computed from
R KZ ) = R4- + R 4- .
an c
The result is stored. Then m is incremented by 1, and the computations
are repeated from eqs. C (la-e).
C.3.3 Upward Infrared Radiation
For each level m > 1 , the following variables are initialized as
T = 1 C (6a)
ID, in
e = 1. C (6b)
c
efl = 1 C(6c)
cl
R+ - 0 C (6d)
R * = 0 C (6e)
C
Then for each level n starting with n = m-1 and ending with n = 1 ,
the following procedure is used. .
Find the first cloud base at a level n < m and set '4 , if any are
found. Then, if the cloud base is at level n , compute
R t = R t + e. 0T1* e e C (7)
c c inac
Then a new e is computed as
c
e = e (1 - e.) ,
c c £ '
and one is substracted from i for the next lower layer of clouds.
149
-------
The path length for aerosols is computed from
i=n
Then the transmlsslvity is computed from
.a_. , N
T =• e d m T(n,m)
n,m
and
a n+l,m n,m
Then
.+ T i1*
R t = R i + 01— — e e . C (8)
a a l2Jac
fn
.-
Then n is decreased by 1, and, if n > 0 , all equations from C(6a-e)
through c (8) are repeated. If n = 0 , the following is computed:
R Hz) <= RQ + + Rt + T. e^ aTj . C(9)
am a c l,m c I
Then, the procedure is repeated for all levels m starting with eqs.
C (6a-e).
C.3.4 Infrared Radiative Cooling Rate
The infrared cooling rates are computed for each layer when three or
more levels have Rt(z) and R4-(z) computed. Then the cooling rates are
interpolated from layered cooling rates to m levels.
The net flux F (z) is computed for each level from
F (z.) - R KZ,) - R i(z.) ; 2 < i < N .
ni ai ai ~~
The cooling rates for the layer are computed for layer i , i >. 3 ,
150
-------
3T
3t
__?__ fW - Fn(«i-l>]
V103 i Pi - Pi-! J '
For level 2, the first level above the interface,
3T
Ut
For levels of i > 2 , the cooling rates are computed from
P±
T
ar
C (10)
Storage will be saved if F and 9T/3t are saved for only levels
and layers of 14-1 and i .
C.4 SOLAR RADIATION
C.4.1 Downward Solar Radiation
The downward solar radiation is based on methods developed by Kondrat'yev
(1969), McDonald (1960), Manabe and Strickler (1964), and Atwater (1970).
Compute the cosine of the zenith angle from
oosZ » sin 6 sin$ + cos5 cos$ cosh .
If cosZ <. 0 , the sun is on or below the horizon. Therefore, in this case,
(s2)t - o ,
a - 1 ,
where 1 <. i <. m , and no further computations for solar radiation are made.
Because of numerical problems as cosZ -»• 0 , a lower limit is set on
151
-------
oosZ and the radiation is reduced by a fraction, V , defined as follows:
V - 1 f .
and oosZ is unchanged if cosZ >_ 0.17365 ;
cosZ
V
0.17365
and cosZ is set to 0.17365 if oosZ < 0.17365 .
The initial flux for a horizontal surface at the top of the atmosphere
is computed from
c = VIQ oosZ .
The flux is then reduced by clouds and aerosols above level N .
* * -
The solar optical path length is computed from
p<0
M • -
1013 cooZ '
The transmission t|> is then computed for all cloud layers above level N
with j being the first cloud layer above N using the equations given in
Table C-2, and using the cloud type as specified on input. Then
If n < j , T^l . TC_ is also the fraction of clear sky above level
c c TD
N .
For each level m , the following are computed. The path length for
water vapor from level m is computed from
ng-1
um * I Auwt • C (11)
i=m
152
-------
TABLE C-2
Transmission equations for clouds.
Cloud type
1.
2.
3.
4.
5.
6-
7.
8.
9.
Fog
St
Sc
Cu
Cb
As
Ac
Ci
Cs
Equation
* \Ji ° 16
1JJ a 26
4» = 36
t|» = 36
* - 23
il/ «> 41
. * = 54
i|i = 87
4» • 90
.26
.84
.58
.58
.63
.30
.56
.17
.55
+ 0.
- 1.
- 1.
- 1.
+ 1.
- 0.
- 2.
- 1.
- 6.
54
01
49
49
45
14
36
79
38
M
M
M
M
M
M
M
M
M
* if) is in percent.
For aerosols, the path length is
u,
>am
ne-l
i=m
Further, T and T are computed for the layers above level m by
cl S
and
Then
m
s r_
m TH
The downward flux at level m is computed starting at level N and
ending at level 1 . For each level m , the following procedure is re-
peated. First compute
153
-------
Y - 1.041 - 0
[0.000949P(z ) .+ ..051~|*
—
CO8Z J
T - .485 + .515 Y
[u 7s
—HT
cosZJ
The values T , T , and T are computed for the region above the layer.
C Si S
Then the variables at level m are computed from
G c T- ,
m v 1 '
and
W «= - c T- ,
m v 2 '
F = A .- (G + W ) .
m m+1 m m
However, when m is equal to N , the variables A Ln are the same as for
m+1
level N . Then the flux at level N is given by
W = Fm • C(12)
and further computations repeat from eq. C (11) if N > 1 after one has
been substracted from m . If N = 1 , the incoming interface flux has
been computed in C (12) , and no further computations are necessary for down-
Ward solar radiation.
For N > 1 and m < N , the computations proceed as follows. If a
cloud base i is at level m , the new transmission for total clouds is
computed from
where i|>. is computed according to Table 2 for the cloud type i . If
154
-------
there is no cloud base at level m ,
Lcm
cm+l
Then compute
m
V
Then the downward flux at level m is computed as
Rs(zm) = T«m (1 ' T«— + T«-}
C(13)
C.4.2 Solar Radiative Heating Rates
The basic assumption made is that only downward radiation flux will
have an influence on radiative heating. The net flux absorbed in a given
layer is
F (z ) - A ... (W ., - W )
n m m+1 m+1 m
m
-V
and
3T
3z
g
m = 103cp
For layers with m <, N - 2 , the radiative heating at the level m+1 is com-
puted from
3T
3z
m
P(zm+2) - p(zm)
[(HI
13ZJ
m+1
3T
155
-------
C..4.3 Upward Solar Radiation
The basic assumptions included here are that:
1) The upward flux of solar radiation does not influence the
radiative heating rates.
2) The upward flux is unaltered by the atmosphere.
The upward solar radiation is composed of a) the radiation reflected
at the surface, and b) the radiation scattered upward by atmospheric levels
below level N (if N > 1 ) . The upward flux at the surf-ace is computed
from
Rt = a R (0) . C (14)
O o
The albedo is computed for oceanic conditions from results of Sverdrup et
al. (1942) given by
a - - 0.0139 + .0467 tcmZ
with limits given by
0.03 i a i 1.00 .
The albedo for inland lakes uses data obtained by Nunez et al. (1971) and
is dependent on the total cloud amount. The total cloud amount is given by
nc
E = 1 - n (1 - e.) . C'(15)
i=l i
When E <. 0.5 , the albedo is computed from
When E > 0.5 , the albedo is computed from
156
-------
-00462 z ; a ^ 1-° •
Over land, the albedo is specified as an input item.
The upward radiation from layer m is computed from
AR t(z )
s ra
A ., (G .. - G ) + 0.9 F (Tc ., - T. )
m+1 ro+1 m m Hn+1 °m
The total upward solar radiation at level N is
N-l
= a R (0) + 7 AR (z ) .
s *•, s m
n=l
C (16)
C.5 TERRAIN EFFECTS ON RADIATION
The angle of the slope is computed from the input gradients of slope.
The angle B is defined as the angle between the horizontal and the slope.
It is also the angle between the normal to the slope and the zenith-point.
The angle B is computed from
B
tan
-i
The azimuth of the alope n is computed from
3 . -l 3H
3* - tan -
If n > 2ir , then substract 2ir from n for a corrected n • Compute the
solar azimuth a from
157
-------
-cos 6 sirih
and check to insure proper quadrant. Then the cosine of the incident angle
is computed from
oosi. = cosB oosZ + sinB sinZ cos (a - n)
The solar radiation is then modified to
The effect on infrared radiation is then computed from
R£ (0) = R+(0) cosB .
Cl
C.6 REFERENCES
Atwater, M.A., 1970: Investigation of the Radiation Balance for Polluted
Layers of the Urban Environment. Ph.D. Thesis, New York University,
New York, N.Y., 116 pp.
Kondrat'yev, K., 1969: Radiation in the Atmosphere. Academic Press, N.Y.,
912 pp.
Kuhn, P.M., 1963: Radiometeorsonde observations of infrared flux emissivity
of water vapor. J. Appl. Meteor.t 2t. 368-378.
Manabe, S. and R.F. Strickler, 1964: Thermal equilibrium of the atmosphere
with a convective adjustment. J. Atmos. Soi.f 21t 361-385.
McDonald, J.E., 1960: Direct absorption of solar radiation by atmospheric
water vapor. J. Meteor., 17t 319-328.
Numez, M., J.A. Davies, and P.S. Robinson, 1971: Solar Radiation and Albedo
at a Lake Ontario Tower Site. Third Rpt., Dept. of Geography, McMaster
University, Hamilton, Ontario, Canada, 92 pp.
Sverdrup, H.V., M.W. Johnson, and R.H. Fleming, 1942: The Oceanst Prentice-
Hall, Inc., N.Y., 1087 pp.
Zdunkowski, W.A. and F.G. Johnson, 1964: Infrared flux divergency calcula-
tions with newly constructed radiation tables. J. Appl. Meteor., 4,
371-377.
158
-------
Appendix D
The Non-Hydrostatic Model
The non-hydrostatic three-dimensional model differs only slightly from the
hydrostatic model described in Section 2.0 of the report. This difference is
manifested in the computation of the thermal wind. In the non-hydrostatic model
the thermal wind equation is derived from the third equation of motion with the
geostrophic assumption and the equation of state for an ideal gas. This deri-
vation yields the equations
.H
8T 1 [ 3
^dz + f J U
= Vg(H) - f^Tl Sfd« + £| £:(3=-]dz
.H
— dz - v
x f
z
where 3T/3x and 3T/3y are computed according to equations (11) and (12) of main
text [also, see equations (95a), (95b), and (95c) in Appendix A for the finite-
difference form]. The above equations replace equations (9) and (10) (see
Section 2.0 of the report). The essential difference between the two models in
the computation of the thermal wind is in the addition of the second integral in
the above equations. This term represents the horizontal gradient of the vertical
acceleration.
In the finite-difference form the integrals are approximated through an
application of the trapezoidal quadrature formula which was used in equations
(94a) and (94b) of Section A.9.1 of Appendix A. The horizontal gradients of the
vertical acceleration were approximated by spatially centered differences (see
Section A.13.2).
Preliminary tests with the non-hydrostatic model revealed that on the coarse
grid (—13 km) used in the Los Angeles simulations, the term involving the gradient
159
-------
of the vertical acceleration would be approximately two orders of magnitude less
than the term involving the temperature gradient, and, therefore, the differences
between the hydrostatic and non-hydrostatic solutions would be small. Subsequent
simulations on the Los Angeles coarse grid have revealed significant differences
between the two models. We have attributed these differences to the topographic
influence on the vertical acceleration. There have been comparisons of the
solutions from the two models under identical input conditions. These compari-
sons have shown a tendency toward better temperature verifications but slightly
worse pollution verifications. Our investigations with the non-hydrostatic model
are continuing.
160
-------
Appendix E
Advection Modifications
Numerical approximations to advection problems generally produce "numer-
ical diffusion" of the solutions. This results from the uncertainty involved
in attaching values to the dependent variables only at grid points. Thus,
trajectories which do not pass directly from one grid point to the next do
not explicitly define values of the dependent variables in regions where they
are passing between grid points. This problem is resolved in Eulerian finite-
difference schemes by various methods of interpolating or extrapolating re-
sults from points of explicit solution to off-trajectory points. Since the
adequacy of an interpolation technique depends on distance between known
points, the accuracy of interpolation (and thus the magnitude of the "numeri-
cal diffusion" effect) depends on the wind azimuth (see Figure El), while
for space- or time-variable winds the diffusion effect is also variable.
Numerical tests of the azimuthal effect were carried out. The results,
reproduced in Figures E2 through E6, are presented in the form of concentra-
tion maps for emission from a single source in a steady, uniform windfield
for each of several azimuths. It may be seen that for winds down the grid
line, no diffusive effect is present, while for windfields with two finite
vector components, every grid point in the quadrant downwind of the source
receives some material. The steepness of the concentration ridge as well
as the location and value of the concentration maxima vary with wind azi-
muth.
161
-------
.Source
.02,3
3,1 O
O3,3
Fig. El. Advectlve transfers for conventional Eulerian scheme
with three wind trajectory azimuths.
Trajectory 1: All material advected to point 2,1. Axial
diffusion only. Point 2,1 is only source for
further transfer. Material is confined to
grid-line trajectory.
Trajectory 2: Most material to point 2,1; some to point 1,2,
Intermediate lateral diffusion. 2,1 and 1,2
become sources for further transfer. Material
spreads to entire quadrant.
Trajectory 3: Material evenly split between points 2,1 and
1,2. Maximum lateral diffusion 2,1 and 1,2
become sources for further transfer. Material
spreads to entire quadrant.
162
-------
-+ 4 4
4 -t- -t- +
4 4 -f -t- +
4 4 4- + +
-t 4 ' 4 -+• 4-
4 4- + -*- 4
4 -t 4- 4 4
-I- -t- -I- 4
-l- -t -+
44- +
+ 4-1-
4 -1- +
+ 4- + H
= .30
Fig. E2. Concentration contour map for constant emission
from a point source using upwind differences on a
Eulerian grid. Wind from ~276° (~W.)
x=.25
Fig. E3. Concentration contour map for constant emission
from a point source using upwind differences on a
Eulerian grid. Wind from ~284° (~WNW by W)
163
-------
Fig.
E4. Concentration contour map
for constant emission from a
point source using upwind
differences on a Eulerian
grid. Wind from ~293° (~WNW)
Fig.
E5. Concentration contour map
for constant emission from a
point source using upwind
differences on a Eulerian
grid. Wind from ~304° (~WNW
by N)
x = .25
.3Q Cj * = .2Q
-t-
x = .20
x=.15
Fig.
E6. Concentration contour map
for constant emission from a
point source using upwind
diferrences on a Eulerian
grid. Wind from 315° (NW)
x=.15
x=.10
x=.05
+ +-+--I--I--I- + -4-
_J i i i I i | i
x=.10
x=.05
x = .05
164
-------
E.I Eulerian Modifications
Several techniques for minimizing these effects have been examined.
One Eulerian technique appears plausible but has not been programmed for
trial. In this scheme, the only change from conventional finite-differencing
techniques would be in allowing advective transfers between diagonally
neighboring points as well as along grid lines. That is, convective (pre-
sumably upwind, for stability) differences could be taken with one or all
grid indices incremented by one or more steps. As shown in Figure E&, if
one-cell diagonal transfers are allowed, numerical diffusion can be reduced
so as to contain material within one half of a quadrant. If transfers are
allowed across two-grid cells, as shown in Figure E8, the spreading can be
reduced to a 2TT/16 sector (22-1/2°).
If it is assumed that winds cannot be defined to closer, than a 16-point
compass, the Eulerian scheme can be modified further to completely eliminate
numerical spreading of an advected plume by allowing finite-difference trans-
fers to only the single downwind neighbor that lies closest to the actual
trajectory (see Figure E9). This will tend to preserve the concentration
magnitudes to be expected downwind and would allow the explicit inclusion
of correct lateral diffusion but would give somewhat less precise definition
of the trajectory and thus of the locations of downwind maxima.
Such Eulerian schemes have equal or greater computing efficiency than
conventional schemes in terms of the actual integration additions required,
but would require considerable additional logical control steps. Time esti-
mates have not been made.
165
-------
1,1 @
Source
01,2
2,2
Ql.3
O 2,3
O 3,3
Fig. E7. Advective transfers for modified Eulerian scheme
allowing transfers to diagonal neighbors. Three wind
trajectory azimuths shown.
Trajectory 1: All material advected to point 2,1. Axial
diffusion only. Point 2,1 is only source for
further transfers. Material is confined to
grid-line trajectory.
Trajectory 2: Material split between points 2,1 and 2,2.
Maximum lateral diffusion. Points 2,1 and 2,2
become sources for further transfer. Material
spreads through 45° sector.
Trajectory 3: All material advected to point 2,2 axial
diffusion only. Point 2,2 is only source for
further transfer. Material is confined to
diagonal trajectory.
166
-------
Source
01,2
O1.3
O 2,3
Ol,4
02,4
O3,4
Fig. E8. Typical allowable advective transfers for a 16-point
(two-cell) modified Eulerian scheme. Five wind trajectory
azimuths shown.
Trajectory 1:
Trajectory 2:
As in Figs El and E7.
Material split between points 2,1 and 3,2 which
become sources for further transfers. Material
confined to 22-1/2° sector.
Trajectory 3: All material advected to point 3,2 which becomes
only source for further transfers. Axial diffusion
only. Material confined to half-diagonal trajectory.
Analogous to trajectory 2. Material confined to
22-1/2° sector.
As in Fig E7.
Trajectory 4:
Trajectory 5:
167
-------
O1.3 Ol,4
O2,4
O3.4
Fig. E9. Three typical trajectories in a 16-point nearest-grid-
point Eulerian advection scheme.
168
-------
E.2 Lagrangian Advection
Lagrangian particle-in-cell techniques are substantially free of numeri-
cal diffusion because positional resolution of advected material depends on
the numerical accuracy (significant figures remaining) maintained in the
computations and not on the grid definition.
Since we are concerned with the solution of a coupled set of differen-
tial equations, and since advection is a significant or even dominant effect
for each dependent variable, a generalized Lagrangian scheme was developed
and tested that
° allows different resolution schemes for each dependent variable
(resolution is basically determined by the definition of the
sources; and mass, momentum and energy sources, or sinks, are
defined on scales consistent with available data and expected
influence on resultant fields);
o is compatible with Eulerian techniques for treating explicitly
included three-dimensional diffusion terms ;
o substantially eliminates numerical diffusion; and
o is easily programmable.
The computer algorithm was developed by treating the conservation equa-
tions in the general form
where d/dt is the total (substantial, or convective) derivative; x^ ^x
the i(th) conserved quantity; 41' is a source term, a function of space or
time that may depend on x^ and/or x^ (i-e., X^'s are other dependent
variables); K. is a three-dimensional diffusion coefficient for i(th)
variable, and V is the gradient operator.
Integrated, in finite-difference form, we have
169
-------
) E(2)
where iji equals (|>'At ; and n is the time step index taken along a tra-
jectory. That is, x is the distribution of x a time step ago when it
was a space step upwind. Similarly, ^n is the source distribution then,
etc. Since X- is» i-n turn, dependent on i|>(x ) , etc., the algorithm
can be reduced to
JO *V + JO
so that x at any time is the sum of "emissions" from one step upwind one
time step ago, plus the "emissions" two time steps ago from two steps upwind,
etc.; (Ji may consist of any number of terms including the diffusive term.
The diffusive term is written separately, however, to indicate its explicit
representation in the model.
In equation E(3), ^ becomes the inventory of sources for the pollu-
tion equation, E(l), the pressure gradient (geostrophic wind) and buoyancy
terms for the momentum equations, the solar and infrared heating for the
temperature equation, E(l), etc. At each time step the diffusion term and
the resultant x fields are computed on an Eulerian grid, but the ty , or
"source", terms are accumulated as sets of advected points with the exact
location of each "particle" computed. Since x is eventually swept out of
the field of interest (e.g., Connecticut SO- eventually moves out to the
Atlantic), the i|> sets are screened to retain only "particles" of continu-
ing interest. To avoid excessive recomputation of "particle" positions,
the "particle" positions are noted relative to the current position of the
first emitted "particle" from the same source. Interpolation is then in-
volved in translating the advected cloud of particles onto the Eulerian,
170
-------
fixed grid. This interpolation, which would produce numerical diffusion
in some techniques does not accumulate positional error in the present
tecniques; thus, the only numerical diffusion present is that for a single
(the latest) time step.
A different Lagrangian scheme (PICK) has been advanced for air pollu-
tion modelling (see ref. Sklarew, et al., 1971 of main text). Differences
between the present scheme and PICK are the following.
o The present scheme is presented as a general scheme applicable to
coupled sets of differential equations. PICK was presented as
a limited algorithm for the solution of an equation for advected
mass.
o In the present scheme, the dependent variables are carried as
continuous variables, not discretized as in PICK; thus, it has
a significant advantage over PICK in computational efficiency
and/or accuracy. In PICK, trajectories must be computed for
each "particle", thus the problem is effectively carried out
on a much finer grid than the answers are available on.
0 Both systems carry location along the trajectories as a con-
tinuous variable, but the advective velocity is defined only
at grid points in the present scheme, whereas it is pseudo-
continuous (by interpolation) in PICK. There seems to be
little advantage in the velocity interpolation since it is dif-
ficult to parameterize significant grid scale variability in
the velocity as turbulence. In practice, windfields should be
defined on grid scales such that the point-to-point variability
is small.
Initial trials of this technique were limited to non-diffusive advec-
tion from a steady point source, to compare with the similar tests with the
Eulerian technique discussed above. The first trials with the Lagrangian
technique used linear interpolation to translate between the moving and the
fixed grids. This proved to be highly non-conservative in the resultant
field because the concentration peaks that occurred between fixed grid
points were not represented anywhere; subsequent trials used nearest grid
point translations, so that each particle in the Lagrangian set was assigned
to a grid point in the resultant fixed grid.
171
-------
These results are shown in Figures E10 through E12. It may be seen
that the system is conservative and as non-diffusive as may be (trajectories
at an angle to grid lines simply cannot be represented by a single straight
line of points; there is necessarily a "staircase" effect). The concentra-
tions do differ for different trajectory azimuths, however, This is because
the linear density of a set of diagonally joined squares is less than that
of a laterally joined set. Over a swath of two or three cells width, which
is about all the resolution that can be expected from any grid solution, the
mean concentrations are identical. It is important to note that this accen-
tuates the difficulty of "verifying" any grid model output with point obser-
vations, even though the model may be more "correct" and/or relevant.
Further tests showed that circular (or spherical) arrays of sources
give off plumes that have essentially no azimuthal variation of center-line
concentration. Similarly, solutions for single source emissions resolved on
a grid defined as a matrix of packed spheres give no azimuthal effect. It
is not presently felt that this effect is important enough for dense arrays
of sources (or pseudo-sources for momentum, energy, etc.) to warrant further
exp er imenta t ion.
Introduction of the Lagrangian advection scheme described above pro-
ceeded by stages. Because of the complexity of converting the full urban
boundary layer model to the Lagrangian form, initial tests were carried out
with the two-dimensional "Connecticut Model".
Several logical problems appeared and were solved. Some of the problems
were compounded by the necessity, because of core space limitations, of
treating the grids by successive segments. Some of these problems were
the following:
172
-------
;Source
f- _—
.Ml M .M* .Ml .Ml .Ml *i
Fig. E10. Concentration contour map for ^Hj=.-=
constant emission from a point source
using Lagrangian advection scheme with ^l
nearest grid-point interpolation - fully ^JT— njn
conservative. Wind from 360° (N) .—i--- —
-— Source —-
Fig. Ell. Concentration contour map for ;"".*"'
constant emission from a point source fH!r
using Lagrangian advection scheme with j~:~
nearest grid-point interpolation - fully ——— ——V
conservative. Wind from 117-1 I1)0 (UNU^ 2Zfl*Z.~~*V
conservative. Wind from 337-1/2° (WNW)
Fig. E12. Concentration contour map for -——
constant emission from a point source ~H;~;r
using Lagrangian advection scheme with -.»..»..-,
nearest grid-point interpolation - fully -"• — — — — •- — —— •-••^«jr'-~ —— — — — — ——
conservative. Wind from 315° (NW) j-.-.-.—~.-r ------- -..-_r.-^:" - - •- •-
173
-------
o masking, at each time step, the source points whose trajectories
terminate outside the resultant grid as well as resultant grid
points which terminate trajectories which originate outside the
source grid;
o handling losses to the upper air in zones of horizontal conver-
gence;
o advecting concentration fields as well as emission fields when
wind fields vary in time;
o increasing computational efficiency; and
o advective instabilities introduced by vertical loss mechanism.
The two-dimensional Lagrangian model has been completed and tested. The
computer output for a test run is presented in Figure E13. The horizontal
numerical diffusion, characteristic of Eulerian schemes, has been found to
be substantially eliminated as expected. The absence of such diffusion led
to very high pollutant values at grid squares with very large point sources
(power plants), but this was reduced to expected levels with the introduction
of upward losses modelled by the thermally generated updrafts which are used
in the mesoscale windfield analysis to compute horizontal winds. Computer
time used is approximately 10m per second per grid point per time step with
an additional 20 seconds for initial computations and 20 seconds per change
of meteorological conditions.
Problem solutions for these two-dimensional tests are applicable to the
three-dimensional coupled meteorology Lagrangian model.
The Lagrangian integration scheme for the complete set of differential
equations of the urban boundary layer model was programmed in ANAL70. The
programming of input-output, boundary conditions, and auxiliary equations
was not done because early test runs with the operating Eulerian model
showed that realistic results were obtainable in spite of the diffusive
properties of the model. In fact, CO verifications were comparable to those
174
-------
West
-—Y-*•*'•—*t —"— '—
South
>tt
Bridgeport
IM .!•• .!•• .HI .!•• .IK .IM .!•• .Ill
Sound -*-*~
^U »ft» «*M .>»» ,C« .IM-^M .11
.•M .»«« .« .Jit »3M .ICCt.lir/.MO .101
W^
North
East
Fig. E13. Excerpt from computed concentration map showing
South-Central Connecticut.
• Lagrangian Advection.
• Southwest wind (1 m/sec) veering to WSW,
• S0» emissions, January.
• 5000 ft grid spacing.
• Topographic and thermal wind disturbances.
175
-------
obtained with the two-dimensional PICK method and a much finer grid (see
Section 3.3.4 of the main text).
176
-------
BIBLIOGRAPHIC DATA
SHEET
1. Report No.
EPA-R4-73-025a
3. Recipient's Accession No.
4. Title and Subtitle
Tests of an Urban Meteorological-Pollutant Model Using CO
Validation Data in the Los Angeles Metropolitan Region
(Volume I)
5' Report Date
May 1973
6.
7. Author(s)
Joseph P. Pandolfo and Clifford A. Jacobs
8. Performing Organization Rept.
No. 49Q-A
9. Performing Organization Name and Address
10. Project/Task/Work Unit No.
4121
The Center for the Environment and Man, Inc.
275 Windsor Street
Hartford, Connecticut 06120
11. Contract/Grant No.
68-02-0223
12. Sponsoring Organization Name and Address
Meteorology Laboratory
National Environmental Research Center
Research Triangle Park, North Carolina 27711
13. Type of Report & Period
Covered
Final 9/71 - 2/73
14.
15. Supplementary Notes Supplement: Volume II, "FORTRAN Program and Input/Output
Specification" by J.A. Sekorski
16. Abstracts The urban boundary-layer model, described in a previous report (Pandolfo et al.,
1971), was modified and used in forty test runs. Many.of the runs varied the meteorological
input about a standard (observed) set. It has, therefore, been demonstrated that an economical,
objective, physically consistent, and precisely specified (though with some arbitrary elements)
procedure has been achieved for obtaining and predicting the 3-dimensional meteorological fields
needed. In several of the runs, the input topography, land-water distribution, and other physical
characteristics of the underlying surface was varied. The results demonstrate that ready gen-
eralization to other regions can be expected.
The modelled region was simulated with relatively coarse (8-mile) grid spacing.
This is in contrast to other models which deal with pollutants only, and which are based on
2-mile grid spacing (Sklarew, et al., 1971; Roth et al., 1971). Nonetheless, the temporal and
spatial variations of air temperature, humidity, and wind, are simulated with an encouraging
degree of realism. Temporal and spatial variations of CO are also simulated fairly realistically,
with somewhat less accuracy than in the model described by Roth, et al. (1971), and with accuracy
equivalent to that shown by Sklarew, et al. (1971). It is reasonable to expect improved simula-
tion accuracy with the finer horizontal resolution used in these other models, and the perform-
ance of such simulations with this model is strongly urged.
17. Key Words and Document Analysis. 17o. Descriptors
Numerical Models
Air Pollution Meteorology
17b. Identifiers/Open-Ended Terms
Urban Boundary-Layer Meteorology
17c. COSATI Field/Group 04/02
18. Availability Statement
19. Security Class (This
Report)
UNCLASSIFIED
20. Security Class (This
Page
UNCLASSIFIED
21. No. of Pages
189
22. Price
FORM NTIS-33 (REV. 3-72)
177
USCOMM-OC MBS2-P72
------- |