EPA-650/4-74-045-0
SEPTEMBER 1974
Environmental Monitoring Series
K*:*
i
at
O
iiiiiiiij^
-------
EPA-650/4-74-045-H
SELECT RESEARCH GROUP
IN AIR POLLUTION METEOROLOGY,
SECOND ANNUAL PROGRESS REPORT:
VOLUME I
by
Select Research Group
Department of Meteorology and
Center for Air Environment Studies
The Pennsylvania State University
University Park , Pennsylvania 16802
Grant No. R-800397
Program Element No. 1AA009
ROAP No. 21 ADO
Task No. 14
Project Officer: Kenneth L . Calder
Meteorology Laboratory
National Environmental Research Center
Research Triangle Park, North Carolina 27711
Prepared for
OFFICE OF RESEARCH AND DEVELOPMENT
ENVIRONMENTAL PROTECTION AGENCY
WASHINGTON, D.C. 20460
September 1974
-------
This report has been reviewed by the Environmental Protection Agency
and approved for publication. Approval does not signify that the con-
tents necessarily reflect the views and policies of the Agency, nor does
mention of trade names or commercial products constitute endorsement
or recommendation for use.
11
-------
Ill
TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing)
1 REPORT NO.
EPA-650/4-74-045-a
3. RECIPIENT'S ACCESSION-NO.
4 TITLE AND SUBTITLE
Select Research Group in Air Pollution Meteorology,
Second Annual Progress Report: Volume I
5. REPORT DATE
September 1974
6. PERFORMING ORGANIZATION CODE
7 AUTHOR(S)
Select Research Group
8. PERFORMING ORGANIZATION REPORT NO
9 PERFORMING ORGANIZATION NAME AND ADDRESS
Department of Meteorology and Center for Air
Environment Studies, The Pennsylvania State
University, University Park, PA 16802
10. PROGRAM ELEMENT NO.
1AA009
11. CONTRACT/GRANT NO.
R-800397
12. SPONSORING AGENCY NAME AND ADDRESS
Environmental Protection Agency
National Environmental Research Center
Meteorology Laboratory
Research Triangle Park, North Carolina
13. TYPE OF REPORT AND PERIOD COVERED
Annual Progress 6/1/73-9/30/74
14. SPONSORING AGENCY CODE
27711
15 SUPPLEMENTARY NOTES
Issued as Volume I of 2 Volumes
16 ABSTRACT
Progress reports are included by the SRG task groups involved in: the
development of mesoscale air pollution related prediction models, modeling of
planetary boundary layer (PBL) turbulence and structure, the analysis of acdar
signals for wind and temperature measurements in the PBL, studies of
atmospheric aerosol properties and aerosol-atmosphere interactions, and airborne
measurements on the urban to mesoscale of atmospheric aerosol, turbulence
and radiation.
17.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.IDENTIFIERS/OPEN ENDED TERMS C. COSATI Field/Group
Mesoscale Prediction Models
Boundary Layer Modeling
Pollutant Removal Processes
Acdar, Acoustic Sounding
Airborne Measurements
18 DISTRIBUTION STATEMENT
Unlimited
19. SECURITY CLASS (This Reportj
Unclassified
21. NO. OF PAGES
20 SECURITY CLASS (This page)
Unclassified
22. PRICE
EPA Form 2220-1 (9-73)
-------
IV
INSTRUCTIONS
1. REPORT NUMBER
Insert the EPA report number as it appears on the cover of the publication.
2. LEAVE BLANK
3. RECIPIENTS ACCESSION NUMBER
Reserved for use by each report recipient.
4. TITLE AND SUBTITLE
Title should indicate clearly and briefly the subject coverage of the report, and be displayed prominently. Set subtitle, if used, in smaller
type or otherwise subordinate it to main title. When a report is prepared in more than one volume, repeat the primary title, add volume
number and include subtitle for the specific title.
5. REPORT DATE
Kach report shall carry a date indicating at least month and year. Indicate the basis on which it was selected (e.g., date of issue, date of
approval, date of preparation, etc.).
6. PERFORMING ORGANIZATION CODE
Leave blank.
7. AUTHOR(SI
Give name(s) in conventional order (John R. Doe, J. Robert Doe, etc.). List author's affiliation if it differs from the performing organi-
zation.
8. PERFORMING ORGANIZATION REPORT NUMBER
Insert if performing organization wishes to assign this number.
9. PERFORMING ORGANIZATION NAME AND ADDRESS
Give name, street, city, state, and ZIP code. List no more than two levels of an organizational hirearchy.
10. PROGRAM ELEMENT NUMBER
Use the program element number under which the report was prepared. Subordinate numbers may be included in parentheses.
11. CONTRACT/G RANT NUMBE R
Insert contract or grant number under which report was prepared.
12. SPONSORING AGENCY NAME AND ADDRESS
Include ZIP code.
13. TYPE OF REPORT AND PERIOD COVERED
Indicate interim final, etc., and if applicable, dates covered.
14. SPONSORING AGENCY CODE
Leave blank.
15. SUPPLEMENTARY NOTES
Enter information not included elsewhere but useful, such as: Prepared in cooperation with, Translation of, Presented at conference of,
To be published in, Supersedes, Supplements, etc.
16. ABSTRACT
Include a brief ^200 words or less) factual summary of the most significant information contained in the report. If the report contains a
significant bibliography or literature survey, mention it here.
17. KEY WORDS AND DOCUMENT ANALYSIS
(a) DESCRIPTORS - Select from the Thesaurus of Engineering and Scientific Terms the proper authorized terms that identify the major
concept of the research and are sufficiently specific and precise to be used as index entries for cataloging.
(b) IDENTIFIERS AND OPEN-ENDED TERMS - Use identifiers for project names, code names, equipment designators, etc. Use open-
ended terms written in descriptor form for those subjects for which no descriptor exists.
(c) COSATI FIELD GROUP - Field and group assignments are to be taken from the 1965 COSATI Subject Category List. Since the ma-
jority of documents are multidisciplinary in nature, the Primary Field/Group assignment(s) will be specific discipline, area of human
endeavor, or type of physical object. The application(s) will be cross-referenced with secondary Field/Group assignments that will follow
the primary posting(s).
18. DISTRIBUTION STATEMENT
Denote releasability to the public or limitation for reasons other than security for example "Release Unlimited." Cite any availability to
the public, with address and price.
19. &20. SECURITY CLASSIFICATION
DO NOT submit classified reports to the National Technical Information service.
21. NUMBER OF PAGES
Insert the total number of pages, including this one and unnumbered pages, but exclude distribution list, if any.
22. PRICE
Insert the price set by the National Technical Information Service or the Government Printing Office, if known.
EPA Form 2220-1 (9-73) (Reverse)
-------
ACKNOWLEDGEMENT
The Select Research Group gratefully acknowledges the
financial support? provided by research grant No. R800397 from
the office of Research and Development, Environmental Protection
Agency. The group also appreciates the financial support and
use of facilities from the Department of Meteorology and the
Center for Air Environment Studies of The Pennsylvania State
University.
An interdisciplinary research program such as the SRG effort
cannot possibly succeed without the contributions of the many
individuals who assisted on this project. The group wishes to
particularly thank the many faculty and staff members and graduate
students for their assistance.
-------
VI
TABLE OF CONTENTS
VOLUME I (of 2 Volumes)
Page
ACKNOWLEDGEMENTS v
CONTENTS OF VOLUME II ix
LIST OF FIGURES xii
I INTRODUCTION AND SCIENTIFIC OBJECTIVES 1
II THE DEVELOPMENT OF MESOSCALE MODELS SUITABLE FOR AIR
POLLUTION STUDIES : 6
ACKNOWLEDGEMENTS 7
1.0 INTRODUCTION 8
1.1 Potential use for regional and urban
dynamical prediction models 8
1.2 Some general considerations of the mesoscale
predictability problem 10
1.3 Overview of mesoscale modeling effort 14
2.0 THE REGIONAL MODEL 16
2.1 The basic equations in sigma coordinates for
a Lambert Conformal map projection 16
2.2 The horizontal and vertical grid structures 20
2.3 Finite difference equations 21
2.4 The two-dimensional analog 24
2.5 Kinetic energy budget equations for 2-D and
3-D models 25
2. 6 Lateral boundary conditions 28
2.6.1 Equations for mean motion over domain 29
2.6.2 Lateral boundary conditions for the 2-D
model 34
2.6.3 Lateral boundary conditions for the 3-D
model 36
2. 7 Initial conditions 37
2.8 Two-dimensional flow across the Appalachian
terrain 38
2.8.1 Specifications of the 2-D experiments 39
2.8.2 Results with geostrophic initial
conditions 43
2.8.3 Initialization of the boundary layer winds
considering the effects of surface friction 49
APPENDIX - CHAPTER 2.0 57
-------
Vll
3.0 PRELIMINARY THREE-DIMENSIONAL EXPERIMENT USING
REAL DATA WITH AND WITHOUT TERRAIN 62
3.1 Synoptic Discussion on 12Z Oct. 16 -
OOZ Oct. 17, 1973 63
3.2 Initialization and verification analyses and
specification of time-dependent boundary conditions. 67
3.3 Specification of parameters 75
3.4 Qualitative discussion of results 77
3.4.1 Low-level results 77
3.4.2 Middle level'results 82
3.4.3 Upper level results 84
3.5 Budget equations for the model domain and the
implications of the lateral boundary conditions 84
3.5.1 Mean kinetic energy budget for the 12-hour
forecast period 87
3.5.2 Time variation of the mean motion 89
4.0 NUMERICAL EXPERIMENTS WITH A TWO-DIMENSIONAL NESTED GRID. 97
4.1 The basic equations 98
4.2 The meshed grid system 99
4.3 Experimental results 103
4.3.1 Background experiment with uniform mesh 103
4.3.2 The treatment of the interface momentum
points in the meshed grid experiments 104
4.3.3 Meshed grid experiments with a mean wind
of 10 ms-1 110
4.3.4 Experiment with.mutually interacting grids
with the fine mesh moving through the
coarse mesh 112
4.4 Mesh grid experiments initialized with
Haurwitz waves 116
4.4.1 Initial conditions and the linear
solutions 117
4.4.2 Quantitative analysis of errors 123
4.4.3 Long-wave results: Exp. 1-5 125
4.4.4 Shortwave results: Exp. 6-8 140
4.4.5 Mixed long and short wave results 149
5.0 INVESTIGATION OF SEMI-IMPLICIT MODELS 158
5.1 Advantages of semi-implicit models over
explicit models 158
5.2 Comparison of one-dimensional explicit and
semi-implicit shallow fluid model 160
5.2.1 Development of explicit model 160
5.2.2 Development of semi-implicit model 161
5.2.3 Initialization of models 165
5.2.4 Results 165
-------
Vlll
5.3 Comparison of two-dimensional explicit and
semi-implicit models 1.TL
5.3.1 Development of the two-dimensional
S.I. model 171
5.3.2 Initialization 175
5.3.3 Results 176
5.4 Conclusions of preliminary tests of semi-
implicit models 186
APPENDIX - CHAPTER 5 187
6.0 DETERMINATION OF INITIAL DATA REQUIREMENTS 190
6.1 Development of Stochastic-Dynamic Equations-.. 192
6.2 Initialization Procedure 203
6.3 Energetics of the model 208
6.4 Interpretation of Predicted Variances. 210
6.5 Pure gravity wave experiments 211
6.6 A Monte Carlo comparison 219
6.7 Synoptic Scale Error Compatability 224
6.8 Summary and Plans for Future Research 229
REFERENCES 232
7.0 EXPERIMENTS WITH SIMPLIFIED SECOND-MOMENT
APPROXIMATIONS FOR USE IN REGIONAL SCALE MODELS 234
7.1 Introduction 234
7.2 The "Poor Man's Method" 236
7. 3 Semicomprehensive Methods 244
ACKNOWLEDGEMENTS 270
REFERENCES 271
-------
IX
TABLE OF CONTENTS
VOLUME II (of 2 Volumes)
Page
ACKNOWLEDGEMENTS v
CONTENTS OF VOLUME I ix
LIST OF FIGURES xii
III TASK 1-C BOUNDARY LAYER MODELING 272
1.0 SUMMARY OF PROGRESS 273
2.0 SHORT-TERM FORECASTS OF TEMPERATURE AND MIXING
HEIGHT ON SUNNY DAYS 275
3.0 ATMOSPHERIC TURBULENCE MODELING. 281
4.0 ATMOSPHERIC BOUNDARY LAYERS AND THE PRESSURE
GRADIENT-VELOCITY CORRELATION MODEL 284
5.0 NUMERICAL MODELING OF TURBULENCE FLOWS 289
6.0 MODELING TURBULENT FLUX OF PASSIVE SCALAR
QUANTITIES IN INHOMOGENEOUS FLOWS 293
7.0 EULERIAN AND LAGRANGIAN TIME MICROSCALES IN ISOTROPIC
TURBULENCE 304
8.0 NOTES ON TURBULENT FLOW IN TWO AND THREE DIMENSIONS 313
IV SRG ON AIR-POLLUTION METEOROLOGY
Part 1
1.0 ATMOSPHERIC EFFECTS ON PARTICIPATE POLLUTANTS 327
1.1 The sampling program 329
1.2 Chemical analysis of particulate matter 331
1.3 The numerical modeling program 334
1.4 Progress in field sampling of agglomeration 337
REFERENCES 341
Part 2
1.0 ATMOSPHERIC REMOVAL PROCESSES FOR AIR POLLUTANTS 349
1.1 Synopsis 350
1. 2 Personnel 351
1.3 Accomplishments 352
1.3.1 Global Emissions and Natural Processes
For Removal of Gaseous Pollutants 352
1.3.2 Rock 416
1.3.3 SO Solubility 419
1.3.4 Rate of S0? absorption by sea water 426
REFERENCES 431
-------
V OBSERVING SYSTEMS FOR URBAN AND REGIONAL ENVIRONMENTS 433
Part 1
Preface 434
1.0 TASK 1-D INTERPRETATION OF ACDAR SOUNDING OBSERVATIONS.. 437
1.1 Introduction 437
2.0 ANALYSIS OF DOPPER-SHIFTED MONOSTATIC ACDAR SIGNALS 441
2.1 System geometry 441
2.2 Surfaces of constant Doppler Shift 444
2.3 Horizontal wind 450
2.4 Antenna weighting 465
2.5 Total spectrum. 468
2.6 With antenna weighting 475
2.7 Wind shear included 483
3.0 MONTE CARLO METHOD FOR EVALUATING ACDAR
SCATTERING VOLUMES AND SYSTEM FUNCTIONS 503
4.0 MEASUREMENTS OF SOUND REFRACTIVELY TRANSMITTED
IN THE PLANETARY BOUNDARY LAYER 510
4.1 Introduction 510
4.2 System description 510
4.3 Measurement of inversion layer temperature
gradients 517
4.4 Fluctuations of signal levels and possible
association with atmospheric gravity wave
motions 520
REFERENCES 526
5.0 TEMPERATURE PROFILE MEASUREMENTS IN INVERSIONS
FROM REFRACTIVE TRANSMISSION OF SOUND 527
6.0 ANALYSIS AND SIMULATION OF PHASE-COHERENT
ACDAR SOUNDING MEASUREMENTS • • 547
Part 2
1.0 AIRBORNE MEASUREMENT SYSTEMS 592
1.1 Introductory Remarks 593
2.0 PSU ISOKINETIC INTAKE FOR AIRBORNE AIR SAMPLING 594
2.1 Design of the PSU Model II Probe 595
2.2 Sampler Performance 600
2.3 Conclusions 602
-------
XI
3.0 INSTRUMENTATION FOR MULTIWAVELENGTH AIRBORNE
PRECISION SPECTRAL RADIOMETER MEASUREMENTS .............. 604
3 . 1 The radiometers ................................... 604
3.2 Radiometer mounting ................................ 610
3.3 Radiometer signal conditioning ..................... 612
3.4 Correction of radiometer outputs ................... 619
3 . 5 Data analysis ...................................... 621
REFERENCES [[[ 622
4.0 AIRBORNE MEASUREMENTS OF TURBULENCE IN THE
PLANETARY BOUNDARY LAYER ................................ 623
5.0 AIRBORNE MEASUREMENTS OF AEROSOL IN THE ST. LOUIS
URBAN AREA .............................................. 633
6.0 SURFACE MEASUREMENTS OF AEROSOL IN A RURAL AREA
USING DIFFERENT METHODS .................... , ............ 641
6 . 1 Introduction ....................................... 641
6 . 2 Instrumentation ................................... 641
6 . 3 Summary of a selected data ......................... 643
7.0 TECHNIQUES FOR THE "MESOSCALE" INTERPRETATION OF
AIRCRAFT MEASUREMENTS ................................... 651
VI OTHER CONTRIBUTIONS .......................................... 659
Part 1
1.0 A GENERAL APPROACH TO DIFFUSION FROM CONTINUOUS
SOURCES ................................................
1 . 1 Theory
1.2 Analysis
Part 2
2.0 THE NIGHT-TIME MIXING DEPTH AT PHILADELPHIA ......... , . . 666
2 . 1 Introduction ...................................... 667
2 . 2 Analysis of observations ......................... 668
2.3 Results ........................................... 668
Part 3
3.0 SO CONCENTRATIONS AT KEYSTONE, PA ..................... 671
3.1 "The purpose of the project ........................ 672
-------
Xll
LIST OF FIGURES
VOLUME I
No. Title Page
Chapter I
I SRG Task Group Organization and Interactions 2
Chapter II _. part 2
2.1 30 x 50 grid on Lambert Conformal map projection.... 18
2.2 Error in geostrophic wind as a function of
horizontal distance for a 1°C temperature error
integrated over 200 mb depths centered at
300, 500, 700 and 900 mb 32
2.3 Magnitude of acceleration due to boundary flux
of momentum as function of wind speed and domain size (L)
The isotachs are labeled in ms~l 35
2.4 Temperature sounding for Washington, D. C.,
12Z Oct. 16, 1973 40
2.5 Time-averaged vertical cross section of u
for Exp. 2D-1 44
2.6 Time-averaged vertical cross section temperature
departure from Exp. 2D-1 45
2.7 Time-averaged vertical cross section of
vertical velocity (cm s~l) for Exp. 2D-1 46
2.8 Kinetic energy budget for Exp. 2D-1 48
2.9 East-west profiles of u-component in boundary
layer at various times in Exp. 2D-1 50
2.10 Hodograph of non-linear solution to simplified
boundary layer equations with geostrophic initial
conditions 53
2.11 Hodograph of non-linear solution of simplified
boundary layer equations with initial conditions
given by linear solution 55
2.12 East-west profiles of u-component in boundary
layer at various times in Exp. 2D-2 56
-------
Xlll
No. Title Page
CHAPTER II - Part 3
3.1 The 850 mb analysis for 12Z Oct. 16, 1973 65
3.2 The 850 mb analysis for OOZ Oct. 17, 1973 66
3.3 The 500 mb analysis for 12Z Oct. 16, 1973 68
3.4 The 500 mb analysis for OOZ Oct. 17, 1973 69
3.5 The 300 mb analysis for 12Z Oct. 16, 1973 70
3.6 The 300 mb analysis for OOZ Oct. 17, 1973 71
3.7 Observed u-component and geostrophic u-component
at level 2 (a = 0.325) at t = 0 73
3.8 Forecast 12-hour surface pressure change '(mb)
for Exp. 1 ending OOZ Oct. 17, 1973 78
3.9 Observed 12-hour surface pressure change (mb)
ending OOZ Oct. 17, 1973 79
3.10 Twelve-hour forecast velocities at level 6
(a = 0.925, p = 939 mb) from Exp. 1. Contours
of terrain are labeled in m 80
3.11 Twelve-hour forecast velocities at level 6
(a = 0.925, p = 929 mb) from Exp. 1A,
no terrain 80
3.12 Verficiation velocities at level 6 (a = 0.925,
p = 929 mb) from interpolation of synoptic
analysis to mesoscale grid 81
3.13 Twelve-hour forecast velocities at level 4
(a = 0.625, p = 709 mb) from Exp. 1 83
3.14 Twelve-hour forecast velocities at level 4
(a = 0.625, p = 709 mb) from Exp. 1, no
terrain 83
3.15 Twelve-hour forecast velocities at le\rel 2
(0 = 0.325, p - 489 mb) from Exp. 1 85
3.16 Twelve hour forecast velocities at level 2
(a = 0.325, p = 489 mb) 85
-------
XIV
No. Title Page
3.17 Twelve-hour forecast temperature change at level 2
(p - 489 mb) from Exp. 1A, no terrain 86
3.18 Twelve-hour forecast temperature change at level 2
(p - 489 mb) from Exp. 1, with terrain 86
3.19 Components of kinetic energet budget from
Exp. 1 and 1A 88
3.20 Solution of mean vector motion over domain for
constant B and V 91
"" *v R
3.21 Components B /f and B /f (mean non-linear accelera-
tions scaledUby f) for Exp. 1A 93
3-22 Hodograph of mean wind (V), mean geostrophic wind
(V ), and -iB/f for various times in Exp. la
atglevel 1 (p = 342 mb) 94
3.23 V at t = 0 and t = 12 hours and hodograph of mean
velocity, V for 12-hour forecast, Exp. 1A 95
CHAPTER II - Part 4
4.1 Boxes representing incremental volumes of mass in
the 4:1 meshed grid system 100
4.2 Boxes representing incremental volumes of momentum
in the 4:1 meshed grid system 100
4.3 Time variation of total energy in uniform CMC
experiment (Exp. 1) 105
4.4 Time variation of total energy in Exp. 2-7 107
4.5 Time variation of total energy and maximum
v-component in Exp. 8, 9, 10 Ill
4.6 Illustration of movement of FMG within CMC 114
4.7 Time variation of total energy and maximum north-
south velocity component in Exp. 11 115
4.8 Initial height and vorticity field on CMC for
Exp. 1-5 127
4.9 Initial u- and v-component fields on CMC for
Exp. 1-5 128
-------
XV
No. Title Page
4.10 Initial height and vorticity field on
CMC f or Exp . 1-5 ................................... 129
4.11 Initial u- and v-component fields on CMC
for Exp. 1-5 ....................................... 130
4.12 Forecast height and vorticity fields on CMC
for Exp. 1 ......................................... 131
4.13 Forecast u- and v-component fields on CMC
for Exp . 1 ......................................... 132
4.14 Forecast relative vorticity and height field in
Exp. 2 (1-way no diffusion) ........................ 134
4.15 Forecast relative vorticity and height field in
Exp. 3 (2-way no diffusion) ........................ 135
4.16 Forecast relative vorticity and height field in
Exp . 4 (1-way with diffusion) ...................... 136
4.17 Forecast relative vorticity and height field in
Exp. 5 (2-way with diffusion) ...................... 137
4.18 Initial height and v-component field on CMC for
for Exp . 6-8 .......................................
4.19 Initial height and v-component field on FMG
for Exp . 6-8 ....................................... 142
4.20 Forecast v-component and height field on FMG at
6 hours for Exp. 7 (1-way with diffusion) .......... 144
4.21 Forecast v-component and height field on FMG at
6 hours for Exp. 8 (2-way with diffusion) .......... 145
4.22 Forecast v-component and height field on FMG at
12 hours for Exp. 7 (1-way with diffusion) ......... 146
4.23 Forecast v-component and height field on FMG at
12 hours for Exp. 8 (2-way with diffusion) ......... 147
4.24 Ratio of variance to MSB for Exp. 6-8 .............. 148
4.25 Initial <(> for Experiments 9-11 (mixed long and
short wave) ........................................ 150
4.26 Initial v-component for Exp. 9-11 (mixed long
and short waves) ................................... 152
-------
XVI
No. Title Page
4.27 Forecast v-component at 13.3 h from Exp. 9
(uniform FMG) 153
4.28 Forecast v-component at 13.3 h from Exp. 10
(1-way meshed system) 154
4.29 Forecast v-component at 13.3 h from Exp. 11
(2-way meshed system) 155
CHAPTER II - Part 5
5.1 Comparison of one dimensional explicit and S.I.
forecasts at 24 hours 167
5.2 Number of iterations required per time step as a
function of a for four different first guesses 170
5.3 Two dimensional staggered grid 174
5.4 Twenty-five hour forecast of h by explicit model... 178
5.5 "Error" height fields for the explicit model at
25 hours 181
5.6 "Error" height fields for the S.I. model at
25 hours 182
5.7 Root-mean square differences between explicit
and S.I. model height forecasts and the linear
solution 183
5.8 Number of iterations required per time step as a
function of a for 2 different first guesses 185
CHAPTER II - Part 6
6.1 Spatial and temporal variation of au, a^, and h for
the pure gravity wave integration. Graph (a)
pertains to 0 minutes; (b) 10 minutes; (c) 20
minutes; and (d) 30 minutes 213
6.2 Spatial and temporal variation of the standard
deviations of u and h during the first five
minutes of the pure gravity wave integration.
Curves are labeled in minutes , 214
-------
XVI1
No. Title
6.3 Time variation of the certain energy modes.
TCE is the total certain energy, CPE is the
certain potential energy and CKE is the certain
kinetic energy 216
6.4 Time variation of the uncertain energy modes.
TUE is the total uncertain energy, UPE is the
uncertain potential energy and UKE is the un-
certain kinetic energy. TCE, the total
certain energy, is included for comparison
with the total uncertain energy 217
6.5 Time variation of the total energy modes. TPE
is the total potential energy, TKE is the total
kinetic energy and TE is the total energy 218
6.6 Monte Carlo solution with 500 trials (dashed line)
and the stochastic-dynamic solution after 20 minutes
of integration for conditions of pure gravity wave
motion 221
6.7 Monte Carlo solution with 500 trials (dashed line)
and the stochastic-dynamic solution after 30 minutes
of integration for conditions of pure gravity
wave motion 222
6.8 Transition of the uncertain kinetic energy of the
v velocity component for various initial, standard
errors of the v component. The curves are labeled
in terms of the magnitude of the standard deviation
of the initial v error, in m s~l. The initial
standard deviation of the equivalent temperature
error was 1°C 228
6.9 Transition of the standard deviation of v in terms
of its spatial average for various initial, standard
errors of the v component. The curves are labeled
in terms of the magnitude of the standard deviation
of the initial v error, in m s 1. The initial
standard deviation of the equivalent temperature
error was 1°C 230
-------
xvin
No. Title Pag<
'CHAPTER II - Part 7
1. Buoyancy Compensation 248
2. Comparison of predicted relation between a /u
.w
*
and z/L with observations from the AFCRL Kansas
1968 Field Program .................................. 263
3. Comparison of predicted relation between l/cf>
(unstable side) or $ (stable side) and z/L with
observations from the AFRCL Kansas 1968
Field Program ....................................... 264
4. Comparison of predicted relation between l/, (stable side) and z/L with
observations from the AFCRL Kansas 1968
Field Program ....................................... 265
5. Comparison of predicted relation between
l^o/T^I and z/L with observations from the
AFCRL Kansas 1968 Field Program ..................... 266
6. Comparison of predicted relation between
R./(z/L) and z/L with observations from the AFCRL
Kansas 1968 Field Program ........................... 267
7. Comparison of predicted relation between
-u9'/w9' and z/L with observations from the
AFCRL Kansas 1968 Field Program ..................... 268
-------
I INTRODUCTION AND SCIENTIFIC OBJECTIVES
The "Select Research Group in Air Pollution Meteorology"
(SRG) sponsored by the U. S. Environmental Protection Agency
(EPA) under the auspices of the Meteorology Laboratory was
established to identify and research long range, inter-
disciplinary problems in air pollution meteorology. Following
national competition, the proposal by The Pennsylvania State
University for a five year program was selected and first funded
in May, 1972.
This report is a progress report from the SRG which
summarizes, in particular, research activities it conducted
and work it completed during the Ilnd year of the grant program.
The primary objective of the SRG program is to develop a
comprehensive air pollution model which can be used to examine
fundamental problems which require solutions in order for
operational air pollution prediction models to be constructed
and used. Clearly, a control agency such as the EPA must have the
capability of quantitatively relating air pollution sources to
their subsequent distribution in an urban or regional environment.
The development, testing, and evaluation of a regional model
by the SRG (Task Group 1A) is, of necessity, an intra and inter-
disciplinary enterprise (see Fig. I). The design and coding of the
basic numerical model requires, firstly, expertise in numerical
weather prediction and dynamic and synoptic meteorology.
-------
o
H
o
w
H
W
a
PL,
W
H O
cd M
o o
•H
60 CO
O 4->
•H a
O Q)
^ 6
O QJ
0) S-l
4J 3
0) CO
S cd
(0
T)
ex
o
c
o
•H
4-1
0)
•H
TJ
I
60 /
£3
s
eu
o
o
I
a)
n
3
CO
C
O
CU
e a)
t-l 4-1
o a)
-rl O
60 (3 H
O 0) O
H iH CO W
O 3 O Cd
H ,0 )^ -H
O ^ (U T3
• • •
CO
a
a
o
Q) O
36
H CO
O
-------
Since many of the important physical processes which determine
meteorological conditions on the mesoscale occur within the
atmosphere's planetary boundary layer, it is essential to have
scientists specializing in atmospheric turbulence and boundary
layer processes (Task Groups IB and 1C) actively contributing to
the modeling effort.
Some of our most powerful continuous measurement techniques for
space-time variations of the meteorological structure in the lower
atmosphere are based on indirect atmospheric probing systems such
as acdar and lidar. Thus, it is also essential to include an in-
direct probing effort (Task ID) in order to provide specialized
observations and advice concerning data assimilation techniques
to the modeling groups.
If air pollutants were simply inert particles or gases carried
about by the wind, it would be possible to develop a straight-forward
dynamical model to predict concentrations of pollutants as a function
of time and location. However, most pollutants are in a state of
constant modification by the atmosphere — through processes such as
photochemical reactions and washout. It is, thus, necessary for the
dynamic modeling group to also interact on a continuing basis with
atmospheric chemists (Tasks 2 and 3). Cross fertilization between
these groups better enables the modelers to develop optimal procedures
-------
for treating those variables which determine, for example,
atmospheric aerosol properties and, also, provides the chemists
with a better, and necessary, appreciation for the characteristics
of the atmospheric "laboratory" environment.
Finally, an instrumented aircraft is, presently, the most
practical and economical sensor platform available capable of
providing a large variety of meteorological and air pollution-
related observations on the mesoscale. Direct "synoptic-type"
observations on a regional scale are both prohibitively expensive
and difficult to effectively utilize. An aircraft can transport
sensors over mesoscale distances so that both horizontal and vertical
profiles of meteorological variables, and/or air pollution parameters,
may be obtained. The airborne measurements group (Task IV), thus, has
the responsibility of providing measurements support to each of the
others. For these purposes, specialized meteorological, turbulence,
aerosol and radiation airborne measurement capabilities have been
developed for use both at PSU and by scientists in the EPA Meteorology
Laboratory.
The principal goal of the Select Research Group is one that guides
nearly every large research group in Air Pollution Meteorology. It is
widely recognized that the ultimate achievement of this goal will not
be reached in any limited period. The objective that will be attained
-------
within the SRG project period is significant progress toward the
development of a regional-scale air hydrodynamic model suitable
for studying mesoscale perturbations to large-scale flows that
are important in the transport of air pollutants. The model will be
applicable to at least a limited number of practical situations
and may be used in conjunction with Lagrangian Puff models to
simulate the modification of pollutant concentrations. Also, the
forecast wind and temperature distributions and the height of the
mixed layer may be used in other more simplified air pollution
models. Of course, the ultimate perfection of air pollution
modelling will require efforts and time far beyond the lead-off
research presently being performed by the SRG scientists.
-------
II THE DEVELOPMENT OF MESOSCALE MODELS SUITABLE
FOR AIR POLLUTION STUDIES
Richard A. Anthes
Nelson Seaman
Joseph Sobel
Thomas T. Warner
-------
ACKNOWLEDGEMENTS
Mr. James Hoke contributed to the solutions of the linear and
non-linear boundary layer equations in Section 2.8.2. Mr. Ed O'Lenic
did the synoptic analyses of Part 3. The three-dimensional model
experiments were run on the NCAR computer. NCAR is supported by
the National Science Foundation.
-------
THE DEVELOPMENT OF MESOSCALE MODELS SUITABLE
FOR AIR POLLUTION STUDIES
Richard A. Anthes, Nelson Seaman,
Joseph Sobel, Thomas T. Warner
1.0 INTRODUCTION
The Pennsylvania State University (PSU) and EPA are in the early
stages of developing time-dependent dynamic models for use in air pollution
studies. PSU is developing a general, hydrostatic model suitable for fore-
casting flows with characteristic horizontal wavelengths of 100-600 km under
a variety of synoptic conditions. EPA is working on a general, non-
hydrostatic model for forecasting smaller scale circulations (horizontal
wavelengths 5 to 30 km). The EPA model is designed for prediction of
flow over and in cities (in particular, St. Louis in support of the RAPS
program) while PSU is modeling perturbations to the synoptic flow induced
by terrain variations, land-sea and land-lake contrasts, etc. The PSU
model is termed a regional model.
1.1 Potential use for regional and urban dynamical prediction models.
The primary goal of EPA in developing the above model lies in utilizing
the models as research tools capable of providing meteorological information
that is necessary for rational decision making over questions concerning air
standards. Specifically, information such as expected flow patterns and
-------
concentrations of pollutants under a variety of weather conditions at
diverse locations is required in the preparation of environmental impact
statements that are now required before proposed actions such as nuclear
power plant construction are approved.
Both the urban and regional models will predict flows for time
periods of 1 to 24 hours, and hence, will be directly suitable for
estimating maximum and average pollutant concentrations over a few hours.
The models will not be directly capable of predicting annual averages, but
stratification of short-range model forecasts according to synoptic
conditions, together with frequency distributions of synoptic regimes, may
make annual estimates possible.
Both the regional and urban models are basically hydrodynamic models
rather than direct pollutant prediction models. Presumably, if the
dynamics are modeled correctly, so that the three-dimensional wind,
temperature, and moisture structure are known, Lagrangian dispersion
models may be utilized with the help of the output from the dynamic models
to predict concentrations associated with individual plumes. Of course,
if the plume size becomes larger than roughly six times the grid size
of the model, explicit prediction of the plume's subsequent behavior is
possible. This latter event is likely to occur often in the urban model,
and may occur occasionally in the regional model under regional air pol-
lution episodes caused by emission of pollutants from multiple sources
over a long period of unfavorable weather conditions.
Although neither the urban nor the regional model is contemplated for
operational use on a daily basis, it is possible that data from realistic
3-D models could be used in lieu of expensive observational networks to
-------
10
provide dense, accurate data for use in statistical air pollution models.
It may also be noted that the interest in and potential benefits of
general three-dimensional models are not restricted to air pollution
studies, but include resolution of regional and urban variability in
cloud cover, rain and snowfall, temperatures, etc. that are important
for effective land use. Finally, the 3-D models are powerful tools for
increasing our understanding of the smaller scales of atmospheric motion.
Eventually, it may be desirable to allow the urban and regional
models to interact. The most obvious mode of interaction is to use the
output from the regional model to provide time-dependent boundary condi-
tions for the urban model. Assuming that the important transfers of
energy are downscale, the regional model could be run for 12 hours in-
dependently of the urban model. The mean values of wind, temperature and
moisture in the vicinity of the city could be saved, and used on the
boundaries of the fine mesh urban model.
1.2 Some general considerations of the mesoscale predictability
problem
A basic problem in the development of both the urban and regional
models is the lack of knowledge concerning the predictability of flows
on these scales. Unlike large-scale, quasi-geostrophic flow, no scaling
of the basic hydrodynamic and thermodynamic equations will produce a
simplified set of equations that will apply under most meteorological
conditions. For example, the release of latent heat by condensation,
a secondary effect for planetary flows over short (1-2 day) time periods,
-------
11
varies from being totally negligible under many conditions to dominant
under others. Similarly, non-hydrostatic effects may not be ignored
under all situations.
The general prediction problem may be discussed by consideration
the equation
F(a,x,t) (1.1)
where a is the variable to be forecast (e.g., wind component, temperature,
humidity, pollutant concentration), a is a climatological or large-scale
->
mean, V is the three-dimensional vector velocity, x is the three-
dimensional position vector, and F is some (complicated) function of
time and space. The first term represents the redistribution of initial
perturbations by the wind field, while the second term represents mod-
ifications of the initial perturbation by local forcing. Integrating in
time, we find
(a-a) (x,t) = (a-a) (x,t ) - V • V (a-a) dt
~ o J
t
o
(1-2)
ft
+ F(a,x,t) dt
J ~
o
that is, the perturbation to be forecast is determined by the initial
distribution of the perturbation, the redistribution of the perturbation by
the wind, and the production of additional perturbations by the local
forcing.
-------
12
The predictability of various scales of flow under different
synoptic conditions can now be discussed in terms of (1.2). In general,
a prognostic model must include:
1. knowledge of the initial conditions,
2. specification of boundary conditions on the horizontal and
vertical boundaries, and
3. representation of the local forcing function by either explicit
or parameterized methods.
Accurate treatment of all three aspects of the problem for the urban and
regional models presents a formidable challenge. It is doubtful whether
simultaneous initial observations will ever be possible on a 30 x 30 x 6
computational mesh. Neither is it likely that a dense network of time-
dependent boundary conditions will ever be realizable, because larger-
scale models are only able to provide mean values along the lateral
boundaries. Finally, the local forcing function will, under various
circumstances, include:
1. complex terrain effects,
2. condensation heating and evaporative cooling,
3. infrared and short wave radiation,
4. conduction of heat at the earth-air interface and turbulent
transfer of heat upward in the boundary layer.
While a general dynamic model demands accurate treatment of the
initial conditions, boundary conditions, and forcing function, under
particular situations the problem may be simplified so that one or more
of the above effects dominates and the other (s) may be treated fairly
simply. For example, the evolution of large-scale flows (wavelengths
2,000-10,000 km) over 1-2 days is determined mainly by the conservation of
-------
13
potential vorticity. Local forcing is weak, as evidenced by the success
of barotropic models. Lateral boundary conditions are eliminated by con-
sidering the entire global domain. Here the specification of the initial
conditions is of supreme importance in the prediction of synoptic-scale
flows. We might then ask whether a comparable set of initial data is
necessary for the regional and urban models. Fortunately, we suspect that
they are not. As the horizontal wavelength of atmospheric disturbances
decreases, the solutions are determined more by local forcing than by
initial conditions, i.e., the atmosphere "forgets" its initial state more
quickly and adjusts to the regional and urban scale forcing. Preliminary
numerical experiments with the regional model indicate that large pertur-
bations to the mean flow develop within about 6 hours as a consequence of
mesoscale terrain variations. Certainly, the adjustment time of flows over
cities is much less, so we hypothesize that in many synoptic situations,
if_ _lpcal forcing is modeled correctly^ the details of the initial perturba-
tions are not particularly important.
The relative importance of the lateral boundary conditions varies
tremendously under differing synoptic flows. For stagnant air situations,
the treatment at the boundaries is relatively straightforward. However,
when even slight flow occurs across the boundaries, it becomes important
to specify at least the mean values of a correctly. Larger scale models
may be used to provide lateral boundary conditions for the regional model,
which in turn could provide boundary conditions to the urban model.
Exceptions to this hypothesis occur when temperature or moisture dis-
continuities such as those associated with fronts or "dry lines" exist.
In these situations the initial location of the discontinuity will be
important.
-------
14
The one aspect of the mesoscale prediction problem that will nearly
always be important is the numerical treatment of the local forcing
function (if the local forcing were unimportant, a mesoscale model
would be unnecessary). Included in the representation of the local
forcing is the parameterization of important physical processes.
Parameterization is the inclusion of the effect of some physical process
in the model in terms of one or more simplified parameters. For example,
the representation of the effects of surface friction by the quadratic
stress law and a drag coefficient is a parameterization. At first, we are
planning to parameterize:
1. the long- and short-wave radiation effects,
2. conduction of heat at the ground, and
3. the transfer of heat and momentum by subgrid-scale motions (the
eddy fluxes).
Eventually, we may want to parameterize condensation and evaporation; however,
since many dynamic and air pollution problems are not affected by phase
changes of water, we will neglect these processes for the time being.
1.3 Overview of mesoscale modeling effort
The next sections summarize the research supported by EPA under
Tasks 1A and IB since the First Annual Project Report (June 1973). First,
the numerical aspects of the three-dimensional regional model and its
two-dimensional analog are presented. Results from several preliminary
experiments with the 3-D and 2-D models are presented. These include a
first attempt at forecasting with real data, obtained from the field ex-
periment in the Middle Atlantic States in October 1973.
-------
15
An important part of the mesoscale modeling problem is the meshing
of the regional model with large-scale models, and perhaps eventually with
small-scale models. Section 4 presents results from one-, two-, and
three-dimensional meshed grid experiments.
Because mesoscale models, meshed or otherwise, are voracious
consumers of computer time, we are investigating the feasibility of semi-
implicit rather than explicit modeling techniques. In Section 5, results
from one- and two-dimensional semi-implicit models are compared with their
explicit counterparts. Indications are that an increase in computational
efficiency by a factor of four may be possible with the semi-implicit
technique.
Finally, a major effort is being directed toward the theoretical investi-
gation of the initialization requirements and the predictability of the meso-
scale flows. Through the use of stochastic-dynamic models, insight into the
initial data requirements in terms of spatial density and error requirements
may be made. The growth of uncertain energy provides estimates for the time
scale of the predictability of the mesoscale.
-------
16
2.0 THE REGIONAL MODEL
The numerical details of the three-dimensional (3-D) model are
described in this section. First the basic equations are presented in
generalized a-coordinates written for Lambert Conformal map coordinates.
Then the structure of the grid, the finite difference equations, and the
lateral boundary conditions are described. Finally, some results from
the 3-D model and its two-dimensional analog are presented.
2.1 The basic equations in sigma coordinates for a Lambert Conformal
map projection
For middle latitude regions, Lambert Conformal map coordinates
produce less distortion of the earth's geometry than either polar
stereographic or Mercator map coordinates (Saucier, 1955). Because the
authors were not able to locate in the literature the derivation of the
equations of motion, thermodynamic equation, and continuity equation
for the Lambert Conformal projection, an abbreviated derivation is
presented in Appendix 2.1. The vertical coordinate is a, defined by
P - P^.
where p is pressure, p is the pressure at the top of the model (assumed
constant) and p is the surface pressure. This coordinate system follows
s
the terrain, and allows the top of the model to be placed at any pressure
surface.
-------
17
The complete equations of motion in a-coordinates, written for
Lambert Conformal map coordinates in flux form are
8p*u _ _ 2 ,5p*u u/m 3p*u v/nu _ 8p* ug
3t " m { 3x 3y ; 8a
-mp*
RT
w (2 cos * y - «> +
IT cl
u v/m 8p*v v/mN
3t ~ " 9x + 3y ) ~ 9a
RT
cos
where x and y are the Lambert Conformal map coordinates (Fig. 2.1), u and
v are the speeds relative to the earth in the x and y directions, w is
the vertical velocity, p* = p (x,y)-p , p is surface pressure, R is the
s t s
gas constant for dry air, T is temperature, f is the Coriolis parameter,
a is the radius of the earth and FU and FV are the frictional forces (to
be discussed later). The quantities m, r, and y are defined by
rtan
[
tan
(sin (b-n) ,-uy . vx.
Y = ~ •*•—— ( J + —)
Y a cos (j) v r r ;
-------
Figure 2.1 30
Lambert Conformal map projection
-------
19
rtanJV2_ n ( .
Ltan i|; /2J (Z'b)
where n = 0.716, and ip = 30°, is latitude and ty is colatitude.
For these choices of n and i/J , the projection is true at 30°N and
60°N. Note that a Mercator projection, true at the equator and
<{>, = 90-1^.., is obtained by setting n = 0, and the polar stereographic,
projection is obtained by setting n to equal 1.
For consistency in the kinetic energy equation when the hydro-
static equation is used instead of the vertical component of the
equation of motion, the terms multiplied by w are neglected in (2.2) and
(2.3). Furthermore, y is at least two orders of magnitude smaller than
f in middle latitudes, and may be neglected. The equations of motion
utilized in the model are therefore given by (2.2) and (2.3) without the
underlined terms.
The complete continuity equation in Lambert Conformal, cr-coordinates
is
/0 ...
(2-7)
* 2 r8u/m , 9v/m-, 9p*a -2wp*
=-m [ + ] • ~ ~~
The last term in (2.7) is small and is neglected in order to conserve
mass. The thermodynamic equation is
9p*T 2 ,--3up*T/m 3vp*T/m, _ 9p*Tq
(2.8)
cpa cp
FT
-------
20
where Q is the diabatic heating per unit mass, c is the specific heat
at constant pressure, and co = dp/dt is related to a by
= p*a + a . t (2.9)
where
The term FT represents the lateral diffusion of heat by subgrid-scale
eddies. In the a-system, the hydrostatic equation may be written
P*RT
P*a
2.2 The horizontal and vertical grid structures
The staggered horizontal grid is identical to the SI grid
described by Anthes (1972) and Anthes and Warner (1974), and is not
illustrated here. The vertical grid is also staggered, with the
velocity components and temperature defined at one set of levels while
•
geopotential height and vertical velocity (a) are defined at adjacent
levels. The horizontal and vertical grid spacings and the domain size
are variable. The model has been tested with 30 x 30 and 30 x 50 grids
of 20 km resolution and a 30 x 30 grid with 5 km resolution.
-------
21
2.3 Finite difference equations
The finite difference equations for the SI grid conserve mass,
momentum, and total energy.
We first define the following finite-difference operators:
_ ai. j + 1/2 ai, j - 1/2
-- - -- — -
x
Ax
a = > - ' (2.12)
y Ay
1/2
and
a = ^ai + 1/2, j " Ui - 1/2, j
where j is the east-west index and i is the north-south index.
We also introduce the following four-point operators:
1. .1 + 2ai. .1 + "i - 1, .1
and (2>13)
a. . . , + 2a. . + a. . _
= i, J + 1 i. J i, J ~ 1
-------
22
For vertical differences and averages, we define
-0 (0tk + 1/2 + \ -
and (2.14)
E (ak + 1/2
The finite-difference equations associated with the staggered grid
for the u- and v-component equations of motion are
' —y —xy —a -xy -—y
9p*u ~ 2r/-x p*u , , ,-y p*v XN , 6a p*u . ,
C ~ - m [(u £— ) + (u-7 £— ) ] T—c mp* 4>
9t L m x m yj oa r x
- m
(2.15)
+p*[(KHu)y]
y y v
and
3p*v ~ 2 r ,—x p*u s ,—y p*v v <5a p*v
£. - - m [(v •*-— ) + (v7 •*-— ) 1 j—*•—
3t tv m 7x m xyj 60
—xy —a —xy —x
- mp*
-------
23
where the finite-difference analogs for the vertical friction terms, FU
and FV are discussed later. In this section the terms (p*u) and (p*v)
represent (up**^) and (vp**), respectively.
The analogs for the continuity and thermodynamic equations are
- *
-* <2-i?>
and
(2.18)
S?L
+ Q + p* [V T ] + p* [K T ]
c a c H xx v H y y
P .P
^-^
The notation T in equation (2.18) signifies that potential tem-
perature, rather than temperature itself, is linearly interpolated
between CT-levels. This provides for a substantial improvement in the
calculation of static stability since potential temperature is much
more nearly linear in a than is temperature. Mintz and Arakawa (see
Langlois and Kwok 1969), have adopted a similar procedure for their
general circulation models. The finite difference expression for
oj is
dt ** " ' "V3t ' ~ * x ' v P y
-------
24
In equations (2.15) and (2.16) the use of the four-point averaging
operator is necessary if the finite-difference equations are to conserve
kinetic energy (Anthes, 1972). For the approximation of the time derivatives,
the centered-time or leapfrog scheme, which has second order accuracy, is
adopted.
2.4 The two-dimensional analog
Because of the complexity and relatively high cost of the three-
dimensional model, it is frequently useful to test numerical aspects
of the model (such as vertical resolution) or the parameterization of
physical processes (such as turbulent transfers of heat and momentum in
the planetary boundary layer) with a two-dimensional analog to the 3-D
model. The 2-D model is designed to be as much like the 3-D in structure
as possible, so that the many properties of the model that are not inherently
three-dimensional may be investigated with this more economical model. The
2-D analog is a vertical, east-west cross section through the atmosphere
in which the y-variation of the forecast quantities are neglected. For
simplicity, the map scale factors are considered to be unity. With these
assumptions, the equations of motion, continuity equation, and thermodynamic
equation for this model corresponding to (2.2), (2.3), (2.7), and (2.8) are
3p*u _ _ 3p*uu _ 3p*u _ RT
3t ~ ~ 3x 3cr (1 + p../p*a) 9x ~ p
(2.20)
+ fp*v + FU
-------
25
9p*v 9p_*uv 9p*vg c, N , T,T, /o TIN
-§— = § t^- f (u - u ) + FV (2.21)
9t ox da g
dp* _ 9p*u 9p*a ,„ .
9t ~ ~ 3x ~ 9a (2>22)
and
9t 9x 9a c a c
P P
where u is the specified geostrophic wind.
o
The equation for to is given by (2.9) and (2.10) after neglecting
(v TT*— ) , and the hydrostatic equation is given by (2.11).
dy
The vertical grid structure is identical to that of the 3-D
model. The horizontal grid is staggered with u and v defined at points
x. = (j-l)Ax and all other variables defined at points x. = (j - 1/2) Ax,
J J
2.5 Kinetic energy budget equations for 2-D and 3-D models
The kinetic energy equation, formed by the multiplication of (2.2)
by u and (2.3) by v and the use of the continuity equation, is
-------
26
+ V • V p* k + = _ p* v . vcj,
~ d° " (2.24)
- RT V • Vp* + u • FU + v • FV ,
2 2
where k = (u + v )/2.
For diagnostic purposes, it is convenient to integrate the kinetic
energy equations over the mesoscale domain in order to obtain the equation
for the time rate of change of the mean kinetic energy,
3p*k 1 I rfyN
3t g K j - (uP*k)E + W dy
a x y a yg
rx *
+ E - (vp*k)N + (vp*k)g dx - 9gP*k - p*v • V<|) (2.25)
LJ .
- RT V • Vp*^ + (u • FU + v • FV)] da
where the subscripts N, E, S, and W denote the northern, eastern,
southern, and western boundaries, respectively, and the ( ) operator
denotes an average over the horizontal domain.
-------
27
For the first two terms on the right side of (2.25) represent the
change in mean kinetic energy in the volume by the net flux of kinetic
energy across the east-west and north-south lateral boundaries respectively.
Unlike their relatively insignificant role in large-scale models, these
terms may be very important in mesoscale models, contributing as much
or more as the linear terms. This fact emphasizes the need for accurate
specification of the velocities on the lateral boundaries. The third term
on the right side of (2.25) is the vertical flux divergence of kinetic
energy, and vanishes if the entire depth of the model is considered,
•
since a = 0 at the upper and lower boundaries. The next two terms represent
the conversion of potential to kinetic energy by cross-isobar flow. Over
variable terrain the terms are individually very large, but are of
opposite sign and tend to cancel with each other. The last term in
(2.25) represents the loss of kinetic energy by frictional dissipation.
While the interpretation of the total kinetic energy over the
3-D mesoscale model domain is straightforward, a total kinetic energy
cannot be assigned to the 2-D domain, which represents a slice through
the atmosphere. Instead, the kinetic energy per unit area, K, is defined
fXF
K = £- ] I kp* dx da = g " | kp* da , (2.26)
ax,, a
where L is the length of the 2-D domain and equal to (x^-x^). Similarly,
the time rate of change of kinetic energy per unit area in the 2-D model
is
-------
28
3K i r (up*k)E - (up*k)w
3t g
a
- u(p* + .. . . . ) + fv u P* (2.27)
3x (1 H- pt/p*a) 9x g
+ (u • FU + v • FV)
where the operator ( ) indicates an average along the entire length of
the domain.
2.6 Lateral boundary conditions
The numerical treatment of the lateral boundaries is a difficult but
very important aspect of a limited area forecast model. The problem on the
interior of the domain is well posed only when the proper set of boundary
conditions is prescribed. The solution on the interior may be completely
different for apparently minor variations in the boundary conditions. Much
discussion has appeared in the literature concerning the proper specifica-
tion of variables on the open boundaries (Moretti, 1969; Shapiro, 1970). As
discussed by Moretti, the only straightforward case occurs when the flow
is supersonic, in which case the value on an inflow boundary may be
arbitrarily specified, and the values on an outflow boundary obtained by a
linear extrapolation outward from the interior of the domain.
In the subsonic case, however, waves may propagate upstream, so
that values at inflow boundaries are partially dependent upon the
flow on the interior of the domain. Thus the arbitrary specification
of either constant or time-dependent boundary values may make the problem
ill-posed. However, in the absence of proven computational methods to
-------
29
compute boundary variables correctly under general conditions, we are
forced to take the pragmatic view that if the lateral boundaries are
located far enough away from the region of interest, the errors introduced
at the boundaries will remain within some acceptable tolerance in the
interior of the domain during the forecast period. We must therefore
search for a set of conditions which minimize the errors generated by
the boundaries and their feedback into the interior.
Although it is normally recognized that a set of boundary conditions
that minimizes the generation of high temporal and spatial frequency
components to the numerical solution is desirable, if not absolutely
necessary, it is less generally recognized that smoothly varying
solutions near the boundary are not sufficient to guarantee accurate
solutions on the interior of the grid. It is apparent that as the size
of the horizontal domain decreases, the specification of the velocity
components and temperature along the boundaries affects the mean (wave
number zero) values of these quantities over the entire domain to an
ever increasing degree. Thus on a mesoscale domain of 600 x 600 km, a
set of boundary conditions may be computationally "stable" and produce
"smooth" results, but even small errors in the treatment of temperature
or velocity may profoundly affect the mean kinetic and internal energy
budgets over the domain.
2.6.1 Equations for mean motion over domain. The importance of
specifying accurate values of temperature, surface pressure, and velocity
components on the lateral boundaries may be shown by deriving the equations
for the mean motion over the domain. For convenience, we will neglect
vertical velocities, friction, and the covariance of p* with the velocity
components, and integrate (2.2) and (2.3) over area, obtaining
-------
30
and
where
(u
-fu
v)E - (uv;
L
X
' + *r°
lj
"Ly ~ 0
i / 2 2.
) ( v — v )
W V N S '
L
y
3t
x y
(2.28)
Ly - Lx
(UVN) " (UVS}
3t L L
x y
(2.29)
, --D c — D
+ f v - fv
f = _ RT
U
_ _
g (p* + Pt/a) 3y 3y
(2.30)
f v =
g (p* + p/a) 3x 3x
and L and L are the lengths of the domains in the east-west and
x y
norht-south directions, respectively.
For simplicity, let us consider only the u-component equation.
The importance of the specification of mass (given by the surface pressure
—D
and temperature) on the lateral boundaries rests in the term fv , which
&
(neglecting the covariance of , ^ j—^ with -^— across the domain)
simply represents the difference in surface pressure and geopotential
(which is determined by the temperature) across the domain. Thus the mass
-------
31
variables on the lateral boundaries determine (together with the mean
Coriolis acceleration) the net linear acceleration of the air in the domain.
Of course this is true for any non-periodic domain, but the larger the
domain, the greater error that may be tolerated. This point is illustrated
in Figure 2.2 which shows the errors in the calculated geostrophic
wind as a function of horizontal distance when RMS temperature
errors of 1°C are hydrostatically integrated over a 200 mb depth centered
at the given pressure level. For synoptic scale domains (L >~ 3000 km),
the error is about one ms at all levels. The error increases rapidly
as the domain size decreases below 1000 km. For a regional scale model,
with a domain size of 1000 x 600 km, a 1°C temperature error on the lateral
boundaries will produce large (-5 ms ) errors in the mean geostrophic
wind across the grid, especially in the upper levels. These errors,
if they persist in time, may produce large erroneous accelerations of
-4 -1 -1
the mean motion over the domain. For example, with f = 10 s , a 5 ms
error in geostrophic wind which persists for six hours would produce an
erroneous change of 11 ms in the mean velocity over the domain.
The preceding analysis indicates that although it might be tempting
for computational stability purposes to calculate the temperature on the
lateral boundaries by some type of outward extrapolation from the interior
of the domain, or one-sided differencing scheme, the extreme sensitivity
of the mean acceleration to small temperature errors leads to us reject
these possibilities. Instead, we conclude that for mesoscale domains,
it is preferable to specify the mass variables at all boundary points in
a realistic and accurate way, in order to insure that the mean geostrophic
-------
32
i
Jt
Figure 2.2 Error in geostrophic wind as a function of horizontal distance
for a 1°C temperature error integrated over 200 mb depths
centered at 300, 500, 700 and 900 mb
-------
33
wind and the associated linear acceleration term remain accurate during
the integration. This specification may come from a large-scale model
(one-way interaction), or in the research mode, from consistent analyses
of observations.
While the specification of surface pressure and temperature on the
lateral boundaries affects the pressure gradient force terms, the importance
of accurately specifying the velocity components on the boundary is
exemplified by the term
2 _ u 2
E W
Ly
= B. in (2.29).
L - k
x
To illustrate the magnitude of this term, let u = u + Au, where
E W
Au may consist of either a physical variation across the domain or an
error. Then, neglecting the variation of u with y, the boundary flux term
is
2
Au + Au ]
\=- L
x
With moderate to strong winds when the non-linear terms are
important, Au is smaller than u so that the boundary term is of order
uT7 A
W u
L
x
Figure 2.3 shows the magnitude of B, (divided by f equal to
-4 -1 —D
10 s for comparison with v discussed earlier) as a function of wind
o
speed and domain size, L, for a value of Au equal to 10% of the mean wind
speed.
-------
34
/ O
Non-linear accelerations greater than 10~ ms~ (corresponding
to a value of B,/f equal to 1 ms~ in Fig. 2.3) become comparable in
magnitude to the linear acceleration term, f(v-v ), especially if they
O
persist over the entire 12-hour forecast. From Fig. 2.3, we see that
this value is exceeded by a considerable margin for wind speeds in
_^
excess of 30 ms and mesoscale domain sizes of the order 1000 km or
less. On the other hand, for domain sizes in excess of 5000 km, the
mean non-linear accelerations are very small, except when the mean wind
approaches 100 ms . In the discussion of the preliminary 3-D experiment
(see Section 3), we will see that in the upper layer where the wind speeds
are greatest, the values of B, /f are of order 5 ms and therefore of
K.
comparable magnitude to the linear acceleration over the domain.
2.6.2. Lateral boundary conditions for the 2-D model. A technique
which has found widespread use in theoretical numerical calculations is the
specification of cyclic or periodic boundary conditions in the direction
of the mean flow. This treatment has the desirable numerical property that
all calculations utilize the same finite difference equations and therefore
have similar truncation errors. Physically, however, these conditions imply
an infinite repetition of solutions downstream. If the domain is not large
enough, disturbances generated in the interior of the domain may propagate
downstream, pass through the cyclic boundary, and eventually return to the
interior of the domain. For this reason, we have chosen to experimentally
determine a set of non-periodic, open conditions that allow gravity waves
to pass out of the domain.
A set of open lateral boundary conditions which appear to give
computationally stable results was found through experimentation. A sug-
gestion that these conditions do not adversely affect the solutions on the
-------
35
10
IOOO
3900
Figure 2.3. Magnitude of acceleration due to boundary flux of momentum
as function of wind speed and domain size (L). The isotachs
are labeled in ms"-"-.
-------
36
interior to a large degree is supported by the experiments in which the
lateral boundaries are moved farther away from the interior of the domain.
(See Anthes and Warner, 1974, p. 52-56).
In adiabatic experiments (Q = 0 in (2.8)), the surface pressure and
temperature on both the inflow and outflow boundaries are specified as a
smoothly varying function of time. This specification prevents any
spurious pressure difference across the domain, which would lead to an
unrealistic net acceleration of the flow.
For the velocity components, the values at the inflow points are
specified in a manner similar to the specification of temperature and
pressure. The values on the outflow boundary are obtained by a Lagrangian
extrapolation outward from the interior. These boundary values are required
only in the computation of the non-linear horizontal momentum flux divergence
terms; they are not required in the computation of the horizontal divergence.
In fact, this is an advantage of the horizontally staggered grid. The first
computation point for pressure utilizes only velocity points interior to the
grid. Since the divergence calculation is very sensitive, this property
appears to be quite important in obtaining smooth solutions near the
boundaries. Also, these boundary conditions allow gravity waves to
propagate out of the domain in both directions, with no reflection at the
boundaries.
2.6.3 Lateral boundary conditions for the 3-D model. The lateral
boundary conditions for the 3-D model are more complex than those for the
2-D model, because of the staggering of the grid in two dimensions. Again,
the surface pressure and temperature on all the boundaries are specified
as smoothly varying functions of time. Boundary values for the velocity
components are needed only in calculating the non-linear horizontal flux
-------
37
divergence terms, such as 8(up*) a/3x where a may be u or v. For evaluation
of these terms in the preliminary experiment to be described in Section 3, the
quantity a is specified on the inflow boundaries, and extrapolated outward on
outflow points. The pressure-weighted flux velocity component normal to
the boundary ((up*) on the east-west boundaries and (vp*) on the north-south
boundaries) is extrapolated outward regardless of the direction of flow.
The above lateral boundary conditions on the velocity components, while
yielding computationally smooth solutions, produced a large source or sink
term in the kinetic energy budget for each layer, as discussed in Section 3.
In order to obtain realistic values for the net boundary flux of kinetic
energy, it is probable that the velocity components should be specified
on all boundaries, even if this amounts to an overspecification. In any
case, we find that because of the large contribution to the kinetic energy
budget by the boundary flux term, great care must be taken in the actual
values of the velocity components used in the specification.
2.7 Initial conditions
The initialization of numerical models with real data is usually
considered to be a problem of obtaining an initial quasi-balance between the
mass and momentum fields so that the amplitude of gravitational oscillations
is minimized. Although it is certainly desirable to minimize the small
scale imbalances, we have found that for mesoscale models with open
domains, the presence of small scale imbalances is less injurious than with
large-scale models in closed or cyclic domains. In short, the mesoscale
model adjusts to initial imbalances very rapidly through the propagation
-------
38
out of the .domain of internal and external gravity waves. More important
than the small-scale imbalances appears to be the large-scale imbalances
induced by the incompatibility of the lateral boundary conditions on mass
and the mean wind over the domain.
These large-scale imbalances produce low-frequency inertial gravity
waves which affect the mean motion over the domain to a large extent.
These somewhat surprising results are discussed in Section 3.5, which
presents preliminary results from a real data experiment.
Prior to the experimentation with real data, fictitious initial
conditions are specified in order to study the effects of topography,
diffusion parameterization, and the numerical aspects of the model. Both
the artificial initial conditions and the real data initial conditions
consist of specifying the sea-level pressure and the three-dimensional
temperature and wind fields. The surface pressure on terrain above sea
level is calculated iteratively from the sea-level pressure and the vertical
temperature profile, using the hydrostatic equation and alternatively
improving the surface pressure and mean temperature estimates until
convergence occurs.
2.8 Two-dimensional flow across the Appalachian terrain
To illustrate the nature of the perturbation flow induced by the
smoothed Appalachian terrain, the importance of the fluxes of kinetic
energy across the boundaries, and some of the aspects of mesoscale model
initialization, the results from two experiments with the 2-D analog are
presented.
-------
39
2.8.1 Specifications of the 2-D experiments. In both 12-hour
experiments, the model was run with six layers (velocity components defined
at a = 0.95, 0.85, 0.7, 0.5, 0.3, and 0.1), p = 0, and the temperature
sounding for Washington, D. C., Oct. 16, 1973 (shown in Fig. 2.4). The
horizontal grid is constant (20 km) over most of the interior of the domain
and stretched at both ends to increase the distance of the lateral boundaries
from the central region of interest. The x-coordinates of the grid points
at which the velocity is defined are given by
0, 120, 220, 300, 360, 400, 420...1360, 1400, 1460, 1540,
1640, 1760 km.
In these experiments the effects of vertical diffusion of momentum
are modeled according to
FU = - g —z— k = 1,2...6 (velocity levels)
where (2.31)
T = p(z) K |U k = 1,2...6 ( levels)
zx z oz
and
T = p C |vJ u k = 7 (surface stress) . (2.32)
ZX S D -vO t)
Similar expressions are written for the v-component.
-------
40
zoo
4*0
t
v^
a.
300
(0*0
-40
Figure 2.4 Temperature sounding for Washington, D. C., I2Z Oct. 16, 1973
-------
41
In (2.31) and (2.32), T is the stress in the x-direction due to
zx
vertical mixing, p(z) and p are constant densities, C is the drag coefficient
-3 -3
(equal to 3 x 10 over land and 1.5 x 10 over water), and K is the
z
2 -1
eddy diffusivity coefficient and is equal to 25m s
The above crude formulation of the PEL is to be utilized while a more
realistic parameterization is developed.
In the 2-D model the horizontal diffusion terms take the form
FUR = p* K^ -^
(2.33)
s2
FV = £* K 9 V
•a P ^u 9 *
n n r\ £•
and are added to the momentum tendencies mainly to prevent the accumulation
of energy in the short wavelengths. In these experiments, ILj varies with
grid size according to the formulation
[5 x 104 + k^
9V
,22-1 /o o/\
) m s , (2.34)
min
where k =0.4 and Ax . equals the constant grid spacing of 20 km over
o mm
the interior of the domain.
The lateral boundary conditions consist of steady state pressures
and temperatures on both the eastern and western boundaries On the
outflow (eastern) boundary, the velocity components are extrapolated from
-------
42
the interior using a Lagrangian extrapolation. On the inflow (western)
boundary, the velocity components are set equal to the next inner value
provided the speed at the inner point is less than the geostrophic wind
speed. This procedure prevents an unrealistically large horizontal
shear from developing in the boundary layer where surface friction reduces
the wind speed to roughly 60% of the geostrophic value.
The only difference between the two experiments is the initialization
of the winds in the boundary layer. In Exp. 2D-1, the u-component is set
_i
equal to a geostrophic value of 30 ms at all points. As surface friction
destroys the geostrophic balance, the boundary layer winds undergo a large
oscillation before reaching quasi-steady values, which is undesirable
for a short-range forecast. Exp. 2D-2 illustrates that by taking surface
friction into account initially, the amplitude of the oscillation during
the adjustment period may be reduced.
As revealed by previous experiments on this scale with the early
version of the regional model, an undisturbed initial flow adjusts to a
quasi-steady state after roughly six hours. Actually, the solution never
reaches a steady-state, because high frequency inertia gravity waves are
continuously propagating through the model domain. Furthermore, for
long integrations, the absence of any effective mechanism (horizontal
and vertical diffusion are too small) for removing the upward propagating
energy near the top of the model, produces a gradual modification of the
solutions by the reflection of energy at the top of the model. Aliasing
of the vertically propagating waves produced by the inadequate vertical
resolution near the top also contributes to the problem for long inte-
grations. However, in most of these experiments with steady forcing, the
-------
43
most rapid adjustment occurs during the first six hours, the changes in the
low-level solution between six and twelve hours are considerably smaller
than between zero and six hours.
In order to minimize the effect of the transient part of the solution
due to traveling gravity waves, and to emphasize the stationary portion
of the perturbed large-scale flow, it is convenient to present time-averaged
cross sections from the last six hours of the forecast. These are obtained
by averaging the solution at one hour intervals from 7 to 12 hours. Be-
cause the high frequency gravity waves affect the vertical velocity field
more strongly than the horizontal velocity or temperature fields, the averaged
cross sections of vertical velocity show a greater departure from the hourly
solutions than do these latter sections.
2.8.2 Results with geostrophic initial conditions: Exp. 2D-1.
Figures 2.5, 2.6, and 2.7 show the time averaged cross sections of u,
temperature departure (from the horizontally isothermal initial state),
and vertical velocity. Even with only six layers, the classic characteristics
of perturbed flow over a wide ridge are evident. (Queney, 1960; Hovermale,
1965). A maximum in the wind velocity over the crest of the mountain is
crudely resolved. In the boundary layer, surface friction reduces the u-
component to about 60% of the geostrophic value.
The maximum cold temperature anomally of 3°C occurs at the top of
the mountain, while warm air tilts upward from the Piedmont region over
the top of the mountain.
The time-averaged vertical velocity cross section indicates a mean
maximum updraft of over 10cm s over the upwind slopes, extending to an
upper pressure level of about 700 mb. This mean upward motion in the
lower troposphere appears to be sufficient to account for the low deck of
-------
44
ID
_ in
to
N
oo
Oi
ci
o
II
b
8
oo
o
o
o
o
CM
o
o
00
8
tO
o
o
8
CO
o
CM
e
o
o
0)
to
en
03
o
rt
o
a)
M
n)
M
a)
n)
01
t-t
•H
PL4
o
o
CVJ
o
o
CO
-------
45
s
00
o>
d
{g •' 8 "S
O' O I"
b b b
8
00
O
O
u>
O
O
O
O
CM
O
O
CO
O
O
(0
O
O
, O
8
O
O
CM
O
O
O
O
O
O
00
ca
ex
0)
0)
M
3
4J
CO
>-J
01
O
Q
n) ON
i
01
V-i
3
00
•H
PL*
O
O
O
-------
46
E oo
* 10
— ro
N
ID
q
N-'
o
o
00
o
o
to
o
o
o
o
CM
o
o
CO
o
o
to
o
o
o
o
CM
•H
o
o
fl)
o
•H
HI
§
•H
4-t
O
ft
rt X
•H O
H M-l
CN
0)
3
60
•H
Pn
o
o
CM
o
o
o
o
CO
o
o
CO
o
o
o
-------
47
stratocumulus clouds that is frequently observed over the Appalachians in
strong northwesterly flow. Mean subsidence of 10cm s occurs over the lee
slopes.
We now consider the components of the kinetic energy budget in Exp. 2D-1,
shown for each hour of the 12-hour forecast in Fig. 2.8. Because the
western boundary of the domain is at a higher elevation than the eastern
boundary, the conversion of potential to kinetic energy by the u-component
2
(labeled A) is positive over the domain (about 20 watts/m ). This generation
is more than offset by the frictional dissipation (E), which represents
mainly the effect of surface friction. The conversion of potential to
kinetic energy by the term (fvu p*) in (2.27) oscillates slowly around zero
o
as the v-component adjusts to the ageostrophic u-component produced by
surface friction and the mountain-induced wave.
The net flux of kinetic energy across the lateral boundaries is
represented by B in Fig. 2.8. This term is positive, indicating that
kinetic energy is being advected into the region. The magnitude of the net
boundary flux exceeds the sum of the generation term and the frictional
dissipation term, and in fact, the total kinetic energy tendency (D)
closely follows the variation of the boundary flux of kinetic energy. This
large flux of kinetic energy across the boundaries, even in this simple
experiment in which there is no variation of the geostrophic wind across
the domain, is one indication of the importance of the boundary term. (In
an identical experiment, in which the velocity at the inflow boundary was
held constant rather than allowed to decrease as friction reduced the
velocities on the interior, the boundary flux was much larger and the wind
speed significantly higher far into the domain.)
-------
48
50 i
KINETIC ENERGY BUDGET 2~D MODEL 13 RUN FEB. 21, 1974
40
30
CVJ
20
10
0
-10
- A
-20
-30
-40
-50J
5678
TIME (H.)
10 II
12
A — GENERATION BY CROSS ISOBAR FLOW (u COMPONENT)
B - NET FLUX OF K ACROSS LATERAL BOUNDARIES
C — GENERATION BY CROSS ISOBAR FLOW (\i COMPONENT)
D — TOTAL KINETIC ENERGY TENDENCY (A + B + C + E)
E — DISSIPATION OF KINETIC ENERGY BY FRICTION
Figure 2.8 Kinetic energy budget for Exp. 2D-1
-------
49
2.8.3 Initialization of the boundary layer winds considering the
effects of surface friction. One of the difficulties with mesoscale
forecasts is that the short time scale of the forecast means that greater
care will have to be taken with the initialization procedure; a 6-hour
"adjustment" period is undesirable in a 12-hour forecast. If only gravity
wave modes, which traverse the entire domain in an hour or so, were
present, the adjustment period would be tolerably short. However, inertial
modes with periods of roughly 15-20 hours in middle latitudes are also
present, and for short-range forecasts, it is desirable to minimize the
amplitude of these waves that are generated by initial imbalances.
One source of initial imbalance is the sudden introduction of surface
friction in a geostrophic current. Fig. 2.9 shows the x-profiles
of the u-component in the boundary layer at 0.67h, 3.33h, and the 7-12h
average in Exp. 2D-1. The u-components decrease rapidly at first from
the geostrophic value of 30 m s , but instead of converging monotonically
toward the 7-12 hour average, they "overshoot" their equilibrium values.
This large imbalance and subsequent overshooting produces a large
amplitude inertial wave which persists throughout the 12-hour forecast.
To illustrate the nature of this fractionally induced oscillation, we
consider the following simplified boundary layer equations in which the non-
linear terms are neglected, and the geostrophic north-south wind component
is zero,
~ = fv - Ku (2.35)
-------
50
o>
>
o
NOUVA313
o
ro
8
o
o
in
CO
3
O
•H
4J
ni
cu
^
ca
13
o
•H
e
0)
o
W iH
3i
•H CN
O •
M CX
CX X
w
4J
en a
a) -H
? CO
4-J CU
cn B
W
60
-------
51
f(u - ug) - Kv (2.36)
where K represents a frictional coefficient. A multiplication of (2.36)
by i = /^T and addition to (2.35) yields
3V
^+ [K + if] V = ifVg, (2.37)
where V = u + iv and V = u + iV . The general solution is
~g g g
V = c& (K+ lf)t + (2.38)
and with the initial conditions that V = V at t = 0, the solution is
~g
KV f . ifV
V =2 e~(K lf)t +-^8-
~ (K + if) e K+if
A separation into the real and imaginary parts yields
u = g v , {1 + e Kt [(f)2 cos ft - | sin ft]} (2.40)
[1 + (f)2] f f
a T^ — ITt" TT
v = S 9 (f) (1 + e ^ [-(f) sin ft - cos ft]} (2.41)
[1 + (f)2] f f
-------
52
The time dependent solutions (2.40) and (2.41) contain a transient part
(an inertial oscillation which damps with time) and a stationary part, given
by
U
u(t + ») g (2.42)
[1 + (f) ]
U
v(t - ") - g g , & . (2.43)
[1 + (f) ]
If the friction is abruptly introduced in a geostrophic current, the transient
part is meaningless as far as the mesoscale forecast is concerned. The time
constant for the damped oscillation is 1/K. In the simple PEL parameterization
in which the surface stress is given by the quadratic stress law (eq. 2.32)
and the stress is assumed to vanish at the top (H) of the PEL, K takes the
form in (2.35) and (2.36)
K(t) = -^ (2.44)
so that K is not strictly a constant. As an approximation, let
|v(t)| = |V | = 20 ms'1, C = 1.5 x 10~3, and H = 1 km, which yields a time
-v -~g U
constant of 9.3 hours. Fig. 2.10 shows the hodograph of the non-linear
solution (obtained by a numerical integration of eqs. (2.35) and (2.36))
for these parameters. The three hourly positions are indicated by the symbol
x. As shown by Fig. 2.10, the u-component decreases rapidly during the first
three hours and overshoots its equilibrium value by about 4 ms . This
behavior is similar to the variation of the boundary layer u-component in
the 2-D model.
-------
53
Figure 2.10 Hodograph of non-linear solution to simplified boundary layer
equations with geostrophic initial conditions
-------
54
If it is desirable to eliminate or at least minimize the amplitude
of the transient part of the boundary layer solution, the model must be
initialized to include the effects of surface friction. In the simple model,
this may be done by initializing the u- and v-components with their linear,
steady state values. Fig. 2.11 shows the time variation of the hodograph
for the non-linear model with the initial conditions given by the solutions
CD 'V I
to the linear model (K = —TT % )< The amplitude of the oscillation is
n
reduced by an order of magnitude.
Of course in the mesoscale model, the non-linear advective terms and the
terrain induced gravity winds interact with the boundary layer winds in a
complicated way, nevertheless the amplitude of the frictionally induced
inertial oscillation may be reduced by initializing the u and v components
according to (2.42) and (2.43). Exp. 2D-2 is identical to Exp. 2D-1 except
that the boundary layer winds are initialized in this way, with u = 20.1 ms and
v = 9.8 ms . Fig. 2.12 shows the time evolution of the x-profiles of u in
Exp. 2D-2. Compared to Exp. 2D-1, the amplitude of the initial oscillation
is significantly reduced. Although the solutions in Exp. 2D-1 and 2D-2
differ considerably during the first six hours, the 7-12 hour averages are
very similar in both experiments since the time constant for C = .003,
|V | =20, and H = 1 km is about 4.6h.
o
-------
55
$•/•/<•/0«i
S.I
Figure 2.11 Hodograph of non-linear solution to simplified boundary layer
equations with initial conditions given by linear solution
-------
56
d>
o
NOHVA313
C\J
en
3
O
•H
(U
^
CO
\J
/
cd
"O
fi
o
c
OJ
c
o
p.
o
o
o
M CM
0) I
•H Q
•H CM
O.
x
w
en
a)
S
l
j-i
en
n3
W
CM
0)
3
00
•H
/
/
o
CJ
o
I
n
-------
57
APPENDIX - CHAPTER 2.0
EQUATIONS OF MOTION FOR THE LAMBERT CONFORMAL PROJECTION
The Lambert Conformal projection is obtained by projecting the image
of the earth's surface on a right circular cone which intersects the
earth's surface at two latitudes, $.. and _. The axis of the cone
coincides with the north pole. The cone is then cut along one side
and flattened. A cartesian (x - y) grid is oriented on the flattened
image, and the latitude and longitude coordinate of any point related
to its cartesian coordinate (x,y).
Normally, the x-axis (y = 0) coincides with 0° longitude,
and the relationships between x and y and latitude (((>) and longitude (X)
are given by
x = r cos X (1)
*
y = r sin X (2)
*
X = n X (3)
rtan \b/2 ,n
[tan /2]
where X is the longitude (positive east, negative west), n is the constant
of the cone, fy is the colatitude of , fy is colatitude (90° - ), and a is the
radius of the earth.
-------
58
The above projection produces an x - y grid for which the y-axis is
oriented north-south at longitude 125.7°W. Since the regional model is
centered near 80°W, the x and y directions on this grid are significantly
different from the east-west and north-south directions on the earth. There
fore, in our model, the projection is rotated 45.70° so that the y-axis
coincides with 80°W longitude. To accomplish this rotation, we define
X"= X-- 45.70 (5)
X' = nX"
(6)
The relations (1), (2), and (4) remain the same. With this rotation,
the point at 35°N 80°W is located on the grid at x = 0, y = -7157 km
(See Fig. 2.1).
To derive the equations of motion for Lambert Conformal coordinate,
we make the following definitions:
x,y,z map coordinates (projected system).
x ,y ,z curvilinear earth coordinates, x positive eastward,
SSS -- - - -- •«
y
positive northward, z positive upward,
U,V,W velocity components relative to earth along the x ,y ,z
directions. s
u,v,w velocity components relative to earth along the x,y,z
directions.
The equations of motion are
£Q _ UV tan cf> + M = _ a |2_ + fv _ 2 QW cos * + Fx (7)
dt r r 9x s
e e s
, VW = _ _ _
dt r r 8y *
e e s
. _ a _ - g + 2 un cos cj) + Fz ,
g
e s
-------
where U = r cos
59
dX
dt
V= reif
dr
dt
dx = r cos d> dA
s e r
dy = r d
s e r
dz = dr
s e (10)
and r is the radius measured from the center of the earth, ft is the
e
angular velocity of the earth, a is specific volume, and Fx and Fy are
S S
the frictional forces in the curvilinear earth coordinates.
Now, from the definitions of x and y in (1) - (2) , it may be
shown that
,dx,. /-cos A1 -sin A' cos d>N /ddK
= m r (- ' ' > (dA}
N /
sin A' cos A' cos *> (dA}
,-cos A' -sin A' cos (JK , s.
= m (-sin A' cos A' cos 4>} (dx }
S
tan
" " (13)
Furthermore, the relationship between the velocity components in the
map coordinates (u,v) and in the curvilinear earth coordinates (U,V) is
,u.. _ /-sin A' -cos A\ ,u\
(v) ~ ( cos A' -sin A'} V (14)
-------
60
which may be inverted to yield
,Us ..-sin A' cos A\ ,u«
V ~ (-cos A' -sin AJ V • (15)
Next, the relationship between the acceleration in the (x,y) and
(x ,y ) coordinates must be found. A differentiation of (14) yields
S S
' \ I \
du\ /dU\
dt I _ /-sin A' -cos A'M ?M dA' /O -1\ ,iu
^/~VcosA' -sin XV I g-/ dt Vl O) V
dt j \ dt/ /
Finally, the relationship between the derivatives in the two coordinates is
*s
9x I -11 -11
s 1 , -sin A cos A
3 = m ( -cos A' -sin A' M 3 j (17)
Now, with the use of (17), (16) and (14), the horizontal component equations
of motion may be written
= - m (u + v ) + v (Y + f) - ma + w 2^ cos $ sin A' -
e
- Fx sin A' - Fy cos A' (18)
-------
61
+ v ?> + (Fxs cos x> - Fys sln r)
(19)
f)u-maj-- w/ (2^ cos cos X ' + — )
-
dw 9 9
The derivation of the equation for -r- is much simpler, since -^ - -^
s
w = W, so that after substitution for u,v,and w in terms of U,V, and W,
we obtain
2 2
dw u + v . on A / uy . vx x 8p . _
—- = h LU cos tp ( - —*- + — ) - a TT^— g + Fz
dt r Yrr3zes
where r is given by (4). For hydrostatic systems, in which (20) is approximated
by
01 l^+ g = ° (21)
12 2
the kinetic energy equation will contain -r- (u + v ) rather than
1222
•=• (u 4- v + w ). Hence, for consistency, the small terms involving w in
the horizontal component equation (18) and (19) are usually neglected.
To obtain the equations given in (2.2) and (2.3), the vertical coordinate
z is transformed to a (for example, see Haltiner, 1971) and the equations
are written in flux form with the aid of the continuity equation.
-------
62
3.0 PRELIMINARY THREE-DIMENSIONAL EXPERIMENT
USING REAL DATA WITH AND WITHOUT TERRAIN
In this section we describe our first attempt at initializing the
3-D regional model with real data. Since the primary purpose was to
test the numerical stability of the model under strong winds and time-
dependent boundary conditions for a 12-hour forecast, a very crude
initialization scheme was tested. This scheme consisted of simply
interpolating data from large-scale, hand analyses to the 20-km grid.
The reaction of the model to this set of unbalanced data has revealed
several important aspects of this model in particular and mesoscale modeling
in general. For example, it has shown (1) that the model is quite stable
numerically, (2) that high frequency, short-wavelength noise does not
cause much difficulty even with unbalanced data, and (3) the use of
unbalanced data apparently produces large linear and non-linear acclerations
of the mean flow over the domain. Because of the crude initialization
scheme, the forecasts do not appear to be very accurate, although lack of
verification data on the scale of the model hampers their evaluation.
Nevertheless, a comparison of 12-hour forecasts with and without the
complex terrain of the Appalachians, but identical in every other
respect, indicates the importance of terrain on this scale of motion.
-------
63
3.1 Synoptic Discussion on 12Z Oct. 16 - OOZ Oct. 17, 1973
As part of the SRG research program, The Pennsylvania State
University conducted a pilot field experiment during the period
Oct 16-17, 1973. One of the purposes of this program was to
obtain data to supplement the conventional surface and upper air
data network and to be used in the verification of the regional
model. Part of the program included the flying of the Penn State
aircraft over the central portion of the regional model's domain.
Although it was realized that one aircraft could not provide an
adequate data base for mesoscale model verification, we hoped to
obtain enough low-level data to partially verify the model and to
aid in the identification of the scales of the perturbations caused
by moderate to strong flow over rough terrain.
From the beginning of the regional model development, we have
concentrated on the modification of synoptic-scale flow by gravity
waves produced when statically stable air flows over rough terrain.
Because the terrain in the regional model's domain contains
appreciable energy in the wavelengths resolvable by the model
(100 - 600 km), variations in terrain should be a major cause of
mesoscale variation in weather (temperature, wind, cloudiness,
rainfall, turbidity) in this area.
The synoptic situation chosen for the first test of the model
was one in which terrain-induced gravity waves were expected to
produce measurable (and hopefully predictable) mesoscale perturbations
to the large scale flow. Upward propagating gravity waves are
favored over mountainous terrain when the basic current is strong,
-------
64
statically stable, perpendicular to the ridge, and nearly constant
in direction with height. A strong northwesterly flow of cool air over
the eastern third of the U. S. satisfies these requirements, since the
mean terrain contours are oriented from the southwest to northeast.
Furthermore, a strong cross-domain flow was desired to test the computa-
tional stability of the model under strong advection conditions, which
are known to be difficult to handle numerically.
The weather situation chosen for the field program and subsequent
modeling tests included the two day period from Oct. 16 to Oct. 17, 1973,
and satisfies the preceding conditions quite well. During this period,
northwesterly flow behind a cold front intensified slowly. At 12Z Oct. 16
the mean wind speed increased with increasing elevation from about
15 ms at 850 mb to 50 ms at 300 mb. During the next 12 hours, a
strongly baroclinic zone and an associated jet maximum of 65 ms drifted
southward through the model domain, providing a good test for the time-
dependent boundary conditions.
The lower tropospheric flow during this period is shown in the 850 mb
maps for 12Z Oct. 16 and OOZ Oct. 17 (Figures 3.1 and 3.2). Although
the temperature gradient is rather weak at this level, cold advection
increases throughout the period. Furthermore, the geostrophic winds
become more northerly during the 12 hours as a broad trough moves east-
ward. As discussed later, this tendency for increasing northerly winds is
correctly predicted by the model.
-------
65
o
O
si
C--1
(-1
o
4-1
0)
•H
Cfl
CO
P
ra
o
LO
oo
0)
H
-------
66
o
o
O
O
M
O
CD
•H
CO
cfl
fi
CO
O
m
oo
0)
CM
00
-------
67
The middle tropospheric flow at 12Z Oct 16 and OOZ Oct 17
is illustrated in the 500 mb maps shown in Figures 3.3 and 3.4.
The upper level front is much sharper at this level, as the
packing of the isotherms over Pennsylvania and New York in
Fig. 3.3 indicates. During the 12-hour forecast period,
this intense baroclinic zone drops southward into the mesoscale
model domain — the temperature at 500 mb along the northern
boundary falls by about 10°C. At 500 mb, as well as at 850 mb,
the winds veer with time, becoming more northerly.
The upper level flow is shown in the 300 mb maps at 12Z
Oct 16 and OOZ Oct 17 (Figures 3.5 and 3.6). The upper level
front reaches the northern part of the mesoscale domain by
OOZ Oct 17. The wind speeds over the model domain increase with time
as the baroclinic zone in the middle and lower levels intensifies.
Thus the synoptic pattern of Oct. 16 and 17 should be quite
favorable for the production of mesoscale perturbation^ to the mean
flow. Northwesterly flow of stable air increased in speed with
height but changed direction only slightly. Strong cold advection
in the middle troposphere was present to check the model's reaction
to time-dependent boundary conditions.
3.2 Initialization and verification analyses and specification of
time-dependent boundary conditions
As mentioned earlier, the initial data for the 30 x 50 grid
pictured in Figure 2.1 were obtained by reading temperature and
-------
68
o
o
-------
69
o
o
N
O
O
O
14-1
co
•H
CO
O
O
in
0)
,e
H
ro
a
•H
-------
70
oo
r»-
CTi
O
O
CM
CO
•rl
ca
a
•§
O
O
en
m
•
CO
01
M
60
•H
-------
71
CO
O
O
O
O
O
m
CO
tfl
s
•i
O
O
CO
0)
M
•H
-------
72
wind data for every tenth grid point from the hand analyses at
the surface, 850, 700, 500, 300 and 200 mb pressure levels.
These data were linearly interpolated to the remaining grid
points on the surface and constant pressure charts. The data
were then linearly interpolated to the a-surfaces defined in
Table 3.1.
Table 3.1 Sigma surfaces for velocity and temperature
in 3-D model experiment (p = 250 mb)
Level Index
1
2
3
4
5
6
P - Pt
Ps - Pt
0.125
0.325
0.475
0.625
0.775
0.925
Mean pressure of
a - surface (mb) at t = 0, Exp 1
342
489
599
709
819
929
Figure 3.7 shows an example of the unsmoothed analyses of u
and the geostrophic u calculated from the analysis of the mass field
on the sigma = 0.325 surface. The observed u-field is fairly
smooth, but errors in the temperature analysis and the errors
introduced by the bi-linear interpolation scheme produce a very
noisy geostrophic wind field. It is noteworthy that even though
imbalances of over 10 ms~ between the mass and momentum fields
-------
73
•: •i
UJ
oa
o
I
8
HI
o
* 'S'Sl '". * '«"£*«*: * *
t": jiri's «sss/s j isis;:; j :=;
~.r.;~. i -.z::: *. °.t'.tf. - -;-.;- - -. •;-•;-
sH:s ; 5;S-~ - sjt/: s ~;=;= ~ ~ s:;:;
,z|t : :§;:« sjssyj JsiSi"; i j=;:;"2
0!
0)
c
01
c
o
o-
6
o
o
i
u
•H
a
o
1-1
4-J
to
O
(U
60
C
cfl
4J O
e
OJ II
c
O *J
o.
6 4J
o «
o
I /->
3 in
Cvl
T3 ro
01 •
> O
M
a) ii
tn
^ to
O ^
a)
M
-------
74
exist initially, the horizontal scale of the imbalances is so small
that the irregular mass field quickly adjusts to the smooth wind analysis.
This rapid adjustment of the mass (temperature) field to the wind field
on the small scale means that fine scale temperature perturbations, whether
real or introduced as a result of analysis errors, will quickly dis-
appear unless the small-scale features appear in the wind field as
well. Thus it is important to realize that a very accurate observation
and analysis of the temperature field will have little value in mesoscale
models unless the winds are initialized with comparable detail and accuracy.
The unsmoothed and unbalanced data were used as initial
conditions for a 12-hour forecast. The same procedure was used
to obtain analyses at OOZ Oct 17. The temperatures, surface
pressures, and wind components from the analyses at these two
times were used to specify the time dependent lateral boundary
conditions.
For temperature and pressure, the boundary values at every
time step in the forecast were specified by a linear interpolation
in time between the 12Z and OOZ analyses. The specification of the
wind components was somewhat more complicated. Based on experiments
with fictitious data, the most stable set of boundary conditions
for the velocity components consisted of an outward extrapolation of the
pressure-weighted velocity components regardless of direction of
flow, but treating the velocity components themselves differently
depending on whether the velocity component normal to the boundary
-------
75
was inward or outward. These conditions, needed only for the evaluation
of the non-linear momentum flux terms, are given in Table 3.2.
Although these boundary conditions produce smooth solutions,
the previous discussion on the importance of the boundary values of the
velocity components and the results to be presented later suggest that
a complete specification of the velocity components on all lateral
boundaries may be desirable, even if this amounts to an overspecification.
The resultant generation of short wavelength noise near the boundary can be
controlled by the eddy viscosity term, without significantly affecting the
mean motion over the domain.
3.3 Specification of parameters
In both experiments (with and without terrain) the upper pressure
surface, p , is 250 mb. The grid size is 20 km and the time step is
30 s. The vertical diffusion of momentum is neglected except in the
lowest layer of the model where the quadratic stress law, discussed in
Section 2.8, is utilized. However, the drag coefficient, C , was in-
advertantly set to rather low value (for this rough terrain) of
_3
1 x 10 , so that surface friction is probably underestimated in these
experiments.
The horizontal diffusion of heat and momentum, included mainly
for the purpose of preventing the accumulation of energy in the
short wavelengths, is modeled according to (2.15), (2.16) and (2.18), where
K^ is given by
4 *2
K = 5.0 x 1CT + -2 (As)' D (3.D
-------
76
0]
01
•H
^1
cd H
13 W
fi W
3 ^E
0
rO
,_,
cd
c
o
Jf
P,
13
a
cd
3
*
P./-V
^ fi
S
"3 W
> rH H
0 &
T3 CJ O
C3 c/2
cd t-i
O
IS
CD O
4J M
C
0) rl .
C 0
0 -H
f\t ^f
S
13
O ^3
O 4J
>, O
M 4-1
cd
13 rH
ti cd
3 3
o o-
.n 0)
rH 4J W
cd CU H
M CO pi
o> o
4-1 CU S5
ca t-i
iJ cd
CM
ro
O) Cri
rH <;
§
M-l
•H
01
•H
U-l
<4-l O T-I O
•H , O
V 01 A
CM p,
3 CM W CM
^^ x-\
n 3 n 3
* *
!-4 Pi 1— 1 Pi
13
CU
-a
cu
CU
£3
O
a
M-l
•H
M-l
•H 13
0) O
iH O -rl
1 , "4-1 V
g A) i-l
g rH 0) V
3 >*j S* «
II S II )S
P5 r^
W /•— \ y^N
II
rH P, iH
3^3
13
CU
Tj
0)
01
52;
4J
55
nj ti i
•H -H
rH O 13
1 01
>< A -H
<£ M-l
S rH -H
R 1 CJ
3 X 01
«rf P,
" s M
^ y^s ^d
« ^ ^
Is ^c IE
R cu R
3^3
rtf
t>>.
o
A
CM
^-x
^
*
ex,
»^/
0
V
rH
1
X
^
s
pH
/-N
>
*
p.
^-^
<4H
•H
CM
>
II
rH
0
v|
CM
"3"
*
P.
13
01
•H
M-l
•H
CJ
cu
p.
CO
II
rH
0
A
CM
3
*
P.
14-1
•H
13
cu
13
cu
cu
0
A|
rH
1
>^
1
x-v
3
*
P.
13
0)
•H
M-l
•H
CJ
01
P.
CQ
II
>
O
V
r-
I
><
^
/-v
•K
P.
v../
13
cu
13
CU
CU
4H
•H
CM
>
II
rH
•O
CU
13
cu
0)
S3
4-1
O
M-l
•H
T3
0)
MH
0 iH
0
V 0)
P.
CM CO
">" II
*
P, rH
O
A
CM
K>
4<
P.
cu
13
cu
cu
S3
*\*
f>
-------
77
where kQ = 0.4, As = 20 km, and D = [ (|^ - |V + (~ + |V]1/2- After
an initial forward time step, the centered-in-time differencing scheme
is utilized without time smoothing. The model was run on the NCAR CDC 7600
system, and requires about 30 min. of computer time for a 12-hour forecast.
3.4 Qualitative discussion of results
3.4.1. Low-level results. Figure 3.8 shows the forecast surface
pressure change in the experiment with the smoothed Appalachian terrain,
which is contoured in Figure 3.10. (Hereafter, this experiment, which
includes terrain, will be denoted as Exp. 1; the identical experiment
except for the setting of the terrain elevations to zero will be designated
Exp. 1A.) The observed 12-hour pressure change is shown in Figure 3.9.
Both the forecast and observed pressures over the mesoscale domain
have risen during the 12-hour forecast; however, the forecast pressures
in the center of the domain have risen 2 mb more than the observed pres-
sures. This erroneous rise cannot be attributed to inaccurate treatment
of the terrain effects, since nearly the same rise is observed in Exp. 1A
when the terrain is absent. Instead this error is probably a manifestation
of the "pillow effect" (Benwell and Bretherton, 1968), and is related to
the inaccurate initialization procedure, and in particular the treatment
of the velocity components. This possibility will be investigated further
after a more sophisticated analysis of the initial data that produces
the time-dependent boundary conditions is tested.
The low-level velocity fields at 12 hours (OOZ Oct. 17) for
Exp. 1 and 1A are shown in Figures 3.10 and 3.11, and the verification
map is shown in Figure 3.12. In both experiments the winds have become
-------
78
« 10 . .- .0
Figure 3.8 Forecast 12-hour surface pressure change (mb) for Exp. 1 ending
OOZ Oct. 17, 1973
-------
79
o
o
Nl
O
O
oo
C
•H
T3
C
01
O
(U
(-1
60
•H
-------
80
Figure 3.10
Twelve-hour forecast velocities at level 6 (a =0,925, p = 939 mb)
from Exp. 1. Contours of terrain are labeled in m.
Figure 3.11
Twelve-hour forecast velocities at level 6 (a =0.925, p = 929 mb)
from Exp. 1A, no terrain
-------
81
i\ll>l
i V" T i
i i i i i i l
i i i i i i t
i i i i i i l l
i t i i t i l l
t i i i i i i l
l l l t l l l I
t t t 1 1 1 L
I I 1 I I I I
Figure 3.12 Verification velocities at level 6 (a =0.925, p = 929 mb) from
interpolation of synoptic analysis to mesoscale grid
-------
82
more northerly in the low levels in agreement with the observed tendency
for stronger northerly winds. However, both experiments overemphasize
the magnitude of the northerly winds, especially near the southern boundary,
where northerly components approach 15 ms . These speeds are much
stronger than the observed surface velocities. This error is related
to the erroneous pressure change discussed earlier; the north-south
pressure gradient is too strong. A second reason may be the weak effect
of surface friction caused by the small value of the drag coefficient. Never-
theless, the low-level winds do illustrate the effect of terrain, with
the terrain (Exp. 1) producing a lee-side trough, especially in the northern
half of the domain.
3.4.2 Middle level results. Under the influence of the vertically
propagating gravity waves, the effect of terrain is strongly evident
at middle as well as low levels. Figures 3.13 and 3.14 show the 12-hour
forecast flow at level 4 (p = 709 mb) in Exps. 1 and LA. Without terrain
the flow is smooth with a broad trough traversing the domain. With
terrain, however, a sharp mesoscale trough appears in the lee of the
Appalachians. This perturbation is so strong that the wind at several
grid points is reduced nearly to zero.
Because of the low vertical resolution there are undoubtedly serious
phase and amplitude errors in the vertically propagating waves. However,
these experiments show that the mesoscale terrain variations are capable
of producing mesoscale perturbations to the flow which can be resolved by the
model in a stable 12-hour forecast. For improved accuracy, the model could,
of course, be run with more layers.
-------
83
Figure 3.13
Twelve-hour forecast velocities at level 4 (o = 0.625, p = 709 mb)
from Exp. 1
Figure 3.14 Twelve-hour forecast velocities at level 4 (a = 0.625, p = 709 mb)
from Exp. 1A, no terrain
-------
84
3.4.3 Upper level results. Finally, the upper level flow (level 2,
p = 489 mb) is shown in Figures 3.15 and 3.16. Even at this level, the
terrain has produced a noticeable distortion in the large-scale flow.
Also notable in Figures 3.15 and 3.16 is the smoothness of the flow in
spite of the strong winds (60 m/s) blowing across the boundaries.
The temperature change field (Figure 3.17) in Exp. 1A at level 2
(p = 489 mb) illustrates the southward propagation of the observed in-
tense baroclinic zone and shows how the time-dependent boundary conditions
affect the temperature. On the whole, the temperature forecast is reasonable,
except near the eastern (outflow) boundary where 2Ax noise is apparent. This
noise is caused by the overspecification of the temperatures on the outflow
boundary and the subsequent upstream propagation of the short-wavelength
advective temperature waves. This noise could, of course, be reduced by
using a higher value for horizontal diffusion or by occasionally smoothing
the temperature forecasts in space. However, this noise is preferable to
the alternative of extrapolating the temperature outward to the boundary
from the interior, in which case the errors in temperature along the
boundary adversely affect the inner motion on the grid through the processes
described in Section 2.6. The effect of terrain on the temperature fore-
cast at all levels is to perturb the large-scale temperature field by
several degrees Celsius. Figure 3.18 illustrates the temperature perturba-
tions for level 2 in Exp. 1 with terrain for comparison with Figure 3.17.
3.5 Budget equations for the model domain and the implications of the
lateral boundary conditions
In the previous section, a sample of the 12-hour forecast results
in experiments with and without terrain were presented. The results,
-------
85
Figure 3.15 Twelve-hour forecast velocities at level 2 (a = 0.325, p = 489 mb)
from Exp. 1
Figure 3.16 Twelve hour forecast velocities at level 2 (a = 0.325, p = 489 mb)
from Exp. 1A, no terrain
-------
86
Figure 3.17 Twelve-hour forecast temperature change at level 2 (p
from Exp. 1A, no terrain
489 mb)
Figure 3.18 Twelve-hour forecast temperature change at level 2 (p 489 mb)
from Exp. 1, with terrain
-------
87
while emphasizing the importance of terrain and the "smoothness" of the
solutions, showed large accelerations in the mean motion at some of the
levels, indicating that large-scale imbalances between the mass and
momentum fields in the initial, and probably the boundary conditions, were
producing a "sloshing" of the winds and temperatures over the domain. In
this section we seek to explain these imbalances by an examination of the
mean kinetic energy budget for the domain and the time-dependent behavior
of the mean motion at various levels.
3.5.1 Mean kinetic energy budget for the 12-hour forecast period. The
mean kinetic energy budget for the entire mesoscale domain provides an
important check on the consistency of the model as well as gives some clues
on the importance of the various processes in determining the mean motion.
The components of the kinetic energy budget (eq. (2.25)) for Exp. 1 and
1A are presented in Figure 3.19. The total kinetic energy tendency evaluated
from the right side of (2.25) is shown for Exp. 1. As a check, a second
budget is computed by mulitplying the left side of (2.15) and (2.16) (the
actual tendencies used in the forecast) by u and v and integrating over
the domain. The total budget from this method is labeled LS in Fig. 3.19.
The difference between these two budgets is the uncertainty in evaluating
the boundary flux of kinetic energy terms in (2.25). As shown by the
close agreement of the two budgets, this uncertainty is small compared
to the magnitude of the boundary flux terms.
Although the presence of terrain greatly affects the components
of the generation of kinetic energy by cross-isobar flow, these two
terms tend to cancel and the net effect is a similar total budget for
-------
88
\i
13
w
T)
>>
to
M
Ol
a
-------
89
both experiments. (The net kinetic energy tendency for Exp. 1A would
be nearly indistinguishable from the tendency plotted for Exp. 1 in
Figure 3.19.) Thus the effect of terrain in this experiment is mainly
to redistribute kinetic energy within the domain rather than to create
or destroy large amounts of kinetic energy.
The striking feature in the kinetic energy budgets of both experiments
is the near cancellation between the large, positive production of
kinetic energy by cross-isobar flow and the large loss of kinetic energy
through the non-linear, boundary flux terms. Frictional dissipation by
horizontal diffusion and surface friction is very small compared to these
12
two components. As discussed in Section 2.6, the accleration of the mean
f\
motion over the domain is strongly dependent on the lateral boundary condi-
tions of mass and momentum. In the kinetic energy budget, we see this point
again; the net loss of kinetic energy through the lateral flux of kinetic
energy appears to be acting like friction, causing the winds to deflect
toward low pressure and resulting in a large, positive generation of kinetic
energy through cross isobaric flow. This point may be further illustrated
by considering the behavior of the mean motion at various levels in the
domain.
3.5.2 Time variation of the mean motion. Before presenting the results
for the time variation of the mean motion at various levels over the domain,
we present the solution for the mean equation of motion (2.28) and (2.29).
With the definitions
-------
90
(UE" UW } (uv]k - (UV)
S
LX Ly
_ [(uv)E- (uvyV (yg - v2g)
v L L~
x y
B = B + i B (3.2)
u v v '
V = u° + i v°
~g g g
the equation of mean motion may be written
9V
«-• = B - if (V - V ) . (3.3)
O t "" —' *^fi
The complete solution to this equation with B and V equal to a constant
is
iB iB
V(t) = V - -- + [V° - V + --] e"1"
This solution, which is pictured in Fig. 3.20, consists of an oscillation
with the inertial period about the stationary solution, V - iB/f. If
-------
91
TJ
(2
cfl
PQ
H-l
O
TO
t-i
60
O
C
cfl
4-1
0)
C
O
O
M
O
14-1
C
•H
ca
o
•a
o
§
•H
4-)
O
M
O
4-J
a
G
a
(3
O
•H
O
C/3
m
0)
3
oo
•H
-------
92
the large scale motion is not in balance with the prescribed boundary
conditions for mass (V ) and momentum (- iB/f), the mean motion over
~S ~
the domain will undergo the oscillation indicated in Fig. 3.20. It is
therefore desirable to minimize the departure from equilibrium, given by
[V° - V + iB/f]
~v rv o **"*
Fig. 3.21 shows the magnitude of the components B /f and B /f
at the various levels as a function of time for Exp. 1A. In the upper
levels this non-linear acceleration is very large, and verifies the
crude analysis of Section 2.6. Only in the lowest two layers
are the f-scaled, non-linear accelerations less than 3 ms . In the
upper layers, these scaled accelerations exceed 15 ms , dominating the
ageostrophic flow at these levels.
Although the boundary conditions for mass and momentum were not
constant in time, so that V and B varied over the 12-hour period, the
"" o ^*
type of oscillation present in the simple solution with constant B and V
is clearly present in the hodographs of the mean motion at the various levels.
For example, Fig. 3.22 shows the values of V and -iB/f, for the 12-hour
""& ""
period and the hodograph of the mean motion, V . During the 12-hour period
the mean wind swings toward low pressure, producing the large generation of
kinetic energy by the cross-isobar flow.
In the low-levels, the term B/f is small, and the mean wind is
observed to chase the time-dependent, mean geostrophic wind, as shown
in Fig. 3.23. We see now why the low-level winds become too strong from
-------
time
Figure 3.21 Components B /f and B /f (mean non-linear accelerations scaled
by f) for Exp. 1A v
-------
94
cat
•H
I
*o
c
CD
00
-o
co en
o II
a> icx
60
C
(0
4)
' 0)
I CO
Cfl
-a
5 K
C H
a
0) C
B -^
W-l CO
o cu
03 T-l
CO CO
VJ 3
00 O
O i-l
•a i-i
O CO
PC >
CN
CN
0)
I-I
00
-------
95
X
•O
a)
CO
CO
.. •*• -
O
•*'
4-1
•H
U
O
a
cfl
60
O
LT|
cfl
CD
P.
X
rH Cfl
n)
II U
(!)
•U S-l
O
T3 U-l
rt t-i
3
O O
.e
n i
60
CM
01
m
i
3
M
•H
-------
96
the north — the mean geostrophic wind, specified through the boundary
conditions on temperature and pressure, becomes strong northerly.
The conclusion from these results is that for mesoscale domains of
order 1,000 km across, the boundary values of mass and momentum may
produce very large accelerations in the mean flow over the domain. Unless
proper care in the initialization and specification of the time-dependent
boundary conditions is taken, the mean motion will undergo large
oscillations which may completely ruin the mesoscale forecast. In con-
trast to large domains, therefore, boundary conditions which simply minimize
the noise generated at the boundaries are not sufficient to guarantee
realistic solutions on the interior of the domain. These results suggest
that it is more important in the initialization and specification of lateral
boundary conditions to make sure that the large-scale motions are balanced
than to eliminate small-scale imbalances within the domain. Provided
that the winds are initialized accurately, the small-scale mass imbalances
will quickly adjust to the wind field. Large-scale imbalances, however, which
are forced by the lateral boundary conditions, may adversely affect the
solution for the entire forecast.
-------
97
4.0 NUMERICAL EXPERIMENTS WITH A TWO-DIMENSIONAL NESTED GRID
The motivation for nesting a fine mesh within a larger mesh is
obvious; the horizontal resolution in a numerical model may be con-
centrated where it is needed the most (e.g., over populated regions or
over atmospheric phenomena which have sharp gradients and fine scales).
Furthermore, by nesting a fine mesh model such as the regional model within
a synoptic-scale model, some of the deleterious effects of the proximity
of lateral boundaries to the center of interest may be avoided. On the
other hand, additional computational difficulties may be expected where
the grid meshes change size.
At least two approaches to the nested grid problem exist. In the
first, and simplest, a large scale model is integrated independently
of the fine mesh model, and the output from the coarse mesh model
utilized to provide time-dependent boundary conditions for the
fine mesh model. In the other approach, the model equations on the two
(or more) meshes are integrated simultaneously so that two-way interactions
are possible. Phillips and Shukla (1973) have argued from a theoretical point
of view that the simultaneous, two-way interaction approach is preferable,
and -HHur arguments tend to be confirmed by the results of Harrison and
Elsberry (1972).
However, the two-way interacting mode is considerably more difficult
to treat within operational constraints, and it is apparent from some
results (Hovermale, 1974) that a one-way interacting meshed model can produce
very good results. Furthermore, in both cases the difference in computational
phase speeds of waves on meshes of varying size produces problems along the
-------
98
interface. In the 2-way system, this noise can eventually contaminate both
grids. In the 1-way system, of course, the noise cannot feed back to the
coarse mesh. The problem of different phase speeds is aggravated as the
ratio of course mesh grid (CMC) to fine mesh grid (FMG) increases. Indeed,
it may turn out that for a ratio CMG/FMG of 2, the two-way interacting system
is preferable, while for larger ratios (say 5 to 10), the one-way interaction
is more desirable.
With the ultimate goal being the meshing of the regional model with
20 km resolution with a larger scale model (say with 100 km or greater
resolution), we first investigate the numerical problems associated with the
horizontal meshing of staggered grids in a simple numerical model - the
barotropic shallow-fluid equations. In spite of the simplicity of this
model, meteorologically relevant solutions corresponding to gravity waves,
gravity inertia waves, Rossby and Haurwitz waves exist, and the behavior
of these waves in the meshed model provides valuable insight into the
meshing of multi-level, baroclinic atmospheric models.
4.1 The basic equations
The basic equations consist of the well-known shallow fluid equations
of motion on an f-plane, written in flux form as
. duuh . duvh . ,9h - , _. r/-\\
•r + -5 + T + gh -r fvh= 0 (4.1)
dt dx dy dx
9vh 3uvh dyvh , 3h . , _ // o\
•~r f- —5— + —5— + gh -5 1- fuh= 0 (4.2)
8t 3x dy ° dy
-------
99
and the continuity equation
Here, the symbols, u, v, h, f, and g have their usual meteorological
meanings. In addition to the model described by these basic equations,
some experiments consider the advection of a passive quantity, A, governed
by the conservation equation
9hA , 3uhA 9vhA „ // /\
+ - + -= ° • (4<4)
4.2 The meshed grid system
Because the prediction variables in the regional model are staggered
in the horizontal (in order to economically reduce truncation errors), the
CMC and the FMG are both staggered, with h (and A) defined at x-points
and the velocity components defined at '-points (Figures 4.1 and 4.2). The
grids are meshed so that the momentum (•) points coincide. In the two-way
interacting experiment, tendencies at the CMC points that lie within the FMG are
not calculated. The FMG consists of 21 x 21 '-points and 20 x 20 x-points. The
ratio of coarse to fine mesh size is 4:1.
The finite difference equations for the two-way interacting model are
written to conserve fluxes exactly across the interfaces of each "box",
which are (in general) centered on a grid point. The boxes, representing
incremental volumes of mass or A, are centered on the x-points as shown
in Fig. 4.1. Because of the geometry of the meshed grids, each mass "box"is a
-------
100
•/ —
Figure 4.1 Boxes representing incremental volumes of mass in the 4:1 meshed
grid system
Q
•
°S
«i
J
«
C
1
•'
«
•
•
i— ,
r-*
—4
•-T
—t
t—t
•
&±*
'I
•IM^'I
• Q
M-.
I'J-
• e
!«
'•J
•
',
f>
.
••:
o
•
•I
• Q
,-i',
H-J
• o
'!'
3
•1-
'"I
— 1
-H
r^
h—^
-1
k— 1
•-^
k— •
!•"!
•
e
•
c
0
•
0
o
•
9
Figure 4.2 Boxes representing incremental volumes of momentum in the
4:1 meshed grid system
-------
101
square. However, with a staggered grid, it is not possible to preserve
square boxes for the momentum points (Fig. 4.2), and therefore, some of
the '-points near the interface represent irregular shaped volumes.
The finite difference equations across the interface of the various
boxes were written following Koss (1971). A difference between Koss's
model and this one is that Koss used the same time step for all grids .
In these experiments, At is varied with the grid size Ax to increase
the computational efficiency. Because of this variation in time step,
special care must be taken to insure that the fluxes of momentum and mass
are conserved across the interface between the CMC and FMG. The procedure
for the conservation of mass is illustrated for the grid volume labeled
i,j in Fig. 4.1. The continuity equation for this volume is written
-"-1-1/2, i- . .
where
,
, (uh)
(4.6)
Here T denotes the time step of the CMC (At), and T' denotes the
time step of the FMG (At/4). Thus, given values of uh on both grids at
time T, the continuity equation is first integrated forward on the FMG to
T + At in a series of four short (At/4) time steps. The fluxes of mass
across the FMG interface are accumulated and used to calculate the flux
divergence of mass at the CMG point. In this manner, at time (T + At),
the flux of mass (or A) into (out of) the CMG equals the flux out of (into)
the four adjacent FMG volumes.
-------
102
The conservation of the momentum fluxes is more complicated because
of the irregularly shaped boxes at the interface, particularly at the corners.
As discussed later, the boxes at the interface can be treated as CMC or FMG
points. The best results (to be discussed later) were obtained when these boxes
were treated as FMG points. If treated as CMC points, the variables
at the momentum points within the FMG are integrated first using
four time steps of At/4 to arrive at values of time T + At. The
fluxes across the interface wall between the FMG and CMC are accumulated
and used in the evaluation of the horizontal flux divergence of momentum
at the circled CMC points only. (Values of mass weighted momentum are
linearly interpolated between the circled --points for the FMG fore-
casts.) If the boundaries are treated as FMG points, it is not
necessary to accumulate fluxes, but it is necessary to assume that the
fluxes through the CMC side of the interface momentum boxes are con-
stant during the four FMG time steps. As in the CMC case, forecasts
are made only at the circled points > with the remaining momentum values
calculated by linear interpolation.
The finite difference equations corresponding to (4.1) - (4.4) at the
grid points away from the interface are
fvh <*.7)
. - )x - )y - - fuh (4.8)
-------
103
T -
- (in? A*) -
-------
104
1300 time steps. With double precision, a similar steady mass loss
-13
occurred, but the total loss after 1300 steps was only -20 x 10 %,
verifying that round-off errors were responsible for the slight loss.
The variation with time of the total energy on the uniform CMC
is illustrated in Fig. 4.3. The total energy, defined as
E=Z [g(h-h)2 + yh(u2 + v2)] (4.11)
shows an oscillation between about +; 5% of the initial energy. In
contrast to the cause of the variation in mass, the oscillations in
total energy are not caused by round-off errors, since an experiment
run in double precision showed the same variations in total energy.
Instead, these small oscillations are probably caused by the artificial
numerical treatment of the northern and southern boundaries. Whatever
the cause, the amplitude of the oscillation does not change with time.
The relatively high percentage change (up to 5%) is a consequence of
the small amount of perturbation energy present initially; in later
experiments with stronger perturbations, the percent changes in total
energy are much smaller.
4.3.2 The treatment of the interface momentum points in the
meshed grid experiments. The experiments in this section are designed
to determine the most stable way of treating the prediction points for
momentum which lie on the interface between the two meshes (Fig. 4.2).
The treatment of the prediction points for mass (Fig. 4.1) is much more
obvious, since each x-point belongs entirely within either the CMC or
the FMG.
-------
i
105
§
W
•H
S-i
CU
CX
o
O
LM
•rl
•S
I • •
S
<0
CU
C
0)
tfl
4-1
O
g
4-1
O
C
O
=•-
•H
M
I
01
0)
M
60
•H
-------
106
The basic choice in the treatment of the interfacial momentum
points is whether to treat them as FMG or CMC points. If treated as
FMG points, forecasts are made with the smaller time step of -7— at the
same time that the other FMG forecasts are made. If treated as CMC
points, the forecasts are made with the larger time step of At before the
FMG forecasts are made.
A second important parameter in determining the relative stability
of the meshed integration is the CMC time step, At. One might expect that
the computations would be linearly stable if the CFL criterion for a
staggered grid (^- < —) were satisfied on both grids. The effect of the
size of At on non-linear instability is less certain, however. Several
experiments with varying At are run to determine its effect on the
overall stability of the meshed system.
In the first set of experiments (Exp. 2 through 7) with the meshed
grid, the initial conditions consist of zero mean wind superimposed with
random perturbations of maximum amplitude 1 ms . Experiments
2 through 7 are identical except for the treatment of the momentum points
on the interface between the two meshes. Some of the important details of
the experiments are summarized in Table 4.1. Note that none of these
experiments utilizes any spatial or temporal damping devices.
The behavior of the six meshed grid experiments may be summarized
by considering the time variation of total energy, which is shown in
Fig. 4.4. In contrast to the uniform mesh experiment (Experiment 1), all
of the meshed grid experiments show an increase of total energy. However,
the rate of increase varies strongly depending on how the interface momentu
-------
W
p.
4. ,
•
e
e
VD *«» «
••v ^
•
& •*«
X e
W » % .
1 ••*.
tt .'o
a °
. tt • »
* • •
B .^ . •
B * °.
p. • • *:
•„ P. . •
* X . «•
• o3 . •
• w • »
_ B « •
B % '• ~ '•.
*B X \"% J
rt *. V •
& ' °n J .' ^'
w C * •
v.^ B - : . .
^ &B ' ** : &
^ n < • td
^ B %'. "
^ ^ 0 f •
N a < •.
N B » .*
\ *• --. i
» B w . »
< « * • •
< '. 8 '
\ * « : i
B - « • •
v " * : 5
r- V. V * . J
««*,«-• N 0B "-;• o° '
i 8:8 8 88 N °o «: ".
^HtHgE \o i?;
1 1 1 1 1 i ' •• \ • i
1 1 1 1 I \ V H >
s e B 5 B . x ^ ^. {
-'/ £ 1
ID >k *
1 u * •
1 P° *'.;
NDD 5
* ... v°o ,-:
^ - » a **•
^ ^« *•
4 Q j..
* - VD i
< a 3
•>«*,!
^<>:
k^.
*S. * . < °
a « . .
NO «: .<
a *4
II' ' ' ' f°\-*
2 S SI'SS Sg
*4 *~
». E
o g
i8 !>
•: < 5
Eff H
u -P
107
o
§
«**
•-•
r--
0 1
O fsl
O ^
i 0
p.
M
KS
w
.s
p^
60
o &-1
o 0)
O ,-t
03 C
QJ
rH
Cfl
4J
^ 0
tn 4J
0 ^
v O
o «
g S C
10 S o
H -H
4J
nj
•rt
(-1
rt
0)
o B
• g -H
•» H
-*
--d*
0)
(^
o 3
e oo
S -H
PM
'•
«
-------
108
in
43
Cd
H
co
4-1
G
5J
I
•H
M
cu
cv
t-1*
X
W
M-l
0
fr
rt
s
9
cn
C
•rl
CO
44
M
1
<*!
(3
O
•rl
4-1
cd
rl
60 ^N
CU CO
1 | fK
4_l ^
(3 CU
M 4-1
CO
MH O
o cy a
0 E
43-9
4-J 4-1
60'^
(3
cu
iJ
CO
a
rH 0
Cd 'rl
, I _L I
•rt 4J
4-1 -rl
•H -a
a i3
H 0
CJ
CO
4-1
13
MH rH -rl
O Cd O
•rl PM
4J CJ
13 cd 0
cu MH 3
0 rl 4J
4-1 CU (3
Cd 4-1 CU
CU (3 0
MHO
H S
as
W
CJ /-N
a co
0^
4-J
<
cj
O 0
W 4?
< ^
•
6
w
-a
CO -H
•H M
CO 60
flj
43 T3
0)
43 43
co co
•rl CU
H 0
43
Cd 43
4J 4-1
CO -rl
CU >
0 13
•U O
CO
O -rl
g JH 4J
ft CU
O O 'rl
MH CU
•rl M ft
C 0 X
D MH CU
O
O
m
H
CU
CO
•rl
o
+ 1 a
T^ 1
000
•a
II II C
cd
3 > M
1
O
o
oo
o
CM
H
0)
60
a
.3
CJ
8*
O
O
rH
*
0)
rH CO
43 ft
cd cu
4-1 4-1
CO CO
§0
>
rH H
M CU
rt 4-i
CU MH
C cd
•rl
J H
0
sfr
^
'0
O
ao
o
CXI
CM
<4H s-\
O •*
•
43 vf /-N
4J <•
>
O 60 <•
M -rl
60 b
60
•a a) ^H
•H CU fn
ft CO ^^
« ^
M W
/— s W
CU sf M-i
M • M-l o
0
CU -rl &O
43 fu O M
4J M 60
CO 60
• CU 13
|5 CO U -H
O <^ O ft
rH t-H cd
to W to (i;
0 00
o o in
-d- CM -a-
rH rH
» ? f
§o o
g §
0 00
o o o
sr -a- oo
o o o
CN CM CM
n -* m
ft
^t
?s
H
(3
•H
(3
cd
43
4-1
M
M-l
O
43
4-1
&
O
M
60
T3
•H
ft
Cd
1-1
CO
CO
cu
hJ
o
o
oo
£
O
o
o
vO
O
CM
vO
&
id
c
•H
(3
cd
43
4J
W
M-l
O
43
4-1
&
O
M
60
•X)
•H
ft
Cd
M
CO
CO
cu
I-J
o
o
oo
£
o
o
o
m
o
CM
r^
O
o
o
H
M
CU
4-1
M-l
cd
cu
CO
•H
O
(3
M-l
O
43 •
4-1 CO
& ft
0 Q)
M 4J
60 CO
CU CU
0 0
0 -H
W 4-1
O
O
CM
,— |
CU
CO
•H
O
+ IS
0 0
H O
T3
II (3
cd
3 M
O
O
O
•tf
O
CM
oo
60
13
•H
O
f3
CU
M
CU
MH
MH
•H
T3
CU T3
0 CU
•H CO
4-1 CO
CU
TJ M
M ft
§ 3
1 a
CJ
cd co
43 iH
M 0)
CU CO
rH -H
3 O
w a
O
O
CM
,_|
O
O
O
>*
O
CM
CTi
T3
CU
rH
O
M
4J
a
o
o
co
•rl
43
4-1
&
0
M
60
4-1
3
43
CU
CO
•rl
O
(3
CU
0
O
oo
o
o
0 vH
CO
(3 CO
•H cd
ft
O
a td
u
MH
13 0
•H
43 -3
4J 3
•H O
*•£
CO
cu cd
1 *
§ o
rH
U rH
a o
E MH
O
O
CM
,— |
cu
CO
•H
o
+ 1 §
o 0
CM O
13
II 13
cd
3 M
CJ
o
o
o-
o
CM
rH
rH
•
4-1
G
•H
4-1
(3
0
O
-------
109
points are treated and the size of the time step. In Experiment 2, in
which the points are treated as CMC points and the CMC time step is
80, the model exhibits a very rapid increase in total energy. This in-
crease resembles a linear instability, even though the linear CFL
criterion is satisfied on both grids (^ = /gH ^- = 0.4
-------
110
4.3.3 Meshed grid experiments with a mean wind of 10 ms~ . The
experiments discussed in the previous sections investigated mainly the
effect of gravity waves passing through the meshed system. Advective
effects were negligible in the absence of a mean velocity. In this
section, three experiments (Experiments 8, 9, and 10) are run with a
mean west wind of 10 ms . This mean wind is initially in geostrophic
balance with a north-south pressure gradient. The Coriolis parameter is
-4 -1
10 s . Two measures of the behavior of the meshed systems are shown,
the percent variation of total energy and the time variation of the
maximum north-south component. Ideally, the current should remain in
approximate geostrophic balance, so that the growth of the maximum
north-south velocity is a measure of the error in the forecast.
Experiment 8 is identical to Experiment 4, except for the initial
conditions, and contains no spatial or temporal damping. Experiments
9 and 10 are identical to Experiment 8, except that an effort is made to
reduce the growth of noise and energy through the damping of high spatial
or temporal frequencies. In experiment 9, the Euler-backward scheme is
used (rather than the leapfrog) to damp the high time frequencies. In
Experiment 10, a constant horizontal eddy diffusivity coefficient of
3 2-1
2.5 x 10 m s is used on both meshes.
The results from Exp. 8-10 are summarized in Fig. 4.5 Without
any smoothing devices in space or time, the meshed system exhibits a slow,
irregular increase in total energy, reaching a value of 0.01% after
1200 FMG time steps of 10 s (3.3h). During this time the maximum v-components
-------
Ill
-» o •
V *
Sx.
«,:
i
g
/o •
°. •
0\ 4
§
1 N
.g
v>
.§
i§
-a
g
c
0)
c
o
o
o
I
i
•rl
X
cd
a
T)
Cfl
S-i
Q)
c
0)
go
n)
•H
> 00
a) .
.a &
H W
o-
0)
00
•H
-¥ H
•I
I ^~
g
-------
112
show a high frequency oscillation superimposed upon a gradual increase.
Values of 0.8 ms are reached late in the forecast.
When the Euler-backward time integration scheme is used (Exp. 9),
the high frequency time oscillation in both the total energy and the
maximum v-component is nearly eliminated. Furthermore, the total in-
creases in energy and the maximum v-component are greatly reduced,
indicating that the high time frequencies are primarily responsible for the
overall growth in Exp. 8. The disadvantage is that the Euler-backward
scheme requires nearly twice as much time.
In Exp. 10, the effect of damping the high spatial frequencies
2 2
through the addition of the terms K^V (hu) and K^V (hv) in equations (4.1)
and (4.2) is investigated. This value of K^ is rather small for the
\At
present grid size; on the FMG, for example, — r— = .001. As shown by
Fig. 4.5, damping the short wavelengths improves the behavior of the model.
The rate of growth of the total energy is slowed and the growth of the
maximum v-component is nearly eliminated. These results suggest that the
meshed grid system can be integrated for reasonably long periods of
time, even in the presence of a substantial advecting velocity, provided
some method is used to prevent the increase of energy in the short
wavelengths or high frequencies.
4.3.4 Experiment with mutually interacting grids with the fine mesh
moving through the coarse mesh. In the numerical modeling of certain
atmospheric phenomena such as hurricanes, it may be desirable to follow a
certain portion of the flow by allowing the fine mesh to move within the
coarse mesh. In this section, the fine mesh is initialized with a "cloud"
-------
113
of a passive contaminant, A, in a mean flow of 20 ms . As the cloud
is advected downwind, the FMG is allowed to move within the CMC in
order to remain centered over the cloud.
When the maximum concentration "of A travels four FMG grid points
downstream, the FMG moves one CMC grid length as illustrated in
Fig. 4.6. New values for mass and momentum on the points on the leading
edge of the FMG that now coincide with old CMC points are set equal to
the old CMC values. All other new FMG points near the leading edge
(Fig. 4.6) are obtained by a linear interpolation between the neighboring
CMC and FMG values. The data on the trailing edge of the FMG are utilized
to define the newly uncovered CMC points.
Fig. 4.7 shows the time variation of total energy and the maximum
north-south velocity component in Exp. 11. The * indicates the times at
which the FMG moves. As can be seen by a comparison of Fig. 4.7 with
4.5, the behavior of the integration on the moving grid is quite comparable
to the behavior of the stationary grid experiments. In fact, the varia-
tions of total energy and maximum v are less when the FMG is moved. This
somewhat surprising result is apparently caused by the initialization
procedure on the new FMG points each time the grid is moved. The
interpolation periodically replaces a developing noisy field on the first
four columns of FMG points with a smoothly varying field. We conclude that
completely interacting meshed systems in which the FMG moves within the CMC
introduce no special stability problems compared to integrations on
stationary grids.
-------
114
BEFORE MOVE
FMC POINTS
CMC
9
* • • ...... •
NEW FM1 POINTS
* ..... o o o •f
I • ••••OOOO
NEW CMP, POINTS • • • » • O P O O
« « « • • O O o o
• • • • • . O O 0 *^
OLD CMH POINTS-
NOW FMO POINTS
AFTER MOVE
DATA AT NEW FMO POINTS (DENOTED BYO)
OBTAINED BY LINEAR INTERPOLATION
BETWEEN OLD CMC VALUES AND FMG
VALUES .
Figure 4.6 Illustration of movement of FMG within CMC
-------
115
•H
a
o
iH
0)
O
CO
(-1
o
a
cd
-a
cfl
0)
a)
cfl
4-1
o
4-1 •
O O-
o
•H a
4J -H
(0
•H 4J
M C
Cfl QJ
> (3
O
•H O
H O
60
-------
116
4.4 Mesh grid experiments initialized with Haurwitz waves
In order to examine the effects of the moving and meshing techniques
on a different physical system, a series of experiments was initialized
with Haurwitz waves of various wavelengths. Haurwitz waves are well
suited for such a test, not only because they are analogous to waves
in the real atmosphere, but also because they have linear solutions
which can be qualitatively and quantitatively compared to the numerical
forecasts. Several of the critical questions these experiments are
designed to explore are:
(1) How do the meshing and moving techniques affect the forecast of
waves that can be adequately forecast on a uniform course mesh
grid?
(2) How do differences in numerical phase speed on the grids of
different size contaminate the forecasts?
(3) How does the nested grid structure handle the non-linear
interaction of disturbances with different wavelengths?
(4) Which of the following meshing techniques is superior?
a. The one-way interacting system in which the boundary values
of , u and v on the FMG are interpolated in time and space
from a previous independent forecast on the CMC. In this
case the grids can interact in only one direction - from the
CMC to the FMG.
b. The two-way interacting system - the boundary values of , u
and v on the FMG are computation points, and mutual interactions
between the CMC and FMG are allowed.
-------
117
Table 4.2 summarizes the Haurwitz wave experiments. The experiments
can be broken into three groups with the primary difference being the
east-west wavelength. Experiments 1 through 5 (Group A) were initialized
with a 3200 km wave (wave number 1 on the CMC) . Experiments 6 through 9
(Group B) are identical except the wavelength is only 800 km
(wave number 4 on the CMC. Experiments 10 through 12 (Group C) include
both wavelengths.
4.4.1. Initial conditions and the linear solutions. A
Haurwitz wave is a Rossby-type wave restricted to a channel of finite
width (D) and subject to the boundary conditions that v (the cross-channel
component of velocity) be a maximum at the center and vanish at both
edges of the channel. In order to derive the linear non-divergent solution
we start with the linearized equations of motion for a rotating fluid
restricted to horizontal non-divergent motion,
(4.12)
(4.13)
'
(4.14)
where cf> is the geopotential.
-------
118
Table 4.2
Summary of Meshed Experiments Initialized with. Haurwitz Wave
Experiment Wavelength of Grid
Group Number Initial Wave Structure
Number of
Time Steps and
Length of Forecast
A
B
C
1
2
3
4
5
6
7
8
9
10
3200 km
ti
ii
ii
ii
800 km
ii
M
Combined 800
and
3200 km
Uniform CMC
1
2
1
2
way
way
way
way
meshed
meshed
meshed
meshed
2
2
Uniform CMC
1
2
way
way
meshed
meshed
Uniform FMG
1
way
meshed
2
2
2
2
8
0
0
0
x 10
x 10
0
x 10
x 10
x 10
x 10
x 10
5
5
5
5
5
on FMG
5 on CMC
350/24
1400/24
1400/24
1400/24
1400/24
350/24
1400/24
1400/24
800/13
800/13
hrs
hrs
hrs
hrs
hrs
hrs
hrs
hrs
.3 hrs
.3 hrs
11
2 way meshed 2 x 10 on FMG 800/13.3 hrs
8 x 105 on CMC
-------
119
We restrict propagation of the wave to the x direction and seek
solutions of the form
^_ Xx)
u(y) v(y)
where
A is amplitude
oj is frequency
t is time,
and K is the wave number in the x direction = T—where X is the wavelength.
A. X
X
Substituting (4.15) into (4.12) - (4.14), we obtain
/s /\ /\
+ KU) u - fv + 1K<(> = 0 (4.16)
-t- KU) v + fu + = 0 (4.17)
ay
+ i == o . (4.18)
With the assumption that the phase speed C = - co/K ^ U, we can
/*> s* /\
eliminate (j> and u to get a second order differential equation for v,
-------
120
9y U - C
- K2) v- 0 , (4.19)
where $ = v~ is a constant.
R 91/9
With the simplification S = 0^-= -- K ) , the general solution
U - C
becomes
v = G e1Sy + Ce- . (4.20)
An expression for u may be found by substituting (4.20) into (4.18),
C1S iSy , C2S -iSy .. ...
u = — e J + —jjr- e J . (4.21)
Next (4.21) and (4.20) are substituted into (4.16) to yield
(4.22)
To evaluate the constants CL and C_, we apply the first boundary condition:
1) at y = 0,
9v
v is a maximum, or -r— = 0,
to (4.20) in order to obtain C.. = C .
-------
121
The second condition,
/••»
2) at y = + D/2, v = 0,
applied to (4.20) yields
cos Sy = 0 . (4.23)
For the non-trivial case C is non-zero but arbitrary and
R , 1/2
cos Sy - cos[(_ p - K ) (+ D/2)]- 0 , (4.24)
U - C
which implies that the phase speed of the wave is given by
() +
(4-25)
Now let L be the wave number in the y direction. It can be shown
that
S = 5= L (4.26)
if the wavelength in the y direction (X ) is 2D.
We can now rewrite (4.20), (4.21), and (4.22)
v = 2C cos Ly (4.27)
-------
122
2C,L
u - - 1 --— sin Ly (4.28)
/\ £* — L* 4 t TTTT \
4> = - i -H1 cosLy + 12 ((° 0KUj LC, sin Ly (4.29)
K K2 1
Now to get the initial fields we use (4.15) with t = 0. Also,
since £ is arbitrary, let C = i/2 and solve for the real parts of u,
v, and .
Re(v) = - A co's Ly sin Kx (4.30)
Re(u) = \ A sin Ly cos Kx (4.31)
Re(cf>) = |- A cos Ly cos Kx - ^ +.KU^ LA sin Ly cos Kx.(4.32)
K K2
Finally, upon substituting to = - CK in (4.32) we get
Re(4>) = f A cosLy cos Kx - -^ | - (4.33)
K (L2 + K2)K
A sin Ly cos Kx
Equations (4.30), (4.31), and (4.33) are used to determine the initial
conditions for the numerical experiments by specifying f, 8, K, L, and A.
Similarly, at any time t the complete analytic solution is
u = ^ A sin Ly cos (ait + Kx) (4.34)
v = - A cos Ly sin (uit + Kx) (4.35)
-------
123
and
g
f 2 2 L
[— A cos Ly cos Kx-(L + K ) — A sin Ly] •
K K
(A.36)
cos (cot + Kx)
This solution can be compared to the numerical forecasts to
quantitatively and qualitatively evaluate the effects of the various grid
structures, provided the non-linear terms in the numerical forecasts remain small,
4.4.2 Quantitative analysis of errors. One simple statistic
for quantitatively comparing the forecast fields to the analytic solutions
is the mean square error (MSB). Another measure of the skill of a fore-
cast is the ratio of the variance of a particular field to its MSE.
A large value of the ratio means the forecast has considerable skill.
When the ratio decreases to 1, we can say the forecast has no skill.
For example, consider the effect of phase errors on this quantitative measure
of forecast skill.
If the analytic solution of v is given by
v = A cos y sin x , (4.37)
then the variance of v when we integrate over the entire forecast domain
is
VARv - . (4.38)
-------
124
Furthermore we define
Z E (v - v )2
MSEv = X y N —
where vf is the forecast value and v is the analytic value of the v-component
I 3.
and N is the total number of grid points.
In the limit as the grid size becomes small, the MSB becomes
f f 2
(v. - v ) dx dy
•L-* — . (4.40)
dx dy
If we consider only phase errors in the x-direction in the numerical
solutions, we can express the forecast solution of v as
v = - A cos y sin (x-y) (4.41)
where y is the phase error.
Then the mean square error is
2
(-A cos y sin (x-y) + A cos y sin x) dx dy
MSB = J-« n (4.42)
V ' ' dx dy
When integrated over the entire domain, we get
,2
MSEV = |- (1 - cos Y) . (4.43)
2
The maximum MSEV (A ) occurs when the analytic and forecast waves
are 180 degrees out of phase, in which case the measure of forecast skill is
-------
125
VAR
R = ^= °'25
A ratio of 1.0 occurs when the forecast and true waves are 60° out
of phase. These examples will be useful in judging the effect of phase
errors in the numerical solutions on the meshed grids.
4.4.3 Long-wave results: Exp. 1-5. This set of experiments was
designed to test the effects of the meshing and moving techniques on a
wave that can be adequately resolved on a uniform course mesh grid.
Equations (4.30), (4.31), and (4.33) with the following parameters were
used to initialize the model :
A = 10
U = 20 mps
AXCMG = - 5
AXFMG -
AT™,_ = 240 seconds
CMG
AT—,., = 60 seconds
—5 -1
f at south end of channel = 7.0 x 10 s
-4 -1
f at north end of channel = 1.0 x 10 s
3 = 1.6 x 10"11 m"1 s""1
X = 20 Ax_..r = 3200 km
X CMCj
Xy - 26 ^ - 416°
-------
126
In the nested grid experiments (2-5) the FMG is initially positioned
along the trough axis in the center of the channel. When the trough moves
four FMG increments downstream, the FMG is moved in order to remain
nearly centered along the axis . The movement is described in
Section 4.3.4.
The numerical wave speed, C , can be computed from
, C AT
CN = KZf ARCSIJ* (~AX~ SIN KAx) » (
where C is the analytic phase speed (Thompson, 1961). Both equation (4.25)
for the analytic phase speed and (4.44) give 17.39 mps indicating that
the uniform CMG should forecast the phase speed of this wave very accurately.
Figures 4.8 and 4.9 show the initial height, velocity and vorticity
fields on the CMG for the experiments of Group A and also indicate where
the FMG is located when a meshed grid structure is appropriate. Figures
4.10 and 4.11 show the same initial fields on the FMG for experiments
2-5. Figures 4.12 - 4.13 show the results of experiment 1 (uniform CMG)
after 24 hours. All the fields are quite smooth with only a slight loss
of amplitude of the initial disturbance. The trough axis has moved nearly
9 grid points which yields a phase speed of about 16.5 ms . This value is
about 5 percent less than the analytic linear phase speed, but still quite
satisfactory considering the effects of the numerical lateral boundary
•
conditions and the non-linear terms in the forecast equations.
-------
127
463.146'.446'.74ftj. 94 6'. .'14V.. 146'.. 1464.9463.4463.£463.34 6A.O" 462.IU6 2.646? .4462 .44< 2 .446?.54t2.7462.1
W, 7. ",<, 7. )',f ft.7V.ft.446f . lit5.F4I 5.6465.5465.S46S.6465.8466.14K .4466. 74f. 7 .04ft 7. 347 7 . 54f,7 .6
77 7 7 ! r
77777/777777777*7777777','7
71 .'!', 7 1 . 7'. 71 . 1470.H47C.1 '. "1 .',;!:•<. :\ >,<.!>.
7777777777777777777777777-
,9 .5470 . 14 70. «4 71 . 1*. 71 . 7471 . 9
;vnjijj.H M i MM M Tn i ' ' IT,'
7.2497.0446 . 6 vHi^^n s . 6 4 r,4 . 44r,4. 44^, 3. ^413. 54 S t. 1441. 349 1.5'-9 3 .4404 .4 49 5 .0 4 9 5 . f^^ff^f- .7497.0497.2
<33 I 33 33 333 33 333)3 33 333 (3 33 3 33111 33VJJ
33333333131313.
"448.749».h44n.444UT2497.«497.5457.2496.9496.7496.6496.441ft./446.4497.2497.5497.8498.249H.4498.6496.7
CUTA SCJIEC PY " O.IOOf (JC, CMNIUUR IMFRVAL'IS
PI'l M Wc V'l-i t If I IV
-0.1 -0.1 -'1.1 -C.I -0.') T?.0 0.1
0.n -0.1 -O.I -0.1 - 0
I'I
-0.7 -'1.7 - l.-i -0. ) -0.1
-n. i -0.3 -o. •> -n. 7 - o
2.6 ?.« 2.8 2.1, J.a
2.7 3/0- 1X0 2.7 ?.
77T77777/
777777777
-3/0 -2.7 -2
7777777777
7777777777
?.6 2.F) ?.8 2.6 2/0
I 111 111tl
1U1I1I1U
111111111
l.H 1.4 f.1
111 1 111 11111 1
_ 1
1.3 1.4 1.4 T;V
11111111111
-0.7 -0.7 -0.5 -0.3 -0.1
9999999599959999999090
9999449999999999999994«
0.1 -0.1 -0.1 -0.1 -0.0
o.fl o.l o.l o.i o.i
O.IOOE 06, C^NTLlljR IN'FRVAL IS
0.0 -0.1 -0.1 -0.1 -0.
Figure 4.8 Initial height and vorticity field on CMC for Exp. 1-5
-------
128
6\15.^ in.8 12.7 12.4 12.7 13.H/1S.5 J17.6
TATA "iC Al FO
£i'3<.)<)OlJ9tK)f1^^^*"^^*i^^t 9 t; c ^'J'iT "3199
*) 99^99''^ i 'i^^^^ * ^"^»^C C 0 q q r^ q q 9.
0.0 -1.1
Q.O -0.4 -0.7 - I
49999,99999949 9994 9S9999
qqqgqgqqqnqq
0.7 1.0 l.l I.? 1.1 L.O 0.7 0.*
o.o o.o o.o o.o o.o o.o o.o o.o o.o o.o qlo o.o o.o o.o o.o o.o o.o o.o o.o o.
TATA Sr*L60 PV 6.1CCE Olt CONTOUR INTF B VAL IS 2 "~ ' ""
Figure 4.9 Initial u- and v-component fields on CMC for Exp. 1-5
-------
129
.1.1 IT.7 IT.', IT.^ IT.* IT.l IT.-) 17.1 IT.} 17.1 17.3 17.5 17.6 I f . •; IT.*, 17. I 17.
9.6 IS.7 1-9.7 1-J.7
Figure 4.10 Initial height and vorticity field on FMG for Exp. 1-5
-------
130
V7«
Figure 4.11 Initial u- and v-component fields on CMC for Exp. 1-5
-------
131
'I'tTPUT AFTER 150. T1MF STTPS CK
IVF V'JRT I T ITY
Coin—«FTFR 23.33HCURS
"API J> Q
0.2 0.1 0\0 -0.0 -0.1 -0.1 -0.2 -C.2 -C.2 -0.2 -0.1 -0.1 -0.0 ^•
9 j 9909999999999991,^ cqqcq99oo
1.7 I 2.4 2.9 /3.1 \2
3
11111 I 3
1.7 I 2.4 2.9l 3.1 12
inn
I.I, 12.3 2.8 2. 0.1 -/-I -0.? -0.4 -O.t> -0.7 -0. fl"-5. H -0. / -0.6 -0.3 -I
09999999 9 9 9
-------
132
P4rT tf 1/6 (**LO «XX X** *tt
11 IT.
1111lnI 11
11111 I till II1llll1 1 I
,-M.n
II111 111 111III 1111
76
P4TA <;r,AtFD FV 0.100F 01, CnNTOUR TNTFPVftl. IS
rnorrA«iT VL AFT^P 350,riMP STPPS
o.o o.d""o.o o.o o.o o.o o.d n.o o.o o.o o."6" o.o" o.o " T .0 o.o o.o o.o o.o o.o o.o
o
TS <;C»LEO it o.ioot oil rnntnup INIPPV«L I •; 7
Figure 4.13 Forecast u- and v-component fields on CMC for Exp. 1
-------
133
Experiment 1 is now used as a standard for comparing the performances
of the different meshing techniques. Because the wave is well-forecast
on the CMC, the addition of the FMG cannot greatly improve the forecast.
In fact, there is no small-scale information in the initial fields.
However, because of the discontinuity in grid spacing and the moving
of the FMG, small-scale noise will be continuously generated at the
interface and may grow and contaminate the forecasts on both grids. Since
we cannot improve on the large-scale forecast, we would hope that the
addition of the movable FMG does not lead to a deterioration of the entire
forecast.
Figures 4.14 and 4.15 show the results after 1400 AT—,., for
rM(j
Exp. 2 and 3 when no numerical smoothing devices were used. The FMG has
moved 8 times. The results from both meshing techniques look similar.
There is some small-scale noise in all fields although it is most
noticeable in the vorticity patterns. In both experiments, the phase
velocity has increased slightly to about 17 mps approaching the linear
phase speed.
Figures 4.16 and 4.17 show the results after 1400 AT from
5 2-1
Exp. 4 and 5 where a horizontal diffusion term with K = 2 x 10 m s was
added to the momentum tendencies. All the fields are noticeably smoother,
especially the vorticity. The phase speed does not change. From inspection
of the velocity and height forecasts, there is little reason to choose one
meshing technique over the other. However, the vorticity fields associated
with the one-way interacting system appear to be somewhat better.
-------
134
V80
rsr-t-ftti.it- py
-'iiumi 33 i in n in ,
nj ni33 i 3*1 j*-il
1.2 4.R ?.-. ?.>) ^.3 t.'i ?.<•> fi.n /,.« «.? } 1.0 ?.^
1333*131 3 33)J313 ^^J^^U^U^UV
- - - - - 3 13J33
.1 1.? ?.3
1 1 1
mi
1 1 1 1 1 n
U 1 llll I
33
13
3331*13 1
33333313 t
3?3333 1
IllJJ 3
1 3J33
1 333
U
1 .7
3333 *
133133
333333
Illll 3
llll 1 I
111 1 1 I 3
mn
i.
m
1 1 1
1 1 1
l.o
U
mi
inn
1.7 2.3 2.* 2.2
111
I 1
1
3.3 2.5 2
111L1 till
133
8 3.1 3.2 3.3 2.1 3.7 2.« 3.3 3.2 2.4
333 1 lit
imimiimii ui iimn
i
i m
. 7
llll
llll
J
1 1 1
1 1 1
1 1 t
1.
1 1 !
11 1
C4TA_ S_C_ai EC PV (KICK).! r_61_rnMf1(,R JMFB VAL. IS J
Figure 4.14 Forecast relative vorticity and height field in Exp. 2 (1-way,
no diffusion)
-------
135
I99999>)')'?<}99'('
V76
uiiiiiiniii
niiuiuiu
iiiiiiuiiu
llllT
111 lllllll 1 1 1 l
.2
-nil 331
11333333333 11 /^~\ /^T»3T^«~O33 1 3 3 *33
.4 3.5 3.f, 2.3 (4.ll 3.5(4.fc 5.4^4.3 3.6
13333131311113 b33/ NOn'X 3 331333
1333331 /*\13 ^"7 U llllTl1- '^31133331333
11
llll
1.1 2.2 I.*- 2.9 l^
11 in ilium urn
iimmim u um
1.5 1.7 1.0 ?.4 1.4
uinim nn inn
iiimiii mi
.5 1.9 1,3
3.4 3.T 4.2 3.1 3.0
mim mil
2.1 1.5 1.5
D»T» SC*LfH it O.IOOE 06. CONTnUR INTEIVAL IS
Figure 4.15 Forecast relative vorticity and height field in Exp. 3 (2-way,
no diffusion)
-------
136
"'/'.I 2.1 ;.* '.•< !.* ?.* '.' ?.T J~ !3 T.F. 7.7 !.7 ?. /, ?.< ,.'. J.' ^.1 ?.4
11/11 ^
X
Figure 4.16 Forecast relative vorticity and height field in Exp. 4 (1-wav
with diffusion)
-------
137
2,7 ?.<. 2.5 2.6 2.6 2.6 2.6 2.6 2.5 ? .ft 2.5 2.* 2.* ?.* 2.1
7.1 2.<. 2.5 2.6 2.7 2.7 2.7 2.T 2.6 2.*, 2.6 2.5 2.5 2.5 2.3
2.5 2.6 ?.T 2.7 2.7. 7.7 7.7 7.7 2.6 2.5 7.
2.1 2.5 7.6 2.7 2.7 2.7X2-V 2
?.3 2.<. 7.6 2.6 2.7 2.T 2.7 2
2.5 2.6 2.6 2.6 2
Figure 4.17 Forecast relative vorticity and height field in Exp. 5 (2-wav
with diffusion) *
-------
138
The statistics for Exp. 1 in Table 4.3 do show a high level of
skill - the MSE are low and R is high. As mentioned previously, the
phase speed in Exp. 1 is about 0.9 mps slower than the analytic value.
If this were the only error in the forecast, after 24 hours, the phase
lag would be between 9 and 10 degrees. Letting y = 9.5° in (4.43)
we get R = 35. Because of the lateral boundary conditions, the non-linear
term, and round-off errors, R in Exp. 1 is not equal to 35. However, it is close
enough to suggest that much of the forecast errors is caused by the
phase lag.
Both meshing techniques show a decrease of skill compared to the
uniform CMC. As a result of the differences in the numerical phase speeds
on the CMC and FMG, the solutions will tend to separate at the interface
and shortwave noise will develop and gradually contaminate the solutions -
most noticeably on the FMG. However, the statistics still indicate a
high level of skill, with the one-way technique appearing to be slightly
superior. This is especially true on the CMC because the noise generated
at the interface cannot feed back to the large scale forecast in the one-
way system.
Although the addition of diffusion (Exp. 4 and 5) does produce noticeably
smoother fields (Fig. 9 and 10), the MSE actually increase
slightly in some cases as the amplitude of the wave is reduced. This is
another indication that most of the errors in the forecast are due to
phase errors which are not reduced by any numerical smoothing devices
such as diffusion.
-------
139
Table 4.3
Measures of Error for Meshed Grid Experiments
Experiment
Number v v v u u u
1
2
3
4
5
2
3
4
5
0.91
0.91
0.94
0.95
1.07
1.57
1.52
1.59
1.58
22.2
22.2
22.1
18.2
17.0
17.3
10.6
14.1
10.8
CMC POINTS
24.3
24.3
23.5
19.2
16.0
FMG POINTS
11.0
10.6
8.9
6.8
0.52
0.52
0.56
0.67
0.62
0.226
0.292
0.291
0.248
15.7
15.7
15.8
13.1
12.1
5.7
5.7
4.6
4.2
30.5
30.5
28.4
19.6
19.6
25.4
19.5
15.8
17.1
-------
140
In summary, for a long wave that can be well resolved on a grid with
uniform size, the addition of a nested fine mesh grid will not
necessarily produce a better forecast. For example, the National
Weather Service's addition of a Limited Fine Mesh (LFM) model - no
matter how the boundaries are treated - will not improve the forecasts
of the planetary scale waves which are well-forecast by the larger scale
model. It is encouraging however, that the meshed forecast is still of
high quality despite the continuous generation of noise at the interface.
4.4.4 Short wave results: Exp. 6-8.
Experiments 6, 7, and 8 are essentially identical to Exp. 1, 4, and 5
except that the wavelength is reduced by a factor of four. The initial height
and v fields (Fig. 4.18 and 4.19) show four 800-km waves on the CMC. The FMG
contains one complete wave. The purpose of these experiments is to test the
meshing techniques under very severe conditions - when the numerical phase speed
of the wave becomes quite different on the two grids. The analytic phase
speed of an 800 km wave is 19.9 ms , nearly equal to the mean wind speed.
The numerical phase speed computed from (4.44) is 19.8 ms on the
FMG but only 14.9 ms on the CMC. Thus, there is nearly a 5 ms difference
in numerical phase speeds between the grids. After six hours, this dif-
ference will result in a 45° phase difference, and after 24 hours the waves
will be approximately 180° out of phase. These errors, of course, will
introduce very large errors at the interface and greatly contaminate the
forecasts on both grids. We should therefore expect the forecasts to rapidly
deteriorate in skill with time.
-------
141
PART '
0.0 0
t
~T.I 0
H3 4 0
7.1\ 0
III
111
~ S.4 In
111
"ill
111
~~l.! 11
11 1
11 1
9.2 b
11 1
n li
~ T.4 la
ml
ml
7.1 JO
111
X.4 0
0.0 0
(
m~~scj
' i /
0 0.0 0.0
999999999999
0 -l.l -0.7
0 -3^4 -2.1
99
0
9<
9<
0
9
9
0
9
q
0
9
0
q
9
0
9
9l
0
99
9
0
9
99
0 -
991;
5.4 -X.3
39
7.1 -4.
9
S.4 -5.
77 < '1
77 9
.2 -5.
77 • )
777 ' 1
.5 -5.'
777 9
777 )
.2 -5.
77 J
77 9
.4 -5.
7 9
>9
7.1 -4 '.
99
J 199
9 9999999
1.4 -2.1
999999999
0 -l.r -0.7
999999999999
999999999999'
0 0.0 0.0
<
FTJTTC "~ ff.ISl
j S
,0.0 0.0 0
O.T 1.1 0
?.\ 3.4 0
/mm
'.4 7. 0
mil i
mm i
'.28. 0
mil i
mi i
' .7 9. 0
mi i
mi i
•i .9 9. 0
llll 1
mi i
1 .7 9. 0
llll
mil i
'.28. 0
iimiii
iiimii
< 4 7.1 0
mill
mm
3 V 5.4 0
2.1 3.4 0
0.7 — la o
9
0.0 0.0 0
>
£ Tirr cwrou
* '
0 0.0 0.0
0 -1. 1 -0.7
0 -3.4 -2.1
0 5.4 \.3
9' 1 V<7<
9 \ 9
'J -8.4 -5.
9177 9
•. 1 77 i
C -9.2 -5.
I 77 '9
9 1 777
1 -9.5 -5.
9 777 9
9 777 ' 1
0 -9.2 -5.
9 77 ' 9
9< 77 '9
L! 8.4 -5.
' 7 9
( 7.1 -4.
y9
5
99
99
(1 -3.4 -2.1
99999S999999
999995999999
o -i.r -~a.~7
999999999999
999999999999
0 0.0 0.0
0
*TRTOi«r n
0.0 0.0 0
}
O.T~ r.i o
2.1 3.4 0
A
1./5.
/llll
mm
mini
mini
.2 1.4
mil i
1111 1
.7 9.2
llll 1
mi i
.9 "3.5
llll 1
llll 1
' . 7 9.2
llll 1
mil i
' .2 e.4
iiinii
mm
4 4 7.1
mm
1 1 1 11
2.1 3.4
9
0.0 0.0
0
0
0
0
0
0
6
0
"0
9
0
'
0 0.0 0.0
99999999999
0 -l.l -0.7
0 -3.4 -2.1
-8.4 -•>. '
77 9
77
-9.7 -5.
77 1
777 T
-9.5 -5.'
777 • 9
777 ' !
-9.2 -5.
77 ' •>
77 ' 9
-8.4 -5.
7 9
^,5
999 9999999
-3.4 -2.1
99999999999
99999999999
-f.r-o.7
99999999999
99999999999
0.0 0.0
C
3 ,
0.0 0.0 0
0.7 1.1 0
7.1 3.4 0
'-/-\ C
/I 111 11
/111 111
.4 7.1 0
Illllll
llll 1 111
.2 8.4 0
mil 11
mi 11
.7 9.7 0
llll 11
llll 11
.9 9.5 0
llll 11
llll 11
llll 11
11111 11
1.7 8.4 0
imim
iiinii
mm
until
3.\ -5.4 0
2.1 3.4 0
0.7" f.r~0
S
9
0.0 0.0 0
«
3
n 0,0 n.c
99 9 9 9999 cj 9 <;
0-1,1 -0, 1
9999999999^
0 -^.4 -2,1
f)99/ \*IQ£. 7475.1
~5999999999qq999l~
qqqqQqqq9qqqq9999999999999999999959q9S999999999Q99999999999999999999999qq9999999999999999999999
?fl 2. 948 2. 9 48 1.648 0.04 8 1.64 8 2. 94 82. 9'. 81. 6't 80. 84 B1.6482 .9482 .948 1.6480.848 1 .6482 .94 W 7.948 1. 6480.H481.6
_ _ __ .- _ imimT
VirV<85.7485.74847>*»1.7i*<<4485. 7485. 7484.~>ta3. 74BXV.85.7465 .7484 .f»t3 •JA**': 4 4 8 5 .7465.74a4/*S«J.7A**. 4
11 m 11 u 1111 imT 11111111111111 m 111 TTfTrtT1111111111111111111 irrii 111111 m 1111 m i m imi 11
1111111 M 111111 m m m MI mil miiiiiiiiiiiii i u,> 111 u 111111111111111111 ILWJJI 11111 m i m i
48a/Tl88.34>%^496.5487.5>«**B7348a.3'»ii. 2486 .54s 7.24>
-------
142
v»r
Figure 4.19 Initial height and v-component field on FMG for Exp. 6-8
-------
143
Figures 4.20 and 4.23 show the height and v-component fields
after six and twelve hours from Exp. 7 and 8. The difference in
phase speeds is obvious in the v fields, where the wave moves
faster in the middle of the FMG compared to points close to the north or
south boundaries which are affected by the slower moving wave on the CMC.
The height fields have become noisy by 12 hours and all fields continue
to deteriorate out to 24 hours. Qualitatively, there again is little
difference between the two meshing techniques.
The ratio of MSE to variance, R, for Exp. 6-8 are shown in Fig. 4.24.
There is considerable skill in all experiments on both the FMG and CMC at
6 hours, although not nearly as much as for the experiments in Group A.
By 18 hours, the difference of phase speeds and the resultant noise
i
generated at the interface has virtually destroyed all skill in the fore-
casts, and R is close to 1.
The statistics in Fig. 4.24 indicate that at least out to 12 hours,
the two-way meshing technique is somewhat better than the one-way technique.
However, this apparent superiority is not consistently shown by the
statistics for the other variables. These experiments indicate that
neither meshing technique will be satisfactory when there is a large
phase speed difference in waves which move along the interface. This difference
occurs when the wave is too short to be adequately resolved on the CMC.
One encouraging result of these experiments is that the meshing of a
grid with finer resolution in a larger domain does produce a better
solution on the FMG than on the CMC. Fig. 4.24 shows that R for Exp. 7 and 8
is consistently higher on the FMG out to 12 hours. The improvement is
observed for all variables. However, it must be noted that eventually
-------
144
Figure 4.20 Forecast v-component and height field on FMG at 6 hours for
Exp. 1 (1-way with diffusion)
-------
145
-•« ->, o
-6.0 -5.8 -5.5 -5.3 -5.d-3.i -0.
77J777777777T77777J777777
777777 77777777777777777
-*.9 t-J^Ol-
7777777777777777777777177
777777J 77777777J77777777
-•>.'. -6.0 -1.9 -•-.'. -'..1
.7^77777777/777
-^.O -5.
77 t77777
177777777
- '. . f - -i .
1111
1 I 1
Hill 111
1 I!1111 111111 111 111
111111111111111111111111111111U11 111111111 n i
1111111 i 111111 U111111111111 U I It 1 \11 1111L111
flif
^ 1 U 111 11111 I L 111 1111 11 1 t • i
111 111 11(1 111 111t 11 11 I 111 1
n u it u 111111111
111 11 11111 111111 111
I 11t11U1t 1 1 1 1 ] 1 I U
1111 1111II 111 111111
111 1
Figure 4.21 Forecast v-component and height field on FHG at 6 hours for
Exp. 8 (2-way with diffusion)
-------
146
VBY
vey
Figure 4.22 Forecast v-component and height field on FMG at 12 hours tor
Exp. 7 (1-way with diffusion)
-------
147
,:&».,.
Figure 4.23 Forecast v-component and height field on FMG at 12 hours for
Exp. 8 (2-way with diffusion)
-------
148
EXPERIMENT 6 *
(UNIFORM CMC)
10 4- EXPERIMENT 7
4 (1-WAY MESHED)
FMG POINTS O
9 T" CMC POINTS •
EXPERIMENT 8
8 + (2-WAY MESHED)
FMG POINTS A
7 J- CMC POINTS A
6 -•
5--
3--
1- •
1 1 1 1
6 HRS 12 HRS 18 HRS 24 HRS
TIME
Figure 4.24 Ratio of variance to MSE for Exp. 6-f
-------
149
the noise generated at the interface will contaminate the forecast on the
FMG to the point where the advantage of the finer resolution is lost.
This critical point in time is undoubtedly a function of the ratio of
the grid sizes, but in this particular case, it occurs at about the
period of the short wave on the CMC (13 hours). We might speculate
that in a meshed, mesoscale model, for a traveling disturbance of the
small cyclone scale that may not be well resolved on the CMC, this
critical time would be on the order of 2 or 3 days. It is encouraging
that this is probably longer than we would want to run a mesoscale
forecast model.
4.4.5 Mixed long and short wave results. The use of nested grids is
warranted when the fine mesh contains energy in shorter wavelengths than
does the coarse mesh. To test the meshing technique under these conditions, a
final set of experiments was run combining the 3200 km wave of Group A and
the 800 km wave of Group B. The long wave covers the entire forecast domain.
The short wave has an 800 km wavelength in the x direction, a 1600 km wavelength
in the y direction, and is totally contained on the FMG. The short wave
trough axis is initially located 5 FMG points upstream of the long wave
trough axis (Fig. 4.25). Since analytically and computationally the short
wave moves faster than the long wave, the waves will tend to phase during
the first 12 hours as the short wave moves into the long wave trough position.
Because of the non-linear interaction between the two waves, it is
not valid to linearly add the analytic solutions at any given time to
produce a control forecast. Instead, a uniform FMG forecast (Exp. 9) over
the entire domain was made and used as the control to compare with the
results of the nested grid experiments.
-------
150-
Ssr-t:;sr«=Ei •
;:z:rz;~Err;r c
s:;::r:::;;-; .
:zr::E:rz-r:-"
i:=l:i*E~i-=- •
:Ezi:E:=i:z:; ;
;EEiEE:±:":-^ "
\ rfHMiiilj;
i i nninmi
' i 'mm
i i illslnll-
\ I =EEJEE-EEf
i flilHltl
i i "lllillj
i ? slfiMs
= iEHH:; = = j i S E
i it:: : ;'
*. : * - ' - '
: i t :
liiiiliilltf li ' i r : : E:;EtE3r
ITis-.'HiPf-*- ~ '-. • -".=:.====??
•:;;-: ! = 5 • • • 5 i
MM
••.: \ • • • •
£ Mill
in miii
:;EEt:Ef=:j j :
::::>:::::: s s
:;::;::::: = !::: 5
""llKUlllli '*
-•-'.•:•
iiil!" fill
=; = £; s t 7 S
z:r=s i j £ i
=* i i i i i
=: s « s J *
:
:
s
:
! i
; i iliUi
i : i •::•:
! s ; •.:::-.
i -. i '•••'.•)
i I I lilli
i i Him
i i iiiili
i \W
• * ?;:::
» ? i*:i=.
! I ilij!
i Uiil!
: t :s:;:
z ! rs:s:
'
• : r :i;:i
; i 2 :x:;:
i = ; liili
i \ i \vM
c :»:;:
!. :=•••
ill!
I lilli
• s
1 '. 5
;' i'iiti \
H isisi t
5 flilj J
i m I
i lilli •"
i lilli '•
S-s'ss- ;
Hill! '
i'i'i ?
01
cfl
o
J3
CO
13
C
cd
bO
C
O
•H
PS
01
•H
M
cu
ft
S-i
o
cu
n
60
•H
FJH
-------
151
The initial geopotential and v-component fields for the
control Exp. 9 are shown in Fig. 4.25 and 4.26. The square outline
indicates the initial location of the FMG for Exp. 10 and 11. In the
one-way case (Exp. 10), the CMC is initialized only with the long wave.
Because the CMC is integrated independently of the FMG in Exp. 11, it never
senses the existence of the short wave. As the results of Group A show,
the CMC adequately forecasts the long wave. However, the boundary points
of the FMG, specified continuously by the CMC and therefore containing
no short wave information, introduce noise on the FMG.
Unlike earlier experiments, a larger diffusion coefficient on the CMC
was found to be necessary in order to control the spread of noise from
the interface of the meshed grids. K = 2 x 10 was still used on the FMG,
but K = 8 x 10 was used on the CMC.
Fig. 4.27 shows the v-field from Exp. 9 after 800 AT_.., or 13.3 h.
FMG
Fig. 4.28 and 4.29 show the v-fields on the FMG from Exp. 10 and 11
at the same time. The velocity fields show quite clearly that the short wave
has moved into the long wave trough position. The 0 ms isotach marks
the trough axis. Initially the short wave trough was 5 AX upstream of
the long wave trough. After 800 time steps, the 0 isotach runs nearly
north-south indicating the waves have phased.
Both the two-way and one-way interacting systems produce solutions
with skill, however, the two-way experiment appears qualitatively more
like the control uniform FMG experiment. In fact the v-field from Exp. 11
(two-way system) looks very similar to that of the control experiment.
-------
152
_".__- -._--_.*._;=- = :-. = = -. = = -H- = r~:
'. ~ ~- ' _\ : ~-:"'~~'~~-
0)
O
42
CD
c
cfl
W)
O
13
0)
X
•H
I :- ; r : .:|=::=:HE*=:"=
ft
X
w
C
01
a
o
a
e
o
o
I
cfl
•H
VO
CN
-*
0)
^1
3
•i|-i!-i!-
"s:':rlj":£7;:":iJs"?5"ir:t"i;T::"ij7i
::::(;::::!;:!!:!!:!:::!: u:i!:i::!>:h:c
..ii.::.>:.>.!:.;:.!;.!::!
-------
153
j;
°Zf'.-I'~~.':z~, :;'.:::& " """" *""
nr.w,:
•••-i-----
mm ..
!:{!:Et:-!i:i!:!
\-\l~\l-\l-li-l
'.'.'.;•.
-li.=!
'il':'
-!|:!i:r
iliillili
iI!
ft
X
a
e
o
1-1
43
CO
ro
4-1
fi
0)
O
CX
O
o
o
u
2
::::;;: 5 : s : a : ; ;;;:•;:::::;
--i.il.. I
: ; : j : .: (
4-1
co
CO
o
a)
CN
-------
154
e
cu
T3
cu
,C
cn
01
cd
is
ex
w
S
o
x
ro
n
01
C
o
ex
o
o
cn
tt)
o
a)
oo
CN
-*
cu
60
-------
155
a
0)
4J
to
>,
(0
T3
I
o
JJ
4-1
•fi
ro
fl
a
01
o
I
o
l
CS1
~3-
01
M
3
o O- ff> O Cf (> OiO- (71 O CMC1 OP-CTOff-ff-off- O"
{TOO- CT-IO tf O> O O* O- O
iniin | in tr» i in in | in in i in in i .n «n
mini in in Mt in in in m in ut **
-------
156
Table 4.4 shows quantitatively the same result. There is considerable
skill in both methods at 13.3 hours, certainly more than the experiments
of Group B where the short wave is very poorly forecast on the CMC. Unlike
the previous experiments, there is a clear choice as to the preferred
methods. In all cases, the two-way MSB's and R's are better than the
comparable one-way results.
Other investigators (Harrison and Elsberry, 1972 - Phillips and
Shukla, 1973) have indicated that the two-way system is the superior technique
to use in a meshed numerical model. The last group of experiments in which
there was a clear, distinct separation of the scales of motion on the two
grids tends to confirm this result. However, the results from the experiments in
which only a single scale of disturbance was present, indicate there was little
difference in the two techniques. It is also possible that as the ratio of
the grid lengths on the CMC to the FMG increases, the one-way system will
prove more viable than the two-way system, because disturbances which
are adequately resolved on the FMG will be badly aliased as they move onto
the CMC. In such a two-way system, the disturbances entering the CMC would
have to be diffused in space to prevent large numerical errors.
-------
157
Table 4.4
Integrated Measures of Errors for Mixed Haurwitz Waves - 13-1/3 hours
o
3
COINCIDENT
CMC
U
J
VAR
u
1 way 10.95
2 way 11.93
1 way 5.99
2 way 12.32
1 way 6.41
2 way 11.94
MSE
u
.89
.66
1.71
.60
1.63
.61
VAR
MSEU
12.30
18.08
3.50
20.53
3.93
19.57
VAR
V
14.32
16.44
13.99
22.65
13.20
20.53
MSE
V
1.06
.75
2.02
.95
1.98
.93
VAR
MSEV
13.51
21.92
6.93
24.62
6.67
22.08
-------
158
5.0 INVESTIGATION OF SEMI-IMPLICIT MODELS
5.1 Advantages of semi-implicit models over explicit models
Because of the small grid sizes demanded by regional-scale models, the
time step must be quite small in order to maintain the CFL criterion
cAt
Or— < 1) which is necessary for the stability of explicit models. For
an external gravity wave speed of 300 m/s and a grid length of 20 km, the
maximum possible time step is 66.7 s for an unstaggered grid. For the
cAt /"~2~
staggered grid used in the 3-D mesoscale model, the criterion is -r— < -r— ,
AX fL
which yields a maximum time step of 47 s for the 20 km grid.
In many meteorological modelling applications, as the horizontal grid
scale decreases, so does the time-scale of the dominant scale of motion.
Thus a thunderstorm model may require a 200 m grid and a time step of 0.67 s,
but the time scale of the thunderstorm is one hour, requiring roughly 5400 time
steps of calculation for the complete forecast. However, the mesoscale in-
cludes phenomena with not only short time scales (as mountain waves or squall
lines), but also longer scales associated with diurnal variation in heating
or inertial effects. Thus, it may be necessary for mesoscale models to
predict on a relatively fine mesh for 24 hours in order to model diurnal
heating effects and inertial oscillations. For these longer integrations,
it would be more desirable than ever to utilize computational methods that
can tolerate larger time steps. For example, the maximum time step of 47 s
discussed above for the 20 km regional model is much smaller than is necessary
to adequately resolve the diurnal heating cycle. If a method could be found
-------
159
which would be stable for a time step of 5 minutes, for example, the part
of the solution determined by the diurnal heating variation would likely
be very nearly the same.
Recently, Kwizak and Robert 0-971), and Gerrity et al, (1973) have utilized
so called "semi-implicit" models to allow the use of greater time steps in
atmospheric models. By writing the finite difference analogs to the
relevant dynamic equations in a form which treats the terms in the equations
that are primarily responsible for the propagation of gravity waves
implicitly rather than explicitly, larger time steps than those specified
by the CFL restriction can be used. The "rule of thumb" appears to be
that the time required to make a forecast can be reduced by a factor of
four to six through the use of "semi-implicit" techniques. One of the costs
of such a gain is an increase in the complexity of the model; another may
be a loss of accuracy or a significant change in the adjustment times of
models. The latter possibility needs further investigation with atmospheric
models, especially on the mesoscale where external and internal gravity
waves play a very important role in the adjustment of the mass and momentum
fields.
In spite of the potential drawbacks, however, the probable gain in
speed by a factor of four is attractive and we are investigating the pos-
sibility of utilizing semi-implicit modeling techniques for the mesoscale
model. In order to develop familiarity with the technique and to isolate
such problems as lateral boundary conditions and staggered horizontal grids,
-------
160
a number of one- and two-dimensional experiments have been conducted
with the "shallow fluid" equations with both explicit and implicit
models. The results from these experiments (presented in the following
sections) have been sufficiently encouraging to warrant the coding of
a 3-D semi-implicit atmospheric model to be analogous as far as pos-
sible to the explicit regional model already under development.
5.2 Comparison of one-dimensional explicit and semi-implicit shallow
fluid model
The initial comparisons of the explicit and semi-implicit (S.I.)
techniques are made with one-dimensional shallow fluid models, which are
governed by the equation of motion and continuity equation,
9hu 9huu , 3h fr. ,.
9x
In analogy with the 3-D atmospheric regional model, a staggered grid
is used with the velocity components and height defined at alternate grid
points. The staggered grid provides for a better resolution of the pres-
sure gradient force (•=—) and the horizontal mass divergence (-5—) than a
grid in which all variables are defined at all grid points.
5.2.1 Development of explicit model. With the notation:
„ . a, + 1/2 - j - 1/2
x Ax
-------
161
—x i + 1/2 i - 1/2
a = ~* ~* fs ^
the finite difference equations for the explicit model may be
written:
S* ^ - gh* h (5.5)
X X
h; = - (hu)x , (5.6)
where (hu) is defined at x. = (j - 1) Ax, j = 1, 2...J ; and h is defined
J IU3.X
at x. = (j - 1/2) Ax for j = 1, 2..., J -1. The decoupled velocity, u.,
j max 3
is obtained from (hu); according to
2(hu)
u = jr r^ r- (5.7)
^ j+1/2 j-1/2'
The finite difference equations (5.5) and (5.6) provide a complete
set of finite difference equations which may be integrated as an initial
value problem. The first time step consists of a forward rather than
centered time difference.
5.2.2 Development of semi-implicit model. The semi-implicit model is
also based upon equations (5.1) and (5.2) in accordance with the semi-
implicit form of primitive equations as used by Kwizak and Robert (1971)
and by Gerrity, McPherson, and Scolnik (1973) . The pressure gradient and
The semi-implicit scheme was also reported by Kurihara (1965) as
"partly-implicit" and by McPherson (1973) as "implicit-backward".
-------
162
divergence terms are treated implicitly; other terms (except temporal
derivatives, of course) explicitly. However, in the flux form of the
equations, the term (hu) is a combination of the advection and divergence
terms of the continuity equation. It was considered to be the divergence
of the mass-weighted momentum, hu, and as such was treated implicitly.
After linearization of the pressure gradient term by the perturbation
method, the semi- implicit finite difference forms of equations (5.1) and
(5.2) become, with
h = H + h? (5.8)
and T denoting the time step (t = tAt) ,
_ 7r _?f T_i -T" T
hu + At gH h = hu - At gh1 h
X
(5.9)
- At (uT huT ) = F
X -L
h2t + At hu2t = hT l = F0 . (5.10)
x 2
In equations (5.9) and (5.10)
(5.u)
where L is any variable.
-------
163
Equations (5.8) - (5.10) now constitute a system of semi-implicit
finite difference equations. They cannot be solved directly, as can the
explicit equations, because the left sides of (5.9) and (5.10) each con-
tain two terms which are unknown at time t = T.
Differentiation of (5.9) by x yields
- At gHh + F]_ , (5.12)
x
and the substitution of (5.12) into (5.10) gives
h2t - At2 gH h2t - At F + F = G . (5.13)
X2x J_ L-
X
Equation (5.13) is a Helmholtz type equation and can be solved for
—2t —2t
h by the relaxation method. If we define h = T and:
E(AX)2V2T (5.14)
then equation (5.13) becomes
- X2T= E (5.15)
where
T = h2t (5.16)
-------
164
>2 _ (Ax)2 ,r ,,v
A = —***"— (5.17)
AtT gH
E = - A2G . (5.18)
The unknown T may be found by relaxation. The quantity E is known at
t = TAt for each grid point. An initial guess, T°, must be made for T at
each time step. After experimentation, the best first guess was found
to be linear extrapolation of the height tendency,
T° = hT + (h1'1 + T7'1) . (5.19)
The quantity T is solved by Leibmann (simultaneous) relaxation. A
residual is calculated and used to correct the current guess at a
particular grid point, followed by an alteration of the value of T
based on the residual. The m iteration is calculated by
aR-
(5.21)
where a is the relaxation coefficient and varies between 0.5 and 2.0 in
these experiments. This iterative procedure is continued until the
maximum residual across the grid is less than some arbitrary value, e.
The choice of e influences the number of iterations (and thus the speed of
computation), and the accuracy of the resulting forecast.
-------
165
—2t —2t
Once a value of T = h is found, hu is obtained by substitution of
fj .
h into equation (5.9). Then the forecast values are easily recovered
according to
hT+1 - 2h2t - h^1 (5.22)
(hu)T+1 = 2 hu2t - hu1'1 , (5.23)
5.2.3. Initialization of models. Both one-dimensional models are
initialized in the same manner. The number of grid points is 96, the
grid size Ax is 50 km, and cyclic boundary conditions are used. The mean
fluid depth is 1000 m and since the model is meant to be isentropic, a
reduced value of gravity g* = 1/4 g, is used (Seaman, 1973) . A perturba-
s\
tion of amplitude h = 25 m and wave number one is superimposed on the fluid.
The initial velocity field is placed in linear balance with the height
field by the relation
(5'24)
With all energy initially in wave number one, numerical and physical
dispersion should be minimized during the forecast process. Both models
use explicit, forward- in- time, centered-in-space finite differences for
the initial time step.
5.2.4. Results. The explicit model was run with At = 300 sec. For
*
the above values of g and mean depth, the linear phase speed of the gravity
wave is 49.44 m/sec. (same for the S.I. model case).
-------
166
cAt
For these choices of parameters, -r — =0.3 which satisfies the linear
computational stability requirement for the centered-in-time, staggered-
in-space finite difference scheme in the explicit model, which is
Several 24-hour forecasts were performed with the S.I. model to "tune"
the model parameters, a, e, and At. The basis for comparison in these
experiments was the 24-hour explicit forecast, in which the wave moved
at a speed of 49.35 ms~ , which is 99.82% of the linear phase speed. The
forecast height profiles from two of the S.I. model forecasts (with
At = 1800 and 3600 s) are compared to the height profile from the explicit
model in Fig. 5.1. In the S.I. experiment in which At = 1800 s (six times
greater than the explicit model time step) , the forecast profiles from the
S.I. and explicit models agree very closely. For a S.I. time step of
3600 s (12 times the explicit model's At), the S.I. wave lags noticeably
the explicit wave.
Sensitivity experiments with the S.I. model revealed the dependency of
the forecast accuracy upon the choice of e. It was necessary to decrease e
as the length of the time step increased, in order to prevent the growth of
short wavelength errors in the forecast (2 Ax noise) . By reducing e from
6.27 to 0.0627, the time step could be increased from 300 sec to 1800 sec.
This resulted in a 24-hour forecast essentially free of short wavelength
noise, as shown in Fig. 5.1.
Although short wavelength noise could be successfully controlled by
reducing e as At increased, other distortions from the explicit results
increased as the time step was lengthened. With a At = 1800 sec, the
-------
167
T)
a
c o
O L-l
• ^
CO II
q
S <
d
o
CN
JJ
Cfl
CO
4-1
Cfl
CS
O
0)
CO
a.
X
0)
Cfl
§
•H
Cfl
I
0)
a
o
o
ti
o
CD
o
CJ
60
•H
-------
168
observed phase speed decreased to 49.02 m/sec or 99.16% of the analytic
phase speed. Also, the amplitude of the wave oscillated by + 1.1 m.
In spite of the above problems, the S.I. forecast using a time step
of 1800 s was still quite acceptable at 24 h. Mass remains constant within
0.007%, and the total energy change is only 0.012%. Furthermore, the
wave remains smooth and moves at nearly the linear phase speed.
Several S.I. experiments were conducted with even longer time steps,
up to 3600 s. A smooth solution was obtained by decreasing e, but other
errors increased significantly. With At equal to 3600 s, the computational
phase speed was only 47.74 m/sec, or 96.6% of the linear phase speed. Also,
the oscillation of the height of the wave crest increased to + 1.8 meters.
For these reasons, a At = 1800 s was chosen as the maximum time step
capable of giving an acceptable forecast.
In addition to accuracy, the second important aspect of the S.I. model
results is the efficiency of the computations. The efficiency or speed
depends not only on the length of the time step, but also on the number
of relaxation iterations required for a 24-hour forecast. The objective
is to use the fewest possible iterations per 24 hour forecast period that
are consistent with accurate results. The number of iterations per time
step is dependent on At, a and e. If e is too large, the resulting "noisy"
solutions required more iterations to satisfy even this coarse tolerance.
On the other hand, if e was too small, the noise was suppressed, but extra
iterations were required to achieve the fine tolerance. Efficiency was
obtained by a compromise between the two extremes.
-------
169
A third factor affecting the rate of convergence is the over-
relaxation coefficient, a. The number of iterations required per time
step was found to be highly dependent upon a, with a = 1.4 being the
best choice for the 1-D model. Fig. 5.2 shows the dependency of the
number of iterations required for a 24-hour forecast (and the average
number of iterations per time step) as a function of a for four different
first guesses.
The sensitivity tests showed that by judicious choice of e, a and
the first guess, the number of iterations per 24 hours could be reduced
to 255, or between 5 and 6 iterations per time step with At = 1800 sec.
When At = 3600 sec, the iterations per time step rose rapidly to between
17 and 18, or 421 per 24 hour forecast. Thus, the 1 hour time step is not
only of poor quality, but is less efficient as well. An efficiency ratio
may be defined as
E - TEX/TS.I. (5'25)
where T and T T are the times required to compute a 24 hour forecast
KX o. J..
with the explicit and S.I. models respectively. The maximum efficiency
obtained with acceptable accuracy with At = 1800 s was found to be 2.35.
Much of the literature, for example Kwizak and Robert (1971), indicates
that implicit models should be able to produce efficiency ratios of roughly
4.0. The efficiency ratio of the S.I. model results discussed above might
therefore appear disappointing. However, it must be remembered that models
-------
170
1200 -
LOGO -.
200 -
1. T°
2. T°
3. T°
4. T°
nT-1
— f^ L j_ i, t~1
y (.n + h
hT + f (h^-1
T^1)
i
15
6.8
1.0
1.2
1.4
1.6
Relaxation Coefficient a
Figure 5.2 Number of iterations required per time step as a function of a for
four different first guesses
-------
171
discussed in the literature are 2- and 3-dimensional. Here, the
extreme simplicity of the one-dimensional explicit model means that
relatively few computational operations are required to make a 24 hour
forecast when compared to the more complex equations and the relaxation
process in the S.I. model. While both implicit and explicit models of
two-dimensional flow are more complex and require more computations to
produce a 24 hour forecast, the relative increase is greater in the
explicit case. Thus the S.I. model should become more efficient compared
to the explicit model as the number of dimensions increases.
5.3 Comparison of two-dimensional explicit and semi-implicit models
In this section, two-dimensional explicit and S.I. solutions of
Haurwitz waves are compared. The basic equations of the explicit model
are given in Section 4; the staggered grid for these experiments is
shown in Fig. 5.3.
5.3.1 Development of the two-dimensional S.I. model. Equations (4.1)
(4.3) were the basis of the S.I. model. The pressure gradient terms,
divergence terms, and tendency terms are treated implicitly, while all
other terms were treated explicitly.
After linearization of the pressure gradient term, (4.1) - (4.3)
are written:
-2t
h 2t+ At [(hu)y + (hv)x] = hT -1 = F (5.26)
x y 4
-------
172
_ _ — 2t
hu + At g H h7 = (hu)T"1 - Atg .(h7
x x
p
- At (i^hu7 )x - At (Iryrhvx ) + At f (hv)T
—2t
hv~2t + At g H h* = (hv)1"1 - At gOT7^7 h*)T
(5.28)
™*^T •«—**»J7 • iitr ^^^^_Y T"
- At (v TiiP ) - At (^"hv ) - At f (hu) = F,
y
Equations (5.26) - (5.28) now constitute a system of semi-implicit
finite difference equations. As in the one-dimensional case, these
equations cannot be solved directly; they must be altered and then
solved by relaxation. First, equation (5.27) is differentiated by x and
equation (5.28) by y in order that the divergence term may be eliminated
in equation (5.26).
2t 2t
hu7 + At g H h77 = F^7 (5.29)
x
2t 2t
hV* + At gH h™ = FT* (5.30)
y yy 6^r
A substitution of equations (5.29) and (5.30) into (5.26) yields
r2t —rrr2t
+ F7
x
(5.31)
li + At [-At g H h77 -At g H h + F
77 x
F6 '
y
-------
173
2t
h2t - At2 g H [hg + h** ] = F4 - At (F| + F^ ) = G (5.32)
xx yy y
Equation (5.32) is a Helmholtz equation and may be solved by the
—2t
relaxation method for h , hereafter referred to as T. Then since Ax - Ay
on the two-dimensional grid, equations (5.14) and (5.15) are valid for the
2-D model as well as the 1-D model.
In order to solve for T, values of E.. must be known at all interior
x grid points (see Fig. 5.3). Then, a Liebmann (simultaneous) relaxation
is carried out given a first guess of T, according to the relations
,m = Tm _m+l m m+1
j """ Vl.j "*" "i.j+l "•" ^.j-l
(5.33)
• 2 x _m _
aR?
o
4 + X
2 22
where X = Ax /(At gH) and the relaxation coefficient, a, satisfies
0.5 < a < 2.0. The relaxation process is considered completed when the
maximum residual across the grid has become less than some arbitrary e.
The choice of £ influences the number of iterations, and thus the speed
of the computation, and the accuracy of the resulting forecast. For
-------
174
— —-X —
o
CM
m
iH ffl
O C
>s 3
u o
0
CM
CM
ON
m
»
2
ON
sr
I-J
M
ON
CM
ON i
y
JZ O
•H
a
13
•H
60
m u-
rH r-
3 >
J3 f,
0)
60
60
cd
co
4-
*
— ]
Cf)
C
O
•H
CO
cu
a
T)
CM
«s
CO
CM
CM
CO
"
O CM
O
CM
O rl
•H ct)
rH T)
O C
>. 3
O O
m
o
IN
CM
O
CM
-------
175
example, if one wishes that the maximum height difference between suc-
cessive relaxations be less than X, then by definition of T we
require
m+1 m
"
4 + A2
< X/2 (5.35)
2
for all i and j, and e must satisfy |e| < ^r .
2 t; 21
Once a value of T is found for each grid points, hu and hv are
easily obtained from equations (5.26) and (5.27). Forecast values of h,
hu, and hv at T + 1 are then recovered by use of equation (5.11).
5.3.2 Initialization. Both of the 2-D models were initialized in
the same manner. The grid consists of 20 x 15 "dot" points and 20 x 14
cross points (Fig. 5.3) with a spacing of 120 km. Cyclic conditions
are assumed to eliminate the east and west boundaries. Free-slip walls
are assumed along the north and south boundaries so that
(hu)- . = (hu),
-••» j ^» j
(hu) = (h»)14fj (5.36)
(hv)1§j = (hv)15>. = 0
The heights on the boundary row of x points are specified from the
analytic linear wave solution for the Haurwitz wave. The initial pertur-
bation of the height surface consists of a Haurwitz wave superimposed on
a north-south height gradient corresponding to a mean geostrophic wind of
_T
20 ms ,
-------
176
h1 = [A ~- cos £y - A M ( -1 ) sin £y] cos kx (5.37)
kg kg 2 _ 2
u' = A • sin £y cos kx (5.38)
v? = A cos £y sin kx , (5.39)
where h', u1, v1 are the perturbations superimposed on the mean fields.
A is an amplitude, k and £ are wave numbers in x and y directions
respectively (see Fig. 5.3), and 3 is the north-south variation of
Coriolis parameter, f. With the initialization given by (5.37) - (5.39),
the analytic linear phase speed of the Haurwitz wave is
C - U - , , . (5.40)
JT + kZ
The same three equations imply a balance between the height and velocity
fields. However because h is specified rather than calculated from the
finite difference equation along the north and south boundaries, there is
some flux of momentum and mass across the lateral boundaries, and total
mass is not conserved exactly. The errors associated with the lateral
boundary conditions also excite high frequency gravity waves, which will
be discussed later. As with the 1-D model, both 2-D explicit and S.I.
models are started with an initial forward-in-time, explicit, time step.
5.3.3 Results. The 2-D explicit model used a 300 second time step.
The linear gravity wave speed for this model with H = 4816 m and g = 9.8 m
-------
177
is 217 ms . Therefore, -r— equals 0.542 s which satisfies the stability
requirement for the 2-D staggered grid (-T— < —) .
£AX £m
The 25 hour forecast of h by the explicit model is shown in Fig. 5.A.
In this experiment, the computational phase speed was 16.48 m/sec, or about
95.4% of the analytic linear phase speed of 17.28 ms . The mass flux
across the lateral boundaries was large enough to cause some distortion
in the height and velocity fields. The distortion does not become severe
enough to destroy the wave pattern, but is evident in the velocity field
(not shown) along the southern boundary. The total mass loss is 0.095%,
which represents an average decrease of height of almost 5 meters. This
mass (height) loss is not uniformly distributed, but tends to be concentrated
in the central channel of the grid where the wave amplitude is greatest.
Inspection of the deviation of the explicit model results from the
linear solution showed that the period of the mass oscillations was be-
tween 1/2 and 2 hours, and apparently caused height errors greater than those
associated with phase speed errors. The dependency of these errors upon
the boundary conditions was demonstrated by varying the lateral boundary
conditions. Use of steady state conditions on h on the northern and
southern boundaries caused total mass fluctuations of about the same range
as those described above, but their distribution in time and space was such
that the distortions of the velocity field were about five times as great
as those associated with the analytic boundary conditions on h. On the
other hand, an explicit calculation of h of the heights along the boundaries,
eliminated the mass fluctuations and the distortion in the velocity fields.
-------
178
§
•H
O
CX
X
0)
cn
cfl
CJ
O
M-l
3
O
-------
179
However, the amplitude of the height wave increased noticeably. The
direct calculation of h on the lateral boundaries might have been the best
choice if it were possible to do the same for the semi-implicit case.
However, since lateral boundary values for height could not be calculated
directly in the S.I. model, but instead had to be specified, the analytic
boundary conditions on h were used for both models.
Based on the 1-D experiments, an optimum over-relaxation coefficient was
quickly found to be 1.4 for an e of 0.149 and a At of 1800 s (X = .094).
This value of e corresponds to a tolerance of O.lm between successive
relaxations for the height field (see eq. (5.35)). A number of experiments
with the S.I. model were carried out and the results compared to those from
the explicit model.
The oscillations in mass associated with the S.I. model were of the
same order of magnitude as those from the explicit model, reaching as much
as 0.103%. Because of the non-uniformity of the mass (height) losses
across the grid, it was desirable to compare implicit and explicit model
results by normalizing the height fields. This normalization consisted
of adding or subtracting the mean mass loss or gain from the heights for
output purposes only. The oscillations of the total mass in the two models
were not in phase because of the unequal time steps.
The 25-hour forecast of height contours from the S.I. model with a
time step of 1800 s, an a of 1.4, and an e of 0.149 was very similar
to the 25-hour explicit model forecast shown in Fig. 5.4. The contours
from the S.I. model would appear identical to those in Fig. 5.4 and are
therefore not shown. The computed phase speed of the wave in the S.I.
-------
180
model was identical to the speed in the explicit model (16.48 ms~ ).
Also, the forecast velocity fields (not shown) were very similar in
both models. Even the distortions in the velocities near the lateral
boundary were similar in both models. Fig. 5.5 and 5.6 show the "error"
fields (linear solution minus the non-linear, numerical solution ) from
the explicit and S.I. models. The error patterns are similar both in
magnitude and in distribution. The greatest mass loss in both models
occurred in advance of the ridge, while the greatest loss occurred in ad-
vance of the trough. The S.I. model showed no greater tendency to develop
high spatial frequency (2Ax) noise.
The root-mean-square (RMS) difference between the numerical fore-
casts and the analytic linear solution is another measure of the quality
of a numerical forecast. The linear solution should be a good approximation
to the non-linear solution because only one "long" wave was present
initially, and therefore non-linear effects should be small. Root-mean
square differences were computed for both models at each time step for
50 hours of forecast time. The results, shown in Fig. 5.7, indicate that
the RMS differences oscillate with periods of about 2-5 hours. These oscil-
lations reflect the gravity waves generated near the boundaries. According
to Kurihara (1965) and McPherson (1973), semi-implicit models damp short
waves (the 4 Ax wave is damped the most), without significantly altering
the longer waves. In Fig. 4.7, the oscillations of RMS differences are
well correlated during the first 24 hours, indicating that the two models
This field is not strictly an error since the analytic non-linear solutioi
is unknown. However, if the non-linear effects are small, the linear
solution will approximate the non-linear solution. In any case, the line;
solution serves as a useful basis for comparison.
-------
181
-T— 1 1 r
o
4-1
oj
I
•H
O
01
r-H
0)
•H
H-l
01
in
01
00
•H
Pu
-------
182
3
O
,C
LO
CN
O
Ol
42
tn
•o
1-1
0)
•H
U-l
•a
M
O
M
(j
W
00
•H
-------
183
\ \
4-1
ao
•H
0)
Hi
T)
I
T3
c
0)
x
0)
n
a)
a
o
a) -H
& -u
01 i-H
a) o
o ra
0) 01
m C
IH -H
QJ X!
h 4-)
n)
3 T3
cr a
tn rt
C tn
to 4-1
01 tn
e t«
I o
4-1 0)
O l-i
o o
SO
•H
fn
snu
-------
184
are generating gravity waves in the same manner. While no obvious
damping with time of the RMS differences in the S.I. model is evident,
it is easily seen that the explicit RMS differences are steadily growing.
Thus, it is possible that the damping of short wavelengths in the S.I.
model might be suppressing growth of RMS differences with time. In summary,
the solutions from the S.I. model are considered to be of nearly equal
accuracy compared to the explicit model results; the major errors in both
models are associated with the lateral boundary conditions.
Finally, it is necessary to consider the efficiency of the 2-D S.I.
model. A number of experiments with different over-relaxation coefficients
were conducted. The results are summarized in Fig. 5.8, which shows the
number of iterations per 25 hour forecast for two methods of specifying
the first guess. Values of a between 1.4 and 1.6 proved optimal. It can
also be seen in Fig. 5.8 that changing the initial guess has only a
minor effect.
An interesting result from the efficiency study was the effect of the
size of the time step on the number of iterations. With At = 300 seconds,
each time step required an average of 3.55 relaxation iterations. With
At = 1800 seconds, the average number of relaxation iterations per time
step increased to 13.6, but the total number of iterations per 25 hour
forecast decreased from 1065 to 688.
The large increase in the number of relaxations required for each
time step as At increases, is apparently related to the greater distance
travelled by the fast moving, short-wavelength gravity waves,and, to a
lesser extent, the slower moving Haurwitz wave. The first guess of T is
-------
185
1SOO -
1200 -
j 1000 -
. 800 -
600
1. T° =
2. T° = ~
1.0
1.2
1.4
..6
1.8
Relaxation Coefficient a
Figure 5.8 Number of iterations required pet time step as a function of
a for 2 different first guesses
-------
186
based on an approximately linear advection of the wave, expressed as a
continuity of the previous time steps height tendency at each grid point.
This first guess was quite good when At equalled 300 seconds, but when
At was increased to 1800s, the high frequency gravity waves travel a
significant distance between time steps. This large distance travelled
produces large changes in the heights, so that the first guess becomes a
rather bad guess. It then requires many more iterations to recover and
reduce the maximum residual to less than e. The implications are that
the speed of the convergence in the relaxation part of the S.I. forecast
would be increased if the amplitude of the high frequency gravity waves
were reduced.
A comparison of the computational speed of the explicit and most
efficient S.I. model (a = 1.5, At = 1800), described by the efficiency ratio
E (in (5.25), shows a significant increase of speed in the S.I. model.
The efficiency, E, is 3.725, which approaches the efficiency ratio of 4.0
described for two-dimensional models in the literature.
5.4 Conclusions of preliminary tests of semi-implicit models
It is apparent that the development of the one and two dimensional
semi-implicit models was reasonably successful. The forecast accuracies
of the S.I. models were comparable to the accuracies of the explicit models,
and the S.I.models were roughly two to four times faster. With this en-
couraging result and with the experience gained in the 1-D and 2-D S.I.
models, it is reasonable to begin the development of a three-dimensional
semi-implicit numerical mesoscale model.
-------
187
APPENDIX - CHAPTER 5
TRIANGULARIZATION OF THE COEFFICIENT MATRIX PRODUCED
BY 3-DIMENSIONAL IMPLICIT MODELS
Sela and Scolnik (1972) describe the transformation of a square
coefficient matrix into an equivalent lower triangular matrix. This
procedure is required for efficient solution of the family of Helmholtz
type equations which is produced by a three-dimensional implicit model. In
the two dimensional case there is one Helmholtz equation for the one layer
considered and the "coefficient matrix" is one by one. As the layers in-
crease, unknown parameters multiply and are contained in a system of
Helmholtz equations as a square matrix. From Gerrity, McPherson, and
Scolnik (1973), these are represented in the matrix equation
V2 CT + ET = I , (1)
-• -
where T is the vector of unknown parameters to be forecast and F is the
vector of known atmospheric parameters at t = T, and t = T - 1. C and E are
square coefficient matrices. The length of the vectors is determined by the
number of layers, so that there is one individual Helmholtz equation in the
system (1) for each layer.
The solution of (1) can be accomplished by simultaneously relaxing
in three dimensions, but this is very inefficient and would destroy any
possibility for an increase in efficiency using S.I. models. The use of the
following matrix algorithm effectively decouples the layers of the model by
-------
188
triangularization of the coefficient matrices. Then, relaxation can be
performed on each layer individually, roughly doubling the efficiency of
the model.
If the original matrix to undergo triangularization is C.., the process,
after Sela and Scolnik (1972) is given by
Cl = L1U1 • '
U1L1 = C2
C2 = L2U2
Uk-lLk-l ' Ck
where C, is now lower traingular. The U matrices are upper triangular
and L matrices are lower triangular. The step ...
C = L U
n n n
is the key. The code for performing this operation numerically was supplied
by NMC and the example given by Sela and Scolnik on a 10 x 10 matrix nas
been duplicated.
-------
189
As a check, the original matrix can be recovered. If
= \-l\-2- - U2U1
then
Cn = B"1 L. U, B = B"1 C. B (4)
1 k k k
and likewise
B C B
-------
190
6.0 DETERMINATION OF INITIAL DATA REQUIREMENTS
In consonance with the historical emphasis placed on the numerical
modeling of large scale motion fields, there has been much research
dealing with the predictability of atmospheric motions of this scale.
That is, given certain limitations on (1) the initial data density
and accuracy, (2) knowledge of the local forcing and (3) the ability
to utilize realistic and mathematically correct lateral boundary values,
over how long a period of time will the forecast errors remain within
some definable limit? Our objective in this section will be to in-
vestigate the effect of the first type of error, which results from
inaccuracies in the initial data, on numerical prediction models of
the regional scale atmosphere.
That it is possible to experimentally seperate the error sources
according to the three ways listed above is not completely obvious
and deserves some discussion. Consider a method of categorizing fluid
flows in terms of whether initial conditions, boundary conditions or
local forcing dominate. Initial conditions almost completely dominate
the solution early in the integration. Local forcing becomes pro-
gressively more important with time. Boundary effects, whether real
or fictitious, permeate the solution in the interior of the domain at
a rate that is directly proportional to the group velocity of the waves
that are generated at the boundary and inversely proportional to the
size of the domain. These same considerations that explain the
-------
191
dominance of one factor over the others in determining the general
solution, also may be used to explain in what sense the forecast is
most susceptible to error. For example, the impact of erroneously
specified initial data will be most severe under circumstances
where the boundary conditions and local forcing are dominated by the
effect of the initial conditions. Therefore, if the predictability
experiments are structured such that the state of the flow at any
point in time is determined almost entirely "by the initial conditions,
all unwanted error sources are esentially precluded.
Conventional predictability studies involve performing a control
integration of a numerical model and using the solution to represent
the true atmosphere. The solution at some arbitrary time is then
perturbed by a specified amount to simulate measurement or analysis
errors and a new forecast is made. The growth rate of root-mean-
square errors between the "true" atmosphere and the forecast reflects
some average measure of predictability of the atmosphere as represented
by the model. This technique has also been used to evaluate the
"compatability" among the errors in the dependent variables. That is,
it is possible for an excessive amount of error in one of the variables
used for asynoptic updating or the initialization of a model to
nullify the informational content of the others. Since there
will always be errors in the initial data used as input to numerical
models of the atmosphere applied at any scale, it is pertinent to
consider the implications of error "compatability". This concept of
-------
192
error compatability will be useful in addressing the problem of
determining initial data requirements for regional scale numerical
models. The modeling technique used, however, will differ completely
from the one discussed above, in that it will embody the concept
of stochastic-dynamic prediction (Epstein, 1969; Fleming, 1970).
6.1 Development of Stochastic-Dynamic Equations
The equations for this particular stochastic-dynamic model
describe the barotropic flow of an inviscid, incompressible, hydro-
static fluid in a rotating coordinate system. Invariance of all
forecast variables in the y (north-south) direction has been assumed.
The equations of motion and continuity equation are
i+ui-fv + 8(s + h) -0 , (t.l)
(6.2)
and
9h . 9h . 3H
,, „.
(6'3)
where s is the elevation of the lower boundary and h is the depth
of the fluid above the lower surface. Note that the mean north-
3H
south height gradient, -r- , is constant and is geostrophically
dy
balanced by a mean east-west velocity, U, where
-------
193
f 7;
The model consists of two layers, the top layer being inert and
serving only to simulate a stratosphere by exerting a buoyant effect
on the lower layer. Prior to considering the influence of the upper
layer, the total depth of the model must be scaled so that the
dynamic behavior of the model will be comparable to the behavior of
the atmosphere. This depth can be considered as the sum of a mean
depth and a variable quantity h' which is a function of time and
space. To develop a realistic magnitude for this mean depth, we
consider the atmospheric scale height, H, where
RT
H - —- , (6.5)
which equals 8 km for a surface temperature, T , of 280°C.
The effect of an upper layer of less density on the motion of
the lower layer may be accounted for in two ways. One method is to
adjust the gravitational forcing in the lower layer so that it is
effectively reduced by a factor proportional to the potential
temperature difference between the upper stratospheric layer and the
lower tropospheric layer. If we define the effective gravity by g*,
it can be shown that
g* - S e T g , (6.6)
T
where 6 and 9 refer to stratospheric and tropospheric means
S J.
-------
194
respectively. Cressman (1958) found that barotropic forecast errors
were minimized when the gravitational forcing was reduced by a
factor of 4 or greater. This factor of 4 is consistent with (6.6)
if appropriate values of 0 are substantiated. An alternative to
using g* in the equations of motion is to reduce the depth of the
lower layer by a factor proportional to the reduction in gravitational
forcing. This is justified by the fact that in the linear solution
for the phase speed of external gravity waves, acceleration of gravity
appears multiplied by the mean depth of the fluid. Application of
either method would then have the same desired effect of
decreasing the phase speed of the external gravity waves to one that
is more characteristic of the more prevalent internal waves. Since
direct use of g* would drastically affect the geostrophic relationship,
the condition for inertial instability, etc., the mean depth will be
reduced by the factor of 4. The time dependent depth can be re-
presented as
a RT
h = aH + h1 = + hf , (6.7)
O
where a is equal to 1/4.
Equations (6.1) - (6.3) will be converted to spectral form by
representing the dependent variables in terms of truncated Fourier
series. Specifically
u(x,t) - £ u (t)eimkx , (6.8)
m < Km
-------
195
v(x,t) - Z V(t)elmkx
|m| < K m
(6.9)
and
h(x,t) - Z H (t)e
m < K m
imkx
(6.10)
where m is an integer wavenumber, k = 2ir/L, L is the domain length,
i = /-I and U (t) equals [U (t)]*. The spectral forms of equations
m —m
(6.1) - (6.3) are
U - - ik
m
m-p
(m-p) U U + f V
x *' m-p p m
< K
(6.11)
- ikgm(S + H )
m m
V - - ik. ;Z (m-p) U V - f U
m
IP
| m-p
1 K
< K
p m-p m
9H .
g 6
m
(6.12)
and
H = - ik
m
m-p
Z (m-p)
1 K
< K
U - V ^
p m 3y
(6.13)
- ik Z (m-p) H U
IP
| m-p
<_ K
< K
p m-p
where
-------
196
s(x) = £ S elmkx (6.14)
|m| < K m
and
6=1 , m = 0
m
6=0 , m ^ 0
m
In the first annual report of the Select Research Group
(hereafter referred to as SRG-1) we reported on some preliminary
experiments using a slightly simpler version of this set of spectral
equations. The purpose of the experiments was to gain a familiarity
with the spectral shallow fluid equations, since they will be used
as a basis for the stochastic-dynamic model. The derivation of the
stochastic-dynamic equation was outlined in SRG-1 and will be re-
produced here for reference. For notational convenience let
A = - ikm ,
m
B = - ik(m-p)
m,p
C = - ikgm ,
m
and
Dl = - 8 97
D --
U2 3y
The spectral equations (6.11) - (6.13) become
-------
197
U=ZB UU +fV+C(S+H) , (6.15)
m p m,p p m-p m mm m
V = Z B UV - fU + D, 6 (6.16)
m p m,p p m-p m 1m
and
H=EB UH +VD0
m p m,p p m-p m 2
(6.17)
+ £ B H U
p m, p p m-p
Summation subscript limits will be unstated for convenience, but
it will be assumed that only those expansion coefficients pertaining
to wavenumbers inside the permitted region will be included in the
summation; for example,
Z implies E
IP
I m-p
1 K
< K
An application of the expected value operator (E) to (6.15) - (6.17)
yields
E(UP Vp> + f E(V
(6.18)
C [E(S ) + E(H )
mm m
E(V = $ Bm,p E -
-------
198
and
E(H ) = £ B [E(U H ) + E(H U )]
m p m,p p m-p p m-p
D E(V )
z m
Using the definition of covariance,
(6.20)
E(x2) (6.21)
and the commutability of E( ) and d( )/dt (Epstein, 1969) we
obtain
E(U ) = £ B [COV(U ,U ) + E(U ) E(U )
m p m,p p' m-p p m-p
fE(V ) + C (E(S ) + E(H )]
m m m m
E(V ) = Z B [COV(U ,V ) + E(U ) E(V )
m p m,p p' m-p p m-p
- fE(V + Dl6m
(6.22)
(6.23)
and
E(H ) = Z B [COV(U ,H ) + E(U ) E(H )
m p m,p p m-p p m-p
+ COV(H ,U ) + E(H ) E(U )] (6.24)
p m-p p m-p
D.E(V )
2 m
-------
199
We now have the prognostic equations for the first moments of
the Fourier coefficients of u, v and h. To obtain prognostic
equations for the second moments, we first obtain, from the
definition of the covariance, the relationship
E(x2)
(6.25)
E(x2)
The expression for the expected value of a triple product will also
be required and is obtained from the definition of a third moment
(T) . Specifically
E(x2) COV(X;L,x3)
COV(x1,x2) + ECx-j^) E(x2) E(x3> (6.26)
Through application of equations (6.15) - (6.17), (6.22) - (6.24)
and the definition (6.26) to the relationship (6.25), the forecast
equations for the 6 covariance matrices may be derived. They are
-------
200
COV(Ul,V.) = I Bljp [E(Up) COV(U._p,V.)
+E(U._p) COV(Up,V.) +T(Up,U._p,V.)]
+E(V._p) COV(Up.D±)
+ f COV(V V.) + C. [COV(S.,V.)
i> J ! !• J
COV(H ,V.)] - f COV(U ,U.
H^.) -Z B1)p [E(Up) COV(H.,U._p)
COV(H.,Up)]
[E(Up) COV(U.,H._p)
E(H. ) COV(U. U )
J-P i» P
COV(Di*UJ-p)
(continued on following page)
(6'27)
-------
201
+ E(U, ) COV(U.,H )] + f COV(H. V )
J~P l P J » i
C± [COV(H , ±)
COV(U ,V.) (6.28 cont'd)
W YP]
COV(D1,D.J) « E B±)p [E(Up) COV(U.,U±_p)
1-p) COV(U.,Up) +
(6.29)
E(U,
P
cov(u ,v ) + c [cov(u ,s )
J J
COV(U ,H±)] + f COV(U±,V
+ C [COV(UlfS ) + COV^.H )]
-------
202
COVV =£Bi,p EE
-------
203
and
(6.32)
COV(H.,H.) = I B.>p [E(Up) COV(H.,H._p)
+E(H._p) COV(H.,Up) +T(H.,Up,H._p)
E(H ) COV(H ,
P j
E(U. ) COV(H.,H ) + T(H.,H ,U. )
i-P J P 3 P 1-P
[E(V
E(H._p) COV(H.,Up)
E(H
+E(U._p) COV(H.,Hp) +T(H.,Hp5U._p)]
D2 COV(H.,V.) + D2 COV(H.,V.)
J J
6.2 Initialization Procedure
The covariance matrices for the complex expansion coefficients
U, V and H contain the uncertainty information about u, v and h.
To define these matrices from the knowledge of the initial data
-------
204
errors, a multiple, linear regression approach will be utilized
(Mood and Graybill, 1963). Let B be a column vector representing
the p coefficients of a set of p expansion functions that are
orthogonal in physical space, and X be a matrix whose elements
X. . represent the values of the jth function at the ith of n.
!»J
gridpoints. Let $ be some measure of the state of the atmosphere
such that
$ = BX (6.33)
where
4>1 XH X12 ••' X1P
( : ) , X = (X21 X22 ''' X2p)
A. -. A. _ ••• •"•
nl n2 np
B = (: ) ,
p
and
(6-34)
If the $, must be obtained, we first must observe $ at some set
of regularly or irregularly spaced points. The following assumptions
will be made about the observations of $ at each i point: the
observations, say Y , of fy., consist of a true value plus random
-------
205
errors; the true values of $. are the expected values of the
observations. The observation errors, e., have a mean of zero
2
and a known variance a , i.e.,
E(Yi) - $± (6.35)
and Y± = $± + £± (6.36)
or Y = $ + e (6.37)
where
£
1
e
n
It will also be assumed that the e. are independent. We now have
Yi = £i * xi + £i (6-38)
or
Y = BX + e . (6.39)
Equation (6.39) corresponds to a functional relationship
type of multiple linear regression problem. The Y and E are
random variables, X is a known physical variable (non-random) and
B represents a set of unknown parameters. Estimates of 3., which
J
/\
will be referred to as 3.» will be obtained. In matrix notation
-------
206
3
B = ( ) (6.40)
which is the value of B such that the sum of the squares
n 2
ill £i
is a minimum. We need not know the density function of Y. It
can be shown (Mood and Graybill, 1963) that
~ -IT
B = S X Y (6.41)
where
s = XTX
The Gauss-Markoff Theorem verifies that B is the best linear
unbiased estimate of B. Further calculations show that
COV(B) = a2 S'1 . (6.43)
The covariance matrix is a function of the data error (a being
independent of x in this case) as well as the form of S, which is
determined by the nature of the expansion function used and the
position of the observations.
As an example of an application of this approach, consider
evenly spaced observations of the u component of the wind velocity
-------
207
The observed values of u, representing the true value with the
superimposed error, as well as the variance of the observation error
of u, are the initial conditions represented in physical space.
Using previous notation, this corresponds to a knowledge of Y and a.
This information must be used to obtain a probability statement
about the initial conditions in Fourier space, i.e., the first and
second moments of the Fourier expansion coefficients of u. The
elements of the matrix X consist of the values of the exponential
expansion function as it is evaluated for all x corresponding to the
location of the observations and for all allowed wavenumbers. Column
vector B represents the expansion coefficients for the unknown true
/\
values of u,and B is the best estimate of B based on error contaminated
data. The expected values (first moments) of the coefficients, U.,
/\
are obtained from equation (6.41) since B = E(B), and the second
moments, in terms of the covariance matrix COV(U.,U.), are calculated
using equation (6.43).
Thus, given the location of the observations of u, v, and h
along with the respective error variance for each variable, we can
calculate COV(H.,H.), COV(U.,U.) and COV(V.,V.) from equation (6.43).
Having no information about covariances among expansion coefficients
for different variables, these are initially set equal to zero.
The specification of statistical moments of higher order than the
second will not be required since, thus far, the assumption that
they remain small and can be neglected has proved to be an acceptable
stochastic closure scheme, especially for the limited duration of time
over which we require integrations.
-------
208
6.3 Energetics of the Model
The energetics of a stochastic-dynamic model are considerably
more complicated than for a deterministic system of equations, since
the conventional modes of energy must now be partitioned because
of the uncertainty in the variables that define the energies. That
is, there is energy associated with the expected values of the
variables as well as uncertain modes of energy that are related
to the errors in the initial data.
The total energy of the analogous deterministic system of
equations is equal to the sum of the kinetic and potential energies
defined by
1 CL 2
KE = -- hu dx (6.44)
0
and
PE = f I h2dx , (6.45)
where L is the length of the domain of integration. In wavenumber
space these become
KE - H U U (6.46)
2 m n m n -m-n
and
-------
209
The stochastic-dynamic energy modes for this model are the
certain and uncertain potential energy and certain and uncertain
kinetic energy, denoted by CPE, UPE, CKE and UKE respectively.
They are defined as
CPE = |t E E(H ) E(H-m) , (6.48)
2m m
UPE - [2 COV(H-m,Hm) + COV(H ,H )] , (6.49)
t. m o o
CKE = \ I I [E(V E(V E + E(V C°V(Vn'V-m-n)
E(Un) COV(Hm,U_^n) + E(U_m_n) COV(Hm,Un)
+ E(Vn) COV(Hm,V_m_n) + E(V_m_n) COV(Hm,Vn)] .
Also, the total energy (TE) of the system is given by the sum of
the certain and uncertain energies,
TE = CPE -4- UPE + CKE + UKE (6.52)
and should be conserved. In utilizing the information in the
energy budgets, it is convenient to subtract the large amount of
potential energy, associated with the mean fluid depth, that is
unavailable for conversion to other modes of energy. Therefore,
-------
210
certain potential energy will hereafter be defined as
CPE = -^ [Z E(Hm) E(H-m) - E2(H )] , (6.53)
2. m o
so that now the total energy is given by
TE = CPE + UPE + CKE + UKE + |^ E2(H ) . (6.54)
i o
The information obtained from a study of the time variation of
these energy components will be both necessary in studying the
effect of initial data errors on mesoscale predictability as well
as helpful in verifying the proper formulation and coding of the
forecast equations.
6.4 Interpretation of Predicted Variances
The stochastic-dynamic model produces time-dependent covariance
matrices that contain the uncertainty information about the expansion
coefficients. For a physical interpretation of these results, the
point variances in physical space must be obtained. Substitution of
h = & H F (6.55)
m m m
into the definition for the variance,
VAR(h) = E(h2) - E(h) E(h) , (6.56)
provides
VAR(h) = £ £ E(H H ) F F
mn mn mn
(6-57)
- z £ E(H ) E(Hn) F F
m n m l m n
-------
211
The quantity F is the basis function for the expansion, in this
irolcx
case e . Using the definition of covariance provided by
equation (6.21), the above expression becomes simply
var(h) = E z COV(H ,H ) F F . (6.58)
v ' m n m' n' m n
Given the covariance matrix for H or any other dependent variable,
the more physically meaningful spatial variations of the standard
deviations may be obtained from equations of the form of (6.58).
6.5 Pure Gravity Wave Experiments
An experiment was performed with the stochastic-dynamic model
for conditions of pure gravity wave motion (no mean advection).
Initial conditions that produce strong nonlinear interactions were
chosen, to provide a non-trivial test of the model. Highly nonlinear
conditions also allow for a relatively rapid transfer of energy from
certain to uncertain modes. Under these conditions we expect a
relatively smooth increase in the amount of uncertain energy, meaning
the forecast is becoming less credible with time.
The initial standard error of each of the depth observations
was 200 meters. The initialization routine used 17 evenly spaced
information or observation points and allowed 16 Fourier components
(8 complex expansion coefficients). There was no uncertainty in
the wind field. The Coriolis parameter and the v wind component were
-------
212
set equal to zero, the domain size was 600 km, the depth of the
unperturbed fluid was one kilometer and the centered time differencing
method was used after a single forward time step. The initial
first moments consisted of a mass perturbation, situated in the
center of the domain, with an amplitude of one kilometer above the
unperturbed level. The initial velocity first moments were zero.
Figure 6.1 shows the initial conditions as well as the
evolution of the second moments and the depth with time. The
uncertainty in both the u and h fields grew most rapidly on the
leading edge of the propagating gravity wave. This is reasonable
in the sense that uncertainties in the wave position produce the
greatest "point uncertainties" where the gradients of the wave are
greatest. Although there is some tendency for large uncertainties
at the trailing edge of the separate waves, this large gradient in
h did not exist until the waves completely separated, so the total
uncertainty growth should probably not be comparable to that on the
forward edge of the wave. The regions of large a are
more highly correlated with the values of high a, rather than with
gradients of u. Figure 6.2 shows the changes in a and a during the
first five minutes of the forecast and reveals a smooth transfer of
uncertainty from the mass field to the wind field. Except in the
center of the domain where transfers of uncertain energy are being
affected by nonlinear interactions among the first moments, the
second moments approach an equilibrium value. Although a detailed
-------
213
200
100
-r 2.0
1^ f "*
.-,20 'w --1.5 I
10 -*-1.0
(a)
200-•
100
-r25
--1.5
10 -KL.O
(b)
-r25
200-•
100
200
100
(c)
X^ "\
I
(d)
Figure 6.1. Spatial and temporal variation of au, a^ and h for the pure
gravity wave integration. Graph (a) pertains to 0 minutes;
(b) 10 minutes; (c) 20 minutes; and (d) 30 minutes.
-------
214
200
150
100
w
e
Figure 6.2. Spatial and temporal variation of the
standard deviations of u and h during
the first five minutes of the pure
gravity wave integration. Curves are
labeled in minutes.
-------
215
discussion of the initial adjustment of a and a will be provided in
a section to follow, we point out the following in regard to this
experiment. Errors in h (a, ^ 0) cause random gravity wave
components which are reflected in the initially zero u field as
error components (a ^ 0). Uncertainty about the fluid depth
breeds uncertainty in the velocities, the mechanism for this
transfer of uncertain energy being the normal modes of energy
exchange characteristic of gravity wave generation and propagation.
These uncertain energy exchanges take place irrespective of the
prevailing first moment fields.
Figures 6.3, 6.4, and 6.5 show the energy transitions for the
certain modes, the uncertain modes and the total energies respectively.
Initially, the accelerations produced along the leading edge of the
height perturbation are reflected in the growth of the certain kinetic
energy. The uncertain component of the height observation is
responsible for the simultaneous growth of the uncertain kinetic
energy. The certain and uncertain potential energies show corresponding
decreases. The total energy is conserved to an acceptable degree
and the changes in the certain energy modes correspond to what one
might expect from two gravity waves propagating in opposite directions
on a domain with cyclic boundaries. For example, at about 30 minutes,
the two waves begin to come into phase again. The total uncertain
energy shows little change during the early period of the integration,
but soon begins to grow rapidly eventually exceeding in magnitude
the total certain energy at 43 minutes.
-------
216
o
VO
O «
4-1 CJ
CU T3
43 C3
CO >>
•H 60
^
/ N .
/
/ .v
/ ••• Vx
/ •• x
/ / Nl
/ -•' 1
/ -•• v
/.-
:f N
'\
1 f
***•* -—•• ^* "
**-i»*^ *^*
^
»•
/
/
/
i
{
1
/ ••
'
/
1
1
rt / /
CJ / '
'
>• 1
O / |
pd
2 W W ! W
w H u \ cj /
X *
a \ /
M NX
;s
& ..--• "-^
w ...-••' ^^^
« s **^-
/ ~~*
-1 . 1 1 f— i • 1 1 1 1 1 • . . 3
I
. O
1 m
m
•
LO
" CO
; o
rn
in
eg
-o
eg
ui
, — 1
O
, — 1
-in
s^
CO
CU
4J
3
a
•H
e
v_/
Cl^
P
H
M CU
u a
H CU
,_!
• n)
CO -H
CU 4J
r^ qj
O CU
S 4J
O
>> o.
00
!-i C
CU -H
C n)
0) 4-1
^
C3 CU
•H 0
cfl
4J CU >-,
M 43 00
0) 4-1 S-i
O CU
CO fi
CU -H CU
J3
4J W O
PM -H
M-l CJ 4-1
O Q)
" C
n >> -H
o oojrf
•H Vj
w cu c
OJ C -H
•H CU tfl
^ 4-1
tfl C U
> -rl CU
tfl CJ
CU 4-1
a n cu
•H CU J2
H O 4J
ro
S-l
oo
•H
-------
217
o
I
g
psJ M
W O
55 E"" *
g
H
Crf
M
y
£>
Jill 1 1
ro
/•
1 ''«
\ /;
\ / v
'x- / v
\ / «•
• / *^
I / ..«•• -^
\ / c^ ^
\ / i^^»
\ / fv .'••
V >'N
A V
A A
/ * * 1
/ \ *• ^
/ 1 •
/ V *
/ \ • ^
/ • **
/ \ ^
/ • f
1 \ l«
/ * t
/ l- fj
/ \ \
1 \ v
/ w ^ v
I ^ \ 0
I 1 V
1 \ ; :
. v'
' -"j
i '/-
• W / '. W '
i ^ .•"*v*-*.?
• -*•* ^-v
r— i I ' \ i "1 i i i
r°
,^°
•m
o
•m
m
.0
*^j
CO
O
•m
•CM
o
CN
m
I-H
o
i — i
m
o
CN i-H O
G m r*T\ TrwTmTrr
"cd 2 M
4J ED M
o cu
4J 13 C •
(3 0) >>
cu n) oo
4-1 ^i *H CU
00 cfl C
CO M 4-1 CU
•H 0) )-i
C cu C
W CU 0 -H
ED rt
H H rH 4J
ca rt j-i
•H 4-1 CU
.4-100
co c 4J a
CU CU 3
T3 4-1 CU
O O JS r-l
S a, 4-1 cd
>~, c « o
00 -H W 4->
t-i nJ o
CU 4J H CU
C >-i -C
cu a> 4-1
/-N CJ •
co (d d! >-> J2
CU -H 3 00 4J
4J nj (-4 >H
3 4-1 CU Q) >
a M js c
•H a) 4-i cu C
0 CJ O
^ c co o co
3 -H -H -H
W 4-1 J-l
§ CU W CU Cfl
M ,c! f^ a CX
^ o
4-1 •> O
O >, G
00 -H J-i
C S-i 03 O
O CU -U M-l
•H a n
4-) 0) 0) T3
03 O CU
•rl C C Tl
H iH 3 3
cd to T— t
> 4J CU O
^ i^l C^
CU QJ 4J -H
3 0
•H C co co
H 3 -H -H
**
cu
VJ
3
00
Z- <7 TI
-------
218
co
cu
•u
iH CO
CB -rl
4J
O W
4J H
CD T3
JS C
4J cfl
CO >,
•H 60
)-l
H 0)
PU 0
H CU
. TT
- v TT
bO
-------
219
6.6 A Monte Carlo Comparison
The results from the stochastic-dynamic model may be compared
with the results from a Monte Carlo procedure that also provides
error growth information about a numerical forecast. In the Monte
Carlo scheme, the deterministic, spectral, model equations (to
which the stochastic-dynamic equations degenerate if second and higher
order moments are assumed to be zero) are applied successively to
members of an ensemble of initial conditions". Each member is obtained
by using a "quasi-random" number generator that simulates the
observation errors by producing deviations to some "known" state
of the fluid. The errors are normally distributed, with a given
variance, and are applied at each point. If there are p distinct
grid points and n Monte Carlo trials, there are n x p random,
independent errors generated. This method does not suffer from
the requirement of an artificial stochastic closure scheme and has
an unlimited potential accuracy as the number of trials is increased.
It is thus a desirable method of verifying solutions obtained from
the integration of a set of stochastic-dynamic equations. The
large amount of machine time required by this method, however,
prevents its routine use.
A Monte Carlo solution to the gravity wave problem,
discussed in the last section, was carried out to check the overall
performance of the stochastic-dynamic model, both in terms of the
-------
220
correctness of the machine coding as well as the appropriateness
of the stochastic closure scheme. The deterministic, spectral
model used in the Monte Carlo runs did, out of necessity, have
machine coding that was not structured identically to that used
in the stochastic-dynamic model. The possibility exists therefore
that differential roundoff and conversion errors may have produced
some unknown degree of difference between the solutions from the
two methods. All other conditions were held identical between the
two models and any significant differences in the solutions are
probably due to the closure scheme utilized.
It is possible to compare the two solutions in terms of the
energetics as well as in terms of the actual spatial variation of the
uncertainties. The latter method was chosen because it seemed more
demanding in the sense that compensating errors could be hidden in
the spatially integrated energy modes. Figure 6.6 shows the Monte
Carlo solution for a and a after 20 minutes of integration time.
The number of trials was 500. Figure 6.7 provides the same information
but corresponds to 30 minutes of integration time. The true solution,
at any point in the time scale, should be symmetric about the location
of the original perturbation in h, however, the 30 minute Monte Carlo
solution for a does show some asymetry in the center of the domain.
In other experiments with 50 and 200 trials, this central region
seemed to have less asymetry and did not exhibit errors of any
greater magnitude. In most other respects, the solutions agree well,
-------
221
200 --
150 --
100
Figure 6.6. Monte Carlo solution with 500 trials
(dashed line) and the stochastic-dynamic
solution after 20 minutes of integration
for conditions of pure gravity wave
motion.
-------
222
250
200 - -\
150 --
100
25
20 -~
w
0
15 --
10
Figure 6.7. Monte Carlo solution with 500 trials (dashed
line) and the stochastic-dynamic solution after
30 minutes of integration for conditions of
pure gravity wave motion.
-------
223
especially at 20 minutes, even though there does appear to be
some difficulty in preserving the amplitude of the peak in CT, .
Fleming (1970) points out that the third moments do not affect
the energy conservation nor are they directly involved in the
conversion from the certain to uncertain energy modes. They are,
however, directly involved in the transfer between uncertain energy
components. Since the conclusions concerning the compatability or
errors in the mass and momentum variables will be strongly dependent
upon the modeled exchange between uncertain kinetic and uncertain
potential energy, the Monte Carlo control experiment should be
examined carefully to determine if the stochastic-dynamic solution
is strongly affected by this closure assumption. In this example,
the results appear to be encouraging. The initial a, was 200 meters
and (J was zero in this particular experiment. The results in
Figure 6.6 show the maximum discrepency between the Monte Carlo and
stochastic-dynamic solution for a to be only 3 or 4% of the total
change in a up to that point in the forecast. The spatially
averaged discrepency is on the order of 1 or 2% of the total change.
These percentages are approximately the same for the a solution
at 30 minutes, shown in Figure 6.7. The average difference of a
at 20 minutes, again in terms of percentage of the total change, is
about 12 percent. By 30 minutes, this percentage has increased
somewhat but it is still less than 20%.
-------
224
Considering the highly nonlinear character of the flow, it
is probable that these relatively small observed discrepencies were
due to the neglect of the third and higher moments. For flow condi-
tions more representative of the atmosphere (i.e., much smaller
amplitude waves) this source of error should be significantly reduced.
Further Monte Carlo control experiments will be performed, however, as
the model is applied to differing atmospheric flow conditions, to insure
that the closure method remains satisfactory.
6.7 Synoptic Scale Error Compatability Experiments
To visualize the process of error or uncertainty exchange
between variables, consider one member of the ensemble of forecasts
that the stochastic-dynamic method represents. Let the model,
defined by equations (6.15-6.17) be initialized with values that
represent the exact state of the fluid plus observation errors
superimposed as perturbations. After the integration begins, the
variables attempt to adjust to a state of dynamic compatability.
This adjustment will occur because of true imbalances in the actual
fluid as well as because of imbalances inherent in the random
observation errors. The scale of the true imbalance is physically
determined, while the scale of the error-related imbalances is
artifically imposed by the grid point or observation point spacing.
The mode of adjustment of the mass and velocity fields to a
dynamically consistent state is scale dependent. Linear geostrophic
adjustment theory provides guidance in understanding the adjustment
-------
225
mechanisms (Washington, 1965). Briefly, small-scale imbalances
are eliminated rapidly through the adjustment of the mass field
by gravity waves. Large-scale imbalances respond slowly to
inertially dominated gravity-inertia waves that attempt to alter
the velocity field. Thus, the velocity field dominates the
adjustment process on the small scale and the mass field dominates
on the large scale. At intermediate scales, the adjustment is
more mutual. In all cases, the process is manifested by damped
inertial oscillations around the equilibrium state, the damping
rate being rapid for small-scale and slow for large-scale disturbances.
A knowledge about the relative dominance of these adjustment
mechanisms in a particular situation can provide information about
the flow of uncertain energy during the adjustment process. For
purposes of illustration, assume no errors exist in the initial
wind field, but that the mass field data have a random component
with prescribed variance. Also, assume no imbalances exist in the
true state of the system. If the adjustment to equilibrium is
mutual, an error in the observation of the mass field will be rapidly
reflected in the velocity field, producing an increase in the
uncertain kinetic energy from its initial value of zero. As the
velocities increase, they approach a condition of geostrophic
equilibrium. These geostrophic velocities represent uncertain
kinetic energy since they resulted solely from errors in the mass
field. Since the adjustment is mutual, the mass field perturbation
has a tendency to disperse and become smooth so as to be in balance
-------
226
with the initially smooth momentum field. The net effect of the
mutual adjustment is that some of the initial uncertain potential
energy is partitioned to uncertain kinetic energy related to the
geostrophic velocities and some is partitioned to the ageostrophic
velocities related to gravity waves. Transfer of uncertain energy
during the adjustment has contaminated one variable with the
uncertainty from another. If random errors are present in both
mass and velocity fields, which is generally the case, the ad-
justment process becomes somewhat more complicated, with uncertain
energies traveling both ways between the mass and velocity fields.
The condition in which neither velocity nor mass related variables
have a net effect of increasing the errors of the other variables
is defined as one of error consistency or compatability.
Thus far, the discussion of uncertain error transfer in terms
of the adjustment process has been confined to a single deterministic
forecast in order to simplify the consideration of the dynamics
involved. Since the stochastic-dynamic method represents an infinite
ensemble of such forecasts, each initialized with different errors
drawn from the same distribution, the modeled transfers between
uncertain energy modes pertain to an infinite number of adjustments
toward equilibrium. The statistics of such an ensemble of forecasts
are produced by the stochastic-dynamic model discussed in Section 6.1,
which has been used to investigate this question of error consistency
for the synoptic scale. The domain size for the initial experiments
was 10,000 km which means that the shortest allowed harmonic was 1250 km.
-------
227
The first moments of the forecast variables defined a pair of
"highs" and "lows" on the grid. A geostrophic u component of
25 m s was specified and the v component was defined by the
spectral form of the geostrophic wind equation
V = i mkH . (6.59)
m f m
The hydrostatic approximation and gas law for an incompressible
fluid can be used to yield a correspondence between depth errors
in the barotropic model and temperature errors in the atmosphere.
The relationship is
R Pt
6h = - (1 - ~)
-------
228
CNI
I
0
H
H
O
rH
s—'
I
w
u
M
H
u --
TIME (hours)
Figure 6.8. Transition of the uncertain kinetic
energy of the v velocity component for
various initial, standard errors of
the v component. The curves are labeled
in terms of the magnitude of the standard
deviation of the initial v error, in m s~
The initial standard deviation of the
equivalent temperature error was 1°C.
-------
229
equivalent temperature error was 1C0. When no uncertainty existed
in the wind field, the uncertain kinetic energy of the v component
increased as the uncertain potential energy decreased (not shown).
For a three m s standard error in v, a large amount of uncertain
energy was transferred to the potential mode. The standard error
in the wind field that corresponds to 1C0 standard error in the
temperature field is shown to be about 1.0 m s . As the standard
error is changed above or below this value of 1.0 m s , the direction
of the flux of uncertain energy between the two modes is changed.
Another way of looking at the same results is in terms of the space
average of the point values of the standard deviations. Figure 6.9
displays the change in a during the same time period that Figure
6.8 results apply. If the error in the initial wind field was
excessive, in terms of its consistency with the mass field error
~x ™1 ~~x
(e.g., an initial a of 3 m s ), the a decreased from its initial
value to the consistent value during the first 2 hours of the
integration. Similarity, if the wind field observations had an
accuracy that was higher than the level of consistency, their accuracy
was degraded by the errors in the observations of the mass field.
In the error consistency experiments thus far performed, the
results have proved to be independent of the synoptic situation
reflected by the first moments.
6.8 Summary and Plans for Future Research
When the stochastic-dynamic model was tested for conditions of
pure gravity wave motion, reasonable solutions were produced in terms
-------
230
ax
v
4 -r
3
2 4-
TIME (hours)
Figure 6.9. Transition of the standard deviation of v
in terms of its spatial average for various
initial, standard errors of the v component.
The curves are labeled in terms of the
magnitude of the standard deivation of the
initial v error, in m s~^-. The initial
standard deviation of the equivalent temperature
error was 1°C.
-------
231
of total uncertain energy Increase, the spatial variability of the
uncertainty and the flux of uncertain energy between the potential
and kinetic modes. The Monte-Carlo analog of this experiment,
which was used as a "control", indicated that the stochastic
closure scheme of dropping third and higher moments was acceptable
for these particular flow conditions and for integrations of up to
30 minutes in this case. The question of error compatability or
consistency was investigated for the synoptic scale. It was
determined that, for a standard error of 1°C in the temperature
observations, the wind field had to be determined to an accuracy
corresponding to a standard error of about 1.0 m s , in order that
neither field contaminate the other with its error. Similar experiments
will be performed for regional scale domains and flow conditions.
The Monte Carlo technique will be used periodically to verify the
adequacy of the stochastic closure technique under new experimental
conditions.
-------
232
REFERENCES
Anthes, R. A. , "Development of Asymmetries in a Three-Dimensional
Numerical Model of the Tropical Cyclone", Monthly Weather
Review, 100(6), 1972, 461-476.
Anthes, R. A., and T. T. Warner, "Prediction of Mesoscale Flows over
Complex Terrain", Research and Development Tech. Report ECOM-
5532, Atmospheric Science Laboratory, White Sands Missle Range,
New Mexico 88002, March 1974, 101 pp.
Benwell, G. R. and F. Bretherton, "A pressure oscillation in a 10-
level atmospheric model", Quarterly Journal of the Royal
Meteorological Society, ^4, (400), April 1968, 123-131.
Epstein, E. S., "Stochastic Dynamic Prediction", Tellus, 21,, 1969,
739-759.
Fleming, R. J., "Concepts and Implications of Stochastic Dynamic
Prediction", NCAR Cooperative Thesis No. 22, The University
of Michigan and Laboratory of Atmospheric Science, NCAR,
1970, 171 pp.
Gerrity, J. P., McPherson, R. D., and S. Scolnik, "A Semi-Implicit
Version of the Shuman-Hovermale Model", NOAA Technical
Memorandom NWS NMC-53, National Meteorological Center,
Suitland, Maryland, July 1973, 44 pp.
Haltiner, G. J., Numerical Weather Prediction, Wiley, New York,
1971, 317 pp.
Harrison, E., and R. Elsberry, A method for incorporating nested
finite grids in the solution of systems of geophysical equations.
J. Atmos. Sci., ^9, 1972, 1235-1245.
Hovermale, John B., 1974, Personal communication.
Koss, W. J., Numerical integration experiments with variable resolution
two-dimensional Cartesian grids using the box method.
Monthly Weather Review, j)9, 10, 1971, p. 725-738.
Kurihara, Y., On the use of implicit and Interactive methods for
the time integration of the wave equation, Monthly Weather Review,
Vol. j)3, No. 1, 1965, pp. 33-46.
Kwizak, M., and A. J. Robert, A semi-implicit scheme for grid point
atmospheric models of the primitive equations, Monthly Weather Review,
Vol. 99, No. 1, 1971, pp. 32-36.
-------
233
REFERENCES (Cont'd)
Langlois, W. R. and H. C. W. Kwok, Description of the Mintz-Arakawa
numerical general circulation model, Numerical Simulation of_
Weather and Climate, Tech. Report, No. 3, Dept. of Meteorology,
U.C.L.A., Feb. 1969, 95 pp.
McPherson, R. D., Damping properties of the implicit-backward integration
method, NOAA Office Note 81, National Meteorological Center,
Suitland, Maryland, April, 1973, 23 pp.
Mood, A. M., and F. A. Graybill, Introduction of the Theory of
Statistics, Second Edition, McGraw-Hill, New York, 1963, 443 pp.
Moretti, G., Importance of boundary conditions in the numerical
treatment of hyperbolic equations, The Physics of Fluids,
Supplement 11, 1969, 13-20.
Phillips, N. A. and A. Shukla, On the Strategy of Combining Coarse
and Fine Grid Meshes in Numerical Weather Prediction,
-J. Appl. Meteor. , 12, 1973, 763-770.
Saucer, W. J., Principles of Meteorological Analysis, University
of Chicago Press, Chicago, Illinois, 1955, 438 pp.
Seaman, N. L., "A Numerical Investigation into Effects Contributing
to Nocturnal Low Level Jet Oscillations, Master of Science
Thesis, Department of Meteorology, Pennsylvania State University,
August 1973, 86 pp.
Sela, J. and S. Scolnik, "Method for Solving Simultaneous Helmholtz
Equations", Monthly Weather Review, Vol. 100, No. 8, 1972, pp. 644-645.
Shapiro, R., "Smoothing, Filtering, and Boundary Effects",
Reviews of Geophysics and Space Physics, J3(2), 1970,
359-387.
Thompson, P. D., Numerical Weather Analysis and Prediction,
The Macmillan Co., 1961, 170 pp.
Washington, W. W., "Initialization of Primitive-Equation Models
for Numerical Weather Prediction", Ph.D. Dissertation,
The Pennsylvania State University, 1964, 84 pp.
-------
234
7-° EXPERIMENTS WITH SIMPLIFIED SECOND-MOMENT
APPROXIMATIONS FOR USE IN REGIONAL SCALE MODELS
*
Alfred K. Blackadar
7.1 Introduction
In connection with the development of the regional air-pollution
model, it has been recognized that there is a need for a very simple
method for taking approximate account of the rates of exchange of
momentum, heat, water-vapor, and other properties that result from the
action of sub-grid scale motions. The method that has been tried up
to this time has been equivalent to a K-type closure with the turbulent
energy constituting the most important parameter for the determination
of the value of K. The method was based on a Prandtl type of argument,
and since such an argument has only a very limited physical basis, it
was not possible to be very confident of its applicability to situations
that are likely to depart rather widely from the familar situations in
the surface layers.
During the past year emphasis has been placed on the use of the 2nd-
order closure schemes to derive and experiment with various sets of simplified
equations that might be incorporated into a regional scale model for the
purpose of calculating the most important fluxes of heat and momentum. As
in earlier studies, it is recognized that any abridgement of the complete
set of equations for the second moments carries with it the likelihood
of inferior results. However, as a practical matter, the demands placed
upon the most advanced modern computers by the regional scale model are
such that operational applications will require the use of simple methods
-------
235
wherever these suffice, and that if the more sophisticated methods are
needed, it may be for only a limited portion of the total space-time
domain of the forecast. In this regard, it appears quite feasible to
operate the regional model in several modes so that various schemes of
varying complexity can be called in response to particular situations.
The second moment equations on which the last year's experiments
have been based are the following (see Lumley and Khajeh-Nouri, 1973)
3u u „ 3u 3U
1 K . T_ 0 1 K. 0 / " -y
—Z h U. -5" U.U = - U.U -r- U.U. -r—- - "5 (u.U.tL )
3t j ox. i k j k 3x. j i 3x 3X i jTc
J J J J
1 f 3p , 3p \ • g / * 5T i s 5T\ /1 \
— (u, -r*— + u. -5*—) + -** (o u 9' +o u 9') (1)
po k\ ^\ e 3i k 3k i
3u
u.e' + u -Ir— u.ef = - u.e1 -^ - u.u
rv,<-*.V/ 'U f\ «.u ",v r\ -U.»-k.r\ r\ « . u. , \
3t i uj 3X< i 3 3x. i j 3x. 3x,. i j
j J J
(2)
'2 _ '2
TWT 3u.0'2
at w ' "i 3x 9' = - 2 u.9' -^ - —A- = - 2e (3)
o L- .J O'^ • 1 O" • OA * O
J J J J
In these equations, U. and u are the components of the mean and
fluctuation velocities, respectively; 9 and 8' the same for potential
temperature; 6 . the Kroenecker tensor; p the pressure fluctuation; and
IK
2
e and efl are the rates of dissipation of turbulent energy and 0' respectiv
-------
236
The limitations of the computers likely to be available for
operational use require that all reasonable simplifying assumptions be
incorporated into to the use of these equations. It is a property of
these models that the horizontal grid spacing is one to two orders of
magnitude larger than the vertical spacing. It is also fairly well
established that there exists a deficiency of energy in the wavelength
spectrum in the neighborhood of the dimension of the horizontal spacing
(20-30 km). For these reasons the horizontal sub-grid scale fluxes have
been neglected and the vertical fluxes have been treated as if the tur-
bulence is horizontally homogeneous in mean properties. It is also
assumed that the time rate of change of the second moment is small com-
pared to other terms in the second-moment equations. This assumption is
less easily justified when diurnal variations are involved, and should
be abandoned if computational limitations permit. In some cases, it may
be simpler to retain time-change terms so as to permit the equations to
be solved prognostically rather than diagnostically.
In writing equation (1) for the stresses we have omitted the
Coriolis term - 2(e..,, K. upui, + S-p ^o u u-) from the right side. Wyngaard
et al (1974) have discussed the possible importance of these terms. Further
study must be made to determine whether these terms need to be retained
in the regional-scale model.
7.2 The "Poor Man's Method"
In accordance with the simplifications proposed in the previous section,
we may obtain the following equations for the cases k = 3, i = 1,2:
-------
237
1 / _iP_ _L l2^ 2 9U 3 2 , g -HT
— (w -r^- + u v*1) =-a TT TT- uw + ^ u6'
p 3x 3z w 3z 3z 5
o o
(4)
1 / IE. -u IE-* 2 8v 8 2 , g -57
— (w -^- + v -r*1) = - a -5 -5— vw + -^ v6'
p dy dz w dz dz 7;
o 0
where a = w' .
w
In accordance with the views of Rotta (1951 a,b) the pressure-
correlation terms on the left side constitute the mechanism for the
restoration of isotropy. It is easily shown that if the turbulent motions
are non-divergent, these terms contribute at most only to the flux divergen
but do not constitute a source of turbulent energy. Accordingly, we put
(5)
~ (v I2 + w |£)
p 3z 3y
where cm is some universal constant. This parameterization is equivalent
to the first order term in Lumley's expansion if the relaxation time T is
set equal to Jl/q where £ is the length scale of the largest sub-grid scale
motions and
? 22
j - 2E - u' + v' +
It has more recently been shown by Lilly (1967) and Crow (1968) that t
pressure-correlation term also contains a portion that is proportional to
the shear. Accordingly, equation (4) takes the form
-------
238
c £ quw = -a a -r— + -** u9' - -r— uw
m m w 3z 5 °z
9 (6)
m w
The value of a for isotropic turbulence is 2/5. Lumley (1973) has
suggested a functional expansion to terms of third order of anisotropic
invariants by which the value of a can be estimated. Difficulties with
m
the use of this expression have not yet been resolved, and it appears
necessary at the present time to treat this factor as one of the quantities
to be derived empirically.
Following a procedure suggested by Mahrt (1973) we neglect the flux-
divergence in our following discussion and introduce the so-called
stress-Richardson number Rs defined by
3U
= a Rs
i-i9
A 1 w dz 1 = 1,2
which, of course, implies that the vector u.9* is parallel to the shear
vector. As we shall see, this assumption is entirely reasonable.
These assumptions result in the following equations for the stresses;
_1 0 ~IT
uw = [c (a - Rs) —] Ho ~
m ^ m q w 3z
_, a „
vw = [c (a - Rs) —] la ~-
m m q w 9z
(8)
These equations are K-type fluxes and are quite similar to Prandtl's
mixing-length result.
-------
239
Although the exact meaning of £ in (5) and (8) is not clear, it
must ultimeately be identifiable with either the grid spacing (which
defines the turbulence to be parameterized) or some length scale of
the turbulence spectrum. It is therefore relevant to investigate the order
of magnitude of the bracketed expression. At the present time, this can
only be done for isotropic turbulence, and we are mostly interested in the
case of free-convection, as this state presumably characterizes the free-
atmosphere whenever convection is occurring. Under these conditions, w
and 6' may be expected to be well correlated and we can find a reasonable
estimate of 9' by using the regression formula
w6'
w
Therefore,
w6' —
u9' = —;~~ uw
% (9)
-5T we'-
ve = —j— vw
a
w
Since uw and vw are proportional to the corresponding shear components,
the assumption underlying (7) is justified. We may now look to the surface
boundary layer for further indication as to the limit approached in free
convection. Introducing
-------
240
— 2
uw = - u.
we obtain for free-convection
w
where y is tne constant in the KEYPS equation, which is most likely close
to 9.0 (Businger, 1971). For free convection, several authors find
i/o
2- 1.8 (-V)173
Accordingly, we estimate that Rs approaches a limit of about -.16 and,
for isotropic turbulence,
2
a
-~- (a - Rs) = 0.19
2. m
q
Subsequent work indicates that c as defined here is of very similar
m
order of magnitude. We thus have reasonably good justification (subject
to the assumptions made previously, for putting
-------
241
K = c.. £q (10)
m 1 M
where c is a constant of order unity.
A detailed analysis of the heat flux will be given later. It can be
shown that, subject to the same simplifications, a K-type transfer takes
place, with 1C of the form,
(11)
It appears that under free convection c~ does approach a limit which is
a little more than twice c.. .
It may be expected that c- and c~ will be functions of stability, and
the treatment of these quantities as constants is an undesirable, but
necessary, feature of the Poor Man's Method. Removal of these restrictions
can be accomplished by going to more complicated schemes as will be shown
later.
To obtain a closure it is necessary to define £ and to obtain an
expression for q. The latter follows from (1) by putting i = k and summing
over k. Under conditions of horizontal homogeneity, this equation becomes
(12)
We finally introduce the flux-gradient expression
-------
242
w (q' + 2 -) = - c3
o
together with (10) and (11). We also assume
3
CF 1
e - -I- (13)
The result is
1 3q2 , 1 fTT 3q2 . TT 3q2 . r,30\2 . ,&fJ2 g „ 30
2 ^T + 2 tU 3x~ + V 3T] = Cl *q "T^ + (^ ' f C2 ^q 31
y
(14)
+ fc =3 <" "1
2
If we assume stationarity and neglect the turbulent flux of q , we
2
can solve this expression for q
This expression was the basis of the experiments reported during the first
year of the SRG grant.
It has been found by E. Mu'ller (1973) , in connection with a much finer
grid-mesh model employing a similar parameterization scheme, that the
omission of the flux-divergence of the energy caused the gradients of K and
K. to become unrealistically large. Whether this problem develops in the
regional scale model is not yet known. If it appears that retention of the
flux divergence is required in (14), then, as pointed out by Hinkelmann (1973),
it will be advantageous to retain also the time derivative. The integration
can then be carried out as an initial value problem with fading memory,
-------
243
rather than as a boundary-value problem that would require relaxation at
every time step.
The necessity to designate i poses the greatest conceptual difficulty
for the Poor Man's Method. When the grid-spacing is smaller than the
length scale of the energy containing turbulence, it is reasonable to
identify H as that length which separates the subgrid scales from the scales
of motion that are predicted by the model. If, on the other hand, the peak
of the subgrid scale turbulence spectrum falls at length scales smaller
than the grid spacing, & must be identified with the spectral peak, and is
probably best defined by (13). This would make Rotta's relaxation time
2 _
proportional to q /£ as assumed by Lumley. There appears to be no simple
scheme for accomplishing a realistic specif iction for H in the Poor Man's
Method. It is suggested that reasonable predictions could be made by
rewriting (14) in the following way:
1 3q2 .1 rTT 3q2 , „ 3q2 „ r ,311,2 ^ -.SV.2, . g 39
2 ^T +2 [U IT + V 3y~] = °1 £q I(3^ + (9^ ] " C2 *q 31
(16)
3 „ 3q CE 3
TT— c_ x,q ^r3 --- T— q
3z 3 M 9z Az
It is then suggested that if the Richardson number is less than zero, or
some other more appropriate value, i be made equal to the vertical grid spacing
2
Az and that otherwise it be made zero. Since negative solutions of q could
conceivably result from errors in the initial values or from parameterization
errors, it will probably be necessary to change negative values to zero
artifically. By this proposed scheme the largest precentage errors will be
confined to conditions when the fluxes are small and therefore relatively
unimportant.
-------
244
7.3 Semicomprehensive Methods
a. Approximations
Between the poor-man's method and the truly comprehensive
method of Lumley and Kahjeh-Nouri lie a number of possibilities for
more or less simplified methods for modeling exchange. Whether such
simplified methods are necessary or sufficient for the needs of the
regional model is not readily apparent, and may not be definitely
known until a truly comprehensive method has been perfected and applied
in a way that can provide comparisons. It is possible, however, to
subject semicomprehensive methods to rigorous tests in the atmospheric
surface layer where the behavior of the various turbulent fluctuations
statistics has been well observed.
Parameterizations similar to those to be described in this section
have been used by Deardorff (1972) and by Wyngaard et al. (1973) to
predict the distribution of fluctuation statistics in the planetary
boundary layer. Wyngaard, and to a lesser extent, Deardorff, have used
surface layer observations to provide values of the constants used in
the integration. Neither author appears to have tested his method in
detail within the surface layer itself, and throughout the planetary
boundary layer, the needed quantities for verification of these methods
are not generally known from observations.
The method to be discussed below is one of several that have been
tried; it is also the most comprehensive one, and the one that was most
-------
245
successful in predicting the surface layer moments. It is incomplete
in that all of the flux-divergence terms, some of which are known to be
important, have been neglected. This shortcoming has been to a large
extent removed in experiments now underway, but these have not yet been
completed.
In addition to the neglect of the flux-divergence terms, it id
assumed that the surface layer is horizontally homogeneous and, in a steady
mean state. The equations (1), (2), and (3) yield the following set,
after eliminating e with the use of the energy equation.
r 2. _i_ p_i ^.gQ, f,-,^
— [w -r*1 + u -^] = - a -~ H ^ u6 (17)
p dx dz w dz 7;
o o
_
p dz 3 dz 3 5
fc)
o
1-6' |E. .^rlH-^li (20)
p dX dz dz
P0 9 3z °w 8z + g 9 (21)
= - w6- f (22)
b. Parameterizations
As in the poor-man's method we shall use a Rotta type
parameterization of the pressure-correlation terms. However, in view of
the presence of distinctly anisotropic modes of turbulence, it is possible
-------
246
that better expressions for the relaxation time can be found.
The regional-scale model being considered is characterized by
a large horizontal grid length. The vertical grid-dimension is, by
contrast, at least an order of magnitude smaller. Since the dimensions
of the horizontal grid-scale eddies, if they exist, are so much larger
than the vertical scale, it does not seem that this scale should be
relevant in determining the relaxation time for the return to isotropy.
At the same time the velocities of the horizontal modes are less likely
to be relevant for the same reason. It seems, therefore, that a
more reasonable measure of the relaxation time relevant to the vertical
transfer should be Ho rather than £q~ , with £ being identified with
the vertical scale of the eddies, which would correspond to a small
multiple of the vertical grid separation, at most.
Accordingly, we introduce the following expressions
1 r 9? . 8p, n-l _ D^2 8U /00\
— lu TT*- + w -^-1 = c x a wu -pa -5— (23)
p oz dz m w m w oz
in which c is some constant to be determined empirically.
m
2
B cf 9U/3z is an expression for the Lilly-Crow compensation term.
m w
Crow has pictured this term as an effect of the stretching of "vortex
springs" by the deformation field of the mean flow. Each vortex filament,
supposed to be randomly oriented in an isotropic turbulence field, is
associated with a negative pressure disturbance. When the fluid is
stretched in one direction, the vortices whose axes lie in this direction
-------
247
suffer a further depression of their pressures and thus give an elastic
resistance to the deformation. The elastic response to this stretching
of "vortex springs" is, at the same time, being destroyed by the
scrambling effect of the turbulence over a time interval proportional
to &q or, in this case, perhaps better, i O . For isotropic
w
turbulence 3 is 3/5, but is value for typical surface-layer anisotropies
is not known at the present time.
In like manner, one obtains
2 w |2- = c Jf1 a (a 2 - i q2) (24)
3zmww3
The shear compensation term for this equation depends only on the
tangential components of the anisotropic stress tensor, and is too
small to be important. Lumley has recently suggested that this pressure-
correlation term contains a part that modifies the buoyancy term very
slightly.
In a similar manner, one obtains
)' I2- = c. Si1 a__ U6' (25)
0. — t-. As W
dx h w
and
(26)
-------
248
a.
The compensation term for the heat flux pressure-correlation may
be pictured as arising in the following way. Consider, in Figure 1,
a single Fourier mode of the vertical temperature distribution. At A
and A' the buoyant forces are upward while at B, B', etc. they are down-
ward. The result is to create positive pressure fluctuations at b,
b', etc. and negative fluctuations
at a, a', etc. Thus negative
temperature fluctuations are
associated with negative values of
9p/9z and vice versa. According to
Deardorff (1973) this effect was
first quantitatively evaluated for
» isotropic turbulence by Taulbee in
—~~~~~~ a private communication. For
a.
isotropic turbulence, 3 is 1/3.
w
Fig. 1 BUOYANCY COMPENSATION.
The dissipation term in (22) may be written, at least for isotropic
fields, as
,.2
(27)
where Afl is the Taylor microscale for temperature fluctuations (see
Tennekes and Lumley, 1973, p. 96). In the absence of a flux-divergence
2
of 9' , this term must be of the same order of magnitude as the gradient-
-------
249
2
production term, which is of order a 0' /£. We may thus arrive at the
w
result
<28>
where cfl is some constant to be determined empirically. Similar
reasoning suggests the relation
q2/£ (29)
where c is another constant to be derived from observations.
ti
The parameterizations described here enable us to arrive at a closed
set of six diagnostic equations for the turbulent moments, applicable
to steady-state, horizontally homogeneous situations. These equations ar<
c A" a wu = -a a + ue' (30)
m w m w dz ;r
0-1 _. ,_ 2 1 2N 2 — 3u , 4 g —TOT /m \
c H a (a - — q)=—wu-5iL + —^w'G' (31)
m ww33dz3o
^r1 a w?7" = - a 2 |£ + OL & 9'2 (33)
h w w 9z h 5
t)
-------
250
(35)
c. Application to regional models
Although the parameterization described above has not yet been
incorporated into the regional model, such incorporation appears to be
quite simple to do. The set can be simplified considerably by algebraic
substitution and reduced to a dependence on a single parameter that
can be furnished at each time step by the model.
Elimination of 0' between (33) and (34) gives
-
-wF7 = K, = -S (1 + - & )- (36)
8z ^ ch cfl a 2 9 9z
9 w
Thus a generalized exchange coefficient of the type suggested by simple
mixing length theory results from these equations. Even though the
relation between a and 39/3z has not yet been determined, it is clear
w
that with lapse conditions K, is larger than with inversion conditions.
The equation suggests the non-dimensional form
i 2s ^ i
K. I to - — (1 + S-)"1 ,„,
V w c c C' (37)
where the parameter s defined by
-------
251
* (38)
a 2 9 9z
w
is calculable at each grid point at the conclusion of any time step.
Using (30) and (32) we proceed similarly to eliminate u0' with
the result
c a i k-iS
m w ,.. s N — 2 , h N 9U , 0.
—„— (1 H ) wu = - a (1 ) -5— (39)
H c c, w c, dz
m h h
Thus the Reynolds stress is proportional to the shear of the mean motion
and we can define a momentum transfer coefficient
Fi-r (1-r1)/(1 + rr-) (40)
w m h h m
In order to obtain a , q must be eliminated between (31) and (35)
w
The result is
3c_c a 2 = (c - 2O k H2 (f
E m w m Em 8z
f
(41)
which suffices for obtaining a since k and k, are already known. The
w m h J
original set of six equations is thus reduced to a set of thee explicit
algebraic equations that can be evaluated sequentially.
-------
252
The most important difficulty in this procedure is that at the
initial time s cannot be determined unless a is given at each point.
w
To overcome this difficulty it may be expedient to retain the time
derivative and advective terms in (41) and to simply guess at the initial
a distribution. It is reasonable to expect that the state of turbulence
will lose its dependence on the initial state over a period comparable to
a 2/e.
w
If other statistics are desired they can be calculated from the
equations. For example, the Richardson number is obtainable from
Ri = (cm - 2cE) kms/[(4cE + cn) V + 3cE cj (42)
We may note that as s ->• °°, Ri approaches the limit
R1c = (Cm - 2cB> (2am ah Ch ' VKSa + 4cE> C9
(43)
+ 6CE Cm ah]
a quantity which observations indicate to be about 1/4.
d. Application to the surface layer
Although the application of the method to predictive models
appears to be quite straightforward, it makes sense first to apply it in
the surface layer for two reasons. First, the wealth of empirical
-------
253
knowledge of this layer provides an opportunity for a critical
test of the parameterization, which has so far been only partially
justified. Secondly, comparison of the derived equations with
emprical data provides values of c , cfi, c, , and c that are needed
m u n. ij
for more general application.
The success of the empirically based similarity theory suggests
that it would be useful to attempt to nondimensionalize the equations
by scaling £ by the Monin length
_ o / 2 '
9 (-wu) (44)
kg w6'
The definition (38) of s gives
where
s = r- x3 £ (45)
x = u/a (46)
u,, = - wu (47)
and
C = A/L (48)
-------
254
It is useful here to treat as the Independent parameter a quantity
defined as
(49)
We then obtain, by transformation of (37), (40), and (41), the set
and
i 2V
k = — (1 - ." ) (50)
h c k c'
c, k a - p
i h m ,.... N
k = - - • - 7: — (51)
m kc.c + p/k, v
h m h
4cE)] (52)
To complete the set it is necessary to adopt an appropriate
assumption for £. Experience in the surface layers has shown that the
largest flux-supporting eddies are proportional in size to the height.
With no loss of generality we can make the coefficient of proportionality
equal to unity since the values of c , c, , CQ, and c_ have yet to be
mat) h,
determined. Accordingly, we have from (49)
f = p/x3 (53)
-------
255
The set of four equations (50) through (53) thus permit k, , k ,
, , ,
and a /u. to be calculated as functions of z/L.
w *
The equation set (30) to (35) permits other functions to be cal-
culated as well. After some reduction one can obtain
where
r •
L* 6i h w
(55)
r> ^r (56)
m
w
and
2 u u £
o ~ 7~ t 2 ~ ko"-' ^57^
u,z CE k a w
* m w
Additional derived quantities may also be calculated, such as
A - = _ .
*m - S" K k a
* m w
-------
256
kz9l=k_!i (59)
T 3z k. a ^ '
J- . W £i l\- W
* n w
(60)
Although p is the independent parameter, there are constraints
imposed on the range of its values by the equations. The upper limit
comes when k, approaches zero, which then implies that k , u., and a
h rr ' ^ m * w
also become zero. This is the critical state. The negative limit of p
is found when x goes to zero implying u. approaching zero and a ^ 0.
* w
We associate this limit with free convection. From (50) and (52) we
thus observe
3k c c kcfi
c< m . , o
c + 4c ±pl2o-
m E h
It is easily observed that the critical state condition is identical to
the earlier condition that s -»• °°. The Richardson number associated with
this state is given by (43).
Let us denote by the subscript L the limiting value under free
convection. Thus
3k c_ c
li m
One easily finds that both k, and k have limiting values for free
convection. From (53) we have
-------
257
for free convection, a result that is in agreement with similarity theory.
-1/3
It follows immediately from (54) that cfQ/T must approach aQ(z/-L)
u * t)
where afi is a predictable coefficient, again in accordance with similarity
theory. Another prediction is that the total turbulent energy should
satisfy
l/o-o/^ n PT
(\ J-/ -* /•" \ ^-/ -' i M -^ //- / \
-PT) (~r) and o " - irr- (64)
2 kc^ v ^Lx ^L7 •""* 2 kc
u* E a
* w
It has been found from experience, however, that the turbulent energy
does not behave in this manner (which is also predicted by similarity
theory); the discrepancy is believed to be due to the vertical transport
of larger scale fluctuations into the surface layer from the free-
atmosphere above.
Under conditions of free-convection, the well-known KEYPS equation
reduces to
This result also follows from (58) provided y satisfies
V - () (66)
-------
258
The KEYPS equation has been largely replaced by an equation
= (1 - 15z/L) suggested by Businger as being more convenient to
use and shown to satisfy the data from the Kansas 1968 Field Program.
This equation gives a different limiting behavior than predicted
from the present theory. It must, however, be kept in mind that the
range of observed values of z/L was quite limited, and that over the
observed range the observations were also well fitted by the KEYPS equation
with y = 9.
Finally, we can verify Priestley's equation for the heat flux with
free-convection
and find
3
-P
kk
h = (—^_) (68)
L
For the determination of the constants c , c, , CQ, and c^ in our
mho .h
experiments, only the neutral state has been utilized. The selection of
the neutral state was made for two reasons. First, under neutral conditions
the neglect of the vertical flux-divergence of u^, a , heat, and energy can
be justified, leading one to have some confidence that the calculated
constants are universal. Secondly, the application of empirical relations
is particularly simple in this case.
-------
259
The following expressions were used. These follow from (50),
(51), (52), and (54) using well known empirical relationships for
neutral conditions.
a a
m , wx /• /• n\
Cm = IT (N (69)
"h
K c
(—) — (70)
X, N a v '
n m
Ic /[-Z^-]2 (71)
u* x* N
CF = c /[3a (-^) + 2] (72)
E m m u* N
The four constants thus depend on three observed quantities, a /u ,
w x
k /k, , and o"Q/T. for neutral conditions, plus a . The last has a value
mho* m
of 2/5 for isotropic conditions, but its value for the anisotropic
conditions in the surface layer is not currently known.
For the experiment described here, the value 2/5 was used for a . An
even better over-all fit was achieved using a value of .55. It may be
noted however, that a value of 1.0, which would be equivalent to
ignoring the vortex-spring stretching completely, led to a totally
unacceptable prediction of the effects of instability on the turbulent
statistics. The values of the four constants are invariant with respect
to a, .
h
-------
260
The following values were adopted from an analysis of the AFCRL
Kansas field experiment data published by Izumi (1971). It was assumed
that k = .40; a different value would have led to somewhat different
values of each of the constants, but would not have affected the
predicted values of the statistics in any way.
a afl Kh
(^) - 1.30 (=£) = 2.60 (—) = 1.18 (73)
* N * N m N
The constants were then found to be
c = 1.30 c = 2.76 CQ = .484 cw = .240 (74)
m n o &
Deardorff (1972) and Wyngaard et al. (1973) have evaluated constants
of a similar nature. However, it is not possible to compare these values
with either of the others. The main reason is that a rather than q was
w
used in the Rotta parameterization.
The value of a, for isotropic conditions should be 2/3; its value
for anisotropic conditions'in the surface boundary layer is not known.
In general the value 2/3 did not work very well and a value of 1.0 was used
instead. This question requires clarification, and it is noted that
Lumley (personal communication) has recently developed functional expansions
for the buoyancy compensation that are expected to be very helpful in
fixing this quantity.
-------
261
The use of these constants in the various equations leads to
quantitative predictions about the limiting behavior of various
fluctuation statistics for free-convection that can be compared with
observations. This comparison is given in Table 1.
Figures 2 to 7 show predicted functions of z/L in comparison with
observations reported by Izumi (1971) for- the AFCRL Kansas Field
Experiment. It is interesting that although $ deviates systematically
from the observed values, the calculated values of cj) are in general
indistinguishable from those calculated from the KEYPS equation
with y = 10.36 to three significant figures over the entire negative
range of z/L.
The agreement of most predictions is acceptable for unstable
2 9
conditions with the exception of -u9'/w9'. The predictions of q /a
w
are also not very realistic - a fact already discussed. The predictions
for the stable side are generally unacceptable. A part of the failure
may be due to the obviously unacceptable association of £ with z under
stable conditions. [The critical Richardson number prediction is independent
of this assumption.] Other causes of the failure may be due to changes of
a and a, as the anistropy tends toward rather extreme states.
m h VJ
-------
262
Table 1
Comparison of Predicted and Observed Constants
Column three values were reported by Businger et al. (1971) and Wyngaard
et al. (1971) for the AFCRL Kansas Field Experiments. Column 4 indicates
values reported by other sources.
Statistics
(tfw/u^)N
(a /T )
9 * N
(Km/Kh)N
R.
1C
a
w
a9
Y
h
r u'9' /ru*,2.
L~" ' \ — / J T
w'9' w L
WL
Predicted Observed
Kansas Other
1 2
1.30 1.15 - 1.4 0.7 - 1.4
2.601 2.9 -2.5
.851 .74 -1.
.48 .21 -1/4
1.82 1.9 1.9 - 2.8
1.13 .95 -1.
10.4 9. 9-15
1.5 0.6 - -2
1.12 l.O3
~2. 2.05
1. Assumed values
2. Higher value is that approached from stable side.
3. Value predicted from the automatic correlation hypothesis, eq. (9).
-------
263
N
O
I
K
UJ
CO
00
O
X
i
_i
Vs
N
O
Z
UJ
UJ
m
z
o
OC.
O
O
IT
Q.
ui
UJ —
tr
o Sg
UJ (O
Q
UJ
O
w
E
z o
o
-------
264
E oo
-©. to
a:
o
co
—. CO
UJ Z
Q <
CO *
a «i
m o
< u-
? <
CO
i y
O
tr
u.
CO
UJ
UJ
o>
to
I- £
< o
Q
UJ
Q
UJ
U. 5
o _. ^
Ul ^
2 Q O
O co O
CO °Z
E iiJ °"
o ^
10
Ul
(T
3
o
-------
265
O)
oc. —
° CO
^ <
UJ
-------
266
N
u.
Z GO
UJ to
UJ O>
UJ
tE
O
UJ
O
Q
UJ
§1
O QQ
O O
in
ui
a:
-------
267
N
o Z
Z <
< cr
o
— o
-i o:
N
z
UJ
UJ
00
o
UJ
o
o <
UJ UJ
I"
n
V)
>
< tE
Q. iii
2 0)
o m
o o
ui
tr
-------
268
(D
00
CM
N
CJ
r
CO
I'
N
2
Z
<
®
o
"S3
l UJ
jji I
JL ~
00 t/>
o 3
5*
_i _i
UJ CC
Q; j_j
U.
0 <
UJ . .
LJ
K
Q.
Z O
O p
52 <
o: >
< oc
Q. uj
o
o o
UJ
(C
10
I
-------
269
It is indeed surprising to find such good agreement between predicted
and observed free-convection statistics when one reflects on the large
effects that the neglected flux-divergence terms might be expected to
have under these conditions. These terms have been reasonably well ob-
served for the energy and vertical velocity variances. There are good
reasons to think that the flux divergence is also important in the u9'
equation. Work is currently in progress to include these terms in the
surface layer verification. At the same time, it is hoped that more
acceptable estimates of a and a can be obtained by application of Lumley's
work.
-------
270
ACKNOWLEDGEMENT
During the past year, the author was on sabbatical leave in
Mainz, Germany, and in Ris^, Denmark. A substantial portion of the
author's support during this period was provided by the Alexander von
Humboldt Foundation and the Danish Atomic Energy Commission; the remainder
was provided out of general funds by The Pennsylvania State University.
The author wishes to thank Professors K. Hinkelmann, F. Wippermann,
and F. Fiedler for many stimulating dialogues about poor-man's methods
of modeling boundary layer processes. These prepared the way for the
fruitful experiments on semiquantitative methods that were conducted over a
three-month period in Ris^ amid the encouragement of and helpful
discussions with Dr. Niels Busch and his associates.
-------
271
REFERENCES
Businger, J., J. C. Wyngaard, Y. Izumi, and E. F. Bradley, 1971: Flux-
profile relationships in the atmospheric surface layer. J^. Atm. Sci. ,
V. _28_, 181-189.
Crow, S. C., 1968: Viscoelastic properties of fine-grained incompressible
turbulence. £. Fl. Mech., V. 33, 1-20.
Deardorff, J. W., 1972: Numerical invesitgation of neutral and unstable
planetary boundary layers. J_. Atm. SCJL. , V. 29, 91-115.
Hinkelmann, J., 1973: Personal communication.
Izumi, Y., 1971: Kansas 1968 Field Program Data. AFCRL Environmental
Research Papers, No. 379. Air Force Cambridge Research Laboratories,
L. G. Hanscom Field, Bedford, Mass.
Lilly, D. K., 1967: The representation of small-scale turbulence in
numerical simulation experiments. Proc. IBM Sci. Comput. Symp.
Environmental Sci., IBM Form No. 320-1951, 195-210.
Lumley, J. L., 1973: Modeling homogeneous deformation of turbulence.
Symposium on Turbulence in the Planetary Boundary Layer, Moscow,
August 28-31, 1973.
Lumley, J. L., and B. Khajeh-Nouri, 1973: Computational modeling of
turbulent transport, pp. 84-133 in Select Research Group in Air
Pollution Meteorology, First Annual Report. Department of Meteorology
and Center for Air Environment Studies, The Pennsylvania State
University, 210 pp.
Mahrt, L. J., 1973: A relationship of the Reynolds stress to local shear
and heat fluxes in the planetary boundary layer. J_. Atm. Sci. , V. J30_,
1577-1583.
Miiller, E., 1973: Ph.D. Dissertation, Universitat Mainz.
Rotta, J., 1951a: Statistische Theorie nichthomogener Turbulenz,
1 Mitteilung. Zeitschrift fur Physik. V. 129, 547-572.
Rotta, J., 1951b: Statische Theorie nichthomogener Turbulenz,
2 Mitteilung. Zeitschrift fur Physik, V. 131, 51-77.
Tennekes, H., and J. L. Lumley, 1973: A first course in turbulence,
MIT Press, 300 pp.
Wyngaard, J. C., S. P. Arya, and 0. R. Cote, 1974: Some aspects of
convective planetary layers. J_> Atm. Sci., V. 31, 747-754.
Wyngaard, J. C., 0. R. Cote, and K. S. Rao, 1973: Modeling the Atmospheric
Boundary Layer. 2nd IUTAM-IUGG Symposium on Turbulent Diffusion
in Environmental Pollution, April 8-14, 1973, Charlottesville, Va.
------- |