November 1971 851
UDC 551.508.313
ON STOCHASTIC DYNAMIC PREDICTION
1. The Energetics of Uncertainty and the Question of Closure
REX J. FLEMING
Air Weather Service, U.S. Air Force
ABSTRACT
I n this, the first part of a two-part study, we examine the rationale of the stochastic dynamic approach to nu-
merical weather prediction. Advantages of the stochastic dynamic method are discussed along with problems as-
sociated with the method. This method deals with the initial uncertainty by considering an infinite ensemble of
initial states in phase space, relative frequencies within the ensemble being proportional to probability densities.
The evolution of this ensemble in time, given by the stochastic dynamic equation set, is based upon the original
deterministic hydrodynamic equation set. One may consider the latter set as a subset of the former. Insight into the
nature of these equations is obtained by deriving the energy transformations associated with them. A simple baro-
clinic model is used to isolate the energy concepts and relations. The energetics yield qualitative and quantitative
information on the nature of the growth of uncertainty. It is found that the baroclinic instability mechanism is
responsible for most of the error growth as would be expected.
Previous predictability studies have considered that the simulation of the forces governing the atmosphere has
been perfect. The effects of imperfect forcing can be viewed with the stochastic dynamic equations by adding another
dimension to phase space for each parameter considered to be uncertain. The effect of the inclusion of this imperfect
forcing is shown by the new energetic relations that result, and by numerical calculation of the changes in the
growth of uncertainty.
The stochastic dynamic equations are faced with the same mathematical problem of “closure” found in ana-
lytical treatments of homogeneous isotroFic turbulence; that is, an approximation concerning higher order moments
must be made to close the system. A number of closure schemes are studied and it is found that the third moments,
which are individually small, should nevertheless be retained. It is shown in the equations and verified by numerical
calculations that the third moments do not affect energy conservation but affect energy conversion between uncertain
components, with the eventual result of altering the forecast of the mean. An eddy-damped third moment scheme
is found to give extremely accurate results when compared to Monte Carlo calculations.
1. INTRODUCTION
There is considerable interest today in the limits of
predictability of the atmosphere. Perhaps a more relevant
concern than the “limit” is what the practical predicta-
bility values will be for the various scales of atmospheric
features in the near future, say, 1, 5, or 10 JTT from now.
What effort must be expended t o achieve what degree
of predictability? This study is concerned with a method
that can deal with the above concern. Moreover, a more
prgfound question that some meteorologists must even-
tually ask themselves is, “What course does nu&e$ical
prediction follow after these practical limits of pr$dicta-
bility for the various features have been reached?”
One answer to the above question is to use the uncer-
tainty of the forecast itself as information. Certainly,
the variance or uncertainty of any meteorological variable
is important information, and may be as valuable as the
forecast of the variable itself. A technique for coping
with the uncertainties that exist in the initial conditions
of a numerical forecast is explicitly given in the method
of stochastic dynamic prediction (Epstein 19693).
In numerical weather prediction, there are several
sources of uncertainty that are part of the mathematical
initial-value problem. Acknowledging these sources, the
stochastic dynamic formulation of that problem forecasts
the expected value and uncertainty of every meteoro-
logical variable at each point in the space domain and at
every instant of time. No assumptions are made concern-
ing the physics of the original problem.
Some of the uncertainties that exist in the numerical
forecasts are due to errors in the initial conditions. A list
(by no means exhaustive) of sources of uncertainty in the
initial conditions of a numerical forecast might be: errors
in the raw data (instrument, transmission, etc.), errors
arising from the heterogeneous data set whose signal only
partially describes the complete inertial-gravity wave and
turbulence phenomena that exist at a particular instant
in the atmosphere, errors caused by the complete absence
of scales of motion unresolved by the paucity of the data,
“random” analysis error (Thompson 1957) arising from
the reconstruction of the grid analysis from a finite num-
ber of data points, errors caused by loss of significant
kinetic energy near the jet stream from smoothing by the
analysis technique, and that due to further smoothing of
the smaller scales that is usually necessary to avoid
amplification (beyond physical reality) of the short-wave-
length components. -
The above are some of the errors that will limit the prac-
tical value of predictability. In addition, the numerical
Now =signed to the National Meteorological Center, National Weather Service,
NOAA, Suitland, Md.
852 MONTHLY WEATHER REVIEW V O l . 99, No. 11
models used in forecasting are not perfectly simulating the
physical processes of the atmosphere. Indeed, some of
these processes have not yet been fully explained. This
obviously decreases the practical value of predictability.
If one could assume that all the above errors could be
eliminated and, in addition, if it could be assumed that
the physics of the atmosphere were “perfectly” known
and/or could be perfectly parameterized for numerical
computation, one would still be faced with a source of
error. The hydrodynamic fields must be represented in a
k i t e manner-in terms of a space mesh or in terms of
orthogonal coefFxients. Thus, the smallest scales of
motion that are unresolved must be ignored or treated
statistically. However, the nonlinearity of the equations
causes statistically treated small scales to give a random
perturbation to the motions of the large scales, causing
these scales to be increasingly unpredictable with time.
This “aliasing effect,” which is summarized by Leith
(1969), increases computational instability. The computa-
tional instability can be dealt with in a number of ways
(using a filtering technique at various points in time, using
a finite-diff erence scheme that inherently smooths, using
an artificial viscosity) ; however, the generation of uncer-
tainty at higher scales still remains. It is this last type of
uncertainty that is generally felt to be responsible for the
limit of predictability of various scales.
The nonlinear interactions between scales (known and
unknown), coupled with physical forces that are not per-
fectly understood, pIoduce uncertainty and make state-
ments of practical predictability only “speculation” when
the usual meteorological equations are used.
The stochastic approach to predicting the atmosphere
appears to be capable of reducing that speculation of
predictability. Consider a model of the atmosphere
described by a certain number of variables. A stochastic
model would specify the complete joint probability dis-
tribution of all the variables at each point in time, and
the whole process conceived as a continuous development
in time would be a stochastic method.
The model equations used in predicting the atmosphere
are assumed to be deterministic; that is, the exact present
state of the dependent variables completely determines
the exact future state of those variables. Since the initial
conditions of the atmosphere are known only in a proba-
bilistic sense, Epstein’s method seems appropriate. As
will he seen later in the paper, the deterministic equations
are but a subset of the stochastic equations and, thus, a
tremendous amount of computer power is required to
utilize the stochastic set. There remain some fundamental
problems and questions to be solved and answered before
these equations can be incorporated on a practical basis.
The purpose of this paper is to resolve those problems
and questions left open in the paper of Epstein (1969b).
The rationale of the stochastic dynamic equations will
be critically evaluated. The physical meaning of the
relationships that exist within these equations will be
discussed for the first time. An energy diagram will be
derived which shows the energetic relationships that
must exist in any prediction or simulation model in which
an account of the uncertainties is made.
The stochastic dynamic equations are actually an in-
finite unclosed set of coupled equations which cannot be
solved until a method of closing them is devised. This
question of “closure” is dealt with and is viewed in
analytical terms.
Freiberger and Qrenander (1965) discuss the same
problem that this paper is concerned with-how to deal
with the uncertainties that limit predictability, rather
than what the limit of predictability is. They propose
stochastic solutions to “limited” problems using an
approach analogous to classical statistical mechanics.
They did not offer the meum to the solution of the com-
plex atmospheric equations, but stated what form that
solution should take. A relevant quote from their paper
is: “Starting from the dynamic equations, whether in
their primitive or in a modified form, we should consider
an ensemble of possible solutions. The uncertainties and
lack of completeness of our initial weather data would
then correspond to a probabilistic superstructure on the
ensemble. . . . The conceptual, analytical and numerical
difficulties inherent in this approach are so overwhelming
that prospects of a completely successful analysis seem
rem0 te.”
It is hoped that the present study is a significant step
toward unraveling those conceptual, analytical. and
numerical difficulties.
I. THE STOCHASTIC DYNAMIC EQUATlONS
To visualize the difference between a deterministic
and stochastic dynamic forecast of a system, consider
that there are N dependent variables describing the
system. The deterministic prognostic equations can be
expressed in the general form
where the X t are the N dependent variables describing
the system. The N variables form an N-dimensional
Euclidean space whose coordinates are X,, . . . , X N .
Such a space is called a phase space, where each point in
the phase space represents a possible instantaneous
state of the system.
The initial state of the atmosphere is then represented
by a single point in the phase space and the determi-
nistic forecast gives the trajectory of that point in phase
space. In figure 1, the N-dimensional phase space is repre-
sented by any two dimensions. In this figure, point S is
the initial position of a single point in phase space and
S‘ is a point on the deterministic trajectory of that same
point at a later time t ,.
The results of all the possible uncertainties in the
final numerical analysis leads t o a great number of possible
initial states in the previously defined N-dimensional
phase space. Each of the different initial states will yield
a different trajectory in phase space even though sub-
jected to the same set of prognostic equations. Slightly
November 1971 Rex j. Fleming a53
FIGURE 1.-Evolution of a point and an ensemble in phase space.
different initial states may yield significantly different
trajectories. This is to be expected in light of Lorenz’
(1963~) discussion of nonperiodic trajectories. He con-
cludes that two atmospheric states differing by imper-
ceptible amounts may eventually evolve into two
considerably different states.
Is there a way to know which initial state is the correct
one? Were the matter of initializing the dependent
variables of the atmosphere merely a question of in-
accurate random measurements, then, in theory, an aver-
age of an inJnite set of measurements would converge to
the “true” state. However, operational predictions
depend upon a single set of data and it is impossible to
know the actual error. Moreover, the determination of the
proper initial conditions of the numerical prediction
problem is more complex than the theory of random
measurements.
Consider an observing net with the same data density as
exists today or even as it might exist 10 yr from now.
Even if each individual observation were a perfect rep-
resentation (Heisenberg’s Principle of Uncertainty ignored
for the moment), such data would not be acceptable
as initial conditions of the numerical prediction model.
Each measured signal would be a superposition of many
atmospheric modes, some for which a complete spatial
description could not be given by the coarse data res-
olution. Indeed, this is why the microscale features,
inertial-gravity wave oscillations, turbulence, etc., are
smoothed and/or modified in the final numerical analysis.
The analysis is made physically and computationally con-
sistent with the numerical model. From the viewpoint of
the model, then, there is no true initial state of the atmos-
phere. However, there ought to exist some state which is
a t least statistically optimal and consistent with the
objectives of the numerical prediction scheme.
The fact that the observations do have errors further
complicates the selection of that optimal initial state
which is consistent with the intentions of the prediction
model. It seems logical, then, to express the initial con-
ditions in terms of a probability distribution. This
would mean considering an infinite ensemble of initial
states in phase space with relative frequencies within the
ensemble proportional to the probability densities. This
initial joint probability distribution will evolve in a non-
linear manner according to the dynamic equations. It
will be seen that the stochastic method only forecasts the
low-order moments of that distribution, and that the
lowest order moments, the means, are by definition the
“expected values” of the dependent variables consistent
with the original analysis, all the uncertainties, and the
nonlinear dynamic processes a t work during the length
of the forecast.
Such an ensemble is shown in figure 1 where E and S
coincide and the ensemble mean, E, evolves to E’. Because
of the nonlinearities, there will eventually be a time
tl where E’ and S’ would differ. This is seen in independent
Monte Carlo and stochastic dynamic calculations pre-
sented in the latter portion of this paper and is also dis-
cussed by Epstein (1 9 6 9 ~~ 1969b).
How one might obtain an approximation of the ini-
tial joint probability distribution or, at least, of the low-
order moments of that distribution will not be the subject
of this paper. For our purposes, we shall consider that
by some method (deterministic or stochastic) the best
possible physically consistent analysis has been made.
Because of all the associated uncertainties, we consider the
analyzed values of the N variables to be expected values
having a variance or measure of uncertainty associated
with each.
There remains the problem of determining the evo-
lution of the ensemble. In terms of the atmospheric pre-
diction problem, this has recently been accomplished
independently and at about the same time by Tatarskii
(1969) and Epstein (1969b). The final result of each author
was the same, but each had a different conceptual basis for
the techniqu,e and each had a different interpretation
of the result. Both authors advocate using the dynamic
equations to find the probability distribution of future
meteorological fields, not the complete probability den-
sity distribution over the entire phase space-this is an
impossible task-but only the first and second moments
of the distribution over phase space. Both drop third and
higher order moments to close the system of equations.
Tatarskii uses a simple nonlinear form of a prediction
equation and from it derives the Liouville equation de-
scribing the evolution of the probability density. Epstein
considers the one-dimensional probability equation of
Gleeson (1966) but extends it to N-dimensions, correspond-
ing to the N-dimensional phase space under consideration.
This equation is exactly the same as Tatarskii’s Liouville
equation. Epstein calls his equations the “stochastic
dynamic” set of equations for the atmosphere. This
nomenclature is adopted here. Epstein presents numerical
results using a simple barotropie. model with only eight
orthogonal functions representing the space domain, where
the amplitudes of these orthogonal functions are functions
of time. He drops third and higher order moments and
calls the resulting equations the “approximate” stochastic
dynamic equations.
Consider writing the hydrodynamic equations in the
general form given by Lorenz (1963~) as
where p and q are dummy indices, (’) refers to a time
derivative, and the a’s, b’s, and e’s are constants describing
8 54 MONTHLY WEATHER REVIEW Vol. 99, No. 11
nonlinear interactions, external forces, and dissipative
mechanisms. The dependent variables, X i , could be grid
point values or amplitudes of orthogonal functions, in
ing the results of the truncated system of equations, is the
fundamental question to be answered.
principle, but the latter form of specifying the spatial
dependence will be used here. From the definition of
expectation, the expected value of f(X) a t time t is
given by
EIKWI=Jrn --m f(-%(Z wz (3)
w
where dX=dXl,, dX2, . . ., dXN and +(z, t ) is the probability
density that satisfies
JJ. . . J+(Z, t)dXi,dXz, . . ., dXN=l
over the phase space. One can then derive the following
relationships for Xf, any variable of the set of variables
XN that makes up the phase space:
E(X)=Pr, (4)
E(Xf Xi) =Pi Pj + uti, (5)
+ ~jut.t+ ~l.turj+ r i j k (6)
and
where :
-% Xk) =PtPj Pk+ Pf Q j k
pr is the mean of Xr,
utj is the variance if i=j, covariance if i #j ,
7 f j k is the instantaneous third moment about the
and where Xi, pi, uu, T f j k are all functions of time. The
equations for the means and the second moments about
the means corresponding to the deterministic set [eq (2)]
have been derived by Epstein (196%) :
mean,
P, !7 P
If there were no uncertainties, then the upq in eq (7) would
be zero, there would be no need for eq (8), and eq (7)
reduces to eq (2).
Solving eq (7) and (8) requires knowledge concerning
the third moments about the mean. One can derive a.n
equation for the third moments (cf. Fleming 1970),
but it involves fourth moments. The complete stochastic
dynamic equations form an infinite unclosed set of coupled
equations. The equations are thus unsolvable until an
algorithm is devised to close the set. The algorithm or
closure scheme chosen by Epstein was just t o drop the
third moment terms in eq (8). The implications of this
technique and other closure methods are discussed below
where the closure problem is studied more fully. It will be
mentioned here, though, that even if the third moments of
the ensemble are initially zero (for example, the variables
are multivariate normally distributed initially) they will
not remain zero. The time required for the third moments
to grow and the extent t o which they will grow, thus alter-
3. APPLlCATlON TO A
Epstein (1969b) used a simple barotropic model of the
atmosphere for his stochastic dynamic equations. In a
barotropic atmosphere, there is no change of wind with
height so that the atmosphere can be characterized by a
single level. Since the early days of numerical fore-
casting, the barotropic model has been used to forecast the
500-mb level of the atmosphere. I n such a model, the
only energy exchange that can take place is between the
kinetic energy of the eddies and the kinetic energy of the
zonal flow. Thus, it is the next logical step t o use a more
complex model of the atmosphere for the stochastic
dynamic equations. On the other hand, there are many
fundamental questions to be answered concerning the
dynamics and energetics of an ensemble in phase space, so
that it is desirable to keep the model simple enough t o
isolate the answers to these questions.
A model which ideally suits the purpose here was used by
Lorenz (1963b) to study different regimes of flow. This
model was extended by increasing the spectral resolution
(Lorenz 1965) in a study of predictability. The ex-
tended model transforms the space domain into orthog-
onal functions which describe three different wave num-
bers in the x-direction and two orthogonal modes in the
y-direction. ([‘Mode” refers to the y-variation.) This
model is the basic model used in this study. It is also
extended to include 15 different waves in the x-direction
for a predictability study to be described in part 11.
The details and equations of this two-level model were
developed by Lorenz (1960, 1963b, 1965) and are only
stated here. The basic equations are the quasi-geostrophic
vorticity equation and the thermal equation, that can be
mitten with 5, y as horizontal coordinates and pressure,
p , as the vertical coordinate. The equations for the model
written in spectral form are
where + and e are nondimensional dependent variables
representing the mean wind and mean potential tempera-
ture; K , K ’ , h, u, and the e; (all described below) are non-
dimensional constants; ai and C,,, are constants depending
upon the choice of the orthogonal functions F, used to
represent the space domain (described below). Quantities
have been made dimensionless by appropriate scale factors
including a length L and a timef’ (fis the Coriolis param-
eter assumed constant in this case). The dot in eq (9) and
(10) denotes differentiation with respect to the dimension-
less time.
November 1971 Rex J. Fleming 855
A simple two-layer model like this is incapable of de-
scribing the complicated boundary-layer physics of the
real atmosphere. Therefore, it is necessary to parameterize
the boundary-layer effects by the use of coefficients of
heating and friction. A frictional drag with coefficient 2~
is introduced at the lower surface, proportional to the
velocity in the lower layer. A frictional drag at the bound-
ary separating the two layers with coefficient K‘ is intro-
duced, proportional to the shear at this surface. There is a
heating of the lower layer, proportional to the difference
between the temperature of the lower layer and a pre-
assigned temperature field e*, characterized by a coefficient
2h. A coefficient for heat exchange between the layers
cancels out when u (a measure of static stability) is con-
sidered constant as discussed below.
In his application of this model to a consideration of
predictability, Lorenz (1965) considered u as a constant in
space and time. This assumption concerning u is used here,
so that predictability results can be compared t o those of
Lorenz. If u is a function of pressure, as in the real atmos-
phere, this leads to products of three variables-which, in
turn, means that at least some third moments must be
retained. If the model was “primitive” instead of quasi-
geostrophic, the equations would contain no triple products.
To solve these equations, we must determine the values
of a, and c i j k . These values are determined by the chosen
set of orthogonal functions F,. The choice of the functions
depends upon the geometry of the space domain in which
the model equations are applied. Following Lorenz (1965),
these functions are given over an infinite channel of width
nL having walls at the surfaces y=O and y=rL, where
the flow in the channel is periodic along the length, with
a specified fundamental wavelength. The functions are
=1,
@om=& cos mylL m=1, 2,
‘P,,=2 cos (mx/L) sin my/L
@;,=2 sin (naz/L) sin my/L
m=1,2; n =l , 2, 3, (11)
‘)71=1,2; n=l, 2, 3
and
where a is an integer chosen to specify the ratio of the
fundamental wavelength in the 2-direction.
It is readily verified that these functions are indeed
orthogonal. For clarification in deriving the energy re-
lationships for the stochastic dynamic form of the above
deterministic spectral equations, the following notation
is adopted for the 14 functions:
and
FIGURE 2.-Initial stream field $ for deterministic and stochastic
dynamic methods.
Note that functions 1 through 7 represent mode 1, and
8 through 14 refer to mode 2. The zonal components
of the flow are functions 1 and 8; the eddy components
are all the other functions.
The method of obtaining the a; and C i J k is given by
IAorenz (19636) and their values are given by Fleming
(1970). There are certain “selection” rules that admit
desired sets of pairs of functions for the nonlinear term on
the right side of the spectral equations; these pairs of
functions are those which interact with the third function
on the left of the spectral equation. For the 28-variable
model (14 orthogonal functions for $ and e), the different
pairs which interact with each of the respective functions
are determined from the “selection” rules and stored in
tables.
The physics of the model to be used can be conveniently
summarized. The quasi-geostrophic equations have been
used to describe thermally forced rotating flow in the
domain of the infinite channel. Sound and gravity waves
have been filtered from the solutions by the hydrostatic
and geostrophic equations. The transport of momentum
by the vertical motions and by the divergent part of the
horizontal motions has been left out of the vorticity
equation. However, the transport of heat by the total
horizontal and vertical motions is present because no such
assumptions were made in the thermal equations. Flow
in the model is slowed by internal and surface friction.
There is no vertical velocity a t the top and bottom of the
model, and no exchange of heat and momentum through
the sides of the channel.
In all the calculations in this study, the finite difference
scheme used for the time derivative in eq (9) and (10) is
the fourth-order Kutta-Simpson method which is described
in many texts (e.g., Hildebrand 1965).
An example of the results of a deterministic and
stochastic dynamic forecast with the model is given below.
The initial stream field for both methods is shown in
figure 2. The linear and constant terms have been set to
zero which means the model is now adiabatic and friction-
less. The mean values of $t and Or: are the same as the
3.i and et for the deterministic forecast. The ensemble
members are normally distributed about the mean for the
stochastic model. The variables have a standard deviation
$56 MONTHLY WEATHER REVIEW Vol. 99, No. 4 4
FIGURE 3.-Stream field after 7 days for deterministic (top) and
stochastic dynamic (bottom) cases. The contour intervals are
not the same.
of 10 percent of the mean, but the initial variables are
uncorrelated-all covariances are zero. The results after
one week are shown in figure 3. We see that the deter-
ministic forecast differs from the ensemble mean forecast.
The features differ in intensity, position, and shape. As
pointed out by Epstein (19693)) this does not mean the
deterministic forecast is wrong, as it is possible that the
deterministic solution would verify since it is a member of
the ensemble. The stochastic solution represents the en-
semble mean and stochastic solutions will be more
representative since they have been designed to minimize
mean square error based upon all the sources of
uncertainty.
It will be seen in later figures presented in this paper
that, with small uncertainties confined to the initial
conditions, the deterministic forecasts of the variables
and the stochastic forecasts of their means will be essen-
tially alike in the early stages of the time integration.
However, the advantage of the stochastic equations is
the additional knowledge given about the uncertanities
of these variables and about how they are dynamically
coupled.
4. ENERGETICS OF UNCERTAINTY
The results shown in the preceding section are similar
to those of Epstein (19693). They give rise to certain
questions. How does the uncerhinty grow? Why are some
physical situations more prone to error growth than others?
0 .5
2 .0
FIGURE 4.-Quasi-geostrophic energy diagram. Values from Wiin-
Nielsen (1968).
What fundamentally determines the evolution of the
hypothetical ensemble of points in phase space? An
attempt will be made to gain insight into the answers to
these questions by viewing the energy processes. It is
realized that complete answers would depend upon ana-
lytical solutions of the originel nonlinear deterministic
set of equations which, of course, are not ovailable.
Many authors have contributed to the understanding of
the atmosphere by considering the energetics of the
deterministic equations. The intent here is to give a
brief summary of the basic energy concepts. This will
provide a foundation upon which to build the concepts of
the energetics of uncertainty.
There are three significant forms of energy in the
atmosphere: internal, potential, and kinetic. I n a hydro-
static atmosphere, the internal and potential energy are
proportional and are conveniently combined to form the
“total potential energy.” In dealing with the exchanges
between these types of energy, Lorenz (1955) found it con-
venient to define the ‘‘available potential energy” as that
portion of the total potential energy which may be avail-
able for conversion into kinetic energy. Lorenz also
derives an approximate expression for available potential
energy which has since been shown t o be energetically
consistent with the quasi-geostrophic equations (cf.
Wiin-Nielsen 1968). Representative values of the gen-
eration, conversion, and dissipation of energy obtained
from the quasi-geostrophic calculations are given by
Wiin-Nielsen (1965, 1968). His values, shown in figure 4,
do not reflect the effects of the Tropics. In this figure,
A and K refer to available potential and kinetic energy;
the subscripts Z and Erefer to zonal and eddy components
of the energy, G(X) refers to a generation of X, D(X) re-
fers to a dissipation of x, @(X,Y) refers to a conversion of
X into Y. The conversion C(AZ,Kz) is taken from
Dutton and Johnson (1967).
A positive generation will occur when there is diabatic
heating where it is warm and cooling where it is cold.
The dissipation of kinetic energy is due to friction effects
in the boundary layer and internal friction within the
atmosphere. The conversion C (Az,&) depends upon trans-
porting heat from a warm latitude to a cold latitude, and
November 1971 Rex J. Fleming 857
alternately, the transport of cooler air from a cold lat-
itude to a warm one. The conversion C(KE,Kz) depends
upon whether the transport of angular momentum is
toward latitudes of higher or lower angular velocity. The
conversions C(A,K) are positive when warm air is rising
and cold air is sinking.
A more complete derivation of the stochastic dynamic
energy relationships obt,ained from the stochastic dynamic
model is given below, but all those deterministic energy
relationships will be used without derbation. The baro-
clinic stochastic dynamic equations and their energetics
will yield, for the first time, substantial insight into why
and by how much a deterministic forecast may differ
from the real atmosphere
5. DERIVATION OF ENERGY EQUATIONS
In the absence of heating and friction, the condition
~a r t q X r X d j q =O on eq (2) implies that F X : is con-
served (cf. Lorenz 19634. The corresponding adiabatic
frictionless model [given by eq (9) and (10) with the linear
and constant terms set equal to zero] conserves the quan-
tity (K+A) where K is kinetic energy and A is Lorenz’
approximated available potential energy. It is easily
verified that
and
S
A=- e:: 2 ~0 i
where s is a scale factor given by P,,LLy2!g (9 is the ac-
celeration of gravity). In the remainder of this paper,
the scale factor will be ignored and u0 will be referred t o
as u.
The stochastic equations applied to the general expres-
sion, eq (2), reveal that the quantity (&+ uii ) is * con-
served. Denoting the conserved quadratic quantities as
energy (they could represent enstrophy-mean squared
vorticity) reveals that the stochastic dynamic equations
provide for a partitioning of the energy as pointed o u t by
Epsteia (1969b). There is “certain” energy associated
with the expected components of the ensemble and “un-
certain” energy associated with the variances. The uncer-
tain energy cannot be specified” (e.g., as t o location)
within the model.
The initial uncertain energy may represent, in part,
some of the actual energy of the atmosphere; that is, it
may be energy lost by smoothing the data, or turbulent
energy, or any small-scale energy unresolved by the
deterministic initial conditions. (In the statistical approach
to analytical theories of turbulence, these variances are
referred to as the “turbulent” energy.) However, the
variances (whatever name they be given) are a direct
measure of the growth of uncertainty in the values of the
variables. The names used here are especially appropriate
in that the growth of uncertainty will be seen to depend
i
((
upon the dynamic coupling between the variables, and
the energy exchange processes associated with these
couplings.
The stochastic dynamic equations applied t o the deter-
ministic equations for the two-level model, eq (9) and
(10) , give the following result:
1
f 2a A =C - [e:+ var (e,)]
where var is variance. The quantity
is the certain kinetic energy;
is the uncertain kinetic energy;
is the certain available potential energy; and
is the uncertain available potential energy.
The quantities +i and e, in the above equations now
refer to the mean values of +, and 8, in the stochastic
dynamic model, w-hereas the same symbols have been
used for the deterministic stream function and potential
temperature There will be no confusion, as only the sto-
chastic dynamic equations will be referred to in the
remainder of this paper.
Taking the time derivative of eq (15) and (16) gives
The zonal component of either energy includes the terms
for which i= 1 and 8. The eddy component is computed by
summing the terms over i =2 ,3 ,. . . , 14; iZ8. Thus,
expressions for A,, Kz, A,, K, and their uncertain counter-
parts U.,, UKz, UA,, UKE can be obtained from the
expressions for i,, it, [vai (+,)I, and [vai (e,)]. By compar-
ing the terms on the right-hand side of the equations for
858 MONTHLY WEATHER REVIEW voi. 99, No. 11
the time derivatives of the various energy quantities, the
exact energy relationships that exist between the quantities
can be determined. Although this is a tedious manipulation,
a procedure is given here since the energy relations which
result are quite important. Instead of writing time
derivatives of the eight different energy quantities, much
space can be saved by “SelectiveIy” grouping terms
and showing that the “grouped” terms are the indicated
energy conversions.
The restriction of adiabatic frictionless flow is removed
and the complete equations for dI, it, [vG ($J ], and
[vai (e,)] are
$,=AMf+AMQi-CMt (19)
= uAS, + B, + u CS t - u Di + Fi - Ei + UASQi + BQt,
(20)
[va; ($JI =~AMU, -~cMu,, (21)
and
[vai (e,)] =2UASU, +2 BU, + 2aCSUi-2uDUf -2Eu,
refer to terms coming from first moments and eventually
define certain energy relations, and M or S refers to a
term’s source as being from the mean flow equation or the
shear flow equation. The terms that come from second
moments about the mean are
1 C S U f =r a;”# cov (e,, #J,
+ 1
where Urefers to terms that eventually define uncertain
energy relations.
The summation sign in the above defhitions implies a
double summation over k and Z with Z>k. Substituting
eq (19) through (22) into the equations for the time rate
of change of kinetic and available potential energy given
by eq (17) and (18) results in energy change equations
that are completely expressed with terms that are known
a t any given instant of time. The resulting equations are
~=C ~~~f ~~f +~Q f -C 2 0 1 ,~+~~~,(~A S i
-uDt+F,--{+uASQt+BQ,) +a?(mUf-C‘mi)
+ U: ( uASU, + B Vi + u CSU~-UDU~- E U t ) ] (26) and
A= E[& (AS + B flu+ GS 1 - D,+ FJu- E& i
+m&i+BQ,Iu) + (ASU<+But/u+ CSu<-oUi- EhJtIu) -
(27)
Selectively regrouping terms gives
November 1971
where
Rex 1. Fleming
Since Cikl=Ckl,=C1(k, eq (37) can be written
859
and the numbered brackets serve only to identify terms.
Equations (28) and (29) can be made dimensional by
multiplying each term by the scale factor (PoL2f3/g).
It can be shown that w ; corresponds exactly to the
expected value of vertical velocity in pressure coordinates
a t the surface between the layers. Thus, the bracketed
terms “l”, that appear in eq (28) and (29) with oppasite
sign, agree with the definition of C(A,K) from quasi-
geostrophic theory. The conversions for the certain A to K
are then given by
C(Az,Kz)=-~Ce,wi i=l and 8 (32) t
and
C(A,,K,)=--Ce,w; i=2,3, . . .,14; i #8 . (33) i
or
Thus, the first term on the right of eq (39) is equal and
opposite in sign to that of eq (36). This also implies that
the sum of all the components in term “9” is identically
zero.
We have shown that term “9” gives C(A,,A,), or
that the negative of term “9” gives C(A,,A,). By the
same method used above, it can be easily verified that
C(KE,Kz) is given by term “5” of eq (28).
In much the same manner as above, the conversions
C(UA;,uAd, C(U&,Az), C(AE,UAZ), and C(AE,UAE)
can be obtained from term “10” in eq (29). The proof of
these conversions is given by Fleming (1970). The method
involves considering the “zonal-eddy” triplets and all
distinct “eddy-only” triplets. Many terms are involved
that must be “tagged” as to their source, but the pro-
cedure is the same as above. The conversions for the uncertain A to K are similarly The remaini+g energy conversions C(UK,,~K,),
given by terms “2” : C(UKE,Kz), C(Ks,UKz), and C(KE,UKE) are obtained
from term “6” of eq (28) just as described in the preceding
i (34) paragraph. A complete summary of the certain and
uncertain energy generation, conversion, and dissipation
C(UA,,UK,)=-~wi’ i =1 and 8
c(uAE,uKE)=-~w;’ i=2,37 . * -,14; i#8* (35) equations is given below. i
(40)
(41)
The generation of available potential energy is pro-
In eq (29), since all the heating effect is contained in
E, and F,, the certain generation is given by terms “7”
and the uncertain generation by term “8”. These genera-
G(Az)=C-8,(8T-Oi) h i=l,8,
portional to the correlation ‘of heating and temperature. i @
G(AE)=x-ei(@-ei) h i=2,3, * * *, 14$#8, i Q
D (&) =E @K[@ - 2##3, + e:] + U: (2 K ’ )e : i= l ,8, (42) tions are more clearly identified at the end of thissection
where the notation is replaced and the complete energy i
equations are summarized.
Since all the friction effects are in CM,, CSi, and Di,
given in eq (28) by the terms “3” and “4”, respectively.
The conversion C(A,, AE) is given by the negative
of term “9” in eq (29). This is not as immediately evident
as the other conversions so a proof is given below. 3’ mce
each triplet of wave vectors capable of interacting down
to some truncated scale, only one of the triplets need be
considered for this proof. Consider the triplet (i,k,Z)
where i is a zonal function and k and I are eddy functions.
Equation (29), considering oniy term “Y’, gives
the certain and uncertain dissipation of kinetic energy is i
i=2,3, - .,14;i#8,
(43)
c(A,,A,) = xx 5 caZ (ek+ J i=l,8, (44) the spectral equations used here completely describe i k,1 Q
c (K E ,~,) =~,~,(~~--a ~)C i k z [~t (~k ~2 +e k e z ) i k,Z
+ ei 2 +#kO J 1 i=1,8, (45)
C (A,,&)= - C&?, i=1,8, (46)
i
(36) Az=Tef(ek+l-$ke,)+ c i k l . . -
C(AE,KE)= -x w :e r i=2,3, . - , 14;i#8, (47) and for the eddies. i
i=1,8, (48) AE=n[e,~k,f(e~$I-~~ei)+e2c,,(ei$k-#,s,)i+ 1 - - .. (37) G(UA,)=-x;var h (e,)
I
860 MONTHLY WEATHER REVIEW VOI. 99, No. f i
i=2,3, * * *, 14;i#8, (49) h - var (e,) G(UA,) = - i a
D(UKZ)=x&[var (#,)--2 cov (#,,e,) +var (e,)] i
+U?(’LK’) var (e,) i=1,8, (50)
and
The energy generations, conversions, and dissipations
of the baroclinic stochastic dynamic model are indicated
in an energy diagram shown in figure 5. The energy
equations would be the same for any number of orthog-
onal functions, as the indices i, k, and I would simply
have more values. If the original model were not a quasi-
nondivergent model but a “primitive” equation model,
the right-hand side of the energy equations would be
different; however, the energy diagram in figure 5 would
still be the same.
The advantage of the stochastic dynamic method and
of these energy equations is that one can quantgy the
growth of uncertainty. Various calculations are included
in this paper to bring this point out. Here, however, we
will only qualify certain features of the energy diagram
shown in figure 5.
Numerical computations indicate that the largest
energy conversion (involving an uncertain energy) is
Az to UAE. This is associated with baroclinic instability
as is the conversion Az to A,. Large-scale baroclinic in-
stability is the prime mechanism for transporting heat
from low latitudes to higher latitudes. The resulting
large-scale eddies are of importance in maintaining the
general circulation of the earth’s atmosphere against
frictional dissipation. It is not surprising that if the
physical variables involved in such a large-scale nonlinear
mechanism are not perfectly known initially, the future
values of those variables (after the instability mechanism
has had time to work) mill be known even less perfectly.
This situation reaches the worst possible stage when the
“true” values of the variables are near a critical condition
where instantaneous development may or may not occur.
November 1971 Rex J. Fleming 861
- N
Y
Y
Y
A
C(+, KE)
G IAE) D IKE)
FIGURE 5.-Stochastic dynamic energy diagram.
The conversions C(Az, UA,) and C(A,, UA,) given by
eq (56) and (57) show the extent of this nonlinear growth
of uncertainty. Both conversions are concerned with the
uncertain transport of an uncertain amount of heat
northward, but the former is by far the most important
as it involves the certain or mean latitudinal temperature
gradient (a relatively large quantity) multiplying un-
- certain eddy components, while the latter conversion
involves certain eddy temperatures (relatively small
quantities) multiplying terms involving the uncertainty
of latitudinal temperature gradient (anti other terms)
which tend to remain small, relatively speaking.
The conversions C(K,, UK.) and C(&, UK’) involve
barotropic instability. These conversions are much
smaller than C(A,,UA,), as would be expected since
barotropic instability plays the lesser role in the real
atmosphere. These conversions are discussed in this
section.
The conversions C(A,, UA,) and C(&, UK,) are
due to the nonlinear exchange of energy between wave
numbers other than zero, that is, the eddies. These con-
versions can and do change sign, but their amplitudes
are small and they will not be referred to in the remainder
of this paper.
The conversions C ( UAz, UKz) and C ( UAE, UK,) in-
volve the variance of the vertical velocity, which is reason-
able, since the conversions C (Az,Kz) and C (A,,KE) involve
the vertical velocity in the deterministic calculation, and
involve the expected value of the vertical velocity in this
stochastic calculation. There are, however, missing terms
in eq (54) and (55) for the uncertain conversions. The
reason is simply that the vertical velocity involves a
Jacobian, which in turn implies that eq (54) and (55)
involve a triple correlation of moments about the mean.
Up to this point, we have used the simple closure scheme
of neglecting third moments. It is, therefore, necessary to
keep third moments if an accurate account of C( UA,, UK,)
anti C(UA,,UK,) is to be made. These additional terms
are described later in this paper when the problem of
closure is discussed.
Up to this point, the uncertain generations and dissipa-
tions are necessarily constrained to act in such a way as
to remove uncertainty; that is, G(UAz) and G(UA,) will
be negative (removing energy from their respective boxes)
and D(UK,) and D(UK,) will be positive (removing
energy from their respective boxes). This occurs because
the physical parameters we have used to simulate the
forces working on the atmosphere have been assumed to
be known perfectly. This is what all deterministic forecast
and simulation models do, also. The implications of this
assumption are discussed in part 11, but the effects on
the energy generation are discussed in a later section of
this paper.
Figure 5 reveals no conversion from AZ to UAz and Kz
to UK,. This is because there is no communication be-
tween the means of the zonal quantities and their respec-
tive variances. This, in turn, is because the deterministic
prognostic equations for the zonal components do not con-
tain other zonal components in the nonlinear terms on
the right hand side of their respective expressions.
ENERGY FLOW IN A BAROTROPIC SYSTEM
The energy flow characterized by figure 5 will be shown
in numerical calculations. Fortunately, there are physical
situations in which the numerical calculations must give
exactly predictable results-this will suggest that the com-
puter progmm is working properly. Numerical calculations
of the complete energy conversions using eq (40) through
(61) were carried out to convince the author that those
conversion equations were correct. By adding from the
conversions the input and output t o each energy “box”
(through a complete Kutta-Simpson time step), we should
obtain net changes that exactly agree with the energy box
values before and after the time steps. Such calculations
were instrumental in discovering initial algebra errors,
but eventually these calculations were exact”, well
within the bounds of roundoff error.
The simplest case to consider is reversion back to a
barotropic model by assuming that all Oc, linear terms,
and constant terms are equal to zero. In this model, the
sum of zonal anti eddy kinetic energy is conserved. The
energy diagram reduces t o that of figure 6.
The total energy in the four boxes must be constant for
all time. Results for a particular case are shown in table 1.
It is seen that the required energy is conserved. The last
digit in the total is in doubt because of the printing format
used, but there is an extremely small “time” truncation
error as well as roundoff error present. Waves 2, 4, and 6
were used in this case. The time step was equivalent to 1.5
hr .
(L
862 MONTHLY WEATHER REVIEW Vol. 99, No. I 1
UKE 3.068
.44 52 .63
161 1121 1181
FXGURE ?’.-Initial energy conditions for stable-unstable barotropic
test. Energy is in kJ.m-2.
KE 998.19
0.013
- - - 188 175 162
161 I121 (181
FIGURE 6.-Barotropic stochastic dynamic energy diagram.
STABLE
TABLE l.-Barotropic energy conservation (kJ m-2)
Time (days) KA Ks UKz UKg Total
0 1040.0 1120.0 10.0 172.0 2342.MMOO
7 728.1 832.0 269.5 512.4 2342.oooOl
14 766.7 403.6 167.8 1003.9 2342.0004
21 821.8 470.8 220.3 829.1 2341.99997
An important point concerning the value of the sto-
chastic dynamic equations is brought out in table 1. There
is not much that can be said about the eddies after two
weeks, as UKE is about two and one-half times as great as
KE. Note that at 3-weeks, however, there is more certain
energy in the system. The stochastic dynamic model does
not allow the uncertaint5- t o grow without limit-and
this is realistic. In a barotropic flow, the energy distribu-
tion oscillates back and forth with more energy alternately
in the zonal component and eddy component. A character-
istic period of 5 to 6 days for this process is shown by
Thompson (1957), Wiin-Nielsen (1961), and others. At 2
weeks, a greater portion of the energy is in the eddies but
there is great uncertainty as to exactly where the eddies
are! The dynamics of the model then dicthte that the
eddies should feed back their energy to the zonal flow, and
at 3 weeks this has happened. The uncertain eddies give
up their energy properly because the stochastic dynamic
equations keep track of the uncertainties “dynamically.”
The flow is actually more predictable at 3 weeks than at 2
weeks. The zonal component is more predictable than any
eddy components; there is one less degree of freedom.
The growth of uncertainty in a barotropic model is
studied with two cases having the same initial energy.
The two cases differ only in the specification of the mode
“2” zonal component amplitude, which results in an
initially stable or unstable barotropic flow. Only the sign
of this component is changed so that the energy is the same
in each case; however, the sign change results in satisfying
instability requirements.
F=l
2.499
F l
1001.24
UNSTABLE
F l
998.54
IJKE 3.07 KE 1WO.82
0.013
-- TJ1 184 176 -
161 1121 1181 16) 112) iiai
FIGURE &-Energy conditions after 3 hr for stable-unstable case
with E in kJ.m-2, E in 10-4 kJ.m-2.s-I.
Figure 7 shows the initial conditions of each case where
appropriate values of +$ have been chosen-to give the
indicated energies, Waves 6, 12, and 18 in the 2-direction
are chosen. Figure 8 shows the energy values and con-
version rates after 3 hr. I n the stable case, the waves
are feeding energy to the zonal flow, while the reverse is
true for the unstable case; where the waves are growing
at the expense of the zonal flow. Note that in both cases
the growth of uncertainty is through C(KE,UKz) and
subsequently C( UKz, UK,). The uncertain energy U&
is acting as a “catalyst” in the early stages of this growth
of uncertainty.
A look at the behavior of the kinetic energy in mode “1”
of waves 6 , 12, and 18 is also given in figure 8. Initially,
each wave had the same amount of certain and uncertain
energy. Since wave 18 is more stable than the other waves,
November 1971 Rex J. Fleming 863
TABLE 2.-Energy check for adiabatic frictionless model with energy i n the units kJ.m*
Weeks A2 AB KZ KE UAz UAE UKz UKE Total
0 10, OOO 3926 1200 1200 100 39 12 12 16,488.300
1 1,628 2680 e34 1886 400 6497 78 3685 16,487.881
2 277 1402 480 1055 1225 7139 380 4627 16,487.791
3 424 3276 1948 444 1883 342 6172 181 16,487.569
4 114 2254 435 1435 832 5954 472 4990 16,481.411
5 142 1195 416 8803 1665 6738 468 5059 16,487.275
6 21 1804 421 963 1027 6852 276 6321 16,487.069
STABLE
1193.4
UNSTABLE
FIGURE 9.-Energy conditions after 3 days for stable-unstable case.
Units are the same as in figure 8.
it gives up its energy faster in the stable case. Wave 6 is
more unstable than the other waves and is gaining energy
faster in the unstable case. However, in both cases, the
growth of uncertainty is greater in the smaller scales or
higher wave numbers. If the energy diagram were further
broken down to show relations within wave numbers, the
cascade of uncertain energy to the higher wave numbers
would be seen. That the smaller-scale errors should grow
faster is nothing new or unnatural. Consider just a linear
displacement of a small-scale feature by the zonal wind or
by a large-scale wave. A small uncertainty in the large
scale has a large effect on the small scale, as discussed by
Lorenz (1 969).
Figure 9 shows the energy and energy transforms at
3 days. The stable case has less total uncertainty. Thus,
3-day forecasts from initially stable situations will have
better verification scores, or, in other words, the forecasts
can be made with greater confidence. The stochastic
dynamic method gives the expected state and a measure
of that confidence.
ENERGY FLOW IN AN ADIABATIC FRICTIONLESS SYSTEM
A more complex model, with more degrees of freedom
than the barotropic, results if the e t are retained. The
linear and constant terms are still zero. This is the adi-
abatic frictionless model. Figure 5 gives the corresponding
energy diagram (with the removal of the generations and
dissipations) and the sum of the energy of the boxes is
constant with time. Numerical results for a typical case
(table 2) show that this energy is conserved even when the
integration in time is carried far past the unpredictable
stage. This model is the least predictable, as is discussed
in part 11; its value here lies in the check it provides on
the energy equations and the model itself. The integration
has been carried out to 6 weeks with waves 2, 4, and 6
and a time step of 3 hr.
CERTAIN-UNCERTAIN ENERGY FLOW IN A BAROCLINIC
SYSTEM
The values of the heating and friction coefficients as well
as the choice of u will determine the type of motion that
will evolve. In a classic study showing the power of this
simple two-level spectral model with but a single wave in
the 2-direction, Lorenz (19633) was able to reproduce the
different regimes of flow found in the rotating-annulus
experiments of Hide (1953), including an irregular non-
periodic flow typical of the atmosphere. In that study,
Lorenz let h=~=2~' so that the controllable external
parameters were 0; and h. Starting with a weak zonal
flow with small superimposed perturbations, the different
regimes of flow occurred for different external parameters.
For his predictability experiment, Lorenz (1965) used
three waves in the 2-direction, and values of u, 0: and h
were chosen to give irregular nonperiodic flow. This study
follows Lorenz in setting h =~=2 ~' . The purpose of
this paper is to study stochastic concepts and not to sim-
ulate the atmosphere (which could hardly be accomplished
by a two-level model). Nevertheless, it is desirable to
choose parameters that give irregular nonperiodic flow
resembling the flow patterns found in the atmosphere.
This has been achieved by using values e:= 9/128,
h= 15/128, and u=5/128. The values of e; and h are slightly
less than those used by Lorenz (1965) for reasons given in
part 11. Only 07 was retained in the linear term, implying
a heating at the Equator and a cooling at the poles.
Changes of these parameters are used in this study to
MONTHLY WEATHER REVIEW VOI. 99, No. d 4
FIQURE 10.-Initial stream field $ of “return to Hadley circulation”
case.
FIGURE 11.-Stream field after 3 weeks of return to Hadley
circulation case.
clarify various concepts. When such changes are made, the
new values are given.
Energy is no longer conserved in the model when forc-
ing and friction are included. However, an additional check
is devised to insure that the stochastic dynamic energy
equations are correctly derived and correctly programmed.
A large value of u is used in a complex initial state of large
eddies. The large u means that the upper level has a much
higher potential temperature than the lower layer-a
very stable situation. With this stable atmosphere con-
taining frictional dissipation, the flow in the model should
return to the steady-state Hadley circulation which has an
analytical solution discussed by Lorenz (1963b).
For this run only, the value of Q is double its commonly
used value, that is, u=10/128. The initial stream field
for this case is shown in figure 10. The resulting flow after
3 weeks is shown in figure 11. The flow is approaching the
equilibrium Hadley circulation with the proper analytical
values for J.l and 8 ,. These values are obtained by solving
eq (9) and (10) for the steady-state zonal flow which
yields :
and
81 = $1 (62)
81 = 8?/(1+4 (63)
The initial values of J., and were 0.0452 and 0.0442,
FIGURE 12.-Energy after 1 day foy return to Hadley circula-
tion case with E in kJ*m-2, E in lO-‘ kJ.m-2-s-I.
FIGURE l3.-Energy after 3 weeks for return t o Hadley circulation
case. Units are the same as in figure 12.
respectively-they converged to the common value of
0.06521 which satisfies eq (63). The energy diagrams in
figures 12 and 13 show that the stochastic dynamic
energy equations we programmed correctly. The eddy
energy and exchanges have nearly vanished, and the
terms G(Az), @(Az, &), and D(Kz) are approaching a
steady-state common va,lue.
November 1971 Rex J. Fleming 865
DAYS
FIGURE 14.-Relative kinetic energy uncertainty as a function of
time in the return to Hadley circulation case. I INITIAL I
Another important point is brought out by this case.
The relative error of the eddies (uK,/KE) is decreasing.
It is natural to expect that the absolute uncertainty of
the eddies will drop because of the dissipation of the
eddies themselves, bubit is not obvious that the relative
uncertaintt would drop. In the barotropic case, it was
seen that when the situation was stable, the relative un-
certainty was greater than in the unstable case. Figure 14
shows the relative uncertainty as a function of time.
The stochastic equations are energetically consistent and
keep track of the uncertainty dynamically. As a result,
the volume of the ensemble of points in phase space will
eventually shrink to a single point-the same point that
the deterministic trajectory will ultimately arrive at-the
steady-state solution of the phase variables.
A particular example of the resulting energy trans-
formation for the barotropic case with irregular flow is
considered next. Many different calculations were made
by varying different initial conditions and forcing param-
eters. Quantitative discussions of these values and how
they change will not be a part of this paper. Moreover,
such a discussion may be of little value in that this
baroclinic model is so severely truncated and contains
the physically inhibiting use of a constant static stability
(cf. Lorenz 1960). However, the qualitative effect of the
growth of uncertainty can be assessed; in particular, the
relative values of different uncertain energy conversions
are typified by the example shown in figures 15 and 16.
Figure 15 (left) shows an initial state of a particular
case and the resulting energy values and conversions
after 1 day (right). Figure 16 shows these values after
1 and 2 weeks. Choosing and maintaining energy values
corresponding to atmospheric values is somewhat difFi-
cult (and perhaps meaningless) with such a severely
truncated system. However, the initial energy values are
close to winter values (cf. Wiin-Nielsen 1967) except for
AZ which is too large by a factor of two. The latter was
necessary to maintain nonperiodic flow in such a trun-
cated system and this will be further discussed in part 11.
Three general comments typical of all calculations and
indicated in figures 15 and 16 are: (1) the conversion
A, to UA, dominates in the growth of uncertainty; (2)
FIQURE 15.-Energy conditions in a case of irregular baroclinic
flow initially and after 1 day with E in kJ.m-2, E in 10-'
kJ.m-2.s-I.
1 WEEK 1 WEEKS
FIGURE 16.-Energy conditions in a case of irregular baroclinic
flow after 1 week and after 2 weeks. Units are the same as in
figure 15.
the resulting uncertain energy values are considerably
smaller for zonal components than they are for eddy
components; and (3) the effects of heating and friction
tend to slow the growth of uncertainty by acting to
reduce the variances of the variables; that is, there is
negative generation and positive dissipation of what has
been called the uncertain energy. This is considered in
greater detail in the following section.
6. EFFECTS O F IMPERFECT FORCING ON ENERGETICS
The atmosphere is primarily forced by net radiative
heating at the Equator and cooling a t the poles. However,
there are secondary forces present which will vary with
latitude and longitude, including complex radiation
processes; release of latent heat and evaporation in precipita-
tion processes; conduction of sensible heat; etc. None of
these atmospheric processes is known perfectly, although
some are known considerably better than others. A
866 MONTHLY WEATHER REVIEW voo. 99, No. 11
prediction model simulates these processes by various
parameters. The question that must be answered is, “How
does uncertainty in the physical parameters affect pre- 21.145
dictability?” The stochastic dynamic equations c a n
provide answers to this question, but to understand these
answers the energetics will be utilized.
Consider, for example, that the et are still constant but
that their values are uncertain. There is then a mean and
variance for each 0: that contributes to the forcing.
Thus, their deterministic tendencies are zero (being in
in the general form of eq (2) with the a’s, b’s, and c’s all
zero) and the stochastic dynamic equations involving
the 0: only become
(&=O (64)
and
I
Thus, the means and covariances involving only the 0:
do not change. However, even though cov (ei, 0:) may be
initially zero, it is readily shown (cf. Fleming 1970) that
correlations between the et and 0; will develop.
The effect of the terms cov (ei, 0:) on the energetics is
considered in the following where previous energy terms
are assumed present. Then eq (22) becomes
FIGURE 17.-Partial energy diagram after 3 days for a typical case
where primary forcing is known “perfectly.” Units are the same
as in figure l5.
cov (e,, e:)= . . . +~F U 2h [vai (e,)]= . . . +- .a:+ 1
(66)
where
Then eq (29) becomes
A= ...
Carrying this term through the (8) equation reveals
that -a:FU completes the definition of the uncertain
vertical velocity. The term [(a?u+l)/u] FU is seen to be
another component of the generation of uncertain avail-
able potential energy. When these are added to their
respective energy equations they become
and h Q(UA)= . . . +- cov (e,, e;). U FIGURE
66.9 I
18.--Partial energy diagram after 3 days for same case as
figure 17, but primary forcing parameter has 2.5 percent standard
deviation. Units are the same as in figure 15. Thus, the inclusion of the notion of uncertainty in
the diabatic heating contributes directly to the generation
of uncertain available potential energy and to the conver-
sion of UA to UK. Indirectly, it eventually contributes to all the terms.
The difference in expressing the physical forcing
“perfectly” versus expressing it “imperfectly” is fully
brought out in part 11. A manifestation of that difference
is revealed in figures 17 and 18 where partial energy
diagrams are shown for two different cases dealing with
primary forcing. In one case, the primary €arcing is
November 1971 Rex J. Fleming 867
assumed to be known perfectly; in the other case, the
0: has a standard deviation of 2 1/2 percent of the mean
value of e:. (The details of the initial conditions are not
given here, as any set of conditions will show the quali-
tative result that is brought out.)
Note that in the “perfect” forcing the G(UA,) is
negative, and when the forcing is “imperfect” the genera-
tion of the zonal component of uncertain available
potential energy is positive. From eq (4 8 ), with no
uncertainty in the heating term, G( UAz) is necessarily
negative. The heating in the model is such that the largest
latitudinal temperatire gradients are reduced and vice
versa. Diabatic heating tends to eliminate uncertainty
and push the system toward equilibrium. But if the heat-
ing is uncertain, temperatures would get higher if heating
is large and vice versa, et and 6: will become positively cor-
related, and G( UA,) may (and does) become positive
[cf. eq (70)].
Many interesting calculations that highlight the growth
of error can be made with the energy equations. The
nature of this paper will only permit a summary here.
We have stated that, by dealing with the uncertainty
dynamically, the expected forecast of an uncertain
ensemble of initial conditions can result and the growth
of uncertainty of the forecast can be viewed through the
energetics. It has been shown that the stochastic dynamic
equations that neglect third and higher order moments
conserve the ensemble energy in those cases where the
energy should be deterministically conserved. Epstein
(1969b) has pointed out that, while dropping the third
moments implies that the time derivatives of the indi-
vidual variances are inexact, the equation for the time
derivative of his kinetic energy did not contain third
moments. This is indeed true and is the reason that the
sum of the energy in the eight boxes of an adiabatic fric-
tionless energy diagram is conserved. It is not assured that
the energy exchange between these boxes is “accurate.”
This last point is emphasized in the following section.
7. THE CLOSURE PROBLEM
Epstein (1969b) has made the solution of the stochastic
dynamic equations tractable by simply disregarding third
and higher ordw moments. He notes that this closure
scheme is necessary in light of the vast amount of com-
putation required. He also notes that the success of the
deterministic equations, which are but a form of the
general stochastic dynamic equations with the variances
implicitly assumed to be zero, suggests that the truncated
set will be more successful than the deterministic, and
should be valid if the time domain is not too large. How-
ever, there is evidence that this closure is not a good one.
That third moments which are initially small or zero will
not significantly alter the forecast of the first moments
until they become large is not an obvious fact. It is pos-
sible that the sum effect of many small third moments
will be significant. There is also the question of their
effect on the second moments. To underscore this last
concern, consider the history of analytical theories of
turbulence.
turbulence
originated from the statistical approach of Taylor (1921).
These theoretical investigations are usually limited to
the idealistic case of homogeneous isotropic turbulence.
Homogeneity means that the statistical properties of
the space domain of turbulence are independent of
position in that domain, and isotropy means that the
properties are independent of direction. The method of
obtaining the statistical moments of the velocity field
from the Navier-Stokes equation leads to the same type
of closure problem faced in stochastic dynamic prediction
of the primitive equations. The assumption of homoge-
neity, however, leads to a tremendous simplification of the
equations, that is, the mean velocity can be considered
zero; After obtaining the Fourier transformation of the
Navier-Stokes equation, the variance of the Fourier
components of the velocity field is then proportional to
the turbulent energy, and the third moments are necessary
for nonlinear transfer of energy between eddies. A closure
technique that was considered for many years was the
quasi-normal theory of Millionshtchikov (1941) and
Proudman and Reid (1954). This scheme retains a
prognostic equation for the third moments but relates
fourth moments to second moments by assuming a
Gaussian distribution. It was pointed out by Ogura
(1963) that negative values of energy are obtained from
this closure scheme after a finite length of time. The
effect of dropping third moments is not as drastic as in
homogeneous turbulence, but its effect must be investi-
gated.
RELATION OF CLOSURE TO ENERGETICS
Insight into the effect of the third moments is gained
from considering their effect on the energetics. The
effect is only derived for the kinetic energy equation in
a barotropic model, but the results for a baroclinic model
can similarly be evaluated and are only stated below.
Considering only the third moment term, one can write
The present day analytical ’ theories of
Again, considering only this term, one obtains
where the right-hand side would fall under the notation
AMU as defined in the original energy equation develop-
ment of the previous section. With the technique used in
discussing term “9” in eq (29), it is readily seen that the
third moment adds to the definition of C(UKE, UKZ).
If the variables et were included, the third moments in
the temperature equation would become part of the term
labeled ASU and would be seen to contribute to C(UA,
868 MONTHLY WEATHER REVIEW V O I . 99, No. 11
UK) and C( UA,, UAE). Consideration of each distinct
triplet also reveals that nonlinear exchange of uncertain
energy between wave numbers is now possible where it
was not possible when the third moments were dropped.
I n the turbulence theory, this latter energy transfer is
the nonlinear turbulent energy transfer between eddies.
The important point here is that even though Epstein’s
simple closure scheme assures conservation of total
energy (certain and uncertain) in those physical systems
where it should be conserved, the neglect of the third
moments yields the result that the correct allocation
of energy to specific “boxes” is not assured. The conse-
quences of that misappropriation must be investigated.
One may summarize the effects of the third moments
by noting that they do not affect the conservation of
energy. They are not explicitly involved with any conver-
sion of certain t o uncertain energy, but they are explicitly
involved in all transfers between uncertain energy com-
ponents. Of course, if a nonlinear transfer between un-
certain modes is incorrect, the subsequent exchange
between those modes and certain modes will be affected.
Evidently, one must consider numerical calculations which
contain third moments to judge their effect on the other
statistical properties.
QUASI-NORMAL CLOSURE AND CONGRUENCY
The derivation of a prognostic equation for third
moments proceeds in the same manner as that used in
deriving such an equation for the second moments (cf.
Fleming 1970). Pollowing the general form for the primi-
tive equations given in eq (2), the prognostic equation for
the third moments is given by
+ik l= [a j p q (P p r.% I g +k 7klP--aPqak 1 +b l p g )
P, P
+a k ~q (k J T j l c +b % T~ l P -‘P ~u 5 1 +x ~ l p q )
+a l P q (& r j k q +& T9kP-flPqu,k+hj.%pg)]
where X refers to fourth moment about the instantaneous
mean. One must now make assumptions about the fourth
moments or derive an equation for them and make
assumptions about the fifth moments. The closure scheme
considered here is analogous to the quasi-normal theory
used in turbulence theory. Assuming that the initial values
are multivariate normally distributed, an expression for
the general fourth moment about the mean, k k t p P , can be
found by generating the first, second, third, and fourth
moments from the moment generating function. The
fourth moment about the mean can then be expressed as
~k I p q =‘k l apq+ a k p O I q +~k q a l P * (74)
In a multivariate normal distribution, the third mo-
ments are identically zero and the fourth moments are
expressible as products of second moments. If this expres-
sion is used in eq (73), the system of equations is closed,
but since it is used in the third moment equation, one can
see how the “quasi” term originates. Inserting eq (74)
into (73) yields
A
+ a k ~q (p ~ 711q+11~ 7 j l p f a j p a t q + ‘Jjq ‘Jip)
+ a l P p (k J rikq+& T j k P + ~5 p ‘k g + a j q a k p
)I -x (b 5 P T k l ,f b.%PT~ l p f b l p T j k p ) * (75) P
Before evaluating this closure scheme we ask, “How
can the merit of any closure scheme be verified?” We are
concerned with more than just the verification of the
original N dependent variables. We are interested in the
evolution of the joint probability distribution of all those
variables; less demanding, we wish to know if the dynamic
coupling between those variables and the measure of un-
certainty associated with them is “representative” based
upon all the associated uncertainties which have been
discussed. Since arriving at the true distribution requires
solving the infinite set of moment equations, we settle for
an approximation of the distribution. This is achieved by
considering a great many samples in a Monte Carlo cal-
culation as used by Epstein (1969b).
In the Monte Carlo method used here, a sample is
selected from the ensemble of points in phase space. The
points are obtained at random from a multivariate normal
distribution with a preassigned mean and variance for
each variable corresponding to the same mean and vari-
ance used in a stochastic dynamic method. The deter-
ministic trajectories in phase space are then determined
for each of the points. The results of each of the deter-
ministic forecasts are accumulated in an array at specific
times and any desired statistical quantity such as the
mean, variance, etc., is computed from that array. When
a great many deterministic trajectories are computed out
to time t , the mean of those values at time t1 is a good
estimate of the expected value. Thus, the collective
behavior of a great many samples in a Monte Carlo
method gives the “norm” by which a closure scheme can
be measured.
Obtaining a good approximation t o higher order
moments may demand an extremely large sample size;
that size will be further compounded by the inclusion of
uncertain forcing parameters. Nevertheless, we shall define
lLcongruencyl’ in an order of moments to mean that they
are “superposable so as to coincide throughout’’ with
Monte Carlo moments of the same order. Ideally, all the
moments of a closure scheme should be congruent with
their Monte Carlo counterparts. As the forecast period is
extended, however, the effect of dropping the higher order
moments d l eventually destroy that conguency-even
November 1971 Rex J. Fleming 869
10 II 12 13 14 15 16 17 18 19
Days
FIGURE 19.-Comparison of quasi-normal closure by showing
cosine components of wave 4.
down to the lowest order moment, the mean. Congruency,
at least in the lowest order, must be maintained in a
stochastic dynamic forecast. It makes little difference how
large the uncertainty becomes because each forecast user
must determine his own limit of uncertainty. There is
much useful information in knowing how uncertain a fore-
cast is-plans can be formed from this information-but
the integrity of the mean must be maintained or that
information may be useless.
The stochastic dynamic calculations that include third
moments must be restricted to a barotropic model because
of the vast amount of storage required. Such a calculation
is shown in figure 19. The +I had values ranging from 0.00
to 0.06. The initial variance of each 4, was l.0X10-5.
Waves 2, 4, and 6 were considered in this case, and the
value of the cosine component of wave 4 is plotted as it
evolved in three types of forecasts: the deterministic,
the stochastic dynamic with quasi-normal closure, and the
Monte Carlo.
It is seen that at 10 days, the quasi-normal closure
scheme is still congruent in the mean, being very close t o
the Monte Carlo; and that both differ from the determi-
nistic solution. After l l )( days, the first negative variance
appeared in the quasi-normal calculation. At 187i days
the calculation begins to “blowup”. The reason for the
blowup is not the negative variance itself, but rather that
energy is flowing from the uncertain boxes to the certain
boxes at an ever-increasing rate. The sum of the energy
is still conserved, but only because the uncertain energy
is becoming more and more negative with the certain
energy becoming more and more positive, eventually
resulting in an overly strong zonal flow that violates
linear computational stability. The initial kinetic energy
(certain and uncertain) was 2,342 units and the maximum
zonal wind in the southern part of the channel was 34 m/s.
After 19 days, the total kinetic energy was still remarkably
conserved, being 2,340.4 units. However, UK, was
-20,261 units with the large positive balance spread over
Kz, K,, UK,, and the zonal wind in the southern part of
the channel had reached 91 m/s.
Orszag (1970) shows that the failure of the quasi-normal
closure in homogeneous isotropic turbulence calculations
can be attributed to the theory’s lack of a dynamically-
Q-NtEddy-Damped ___
9-N+Eddy-Damped ____
l1140l
1.2 (1/20)
n 0.4
0.2
-
-
DAYS
FIGURE 20.-Plot of moment coefficient of skewness of $8.
200 I I I I I I I I
Monte Carlo ____ ( 500)
Deferministic -
IO I I 12 13 14 15 16 17 18 19
DAYS
FIGURE 21.-Comparison of stochastic dynamic closure with no
third moments by showing cosine components of wave 4.
determined relaxation time-so that no slightly perturbed
steady state can be approached. He shows numerical
calculations of a weakly driven model system that begins
with an initial equilibrium ensemble which should relax to a
new stationary state. However, the system overshoots its
new stationary state, and the energy density of some of
the wave numbers is seen t o oscillate with increasing
amplitude as time increases. A very similar effect is seen
in figure 20 where the “moment coefficient of skewness”
(m.c.s.) of one of the components is shown. This computa-
tion is based upon the results of the same case shown in
figure 19. One sees that the moment coefficient of skewness
increases without bound.
The results of the simpler scheme of dropping third
moments for the same case are shown in figure 21. The
figure shows that this scheme is realizable; in fact, other
cases showed no tendency for negative energy out to 6
weeks. But this simple closure scheme does not maintain
congruency for as long as might be desired. Deviations
870 MONTHLY WEATHER REVIEW V O l . 99, No. 11
IO II 12 13 14 15 16 17 18 19
DAYS
FIGURE 22.-Comparison of eddy-damped closure by showing
cosine components of wave 4.
from the Monte Carlo solution become important after
13 days.
EDDY-DAMPED CLOSURE
In dealing with analytical theories of homogeneous
isotropic turbulence, a number of authors have tried to
incorporate a damping effect on the third moments while
retaining the quasi-normal assumption. Early work in this
area was done by Edwards (1964) and more recently by
Orszag (1970) and Leith (1971). The growth of the third
moments is thus stymied by a damping term that is
intended to simulate the “natural” damping of triple
correlations that would take place if all the higher ordered
moments were included. What this damping factor should
be, or what it should be a function of, has not been analyt-
ically determined.
It will be important to know if the simplest type of
damping will work in the stochastic calculations, because
if so, then a more complicated form may ultimately be
avallable. With this motivation, consider the eddy-
damped form of the quasi-normal forecast equation of the
third moments, eq (75), to be
;jkl= . . . +[quashorma1 terms] . . . --krjkl (76)
where k is a positive constant. The damping coefficient, k,
could be a function of time, energy density, other statis-
tical quantities, or any number of things but for our
purpose an absolute constant is preferred. One is not
completely in the dark as to the numerical value of k since
it has dimensions of (time)-’, and one would expect that
the time scale of the dynamic features should be reflected
in its value.
Results of calculations using eq (76) are shown in figures
22 and 23. The initial conditions were exactly those used in
discussing the quasi-normal approximation, the Monte
Carlo method, and the simple stochastic dynamic closure
method earlier in this section. I n figure 22, the damping
coefficient is 0.1. Clearly, the eddy-damped theory is more
congruent than the simpler closure scheme. DifFerent
values of the damping coefficient are used and the results
compared in figure 23. The calculations are in very close
Monte Corb 1500) ----- 0-N f Eddy- Domprd IIIBOI . . . . . ... .
0-N t Eddy-Dompsd 11/40) D o e o o
0-N + Eddy-Dompsd 11/201 -
150
50
6
2 0
- 50
-100
IO I 1 12 13 14 15 16 17 I8 19 20
DAYS
FIGURE 23.-Comparison of various eddy-damped closures by
showing cosine components of wave 4.
agreement d t h the Monte Carlo results out to 18 days.
The extremely good results obtained with this simple
damping factor certainly indicate that the technique is
applicable to stochastic dynamic predictions. The addi-
tional calculations required for this method may not give
sufficient improvement in the early stages of a forecast to
warrant its use. However, recent turbulence calculations
by Leith (1971) indicate that the effect of third moments
being damped can be incorporated without calculating
third moments at every time step. This approach, called
the eddy-damped Markovian approximation” is also
discussed by Qrszag (1971). The saving of computer time
resulting from this method may make the eddy-damped
quasi-normal closure approach to the stochastic dynamic
equations economically feasible.
L l
CLOSURE ON A DISPERSIVE SKSTEM
In his study of predictability, Lorenz (1965) assumedf,
the Coriolis parameter, to be constant. Actually, on a
global scale, f is a function of y and this dependence can
be expressed as /3=dfldy. Prequently, /3 is considered a
constant with a numerical value appropriate for 45’ lati-
tude. The inclusion of the 8-effect causes a significant
change in the direction and amplitude of the phase speed
of the ultra-long waves; however, the calculations pre-
sented in part II indicate that Lorenz’ assumption was
correct; that is, the /.%effect is not important to predict-
ability.
The &effect gives rise to vortex foTces or restoring
forces which act to restore displaced fliiid colbmns and
move waves tvestmard. The waves (Rosshy waves) are
dispersive; that is, the phase speeds are a function of the
wavelength. A simple perturbation analysis yields the
phase speed of these waves as u-p/nz where u is the
constant mean flow and n is the wave number. This shows
that Rossby waves always travel westward relative to
the mean flow and are highly dispersive, but for the
higher wave numbers the 8-term is less effective.
An interesting paper by Benney and Newel1 (1969)
discusses the closure problem that results when trying t o
solve the statistical initial value problem for a system of
waves. The system is assumed to be spatially homogeneous
and the problem is to close the hierarchy of equations for
the statistical moments. These authors claim that in a
November 1971 Rex J. Fleming 871
I l l
FIGURE 24.-Comparison of closure schemes with &effect and
waves 2, 4, and 6 by showing cosine components of wave 4.
100
50
72 O
X
* -50
-100
-150
I I II
0 1 2 3 4 5 6 7 8 9
DAYS
FIGURE 25.-Comparison of closure schemes with &effect and
waves 4, 8, and 12 by showing cosine components of wave.8.
weakly nonlinear system of dispersive waves there is a
“natural” closure. The dispersive nature of the waves
causes a decoupling of triple correlations and there is an
approach to the Gaussian state, that is, a multivariate
normal probability distribution. Benney and Newell pos-
tulate that a second process is at work regenerating higher
moments by nonlinear couplings (due to the nonlinear
inertial forces transferring energy between wave numbers) ,
but that this second process occurs over longer time scales.
It is of interest to test these ideas on a stochastic dynamic
system with dispersive waves, but where homogeneity is
not assumed.
The p-effect, introduced with the same initial conditions
that were used in the quasi-normal calculation (cf. fig.
19), are used to give the results shown in figure 24. No
negative energy occurs and figure 24 clearly shows that
the quasi-normal closure method does not blow up. More-
over, the same calculation using the eddy-damped method
with coefficient k=O. 1 agrees with the quasi-normal calcu-
lation almost exactly out to 16 days. Thus, the damping
term is damping “nothing” or at least a very small
quantity. A comparison of figures 19 and 24 reveals that
the P-effect has introduced a fast time scale into the
system, as expected with no artificial divergence.
The above results are for waves 2 , 4, and 6 which are
highly dependent upon 8. Higher wave numbers are not
so dependent. Figure 25 shows calculations with the same
initial values as before but for the waves 4. 8. and 12
The amplitude of the cosine term of wave 8 is shown. Here,
one sees that the quasi-normal theory blows up after 8
days-the p-effect is not strong enough in these waves to
give a natural closure. Note also that the results from the
quasi-normal closure and the eddy-damped closure cliff er
significantly from the deterministic solution after 6
days-both in amplitude and phase.
These calculations show that the discussion of Benney
and Newell applies equally well to the stochastic dynamic
prediction of the ultra-long waves. However, the practical
application of numerical weather prediction must include
much smaller scales or higher wave numbers than those
where the P-effect is pronounced, so that the discussion of
Benney and Newell is irrelevant for our purposes.
8. CONCLUSIONS
This paper has merely been an introduction into a few
of the important concepts associated with the stochastic
dynamic method. The formally exact stochastic approach
to a nonlinear problem (with associated uncertainties)
leads to an infinite unclosed set of moment equations.
Epstein (1969b) has closed this system of equations in a
very simple way. Calculations given here suggest that this
closure may be sufficient for some applications. However,
evidence has been presented that a better closure scheme
should be used in extended time integrations.
The stochastic dynamic equations which neglect third
and higher order moments conserve the ensemble energy
in those cases where the energy should be deterministically
conserved. However, the neglect of third moments was
shown to lead to a misappropriation of energy between
the uncertain energy components. The quasi-normal
closure scheme fails to give realizable results in the stochas-
tic dynamic equations just as it fails in homogeneous iso-
tropic turbulence theory. The eddy-damped closure
method gave extremely good results out to 18 days as
compared with Monte Carlo calculations.
The stochastic dynamic approach treats the question of
predictability or growth of uncertainty in an explicit
manner. The use of the equations in considering the
predictability of various scales of features is given in part
11. We emphasize here that more can be provided than the
usual notion of an average value of predictability. More
important, the believability of the atmospheric variables
(which depends upon all the dynamic processes a t work and
all the associated uncertainties) can ,be localized in space
and time. In addition, the nature and amount of this
growth can be viewed through the energetics of the
stochastic dynamic equation set. Thus, a full energy
description of the growth of error has been presented, if
only in simple models.
There are two important implications concerning the
stochastic approach that have been brought into focus in
this study: (1) the method has value, and (2) the method
is possible. We amplify on these in reverse order.
The stochastic equations can be solved. The computer
power needed is monumental but not infinite. The infinite
set of moments required in principle is not necessary in
~, -, .____ __. practice. It appears that the dynamical equations them-
872 MONTHLY WEATHER REVIEW Voi. 99, No. 11
selves damp the higher order moments. This was seen in
the extremely accurate results obtained by a simple damp-
ing of third moments. This is obviously related to the fact
that the original nonlinear equations were quadratic. The
author has derived the stochastic equations for a general
cubic deterministic equation set and the closure problem is
more complex. The instability and energy transfer mech-
anisms in the quadratic hydrodynamic set involve triple
correlations of deterministic variables. One might specu-
late that the absence of any known significant interaction
of four OT more variables accounts for the success of the
closure used here.
The value of the stochastic approach lies not in its fore-
cast of the mean of the dependent variables. In fact, in
the early stages of the forecast, there will not be too much
difference between these and the deterministically com-
puted values of these variables. The value of the method
shown here is its ability to forecast the uncertainty of the
variables. We have shown that the stochastic equation
set including third moments is energetically consistent.
This means that the error growth is realistic and that the
believability of each variable will evolve according to the
dynamical situation at hand.
ACKNOWLEDGMENTS
The author is grateful to Professor Edward S. Epstein and Dr.
C. E. Leith for their helpful suggestions.
Financial support was provided by the United States Air Force
and by National Science Foundation grants GP-4385 and GS-19248.
Computations were performed a t the National Center for Atmos-
pheric Research. The author thanks Mrs. Mary Daigle of the
National Meteorological Center, National Weather Service, NOAA,
for typing the manuscript.
REFERENCES
Benny, D. J., and Newell, Alan C., “Random Wave Closures,”
Studies in Applied Mathematics, Vol. 48, No. 1, Massachusetts
Institute of Technology Press, Cambridge, Mar. 1969, pp. 29-53.
Dutton, John A., and Johnson, Donald R., “The Theory of Available
Potential Energy and a Variational Approach to Atmospheric
Energetics,’’ Advances in Geophysics, Vol. 12, Academic Press,
New York, N.Y., 1967, pp. 333-436.
Edwards, S. F., “The Statistical Dynamics of Homogeneous
Turbulence,” Journal of Fluid Mechanics, Vol. 18, No. 2, Cam-
bridge University Press, London, England, Feb. 1964, pp. 239-
273.
Epstein, Edward S., “The Role of Initial Uncertainties in Pre-
diction,” Journal of Applied Meteorology, Vol. 8, No. 2, Apr.
1 9 6 9 ~~ pp. 190-198.
Epstein, Edward S., “Stochastic Dynamic Prediction,” Tellus,
Vol. 21, No. 6, Stockholm, Sweden, 1969b, pp. 739-757.
Fleming, Rex J., “Concepts and Implications of Stochastic Dynamic
Prediction,” N C A R Cooperative Thesis No. 22, University of
Michigan, Ann Arbor, and Laboratory of Atmospheric Science,
National Center for Atmospheric Research, Boulder, Colo.,
Frieberger, Walter, and Grenander, Ulf, “On the Formulation of
Statistical Meteorology,” Review of the International Statistical
Institute, Vol. 33, No. 1, Institute of Mathematical Statistics,
University of Stockholm, Sweden, 1965, pp. 59-86.
Gleeson, Thomas A., “A Causal Relation for Probabilities in
Synoptic Meteorology,’J Journal of Applied Meteorology, Vol. 5,
No. 3, June 1966, pp. 365-368.
n m , 171 pp.
Hide, Raymond, “Some Experiments on Thermal Convection in a
Rotating Liquid,” Quarterly Journal of the Royal Meteorological
Society, Vol. 79, No. 339, London, England, Jan. 119513, p. 161.
Hildebrand, Francis B., Advanced Calculusfor Applications, Prentice-
Hall, Inc., Englewood Cliffs, N.J., 1965, 646 pp.
Leith, C. E., “Numerical Simulation of Turbulent Flow,” Proper-
ties of Matter Under Unusual Conditions, John Wiley & Sons,
Inc., New York, N.Y., 1969, pp. 267-271.
Leith, C. E., “Atmospheyic Predictability and Two-Dipension81
Turbulence,” Journal of the Atmospheric Sciences, Vol. 28, No. 2,
Mar. nW1, pp. 145-161.
Lorenz, Edward N., “Energy and Numerical Weather Prediction,”
Tellus, Vol. 12, No. 4, Stockholm, Sweden, Nov. 1960, pp. 364-373.
Lorenz, Edward N., “Deterministic Nonperiodic Flow,” Journal of
the Atmospheric Sciences, Vol. 20, No. 2, Mar. n963a, pp. 130-141.
Lorenz, Edward N., “The Mechanics of Vacillation,’’ Journal of the
Atmospheric Sciences, Vol. 20, No. 5, Sept. B963b, pp. 448-464.
Lorenz, Edward N., “A Study of the Predictability of a 28-Variable
Atmospheric Model,” Tellus, Vol. 17, No. 3, Aug. 1965, pp. 321-
333.
Lorenz, Edward N., “The Predictability of a Flow Which Possesses
Many Scales of Motion,” Tellus, Vol. 21, No. 3, 1989, pp. 289-307.
Millionshtchikov, M., “On the Theory of Homogeneous Isotropic
Turbulence,” Compte Rendus (Doklady) de l’dcademie Des Sciences
de l’U.R.S.S., Vol. 32, Moscow, U.S.S.R., 1941, pp. 615-618.
Ogura, Yoshimitsu, “A Consequcnce of the Zero-Fourth-Cumulant
Approximation in the Decay of Isotropic Turbulence,” Journal
of Fluid Mechanics, Vol. 16, No. 1, Cambridge University Press,
London, England, May 1963, pp. 33-40.
Orszag, Steven A., “Analytical Theories of Turbulence,” Journal
of Fluid Mechanics, Vol. 41, No. 2, Cambridge University Press
London, England, Apr. 1970, pp. 363-386.
Proudman, Ian, and Reid, W. H., “On the Decay of a Normally
Distributed and Homogeneous Turbulent Velocity Field, ”
Philosophical Transactions of the Royal Society of London, Ser.
A, Vol. 247, No. 926, England, Nov. 4, 1954, pp. 163-189.
Tatarskii, V. I., “The Use of the Dynamical Equations in the Proba-
bility Forecast of the Pressure Field,” Atmospheric and Oceanic
Physics, Vol. 5, No. 3, The American Geophysical Union, Wash-
ington, D.C. Mar. 1969, pp. 162-164 (translation of ‘%pol’-
sovanie Dinamicheskikh Uravnenii Pri Veroiatnostnom Prognoze
Baricheskogo Polia,” Akademiia Nauk SSSR, Izvestiia, Fizika
Atmosfery i Okeana, Vol. 5, No. 3, Leningrad, U.S.S.R. Mar.
Taylor, Goeffrey Ingram, “Diffusion by Continuous Movement,”
Proceedings of the London Mathematical Society, Vol. 20, London
Mathematical Society, London, England, 1921, pp. 196-212.
Thompson, Philip Duncan, “Uncertainty of Initial State as a
Factor in the Predictability of Large-Scale Atmospheric Flow
Patterns,” Tellus, Vol. 9, No. 3, Stockholm, Sweden, Aug. 1957,
Wiin-Nielsen, Aksel, “On Short- and Long-Term Variations in
Quasi-Barotropic Flow,” Monthly Weather Review, Vol. 89,
Wiin-Nielsen, Aksel, “Some New Observational Studies of Energy
and Energy Transformations in the Atmosphere,” WorM: Mete-
orological Organization Technical Note No. 66, Geneva, Switzer-
land, 1966, pp. 177-202.
Wiin-Nielsen, Aksel, “On the Annual Variation and Spectral Dis-
tribution of Atmospheric Energy,” Tellus, Vol. 19, No. 4, Stock-
holm, Sweden, Nov. 1967, pp. 540-559.
Wiin-Nielsen, Aksel, “On the Intensity of the General Circulation
of the Atmosphere,” Reviews of Geophysics, Vol. 6, No. 4, NOV.
1968, pp. 559-579.
1969, pp. 293-297.
pp, 275-295.
NO. 11, NO^. 1961, pp. 461-476.
[Received January 7 , 1971; revised May 24, 19711