T H E UN I V E R S IT Y O F M I C H I G AN COLLEGE OF ENGINEERING Department of Atmospheric and Oceanic Science Technical Report STOCHASTIC DYNAMIC PREDICTION USING ATMOSPHERIC DATA Eric J. Pitcher Edward S. Epstein Stanley J. Jacobs Project Directors DRDA Project 011554 supported by: NATIONAL SCIENCE FOUNDATION GRANT NO. GA-35771X1. WASHINGTON, D.C. administered through: DIVISION OF RESEARCH DEVELOPMENT AND ADMINISTRATION ANN ARBOR July 1974

ACKNOWLEDGMENTS I wish to thank Professors Edward S. Epstein, Karl T. Hecht, Aksel C. Wiin-Nielsen and Dr. Thomas W. Schlatter for serving as members of the Doctoral Committee, It is indeed a pleasure to acknowledge the encouragement and constructive criticism of Professor Epstein, and I have derived considerable benefit from many discussions with him. The contribution of Professor Wiin-Nielsen to my education through excellent teaching as well as informal conversations has been substantial, and his interest in this work is very much appreciated. I am indebted to Dr. Schlatter and the NCAR staff who made available the necessary support facilities. I wish to thank my wife, Peggie, and Ms. Debra Lurkins for the hours and effort they spent in typing this manuscript. This study has received support from the National Science Foundation under Grant GA-35771X1. Acknowledgment is made to the National Center for Atmospheric Research, which is sponsored by the National Science Foundation, for computer time in this research. iii

TABLE OF CONTENTS ACKNOWLEDGMENTS i i i LIST OF TABLES vi LIST OF FIGURES vii ABSTRACT x CHAPTER I. INTRODUCTION 1.1 Preliminary Remarks 1 1.2 Formulation of Stochastic Dynamics 8 1.2.1 Stochastic Dynamic Equations 8 1.2.2 Moment Solutions via a Lagrangian Technique 13 II. APPLICATION TO A SIMPLE SYSTEM 2.1 The Minimum Equations 16 2.2 Elementary Properties of Elliptic Functions 24 2.3 Solutions of the Minimum Equations 26 2.4 The Stochastic Equations 33 2.5 Results 35 III. FORMULATION OF STOCHASTIC FORECAST MODEL 3.1 Motivation for Stochastic Dynamic Forecasting 48 3.1.1 Determination of an Atmospheric State from Observations 54 3.1.2 Stochastic Analysis 59 3.2 External Model Uncertainty 62 3.3 Equivalent Barotropic Model 66 3.4 Spectral Formulation 70 3.5 The Question of Model Uncertainty 73 3.5.1 Parameterization of Model Uncertainty 74 3.5.2 Time Dependence of all) and (2) 76 3.6 Stochastic Dynamic System 80 IV. RESULTS OF FORECAST MODEL 87 4.1 Neglect of External Error Growth 90 4.2 Random Forcing of the First Kind 107 4.3 Random Forcing of the Second Kind 119 iv

TABLE OF CONTENTS (concluded) V. SUMMARY AND OUTLOOK 127 APPENDIX A. PROOF OF STOCHASTIC ANALYSIS ALGORITHM 131 B. EVALUATION OF INTERACTION COEFFICIENTS 136 C. RELATIONS BETWEEN REAL AND COMPLEX QUANTITIES 138 D. NUMERICAL METHOD FOR TIME INTEGRATION 141 E. FORMULATION FOR RANDOM FORCING OF THE SECOND KIND 143 REFERENCES 147

LIST OF TABLES Table Page 2.1 Maximum Value of Each Component 25 2.2 Relevant Constants Appearing in Elliptic 31 Function Solutions vi

LIST OF FIGURES Figure Page 2.1 Comparison of Liouville solution for P2 and p3 with stochastic dynamic neglecting third moments. 38 2.2 Comparison of Liouville solution for o12 with stochastic dynamic neglecting third moments. 40 2.3 Comparison of Liouville solution for a 12 with stochastic dynamic neglecting third moments. 41 2.4 Normalized triple moment T 123 43 2.5 Normalized triple moment T 333. 46 4.1 Typical distribution of station locations whose latitudes exceed 200N. 88 4.2 Standard deviation of least squares estimation of 500 mb height field, Day 0. Units are in meters. 91 4.3 Same as Fig. 4.2 but for Day 1. 91 4.4 Same as Fig. 4.2 but for Day 2. 91 4.5 Same as Fig. 4.2 but for Day 3. 91 4.6 Least squares estimation, Day 0. Units are in meters. 94 4.7 Forecast expected state of height field, Day 1. No random forcing included. 94 4.8 Least squares estimation, Day 1. 94 4.9 Stochastic analysis of expected values,Day 1. 94 4.10 Standard deviation of least squares estimation, Day 0. 98 4.11 Standard deviation of forecast expected state, Day 1. 98 4.12 Same as Fig. 4.10 but for Day 1. 98 vii

LIST OF FIGURES (continued) Figure Page 4.13 Standard deviation of analyzed expected state, Day 1. 98 4.14 Standard deviation of forecast expected state, Day 2. No random forcing included. 101 4.15 Standard deviation of least squares estimation, Day 2. 101 4.16 Standard deviation of analyzed expected state, Day 2. 101 4.17 Time sequence of A6. The dashed vertical lines represent the least squares estimates (expected value plus and minus one standard deviation) from observations on the day in question. The solid vertical lines represent the 24 hr stochastic forecasts, immediately to the left, and stochastic analyses, immediately to the right of the least squares estimates. 103 4.18 Time evolution of A 2. Notation is the same as that in Fig. 4.17. 103 4.19 Scatter diagram of apparent forecast error e vs a, the standard error as inferred from observational and forecast uncertainty; a = (cf + a) 1/2. Tne apparent error is defined as the magnitude of the difference between the least squares estimate and the forecast value. 106 4.20 Forecast expected state of height field, Day 1. Random forcing of the first kind. 112 4.21 Standard deviation of forecast expected state, Day 1. Random forcing of the first kind. 112 4.22 Standard deviation of analyzed expected state, Day 1. 112 4.23 Forecast expected state of height field, Day 2. Random forcing of the first kind. 114 viii

LIST OF FIGURES (concluded) Figure Page 4.24 Standard deviation of forecast expected state, Day 2. Random forcing of the first kind. 114 4.25 Standard deviation of analyzed expected state, Day 2. 114 4.26 Scatter diagram of apparent forecast error e vs a, the standard error as inferred from observational and forecast uncertainty; a = (CT + ) 12, with symbols defined in section 4.1. The apparent error is defined as the magnitude of the difference between the least squares estimate and forecast value at Day 2. Random forcing of the first kind. 115 4.27 Time evolution of AO. Random forcing of the first kind. Notation is the same as that in Fig. 4.17. 118 4.28 Time evolution of Al. Random forcing of the first kind. Notation is the same as that in Fig. 4.17. 118 4.29 Time evolution of A>. Random forcing of the second kind. Notation is the same as that in Fig. 4.17. 121 4.30 Scatter diagram of apparent forecast error e vs a, the standard error as inferred from observational and forecast uncertainty; a = (ao + a2)112, with symbols defined in section 4.1. The apparent error is defined as the magnitude of the difference between the least squares estimate and forecast value at Day 2. Random forcing of the second kind. 123 4.31 Stochastic analysis of expected height field, Day 3. 125 4.32 Forecast expected state of height field, Day 4. Random forcing of the second kind. 125 4.33 Least squares estimation, Day 4. 125 4.34 Stochastic analysis of expected height field, Day 4. 125 ix

ABSTRACT STOCHASTIC DYNAMIC PREDICTION USING ATMOSPHERIC DATA by Eric John Pitcher Chairman: Edward S. Epstein A procedure is presented whereby an evolving ensemble of possible fluid states may be predicted to any desired level of accuracy for a simple system of one wave interacting with a zonal flow. This is accomplished by utilizing the Lagrangian form of the Liouville equation for the conservation of probability. The lowest order moments of the solution are compared with the solution obtained by solving the corresponding moment equations, referred to as the stochastic dynamic equations, subject to the third moment discard approximation. This closure scheme is concluded to be a satisfactory one for stochastic dynamic forecasting on the order of a few days. A set of stochastic dynamic equations is derived for an equivalent barotropic model possessing 106 degrees of freedom. This model is used to make forecasts from atmospheric data. The output of each forecast is an estimate of the expected state of the atmosphere, and also the uncertainty associated with that estimate as measured by the variance information. In addition to the forecast of

the variance field, all of the covariances among the model parameters are predicted. This enables one to make use of Bayes' theorem from probability theory and utilize the covariance information associated with both the forecast and observed fields to arrive at an analysis. Because the basic physical model is not an exact atmospheric analogue, one must incorporate by some parameterization the growth of error during the forecast stage attributable to the model deficiencies. This aspect of the problem is rather difficult from both a mathematical and practical standpoint, the latter implied by the fact that the total computational problem must remain manageable. Two parameterizations are investigated and their relative merits assessed. The forecasts made are credible and the uncertainty fields revealing. It is true that those regions of maximum uncertainty in the forecast are similar to the initial uncertainty pattern, but there are significant changes in response to model dynamics, and the results indicate the feasibility of the approach. Implications for more complicated models are considered. The essential departure from tradition hinges around the fact that the stochastic framework allows us to investigate the evolution of an ensemble of solutions rather than just one. This is especially useful in predicting physical systems such as the atmosphere, where the initial state is imprecise and the dynamical model limited. xi

CHAPTER I INTRODUCTION 1.1 PRELIMINARY REMARKS One of the major goals of meteorology may be stated in rather broad but simple terms: to provide reasonably accurate predictions of weather events and to extend the prediction time over which these prognoses are useful. Several obstacles have appeared in the path leading toward a realization of this goal. Until relatively recent times, not the least of these difficulties was the absence of a device which could carry out the vast number of computations required in order to make a timely forecast. Yet, it was recognized that a satisfactory solution would not be possible unless the atmosphere were viewed as a system subject to physical laws which could be expressed in terms of a set of mathematical relations. An example of this is the pioneering effort of Richardson (1922). With the development of the earliest computers in the middle 1940's, the prospect for a viable solution to the meteorological prediction problem seemed brighter. The physical complexity of the models was of course limited by the computing capability, and this was one of the factors which prevented a rapid increase in forecast skill. Another very important source of forecast error —one which still persists and which will be of major interest in this study —was the error present in the initial state.

The latter was always recognized as a real hindrance to weather forecasting, be it by subjective or objective means. Analysis error was pointed out by Charney and Eliassen (1949), and considered explicitly by Thompson (1957) in one of the early "predictability" studies. In that paper an estimate of two days was obtained for the doubling time of mean square wind error (equivalently, error kinetic energy), which would arise from analysis error alone. The latter is closely related to the distance between observing stations. We shall have further to say about this study later in this section. Few would argue that our ability to observe the atmosphere must keep apace with our ability to devise and solve sophisticated and realistic models of the atmosphere, since the ultimate success of the latter hinges on the former. In recent years an ever growing interest in this aspect of the problem has attracted the attention of many investigators. The effect of initial error, which may be considered observational and/or analysis error, on prediction error generally falls within the purview of predictability studies. These experiments typically study the evolution of two states of an atmospheric model, where one state differs slightly from the other initially. The subsequent "error" growth is investigated in terms of various criteria: dependence on grid spacing, whether initial error is present in (say) the temperature or wind

field, size of initial errors, to name but a few. Examples of some of these experiments are Smagorinsky and Miyakoda (1969), Jastrow and Halem (1970), Williamson and Kasahara (1971), Wellck et al. (1971) and Williamson (1973). These investigations yield a doubling time for error growth of about two or three days. Of equal importance is the aspect of updating the numerical model with new data at regular or irregular time intervals, and the question of four-dimensional data assimilation has spurred a flurry of activity. For various reasons, not the least of which is economic, future atmospheric observing systems will be composed of several subsystems (Kasahara, 1972): satellites, constant-level balloons, ocean buoys, conventional methods, etc. The evaluation of these subsystems has lead to several studies commonly referred to as Observing Systems Simulation Experiments (OSSE). Demanding significant attention has been the procedure by which new data are inserted into the numerical model. These experiments usually proceed in the following way. A model integration is made and this is considered, by definition, to be the "true atmosphere." Errors are added to the initial true state and another integration performed, only this time observations from the true atmosphere are inserted occasionally into what now may be considered a prediction model. (These observations may

or may not be subject to errors depending on the purpose of the particular experiment.) These investigations serve to assess forecast error in relation to the frequency of updates and the accuracy of the observations. The methods of inserting new data may be divided into two broad classifications. In the case where updating involves temperature and/or wind information, then direct substitution is certainly the simplest method. A predicted quantity is replaced by an observed value from the control model atmosphere (see, for example, Charney et al., 1969; Williamson and Kasahara, loc. cit.). It would seem reasonable that the new information at the "updating" point could also be used advantageously in adjusting values near this point. This has been the subject of optimum data analysis, and various approaches have been explored. The method of "optimum" interpolation (Gandin, 1963) utilizes statistics of climatological anomalies to interpolate meteorological variables. The calculation is made subject to a minimization of mean square error of interpolation. Because short-range forecasts have errors not too much larger than typical analysis errors, it would seem appropriate that forecasts should play a role in analysis other than just serving to fill in gaps where there are no data. Rutherford (1972) has proposed such a procedure in which climatological statistics of apparent forecast error are utilized in data assimilation. Tadjbakhsh (1969), Sasaki

(1970) and Petersen (1968) discuss various aspects of optimum data assimilation in which a dynamic equation also enters as an essential ingredient. A recurring aspect of the whole problem of data assimilation and prediction is the uncertainties which are ever present. Observing the atmosphere at a finite number of points precludes a complete specification at some instant. At best, we can only make a judgment about the true state —a judgment which is subject to uncertainty. That only a probabilistic statement can be made was clearly demonstrated by Thompson (1957) in a related discussion on the extent to which the atmosphere is deterministic. Therein the argument was made that, although the deterministic hypothesis is a plausible one, in which two identical atmospheric states would evolve in exactly the same manner given equivalent absorption of heat energy, it is not possible to know if two states are in fact identical because of the coarseness of the observing system. Thus two ostensibly "equal" states may subsequently diverge only because the real differences were initially undetected by the limited observing capability. Consequently, the essential limit of predictability depends upon the average separation between observations. Further elucidation of this point has been given in the theoretical and observational studies of Lorenz (1969a, 1969b). Not only the recognition that imperfect observations

and forecasts are a fact, but the desire to come to grips with this problem in mathematical terms has led to formal procedures for incorporating these uncertainties in the forecast-analysis problem. One such approach has been provided by Epstein (1969b). This method treats the atmosphere as a physical system subject to the well known laws of hydrodynamics, but at the same time subject to uncertainty in its specification from observations. The method, referred to as stochastic dynamic prediction, provides a forecast of the "best" state in the sense of least mean square error, and the uncertainty associated with that state. In this approach then, the total forecast information includes a prediction of the degree of confidence that one can place in the best estimates of the meteorological variables. In the light of previous remarks, this feature has obvious implications for optimum data assimilation. The stochastic dynamic method is rather general and involves the solution of inhomogeneous anisotropic equations of statistical hydrodynamics. Because of the nonlinear nature of the problem, a closure approximation (cf. section 1.2) is necessary as with turbulence models. This closure problem has been studied by Fleming (1971a), with the result that short-term integrations appear possible without the necessity for the retention of statistical moments beyond the second. The method has been used for the study of predictability (Fleming, 1971b), and other

investigations have demonstrated its utility (see, for example, Epstein, 1971; Epstein and Fleming, 1971). Tests of the stochastic dynamic method to date have been based primarily on synthetic data. These simulation experiments have contributed to an understanding of many aspects of the method and demonstrated its feasibility. The potential of the method suggests further experiments, but now applied to real atmospheric data. One such recent investigation is that of Knudsen (1973). In that study a relatively simple model employs atmospheric data and deals with some of the numerical aspects peculiar to the method. To conclude this section we outline the scope of the present study. The remainder of this chapter provides much of the mathematical formalism which will be required. An analytical study is presented in Chapter II. This is a study of the simple three-component system given by Lorenz (1960), and having the pleasing feature of analytic solutions. As will become apparent, this will be an aid in obtaining closed-form solutions of the corresponding stochastic equations. It is then possible to compare directly these solutions and the solutions obtained via the usual closure approximation. The major contribution of this study, however, is contained in succeeding chapters in which our goal is to make hemispheric stochastic forecasts of the 500 mb flow using real atmospheric data.

We shall make use of the vast amount of information in the stochastic forecast and utilize a procedure for optimum data assimilation in the manner proposed by Epstein and Pitcher (1972). It might be useful to point out that, with section 1.2 as backgound, it is possible to read Chapter III essentially independent of Chapter II. The results from the latter required in the former are simply quoted. 1.2 FORMULATION OF STOCHASTIC DYNAMICS We have already given, in somewhat general terms, the rationale behind the applicability of the stochastic dynamic method to atmospheric prediction. Further amplification of the points discussed earlier is postponed until Chapter III, where we shall consider more fully the relevance of the method to the problem of weather forecasting. This section discusses the mathematical details of the stochastic approach as applied to forced dissipative fluid systems in general. 1.2.1 Stochastic Dynamic Equations Much of the formalism which we shall require has already appeared in Epstein (1969b) and Fleming (1970, 1971a) We shall review this here for the sake of completeness. Consider a physical system specified by N parameters, xi. (i = 1,2,....N). Define a corresponding N-dimensional Euclidean space with coordinates xi. Each point in this

phase space defines a possible state of the system. The evolution of the system from one state to another is presumed to follow a dynamical law which may be reduced to a set of N differential equations of the form: i JtojKXjXK-I~ X+C;, t(1.1) where the aijk, b.i and c. are constants. The dot denotes differentiation with respect to time and the first summation on the right-hand side is carried out over all j and k. Let us postulate, for the moment, that our knowledge of the system at time to is given by a probability density p(X, to), where X is a vector composed of the xi. We need not be concerned with the details of ~ at present, noting only that the initial specification of the system is in terms of a cloud of phase points, with the density of points in a particular region proportional to the value of ~ there. By definition, (1.2) where dX = dxldx2....dxN and the integration is performed over all X. The fundamental problem is the prediction of ~(X,t)

10 for t>to, given 4(X,to). From (1.2) the total probability must always be one; equivalently, no members of the cloud of phase points may be created or destroyed. This is simply the conservation law for probability and may be expressed analytically by the Liouville* equation: ____A O ~ (1.3) at:' )xX This is the N-dimensional analogue of the ordinary threedimensional continuity equation expressing conservation of mass. In the present context xi is the ith velocity component of a phase point and is given by (1.1). The excessive computation required to solve (1.3) directly on (say) a grid of points in phase space for even a small value of N(10) is prohibitive. We seek a compromise. In short, the course to be adopted is to solve for the lowest order moments. Define the expectation operator for a quantity f(X) as follows: g [F{t)1 If ~4(rl(x~) t(X) 1(1.4) *Because of the canonical simplicity of the equations encountered in particle dynamics, ZExi/axi= 0, and the usual form of the time dependent Liouville equation is somewhat simpler than (1.3) (see, for example, Tolman, 1938). This simplification is deceptive however because a system of N particles defines a 6N-dimensional phase space, both positions and momenta being required in order to specify the state of the system completely.

11 We now introduce the following definitions:,- E (x;), (1. 5a) T [ — E[-s ] - E(zx. j) -/i.,q (1.5b) FZ~'ij_ (E(xsA;)(xK E (xi xj XV) - At q W S ~ ~tj~ tA; SK (1. 5c) where pi is the mean, aij the variance (i = j) and covariance (i Z j) and Tijk a third moment quantity, likewise centered about the mean. Taking the time derivative of (1.4), making use of (1.3) and the condition that xi -- 0O as xi —+ a, one can easily show that itE lX) F.E [1bl.(1.6) Consequently, pi= E(xi), and from (1.1) we may write ago We =j E ( Xj~) (xS~ e* c, i4C(1.7) or using (1.5b), Av,'Za;jK(.jA jik)iZbLisaiC (1.8)

12 A predictive equation for oj may be obtained by taking the time derivative of (1.5b), i.e., &i=E(*iXj iX ij)MHMj-AiJ. (1.9) By utilizing (1.1), (1.8), (1.5b) and (1.5c), equation (1.9) may be written in the following form: CT~ r - C6. ([a AKW(K.t +Atrit irA ZK. + C 4c KA A1k ZK)l Eblk riK + b eK (1. 10) The equations governing the first moments (means) involve second moments, while the second moment equations (1.10) require knowledge of third moments. This is an infinite sequence which results whenever the basic deterministic equations are nonlinear as regards statistical quantities. For the sake of argument, let us presume that the means and covariances are given, and, in addition, the Tik are zero at to. They will not, in general, remain ijk zero. To establish this fact, one need simply derive an equation for ~ijk' but this of course will introduce fourth moments. Clearly, in order to obtain a solution, certain assumptions about the nature of higher moments must be made —in effect, a closure approximation. Epstein (1969b) argued the case for dropping third

13 moments, and with a simple barotropic model showed that indeed this does give valid results for integrations over a few days. Fleming (1971a) has studied various closure schemes. His results indicate that third moment equations must be retained for long-term integrations ("10 days) in order that the integrity of the mean be maintained. In the present study we shall be primarily interested in short-range predictions, and, based in part on the work of the previous studies, shall adopt the closure scheme proposed by Epstein. Calculations supporting this approximation will be presented in Chapter II. 1.2.2 Moment Solutions via a Lagrangian Technique Previous considerations have led to a sequence of prediction equations for the moments of ~. We have seen that, in principle, a solution is possible only after invoking some closure approximation. In this subsection we propose an alternate method of solution which is conceptually simple, but yet does not require such a closure assumption. Even though computational requirements limit the method to fairly simple systems, one obvious advantage is that it allows a critical evaluation of various closures. Working with the continuity equation for probability in the form given by (1.3) is not satisfactory because of the vast computational requirement. The first step is to derive an alternate, but equivalent, form for the Liouville

14 equation. Consider a volume VO in phase space at time to. At a later time this entity will be transformed into a new volume V. Because each point of VO evolves according to the dynamical equations (1.1), the total number of phase points in V is invariant with time. Equivalently,the probability measure associated with V is constant, i.e., Byvqx0 - ae(1.11) By a fundamental theorem on differentials, AXL -~(x, / Xj E )dx;/ S) 1 v(1.12) or, in general, dX = JdXo where J is the determinant of the Jacobian matrix of the transformation, i.e, r3lya. [)l,(1.13) With this substitution we have ok~0X (cC)q~k I *(1.14) As (1.14) must be satisfied for an arbitrary Vo, we have as a conservation requirement, (1.15) This is the Lagrangian form of the Liouville equation, and gives the value of ~ along a given trajectory in phase space in terms of its value at to. With the previous result it is now possible to express (1.4) in the equivalent form:

15 EL;V~~X) g( o (1.16) The significance of the foregoing is that the integration, to be performed over X at time t, is transformed into an integration over X at to. This of itself is an important result because it obviates the necessity for obtaining a solution for c over a dense grid encompassing a large portion of phase space. Equation (1.16) may be solved by a numerical quadrature without ever computing O explicitly for t>to. Note, however, that E[f(X)] may be computed for an arbitrary time because of the presence of the time dependence in f(X). We intend to make use of (1.16) in the case where analytic solutions for X exist. By making a suitable choice for ~o (e.g., multivariate normal), one can utilize a particular quadrature formula ( to be discussed later), and evaluate (1.16) to any desired accuracy. It is not being advocated that the procedure just outlined is a viable alternative to the solution of the stochastic dynamic equations for moderately large systems. Rather, it provides a means for obtaining accurate solutions to systems with a few degrees of freedom, without having to face the closure problem. The next chapter illustrates the application of both procedures on such a simple system.

CHAPTER II APPLICATION TO A SIMPLE SYSTEM 2.1 THE MINIMUM EQUATIONS In subsection 1.2.2 we presented a procedure whereby, in principle, one may compute the evolution of all the statistical quantities of interest for an ensemble whose individual members are governed by some dynamical equation. Therein, a cautionary remark was made that such a calculation has been found to be feasible only when the number of degrees of freedom is not too large. In the case of a model required to make weather forecasts, such a calculation is computationally unwieldy and recourse must be made to a simpler algorithm such as the stochastic dynamic equations. (Even the latter in their present form, applied to a moderately sophisticated forecast model, present a very severe challenge to existing computers.) Nevertheless, with a sufficiently simple model, ideally one for which analytic solutions are readily available, the moments of an evolving ensemble may be computed with relative ease via the probability integral (1.16). These are then the correct solutions against which various approximate solutions, such as the stochastic dynamic with the neglect of third moment quantities, may be compared. One convenient simple system which represents the nonlinear energy exchange between a wave and zonal flow is the "minimum hydrodynamic equations" of Lorenz (1960). These 16

17 equations govern the evolution of a two-dimensional homogeneous incompressible inviscid fluid for which the vorticity field is constrained to take on the following functional representation: 5=9'lyaft~t,.4 X R+ O Aev s+ AQ81 X, (2.1) where C is the vorticity and 4 the stream function of the flow. Al, A2, and A3 are time dependent coefficients. If Lx and Ly are the disturbance wavelength in the x and y direction respectively, then k = 27/Lx and Z = 27/Ly. The type of flow considered is governed by the barotropic vorticity equation, I V2 —h — X * v'X( ) (2.2) Substituting (2.1) into (2.2) and utilizing the orthogonality properties of the trigonometric functions, we may isolate the time dependence of each spectral coefficient. The following are the "minimum equations": wr c l IA( 3 3od Al Ad RP (2.3) where c1= -[(2+)1 c2 3(2+1)1, c3= ( 2 1)/2 and a = k/Z. The dots denote differentiation with respect

18 to time. The above equations conserve the mean kinetic energy E and the mean squared vorticity V: E 432(F4A)' t (2.4) V sil la A23). (2.5) I AzI i+ Al R 3 The solutions of the above set of equations may be conveniently expressed analytically with the use of elliptic integrals or the associated elliptic functions. The latter will be used for notational convenience. The particular choice of elliptic functions depends on the initial conditions and a. Before displaying the solutions we need to establish expressions for the maximum values taken on by A1, A2, and A3. In order to do this we make use of (2.3)-(2.5). It should be noted that (2.4) and (2.5) hold for any arbitrary time, In the sequel A10, A20, and A30 denote the initial conditions. In the following, the absolute values of extrema will be denoted by A1, A and A.When A1 reaches an extremum, for example, then A1= 0 and 2 = A2. From (2.3) we have'1 =k F1 1 rom two possibilities (the steady state solution A2= A3= 0 emerges as a special case): 2 *2 1) A2= 0 A3= 0 and A3= A3, thus wtbe 0+2OfA;0e A30) - A` i-2. E = 4,I.< Aio I,- +

19 from which to + 4 Ao i>'" (2.6) As A30 + A1 0 > 2) A3=O 0 and A2 A22, thus X( + A'- zt0) 4(A*'%+ A.) from which ItL (w El~~~~~~ E(2.7) AjoD+ > tip E The restrictions on a in the above and succeeding expressions ensure that the maximum value of the squared amplitude is computed and not the minimum. The above expressions were derived by seeking out extrema and not necessarily points of maximum. As an illustration, consider, 2 the two relations given for A,, and assume for the moment that a>l. The one given in (2.6) is obviously a maximum, whereas the one appearing in (2.7) would, in this instance, be a minimum. However, if a<l, then it is possible that the latter would be the maximum value taken on by A2. The apparent conflict is resolved through examina

20 tion of the ratio a2V/E or V/E as is necessary in some cases. For the case of A2 reaching an extremum, then A2= 0 and we have the two possibilities (again, A1= A3= 0 is a special case): 1) A=0 A3= 0 and A2= A 2 thus 3 13 3 E2)o t A3 0 a A A0t + A, thus =(to a o 1 3 Xo C ( + 2 a' V." -L(A * *t A A*0++ a. o from which A ), 4- 0 Io - ~ (2.98) 2) A,= 0 =4> A1= 0 and A2= A 2, thus V* 410 + 4 A + 1 0) Ft I + A ), V= I Pk+ AtI+ I Ao AR t* At from which Al- Ao- 4<, Io J<1s E (2.9) A,, o~, 1'OIL VI,

21 A similar procedure may be followed for the case of A3 reaching an extremum at which time A3= 0. The same line of argument as above leads to the following: A~~t - W4o I EV A}' AtsJ A 1 (2.10) 1Vs = A*o + i t ol The special case a = 1 implies from (2.3) that A= A30 for all time, and from (2.4) and (2.5) that V/E > 2k2 according as JA30 > 0. If A0 > 0, then a Ia a~ V ~a a As = a +' Alo and the solutions for Al and A2 are given by trigonometric functions. For the case a = 1 and A3= 0, then Al= Alo and A2= A20 for all time. The significance of the magnitude of a becomes apparent if we investigate the behavior of a small disturbance on a time independent zonal flow, the latter being a steady state solution of (2.3) characterized by Al= A0 and A2= A3= 0. If we perturb the flow slightly so that 3he significan,: f'

A 2 A3 0s3 where a2 and a3 are small in relation to A10, then neglecting products of perturbation quantities we have from (2.3): xl busb W t + ANewt s where p, q, r, s are arbitrary constants of integration, and -r I I4 A similar set of solutions holds for a3. It is obvious that for a>l this linear analysis would indicate that small perturbations would remain bounded whereas a<l predicts unbounded growth. The former implies a stable zonal flow while the latter an unstable one. Recall that a = k/k = Ly/Lx, so that stability exists if the disturbance is elongated transverse to the basic flow. Let us now briefly note the significance of the ratio V/E for a>1. If V/E < 2k2, then from (2.4) and (2.5)

23 a j A,t C- l 4l < z A, + + I 3 A, or A < i A,. (2.11) In order that the above inequality hold in general, Al must never change sign, for in doing so, it would become zero —violating the above condition. On the other hand, if V/E > 2k2 then A3 A A (2.12) and it is possible for Al to go through zero resulting physically in a reversal of the zonal flow. This is a nonlinear phenomenon in the sense that a large disturbance as determined by (2.12) can bring about the reversal of a zonal flow, otherwise stable to small perturbations. (It turns out that V/E = 2k2 is a discontinuity with respect to the total fluctuation of Al. If (2.11) holds, then the total fluctuation of A1 is less than A*, whereas (2.12) permits twice A*.) Similar inequalities may be deduced from the relation a2V/E < 2k2. In this situation the relevant inequalities relate A2 and A2 These serve as an aid in obtaining the analytic solutions and this will become apparent later. Before proceeding further let us summarize the values

24 taken on by A2 AA2 and A2 subject to the restrictions 1, 2 A 3 on a as well as V/E or a2V/E. Table 2.1 contains such a summary. 2.2 ELEMENTARY PROPERTIES OF ELLIPTIC FUNCTIONS The following is a brief outline of the defining equations for the elliptic functions of Jacobi as well as a listing of some of the simple properties with which we shall be concerned. For further details the reader should consult Davis (1962). If we let A i A~-0 where ko is referred to as the modulus*, then we define the elliptic functions of Jacobi as follows: m (l X nee) ntt b e' (2.13) dmA () t, Al,- s oR A *Formally k2 need not be restricted to the closed interval [0,1]. However in the present application u will be found to play the role of time, and the reality of time imposes the constraint that k2 not be greater than 1. Negative values of k2 are similarly ruled out on the physical grounds that each of A, A2 and is real.

TABLE 2.1 Maximum Value of Each Component Limits on a Initial A1 A A3 Conditions V/E < 2k2 A 2+ a- A2 A2 + 2(1-a-4)-1A2 A2 + 4l(1 — )A2 1o 20 20 30 30 2 20 >1 V/E >/ 2k2 A2 + a-A2A A2 + a4A2 A2 + (l- )A2 10 20 20 10 30 2 20 2V/E 2k2 A + 2(2 A2 + (l_4A)-2 A2 + 1 ( )A 2 10 30 20 10 30 2 10 < 1 a2V/E >, 2k2 A2 + a-4A2 A2 + a4A2 A2 + 1 (1-a4)A2 10 20 20 10o 30 2 10 A2 > 0 A2 + A2 A2 + A2 A2 30 10o 20 20 10 30 A = A2 A2 0 30 10 20

26 From the above definitions it is clear that sn(O,ko) = 0 and cn(0,ko) = dn(0,ko) = 1. We also note that sn(u,ko) and cn(u,ko) may change sign, while dn(u,ko) > 0 for all admissible u and ko. This property will be useful in writing the solutions of (2.3). Using the above definitions the following additional relations may be deduced: sn(u,O) = sin u, cn(u,0) = cos u, dn(u,0)= 1; and sn(u,l) = tanh u, cn(u,l) = sech u, dn(u,l) = sech u. The derivatives of the elliptic functions with respect to u are easily evaluated:' (2.14) iaA( t From the above we observe that the derivative of sn(u,ko), cn(u,k0) or dn(u,ko) involves the product of the remaining two functions, a property common to (2.3). 2.3 SOLUTIONS OF THE MINIMUM EQUATIONS We have seen that the solutions of (2.3) will fall into three basic classes according to whether c is greater

27 than, less than, or equal to 1. For a > 1, then V/E < 2k2 implies as we have already noted that Al cannot change sign. We would thus expect its time behavior to be described by the dn function. We have also shown in section 2.1 that, for a > 1, Al reaches its maximum when A2= 0. If we let t* be this time, then from the properties of the elliptic functions the appropriate set of solutions is:,A- s, Id o (Xct- tAl 4) Al a j(t-t) (2.15) where h, t* and ko are to be determined. The constant s= + 1 according as A,10 the initial condition for Al, is positive or negative respectively. Because the dn function is always positive, the introduction of sl is equivalent to the use of -dn to describe the time dependence of Al if A < 0. This is just a mathematical convenience so 10 that A* * and A3 are always taken to be positive. 1 2 3 Differentiation of (2.15) with respect to time yields: A,-shA AA a A *As A3 —; ~ hA { (2.16) _A: A Al 3 *4 Rz A*A aA.

28 Comparison of (2.16) with (2.3) implies: A A.)3, ^ (2 +r)1-1 (A,,1 A/ = o (21-)/ 7) From the above it is obvious that the sign of h is dictated by sl, and that ko must be positive. Solving for h and ko we get: (2.18) To determine t* we note that at t = 0, wAo Ar byg de) fi n io, Awe3 A* de (u i io, where by definition,

29 (2.19) and wosideOr A, t fand Aas) tcmp o. (2.20) If we consider Al, A2 and A as the components of a position vector in three-dimensional phase space, we can now compute analytically for a > 1 and V/E < 2k2 the paths traced out by the endpoint of that vector for various initial conditions* The special case V/E = 2k2 does not need a special analysis, but the nature of the solution is worth some comment. In this case, Al and A3 are related for all t by the expression: A\i- A (2.21) With the above constraint the solution of (2.3) will approach a steady state asymptotically, irrespective of the initial value of A2. One can show directly from (2.3) that the time dependence of Al (and of course A3) is given by the hyperbolic secant, and that of A2 by the hyperbolic *This particular set of solutions was obtained by Lorenz (1960). The results of that paper, however, were based on a numerical solution of (2.3). Note that the expressions given therein for kg and ~o contain misprints.

30 tangent. This solution emerges as a special case of (2.15) because the condition V/E = 2k2 implies k0 = 1. It was pointed out in the previous section that for this particular value of the modulus the elliptic functions reduce to the hyperbolic functions of the type just mentioned. For V/E > 2k2 the character of the solutions changes because from (2.12) we note that A3 is not permitted to vanish. Consequently, its time behavior must be described by the dn function. Moreover, it can also be shown (analysis leading to (2.10)) that for a > 1, A3 reaches its maximum when A = 0. If we let t* be this time, then the appropriate set of solutions is: A, * cA 1( t -- t*), a) ) A: Al r(AX(tt*),A) i, (2.22) 3 5383 3 (\t-t )8A') where new relations are to be determined for h, t* and ko. The factor s3 = + 1 again according as A30 is positive or negative. The expressions for the relevant parameters are given in Table 2.2. The method of computation for each case listed in that table proceeds along the same lines leading to (2.18)-(2.20). Just as the class of solutions considered in the foregoing has two subdivisions, so does the class of solutions

TABLE 2.2 Relevant Constants Appearing in Elliptic Function Solutions Limits on Initial h kg Conditions Ot 2- (Ot 2 A A ta' ~ a 1 A 2 2 s 1 1 2(02 lA* 2 A23 c 1 V/E < 2k 2 t 1 1 cX4(cx4- 1) A*2 ( 2 O4 A 2+~~~~~~~~~~~~~ 30A >1 L V/E > 2k2 53 X A* 4 1 A1A2 __ Ot2 + 3 2 *2 tan-'( 2 -' 2~ A*3 3 1 0 _ _ _ _ _AIA3 a~2V/E'< 2k2- 2 1~ A* 28 *2 / 1 2 2 02(l + a2 2 1 4 A 2 tanA <1 a ak 4 A*A* A0 (x2V/E > 2k2 53 A 1 A tan-'A2 1 +ca2 3 2c O2 A * 2 t 3 2 0

32 for a < 1. It is easy to show that the condition a2V/E < 2k2 implies that A2 cannot change sign and so is described by the dn function. The solutions for this subdivision are A, — A * (i(t- ts), o ) noA SX A(~ X(t- t*" 0n,) (2.23) tab Al Corn( 0(t)t*),h, where we have used the fact that if A2 has reached its maximum at t* then A; must be zero. The role of s2 is similar to that of s3 in (2.22). Again a new set of relations expressing the constants h, t* and ko in terms of the initial conditions has to be determined. The second subdivision, namely, a2V/E > 2k2, results in the following permutation of elliptic functions: A,- A n (k(t- t{), Aa ) A, iAl C#t( (2.24) 3 SI -As( - Two subdivisions can be identified for a third class of solutions, those for which a = 1. It is easy to show directly for this class that A30 = O implies A1 and A2

33 suffer no change with time. On the other hand, if A2 > 0, 30 then Al and A2 have time behavior which is governed by trigonometric functions. This turns out to be a special case of (2.24) if ko = 0, or equivalently when a = 1. This third class of solutions then can be combined with the second, subject to the provision that a = 1 and A = 0 30 imply no time change for Al and A2. Table 2.2 provides a summary of the expressions obtained for h, kg and.0- It is now possible to compute analytically the solutions of the minimum equations for any value of a and arbitrary set of initial conditions. 2.4 THE STOCHASTIC EQUATIONS Calculation of the statistical moments of an evolving ensemble via (1.16) poses no unusual difficulties. We shall be interested in the case where the initial values are uncorrelated and distributed according to a multivariate normal density, i.e., (rtjw oaja,,o X'Al It A( do Alto ) ] (2.25) where Pio and aio are the respective means and standard deviations of the Aio. With the following change of variables,

34 (1.16) may be written as (As (gl t59;(A(9t 2 )I~~ t (2.27) where AA A3 Equation (2.27) is now in a convenient form for the use of an Hermite quadrature. It can be shown (see, for example, Salzer et al., 1952) that ~II t s(rl. w3(X ) + P X) (2.28) where the weights and remainder term are given by m2* v k LN n! J i;,oo.2(2n)! and xi is the ith zero of the Hermite polynomial Hn(x). (Note that Rn = 0 and (2.28) is exact if f(x) is a polynomial of degree less than or equal to 2n - 1.) The evaluation of (2.27) is accomplished by a three-fold application of (2.28). All of the integrations to be reported here use a value of n = 17. The accuracy of the results is typically greater than five significant figures. Of interest will be the comparison of the results obtained using the procedure outlined above and those

35 obtained by solving the respective set of stochastic dynamic equations. Using (1.8) and (1.10) we may write the stochastic equations corresponding to (2.3) as follows: 0,4 c,(UIhj 43t< 3) (2.29) UMt': c,(q,,q3i 6') (2.30) (2.31).l:2C,(8^613+# sI^) + tC|^ 6(2.32), +:Ac,(, g,,is3,,) + U2t,3 (2.33) 6iSe2Cj(MNItiP^a3)i ic3t123 (2.34) (2.35) fMtI CI (I A t A t 2 C,a ). t t I ) + 2 5x,a 3(', C r(l +G T 3A) (,q, io ) C, + C)C,?$CT t,,L (2.36) "i3,L(O.63,044),33+43 r, el(z,,), ~,~,)3,,,, (2 37) The third moment terms have been included in the above set for the sake of completeness. However, all integrations to be reported in the next section will be made by simply neglecting these terms. 2.5 RESULTS This section describes the results obtained for the evolution of an ensemble of initial states where each member of the ensemble evolves according to (2.3). In particular we shall compare the moment solutions as computed

36 from the stochastic dynamic equations, in which third moments have been neglected, with the solutions via the probability integral (2.27). The latter will occasionally be referred to as the Liouville solution. We shall choose a = 2 and the following initial conditions: P10=.12(3 hr)-1; P20=.24(3 hr)-l; p30= 0. (3 hr)-; aij = 0.(3 hr)-2, i ji; aii = 10-4(3 hr)Y2. Notice that we have taken the components to be completely uncorrelated initially, with standard deviations of 10-2(3 hr)-l. With a scaling of k = (2T/5000)km-l, the above values give a maximum expected value (ensemble mean) of the zonal (westeast) wind of about 64 km/hr with a standard deviation of 5 km/hr. Making use of an earlier result, we note the choice for a implies that the zonal flow is stable to small perturbations. The above set of initial conditions indicates that the initial disturbance is rather large in relation to the zonal flow. In fact, the energy associated initially with the ensemble means is partitioned equally between the zonal flow and the single wave. Nevertheless, this stability characteristic is reflected in the complete nonlinear solution in which the fluctuation of the zonal flow is slight. Discussion will, therefore, be centered on the time variation of P2 and P3' Figs. 2.1 to 2.3 present the Liouville solutions for the ensemble means and standard deviations, as well as the

37 results from the stochastic dynamic equations. The solutions of the latter were obtained by using a fourth order Runge-Kutta scheme with a time step of 1 hr. The agreement between the two calculations for the means is quite good through about 12 days. The reader will recall that the basic deterministic model (aii = 0) is periodic. The solutions for p2 and p3, however, show a reduction in the magnitude of successive extrema. Mathematically this makes sense because the initial ensemble is in that region of phase space where the evolution of deterministic trajectories implies that A2 and A3 oscillate about zero. Adjacent trajectories, however, have differing periods of motion resulting in a phase decoupling among them. It is easy to see that this will lead to a reduction in the maximum variability of the ensemble means about zero. Physically one can think of Fig. 2.1 as displaying the overall reduction in certain energy (see, for example, Fleming, 1971a) associated with the means. This energy is still present in the system but its specification is uncertain. This is reflected in Figs. 2.2 and 2.3 where we observe a general increase in the standard deviation of each component from its initial value. For this simple model we gain further insight into the optimal nature of stochastic forecasts. As the averaging time is increased, the time average of A2 and A3 remains near zero. This is the climatological value for these components. Note the behavior of the stochastic

0.30 Liouville Stochastic 0.20 0.10 0~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~0 -0.10 -0.20 - Q30 0 2 4 6 8 10 12 14 16 18 20 DAYS Fig. 2.1. Comparison of Liouville solution for 12 and 13 with stochastic dynamic neglecting third moments.

39 dyanmic, or better, the Liouville solutions, whose fluctuations about the climate mean are diminished for large time. Even this simple model displays a predictability limit characteristic of models used to make weather forecasts. This result is in agreement with our experience that for large times a deterministic forecast differs randomly from the observed state, and that the best forecast, in a least mean square error sense, is simply climatology. Because ensemble members may alternately converge or diverge, we do not expect that the spread of the probability distribution, as measured by the standard deviations, should increase monotonically. This accounts for the oscillatory nature of the curves in Figs. 2.2 and 2.3. It is quite apparent from the latter two figures that the stochastic dynamic equations with the neglect of third and higher moments do not give acceptable solutions beyond five days for second moments. This is less than one half the time for which there is congruency (agreement with the true solution) in the means. In order to obtain forecast variance information for longer times by way of moment equations, it is clear that there must be some treatment of the third moments. This has been discussed by Fleming (1971a) wherein a quasi-normal eddy damped closure was found to give best results. In that study a Monte Carlo calculation was taken as the reference solution. It should be pointed out that the calculations presented in this

0.25 Liouville Stochastic 0.20 I''~~0.15 0.10 0.05 0 2 4 6 8 10 12 14 16 18 20 DAYS Fig. 2.2. Comparison of Liouville solution for a 22 with stochastic dynamic neglecting third moments.

0.20 Liouville Stochastic 0.,15 -. 0.05 I L \!J' I 9l'JJ 0 2 4 6 8 10 12 14 16 18 20 DAYS Fig. 2.3. Comparison of Liouville solution for ~3a with stochastic dynamic neglecting third moments.

42 section were also carried out by Epstein (1969b), but for just the first six days. Therein a Monte Carlo solution, which itself is subject to sampling error, was taken as the standard. For example, the errors in the means are inversely proportional to the square root of the sample size. Because the integral solution permits a ready evaluation of any order moment, it would be of interest to examine the behavior of representative third moment quantities. The only term neglected in each of the stochastic equations for the variances is the triple moment T123. With an appropriate normalization Fig. 2.4 gives the time behavior of this quantity. Initially zero by assumption, T 123 remains relatively small for about 2~ days for this particular set of initial conditions. The most striking feature is the rather intense oscillations beyond the third day. As the ensemble becomes more diffuse the importance of T 123 in the moment equations would be diminished, but by this time the fidelity of the stochastic dynamic solutions has been lost. Let us now consider a model possessing more degrees of freedom. What are the implications that might be inferred from the foregoing simple system? The answer to this question must necessarily be qualitative, but I believe something can be said. It is safe to say that the

12 1.0 0.8 0.6 bt 0.4 bF 0.2 b= 0.0 -0.2 — 0.4 -0.6 -0.8 -I.0 -1.2 0 2 4 6 8 10 12 14 16 18 20 DAYS Fig. 2.4. Normalized triple moment T123'

44 discard of third moments will be detrimental to stochastic forecasts for "large" times. Also, in models containing many degrees of freedom, the temporal limit beyond which the first and second moments lose congruency probably will vary with each space and/or time scale. Fig. 2.1 shows that the stochastic solutions preserve the congruency of the means of one interacting wave triad for about two system oscillations. Second moment information becomes unreliable within one half that time period. In the author's opinion there is some reason to suspect that one might expect some improvement in second moment statistics in a multidimensional model. Consider the following two situations: 1) an interacting wave triad taken by itself as in the case of the three-component model, and 2) the same triad embedded in the space of many interacting triads of a multidimensional model. In the first situation the certain-uncertain energy exchanges for the simple system are severely limited. After the energy flow has become established, and before the ensemble has substantially diffused throughout phase space, one would expect rather large triple correlations to develop, reflecting the interdependence of the triad members. In the second situation, a given triad member is also a member of many other interacting triads. Qualitatively it is reasonable to expect that in this setting, triple correlations would not gain the same

45 prominence as in a simple three-component model. This being the case, one would expect congruency in the forecasts of second moments for a longer time interval. It should be stressed that the foregoing argument is speculative and a definitive answer would require the solutions of models with increased degrees of freedom. This will not be done here because the forecast model to be described in the next chapter will be used to make 24 hr forecasts. Evidence presented in the preceding discussion and in earlier studies (cf. Fleming, 1971a) would suggest that the third moment discard approximation is quite satisfactory for short-range forecasting. A system with four degrees of freedom has been used by Fleming (1973) to study the adjustment problem in a divergent barotropic model. This particular set of equations gave rise to solutions which were highly non-Gaussian in nature, and the retention of third moments was necessary to obtain satisfactory results for the means beyond even a few days. One particular case studied showed how one of the coefficients defining the height field could become highly positively skewed (Fig. 6 in that paper). In that case, third moments were important from the outset. On the other hand, Fleming pointed out the rather pathological nature of the initial conditions chosen for the previous example, and speculated that the encounter of such an extreme case in a model with many more degrees of freedom

1.6 1.4 1.2 1.0 0.8 0.6 0.4 N to 0.2 b 0 tC -0.2 -0.4 -0.6 -0.8 -I.0O -1.2 -1.4 0 2 4 6 8 10 12 14 16 18 20 DAYS Fig. 2.5. Normalized triple moment T333.

47 would be unlikely. Finally, we present the evolution of a third moment quantity which measures the skewness of the distribution. Fig. 2.5 is a plot of T333 appropriately normalized. Its behavior is also typical of the fluctuations exhibited by Till and T222' Unlike the rather rapid fluctuations which are present in T123 (cf. Fig. 2.4), the oscillations in T333 are approximately one half as frequent. The results from this model would indicate that the distribution does not become systematically skewed.

CHAPTER III FORMULATION OF STOCHASTIC FORECAST MODEL 3.1 MOTIVATION FOR STOCHASTIC DYNAMIC FORECASTING As pointed out in Chapter I, the goal of this study is to predict the evolution in time of flow patterns derived from observations of the 500 mb height field. Traditionally this has proceeded in the following way. Observations made at meteorological stations are interpolated to some ordered grid system. The analysis procedures vary at present from one weather service to the next. Some examples would be: 1) polynomial interpolation (Cressman, 1959), 2) the method of optimum interpolation in which the climatological statistics of the observations are used (Gandin, 1963), or 3) a variation of the previous method in which rational use is also made of the climatological statistics of the fields of apparent forecast error (Eddy, 1964, 1967; Kruger, 1969; Rutherford,1972). Each method of "objective" analysis will generally lead to a different initialization and it is not possible to say which analysis is best on a given day. After the analysis has been made, a set of meteorological equations is used to make a forecast. Obviously, differing initial conditions will subsequently lead to differences among the forecasts for any given prediction model. The variability of the forecasts computable from these different initial states will first of all depend 48

49 upon how much spread there is in the initial states. Secondly, the spread associated with the forecasts will be related in some way to the sensitivity of the model to changes in the initial conditions. This sensitivity can be quite significant in those regions of rapid development or breakdown of a large-scale circulation pattern for example, whereas in other areas it could be relatively innocuous. This study recognizes that observational and analysis errors not only do exist but that a framework can be constructed whereby the analysis and forecast problem can be dealt with in a probabilistic or stochastic setting. In order to utilize a prediction model of the atmosphere, we must first obtain the initial conditions in a suitable form. Let the vector 8o constitute such a convenient set of N quantities. The latter may be values of various meteorological parameters specified either at grid points or in terms of the coefficients of some function set. Because we are always limited to an approximate description of the atmosphere, N must be finite but can be quite large. It will be convenient in the sequel to regard these N quantities as the coordinates of a point in a corresponding N-dimensional phase space. The notion of phase space has proved quite useful in statistical mechanics for some time (Gibbs, 1902), and has also been used by a few authors in meteorological contexts (e.g. Lorenz, 1963, 1965; Epstein, 1969b; Gleeson, 1970). A

50 weather forecast then is the prediction of the trajectory through phase space of the end point of the vector X, whose value at some instant t = t is Lo, Observational errors and lack of detail in the observing network, giving rise to analysis errors, lead us to conclude that an exact determination of SO is impossible at present. Thus, it would seem appropriate from a conceptual viewpoint to think of the various single time atmospheric states as comprising an ensemble of phase points. As used here the term ensemble may be defined as an infinite number of states, or ensemble members, with the density of these states in phase space given by some assigned probability distribution. In practical terms each ensemble member could be thought of as an analysis produced by a different meteorologist. There are many such analyses (ensemble members) which are possible and still consistent with the available data. We are thus led perforce to the consideration of an initial ensemble because of the uncertain nature of the initial conditions. The forecast problem could then be envisioned as the study of the evolution of this ensemble through phase space, or more specifically, the prediction of the probability distribution characterizing the ensemble. For models possessing atmospheric realism this appears to be a much too ambitious undertaking with present day computers. Even if this task were possible, it would

51 probably be wasteful as the entire probability distribution is usually not of interest. Instead, a description of the lowest order moments will generally suffice. With the initial conditions more appropriately viewed as an ensemble of equally likely states*, we shall define the "best" state as the ensemble mean B0. This should not be confused with the true state which of course is unknown. The question that arises is whether the "best" forecast is obtained by using the conventional deterministic equations with SO as initial conditions. Epstein (1969a, 1969b) has shown that this is clearly not the case whenever a governing equation contains nonlinearities as regards statistical quantities. This has already been demonstrated in the previous chapter for a very simple model. The method of forecasting then must be modified if we are to achieve the best forecast in the sense described above. The simplicity of the model considered in Chapter II allowed us to obtain analytic solutions describing the time dependence of each ensemble member. The moment solutions sought could then be expressed in integral form and the results evaluated to any desired accuracy. From the standpoint of computing economy, this type of calculation is not feasible for a physical model which would be required to produce realistic forecasts of the medium to large-scale *Each member of the ensemble is equally likely, but the density of states in a particular region of phase space is given by the probability distribution function.

52 atmospheric motions. However, we have already seen (cf. section 1.2) that it is not possible to derive, in closed form, a set of equations governing just those moments of interest. A differential equation describing the time dependence of a moment of order m involves moments of order m + 1. For m = 1 then, we have the ensemble average of the deterministic equation, the latter assumed to govern the evolution of each ensemble member. This averaged equation gives the instantaneous velocity, B, of the ensemble mean in phase space; 8 is found to depend upon the current position of the center of the ensemble, B, as well as the size and shape of the ensemble as measured by the covariance matrix. (For the purposes of argument we can consider that these two moments are known* at some instant to). Clearly, as the ensemble evolves it will undergo deformations, and we must derive a set of equations for the time behavior of the covariance matrix, in the manner outlined in Chapter I. This is surely a plausible way in which to proceed. Having already argued for the interpretation of the initial conditions in terms of an ensemble in phase space, we are now suggesting, as a logical extension, that a forecast of the covariance matrix is of vital importance to the whole forecast procedure. We have *A procedure for estimating these quantities from observations is outlined later in this section.

53 made mention above that part of the information contained in the covariance matrix, namely the variances, is a measure of the size of the ensemble. This,of course, becomes a prognosis of the degree of confidence that we can place in the forecast of the ensemble mean, the best estimate of the true state. It turns out that the prediction of the uncertainty associated with 3 is the most significant contribution of stochastic dynamic theory to the forecast problem. In order to close the system, a statement must be made about the third moment statistics which have appeared in the differential equations for second moments. Results of the previous chapter and the work of Fleming (1971a) on a barotropic model indicate that the neglect of third moments will yield a satisfactory practical solution for the forecast problem on the order of a few days. Indeed, as pointed out by Epstein (loc. cit.), the success of the predictions of B by the deterministic equations, which implicitly neglect second moment quantities, would tend to indicate that an improved prediction for B is possible with a stochastic set of higher order approximation. By this it should not be inferred that the same measure of success is possible for second moment statistics. We have seen in section 2.5 that this is not always the case. However, recent investigations, referred to in the foregoing, would indicate that the neglect of third moments is a justifiable closure, and gives valid results for integrations over a few days.

54 In short, then, our observational and analysis capability force a probabilistic interpretation of the initial conditions. This fact leads to a probabilistic method for prediction. In a nutshell, this is the motivation and justification for utilizing stochastic dynamic procedures. 3.1.1 Determination of an Atmospheric State from Observations In the foregoing we saw that it is more appropriate to treat the initial conditions as an ensemble. The exploitation of the stochastic method depends upon the estimation of the first and second moments of the probability density function characterizing the ensemble. We shall describe a procedure for accomplishing this. As before, let S0 be a true, but unknown, description of the atmosphere at t = to. Because we shall be working with a forecast model in spectral form, the components of S0 will be the coefficients of a set of functions, spherical harmonics in this study. Let Q be the true value of the meteorological field at P points. Then, (3.1) A typical element xij of the P x N matrix X is the value of the jth function at the ith point. The matrix product is simply the sum of the products of each function value and corresponding coefficient —the calculation repeated for each of the P points. Although SO has been defined

55 as the true value of the N coefficients, the dimensionality of this vector is finite, and consequently X _0 differs from the true value of the meteorological field at the given points. The difference LT is termed the error due to spectral truncation in estimating the field. We shall consider the case where observations Y are made at each of the P points. The expected value of each is equal to its true value, but every observation is subject to an error, i.e., 2' A + Ev e elof., (3.2) The elements of 6,o or observational errors, are assumed to be independent random variables with identical normal density functions. The absence of correlations in the observational errors would appear to be quite a reasonable assumption, as these are essentially observer and/or instrument errors. The implication of identical normal densities is that the mean and variance of each member of so are invariant from one observing point to the next. In principle it is possible to assign different values to these statistics for each point. Choosing E(_o) = O suggests that systematic errors may be neglected. Moreover, we shall consider the case where each

56 observation has the same uncertainty, as measured by its standard deviation. In general, if it is known that the reliability of some observations is greater than others, then it is possible to reflect this in the standard errors assigned to the respective observations. This would be especially useful if different types of data were made available, each having its own degree of precision: radiosonde, satellite, constant-level balloons, etc. In this study we shall make use of one data source,namely radiosonde observations, and shall assign to each observation the same uncertainty. The statistical nature of the error LT arising from the neglect of small scales is not so straightforward. For convenience we assume that LT may be described by a multivariate normal density. The neglected scales will be on the order of 103km, which is certainly greater than the average separation between data points (cf. Fig.4.1), and consequently we would expect these errors to possess correlations. The computation of the appropriate covariance statistics is in itself a rather large problem. Once these quantities are available in a suitable form, their effects may be incorporated without any mathematical difficulty, but the numerical difficulties, notably the inversion of a complex P x p matrix with P ~ 450, would be substantial. Undoubtedly certain simplifications would be possible which would reduce this matrix to a diagonally dominant form. This however would still necessitate

57 certain assumptions of a somewhat arbitrary nature. It was felt the rather extensive computations involved would not alter the major thrust of this study, which is to demonstrate how stochastic forecasts may be made from real data and to examine some of the striking features. Therefore, it was decided to consider the elements of A to be uncorrelated. This problem is worth reexamination in further work and the qualitative effect would be to "reduce" the effective number of independent observations. For the purposes of this study then, we shall consider the elements of s to be uncorrelated and normally distributed, with E(_) = 0, and covariance matrix E(E E') = ao2I. The estimation of a is discussed later in this subsection. If X is an unbiased estimate of 0o, then E(S) = _o where E denotes the expected value. We shall obtain such an estimate by the well known "least squares" procedure. If B is an estimate of To and s the sum of the squares of the residuals, then S, (y - X _ ) X (3.3) Yl'y-2B'X' +'X'2_, where the primes denote transposes of a vector or matrix. Differentiation of the scalar s with respect to B gives the following vector equation: bs- I_ Y. &. _

58 Setting the above equal to zero for a minimum defines the least squares estimate* X, i.e., X X E-XX (3.4) The above are referred to as the normal equations and have the unique solution, ( X)- X (3.5) if X'X is of full rank (P>N). Even if _ is not unique, X' Y certainly is, and the significance of this quantity will become apparent in the next subsection. It is very easy to show that S is unbiased because x E() E or E(i) o When (3.5) exists, then the covariability among the estimates is given by, ig 44(f- 0)P A0Y] cZ ( X X)Y' |,(3.6) - X *Note that the least squares estimate P is irrespective of the distribution properties of c. Under the assumption of normality, x is also the maximum likelihood estimate,and this in part is a justification for the least squares approach (see, for example, Draper and Smith, 1966).

59 To make use of (3.6) we need to know a2. An unbiased estimate of a2 is given by, l2 S E (3.7) where s is computed by using B, and P>N. Note that this estimate for a2 embodies observational errors and errors due to the omission of smaller scales. With the restriction that P>N, we now have a method for estimating the first and second moments of the probability density function defining the ensemble of initial states. 3.1.2 Stochastic Analysis In the meteorological context, analysis is the process by which use is made of all available information to generate the best possible description of the atmosphere at some instant. Output from an analysis procedure becomes the initial conditions for a prediction model. Previous considerations have led us to construct a probabilistic framework, in which the observations infer an ensemble of phase points, and the stochastic forecast becomes the prediction of the evolution of this ensemble. In actual practice, of course, we do not deal directly with the ensemble, or even its probability function, but only the means and covariances of that density function. Consistent with the foregoing, stochastic analysis reflects our best judgment of these lowest order moments at some instant.

60 The analysis procedures to be used are those given in Epstein and Pitcher (1972) and the development outlined here parallels theirs. Before a new set of observations is available the best estimate of the true state of the atmosphere is given by the forecast of the ensemble mean of 3, while the covariance matrix measures its uncertainty. By definition, b -E() (3.8) These are referred to as the prior mean and covariance matrix for B. Our goal is to make a revised judgment about B by utilizing the information contained in the observations. This revision yields a posterior mean b and covariance matrix S. We shall make use of a result in probability theory called Bayes' theorem (see, for example, Winkler, 1972; Raiffa and Schlaifer, 1972; Epstein, 1962) that will allow us to make use of the vast amount of information in a stochastic forecast. Before we can do this, however, some assumptions must be made about the probability density of the forecast, referred to simply as the prior. In the present application, all of the moments of the prior are required to define just the first two moments of the posterior distribution. (In the following, the probability distribution for the observations is presumed

61 known since we have assumed the random vector E is normally distributed.) We are at liberty to choose any prior having the first two moments defined by (3.8), but the following three reasons justify the choice of the multivariate normal: 1) The multivariate normal distribution is "conjugate" to the likelihood (probability density) function of the data (observations), and this facilitates the mathematics involved. 2) Ericson (1969) has shown that the results for the first two moments of the posterior distribution are invariant for a rather broad class of prior distributions, indicating that the results are not too sensitive with respect to the prior chosen. 3) The multivariate normal is the proper prior to choose if only the first two moments are available (Tribus,1969). This avoids the introduction of a bias in the results. For further elaboration on these points, the reader is referred to the previously cited paper. Tests of the formulas obtained on the basis of the foregoing have been made by Epstein and Pitcher (loc.cit.) in a simulation study,and the analyses were found to give valid results*. Taking the prior to be multivariate normal leads to *This is a further justification for the closure scheme adopted, in that it works reasonably well in some tests.

62 the following for the posterior mean and covariance matrix: A A (3.9) S s + r xx) " For an outline of the proof of (3.9) see Appendix A. The prior covariance matrix S is positive definite by definition, and thus S-1 exists. The above formulas permit an analysis even though there may not be a sufficient number of observations to define 3 uniquely (i.e., P<N). The implication for asynoptic analysis is obvious. However, if P>N then S exists and we have A A (3.10) 6- S " + s_~' ) The inverses of the covariance matrices in (3.10) are direct measures of the precision associated with b and 3. In stochastic analysis, these matrices become the weighting factors which are determined by the forecast and observing network. If no stochastic forecast is available, then the prior can be considered vague, or non-informative, and S-1 taken as the zero matrix. In that case we must have A A P>N, and from (3.9), b and S reduce to the ordinary least squares solution. 3.2 EXTERNAL MODEL UNCERTAINTY So far, we have presented the need for stochastic dynamic theory by emphasizing the stochastic aspects of

63 the initialization. Quite apart from forecast error arising as a result of uncertainty in the initial state, there will be prediction error resulting from the inadequacies of any particular dynamic model used. The concept that there is a finite range of predictability for all finite models implies that errors will be encountered in all deterministic predictions. The model is simply incapable of describing certain physical processes or forecasting atmospheric motions smaller than some limited scale. The foregoing is simply a statement of the fact that our prediction models are not perfect. In other words, the trajectory of each ensemble member through phase space cannot be calculated with certainty. Disregarding this external model error, the forecast of internal error growth, as implied by the initial uncertainty, is necessarily an understatement of the actual uncertainty. To emphasize this point, consider a simple analogy. Suppose we have a particle moving in a force field according to some known dynamical law —for example, a projectile in a gravitational field. Given the appropriate initial conditions we can calculate the trajectory exactly*. Imagine that we wish to perform a similar calculation, but this time would also like to consider the buffeting action of air motions on the projectile. Although the primary force might be that of gravity, without a law *Fundamental uncertainty on the order of atomic dimensions is of no consequence here.

64 describing the effects of extraneous influences in detail, our model has lost its determinism and has acquired the character of a stochastic process. Even though the detailed dynamics are different, we may liken the projectile to a phase point travelling in phase space. In the latter case the extraneous influences are all those atmospheric processes which are not included in the dynamic model. This question has been considered previously by Epstein and Pitcher (loc.cit.) in which stochastic analysis was tested on a 28 component two-level quasi-geostrophic model. In one of the experiments discussed therein, a deterministic calculation was made and taken to be the "true" atmosphere. This, by definition, defines a perfect model. "Observations" were made by the addition of random errors to the true state, and forecasts carried out stochastically. The method of stochastic analysis as outlined in the previous section was performed every 24 hours on these imperfect observations and forecasts. After six days the expected values of even the smallest scales in the model had converged toward the true atmosphere, the actual differences being negligible as compared to a typical observational error. This result was possible because the basic model was completely deterministic and the stochastic analysis permitted the information content of the observations to be additive from day to day. The situation in the real atmosphere differs from the above simulation experiment in that there is a lack of

65 large-scale determinacy in the atmosphere. This, in effect, implies a decay rate for the information content of the observations. These model deficiencies give rise to forecasts which may have some systematic error* (see, for example, Rinne, 1970) in addition to a sensibly random component. We must concern ourselves with this aspect of uncertainty when making forecasts for at least two reasons. If for the moment we visualize the prediction problem as the extrapolation in time of the cloud of points defining the ensemble in phase space, then we would expect the size of this cloud, as measured by its covariance matrix, to increase. At the verification time the forecast ensemble and the ensemble defining the new set of observations should display a certain degree of "overlap." In other words, the forecast ensemble should be consistent with the new data, in the sense that their occurrence is admitted by the forecast with a reasonable likelihood. Because our forecast models are only approximations to the real atmosphere, this consistency between forecast and observation is only possible if we admit a mechanism to allow for this extra growth of uncertainty. Otherwise, the stochastic forecast will tend to underestimate the variances (equivalently the amount of spread in the ensemble). *For the purposes of this study, we shall assume systematic errors to be negligible.

66 This leads us to the more significant reason for consideration of dynamic model error. Stochastic analysis, which makes unprejudiced use of both the forecast and observed fields, requries that the error measure associated with each field be a realistic one. If the forecast uncertainty is unrealistically too small for each model parameter, then the new observations will not be accorded their due in the data assimilation or analysis procedure. There must be a balance between the losses due to decay of the information content of older observations, and the gains due to new ones. We shall return to this problem later and consider ways of parameterizing this effect. 3.3 EQUIVALENT BAROTROPIC MODEL Previous considerations in this chapter have focused primarily on the stochastic aspects of the forecast problem. Of equal importance, of course, is the physical model that is used as an atmospheric simulator. The efficacy of a forecast model in detailing atmospheric events is obviously central to any prediction scheme, but especially so in the stochastic method, where the uncertainty in the best estimate of the true state is a prognostic quantity. However, we have to pay the price of increased computer time for this additional information. If a deterministic model is N-dimensional, in the sense referred to in section 3.1, then the corresponding stochastic model

67 predicts N(N + 3)/2 quantities. For large values of N, it is not yet feasible to make stochastic integrations operationally, but this need not prevent us from investigating some of the aspects of the stochastic dynamic method using simpler models. It would be reasonable to assume that many of the effects present in a simple model would also exist, in varying degrees, in a more sophisticated model. It was anticipated that much useful information could be obtained about the techniques by experimentation with a model which traditionally has a reasonable degree of success, but at the same time is manageable. One of the results could conceivably be the feasibility of an extension to more realistic models. In the light of the foregoing argument, it was decided to test the method of stochastic dynamic prediction, on atmospheric data, using as a basic physical model the equivalent barotropic. This model may be derived in different ways but the net result is the incorporation of a divergence effect on the vorticity tendency. The derivations are quite straightforward, and one such procedure is outlined here for completeness. We start with the well known quasigeostrophic vorticity equation:;)t S vI~oQ tvlV~ a0|), (3.11) where, is the local vertical component of the relative vorticity, V/ the non-divergent wind, f = 2Qcose the

68 Coriolis parameter with 0 the colatitude, w = d and p the pressure. V/ may be expressed in terms of a stream function i which is evaluated geostrophically. The basic assumption to be made is that the streamlines and isotherms remain parallel in each pressure surface. The implication here, of course, is that all baroclinic effects are severely reduced. The thermal wind equation reveals that, with this constraint, the wind direction does not change with height; thereby we may write the dependent variables as follows: fWW {(3.12) where ( ) = p1 foP( )dp with po the surface pressure. Po Substitution of (3.12) into (3.11) and averaging over pressure yields J~t +q ^ V/*9~ +ty.Q s___(3.13) noting that A(p)=l and w(O)=O. Expand w, the vertical velocity, on a constant pressure surface: W 1. + v/. vq + ~ Assuming that V/ is geostrophic, and using the hydrostatic condition, we obtain at p,

where po is the surface density.and g the acceleration due to gravity. For an inviscid atmosphere moving over flat terrain wo= 0. Frictional effects will be neglected here as they have very little bearing on forecasts of the order of a day, but the influence of terrain height on the surface vertical velocity will be considered. Letting h be the smoothed elevation, then wo= V/o0 Vh and the lower boundary condition becomes where c0= g zo is the geostrophic stream function. The fo above may be substituted into (3.13) giving smog + EA 21 t Oh V 1 t t V" O - It "6 Ab WO (3.14);t PToaRTo. TO RT where Ao= A(po), To the surface temperature and R the gas constant for air. Choose a level p, (the equivalent barotropic level*) such that A(p*) = A2. Denoting the values at this level by a star subscript, we obtain from (3.14) t + Q1'v4Y I A RT *. ) RT, Wt*-1 a,k (3.15) *In earlier papers (e.g., Charney and Eliassen, 1949), the equivalent barotropic level Pe was defined as the level at which A(pe )=l, whereas A(p,)= 1.25. In general, this implies a greater altitude for p, than Pe'

70 The equivalent barotropic level is near 500 mb (Thompson, 1961), and limitations on data lead us to apply (3.15) at that level. Dropping subscripts, we have at the equivalent barotropic level: at( - bn)= - w. v( + C ) (3.16) where b ~-o Ao t, c $C. As. RTo RTr Because of the various assumptions which have been made in arriving at this equation, investigators, notably Cressman (1958), have found an improvement in forecasts if the parameter b is determined empirically. We shall do likewise and the following numerical values will be used: - 13 -2 -9 -1 — t b =7.8 x 10 m, c = 2.4 x 10 m sec Despite its physical simplicity (3.16) has been used for several years with remarkable success in forecasting the 500 mb height contours for 24 hours, and with moderate skill, out to about three days. We shall work with the spectral form of (3.16) and this transformation is given in the next section. 3.4 SPECTRAL FORMULATION In terms of ~, V/= k x Va and the velocity components may be written,

71 where u and v are the eastward and northward wind speeds, and X the longitude. (Recall that 0 is colatitude.) The stream function and terrain height may be expanded in terms of a convenient set of basis functions. Choosing spherical harmonics as such a set we may write: L (@kt E t)Y*0 (3.17a) ~h(e, X)-E He Ye () X), (3.17b) where the set of 4m(t) is time dependent and the Hm are n n constants. Here a is the radius of the earth, n the set of all non-negative integers, while the Ym(X,e) are non-zero for -n < m < n, wherein the m take on integral values. By definition, where the normalized associated Legendre functions, Pm(p), are given by Rodrigues' formula, P'~hg) (ni i.) (Ii)a) AJ+Y

72 In the above j = cos 0. From (3.17a) the local vertical component of the relative vorticity is easily obtained: where cn = n(n+l). The condition of reality implies where the star denotes complex conjugation. As defined, the Ym satisfy the following orthogonality condition: n ~ n(Y) zYS~m el bL.,"pt t. (3.18) The procedure of transformation of the equations of motion into spectral form is widely known and only the result will be given here. Several accounts are available (Silberman, 1954; Platzman, 1960; Baer, 1964; Merilees, 1968; Eliasen et al., 1970). In short, the transformation is carried out by substituting (3.17a) and (3.17b) into (3.16) and multiplying the resulting equation by the integral operator, lW into )( Yr SAA M3 a, for a fixed m and n. The latter operation isolates the following differential equation for mn: /.s (3.19)

73 where, L(C,,. + L% -AS aOPn( *Pt l) In (3.19) the summation extends over all q, s, p and r subject to the restriction noted. The procedure for the calculation of the interaction matrix Lpmr is given in qns Appendix B. 3.5 THE QUESTION OF MODEL UNCERTAINTY External model error, hereafter referred to as model uncertainty, has been discussed in general terms in section 3.2. In stochastic dynamic forecasting this problem is a perplexing one, but one which the investigator must face if the error information is to be useable. In terms of short-range forecasting, the situation for the future is not without prospect however. Improvement in the basic dynamic model would lead to a better forecast for the ensemble means, and the forecast uncertainty would then be due in an increasing measure to uncertainty in the initial conditions. Growth of error as a result of model uncertainty would play an increasingly minor role. However, as noted in section 3.2, this problem will always be present if the predictability is finite. For the present study, however, we must explicitly consider the effect of model error, and the remainder of this section

74 outlines a procedure whereby this effect may be incorporated. 3.5.1 Parameterization of Model Uncertainty The problem is how to parameterize the model uncertainty. Undoubtedly there are several ways in which one can prescribe this "generation"or "growth" of uncertainty. On the other hand, if the forecasts of the ensemble means are sensitive to the kind of parameterization adopted, then the method is suspect. The explicit procedure used here is to make corrections to the right-hand side of (3.19). At any given time the appropriate correction is unknown to us. As a result we adopt an approach often used in the treatment of similar stochastic problems (Srinivasan and Vasudevan, 1971; Soong, 1973). We simply make corrections which are random. Alternatively, it is akin to parameterizing, in one term, all those processes and scales not included in the deterministic terms. This additional term is random because it is not deterministic. We do not know a priori what this random forcing term should look like, and some possibilities are listed below: C 2)) Ct) inu' (3.20)

75 where a, t (2i, a31, a(4) and a(5) are random functions of nm nm time. It is not possible to support a particular choice of random forcing simply because the purpose is not to model a specific physical process, but rather provide a mechanism whereby the growth of model uncertainty may be incorporated in a rational way. This point is very important and the reader should view this purely mathematical device as an attempt to incorporate error growth which is external to the basic model. The last three choices given above are scale dependent. Of these, the last two symbolically display random forcing functions which would take on, for a fixed time, different values for each spectral coefficient. We have every reason to believe that this would be the best way to proceed. However, as will become evident later (cf. section 3.6), this would require an increase in the dimensionality of the system equal to the total number of model parameters. Perhaps a compromise could be reached in which certain scales would be grouped together to form a number of classes. Within each class the random variables a(i) nm would still be time dependent but not scale dependent. It might be useful to pursue this course of action in subsequent investigations. This appears to be an unnecessary complication at the initial stages of experimentation, and so the last two alternatives will not be considered.

76 The third choice, reminiscent of diffusive dissipation for a positive a (3) is approximately equivalent to the addition of a(3)V2t to (3.16). The wave number truncation in this study (O < m< 10, 1 < n-m < 9: odd harmonics) raises serious doubts about the presence of such a mechanism. This procedure would seem less objectionable and perhaps more appropriate on the grounds of two-dimensional turbulence theory if the model encompassed a much wider spectrum of space scales. In this study we shall confine our attention to the first two choices listed. Computations will be performed using separately both parameterizations and their relative merits will be assessed. The range of values taken on by a(1) and a(2) will be given later. 3.5.2 Time Dependence of a(l'and a(2) In a later development it will be necessary to make a statement about the time dependence of a, or more significantly, concerning the time behavior of its ensemble mean and variance. (The superscript on a has been dropped for notational convenience as the subsequent discussion is applicable to both a(1)and a(2).) Anticipating for the moment a subsequent result, consider that we have an estimate of the ensemble mean and variance for a at some specific instant of time. Since a is a random process, it is certainly reasonable to expect that these will not be the best values to use for the entirety of the succeeding

77 forecast period. Some time later the information about the a-ensemble will be obsolete, and our best estimate of a can be no better than its long term temporal mean, denoted by a, where I T (3.21) and T is on the order of a seasonal length. We shall postulate the following rather well known* first order stochastic process for a: t4~bt + ( Z):Z(*) (3.22) where T is a time constant and Z(t) a realization of white noise with mean zero and autocovariance function y (u) = ZZ a26 (u) Here u is the lag and 6 (u) the Dirac delta function (Jenkins and Watts, 1968). zz (u) is formally defined as follows: GO) -OWWOM I trk)t't) it. (3.23) (White noise, by definition, possesses a constant power spectrum, which leads mathematically to infinite variance for Z(t). The reader will note that the power spectrum of Z(t), given by the Fourier transform of zz (u), is constant. For a non-zero T the output a remains bounded.) *For an arbitrary random function Z(t), (3.22) is often referred to as the Langevin equation. Uhlenbeck and Ornstein (1930) review the application of this equation to the explanation of Brownian motion.

78 The assumed process above implies that, whatever its initial value, a will tend toward its long term mean a in a time interval of order T, and subsequently fluctuate about that value, but will do so in an unpredictable fashion. The formal solution to (3.22) can readily be written as 5tAu + N -.) lvt/*t)l (3.24.) where the impulse response or weighting function h(~) which appears in the particular integral has the simple form: Also, by definition Also, by definition, ao= a(O). It remains to deduce from (3.24) the time behavior of a, the ensemble mean of a, as well as the associated variance. Define E as the expectation or averaging operator over the ensemble. Applying E to (3.24) we get ^ —Eta et S;tU.- r a) -tlj t) (3.25) where by definition E[Z(t)] = 0 and ao= E [ao] To be explicit, a0 is the ensemble mean of a at t=O. The calculation of the variance of a is likewise fairly straightforward. By the definition of variance we have,

79 V" E: - - Making use of (3.24) and assuming ergodicity* for Z(t) we get the following: Z<Ek~ltno. c-,,),l,(,.~),l,(-t,., %z,,c- ltt), V 1A, a E miux &aivv a(tv)(t)& + (Io0-2/ =2[I- 4-4< 2+/t) t.i-&t/t) (3.26) where oa is the variance of a at t= O. Note that in arriv0 ing at (3.26) we have made use of the fact that ao and Z(t) are independent. Equation (3.26) indicates that for t > 0 a3/2T is an asymptote. We may interpret this asymptotic value as the long-term variance of a, i.e., *The ensemble average aensemble average and taverage as defined by (3.23) are equivalent.

80 Relations (3.25) and (3.26) provide us with analytic expressions for the time behavior of the ensemble mean and variance of a. These will be used in the subsequent development of the stochastic dynamic equations. The assignment of appropriate values for the statistical quantities appearing in the foregoing will be discussed in sections 4.2 and 4.3. 3.6 STOCHASTIC DYNAMIC SYSTEM In this section we shall derive the stochastic dynamic equations corresponding to the deterministic system given by (3.19). The latter must now be modified owing to the considerations of section 3.5. The following is the modified form of (3.19) for the random forcing of the first kind: W^% Axl t(yjrs i-C )LZ * Fx<] (3.27) where and the dot denotes differentiation with respect to time. The superscript on a has been dropped for convenience. Before setting down the moment equations, a summary of the notation to be used is in order. From the stochastic viewpoint, the pm and of course a are interpreted as random functions of time. This totality of time dependent random variables forms a multidimensional phase space on

81 which is defined a probability density function, the latter also being time dependent. At a given time, atmospheric data enable us to define, at least conceptually, an infinity of phase points, the density of which is given by this probability function. All of the information that we would require is contained in the probability function; but we have already indicated that its prediction, via the Liouville equation, although a relatively straightforward computation in principle, is prohibitive with present day computers. Having noted that the moments of a distribution which are of ultimate interest are usually the first two, we shall derive from (3.27) the appropriate equations governing the time behavior of the ensemble means and covariances. Recall that E was defined as the expectation operator over the whole ensemble. We shall denote the ensemble average m m of mn by Pn' i.e., (3.28) m u The covariance of mn with Uv' centered about the ensemble averages, is defined: " "- E(f I AI (3.29)

82 Note that amu is a complex quantity in the present formunv lation. Care must be exercised in computing the standard deviations of the real and imaginary parts of i m for example. Appendix C presents some of the relationships between real and complex quantities. Although we shall not compute any third moment statistics, they will of course appear in the development of the stochastic equations. Consistent with (3.29) we shall define the centered third moment as follows: V" t )(I; Hv) Y) or by expansion, VtA Lf,^ t r P fu WmQ (p3.30) The following definitions and identities will also be useful: 1) The ensemble average of Fmn (3.31) 2) The covariance between Fmn and u: mn v v r( IF XC) = EII( FO- Iv )]*, Dn Cor( r.(,'y)or(

83 3) The covariance between a and Fmn: cov(ask, (< -n (C) f) lom D I n~i (3.33) 4) The triple moment quantity: tomr3 ( Fjft ) I "*) >'I ctT ( (4, p ) (3 34) By taking the ensemble average of (3.27) we may write immediately, using (3.29) and (3.32),the equation governing the time evolution of pm, i.e.: As expecte nonlinrits y + C Ns th hn equa 3'c ] (3.35) As expected, the nonlinearity of the prediction equation gives rise to the second moment statistics oPr and qs cov(aMn). Time differentiation of (3.29) yields a set of equations for the behavior of the covariances between each pair of spectral coefficients in terms of known quantities. Performing the indicated operation and utilizing (3.27) and (3.35), the following sequence emerges:

84,V,,,,, +' ",'rV",) -n ct + Dv Cf C E(tS t-) et Ct E ('t ) L.+E(t', Pt'i"=.w ~-' rA c. PIWIV DLIt Cq f,My K% from which we obtain after some cancellation: rt WI% v ++ r(. H,,t 4)+D, nr +b,, CIC,(8{"SL +'I-s, +T )l c Hs tL,

85 v+ T;V,+ to ) 3<D;D c,(,<rv. (3.36) Following exactly the same line of argument leading to (3.36), we may derive an equation governing the time dependence of cov(a, m). Starting with the time derivative of the definition of cov(a,4m), utilizing (3.22) in addition to (3.27) and (3.35), we may write: =D, [ c,{e E ( J +,) + ~C E. (i, ) + E(i. F,,,, V _P,r, t pr IcW(.K~1%~s + 0% ~r,, + c,{ $L, +4st TC~ O -

86 and upon simplification, tCRs cavk(os 7)}Lt s f ter(i, ) tax OX.0 (1 n &-t (oy). (3.37) To accomplish a formal integration of (3.35) to (3.37) a statement must be made about third moment quantities, in effect a closure approximation as referred to in section 3.1. At that point it was indicated that barotropic integrations would suggest a discard of third moments as a good approximation for short-term forecasts. This approximation is reaffirmed here, and all integrations to be reported have been made subject to this closure scheme. In order to carry out an actual integration of the aforementioned set of equations, calculation of the first and second moment statistics of the initial ensemble is necessary. The algorithms for accomplishing this have already been given in section 3.1. The details of the time integration are outlined in Appendix D.

CHAPTER IV RESULTS OF FORECAST MODEL Several integrations of the stochastic dynamic equations have been made, and the purpose of this chapter is to summarize the essential features of some of these numerical experiments. Whereas a stochastic dynamic forecast contains a vast amount of information, in addition to a prediction of the expected state of the atmosphere, the experimenter must necessarily describe those aspects of the results which appear to be the most relevant to the particular problem considered. The reporting of stochastic dynamic calculations becomes especially troublesome, mainly because of the large number of covariance statistics which are present and collectively influence each analysis and forecast. A detailed discussion of these quantities would provide for a cumbersome and rather dreary account. For this reason, we shall focus our attention on the expected state of the atmosphere and its associated uncertainty field, bearing in mind that these are dependent upon all the statistical parameters in the model. The data used in this study are the heights of the 500 mb surface as obtained by approximately 450 radiosonde stations in the Northern Hemisphere for the period December 1 to 6, 1969, and herein referred to as Day 0 through Day 5. A typical distribution of stations is plotted on a stereographic projection in Fig. 4.1. This is representative 87

i x X X' X X XX X CX X x x~x x 0 X.1~~~~~~~~~~~-I'.jr) XX X x x U) X X X XX~~X XX X XX 0 X X XX XX X "x x X XX X X X XX X X X X X xx x XXX Xx X w ~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ rJ' X( x.x Xx r XXr x X X XX XX X X X X X X X XXX CO ~~~~~~ ~~~~~ ~ ~~X X X X X X X X X X X X X X XXXX XXX X X XX — X —~~~~~~~~~~~~~~~~~ XXX X X X XXX X XX X X X X X X X X XX0 00 X ~~~~~~~~~~~~~~X X Xx X XX X XX XX - 00 X X ~ ~ ~~~~~~~~X X X XX X XX X X X X X X X 4 o XX x x~~ ~~~ ~ x( X XX XX X X XX X X X X X X xx x x x x:~~~~~~~ x~~~~~~~~~~~~ x GO L x x x X ~~~~~+x x " x. x x rcx: ) " x ~~ ~ ~ ~ ~ ~ ~~~~~ ~~~X x x X X X xi a)~ XX X X X X X X Xx x X x x x~~~~~~~~~~~~~~~~ 00 x xx:::So oo 65 ~~~~ ~~ ~ ~~~~~~~x. x,(x X X x xx x x X X X E~~~~~~~~~~~~~~~~~~~~~-4 -) ~~~~~-x x x x ~F ~~~~~~~~xxx x xe X U) \ MO X r-I4 x r r ~x,,,X' x / X X X t xx x ~x ~ ~ ~ x x~ x~ ~~~~~~~~~~~~~~~~~ r~ x r~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~rs x x x~~~~~~~~~~~~~~~~~~~~~~~~~~~'x x~~~~~X x ~~~~~~~~~~~~~~~~~~~ ~J x ~~~~~ ~-r- ~:

89 of the irregular distribution of stations north of the 20 N latitude circle, and is presented simply to emphasize that we are dealing with actual station data, notwithstanding those eliminations for obvious transmissions errors. Computation of the 105 spectral coefficients and variance-covariance statistics every 24 hours is relatively straightforward in the manner given by (3.5) and (3.6). In the sequel, reference will not be made to the observations themselves, but to their least squares estimates. These are the best estimates of the true state based on the data available. Less than 10% of the data is located below 20 N. In order to stabilize the least squares fit in these equatorial regions, it was found necessary to insert artificial data there. Approximately 20 points were added and assigned climatological values appropriate for the month of December. For our purposes this imposes no real limitation because the model is ineffectual in low latitudes in any case. In addition, our primary interest in this study concerns midlatitude forecasting on a time scale less than the time required for spurious influences to be felt from the boundary. (It should also be added that these climatological data were assigned the same uncertainty as all the other observations. This uncertainty was greater than a typical day to day variation in the tropical height field, and thus the uncertainty in the estimation of the tropical height field was not made artifically optimistic by the intro

90 duction of climatological values.) Equations (3.35) to (3.37) were derived under the assumption that the external model error could be parameterized through a stochastic forcing term, the latter given by the first possibility listed in (3.20). The need for this process is presented in section 4.1, in which the results were obtained without the incorporation of this effect. The influence of the random forcing of the first kind is described in section 4.2. Its limitations are noted, and section 4.3 is a summary of the calculations made with random forcing of the second kind. 4.1 NEGLECT OF EXTERNAL ERROR GROWTH A useful feature of the least squares procedure is that knowledge of the covariance matrix S^ enables one to construct a corresponding field of standard deviation. This is a graphic illustration of the relative confidence that can be placed in the field recovered from 3. Figs. 4.2 to 4.5 give the standard errors of estimate of the observed height fields on Day 0 through Day 3, respectively. It is clear that the general patterns are similar from day to day. Largest uncertainties appear over oceanic regions reflecting the paucity of observations there. The details, however, do change from day to day. The structure of the uncertainty field depends upon the number and distribution of stations in the observing network, and also upon the value of each function at the individual stations.

.127 Fig. 4.2. standard deviation of least~ squares estimation of 500 b height Fig. 43. Sae as Fig. 4.2 but for Day 1... field, Day 0. Units are in rreters. ~~~~~~~~~~ ~~ ~ ~~~~~~~~~~~~~, L. o..r ),~~~~~~~~~~~~~~~~~~~ Fig. 4.4. Same as Fig. 4.2 but for Day 2. Fig. 4.5. Sane as Fig. 4.2 but for Day 3. ~~~~~~~~.-t,..'I''I\I Fig. 4.2. Standard deviation of least squares estimation of 500 mnb height Fig. 4.3. Same as Fig. 4.2 but for Day 1. field, Day 0. Units are in mesters. I I ~ i..:., ~''i Fig.~~ 4S.4. Sam a s Fi.4.2 butforDy2 tt )') I~~~~~~l~~Fig. 4.5n~. Sam~ as Fi. 4;.2~~ bu fo_~rt~ Day 3.

92 Missing observations and the addition of ones not present 24 hr previous give rise to a restructuring of the uncertainty field, as is quite evident in the Pacific on Days 1 and 3. The procedure for the calculation of these fields of standard deviation is quite straightforward. If xi is the vector of ordered function values at an ith arbitrary point, then the best estimate of the height field at that point is given by, (4.1) The variance of this estimate is easily obtained in matrix notation, i.e., -V E[ )(_x E E X1'* i t^ X;(4.2) where S^ = -2 (XX)-1 from (3.6). Computation of an estimate for a via (3.7), for each of the six days, revealed that a value of 40 meters was representative. Recall that a includes observational errors as well as errors resulting from the neglect of smaller scales. Its value depends upon the number of functions chosen and the density of stations available. We have already discussed in subsection 3.1.1 possible refinements, in which different locations may be

93 assigned different values of a, and correlations between points may be included. The square root of the values obtained by repeated application of (4.2) is the contoured quantity in Figs. 4.2 to 4.5. Fig. 4.6 is the least squares estimation at Day 0 of the height field whose uncertainty field has been given in Fig. 4.2. The expected values of the 105 model parameters at Day 0 along with the variance-covariance information, for a total of 5670 numbers, form the initial conditions for a stochastic prediction. Fig. 4.7 is the expected state from the 24 hr stochastic forecast valid at Day 1. (All forecasts to be reported here are made for a 24 hr period.) However, this forecast was made without incorporating the parameterization of model uncertainty. What this means is that all terms involving the statistics of a in (3.35) and (3.36) have been ignored, and this automatically excludes (3.37) from the set of equations integrated. We intend to confirm that the absence of these terms is mute testimony to the need for a treatment of error growth attributable to the deficiencies of the spectrally truncated equivalent barotropic model. Comparison of Figs. 4.7 and 4.8 shows that there is reasonable agreement between the forecast expected state and that as derived from the observations on Day 1. Notice that the intensification of the pressure ridge over British Columbia and Alaska is predicted with some skill, and the southward push of the cut-off low over eastern Canada is

Fig. 4.6. Least squares estimation, Day 0. Fig. 4.7. Forecast expected state of height Units are in meters field, Day 1. No random forcing included. Fig. 4.8. Least squares estimation, Day 1. Fig. 4.9. St~chastic analysis of expected values, Day 1.

95 detected. We are now in a position to carry out a stochastic analysis according to the procedure outlined in subsection 3.1.2. In this case we wish to meld the spectral coefficients characterizing the forecast field given in Fig. 4.7, with the corresponding coefficients derived from a least squares fit of the data and used to plot Fig. 4.8. The posterior (analyzed) values of the ensemble means are used to construct a height field which becomes an updated or revised best estimate of the expected state of the atmosphere. This analyzed field is given in Fig. 4.9. Under the present circumstances it is expected that Fig. 4.9 will bear a fairly close resemblance to both Figs. 4.7 and 4.8. This is indeed the case although the analysis does differ in detail from the forecast and "observed" fields. Further examination of Fig. 4.9 reveals that on this particular day the analysis tends to "favour" the least squares estimation. Examples of this are the analyzed positions of the trough in the neighborhood of the Aleutian Islands, and the one in the coastal regions of northern Europe. In these, as well as other areas, there is greater similarity between Figs. 4.8 and 4.9 than between Figs. 4.7 and 4.9. This of course does not necessarily mean that Fig. 4.8 is the correct picture, but only that this estimate has a smaller uncertainty than the expected state as given by the forecast. There are many factors, such as correlations between

96 model parameters, which bear on an individual analysis, but which cannot be presented conveniently in one diagram. One important input to the analysis algorithm, however, is the quantitative degree of confidence placed in both the forecast expected state, and that as derived solely from the observations. We have already seen that fields of standard deviation are easily displayed and readily interpreted, and it should be obvious, intuitively at least, what the role of these relative measures of uncertainty would be in stochastic analysis. This is given further consideration below. The reader should be aware, however, that the fields of standard deviation do not give the total picture as regards the uncertainty information. For example, large standard deviations in a particular region do not necessarily mean that the shape of the contours (e.g., a ridge versus a trough) is of no certainty. There may be high spacial correlations that would indicate the presence of a particular pattern. The actual point values may be rather uncertain however. Figs. 4.10 to 4.13 are the uncertainty fields corresponding to the meteorological fields given by Figs. 4.6 to 4.9 respectively. Fig. 4.10 is the standard deviation of the least squares estimation at Day 0 and has been presented earlier as Fig. 4.2. (Figures will be repeated occasionally to facilitate comparisons.) The uncertainty field of the forecast expected state is given by Fig. 4.11. For the most part, the latter shows an increase over the

97 former in the level of uncertainty. We would anticipate this result, and it is simply a reflection of the internal error growth as determined by the model dynamics. It is interesting to note that an appropriate flow regime, such as the southwesterly current impinging on the British Columbia coast, can cause a further increase in the standard deviations locally. This could be thought of as "advection" of uncertainty. As the center of maximum uncertainty moves eastward across the central Pacific, the stochastic forecast predicts a reduction in the highest value, from 85.6 m to 75.4 m. This result is in part a manifestation of the detailed influence of the correlations between the spectral amplitudes in the model. We have already seen in section 2.5, albeit in the case of a very simple model, that the uncertainty of each model parameter is not necessarily a monotonically increasing function of time throughout the forecast period. Given in Fig. 4.12 is the standard deviation field of the least squares estimation, valid for Day 1. With the exception of a considerable portion of the Pacific, in which fewer observations were available on Day 1 than Day 0, the standard error of the forecast (Fig. 4.11) is everywhere greater than the corresponding uncertainty associated with the observed expected state. Qualitatively, this would indicate that the forecast expected state would be given less weight than the least squares estimation in

i,' 1 Fig. 4.10. Standard deviation of least Fig. 4.11. Standard deviation of forecast squares estimation, Day 0. expected state, Day 1. ri~~~~~~~~~~~~~~% _____ ~~~~~~~~~~~~~~00 7.a I~~~~~~~~~~~~~~~~~~~~~~~~~~~~I7,L 3 7~. ~~~~~~~J S~~~~~~. ~~~~~~YF~~~~k\U34: t~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~is Fig. 4.13. Standard deviation of analyzed Fig. 4.12. Same as Fig. 4.10 but for Day 1. expected state, Day 1.

99 arriving at an analyzed expected state. This is indeed the case and corroborates the remarks made earlier concerning Fig. 4.9. Fig. 4.13 is the analyzed standard deviation field. The question that naturally arises is whether Fig. 4.13 is a realistic measure of the confidence that one may place in the analyzed expected state as given in Fig. 4.9. Because the true state is never known to us, a test of the realism of Fig. 4.13 by making a direct comparison of Fig. 4.9 with the true state is not possible. Clearly, however, the validity of Fig. 4.13 depends upon the extent to which the uncertainty in the observed and forecast expected states is faithfully represented. In the case of the uncertainty associated with the least squares estimates, we shall take this as being a realistic estimate of the actual uncertainty. As pointed out in subsection 3.1.1, correlations in the "observational" errors have been neglected, and their inclusion will undoubtedly have some effect. Although the latter should be considered in further work, it is not anticipated that this effect would alter substantially the following line of argument. The forecast uncertainty is subject to scrutiny however. During the discussion of external model uncertainty in section 3.2, the suggestion was advanced that, at the very least, there should exist a consistency relationship between the stochastic forecast and corresponding least

100 squares estimation. There, it was indicated that the respective ensembles should display some overlap. This can be tested from knowledge of the position of the expected state in phase space together with the covariance matrix. We shall return to this later. In an earlier investigation by Epstein and Pitcher (1972), it was demonstrated that, whenever the evolution of a fluid system cannot be described by a purely deterministic model, it is necessary to account for this lack of determinacy when making stochastic analyses and forecasts. In the present context the fluid system under consideration is one half the global atmosphere, and the deterministic model the equivalent barotropic. The intention here is to confirm and validate the above conclusion. The subsequent three figures are presented to show the relevant standard deviation fields 24 hr later at Day 2. From the analysis at Day 1 a stochastic forecast is made, again without the parameterization of external model error, and the uncertainty field plotted in Fig. 4.14. By contrast, the corresponding standard deviation of the least squares estimation is given in Fig. 4.15. It would now appear that the forecast expected state possesses less uncertainty than that as derived from the observations. As one might have expected, it turned out that in this circumstance greater credence was placed in the forecast expected state in arriving at the subsequent analysis. Fig. 4.16 is the standard deviation of the analyzed expected state.

Fig. 4.14. Standard deviation of forecast expected state, Day 2. No random forcing included. 26.3 Ii`5':./-./ k. Lsa:': i'~;...''k...8 2s.5.' Fig. 4.15. Standard deviation of least Fig. 4.16. Standard deviation of analyzed squares estimation, Day 2. expected state, Day 2.

102 The reader will note an even further reduction in the level of uncertainty from that as given by the analysis on Day 1 (cf. Fig. 4.13). It should be clear that realistic values for the predicted uncertainty are crucial to the determination of meaningful uncertainty in the analysis. This is perhaps best illustrated through an examination of the time behavior of some of the spectral coefficients. Recall that pm was defined as the instantaneous m m mean of Pn The real and imaginary parts of pn are defined as follows: Hi'f He + V W. The time sequence of forecasts, least squares estimates and analyses of A6 is given in Fig. 4.17. At Day 0 the analysis is identical with the least squares estimate because the latter is the only information available. Although the forecast and observed values of A6 do not agree 7 at Day 1, they are nevertheless consistent in that they are not widely separated as indicated by their respective standard deviations. This consistency is evident at Day 2 as well. All spectral coefficients do not display the agreement exemplified by the foregoing. Fig. 4.18 gives the time evolution of A7. By Day 2 the consistency between 12 the forecast expected value and that as deduced from observations is not obvious.

103 0.0 7 1.0 I'0 -2.0 Oftwo ~ -3.0 -4.0 1 00 i 2 2 TIME (DAYS) Fig. 4.17. Time sequence of A6 The dashed vertical lines represent the least squares estimates (expected value plus and minus one standard deviation) from observations on the day in question. The solid vertical lines represent the 24 hr stochastic forecasts, immediately to the left, and stochastic analyses, immediately to the right of the least squares estimates. 3.0 - - I II I'g 2.0OI. L I i'0 ~ II"" 0.0 0 I I 2 2 TIME (DAYS) Fig. 4.18. Time evolution of A7. Notation is the same as that in Fig. 4.17.

104 Even though the true values of the spectral coefficients are unknown, a satisfactory check on the validity of the forecast standard deviations is still possible. If D denotes the vector difference between a forecast and corresponding observed state, then the ensemble average, denoted by D, is simply bf-H A (4.3) The associated covariance statistics C of the random vector D are easily shown to be: C E(q — _ -], _+r _ t (4.4) Recall that b was defined as the forecast expected state, while S and S^ are the covariance matrices of the forecast and observed states respectively. Hence the variance of an element of D is simply the sum of the corresponding forecast variance of, and 2t, the variance of the least squares estimate. Moreover, CD is a measure of the total size of the combined forecast and observed ensembles. We shall consider this to be the gross uncertainty of the true state as inferred from knowledge of the forecast and observed uncertainties. Of interest will be the actual difference between the forecast and observed means —this is given by (4.3) —versus this inferred uncertainty.

105 Because of the correlations between the various elements of D, the calculations required to establish the foregoing statistical comparison in detail would be massive. Suffice it to say that it would be necessary to extend the concept of the credible ellipse (Epstein and Fleming,1971) to a credible hyper-ellipsoid in 105dimensional phase space. For our purposes this calculation is not really necessary, and we shall settle for a less precise but more direct solution. The results of this calculation are summarized in Fig. 4.19. This is a scatter diagram of apparent forecast error —defined as the absolute value of the difference between the means of the stochastic forecast and least squares estimates —versus the standard deviations inferred from forecast and observational uncertainty. Taking the normal distribution as the standard, one would expect about half of the points to lie below 0.675a, provided the elements of D were uncorrelated. The latter of course is not true in general, so that Fig. 4.19 cannot be considered as definitive. The figure does indicate the beginning of a trend however. The compatability between the forecast and new observations is poorer on Day 2. One can arrive at this conclusion by noting that the points for Day 2 are shifted to the left, in relation to those for Day 1, without a discernible reduction in the apparent errors. The overall reduction in the values of a is entirely accounted for by the reduction in forecast

106 4.0 DAY 1I DAY 2 3.5 3.0 e=2oe 2o 2.5 0 2.0 OD 0 1.5 e= o — e=oc I.0 e=.675&b 1~ e.:.675 0.5 0 0 0.5 1.0 1.5 0 0.5 1.0 1.5 u (10-8 Sec'1 ) Fig. 4.19. Scatter diagram of apparent forecast error e vs a, the standard error as inferred from observational and forecast uncertainty; a = (a + a)1/2. The apparent error is defined as the magnitude of the difference between the least squares estimate and forecast value.

107 uncertainty. The consequence of this has already been indicated in the discussion of Fig. 4.14: the ensuing analyzed expected height field bears a closer resemblance to the forecast than to the observed expected state. If the results of forecasts for subsequent periods were plotted in a form similar to Fig. 4.19, then the trend to smaller uncertainties in the forecasts would persist. Because the new observations effectively would be ignored, the subsequent forecast would tend to "drift" away from reality. This is a manifestation of the physical simplicity and deterministic nature of the basic model, and has been discussed in section 3.2. At this stage we have indicated, consistent with an earlier investigation (Epstein and Pitcher, loc. cit.), that external error growth must be considered when applying stochastic dynamics to the real atmosphere. The rest of this chapter will focus on this problem and present the main results for the two parameterizations of external error growth considered in this study. 4.2 RANDOM FORCING OF THE FIRST KIND Extension of the previous calculations to allow for an extra error growth, which is attributable to the deficiencies of the physical model, is relatively straightforward. In principle, what we are attempting to do is permit a greater spread among the evolving ensemble members. The net result is a further smearing out of

108 the ensemble as it moves through phase space. Although the selection of a mechanism which generates additional uncertainty is somewhat arbitrary, the subsequent redistribution of this uncertainty among the model parameters is not. This of course takes place according to the specific internal model dynamics. In section 3.5 we presented examples of various stochastic processes which can serve as appropriate parameterizations. The first one listed there involved the addition of a term of the form a(l)(t)mn to each deterministic equation governing the evolution of the individual ensemble members. By assuming a particular process for a(ll, we have seen that this effect may be incorporated into the stochastic dynamic equations, with the behavior of the a(')-ensemble given analytically by (3.25) and (3.26). The covariance matrix must now be expanded to accommodate the covariances between a(1) and Wm These are taken to be zero initially, and develop according to (3.37)* as a result of a non-zero variance assigned to a(1). This brings us to the assignment of appropriate values for the statistics of atl) which appear in equations (3.25) and (3.26). In the absence of any evidence to the contrary, we shall take the long-term mean 1tl)= O. This simply says that the physical model does not have systematic errors *Recall that the superscript on atl) was dropped in (3.35) to (3.37) for notational convenience.

109 which can be accounted for through the addition of a forcing term of the form a(1) mstart of the integration period, our knowledge concerning the ensemble mean ct1) is vague. We can do no better than use the longterm mean. Hence we shall take a(1) = 1)}-= 0. Thus, for the first 24 hr forecast a(1)=0. Parenthetically it can be pointed out that the presence of the covariance statistics, cov((al),m), permits an analysis for (1i) at Day 1 n 0 and thereafter, even though we have no way of directly observing this quantity. Should the analyses for a(1) show 0 a tendency toward one sign over a long period, then our judgment about a(1) could be revised accordingly. Turning to the consideration of the parameters appearing in (3.26), we find that their specification is not so straightforward. The long-term variance of al1), identified previously as az/2T, is in doubt for we do not know a priori what an appropriate range of a(1) would be. An order of magnitude estimate is possible however, and it is found a postiori that this will suffice. This can be done in the following way. Under the assumption that the order of magnitude of 0Inl is computed correctly from (3.27) with a(1)= 0, then the value of a(l) required to make I|(Il)I one order of magnitude smaller than a typical derivative I|nl is taken as a reasonable estimate of the long term standard deviation of ay(). Although this procedure appears somewhat ad hoc, it is not too restrictive because we are not specifying what a(l) should be - it, after all, is

110 a random function —but only what a likely range of values might be. Alternatively, the foregoing could be thought of as a specification of the width of the climatological a(l)-ensemble. In particular we shall take az/2T =8.87 x 10-3sec2. The same value will be assigned to a2 for the first 24 hr period. On subsequent days an analysis for the variance of a(o) is possible, in accordance with earlier remarks, and this becomes the new value for a2 in the next prediction. The time constant of the random forcing must also be specified. We shall consider the meteorological noise to be the order of synoptic and sub-synoptic scales and consequently of still relatively low frequency. The value chosen for T was 0.5 day or 4.32 x 104sec. A few experiments were performed to test the sensitivity of the model to changes in these parameters. Hardly any discernible difference between the various calculations was noticeable in the plotted fields. Virtually no variation was observed in the spectral coefficients from one run to the next. Second moments were likewise rather insensitive to changes in the a(l)-statistics, with differences between runs of about 10% for a 100% alteration in some of the random forcing parameters. In reviewing the results, our prime interest will be the manner in which these calculations differ from the ones presented in the previous section. With this goal in mind we shall examine the stochastic forecasts and analyses

111 starting from Day 0. Fig. 4.20 presents the expected state of the stochastic forecast, valid for Day 1. This is directly comparable with a similar calculation made previously, but without the random forcing effect, and given in Fig. 4.7. There is essentially no difference whatsoever between the two fields. This simply says that for a time interval of at least 24 hr the course of the ensemble mean through phase space is not materially affected, even though the individual ensemble members are subject to random perturbations —an encouraging result. Because of the nonlinearity in the governing equations, we would not expect this situation to persist indefinitely however. On the other hand, comparison of the respective uncertainty fields reveals that the stochastic forecast which includes external error growth produces an overall increase in the field of standard deviation, as compared with the earlier forecast in which this effect was not included. Fig. 4.21 provides evidence of this. The greatest increase occurs in the vicinity of the cut-off low centered over eastern Canada. Here the uncertainty has risen about 20 m over its value in Fig. 4.11, the latter being the uncertainty field from the stochastic forecast computed with the neglect of external error growth. This uncertainty maximum is approximately coincident with the expected position of the trough. Fig. 4.22 is the standard deviation of the subsequent analyzed

Fig. 4.20. Forecast expected state of height field, Day 1. Random forcing of the first kind. Fig. 4.21. Standard deviation of forecast Fig. 4.22. Standard deviation of analyzed expected state, Day 1. Random forcing of expected state, Day 1. the first kind.

113 expected state and differs little from Fig. 4.13. Figs. 4.23 to 4.25 give the forecast expected state, its associated uncertainty, and the uncertainty of the the analyzed expected state, all valid for Day 2. Again, a maximum has developed in the forecast standard deviation field (Fig. 4.24). This maximum is near the trough referred to above. Because the basic physical mechanism in the present model is vorticity advection, it would seem plausible that those regions in which advection is known rather imprecisely would also be those areas in which the predicted meteorological field becomes uncertain the most rapidly. Under this assumption, it would appear that vorticity advection is of some importance in determining the motion of the quasi-stationary trough we have been considering. We could not account for the uncertainty maximum otherwise. In other words, it seems reasonable to expect that the extra uncertainty imparted to the system would be directed to those regions where the flow pattern is the most sensitive to fluctuations in the amplitude and/or phase of its constituent waves. Fig. 4.26 provides a check on the degree of compatability between the forecast expected values and least squares estimates on Day 2. This figure is strikingly similar to the one given in Fig. 4.19, also valid for Day 2, but for the case of no random forcing. Although the present parameterization of model

2 Fig. 4.23. Forecast expected state of height field, Day 2. Random forcing of the first kind. I......'"','':...... Fig. 4.24. Standard deviation of forecast expected state, Day 2. Random forcing of Fig 4.25 Standard deviation of analyzed the first kind. expected state, Day 2.

115 3.5 3.0 e=2o 2.5 2.0 00 o 1.5 e~o'.0' e:.675c 0 0.5 1.0 1.5 ao(108Sec' ) Fig. 4.26. Scatter diagram of apparent forecast error e vs a, the standard error as inferred from observational and forecast uncertainty; a=(oa + cg)1I2, with symbols defined in section 4.1. The apparent error is defined as the magnitude of the difference between the least squares estimate and forecast value at Day 2. Random forcing of the first kind.

116 uncertainty is behaved and does work in the right direction —enhancement of the uncertainty associated with the forecast expected state —further integrations lead to reduced standard deviations for nearly all the spectral coefficients. The reason, of course, is that the decay of the information content of the observations is not being properly modeled. So, as with the case of no external error growth during the forecast stage, the new observations are relegated to a secondary role in the subsequent updating or analysis step. Increasing the long-term variance of ac(l, in order to "feed" more uncertainty into the system and consequently raise the uncertainty in the forecast expected state, is a possible solution, but not a satisfactory one. The reason for this relates to the rather peculiar and unphysical nature of this particular parameterization. It should be evident to the reader that allowing the same range of values for al'c in each spectral equation will, in general, lead to the most rapid error growth in the largest spectral coefficients. (The forcing term after all is given by the product a(l)pm.*) These, for the most part, n represent the largest scales of motion. The contribution of the random forcing terms to the rates of change of the smaller scales will be less pronounced, and this will imply a smaller error growth in these scales. Therefore, while the decay of the information content of the observa

117 tions may be correctly modeled on the largest scales, the treatment of the smallest scales is clearly unsatisfactory. A most striking example of the foregoing is illustrated in the next two figures. Fig. 4.27 gives the time behavior of A0. Note the rather rapid growth of its standard deviation during the forecast period. By contrast, Fig. 4.28 is a similar plot for A7 which is 12 approximately one order of magnitude smaller than the former. In this case, the random forcing has produced no significant error growth. These figures provide an extreme example, but they are nevertheless indicative of the selective nature of this particular parameterization. This result, as well as additional evidence of the decreasing impact of new data, necessitates a reexamination of the forcing function. In summary, it is possible to say that the first choice for the parameterization of external error growth has been partially successful. This has meant a balance between the error growth in the model and the impact of new observations, at least for the large scales. The forecast uncertainty in the smaller scales has not been substantially affected however. Previous studies suggest that these scales would suffer a loss of predictability before the larger ones. A further improvement in the results is discussed in the next section wherein we adopt a different parameterization for external error growth.

118 -3.8 -r.... cF-4.4 o; i -4.6 T -4.0 - x I -O I -1 2 2 3 -3 TIME (DAYS) Fig. 4.27. Time evolution of A. Random forcing of the O. e TIME (DAYS) Fig. 4.2. Time evolution of AO. Random forcing of the first kind. Notation is the same as that in Fig. 4.17.

119 4.3 RANDOM FORCING OF THE SECOND KIND We turn now to a consideration of the second parameterization of external error growth employed in this study, and in particular investigate whether the short-comings of the previous formulation may be averted. The alternate forcing mechanism is listed as the second choice in (3.20) and involves the addition of a random function a(2)(t) to each spectral equation (3.19). Because the nm are complex quantities for m > 0, a(2) is n likewise complex. Let its real and imaginary parts be denoted by ar and ai respectively. The change in the random forcing necessitates certain modifications to the stochastic dynamic equations, and these changes are outlined in Appendix E. The calculations to be discussed in this section are based on the set of prognostic equations given therein. Before we can carry out an integration of these equations, a specification of the a(2)-statistics is necessary. In this regard the same guidelines were followed as set out in section 4.2 in which similar parameters were assigned to at1). Thus for the first 24 hr the ensemble mean a(2)= 0. Once this initial forecast has been made, the covariability between a(2) and lnm permits an analysis for (o2)as well as the variances of ar and ai. Numerical experimentation is necessary as an aid in the selection of appropriate long-term variances for er and ai. These statistics, along with the time constant

120 T, completely define the nature of the forcing parameters. Choosing T =0.5 day, as in the earlier formulation, leads to an undesirable result; the spectral parameters in the model become highly correlated and the covariance matrix tends to lose dimensionality or become singular as a consequence. This value of the time constant implies a large autocorrelation in the a(2)-process. Under this constraint, a(2) is essentially a red-noise (low frequency) process. The results suggest that such a low-frequency forcing of the system produces an organized statistical dependence among the spectral components of the model. Since we do not observe a similar outcome in the absence of random forcing, we must conclude that this interdependence between the model parameters is an artificial one, which is attributable to the particular nature of the parameterization chosen. The preceding difficulty may be avoided by assigning a smaller value to T. The implication of course is that contiguous values of a(2) are now less correlated than before,or in the present context the a(2)-process has been "whitened." The value chosen for T was 10 min. Experimentation led to a choice for the long-term variance of -24 -4 both ar and ai of 4 x 10 sec Comparison of Fig. 4.28 with Fig. 4.29 typifies the major difference between the present formulation and the one upon which the results of the previous section were based. The latter figure, which is again a plot of A1,

121 3.0 TT T OT Q T a,.~ 0 0 I 12 2 3 3 TIME (DAYS) Fig. 4.29. Time evolution of A'2. Random forcing of the second kind. Notation is the same as that in Fig. 4.17.

122 shows a more substantial increase in the standard deviation of the expected value throughout the forecast period. To the extent that all scales now suffer a reduced predictability limit, the second parameterization of external error growth represents a decided improvement over the first. Further support for the superiority of this procedure comes from examination of Fig. 4.30, a scatter diagram of apparent error versus standard error as inferred from forecast and observational uncertainty. The points are shifted generally to the right from those as given in Fig. 4.26. Although this outcome is not definitive by itself, the inference which may be drawn is that the forecast uncertainty is somewhat more reflective of the actual uncertainty. In equivalent terms, the ensemble of states as described by the stochastic forecast enjoys a greater consistency with the corresponding ensemble defined by the new observations. Several forecasts have been made and the effects of this particular prescription of external error growth may be conveniently summarized. As pointed out previously, the present formulation does admit external error growth in all the wave components retained in the model. Moreover, due to the nature of the random forcing, which provides for a uniform level of noise input on all scales, the smaller scales become uncertain more rapidly than

123 3.5e 3.0 2.5 u 2.0 I /S i.5 e=.675o1.0 0.. 0 0.5 1.0 1.5 2.0 o (10-8 Sec' ) Fig. 4.30. Scatter diagram of apparent forecast error e vs a, the standard error as inferred from observational and forecast uncertainty; o=(o2 + ao)1/, with symbols defined in section 4.1. The apparent error is defined as the magnitude of the difference between the least squares estimate and forecast value at Day 2. Random forcing of the second kind.

124 the larger ones. This is in general agreement with experience, and is reassuring. Random forcing of the second kind, while a definite improvement over the previous method, cannot be considered as the answer to the problem of parameterization of external model error. Because of the inherent simplicity of the approach, the forecast uncertainty of some model parameters does not reflect the actual uncertainty. We have noted on previous occasions that this state of affairs is detrimental to a valid analysis. Subsequent forecasts tend to drift away from reality. This is not a fault of the stochastic dynamic method, but arises from the simplicity of the basic physical model and the inability of a random forcing function, of the rather simple form exploited here, to account for the additional forecast uncertainty. An example of this is provided by the following four figures. The analyzed expected state on Day 3 is given in Fig. 4.31. Comparison of the subsequent forecast expected state (Fig. 4.32) with the least squares estimation at Day 4 (Fig. 4.33) reveals certain discrepancies not accounted for by the forecast uncertainty field (not shown here). For example, the forecast expected position of the trough-ridge system over the western portion of North America lags behind the best estimate as deduced from observations. Because these discrepencies are not reflected

Fig. 4.32. Forecast expected state of Fig. 4.31. Stochastic analysis of expected height field, Day 4. Random forcing height field, Day 3. of the second kind. Fig. 4.34. Stochastic analysis of expected Fig. 4.33. Teast squares estimation, Day 4. height field, Day 4.

126 in the forecast uncertainty, the new observations at Day 4 are effectively ignored in arriving at an analysis (Fig. 4.34). It should be evident that a better parameterization of external model error is necessary in order to take optimum advantage of the stochastic dynamic method with the present model. Indeed, one of the major drawbacks in this particular study has been the rather simple physical model used as a basis for prediction. Nevertheless, it would seem that in these initial stages the use of a simple model has permitted a wider range of experimentation and has led to a greater insight into the method itself.

CHAPTER V SUMMARY AND OUTLOOK The goal of the present study has been the furtherance of stochastic dynamic prediction, especially toward the problem of weather forecasting. Recognizing that atmospheric observations define the initial conditions in a probabilistic sense leads to a probabilistic method for prediction. Moreover, this method provides the best forecast from the standpoint of least mean square error. In fact, the most general formulation requires that we forecast an ensemble of initial states, or more specifically, the probability density function characterizing the ensemble. This task is not feasible, nor, as it turns out, really necessary. We simplify the problem by deriving a set of prediction equations for the lowest order moments, the stochastic dynamic equations. In order to obtain a solution of these equations, we have seen that a closure approximation must be made with respect to higher order statistical moments. Previous investigations (Epstein, 1969b; Fleming, 1971a) indicate that, for short-range predictions on the order of a few days, the neglect of third moment quantities yields valid solutions. We have adopted this third moment discard approximation and have presented calculations in Chapter II supporting this assumption. The next phase of the investigation deals with a test 127

128 of the method on real atmospheric data. The physical model employed is the equivalent barotropic, and the initialization of the corresponding stochastic model requires knowledge of the ensemble means of the spectral coefficients as well as the variance-covariance statistics. These moments are obtained by a least squares fit of the 105 spherical harmonics to the data. In this connection, we have already mentioned a couple of refinements which may prove profitable in future investigations. These include more realistic statements of the variances and, especially, the covariances between the observations themselves. Having given an appropriate specification of the initial conditions, we are now able to carry out the stochastic forecast. This yields an expected state and covariance matrix, which provides a measure of the uncertainty associated with that state. Because of the various limitations in the physical model, we have seen that during this phase of the calculations a parameterization of the growth of uncertainty due to these inadequacies must be incorporated. We have utilized two parameterizations. The second gives the more realistic results, but is only partially successful. This is due, in part, to the simplicity of the forcing functions adopted. The limitations of the barotropic model itself combine to make a solution of this problem a rather difficult one.

129 In section 3.5 we hinted at the implications of external model error with respect to forecast models in general. The more sophisticated the physical model, the more accurate will be the prediction of the ensemble mean. It would also appear reasonable to speculate that the increased number of error interactions would contribute to a further spreading of the ensemble members. This would be reflected in terms of increased forecast uncertainty. Both of these effects —better prediction of the ensemble mean and enhanced error growth —would conjoin to make for a greater consistency between the stochastic forecast and the new observations. Nevertheless, the question of external model uncertainty will always be important, although to varying degrees, and further work is required in this area. Although the computational requirement for the stochastic dynamic method is significantly greater than in the case of conventional procedures, there are some indications that extensions to more sophisticated models may be possible with the present computing capabilities. For example, it is possible that time derivatives of every covariance quantity might not have to be computed every time step. This would mean a substantial saving. Alternatively, the use of Monte Carlo methods may prove to be the best approach in the long run, if the sample size required is not too large. A recent theoretical study, based on a two-dimensional homogeneous isotropic

130 turbulence model (Leith, 1974), suggests that a sample size of 8 will give acceptable results in the forecast of the mean for that model. Whether this would be acceptable in an atmospheric model awaits further investigation. Moreover, such a small sample size might not be suitable in order to yield estimates of covariance statistics which would be of sufficient accuracy for optimum data assimilation. We conclude by noting that several problems remain to be solved before stochastic prediction can be considered as a viable alternative to existing methods, but this study has established the feasibility of the approach.

APPENDIX A PROOF OF STOCHASTIC ANALYSIS ALGORITHM This appendix demonstrates how one may deduce the analysis algorithm (3.9) by application of a result in probability theory called Bayes' theorem. It is not our purpose to give a complete treatment here since this may be found elsewhere, notably Raiffa and Schlaifer (1972), but instead to provide the major steps in the argument. The reader desiring fuller detail is referred to the rather exhaustive work cited above. Consider a random vector B whose elements are distributed according to a multivariate density function D1(B), herein referred to as the prior. Consider next a set of observations Y which bear some statistical relationship to B. Then from Bayes' theorem we can deduce D2(BIY), the posterior density function for 5, given the new observations: (fly) = (y1 D (p() R. (Al) In the above, R is simply a normalization factor and is determined by the condition, n 4( ) I'V N te rZ E(,VI a) | (A2) where the integrations are performed over all _. We de131

132 note by Z(YIB) the density function which expresses the conditional probability of Y given S. It is at this point that the relationship between Y and B should be specified. As discussed in Chapter III, we have chosen the regression model, Y- X p. r (A3) where in this application the elements of e are taken to be distributed according to a multivariate normal density with mean zero and covariance matrix a2I. The vectors Y and c have P elements corresponding to P observations. Having postulated the statistical nature of ~, we may write its probability density as follows:? -2 {(e)._~~~~~ tA ay _(A4) where the prime denotes transposition. But from (A3) i(Y!S) is simply (A5) As suggested by the notation Z(Yj_1) is merely the likelihood function for the data. Maximizing this function with respect to B gives the least squares solution _. To utilize (Al) we must now specify the nature of the

133 prior D1(i). This point was discussed in subsection 3.1.2, and on the basis of that argument D1(B) is chosen to be the normal density. The moments required to define this function are supplied by the stochastic forecast. Let b and S denote the prior mean and covariance respectively. If b contains N parameters then the prior is given by'D 12)7 2 s'lT X 2 b)'p S-((A6) where IS-11 is the determinant of S-1. Substitution of (A5) and (A6) into (Al) yields the posterior density: (A7) where r- C'~t y- _ 1)'( X- X )+ ~(- b)' s-'( ) -b)(A8) Expansion of the above and rearranging yields -_ —' ('_-' x_' _- Y._' _'_I _X )

134 - go" Yt + b' S' 6. (A9) The next step is to factorize: T-0"(S'o-d':x'. ) [}' (Sf't+'x''xr'(s",- x._ YSWI] [( S"+ 1' XX" x) ('+.".X'.Y)]'( S " -+~' X)' B' l - (s'+6 rxY)' S I+xll X)Y'111- (S"'W 2x4x)'ts1 1 ateIIY)]a(S-' e~ x x) [_Sp- ( s-+, X)"(s' b + a'" x'Y] +6 y'y t b' s" b -(s'X b + ~t',Y) (- + x'w)' ('b'+6x' This rather cumbersome expression may be written in the following relatively simple form, T, T + T, where - T. - 0 b) "$ - ), [ ~ 6.' Y'Y b S-' -6')]''

135 and 6 — s(s' b+ 2x'Y) ~.~m~6Z xy)(A10) Note that the random variable B appears only in T1, and that the latter is in the standard canonical form for the normal density. Thus D2(IY) is multivariate normal with posterior means and covariances given by (A10), which is the required result.

APPENDIX B EVALUATION OF INTERACTION COEFFICIENTS In section 3.4 the elements of the interaction matrix were defined as follows: PL', tM = c fe t1V Qp Jb *r P1,PSe) ae, (B1) in which p + r = m. Analytical expressions for the evaluation of (B1) have appeared in the literature (Silberman, 1954; Thiebaux, 1971). With a moderately large spectral truncation, these formulas are often not the most convenient to use. We shall evaluate (B1) by Gaussian quadrature. An exact evaluation, apart from round-off errors, is possible if a sufficiently high order quadrature is employed. This is possible because the integrand in (B1) can be written as a polynomial in P, the cosine of colatitude. The derivatives which appear in the above may be expressed in terms of the associated Legendre functions via a recursion formula. This is relatively straightforward by noting that the pm n may be defined in terms of the ordinary Legendre polynomials (normalized), i.e. p,' (~-L )l) (/I (B2) 136

137 Differentiation of the above with respect to 0 yields, OL P.- t qx. PY A {"I (fAn V)(n OrtW +a Y (B3) Making the appropriate substitutions in (B1) and simplifying gives, e E -~ flo Se(B4) where and similarly, E Psw Ps r+. A 32 point Gaussian quadrature was used to make a one-time evaluation of (B4). The majority of the elements Lpmr qns are zero with the necessary conditions for non-zero values given by selection rules (Orszag, 1970).

APPENDIX C RELATIONS BETWEEN REAL AND COMPLEX QUANTITIES In the sequel we shall present expressions which relate real and complex covariances. Complex covariances have been formally defined by (3.29) and are useful in the present context from both a mathematical and numerical standpoint. Of ultimate interest, however, will be the corresponding covariance information for real numbers. Let the following equation define the real and imaginary parts of m tes~ > +; be. ~~~(C1) Again, by definition, E n =pn = Am + iBn where as before pm is the ensemble average of m. In addition, the realn n ity of the streamfunction field implies that (C2) where the star denotes complex conjugation. From (3.29) we have ~vA -- Ev) WeI v 138

139 E( k%, zh V - E( b v) - Be RV n r % V 1, (C3)' ({ E ( ~, 6v ). Bv t+ E (ov b, ) - l I,,C from which results the following identity: Equation (C4) implies directly that VA, -*t.-,,v ), (C5) AL~) V ) 6: or(o ) F u( ) * (C6) Utilizing (C2) we may derive the following relation analogous to (C4): 6 tv + tnv'r,)( + CaO wi''I "1 V m (C7) where J = covy (1n, u). Likewise the foregoing yields the two additional identities:

(113) v im Ws - h Mt T ( )'ppo n (Z J.,A(013) Au?j 1 as ( u:sessD o; aoqt lxo; (63) pus' (8:))'(9D)' (So) mIox; poonpop eq X~UI suo;e.IO;u~AoIO2 OT, (63) TO )d )mb() 0t"1 Ih "k T VAti

APPENDIX D NUMERICAL METHOD FOR TIME INTEGRATION Interchanging summation indices in (3.35) and noting the very simple property, L rmp LPmr we obtain the snq qns' following symmetric form: Utilizing the above makes for increased computing economy. Each member of the system of equations to be integrated may be written schematically as x~R +~FX _ g (D2) where R represents all inhomogeneous terms and fx the linear ones. Following Baer and Simons (1970), we note that a reduction in truncation error is possible if we account for the linear terms exactly. We shall therefore deal with the following equation: t (X (- )) R141 (D3) 141

142 The total number of equations to be integrated is quite large, and this dictates the use of a relatively simple time differencing scheme. For this reason, the Adams-Bashforth scheme was chosen to carry out the time integrations. This scheme is also free of any computational mode (Young, 1968). With this approximation the finite difference form of (D3) is X- 3t'i) — t~.,.....L................ i It or equivalently t+bt %~tf Stil t t-bt (D4) In the above, the time step At was chosen as 0.5 hr. This assured that accumulated truncation errors were kept within tolerable limits. For the first time step an implicit forward difference was used. The relevant difference equation is rbt [ttj btt +b t -ttl|r( (D5) The preceding was solved by an iterative procedure employing an ordinary forward time step to start the iteration. Rapid convergence was obtained after about three iterations.

APPENDIX E FORMULATION FOR RANDOM FORCING OF THE SECOND KIND Parameterization of model uncertainty by random forcing of the second kind, as defined in subsection 3.5.1, yields the following modified form of (3.19): *M. le perts )Li (El) where and where it is understood that a stands for e(2). In this developmenta is a complex random function of time for m > 0, whereas the imaginary part of a is identically zero for m = 0. This constraint is simply a consequence of the fact that n0 is real for all n. The time evolution of the ensemble average Pm is obtained as before by taking the ensemble average of (El): # tX s [.|( t ~ ),s ~!, s4 L, $ + F, +., (E2) In the above and succeeding expressions the notation is the same as that employed in section 3.6. By contrast with (3.35) no second moment statistics involving a appear in (E2).'The 143

144 reason for this of course is because the random forcing enters in a statistically linear fashion, unlike the previous parameterization which is nonlinear in a statistical sense. In the manner leading to (3.36) we may derive another set of equations, each describing the time behavior of the covariance between a pair of spectral coefficients. The following results: WVF. CV $ S +DYW ((V) vL ~(L, ( (3 O"'D. v $ o'sts t+q, <r + t:s,) - C s L;,,, + F% r',V + D tr(Z n)e (E3) Obviously the only difference between (E3) and (3.36) relates to those terms involving the statistics of a. However, for the reason given in the previous paragraph, a fundamental difference between these two equations is the absence of third moment quantities involving a. Because a is real whenever it appears in an equation governing a zonal coefficient, caution must be exercised in u m programming (E3). The form of cov(a'v ) or cov(a'n may change according as m and/or uequal zero. Specifically, if a = ar +iei, where ar is the real and ai. the imaginary

145 part of a, then it follows that Ccnr t n) - Jru(it} t ) bt(di, + ( nC) where individually cov(ar, xn) and cov(ai, 4nm) may be complex. The second term on the right-hand side of (E4) is present in (E3) only when u > 0. Similarly the derivation of the equations governing the time dependence of cov(ar, 4nm) and cov(ai, nm) is quite straightforward. As in the foregoing, we quote the results: +om'((cKi Lt Ev Co.)od) Cot 0 tre cv(t, cxj) does Cnotapr in ES), n or t do (,va ap.-i (E6,)or, ( + D~ mm t c,,j" t), cr4 (E5) does ~var li appear in (E6).

146 Under the assumption that each of ar and ai obeys the simple stochastic process given by (3.22), then, as in subsection 3.5.2, we may obtain the following expressions: o, E- t + (tIo-. ) o(-t/t) ~~ma.~t- 0,2Z~~~~~~ (E7) 0, - m ( t- ij(-t/rt)) i+ O.4(-lt/t) w,here,= cnir 0 (-I t ) ii where B0= cov(r, ai)lt=O In the above 2l(o2l) is the long term variance of ar(ai). We shall take or = 1 The remaining symbols have been defined in subsection 3.5.2.

REFERENCES Baer, F., 1964: Integration with the Spectral Vorticity Equation. J. Atmos. Sci., 21, 260-276. and T.J. Simons, 1970: Computational Stability and Time Truncation of Coupled Nonlinear Equations with Exact Solutions. Mon. Wea. Rev., 98, 665-679. Charney, J.G. and A. Eliassen, 1949: A Numerical Method for Predicting the Perturbations of the Middle Latitude Westerlies. Tellus, 1, 38-54., M. Halem and R. Jastrow, 1969: Use of Incomplete Historical Data to Infer the Present State of the Atmosphere. J. Atmos. Sci., 26, 1160-1163. Cressman, G.P., 1958: Barotropic Divergence and Very Long Atmospheric Waves. Mon. Wea. Rev., 86, 293-297., 1959: An Operational Objective Analysis System. Mon. Wea. Rev., 87, 367-374. Davis, H.T., 1962: Introduction to Nonlinear Differential and Integral Equations. New York, Dover Publications, 566 pp. Draper, N.R. and H. Smith, 1966: Applied Regression Analysis. New York, John Wiley & Sons, Inc., 407 pp. Eddy, A., 1964: The Objective Analysis of Horizontal Wind Divergence Fields. Quart. J. Roy. Meteor. Soc., 90, 424-440., 1967: The Statistical Objective Analysis of Scaler Data Fields. J. Appl. Meteor., 6, 597-609. Eliasen, E., B. Machenhauer and E. Rasmussen, 1970: On a Numerical Method for Integration of the Hydrodynamical Equations with a Spectral Representation of the Horizontal Fields. Report No. 2, Institute for Theoretical Meteorology, University of Copenhagen, 35 pp. Epstein, E.S., 1962: A Bayesian Approach to Decision Making in Applied Meteorology. J. Appl. Meteor., 1, 169-177., 1969a: The Role of Initial Uncertainties in Prediction. J. Appl. Meteor., 8, 190-198., 1969b: Stochastic Dynamic Prediction. Tellus, 21, 739-759. 147

REFERENCES (continued), 1971: Stochastic Prediction of Deterministic Models. Technical Report 037430-2-T, Dept. of Atmos. and Oceanic Sci., The University of Michigan, 22 pp. and R.J. Fleming, 1971: Depicting Stochastic Dynamic Forecasts. J. Atmos. Sci., 28, 500-511. and E.J. Pitcher, 1972: Stochastic Analysis of Meteorological Fields. J. Atmos. Sci., 29, 244-257. Ericson, W.A., 1969: A Note on the Posterior Mean and Variance of Multiple Regression Coefficients. Statistical Res. Lab. Tech. Note (unpubl. ms.), The University of Michigan, 3 pp. Fleming, R.J., 1970: Concepts and Implications of Stochastic Dynamic Prediction, NCAR Cooperative Thesis No. 22, The University of Michigan and the Laboratory of Atmospheric Science, NCAR, 171 pp., 1971a: On Stochastic Dynamic Prediction, Part I, The Energetics of Uncertainty and the Question of Closure. Mon. Wea. Rev., 99, 851-872., 1971b: On Stochastic Dynamic Prediction, Part II, Predictability and Utility. Mon. Wea. Rev., 99, 927-938., 1973: Gravitational Adjustment in Phase Space. J. Appl. Meteor., 12, 1114-1128. Gandin, L.S., 1963: Objective Analysis of Meteorological Fields.Leningrad, Gidromet.; Jerusalem Israel Program for Scientific Translations; available from U.S.Dept. of Commerce, Clearinghouse for Federal Scientific & Technical Info., Springfield, Va., 1965, 242 pp. Gibbs, J.W., 1902: Elementary Principles in Statistical Mechanics. New York, Scribner, 207 pp. Gleeson, T.A., 1970: Statistical-Dynamical Predictions. J. Appl. Meteor., 9, 333-344. Jastrow, R. and M. Halem, 1970: Simulation Studies Related to GARP. Bull. Amer. Meteor. Soc., 51, 490-513. Jenkins, G.M. and D.G. Watts, 1968: Spectral Analysis and its Applications. San Francisco, Holden-Day, p. 157.

149 REFERENCES (continued) Kasahara, A., 1972: Simulation Experiments for Meteorological Observing Systems for GARP. Bull. Amer. Meteor. Soc., 53, 252-264. Knudsen, J.H., 1973: Prediction of Second-Moment Properties in Spectral Form for Non-Divergent Barotropic Flow. Geofysiske Publikasjoner, 30, No. 3, 1-14. Kruger, H.B., 1969: General and Special Approaches to the Problem of Objective Analysis of Meteorological Variables. Quart. J. Roy. Meteor. Soc., 95, 21-39. Leith, C.E., 1974: Theoretical Skill of Monte Carlo Forecasts. To be published. Lorenz, E.N., 1960: Maximum Simplification of the Dynamic Equations. Tellus, 12, 243-254., 1963: Deterministic Nonperiodic Flow. J. Atmos. Sci., 20, 130-141., 1965: A Study of the Predictability of a 28-Variable Atmospheric Model, Tellus, 17, 321-333. 1969a: Atmospheric Predictability as Revealed by Naturally Occurring Analogues. J. Atmos. Sci., 26, 636-646., 1969b: The Predictability of a Flow which Possesses Many Scales of Motion. Tellus, 21, 289-307. Merilees, P.E., 1968: The Equations of Motion in Spectral Form. J. Atmos. Sci., 25, 736-743. Orszag, S.A., 1970: Transform Method for the Calculation of Vector-Coupled Sums: Application to the Spectral Form of the Vorticity Equation. J. Atmos. Sci., 27, 890-895. Petersen, D.P., 1968: On the Concept and Implementation of Sequential Analysis for Linear Random Fields. Tellus, 20, 673-686. Platzman, G.W., 1960: The Spectral Form of the Vorticity Equation. J. Meteor., 17, 635-644. Raiffa, H. and R. Schlaifer, 1972: Applied Statistical Decision Theory (Student ed. ). Cambridge, Mass., The MIT Press Paperback Edition, 356 pp.

150 REFERENCES (continued) Richardson, L.F., 1922: Weather Prediction bL Numerical Process. London, Cambridge University Press, 236 pp. (Reprinted by Dover.) Rinne, J., 1970: Investigation of the Forecasting Error of a Simple Barotropic Model with the Aid of Emperical Orthogonal Functions, Parts I and II. Geophysica, 185-213, 215-239. Rutherford, I.D., 1972: Data Assimilation by Statistical Interpolation of Forecast Error Fields. J. Atmos. Sci., 29, 809-815. Salzer, H. E., R. Zucker and R. Capuano, 1952: Table of Zeros and Weight Factors of the First Twenty Hermite Polynomials. Journal of Research of the National Bureau of Standards, 48, 111-116. Sasaki, Y., 1970: Numerical Variational Analysis Formulated Under Constraints as Determined by Longwave Equations and a Low-Pass Filter. Mon. Wea. Rev., 98, 884-898. Silberman, I., 1954: Planetary Waves in the Atmosphere. J. Meteor., 11, 27-34. Smagorinsky, J.S. and K. Miyakoda, 1969: The Relative Importance of Variables in Initial Conditions for Numerical Predictions. Proc. WMO/IUGG Symp. Numerical Weather Prediction, Tokyo, 1968. Soong, T.T., 1973: Random Differential Equations in Science and Engineering. New York, Academic Press, 327 pp. Srinivasan, S.K. and R. Vasudevan, 1971: Introduction to Random Differential Equations and Their Applications. New York, American Elsevier Publishing Company,Inc., 166 pp. Tadjbakhsh, I.G., 1969: Utilization of Time-Dependent Data in Running Solution of Initial Value Problems. J. Appl. Meteor., 8, 389-391. Thiebaux, M.L., 1971: On the Structure of Interaction Coefficients in the Spectral Equations for Planetary Waves. J. Atmos. Sci., 28, 1294-1295.

151 REFERENCES (concluded) Thompson, P.D., 1957: Uncertainty of Initial State as a Factor in the Predictability of Large Scale Atmospheric Flow Patterns. Tellus, 9, 275-295., 1961: Numerical Weather Analysis and Prediction. New York, The Macmillan Company, p. 87. Tolman, R.C., 1938: The Principles of Statistical Mechanics. Oxford, The Clarendon Press, 660 pp. Tribus, M., 1969: Rational Descriptions, Decisions and Designs. New York, Pergamon Press, 478 pp. Uhlenbeck, G.E. and L.S. Ornstein, 1930: On the Theory of the Brownian Motion. Phys. Rev., 36, 823-841. Wellck, R.E., A. Kasahara, W.M. Washington and G. DeSanto, 1971: The Effect of Horizontal Resolution in a Finite-Difference Model of the General Circulation. Mon. Wea. Rev., 99, 673-683. Williamson, D.L., 1973: The Effect of Forecast Error Accumulation on Four-Dimensional Data Assimilation. J. Atmos. Sci., 30, 537-543. and A. Kasahara, 1971: Adaptation of Meteorological Variables Forced by Updating. J. Atmos. Sci., 28, 1313-1324. Winkler, R.L., 1972: Introduction to Bayesian Inference and Decision. New York, Holt, Rinehart and Winston, 563 pp. Young, J.A., 1968: Comparative Properties of Some Time Differencing Schemes for Linear and Nonlinear Oscillations. Mon. Wea. Rev., 96, 357-364.

~o-t ~Io, C ce