THE UNIVERSITY OF MICHIGAN COLLEGE OF ENGINEERING DEPARTMENT OF METEOROLOGY AND OCEANOGRAPHY Technical Report ON POWER LAWS AND NONLINEAR CASCADES IN LARGE SCALE ATMOSPHERIC FLOW I. Leonard Steinberg Aksel C. Wiin-Nielsen Project Director ORA Project 002650 supported by: NATIONAL SCIENCE FOUNDATION GRANT NO. GA-16166 WASHINGTON, D.C. administered through: OFFICE OF RESEARCH ADMINISTRATION ANN ARBOR March 1971

ACKNOWLEDGMENTS I thank Professors A. C. Wiin-Nielsen, E. S. Epstein, W. P. Graebel, and Adjunct Associate Professor W. M. Washington for participating in my dissertation as committee members. In particular, I am grateful to Professor W. M. Washington for his assistance in obtaining computer time at N.C.A.R. and his helpful discussions during my stay in Boulder and Professor E. S. Epstein for useful discussions and comments on the preliminary drafts of the manuscript. I am very much indebted to Professor A. C. Wiin-Nielsen for his guidance and many suggestions throughout the entire duration of this study in his dual role as Committee Co-Chairman and Project Director. Acknowledgments and thanks are due to James Pfaendtner who programmed the calculations reported on in Chapter III and was most helpful with the debugging of the model, and to Margaret Drake of N.C.A.R. whose help in making the programs compatible with the N.C.A.R. computer was invaluable. My thanks are due to the secretaries of the Department of Meteorology and Oceanography for the preliminary drafts of this report and to the Office of Research Administration for its final production, and in particular Mary L. Falcone. Acknowledgment is made to the National Science Foundation for its sponsorship of this project as well as the National Center of Atmospheric Research, which is sponsored by the National Science Foundation, for use of its Control Data 6600 computer. ii

TABLE OF CONTENTS Page LIST OF TABLES v LIST OF FIGURES vii LIST OF APPENDICES x ABSTRACT xi I. INTRODUCTION 1 1.1 Previous Theoretical Studies1 1.2 Previous Observational Studies6 1.5 Previous Numerical Studies 7 1.4 Outline of the Study 8 II. FORMULATIONS AND DERIVATIONS 14 2.1 Formulation of the Computations 14 2.2 Decomposition of Exchanges and Fluxes into Temporal Modes 23 2.3 Computational Procedures 26 III. DISCUSSION OF RESULTS 29 3.1 Introduction 29 3.2 Three Monthly Means of Exchanges and Fluxes 29 3.3 Results in the Time Domain 37 3.4 Exchanges and Fluxes as a Function of Pressure 50 3.5 Enstrophy and Potential Enstrophy Budgets 55 3.6 Spectral Shape 58 3.7 Coefficient of Eddy Viscosity and Related Quantities 65 3.8 Transfer Coefficient and Potential Vorticity Transport 71 IV. THE NUMERICAL MODEL 76 4.1 Introduction 76 4.2 Formulation of the Model 78 4.3 Numerical Methods 8 4.4 Energetics of the Model 88 4.4.1 Formulation 88 4.4.2 The Diagnostic Calculations 92 4.4.3 The Energy Cycle and Nonlinear Transfers 97 4.4.4'Power Laws' 105 iii

TABLE OF CONTENTS (Concluded) Page V. CONCLUDING REMARKS 111 5.1 Summary and Conclusions 111 5.2 Suggestions for Future Investigation 116 BIBLIOGRAPHY 118 iv

LIST OF TABLES Table Page 1. Static Stability, a(p) (T2 m sec ) 27 2A. Exchanges and Fluxes for the Enstrophy for the Wave Groups in the Time Domain for the Atmosphere (x 10-18 sec-3) 38 2B. Exchanges and Fluxes for the Potential Enstrophy for the Wave Groups in the Time Domain for the Atmosphere (x 10-18 sec-3) 38 3A. Enstrophy (corrected), the Time Mean Change and Residue 56 3B. Potential Enstrophy (corrected), the Time Mean Change and Residue 57 4A. c3 for [ ]S[ [cp]s[p]s and [pc]s(Wave Numbers 1 to 7) 60 4B. c3 for [I]S, [ic]S' [tp]S and [tpc]S (Wave Numbers 8 to 15) 60 5A. c for the Transient enstrophy and potential enstrophy for the spectral region 1 < m < 7 (x 10-1) 63 5B. c for Wave Numbers 8 to 15 for the Transient Enstrophy and Potential Enstrophy (corrected) 63 -25 -2 6A. (V7)2 The Square of the Vorticity Gradient [x 1025 cm sec-2] 67 6B. (V~p)2 The Square of the Potential Vorticity Gradient [x 10-25 cm-2 sec-2] 68 6C. (Potential) Enstrophy Cascade 68 8 2 -i 7. Coefficient of Eddy Viscosity [x 10 cm sec ] 69 8. Heating Amplitude and Coefficient of Eddy Viscosity 93 9. Changes in Total Energy Compared with Residue of Energy Budget Calculations, (x 109 cm2 sec-2) 94 Oat 10. Evaluation of C and a in Ke = C e, 600 < t < 1200, Using Ke at t = 600 and 1200 99 V

LIST OF TABLES (Concluded) Table Page 11. Kinetic Energy Spectral Slopes for 8 < m < 16 as a Function of Time 105 Al. Quantities Used in Potential Vorticity Calculations 124 C1. Means and Standard Deviations of Processing Affecting Enstrophy (x 10-18 sec-3) 127 C2. Means and Standard Deviations of Processing Affecting Potential Enstrophy (x 10-18 sec-3) 128 vi

LIST OF FIGURES Figure Page 1. Vertical resolution of the observational calculations. 28 2. Plotting model for Figures 3a and 3b. 30 3a. Three monthly means of rates of exchanges and fluxes of enstrophy [x 10-18 sec-3] in the atmosphere (0 ~ 100 cb). 31 3b. Three monthly means of rates of exchanges and fluxes of potential enstrophy [x 10-18 sec-3] in the atmosphere 32 (12.5 cb < p < 92.5 cb). -18 -3 4. I(min,A) [x 10 sec ] in the time domain. [Legend is used for Figures 4-13.] Heavy line indicates net. 39 5. I (mjn,k) [x 10-18 sec-3] in the time domain. 40 6. F(mjn,2) [x 10-18 sec-3] in the time domain. 42 7. F (mln,2) [x 10-18 sec-3] in the time domain. 45 8. I(m(o) [x 10-18 sec-3] in the time domain. 44 9. I (mlo) [x 10-18 sec53] in the time domain. 45 10. F(mjo) [x 10-18 sec-3] in the time domain. 46 11. F (mlo) [x 1018 sec-3] in the time domain. 47 12. P(m) [x 10-18 sec-3] in the time domain. 48 13. P (m) [x 10-18 sec-3] in the time domain. 49 14. (a) Vertical variation of F(min,) [x 10-18 sec-3]. (b) Vertical variation of I(mlo) [x 10-18 sec-3]. (c) Vertical variation of F(mto) [x 10-18 sec-3]. (d) Vertical cariation of P(m) [x 10-18 sec-3]. 51 15. (a) Vertical variation of Fp(m n,l) [x 10-18 sec-3]. (b) Vertical variation of Ip(mlo) [x 10-18 sec-3]. (c) Vertical variation of Fp(mlo) [x 1018 sec-S]. (d) Vertical variation of Pp(m) [x 10-18 sec-3]. 52 vii

LIST OF FIGURES (Continued) Figure Page 16. (a) Vertical variation of sum of all processes affecting the enstrophy [x 10-18 sec-3]. (b) Vertical variation of the sum of all processes affecting the potential enstrophy Ix 10-18 sec-3]. 53 16. (c) Vertical variation of the three monthly means of enstrophy and potential enstrophy [x 10-18 sec-3]. 54 17. (a) Spectrum of [~c]s. (b) Spectrum of [1Pc]s. 61 18. (a) Spectrum of the three monthly mean transient enstrophy. (b) Spectrum of the three monthly mean transient potential enstrophy. 64 19. [v pl],s [x 1011 cm sec-'] mean of February, March and April, 1963. 74 20. Schematic of model's vertical resolution, equations, and vertical boundary conditions. 79 21. Illustration of cyclic lateral boundary conditions on a 5 x 5 grid. 85 22. Energy flow diagram for the model. 88 25. Ke and Ae as a function of time for Case A. 98 24. Mean values of CA(mIn,1), -C(Ae,Ke), and D(Ae) for 300 < t < 1200, Case A. 101 25. Mean values for CK(m|n,Q), C(Ae,Ke), and D(Ke) for 300 < t < 1200, Case A. 102 26. Mean value of I(mln,2), 300 < t < 1200, Case A. 104 27. The kinetic energy of wave numbers 8, 12, 16, and 24 as a function of time, Case A. 106 28. The kinetic energy spectrum number at t = 600 and 1200, Case A. 107 29. The kinetic energy spectra at t = 1200 for Cases A, B, and C. 109 viii

LIST OF FIGURES (Concluded) Figure Page 30. The spectra of Ae at t = 1200 for Cases A, B, and C. 110 E.1. Ke, Ae as a function of time for Case B. 132 E.2. Ke, Ae as a function of time for Case C. 133 E.3. Mean values of CA(mln,e), -C(Ae,Ke), and D(Ae) for 300 < t < 1500, Case B. 134 E.4. Mean values of CA(mln,Q), -C(Ae,Ke), and D(Ae) for 300 < t < 1500, Case C. 135 E.5. Mean values of C(Ae,Ke), CK(mln,), and D(Ke) for 300 < t < 1500, Case B. 136 E.6. Mean values of C(Ae,Ke), CK(mIn,l), and D(Ke) for 300 < t < 1500, Case C. 137 E.7. Mean values of I(mJn,A), 300 < t < 1200, Case B. 138 E.8. Mean values of I(mJn,e), 300 < t < 1200, Case C. 139 E.9. The kinetic energy of wave numbers 8, 12, 16, and 24 as a function of time, Case B. 140 E.10. The kinetic energy of wave numbers 8, 12, 16, and 24 as a function of time, Case C. 141 E.11. Ke spectra at t = 900, 1500, Case B. 142 E.12. Ke spectra at t = 900, 1500, Case C. 143 ix

LIST OF APPENDICES Page A. VERTICAL FINITE DIFFERENCE SCHEME FOR THE POTENTIAL VORTICITY CALCULATIONS 123 B. RELATIONSHIP BETWEEN KINETIC ENERGY AND ENTROPHY IN WAVE NUMBER SPACE 125 C. MONTHLY MEANS AND STANDARD DEVIATIONS OF PROCESSES AFFECTING THE ENSTROPHY AND POTENTIAL ENSTROPHY 126 D. DERIVATION OF THE STREAM FUNCTION FROM THE POTENTIAL VORTICITY 129 E. ILLUSTRATION OF RESULTS FOR CASES B AND C 131 x

ABSTRACT Calculations of enstrophy and potential enstrophy changes are carried out in wave number space. The processes of interest are the interactions among waves, the interactions between each of the waves and the zonal flow, and the boundary fluxes. The framework used in the calculations is supplied by the quasi-geostrophic approximation to the equations of motion and the first law of thermodynamics. It is found that the enstrophy is cascaded from small wave numbers [1-7] toward large wave numbers [12-15], with relatively'little gain or loss in the middle range of wave numbers [8-11]. Verification of the existence of an inertial subrange is inconclusive from the calculations based on the atmoshperic data. Power spectrum analyses of the enstrophy and potential enstrophy yield slopes very close to -1, thus indicating a -3 slope for the associated kinetic energy spectrum. 8 2 -1 A value of v - 3 x 10 cm sec is obtained using Leit.h's formulation for the coefficient of eddy viscosity. The cut-off wavelength for the -3 range as formulated by Kraichnan is found to be 70 km. A two-level quasi-geostrophic model (similar to Phillips, 1956) is integrated for two values of the amplitude of the heating function and two values of the eddy viscosity coefficient. Quasi-stationary spectral distributions are obtained in each case. The'power law' xi

characterizing the spectra, varies somewhat from case to case. The nonlinear cascade of available potential energy results in gains at large wave numbers; thus corresponding to atmospheric observations. The nonlinear cascade of kinetic energy results in gains for the small wave number [1-7], and losses for 8 < m < 14 and relatively small gains for m > 14. Gains of enstrophy are observed to be fairly constant for 14 < m < 28. Relatively small gains are recorded for m < 7 while the spectral region 8 < m < 14 acts as the source. In view of this result it is concluded that studies involving enstrophy should extend well past m = 15 and the transition zone noted in the section on theobservational results would most likely not be extended even if more components were accounted for in the calculations. It is concluded that the two-dimensional turbulence theories (Fj~rtoft, 1953; Kraichnan, 1967; and Leith, 1968a) pertain to the atmosphere inasmuch as they predict the dominant direction of the nonlinear cascades of kinetic energy and enstrophy. Inertial subranges in which the cascade takes place resulting in zero gains or losses is not observed in either of the results of the observational or numerical studies. In view of these results it is concluded that the closeness with which atmospheric kinetic energy spectra follow a'minus three power law' should not be explained by arguments postulating inertial subranges. The quasi-stationary kinetic energy spectral slope seems to xii

be affected mainly by the relative intensity of energy sinks and sources at each wave number, and, to a lesser degree, by the nonlinear affects. xiii

CHAPTER I INTRODUCTION In recent years much effort has been put into the investigation of the nonlinear effects in the oceans and the atmosphere.* These processes are important in the redistribution, both spatially and in wave number space, of transport processes and various forms of energy. (Wave number may alternately be defined as the number of wave cycles per unit distance or wave cycles in the interval 0 < X < 2t, where \ is longitude. These quantities will be denoted by k and m, respectively.) Some particular phenomena in which these processes are intimately involved are the maintenance of the Jet Streamns in the atmosphere and the Gulf Stream in the ocean, frontogenesis, poleward heat transport in the troposphere and the limiting of the growth of the barocl. 1Bially unstable waves. 1.1 PREVIOUS THEORETICAL STUDIES The theory of three-dimensional turbulence predicts a universal shape for the spectrum of kinetic energy over those scales where there is statistical equilibrium, that is, there is no change in the kinetic energy of a particular scale relative to adjacent scales. A spectrum may be thought of as a series of Fourier amplitudes distributed among *A summary and extensive bibliography of work in this field is given by Reiter (1969a). 1

2 wave numbers. The shape of the turbulence spectrum is deduced from the assumption that the nonlinear energy transfer toward larger wave numbers is large as compared with the rate of change of energy in each wave number. The universal range lies between the scales at which turbulence is generated and those at which the viscous dissipation is dominant. The spectrum's universality implies that the form of the spectrum will be independent of external conditions. The spectrum is assumed to be a function of E, the assumed constant energy cascade rate, and independent of v, the viscosity. From dimensional arguments the form of the kinetic energy spectrum, K(k, ), the kinetic energy per unit mass and wave number, is deduced to be given by K(k(E) = 2/3 k-5/3 where C) is a constant. This is the Kolmogorov inertial subrange. In atmospheric flows the horizontal scales of the planetary and synoptic waves are about two to three orders of magnitude larger than the vertical scale. Away from the boundary layer the vertical component of the wind and its variations may in general be neglected when compared with the horizontal components. This idea is reflected by the use of the hydrostatic equation. Large-scale atmospheric flow may, therefore, to a first approximation be thought of as being quasitwo-dimensional.

In a two-dimensional flow of a nondivergent homogeneous fluid two conservative properties may be identified in the absence of frictional and generating effects. These quantities are the kinetic energy and any power of the vorticity. The kinetic energy is given by K = V dm (1.1) M where V is the horizontal velocity vector and dm is an elemental mass. The enstrophy is one half the mean square vorticity and has a convenient spectral form. It is given by 21 2 f = M ddm (1.2) M 2 M where M is total mass, and the vorticity is % = V1, where * is the stream function. The second power of the vorticity is used in the theories and observational calculations. In a relatively early study, Fj^rtoft (1953) deduced that in twodimensional turbulent flow, the kinetic energy flow in wave number space must be in one of two ways (1) kl + k2 -+ k or (2) k1 k2 k

4 where k < k < k3, in order for it to be conserved. ~ 2 5 More recent and complete theoretical considerations (Kraichnan, 1967; Leith, l968a,b) predict the power laws which should govern the kinetic energy spectrum in the regions of inertial subranges for two-dimensional turbulence. A "minus three" law governs the spectral region for which k > ki, where ki is the wave number at which energy is being fed into the system. This spectral region is characterized in the theories by a constant enstrophy cascade, Ar, toward large k and zero kinetic energy cascade. The part of the spectrum for which k < ki is characterized by a constant kinetic energy flow toward small wave numbers and zero enstrophy flow. A "minus 5/3" law is predicted to govern this spectral region. The following paragraphs contain some of the assumptions and deductions of the theories as presented by Kraichnan and Leith. If the turbulent fluid flow is quasi-two-dimensional and is to have a spectral region whose shape is preserved in all realizations, then the energy-producing mechanism and the dissipation must be on scales relatively remote from this region. Let L and L. represent the d wavelengths at which the dissipation and energy input, respectively, are taking place. Then, if L is some characteristic wavelength in the s shape-preserving spectral region, a condition for the existence of this region is Ld < L <K L.. ad s i The remaining effect is the nonlinear inertial transfer of kinetic energy or of enstrophy in wave number space. This transfer must be

5 constant with respect to wave number so that the various scales may retain their relative magnitudes. Suppose there exists a spectral range in which the nonlinear enstrophy cascade, rj, is the dominant process. Then the kinetic energy per wave number would depend on r and wave number, k. The kinetic energy per unit mass and wave number may then be written in the form K(k) = C2 T k (1.3a) where C2 is a nondimensional constant. The dimensions of K(k), j, and k are (length)3 (time), (time)-, and (length), respectively. The dimensional equation may be written aS (length) (time)2 = (time) (length)In order for (1.3a) to be dimensionally consistent, a = 2/3 and P = -3. Therefore the kinetic energy must be represented by K(k) = C2 2 k- (1.3b) -2 and the energy cascade, E, would have the form ~ = C ik2 where C3 is a constant. However, since the energy cascade cannot be a function of wave

6 number, C = O. Hence, a spectral range dominated by a constant enstrophy cascade should have the kinetic energy proportional to the -3 power of wave number and should exhibit zero kinetic energy cascade in wave number space. Analogously, a spectral range which is characterized by constant energy flux, e, will have kinetic energy dependent on the -5/3 power of wave number and zero enstrophy cascade. 1.2 PREVIOUS OBSERVATIONAL STUDIES Various observational studies have been carried out in order to ascertain the validity of the assumptions and conclusions of the theories. Ogura (1958), Horn and Bryson (1963), and Wiin-Nielsen (1967) have calculated the power law relating kinetic energy to wave number for various latitudes, pressure levels and times. It was generally found that the spectrum between wave numbers 7 and 15 tended towards a power law with exponent varying between -2 and -3, with mean values closer to -3 than to -2. At some latitudes and pressure levels the exponent is -5. Yang (1967) has made calculations of the nonlinear terms in the kinetic energy equation. He verified, in part, the conclusions of Fj~rtoft. It was observed that the kinetic energy, for the most part, flowed out of the middle waves (5 < m < 10) and into both the long and short waves. An indication from Kraichnan's theory that energy should

7 not pass through to larger wave numbers was also given support by Yang's (loc. cit.) result. Although some kinetic energy was transferred toward large wave numbers, part of this energy was, in turn, being given up to the zonal mode. Further, the gain in waves, m = 13, 14, and 15, is forced in part by the closure condition on the system of equations which was truncated at wave number 15. Julian et al. (1970), have studied kinetic energy spectra from a variety of sources and have concluded that "...for wavelengths < 4000 km a power-law representation would require a constant coefficient of between -2.7 and -35.0.... This coefficient apparently does not depend upon whether the spectrum is computed from geostrophic or balanced winds, or from observed winds." 1.5 PREVIOUS NUMERICAL STUDIES Lilly (1969) conducted a numerical experiment simulating twodimensional turbulence. He used a 64 x 64 point grid with cyclic boundary conditions in both the x- and y-directions. His forcing function was made up of all the two-dimensional components whose larger wave number equaled m.. The total amplitude was kept constant, but the contributions from the individual components were randomly chosen from a Gaussian distribution. He found good agreement with the power laws predicted by the theories. After a time energy began piling up in wave number one, since that is the largest mode allowed in the system. This resulted in

8 deviations from the -5/3 law for m < m.. Aliasing served to distort the short wave end of the spectrum. Manabe et al. (1970), have calculated the kinetic energy spectrum from their general circulation model and obtained a "slope (of)... approximately -3 power." Of the nonlinear effects on kinetic energy they conclude that "...most of the energy, which is converted by baroclinic instability in the medium wave number range is transferred to the low wave number range. 1.4 OUTLINE OF THE STUDY The major contribution of this study is the evaluation of the nonlinear cascade of (potential) enstrophy in wave number space, both from atmospheric observations and a numerical simulation of quasi-twodimensional turbulence. The enstrophy is defined by (1.2), while potential enstrophy is defined through the quasi-geostrophic equations (Phillips, 1963) to be e 1 J 1 2 + ff2 g 0 dM (1.4) where f is the Coriolis parameter (= 22 sin 0 ) at 0 = 45~, a is the static stability taken as a function of pressure only and p is pressure, the vertical coordinate. Taking a = o(p) and f = f, a constant, is consistent with the quasi-geostrophic equations and their integral constraints (WiinNielsen, 1959).

9 The nonlinear effects are, in the open system in which the observations are available, separated into the interactions taking place within the region and the boundary effects. Both processes are evaluated by the summation over allowable triads and are denoted by T(mln,~) (gain rate by wave number m due to interaction of triads composed of the allowed values of wave numbers n and ~) and F(mln,2) (gain rate by wave number m by boundary flux due to the permitted triads containing wave number m), respectively. A triad is a group of three wave numbers in which the sum of two must equal the third. The corresponding terms for the potential enstrophy are denoted by I (m|n,e) and F (mln,e). P P The theories already mentioned indicate that enstrophy transfer should be toward large wave numbers and through a spectral region in which the transfer is not a function of wave number, i.e., it is constant. Other processes involved in changing the (potential) enstrophy were calculated from the atmospheric data. They are I(mlO) and F(mlO) the gain rates of wave number m due to interaction with the zonal mode within the region and at the boundary, respectively. The P-effect (which reduced to a boundary flux effect in the case of enstrophy) is denoted by:(m). The corresponding quantities for the potential enstrophy are denoted by I (miO), F (m O) and P (m). Chapter II contains the formulation of the calculations involving the atmospheric data as well as some discussion of their physical meaning. The derivations used for the decomposition of the results in

10 the time domain is also contained in this chapter. The results of the computations using the atmospheric data are presented in Chapter III. The spectral shape of the (potential) enstrophy is discussed in terms of the power law, with respect to wave number, which best fits the data. If the kinetic energy is governed by a -3 power-law, then the enstrophy should obey a -1 power-law, as shown later. Since the potential vorticity is calculated for the computations mentioned above we have the opportunity to evaluate Green's (1970) suggestion concerning the transfer coefficient for this quantity. From these calculations we may deduce whether or not the action of the eddies may be represented in the following manner a (o)+f) [vp I -K ( 1-5). rp a, 6, where C (0) is the zonal potential vorticity Kp is the appropriate transfer coefficient v is the south-north velocity component a is the earth's radius 0 is latitude [{} ) indicates (1/2j) j (}d2 and () indicates the deviation of f) from [((]. (Square brackets with subscripts indicate a mean with respect to the subscripts while parentheses with the subscripts indicate deviations from the mean. This is a modified form of Reiter's (1969b) notation. ) For example any variable X may be written as follows

11 X(OX) = [X] X+ [(x) ] + [(x) ] + () (1.6) It is known that the transport of momentum does not obey a classical diffusion law of this kind and for the most part the flux is counter gradient. However, the momentum transport may be written as a combination of the heat and potential vorticity transports. The former transport has been shown to be generally of a classical nature (WiinNielsen et al. (1963, 1964)) in the troposphere i.e., heat is transported from a warmer region to a cooler region. If it can be shown that the potential vorticity also obeys this law then the action of the eddies in the momentum equation may be modelled without going into the complete calculations. Green (1970) gives reason to believe that this procedure will be successful. Wiin-Nielsen and Sela (1971) show that when monthly mean momentum and heat transports are used the transfer equation can be applied to the potential vorticity transport with K. (0,t) positive almost everywhere. We can compare the results p obtained here which are based on daily calculations with those of WiinNielsen and Sela (19'71) which are based on monthly and annual mean data. Further, in Chapter III we report on evaluations of various quantities such as the coefficient of eddy viscosity, v, the viscous dissipation scale, kd, the energy input and dissipation which correspond to the calculated value of I, the enstrophy cascade rate. Chapter IV contains the formulation of the numerical simulation

12 model, the equations, approximations and finite-difference methods. The numerical model, using essentially the framework described by Phillips (1956) and incorporating a heating field akin to the forcing function used by Lilly (1969), will be used to simulate disturbances by a quasi-monochromatic heating field. The disturbances will be confined to one of an infinite number of identical square regions. The lateral boundary conditions are, therefore, cyclic in both directions. There will be no zonal mode and so the generation will be that of eddy available potential energy. Kinetic energy will simultaneously be present since the existence of amplitude in the thermal field implies amplitude in the stream function field of at least one of the levels and most likely both. Charney (1966) discusses the energy cascade in three-dimensional flow under adiabatic, inviscid, and quasi-geostrophic conditions. He points out that if the potential temperature at the surface is a constant then an analogy may be made between the flow mentioned above and two-dimensional incompressible inviscid flow. The two conservative properties in the three-dimensional flow are the potential + kinetic energies and the potential enstrophy. These quantities are present in the model used in this study. Results for runs with different values of the amplitude of the heating function as well as the coefficient of viscosity are presented. It is hoped that such an investigation will help to bridge the gap between the two-dimensional turbulence

13 simulated by Lilly's model and the atmospheric turbulence. In Chapter V the major results are summarized and suggestions for future work are made.

CHAPTER II FORMULATIONS AND DERIVATIONS 2.1 FORMULATION OF THE COMPUTATIONS The purpose of these calculations is to determine the nonlinear transfer of vorticity and potential vorticity in quasi-nondivergent flows. The equations are formulated in a spherical geometry using Fourier series along lines of constant latitude and finite differences in the north-south direction. The stream function is specified by M *(k,Cp,pt) = Z A(m;c,p,t) e (2.1) m=-M where X is the longitude cp is the latitude p is the pressure t is the time m is the wave number (the number of cycles in the interval O < \ < 2r) M is a cutoff wave number = 15 i = a. A(m;cp,p,t) is a complex number. (We shall omit cp, p and t from the parentheses for brevity.) 14

15 Let A(m) = - (B(m) - i C(m)) (2.2) 2 and A(-m) = - (B(m) + i C(m)) (2.5) 2 Considering the contributions in (2.1) from m = m1 and m = -ml we get im l -iml A(m!) e 1 + A(-ml) e = B(ml) cos(m1)) + C(ml) sin(mlk) (2.4) which assures that j is always real. The relative vorticity may be written as = oVS 2 (2.5) S where V2 = 1 1 + 1 _ cos l aS 2 2- 2 cos4 6 L J a cos o a6% and a is the earth's radius. Using p = sin cp and 8/0p = cos cp 8/8p we get 1.21- 2 8k2 arbpaclJ = 2 _ a, + K (1 J l )j (2.6) Substituting from (2.1) we get = L{2 i J 1 1j e~ a (m) m i(27) m=-M L- _

16 or M; = Z F(m) e (2.8) m=-M where F(m) = 1 ) - A(m) (2.9) a, 2 21-~ and F(m) = ^ [G(m) - iH(m)]; F(-m) = [G(m) + iH(m)] (2.10) 2 2 Thus: G(m) = | a (1-L2) B(m (m (2.11) a.2 6p ( 21) and aH(m) = -1 a.(l-2) 2 C(m) 2 a, 6 ~-1 2 is give n in a latitude-longitude grid with grid distance AN = A~ = 2.5 = /72 radians. Therefore, 144 B(o) = l4 i 144 i=l B(m) = 72- [ti cos (mni)0 m O 1144 C(m) = 72 [=i sin (mXi)] i=l To calculate G(m) and H(m) at latitude < =. we use (2.11) in finite difference form as follows.

17 B jl(m) + Bj l(m) - 2B (m) G.(rnm) = j+1 j1 j^ —1 G (m) = i o a2 2 B (m) - B (m) m2 B.(m -tan j 2 Aa... (2.12) j 2 M) 2 cos j H(m) may be calculated from (2.12) when C(m)'s replace the B(m)'s. The lowest latitude for which G(m) and H(m) can be computed is p = 22.5~N and the highest latitude is cp = 85~N, since B(m) and C(m) are available for 20~N < cp < 87.5~N. Consider the nondivergent vorticity equation::= -V.V - V.Vf, = k xV7 (2.13) where V is the horizontal nondivergent wind f is the Coriolis parameter (= 2, cos cp, where Q is the earth angular velocity) and k is a, vertical unit vector. This may then be written as. _ 1 a a: 1 af a: at a, ap a cos cp a a cos cp a acp a d(2n sin cp) a cos cp \ a. ac or -1 -1 2 (2.14) at 2 Wax ] 2 a (2214) Substituting for' and ~ from (2.1) and (2.8), respectively, into -(2.14), multiplying by e, and integrating with respect to rom (2.14), multiplying by e, and integrating with respect to X from 0

18 to 23T it may be seen that only those terms for which n + I = m on the right-hand side will give a non-zero contribution. n and I are the wave numbers corresponding to f and 5, respectively. Carrying these operations out we get: aF(m) at 1 M FJrn-n) \A(n) - F(m-n) A ()f -i 22 m A(m) 2 i m-n) F(m-n) -A(n) - 2 a n=-M a E [mF(m-n) -n [ (F(m-n) A(n) - a 2Q m A(m) a n =-M a (2.15) Substituting (2.8) into (1.2) we get I' M - = L[, F(m) F(-m) = L F(m) F(-m) + F 2 m=_M 2 2 j=-M C m=l The rate of change of enstrophy may be derived as follows. Neglecting the zonal component in the left-hand side. d_ =, a F(m) F(-m) m=-M't = rF(-m) F(m) + F(m) 6F(-m) Now, 6F(-m) ot a n:-M J 2 F(-m-n) al) n a Fm-)A )i a2 2m A(-m) (2.16)

19 Thus (F(m) F(-m)} M cA(n) - Z m (F(-m) F(m-n) - F(m) F(-m-n)) n=-M Mi - 2 Z F(-m) (F(m-n) A(n)) a n=-M + F(m) (F(-m-n) A(n)) - 2 2m (F(-m) A(m) - F(m) A(-m)) (2.17) a In order to bring (2.17) into the real domain we note that'(m) F(-m) = 1/4 (G(m) + H(m)2]. Substituting from (2.2), (2.5), and (2.10) into (2.17) we end up with at [G2(m) + H2(m)] = (2.8) M-i 1 -i 1 - Z m ( (B(n)[G(m) H(m-n) - G(m-n) H(m-n)]) - B(n) [ 4a2 n=l o' 2 o 2 + a (C(n)[G(m) G(m-n) + H(m) H(m-n)]) - C(n) - [ I + - (B(n)[G(m) H(m+n) - H(m) G(m+n)]) - B(n)- [ (C 4 4 4 - -(C(n)[G(m) G(m+n) + H(m) H(m+n)]) + C(n) - [ ] aH a6|M 2+ 2 E m(- (B(n)[G(m) H(n-m) + H(m) G(n-m) ) + B(n) [ ] 4a n=m+l *The numbers 1 to 22 in (2.18) and (2.19) refer to the term enclosed by the square brackets. This saves rewriting each term.

20 6 6 + a (C(n)[G(m) G(n-m) - H(m) H(n-m)]) - C(n) [ ] + a (B(n)[G(m) H(n+m) - H(m) G(n+m)]) - B(n) a [ d 8 8 (C(n)[G(m) G(n+m) + H(m) H(n+m)]) + C(n) [ am-1 a I 2 ~G(m) [ + n(-a (G(m)[G(m-n) C(n) + H(m-n) B(n)]) + G( [ ] k2 {-a adJ. 4a n=l 6 10 d 10 + - (H(m)[G(m-n) B(n) - H(m-n) C(n)]) H(m) [ ] 3 11 d6G(m) 11 - - (G(m)[G(m+n) C(n) - H(m+n) B(n)]) + G [ - (H(m)[G(m+n) B(n) + H(m+n) C(n)) + H() [ ] dp~1. adp 1 M 13 __G _ - (G(m)[G(n-m) C(n) - H(n-m) B(n)]) + G(m) 2 () -p 4a n=m+l 4 6H(m) 1[ + - (H(m)[G(n-m) B(n) +H(n-m) C(n)]) - H( -[ (G(m)[G(n+m) C(n)- H(n+m) B(n G() [ 16 166 - (H(m)[G(n+mn) B(n) + H(n+m) C(n)]) + H(m)[ ] m+ --- (B(m)[H(2m) (m) - G(2m) H(m) - B(m) - [L 2a ~ (C G18) 18 - (C(m)[G(m) G(2n) + H(m) H(2m)]) + C(m)) E [ ]]m < M.Ln 6 19 6 19 4- ([ (G(2m)[G(m) C(m) + H(m) B(m)])- G(2m) [ 4a2 ap *(H(2m)[G(m) B(m) - H(m) 2C(m)) - H(2m) 2 0o < + -a( (G(0)[B(m) H(m) - G(m) C(mr)) - G(0) a 21] 2a' dp d

21 + - [B(m) H(m) - C(m) G(m)] a2 In (2.18) the first 40 terms on the r.h.s. represent the changes, due to nonlinear processes, of enstrophy in wave number m. Terms 41 and 42 contain the effect of wave-zonal mode interaction while the last term is due to the p-effect. When (2.18) is integrated with respect to [i the odd terms result in boundary fluxes while the even terms represent the interactions taking place within the region. The first 20 odd terms are denoted by F(m|n,f) while the associated even terms are contained in I(mln,I). Terms 41 and 42 are denoted by F(mlO) and I(mIO), respectively, while the P effect is signified by P(m). The rate of change of zonal enstrophy is expressed by 5(o) G(o0) 1 M 22 = 2- jZ m G(O)[H(m) B(m) - G(m) C(m)] 2a2 Im= 22 cG(O) - [ 1] aGl (2.19) In the calculations M = 15 is taken as the cutoff wave number. This value has been used in the past for reasons of (1) economy, (2) lack ofdata allowing greater field specification, and (3) perhaps most importantly because many studies have shown that energy falls off quickly with wave number. This can be clearly seen in that kinetic energy seems to obey a power law in the region of -3. Enstrophy,

22 however, obeys something like a -1 power law. This indicates that in future studies of vorticity and its exchanges attempts should be made to attain greater resolution. The integration of (2.18) with respect to latitude from c = pS to (p = pN is performed by the trapezoidal rule.* For example, if T denotes the right-hand side of (2.18) except for the P-term then the areal mean contribution may be written as I3(m) = 2i - f T a2 cos dp (2.20) ra (sp -sin cpin S The n-term may be written as CD 1 2Q P(m) = Mg i 2 [ m[B(m) H(m) - C(m) G(m)]coscp dcp psi~np~ S S a ^~(2.21) and by substituting for H(m) and G(m) from (2.11) becomes:(m) = 4 2m B(m) ac= (m) - cosp N P(m) M) B (m) cos a (sincp -sincp ) L B Jm (2.22) and so contributes only to the boundary flux. Analogous calculations were made for the geostrophic potential vorticity, given by Mim a (2.23) =U(0) + U(m) e = V(2.23) m=-M P () (p2 The finite difference scheme for the vertical derivative in (2.23) is shown in Appendix A. *Higher order integration schemes were attempted but the fields proved to be quite smooth and little difference was observed.

23 I(mi n, ) (I (ml n,)) represents the correlation between the gradient of the (potential) vorticity field of wave number mwith the transport of the vorticity field of wave number n by the velocity field of wave number I, where 0 < m.,,n < 15 and the sum of two equals the third. F(mln,e) (F (min,~)) consists of the differences of the correlation of (potential) vorticity fields of wave numbers m, and n with the northward velocity component of wave number I evaluated at cp = 25 ~N and p = 82.50N. I(mlo) (Ip(mlo)) denotes the product of (potential) vorticity transport at wave number m and the meridional gradient of the zonal (potential) vorticity. F(mlo) (F (mlo)) is the boundary flux due to p the correlation of the (potential) vorticity transport at wave number m with the zonal (potential)vorticity. 2.2 DECOMPOSITION OF EXCHANGES AND FLUXES INTO TEMPORAL MODES Many studies (e.g., Bradley and Wiin-Nielsen 1968; Wiin-Nielsen, 1961a, 1961b; and many others) have dealt with some temporal modes of the atmosphere. Without making a detailed frequency breakdown, as done by some (Kao 1968, 1970), three time scales may be identified (Bradley, 1967). These are the stationary or time mean mode, the slow moving or forced mode, and the faster moving, transient or free mode. Any spatial mode is composed of a combination of these temporal modes to varying degrees. Julian et al. (1970) point out that such a separation must be

24 made if comparison of the diagnostic results with the theories are to be made. The turbulence to be described by the power laws must be dominantly of the transient mode since the inertial subrange is defined, in effect, so as not to contain the forced modes. Denoting a mean with respect to a time interval by [ ]t, and the deviation from that mean by ( )t (Reiter, 1969b)wemay write the product of any three quantities F, G, and H in the following manner; [F*G*H] [FG]t[G] t[H]t (l) (2) [F] -[(G)t (H)t]t + [G]t [(H)t(F)t]t (3) + [H] ) [(F)t (G)t]t [(F).(G).(H)] (3.4) The first term on the right-hand side of (3.4) represents the stationary part of the product while the remaining four terms represent the transient part. We can also write [ t = [ ] +([ ]) where s > t and so from (3.4) we get Term (2) = [F] [G]'[H] (4) +[F] ([G]t)s ([H]t) [G]s([H]t).(tF (5) + [H]([F])([G]) s + ([F]t)s ([G]t)s ([H]t (5)

25 Term (4) now represents the stationary mode of the triple product in the time interval s. The remaining four terms on the right-hand side of (3.5) represent the slow transient mode and term (3) in (3.4) is the fast transient mode. Substituting for term (2) from (3.5) into (3.4) and averaging with respect to s yields [F-G.H] = [F] -[G] *[H] + [term 5] + [term 3] (3.6) s s s s s s In this study s = 3 months and t = 5.5 days. The value of 3 months for s was arbitrarily chosen as a compromise between economics and the desire to have a record length which would yield characteristic results and hopefully useful information. Bradley (1967) suggests the value of 5.5 days for t, the averaging time for the running mean filter. This filter reduces all waves with periods of less than 5.5 days to about li/o or less of their original amplitude. Griffith et al. (1956) shows, from power-spectrum analyses of temperatures at University Park, Pa., that maxima exist at periods of 4 and 12 days. The 5.5 day running mean filter, therefore, separates the fluctuations due to cyclonic activity from those of longer time scale. A value of 5 < t < 10 days probably would not alter, qualitatively, the results of this study. The quantities I(mJn,~), I(m|o), F(m|n,e), F(m|o) and:(m), and the corresponding quantities for potential enstrophy, were subjected to time partition as described above. The results are discussed in the

26 next chapter. 2.3 COMPUTATIONAL PROCEDURES The nondivergent flow is obtained from the geostrophic stream function which is defined as the solution of the differential equation V2 =(g/f) V2z - gV Vz (2.28) f where z is the geopotential height, subject to the boundary condition that the wind is geostrophic and the normal component vanishes at the boundary of the octagonal region covered by the National Meteorological Center 1977 point grid. The values of the Fourier coefficients up to M = 15 of this stream function evaluated at latitudes between 20~N and 87.5~N with a 2.5~ interval on each of eight isobaric surfaces, namely, 100-, 85-, 70-, 50-, 30-, 20-, 15-, and 10-cb have been provided by the National Center for Atmospheric Research as the input to the subsequent computations. Differentiations and integrations with respect to both latitude and pressure are replaced by centered finite differences and sums, respectively. Thus, a set of the Fourier components of each of the stream function and vorticity are defined at every latitude from 25 N to 82.5~N with a 2.5~ interval, with the set at latitude ~ representing the mean in the belt between f + 1.25~ and f - 1.25~ on each of the eight isobaric surfaces. The interactions are defined on latitudes 25.0~N to 82.5~N inclusive and the fluxes are evaluated at these

27 limiting latitudes. In the calculations of the interactions and fluxes the values on each isobaric surface are assumed to represent the mean values of the layer which extends halfway to both of the adjacent levels where observations are available, and are weighted accordingly by the thickness of the layer to give the contribution in that layer. Figure 1 shows the vertical resolution of the calculations. Values of static stability are required in the potential enstrophy calculations. The values used are the same as those used by Yang (1967) and given by Gates (1960). The winter values were used for February while the spring values were used for March and April. Table Mlnn 1 lists the values used. Tne static stability is given by aC = a where ia is specific volume, 0 is potential temperature and p is pressure. TABLE 1 STATIC STABILITY, (p) (T m sec2) Layer Layer Winter Spring (cb) 10-15 91.4 95.7 15-20 54.5 42.8 20-30 17.9 13.3 30-50 3.73 3.70 50-70 2.26 2.23 70-85 1.88 1.77 85-100 1.88 1.75

28 10 l 0 I1 20 --- - 20 0 -I4 30 40 40 - -_ - 50 P P (cb)50 (cb) 60 --—. 70Io i6.... 70 80 1- - 85 90 100 100 [ I 8 --'........... —---- 100 (levels of observation) Figure 1. Vertical resolution of the observational calculations.

CHAPTER III DISCUSSION OF RESULTS 3.1 INTRODUCTION The interactions and boundary fluxes derived in section 2.1 are evaluated twice daily, at OOZ and 12Z, for a three-month period extending from February 1, 1963 to April 30, 1963, for each of the 15 components in each layer. These quantities were partitioned into three time modes-stationary, and slow and fast transient modes. Only the values of the three-monthly means of the effects of the various time modes are discussed since the quantities have large daily variations as seen by their standard deviations tabulated in Appendix B. From a preliminary study of the values of I(mln,~) it was decided to form three wave number groups, the large-scale components with 0 < m < 7, the medium scale components with 8 < m < 11 and the small scales with 12 < m < 15. 3.2 THREE-MONTHLY MEANS OF EXCHANGES AND FLUXES The net effect of the various processes discussed below when summed over wave number is a three-monthly mean gain rate for enstrophy -18 -5 -18 -3 (potential enstrophy) of -114 sec (27 x 10 sec ). Figures 3a and 3b contain these values as a, function of wave number in -18 -3 units of 10 sec. Figure 2 illustrates the plotting model for Figures 3a and 3b. 29

50 Gain by wave number m I(m n,) numbr I(mlO) F(ml n, )- Net M(m) - F(mlO) Gain by zonal mode I(OIm) - m — Net F(mlO) Figure 2. Plotting model for Figures 3a and 3b. Arrows indicate direction of positive change. (i) I(mln,l). This term represents the gain of enstrophy by wave number m due to the interactions of the triads composed of wave numbers m, n, and i where n and I attain all possible combinations subject to the restrictions that the sum of any two of Q, m, and n equals the third and 1 < mn, < 15. The sum of I(mln,~) over m equals zero since every term appears twice, with opposite signs. This conservative property serves as a check on the computations and was observed to hold true to the number of digits used. In a closed system with zero mean flow this nonlinear exchange among waves is the only process capable of redistributing enstrophy in wave number space. The results indicate the existence of three wave groups; 1 < m < 7, 8 < m < 11, and 12 _ m < 15 which behave differently. With the exception of m = 4 all of the constituents of the first group lose

31 1 / 1 4-33 — 172 - = 15 36 -- — ~'"6'~"'" —"2'2 16O5 ~2I K 3 62- I-l]_^ 28 3 -C- 38 33 2 — ^ I o 20 1-7 L4 17 2 --- 15 - - 5 i~i4 LII1 10^ e 160' - -- 22 ~!c — _-201-l" 22 --— ~L 1 il s1 z \ ^ ^ 412 0 Figure 5a. Three monthly means of rates of exchanges and fluxes of enstrophy [x 1o-18 sec5] in the atmosphere (0 -00 cb). 10 -1 --— 1 12 13 12-1 - 51 0 01b -- 1 Figure 3a. Three monthly means of rates of exchanges and fluxes of enstrophy [x 10-18 sec-3] in the atmosphere (O - 100 cb).

32 52 10\' - — 52 1 3 2 44 7- 46 253 2 ---- -, 138 m -7 1 8 —'U4-' ~739 24- ( 10 -11 -8 24907 83 1 38-1 - 6-8^-! 6! —4 100 2 3 2 7c —-— 1 —5, —-.2 V 619-I -89102 0 "4* 229 — a-1 17 0.32 —--- 8.,i1193 92.5n cb1). 26 * 2 7J* 4190 - 2 290-. 1711 6 8 "]-~7lo-. 9 -0 —8 - 18-11 1 23 2 17 —=1 -- 7 1 6 12 2 _'" 13 53 -— 3 --- 229 1 2-15 1 3._~-1_ 631_ 3 1 --- 4 - 14 2 — 62 16~1 2 15 1 19-b 831 __ Figure 3b. Three monthly means of rates of exchanges and fluxes of potential enstrophy [x o10-18 sec-3] in the atmosphere (12.5 cb < p < 92.5 cb).

33 enstrophy due to the wave-wave interactions. The loss rate of 51 x -18 -3 10 sec by m = 2 represents the largest contribution of any single -18 -3 wave. The enstrophy loss rate by the group is q = 172 x 10 sec The intermediate group is transitional between the other two groups and may indicate the existence of an inertial subrange. -17 -3 Enstrophy cascades through this wave group with only 10 sec units being absorbed. This represents just under 60 of the rate of enstrophy loss by the group with m < 7. Because of the limited number of waves included in our calculations this transitional group of waves can not be identified as an inertial subrange. The group composed of the smaller scales gains enstrophy at the -18 -3 rate of 161 x 10 8sec. The gain is spread fairly evenly among the components of this group. Due to the closure condition the waves are not linearly independent. In order to examine how local the exchanges are in wave number space the individual terms containing the wave triads would have to be inspected. Although it would not be possible to identify the contributor in any given triad, it would be of interest to calculate, for example, the relative action of waves 3 and 4 vs. waves 20 and 27 on wave number 7. Such a calculation may yield information about the importance of determining the presence of a gap in the wave spectrum. If waves with m = 20, 27 have a significant effect on wave 7 then errors such as those discussed by Lorenz (1969) may affect the large scales even in the absence of an intermediate band of waves.

34 The present results show that the enstrophy cascade is in the predicted sense, that is, from small wave numbers (m < 7) toward the larger wave numbers (m > 7). The existence of the inertial subrange is not verified and more extensive calculations will be required in order to decide the issue. (ii) I (min,). The sense of the potential enstrophy cascade is the same as that of the enstrophy (see Figures 3a and 3b). The magnitudes are larger. Potential enstrophy leaves the first group at the -18 -3 rate of p = 252 x 10 sec. The major source of potential enstrophy for the smaller waves is m = 2, losing at the rate of 65 x -18 -3 10 sec -18 -3 The intermediate group has a gain rate of 23 x 10 sec, which is somewhat over 9% of q. The gain rate of the third group is fairly -18 -3 even with respect to wave number and is 229 x 10 sec. Conclusions arrived at in the light of the above results are: (1) both enstrophy and potential enstrophy cascade from small m toward large m as predicted by Kraichnan (1967), however, (2) there is little verification for the existence of an inertial subrange through which the cascade takes place resulting in little or no gains. (iii) I(m|o). This term represents the gain rate of enstrophy by the waves through their interaction with the zonal flow. It is the correlation of the vorticity advection on a given scale with the zonal vorticity. Comparing (2.18) and (2.19) it is noted that the gain by the waves due to their interactions with the zonal flow is not equal to

35 the loss of the zonal mode due to its interaction with the waves. This is so because the system is not closed, but has boundaries through which there may be advection. This may be written as M 6I(m/O) 6 I(O/m) i/o + (/) 0 but = F(mO) (3.2) at 6t m=l In a closed system, e.g., fluid flow on a sphere, the boundary flux term would disappear and the gain of the waves would be the negative of the gain of the zonal mode due to their mutual interactions. A similar problem concerning energy conversions in an open domain is discussed by Smith (1970). -l8 The net effect of this interaction is a gain rate of 40 x 10 -3 -18 -18 -3 sec. Waves with m = 1,2 lose 8 x 10 and 22 x 10 sec, respectively, while all other waves in group 1 show gains. The middle -18 -3 wave group shows a loss rate of 5 x 10 sec and the short wave -18 -7 group a loss rate of 4 x 10 sec 3 It will be seen in a later section the waves interact with the zonal flow mainly through their slow transient modes. For waves on the scale of the earth's circumference the slow transient and stationary modes are greater than the fast transient mode. For scales with m > 7 the fast transient mode usually dominates. These reasons, in large part, explain the lack of interaction between the zonal flow and waves with m > 6. (iv) I (m|O). The long wave group show a potential enstrophy loss rate of 182 x 10-18 -3 loss rate of 182 x 10 sec while the middle and short wave groups

36 -18 -3 -18 -3 show gain rates of 15 x 10 sec and -4 x 10 sec. Thus the -18 -3 waves have a net loss rate of 171 x 10 sec. The largest loss rate -18 -3 of 101 x 10 sec was by m = 2. (v) F(mln,e) and F (mln,~). The flux terms were evaluated at both northern and southern boundaries separately. The major contribution comes from the southern boundary. Waves with m = 4, 6 and 8 to 12 show losses due to F(mln,~) with the remaining waves experiencing small gains. There is no distinct spectral pattern to this effect. F(mln,2) and F (m|n,~) result in loss p -18 -3 -18 -3 rates of 25 x 10 sec and 24 x 10 sec, respectively, although their spectral distributions are somewhat different. (vi) F(m|O) and F (mlO). The flux due to the interaction bep tween the waves and the zonal mode result in gain rates of 45 x 10 -3 -18 -3 sec and 290 x 10 sec for the enstrophy and potential enstrophy, respectively. The spectral distribution of the gain rates are similar in the two cases with waves 1 < m < 4 accounting for almost all of the change. (vii) P(m) and 3 (m). The rate of change of (potential) entrophy due to the term is (-68 x 18 se -174 x 8 -3 strophy due to the n-term is (-68 x 10 sec ) -174 x 10 sec In both the case of the potential enstrophy and enstrophy the losses are confined mainly to the waves with 1 < m < 7.

37 3.3 RESULTS IN THE TIME DOMAIN In section 2.2 the equations used to decompose the interactions and boundary fluxes in the time domain were derived. The three-monthly means of the temporal modes of the exchanges and fluxes are shown as a function of wave number in Figures 4 through 13. In each of these figures the clear portion of the column represents the contribution of the fast transient mode and the lightly and heavily hashed portions represent the slow-transient and stationary modes, respectively. The solid line indicates the net effect of the three temporal modes on each component. Tables 2A and 2B contain the results in terms of wave groups and we shall for the most part confine our discussion to these two figures. Since, for some exchanges and fluxes the relative effects of the various temporal modes are similar for both enstrophy and potential enstrophy, we shall discuss them in pairs. Figures 4 and 5 contain the temporal decomposition of I(mln,e) and I (mln,e), respectively. For both quantities and m > 12 the fast transient mode accomplishes the major portion of the exchange. The stationary mode is of least importance in the large wave numbers. Wave numbers 2 and 6 experience the largest loss rates in the large scale wave group. In wave number 2 the slow transient mode is the dominant temporal mode while for wave number 6 the fast transient mode dominates. From Tables 2A and 2B it can be seen that, with the exception noted above for m = 2, the exchanges among wave groups take place

TABLE 2A EXCHANGES AND FLUXES FOR THE ENSTROPHY FOR THE WAVE GROUPS IN THE TIME DOMAIN FOR THE ATMOSPHERE [x 10-18 sec-3] Function: I(m|n,~) F(m n,e) I(m 0) F(m O) 5(m) Wave group: L M S L M S L M S L M S L M S Fast ranset -156 10 144 8 -30 13 15 -1 -3 26 -1 1 - 51 -10 -1 transient S slow - 17 7 12 -7 -10 6 45 -5 2 11 4 1 - 82 0 1 transient Stationary 1 -7 6 -7 - 2 5 -12 1 -2 4 -1 -1 -27 -1 -2 Total -172 10 161 -6 -42 24 48 -5 -3 41 2 1 -160 -11 -2 TABLE 2B EXCHANGES AND FLUXES FOR THE POTENTIAL ENSTROPHY FOR THE WAVE GROUPS IN THE TIME DOMAIN FOR THE ATMOSPHERE [x 10-18 sec53] Function: Ip(m| n, ) Fp(ml n, ) Ip(m O) Fp(m 0) 5p(m) Wave group: L M S L M S L M S L M S L M S Fast tras t -244 26 199 13 -27 5 - 62 1 -2 107 7 4 -24 8 0 transient Slow Slow - 27 4 24 -12 - 8 7 -121 9 0 174 0 2 -42 5 1 transient Stationary 2 -8 5 - 4 -4 4 1 5 -3 - 2 - 1 -13 -1 0 Total -249 22 228 - 3 -39 16 -182 15 -4 279 4 7 -79 12 1

Legend Fast transient Slow transient 40 - S tationary FiO 20 - -60 Figure 4. I(mln,e) [x 10-18 sec-3] in the time domain. Legend is used for Figures 4-13. Heavy line indicates net.

604om=1 20 - 010 35 8 c,o "^"1 1 21 134 11 -co0 Poo Hr- 4 T 9 15 ~ -20 - -4o - -6O - -80 Figure 5. I n(mln,~) [x 10-18 sec-3] in the time domain.

41 almost wholly within the fast transient mode. This stems from the fact that the conservation property is applied to each temporal mode alone and the short waves are most active in the fast transient mode. The (potential) enstrophy boundary flux due to the interactions amongwaves, F(mln,~) (F (mln,i)), is depicted in Figures 6 (7). The loss takes place for the most part through the fast transient mode of the middle wave group. The action of short wave group results in gains for all the temporal modes while the long wave loss through the slow transient and stationary modes is partly offset by the gain of the fast transient mode. In the exchanges I(m O) (Figure 8), and I (mjO) (Figure 9), and flux F (mlO) (Figure 11), due to wave-zonal mode interactions the slow transient mode of the long wave group dominates the action. For the flux, F(m|O) (Figure 10), the fast transient mode of the long wave group is dominant. It is noteworthy that, in the main, the slow transient mode is of greatest importance in the interactions between waves and the zonal flow while the fast transient temporal mode dominates in the interactions among waves. The stationary mode is of little or no importance in virtually all the processes evaluated here. However, it is of importance in some processes for specific wave numbers, e.g., m = 3, 4 for F (m 0). The losses both to the enstrophy and potential enstrophy due to the B-effect, P(m) (Figure 12), and Bp(m) (Figure 15), is the result

ho 20 a, ~~~~~~~~~~~~~10 ~~~~~~1 1~~~~~~~~~~~12 m 4 6 8 gV 16 UU:~" co 0?i Ez7yV ^ 7^ 1 -H2 3 -20 -40 Figure 6. F(mln,e) [x 10-18 sec-3] in the time domain.

4020 Figure 7. F(rm|n,A) [x iO18 sec 10] in the time domain. 14 10 -20 -40 Figure 7. Fp(man,~) [x 10-18 sec-3] in the time domain.

-60 -40 -20 H O H 0 4fi 5 6 7 lo l.12 13 14 15 m -20 -40 Figure 8. I(mlo) [x 10-18 sec-3] in the time domain.

20 - m=l 2 613 1.4 ^ V 8 10 7' 9 -20 I OU2~~~~~~~~~~~~~~~~~~~~-62 O -40 -6o -16 -80 - -100 - Figure 9. Ip(mlo) [x 10-18 sec-3] in the time domain.

ho 400 -co 2 -20 Figure 10. F(m o) [x 10-18 sec-3] in the time domain.

120 - 6296 100 1i 80 6)60 -122 O 2 5 Co rIl 40 20 0 3 4 \\\\V//~\\\\~,~l~\\au/I il uuu13 54 1 8 9 ~~~12 M 10 1 -20 Figure 11. F (m o) Ex 10O18 sec-3] in the time domain.

20 m=l 2 3 4 7 8 9 1012 12 13 14 15 E rc-~~~~~~~~~~~~~~~~~~~~~~ ~ (^^-82 -60146 -40 -60-174 -80 Figure 12. P(m) [x 10-18 sec-3] in the time domain.

20 22 i 0'b 1_-20 -J a 1]-66 Figure 13. p (m) [x 10-18 sec-3] in the time domain.

50 for the most part of the action of the slow transient mode of the long wave group. The medium and short wave groups have little effect in any of the temporal modes. 3.4 EXCHANGES AND FLUXES AS A FUNCTION OF PRESSURE The processes, excluding I(mJn,~) and I (mln,e) (which are identically zeros when summed over m) are summed over m and depicted in Figures 14 and 15 as a function of pressure. Figures 16a and 16b show the effect of all the calculated processes as a function of pressure. The vertical dashed lines indicate the pressure-weighted mean values. Tt is seen in Figures 14 and 15 that the intensity maxima for the exchanges and fluxes for both the enstrophy and potential enstrophy occur in the upper troposphere. Yang (1967) shows similar results in the exchanges of kinetic and potential energies due to interactions among waves. These results reflect the presence of maximum amplitudes for the eddy velocity at the level of the Jet Stream. Both enstrophy and potential enstrophy have maxima in the 20 cb to 40 cb layers as shown in Figure 16c. The evaluated processes result in a negative rate of change of enstrophy in all layers of the atmosphere except for 92.5 < p < 100 cb which experiences a small gain. The maximum rate of change of -261 x 10 sec occurs in the layer centered at 20 cb. The pressureweighted mean loss rate is -144 x 10 sec

51 (a) 20 2 - 040 o Z 60 80 100 — -200 0 +200 0- ( ) 20 c, 40o f2 60 80 100 -,, -200 0 +200 ~~~~~~~0 ~(c) 20 460 80o 100, -200 0 +200 of 0 Im o [ 13 (d V t a v i of(d ) 20 40 C) 60 8o 100 -200 +200 Figure 14. (a) Vertical variation of F(mn,e) ['x 10-18 sec-35]. (b) Vertical variation of I(mlo) [x 1018 sec-5]. (c) Vertical variation of F(mlo) [x 10-18 sec35]. (d) Vertical variation of ~(m) [x 10-18 sec-5].

52 (a) 0 pZ 60 100 l -400 -200 0 +200 o0, (b) 20 2 40 C) A 60o 80 100.,' -400 -200 0 +200 0 (C) 20 40 m 0 60 "C. 80 100 - -200 0 +200 +400 0 (d) 20 - ZF 6o 80 100.,, -400 -200 0 +200 Figure 15. (a) Vertical variation of Fp(m|n,) [x 10-1 sec53]. (b) Vertical variation of Ip(m|o) [x 10-18 sec~5]. (c) Vertical variation of Fp(m|o) [x 10-18 sec-3]. (d) Vertical variation of,(m) [x 10-18 sec35].

53 T0 4 |(a) 20 L 00'L 60 80 100., -200 0 +200 ~~~~~~0 (b) 20 80 100 -200 0 +200 Figure 16. (a) Vertigal variation of sum of all processes affecting the enstrophy [x 10" 1sec5-]. (b) Vertical variation of sum of all processes affecting the potential enstropny [x 10 1 sec-3].

54 20 20 /// 4-0 / ///I 6o // i// I 80 / 100 2 4 6 8 10 12 14 Figure 16. (c) Vertical variation of the three monthly means of enstrophy and potential enstrophy [x 10-18 sec53].

55 The potential enstrophy on the other hand experiences losses in the upper and lower troposphere but gains in the middle troposphere. The gains are effected by the action of F (mlo) and (p(m). The p p -18 -3 pressure-weighted mean effect is a gain rate of 29 x 10 sec 3.5 ENSTROPHY AND POTENTIAL ENSTROPHY BUDGETS Our calculations in no way account for the gain or loss of (potential) enstrophy due to production, dissipation, or exchanges involving waves with m > 15. We may, however, make some comments on these effects by calculating the residual changes in (potential) enstrophy allowing for the processes evaluated in this study. We write the equation for (potential) enstrophy change as follows -t (m) = I(mln,i) +I(o|m) +F(mln,i) +F(mlo) +p(m) +residue (3.7) The l.h.s. of (3.7) was evaluated as follows [t(m)] - [(m)] (m) =3.8) = at where r and q represent the first and last ten time periods in the three months. The residues thus calculated are listed as a function of wave number in Tables 3A and 3B. A possible explanation for the residues obtained is as follows. Positive values for m = 1, 2 may be due to vorticity production at the scale of continents and oceans due to differential heating at coast lines. For 6 < m < 9 production of

TABLE 5A ENSTROPHY (CORRECTED*), THE TIME MEAN CHANGE AND RESIDUE Sum of c At/At _ Evaluated Residue, m -12 -2 -14 -2 -l -18 Processes,18 - m~~~~- ---- -- -----— * —----- x 10 sec _ x 10~ sec day _ X 10 sec x 1018 sec-3 x 10 sec 1 50 - 6 -.7 -55 52 2 54 - 68 - 7.8 -102 94 5 51 - 18 - 2.1 - 1 - 1 4 48.4.1 - 15 15 5 48 - 17 - 2.0 - 6 4 6 49 - 24 -2.7 - 74 71 7 46 15 1.7 - 20 22 8 46 9 1.0 - 17 18 9 59 -10 - 1.1 - 21 20 10 41 2.2 8 -8 11 52 - 10 - 1.1 - 16 15 12 29 - 6 -.6 26 -27 15 28 - 14 - 1.7 47 -49 14 24 - 6 -.7 49 -50 15 24 - 10 - 1.2 59 -60 Total 607 -162 -18.8 -114 94 *Correction factor is discussed in next section.

TABLE 3B POTENTIAL ENSTROPHY (CORRECTED), THE TIME MEAN CHANGE AND RESIDUE Sum of _pc, A_^p/At Evaluated Residue, m -12 -2 -14 -2 1 18 - Processes, -18 -3 x 10 sec x 10 sec day 1 x 10 secs x 10 c 1 sec 1 66 - 13 - 1.4 - 49 48 2 71 - 86 - 7.9 -129 119 3 66 - 19 - 2.2 18 -20 4 61 1.1 7 - 7 5 61 - 19 - 2.2 10 -12 6 61 - 29 -.4 - 89 86 7 56 26 3.0 - 7 10 8 57 16 1.9 4 - 2 9 46 - 7 -.8- 1 0 10 49 6.7 21 - 20 11 38 - 7 -.8 - 9 8 12 34 - 4 -.4 41 - 41 13 32 - 13 - 1.5 63 -65 14 29 5 -.5 62 - 63 15 30 - 9 - 1.1 83 -84 Total 755 -162 -18.5 27 - 43

58 vorticity associated with cyclonic generation due to differential temperature advection near frontal zones may account for the residue. Negative residues for 12 < m < 15 may represent dissipation and cascade to waves with m > 15. 3.6 SPECTRAL SHAPE It has been shown (Horn and Bryson, 1963; Wiin-Nielsen, 1967; Julian et al., 1970) that the atmospheric kinetic energy may be written in the form K(m) = b m 7 < m < 15 (3.9) where b and c are constants. b is a measure of intensity while c determines the spectral shape. The value of c has been found to vary between 2.5 and 3 with a tendency toward the latter value, for the spectral region given by 8 < m < 16. The variation with latitude was found to be small while a systematic increase with decreasing pressure, up to about 20 cb was also observed. 2 In Appendix C we show that g(m ) = m K(m ) and therefore c 1 e e should result from a calculation using (3.9) in which 5(m) replaces K(m). The spectra of enstrophy and potential enstrophy were calculated. Their three-monthly mean values were evaluated as a function of wave number and pressure and a least-square regression was used to calculate

59 the values of c. For these calculations the spectrum was divided into two spectral groups. The first contained wave numbers 1 to 7 and the other wave numbers 8 to 15. The results are tabulated as a function of pressure in Tables 4A and 4B. The pressure-weighted mean spectra are contained in Figures 17a and 17b. Also included in Tables 4A and 4B are the values of c after a correction is made to the enstrophy (potential enstrophy) in order to counter the effect of smoothing on the short waves. Using WiinNielsen's (1967) argument, it is found that 2 -2 c(m) = 1- - (m) (5.10) m where 4 (m) is the corrected or "true" enstrophy, c m. is an arbitrary cut-off wave number taken equal to 36 and S(m) is the calculated enstrophy per unit wave number. For spectra of c and p,.74 < c < 1.20 and.69 < c < 1.29, respectively. The slope is maximum at 50 cb for both quantities, indicating relatively greater amounts of enstrophy and potential enstrophy in m = 8, 9 than in m = 14, 15 at this level. However, maxima in all scales are found at this level. The values of c imply kinetic energy spectrum governed by a -d power law where 2.7 < d < 3.2. These results are similar to those of Wiin-Nielsen (1967). The pressure-weighted mean values of c are 1.085 for enstrophy and 1.145 for potential enstrophy. A slope of -5.085 would, therefore, govern the corresponding kinetic energy spectrum.

60 TABLE 4A c FOR [] s [C]scor., [ p] AND [c]or. (Wave numbers 1 to 7) P-Weighted P(cb) 100 85 70 50 30 20 15 10 Mean.0.1.1.1.1.0.1.1.1 cor..0.0.1.1.0.0.0.1 ~p -.-.1.2.2.1.0.1 -.1 pcor. --.1.1.2.1.0.0 --.1 Note: s = 3 months. TABLE 4B C FOR [f], [ c] cor., [ p] AND [pt cor. (Wave numbers 8 to 15) P(cb) 100 85 70 50 30 20 15 10 Pe e Mean e5 1.2 1.3 1.6 1.6 1.6 1.6 1.4 1.3 1.5 t cor..7.9 1.2 1.1 1.2 1.1 1.0.8 1.1 p --- 1.2 1.5 1.6 1.7 1.6 1.4 --- 1.6 | cor. ---.7 1.1 1.2 1.3 1.2.9 --- 1.1

61 8 c =-1.084 CV M c =.048 o 0 I 4 ^ 5. \. 3~~~~~~~~c= 1.084 I I I I I I I I I I 2 4 6 8 10 12 14 16 m - wave number 7- c =.089 OI 1.145 0-C P Ii o 5 0 H 0 O H PL4 I I I I I I I I i I 1 2 4 6 8 10 12 1416 1 m - wave number Figure 17. (a) Spectrum of [S],s' (b) Spectrum of [ pc ]s

62 The slopes for 1 < m < 7 for both the enstrophy and potential enstrophy are about -.1 (Table 3A) indicating that these quantities have no trend in this wave number range. The maximum for both and p occur at m = 2, 3. This is due to the influence of the stationary mode since these maxima shift toward larger wave numbers when the stationary mode is excluded. Julian et al. (1970) suggest that only the transient modes should be subjected to the power law analysis. This was also done. From Tables 2A and 2B we can see that the stationary mode contributes little to the cascade of | and t in wave number space. The theories, if valid at all, should be expected to apply only to the transient modes. The transient enstrophy is given by () =,1 (G)2 + (H)) s = 3 months (3.11) The slope was calculated using the three-monthly mean values after they were subjected to the correction indicated by (3.10). We will here confine our discussion to the results as a function of pressure. The enstrophy of the transient mode for waves with 8 < m < 15 exhibits slopes ranging from -.8 at the 100 cb and 10 cb surfaces to -1.2 at the 30 cb surface. Tables 5A and 5B contain the slopes for the corrected enstrophy and potential enstrophy. Figures 18a and 18b show the mean values of these quantities as a function of wave number. In the wave group with 1 K m < 7 the slopes at 30 cb, 20 cb, and 15 cb are positive showing that a maximum of the quantity is obtained

65 TABLE 5A c FOR THE TRANSIENT ENSTROPHY AND POTENTIAL ENSTROPHY FOR THE SPECTRAL REGION 1 < m < 7 (x 10-1) P-Weighted P(cb) 100 85 70 50 30 20 15 10 Mean -.6 -.2..2 -.5 -.7 -..4 -.2 ---.2.6.6 -.3 -.8 -.4 --.0 p TABLE 5B c FOR WAVE NUMBERS 8 TO 15 FOR TRANSIENT ENSTROPHY AND POTENTIAL ENSTROPHY (CORRECTED) P-Weighted P(cb) 100 85 70 5 0 20 15 10 -Mea Mean.8.9 1.0 1.I 1. 1.1.0.8 1.1 --.6 1.0 1.0 1.1 1.1.8 -- 1.0.p

64 o2 I - C c = -.019:0 4 -.r) —;~ —! 0 0 s'' 2 c = 1.067 I,I I I I I I IO I 2 4 6 8 10 12 14 16 m wave number 8 0 _ c =.002 k 7 2 4 *r -1 4 0z 0 I c =*997 *r - - 2 1 2 4 6 8 10 12 14 16 8 20 m wavenumber Figure 18. (a) Spectrum of the three monthly mean transient enstrophy. (b) Spectrum of the three monthly mean transient potential enstrophy.

65 for m - 7 when the stationary component is omitted. In mid-latitudes on the 20 cb surface, slopes of about.25 are obtained. A value of.33 would correspond to a "-5/3 law" for the kinetic energy spectrum. 3.7 COEFFICIENT OF EDDY VISCOSITY AND RELATED QUANTITIES Leith (1968) deduced the local enstrophy dissipation r to be given by 2 = 7vl| 2 (3512) where v is the coefficient of eddy viscosity. From this equation we may evaluate v, the eddy viscosity coefficient. Using M imk = Z F(m) e -M f0 we find that M1 M M V^ = (+ - a + A 1 a) ~ F(m) e 1, aF(m) im +A im ilf1 e + k, - - _ F(m) e (3.13) a a~p cos -M -M _0 70 and so [(v0 i ) ] = [v72]^ " F(-m) M 2 2 F(m) t F(-m) + C m F(m) F(-m (.14) a2 _ 6 2

66 where 4 and X are unit vectors in the 0 and X directions, respectively. Using the definitions of G(m) and H(m), we find 2 2FM 2 ()2 (m)2 2 (O2 2 - l1/4 + - + (G(m)2 + H(m)2) a2 / cos 42 (3. 15) Here, T is defined as -Z I(mln,~) and is shown in Table 6C as a funcm=l tion of pressure. For the 100, 85, and 10 cb surfaces r is either very small or slightly negative and so these surfaces were omitted from this calculation. This indicates that any similarity theory dependent upon r1 which may be found to hold in the middle and upper troposphere should not be expected to pertain to lower layers where ground effects may be present nor in the stratosphere. 2 2 First (V7)2 and (V7 ) were calculated as a function of latitude and pressure for 15 cb < p < 70 cb and their time means obtained. These values are contained in Tables 6A and 6B. It is seen that the values increase to a maximum at the 30 cb surface. Values of [v] were obtained from s N ( (p) V]sX (p) = (3.16) sA' It must be pointed out that in our truncated system the values of 2 both (V:) and q are underestimated. The net effect on v due to these shortcomings cannot be evaluated a priori. Kraichnan (1967) states "...-3 power law can be expected to extend up to the wave number of the

67 TABLE 6A (V) 2 —THE SQUARE OF THE VORTICITY GRADIENT [x 10-25 cm-2 sec-2] P(Cb) 7~ 30 3~ 20 13 P-Weighted 70 50 30 20 15 cp( ON) _____________Mean 25 2.6 4.6 8.0 8.6 7.8 5.6 27.5 2.9 5.2 9.7 10.8 9.3 6.6 30 3.1 6.2 11.6 13.5 11.1 7.8 32.5 3.3 6.8 12.4 15.0 11.8 8.5 55 3.6 7.6 14.1 15.4 11.7 9.2 37.5 3.9 8.6 16.0 15.3 11.4 10.0 40 4.4 9.1 17.5 15.8 11.5 10.8 42.5 4.9 10.3 18.6 15.7 10.6 11.4 45 5.2 10.9 18.9 14.0 9.0 11.4 47.5 5.4 10.9 18.2 12.8 7.9 11.1 50 5.7 11.4 18.1 11.8 7.0 11.1 52.5 6.0 12.0 18.8 10.9 6.1 11.4 55 6.0 12.2 18.6 10.4 5.7 11.3 57.5 5.9 12.5 18.9 9.8 5.3 11.3 60 5.7 12.3 18.8 9.3 5.0 11.1 62.5 5.7 12.2 18.1 9.2 5.2 10.9 65 5.7 12.2 17.5 9.2 5.3 10.8 67.5 5.6 11.8 16.1 9.2 5.4 10.3 70 5.6 11.5 14.5 8.9 5.2 9.8 72.5 5.5 11.2 12.6 8.2 5.0 9.2 75 5.4 11.1 11.3 7.9 5.2 8.8 77.5 5.4 11.0 11.0 7.8 5.3 8.7 80 5.8 11.6 11.6 8.5 5.9 9.2 82.5 9.2 19.5 17.9 14.9 10.5 15.1 Area 4.7 9.6 15.3 12.0 8.4 9.6 Mean

68 TABLE 6B (V: )2 -THE SQUARE OFTHE POTENTIAL VORTICITY GRADIENT [x 1025 cm2 sec ] P(cb) 20 15 P-Weighted 70 T 0 30 20 15 cp(0N) 70 _ _ Mean 25 3.2 5.3 8.8 9.2 8.2 6.2 27.5 3.5 6.0 10.7 11.6 9.8 7.4 30 3.7 7.2 12.8 14.6 11.7 8.7 32.5 3.9 7.8 13.8 16.2 12.4 9.4 355 4.2 8.7 15.6 16.5 12.3 10.3 37.5 4.4 9.8 17.6 16.3 11.8 11.1 40 4.9 10.4 19.3 16.8 11.8 11.8 42.5 5.2 11.6 20.7 16.7 10.8 12.5 45 5.6 12.3 21.1 14.9 9.1 12.6 47.5 5.9 12.2 20.3 13.6 8.0 12.2 50 6.2 12.8 20.2 12.5 7.1 12.2 52.5 6.5 13.3 20.9 11.6 6.1 12.5 55 6.5 13.6 20.7 11.0 5.8 12.4 57.5 6.4 14.0 21.0 10.2 5.3 12.4 60 6.2 13.8 20.9 9.6 5.1 12.3 62.5 6.1 13.7 20.0 9.5 5.3 12.0 65 6.2 13.8 19.3 9.6 5.5 11.9 67.5 6.1 13.4 17.7 9.6 5.6 11.4 70 6.1 13.1 15.9 9.3 5.4 10.8 72.5 6.0 12.7 13.8 8.6 5.1 10.1 75 6.0 12.7 12.3 8.2 5.4 9.7 77.5 5.9 12.6 11.8 8.2 5.5 9.6 80 6.3 13.1 12.4 8.8 6.1 10.1 82.5 9.7 21.0 18.6 15.4 10.7 16.0 Area 5.2 10.8 16.9 12.8 8.7 10.8 Mean TABLE 6C (POTENTIAL) ENSTROPHY CASCADE [x 101 sec ] P-Weighted P(cb) 70 50 30 20 15 Mean ( sec-3) 49 227 540 310 148 172,p(sec 3) 88 298 579 460 136 252

69 TABLE 7 COEFFICIENT OF EDDY VISCOSITY [x 108 cm2 sec1]'oP(cb) =~ n~ 30 20 15 P-Weighted P(cb) 50 0 20 1 Mean rMean. v 1.1 2.6 3.8 3.6 1.9 2.5 v 2.2 3.8 4.0 4.1 1.7 3.3 order of k - k = (rQ/v) 1/6 Using the values of q and v calculated -2 -1 herein it is found that k 1.5 x 10 km. The corresponding wavelength is Xd 70 km. This result is dependent on the assumption that the inertial cascade process is dominant over other effects for m > m. where m. is the wave number of the energy source. It is most probable that there are other scales between wave numbers m = 8 and m = 400 at which kinetic energy is produced in the atmosphere (Pinus et al., 1967). Therefore, it should not be expected that the real atmosphere will exhibit a "-3" kinetic energy spectrum down to wave lengths of about 100 km. It is of interest to note, however, that d is of the order of the largest convective systems in which the vertical dimension can certainly no longer be ignored. The corresponding viscous dissipation of kinetic energy is given, 2 -1 by Kraichnan, to be Kd ~ (k ) which is calculated to be of the d d -2 2 -5 -2 -1 order of 10 cm sec or 10 erg cm sec. This is negligible when compared with the observational estimates of dissipation made by various authors and reported on by Oort (1964). They estimate the

70 5 -2 -1 dissipation to be of the order of 10 erg cm sec. Kung (1967), using data covering a five-year period, obtains a time mean dissipation 3 -2 -1 of kinetic energy of about 4 x 10 erg cm sect The dissipation of kinetic energy, therefore, does not take place at relatively small scales in large scale atmospheric flow as it does in microscale threedimensional turbulence. Manabe et al. (1970) show that in a numerical simulation of the large scale atmospheric flow maximum dissipation occurs at m = 2, 3 and is negligible for m > 14. h4 -2 -1 to be 8, K. is calculated to be 10 erg cm sec which is larger than estimates of C(A,K ) (Wiin-Nielsen, 1965; Oort, 1964) by a e e factor of about 4. C(A,Ke) is the transformation of eddy available potential energy to eddy kinetic energy. Using the formulation of Leith and Kraichnan for the quantities evaluated in this section we find (a) that the dissipation to be maximum at 30 cb and in midlatitudes. The variation with respect to pressure is in agreement with the calculations of Kung (1967). The eddy viscosity coefficient is 8 2 -1 evaluated to be about 3.7 x 10 cm sec in the layer of maximum dissipation; (b) that in the absence of intermediate scales of energy production the -5 power-law can be expected ideally to extend to scales on

71 the order of 70 km, which, although reasonable as a limit should not be expected to hold in the real atmosphere; (c) that the energy dissipation at kd corresponding to r is on the -2 -1 order of 10 erg cm sec and is negligible when compared to dissipation at the synoptic scales; and (d) that the equivalent energy input Ki, if assumed to be acting at m. = 8 is about a factor of 4 larger than the various estimates of C(Ae,K ) reported on by Wiin-Nielsen (1965) and Oort (1964). 3.8 TRANSFER COEFFICIENT AND POTENTIAL VORTICITY TRANSPORT It is well known that momentum transfer does not obey a diffusion law of the form 3[u], [uv]3 = -K() a (3.17) In the region of large transports such a formulation would result in negative values of K. It was suggested by Green (1970) that only conservative properties of the fluid should be expected to obey this classical diffusion law. In the atmosphere there are no truly conservative quantities. However, quantities such as potential vorticity and potential temperature are considered to be quasi-conservative. It is already known that the meridional transfer of sensible heat can be approximated in a qualitatively correct manner by a large scale diffusion coefficient in most regions.

72 It can be shown that the momentum transport may be written in terms of the potential vorticity and potential temperature transport as follows _ 1! M cos2 o /f R N /,o [vC = a_ Co2 6 6a- p N) a(3.18) a cos < (a) (b) where N = [vT]k and M = [uv]. Therefore, if it can be shown that the potential vorticity transport obeys a diffusion equation then the effects of the eddies in transporting momentum across lines of latitudes may be parameterized through the transfer coefficients for sensible heat and potential vorticity. The potential vorticity transport was calculated using (2.1) and (2.7) as follows M [v]x = 2 cos m[C(m) V(m) - B(m) W(m)] (5.19) im=l The corresponding transfer coefficient may then be evaluated by [v p] K = 6( p(0) + f) a; a5 where M (O) = U(O), as defined by Eq. (2.23) Daily values of the meridional potential vorticity transport and K the transfer coefficient, were calculated as a function of 4 and p. Calculations of the potential vorticity transport using the monthly and annual means of N and M in (3.18) were carried out by

73 Wiin-Nielsen and Sela (1971). The transport was found to be negative everywhere except for the 850 mb surface and a small area around 40~N at 20 cb. The magnitude of the transport was found to be ~.8 x 10 -2 cm sec for the annual mean. In our calculations we find a band extending from 35~N to 55~N in which positive values of the potential vorticity transport are found at all pressure surfaces. Outside of this band negative values are obtained everywhere. Figure 19 shows the distribution of [v: ] as a function of p X,,s and p. Maximum values, both positive and negative, are found at the upper layers of the troposphere (40-17.5 cb) a typical value in the -2 -2 -2 mid-latitude land -.4 x 10 cm sec while values around -.25 x 10 -2 cm sec occur in the northern and subtropical regions. K, was calculated using (3.20). The denominator, the latitudinal tp gradient of absolute potential vorticity, is generally expected to be positive everywhere due to the dominance of the P-term. This is so in the time average. However, there are a sufficient number of daily occurrences of negative values to yield values of Kp the same sign as Up [vp ]. The influence of these cases with negative values are exaggerated in their effect on [K ] because they are usually small. This Up s observation was made by a cursory examination of some 40 time periods. 11 2 -1 Values of Kp are of the order of magnitude of 10 cm sec These values may be compared with those obtained by Wiin-Nielsen and Sela (197l) using the annual mean values of [vLp] and c^(o)/aa4. In

20 20 \'.: \\ / x3 I o / \ \ I /' " -30) 3J -30,,^^ / / 0, \ - /' I % 0 / 11 1 / II 60 - / | I I / / I ( \ -'\ _ — - -'~",/^ ^ - - - - -~ —z~ --,1~~ o-\ / - 60 80 75 70 65 60 5 50 45 40 55 50 25 4(~N) latitude Figure 19. [v p 3xs [x 1011 cm sec-1] mean of February, March, and April 1965.

75 10 2 his study Kp is positive everywhere and of the order of 10 cm -1 sec It is to be expected that calculations based on daily values would result in larger variations with respect to latitude and pressure than calculations based on annual mean values. The calculations herein result in a mean value of Kp, while that in Wiin-Nielsen and Sela result in a value of Kp for the time mean of p. Oort (1964) disC~p P~~p cusses the differences encountered upon comparing results from the mixed time-space domain. Terms a and b in (3.18) were calculated separately to evaluate the relative contributions to the transport. Term b was negative almost everywhere except on the 85 cb surface and between 42.5~N and 62.5~N on the 20 cb surfaces. Term a dominated the result. However, it has been noted (Holopainen, 1967) that this data yields momentum transports too large by about a factor of 2. A correction of this magnitude might result in term b dominating the result and yielding potential vorticity transports of negative sign in a greater number of locations. Whereas Wiin-Nielsen and Sela (1971) have shown that over a long time period the potential vorticity transport can be approximated by (3.20) with positive Kp the present results indicate this not to be the case in the short term. Use of Green's (1970) suggestion and (3.18) to parameterize the eddy momentum transfer should be limited to the flow on an annual or at best monthly time scale.

CHAPTER IV THE NUMERICAL MODEL 4.1 INTRODUCTION Since the attempt by Phillips (1956) to simulate the atmospheric circulation by integration of a simplified set of the Navier-Stokes equations there have been a host of general circulation models of varying complexity constructed (Kasahara and Washington, 1967; Smagorinsky, 1963; Manabe et al., 1970 and others). Some studies are limited to particular kinds of fluid flow; specifically, Ogura (1962a, 1962b, Lilly (1969), and Orzag, (1969) have investigated the evolution of turbulent large scale flows. Ogura performed his computations in wave number space while Lilly's are in the real domain. Orzag discusses the benefits and drawbacks of the two methods. Lilly constructed a two-dimensional model of turbulent flow and was able to approximately verify the "minus five-thirds" and "minus three" spectral ranges. His model was described on one of an infinite number of contiguous identical squares. Thus, the lateral boundary conditions are that the flow is cyclic in both directions. Lilly's initial condition is a spike in the kinetic energy spectrum at m. = 8 (in one run m. =4) with noise several orders of magnitude below the spike for the other spectral components. 76

77 In a study by Bray, reported on by Batchelor (1969) the integration commenced from a variety of spectral distributions. An intermediate spectral region between the wave numbers which initially contained the enstrophy and those at which dissipation was important was found to adhere fairly closely to a "minus one" power law. This study attempts to extend the previous studies by using a combination of Lilly's forcing functions and boundary conditions in the framework of Phillips' (1956) quasi-geostrophic model. The disturbances will be isotropic and homogeneous, and cyclic in the x- and ydirections with zero mean. In the absence of generation and dissipation of energy and enstrophy there exist two conservative quantities, namely, the potential enstrophy and the sum of the kinetic and available potential energies. (These quantities will be defined below.) Phillips (1956) has shown that the two-level quasi-geostrophic model exhibits many of the characteristics of the atmosphere. The conversion of energy from zonal available potential energy, AZ, to Ae, the eddy mode, C(AZ,A ) the conversion from A to K, C(A,K ) e Zi e e e e e and C(K,Kz), the conversion from the eddy to zonal kinetic energy are eZ all qualitatively well represented. K and K denote the eddy and e Z zonal forms of kinetic energy. The model used in this study is a simplification of Phillips' in that we omit the zonal modes of energy. The primary concern here is the interactions among waves. It is known that many of the assumptions (isotropy, homogeneity)

78 entering into two-dimensional turbulence theory (Kraichnan, 1967) are not applicable to the atmosphere. The necessary condition for the existance of an inertial subrange, that the wave numbers at which energy is put into the spectrum is well removed from the range in which dissipation becomes important, is not observed in general circulation models (Manabe et al.,1970). Yet, in the atmosphere (Wiin-Nielsen, 1967; Julian, et al., 1970) as well as in general circulation studies (Manabe, et al., 1970) the spectral slope of the kinetic energy is very close to -3 on a log-log scale. With the aid of a simple numerical model we shall show that a -3 slope is not unique to the set of equations representing the atmosphere but that other quasi-stationary slopes may be obtained by the manipulation of the intensity of the circulation through the amplitude of the heating field and the eddy viscosity coefficient. 4,2 FORMULATION OF THE MODEL The equations to be used for the model are taken, except for the omission of the surface friction effect, from Phillips (1956) and are 0u uau vau -'p PV1 ~t + a~ + - fv =: + v ( 4.1) 5t dx dy o o x / av u~v v~v 2 at + ax y +fu o - + vu (4.2) Cu av aOD 1+ + p ~ (4.3) 1 dQ d cpT dt dt= ne

79 p = pRT (4.5a) and = (h.5b where t is the geopotential height c is the vertical motion when pressures is used as vertical codp ordinate (= d) dt dQ is the nonadiabatic heating rate per unit mass. dt 9 is the potential temperature (~K) p is the density T is the temperature (~K) c is the specific heat of air at constant pressure, and p f = 2 Q sin (45~)(=10 sec ) andQ is the earth angular velocity. 0 Equations (4.1) and (4.2) are to be used at levels 1 and 3 (see Figure 20). The subscripts will indicate the level at which the quantity is evaluated. 0 P = ocb __ =O 0 ~1 25 q = -VVq + vV q + H 11 251ilv 1 2___ _ at 3 -V3 Vq3 + vV2q3 - H 2 50_____ fo~~2= k2 (t+ V5 f V)T + H/2- vv2T 3 75 =..\ Vq5+vy2q5 H 4 100__________4 = o Figure 20. Schematic of model's vertical resolution, equations, and vertical boundary conditions.

8o The vorticity equations on surfaces 1 and 5 are obtained by cross differentiation of (4.1) and (4.2) and then forming the difference (4.2) - (4.1). The geostrophic assumption is introduced for all the velocities in the resulting equations except in terms containing V * V, the horizontal divergence. For this term we substitute from (4.3), the continuity equation. Using the boundary conditions of c = 0 at p = 0 and 100 cb, and evaluating -, at p = 25 cb and p = 75 cb by centered finite differences we get v * = -V = -%2/P2 (4.6) The vorticity equations are (t V1 V) = p + vV 3 (4.7a) 2 v) = +V (4.7b) ( dt 3 3 P2 3 The Coriolis parameter, fo, is taken to be a constant even when differentiated for two reasons. Firstly, we do not wish to disturb the isotropy of the model and so we must treat the two lateral directions alike. Secondly, if (4.7a) and (4.7b) are integrated over the domain df of our model and the approximation of - = a, a constant, is made then dy the P-effect in the integral sense is zero in the light of our boundary conditions and the nondivergent nature of the flow. Since we are interested in integral quantities the retention of the M-terms in (4.7a) and (4.7b) are unnecessary.

81 The thermodynamic equation (4.4) with the use of the equation of state, the hydrostatic equation and a geostrophic stream function defined by ft = 0 (4.8) becomes o 2 2( R dQ (4.9) op'2=.[(- + 2 V)T - c f dt] 2 D o where T =' - V 3, a measure of the thermal field. (Note change in definition of T) and V, =(V1 + V) Finite differences have replaced the vertical derivatives in (4.9) and fG 2 [0 - 0 )( - -)], a constant o2 3 1 3 1 Only two of the nonadiabatic heating effects are included in the last term in (4.9). These are a heating term, which provides a means of putting energy into the system, and a representation of lateral eddy diffusion. We may then write dQ dQh dQd dt dt dt where

82 dQ f c H h - o p dt 2 dQd c f _d p= v2T dt R and R is the universal gas constant. Both of these effects, while being defined at P2 = 50 cb may be thought of as being the average value for an entire column of the atmosphere. Phillips (1956) formulated the dissipative mechanism to include the effects of horizontal diffusion and surface drag. Charney (1959) contends that the vertical diffusive effect should be included while the effects of the horizontal diffusion may be neglected. Smagorinsky (1963) includes all three effects in his atmospheric simulation and finds that the horizontal diffusive effects on the eddy kinetic energy are of the same order of magnitude as the effect of vertical diffusion combined with the surface friction. Manabe et al. (1970) confirm this finding and also show the vertical and horizontal effects to have similar spectral distributions especially for m > 8. Furthermore, the vertical diffusion effects are shown to be limited to the atmospheric layer beneath 85 cb. In this study we do not wish to include surface effects since we certainly cannot expect the two-dimensional turbulence theory to be applicable in that region.

853 We shall also omit the vertical diffusion effect. It may be thought of as being included in the horizontal diffusive effect by an adjustment of the eddy viscosity coefficient. The characteristic of the dissipative mechanism which must be retained is its spectral distribution and this is done. The magnitude of the eddy viscosity coefficient is a parameter in this study which will be varied. Eliminating ao between (4.6) and (4.9), and between (4.7) and (4.9) yields the prognostic equations of the model which are V Vq + vVq + H W "v1 v% + vv0q1 + - -v Vq +v2q3- H (4.10) where ql = V2 - T 2 2 q3 = V*23 + xT (4.11) H, the heating term in Equation (4.10), acts as the energy input to the system. If the dominant wave number of H is taken to be mi then all the two-dimensional Fourier components having max(mym ) = m will be included. The heating, therefore, will not be truly monochromatic. The heating field at time t + 1-is related to the previous heating field by the following equation Ht+l R Ht + [ - R2]1/2 t+l

84 R is a correlation coefficient such that 0 < R < 1 and is a measure of the continuity of the heating field with time. H is the varying part of H and the amplitudes of its spectral components are chosen randomly from a Gaussian distribution. The total amplitude of H is kept constant by normalization and so only the phase relationships change. Since the thermal field at any time is a function of the heating field (as well as other variables) a value of R close to zero implies little correlation between the temperature field and the newly composed heating field. Inhibition of the generation of available potential energy results. In the experiments conducted in this study m. = 8. This corresponds to the spectral region of maximum baroclinic instability as given by theoretical considerations (Derome and Wiin-Nielsen, 1966). The observation studies (Saltzman and Fleisher, 1960, 1961; Wiin-Nielsen, 1959) of the conversion of eddy available potential energy to eddy kinetic energy indicate a maximum at m - 6,7. By using m. = 8 the conversion process is forced to have its maximum at m = 8, thereby maintaining correspondence with the atmosphere. As already mentioned the region of concern will be one of an infinite number of identical contiguous volumes containing five pressure surfaces. Figure 20 contains a schematic of the model's vertical resolution.

85 The lateral boundary conditions are that the flow is cyclic in both the x and y directions. Figure 21 illustrates these conditions on a 5 x 5 grid. (5;5) (5;1) (5;2) (5;3) (5;4) (5;5) (5;1) (1;5), (1;1) (1;2) (1;3) (1;4) (1;5) (1;1) (2;5) I (2;1) (2;5) (2;1) i (3;5) ( (3;l) (3;5) 1 (3;1) (4;5) } (4;1) (4;5) I (4;1) l l (5;5) 1 (5;1) (5;2) (5;) (4) (5;5) (5;1) L _ -__. _ _ -_ _ - - - _ - (1;5) (1;1) (1;2) (1;3) (1;4) (1;5) (1;1) Figure 21. Illustration of cyclic boundary conditions on a 5 x 5 grid. 4.3 NUMERICAL METHODS From the initial conditions a forward difference time step was used to forecast q and q at t=1. This is an unstable differencing scheme and was only used for the first step since the Adams-Bashfort scheme requires that the fields of ql and q be known for two consecutive time steps. For this first step friction and nonlinear terms are neglected. The forward difference time scheme may be written as q = ql + H At q = q - H At (4.1) 3( 3

86 where the superscript denotes time step and At is the time increment. The integration proceeds using the Adams-Bashfort time differencing scheme as follows: t+l t t t-l q+ tq + [3/2 F - 1/2 Ft-/t (4.14) where F denotes the r.h.s. of (4.10). This scheme is slightly unstable (Lilly, 1965) and the amplification factor in the case of a periodic solution is given by 4 1/2 A = (1 + P/4 +... ) (4.15) where P = cAt and c is a characteristic wave speed. If we take 30 m/sec to be an upper bound to the wave speed and Ax to be 475 km, then to insure that a wave will move through Ax in t > 4At we must have At < 3900 sec. We have taken At = 3600 sec. This results in c = 30m/sec2jT/L and P = 2.2 x 102 which yields A - (1 + 10-7). This time differencing scheme has been used by Lilly (1969) and Orzag (1969) with satisfactory results. The Jacobians, -V * VT and -V ~ Vq, are evaluated using the Arakawa (1966) second-order scheme. It conserves kinetic energy and vorticity. Lilly (1969) noted that usage of a fourth-order scheme also having these conservative properties increased the aliasing in the short wavelength portion of the spectrum and otherwise little affected the results.

87 The Jacobian is written as follows. (-v-vq) -i [j. +q,*x' qi, 12(X)2 [(i, jm ip, jm ijp ip, jp qip i +( im, jm ti, jm im, jp-~i,jp)( i, j qim, j +( ip, j lP, jp - im, j im,j jp)(ijpi,j) Lt j imp,,jm im,j)(qi,j qijm) +( *ipj- i,jp) qip,jp-( i, jm- im,j)qim, jm +( tijp jim, j p( i ip, j )i, jm)qip,jm = J(qy)i (4.16) where ip = i+l, im = i-l, jp = j+l and jm = j-l. The linear spectrum is formed by summing the contributions of components whose two-dimensional wave number vector (m,n) describes a square about the origin with sides equal to 2., where m,n and Q are integers. Therefore, we may write 2 2-1 K () = Ke(m,2) + Z K (2,n) (4.17) m=l n=l

88 4.4 ENERGETICS OF THE MODEL 4.4.1 Formulation The formulation for the calculation of the energies and energy transformations for the model is very similar to that of Phillips (1956). The absence of a zonal flow and area means of the stream functions, 41 and. 3' result in the existence of only the eddy form of available potential energy and kinetic energy. The energy diagram pertinent to our model is depicted in Figure 22. G(A )_ Nk A C(A,K) K e e>~ e ee e e e (Ae) D(Ke) Figure 22. Energy flow diagram for the model. The equation for the time rate of change of available potential energy per unit mass is derived by multiplying (4.9) by (1 - if-) which results in a, _ 1 (\L a 2L L 2 A_ - ~ n r (T)2dxdy A - at e 2 2 M t o o 2 f X2 L L o LL - 2 o fo T(v VT)dxdy + 2 Jf f o2 Tdxdy L L P 2 v L L 2 1 L L -\2 - ( W 2 dxdy + 12o f (HT) dxdy L L (4.18) The first term an the r.h.s. of (4.18) represents the nonlinear cascade of available potential energy in wave number space. This term vanishes when summed over the domain if the boundary conditions restrict

89 the net transport of any quantity out of the domain to be zero. Therefore, the integral value of this quantity is not of importance, but its decomposition in wave number space is. The second term on the r.h.s. of (4.18) represents the conversion of energy between the available potential energy and the kinetic energy, and appears with opposite sign in (4.19). The last two terms in (4.18) may be thought of as contributing to the generation of the available potential energy. The diffusion term tends to smooth out the temperature gradients and always yields a negative contribution. The last term contains the product of the heating and the temperature fields. When heating occurs at relatively higher temperatures than the cooling a positive contribution to the time rate of change of available potential energy will result. The equation for the time rate of change of the kinetic energy per unit mass is obtained by multiplying (4.7a) by -Gu,l (4.7b) by -*f, adding and integrating. These manipulations yield* 1 L L1 2 2 -K - -f fL 1 [(\7V)2 + (V )2 dxdy] t t o eo 2 1(V ) + (v ) dxdy0 2 L o L L ) + ( dxdydy - 2 Jo Jo [/l(Vl vl( 3 f 2 L JL (o 2T) dxdy - V fLL ( 2 + 2 2 0 0 2 L2 o' ( 3 L P L (4. 19) *(4.18) and (4.19) may be multiplied by Po/g ~ 105 gm cm-2 to obtain units of erg cm-2 for en egies and erg cm2 sec- for energy transfers and conversions.

90 The first term on the r.h.s. of (4.19) represents the nonlinear cascade of kinetic energy in wave number space and has zero contribution to the time rate of change of total kinetic energy when summed over wave number. The second term on the r.h.s. of (4.19) has already been noted as being identical except for sign as a term in the equation for the time rate of change of available potential energy and represents the conversion between the two energy forms. The last term in (4.19) denotes the dissipation of kinetic energy due to lateral eddy viscosity and is always a loss. We also evaluate the enstrophy and the nonlinear cascade of this quantity in wave number space. The enstrophy is given by 1 L L1 2 2 1 L fL 1 [(V2 2) + (VA 3) ] dxdy (4.20) The rate of enstrophy gain due to the nonlinear cascade is given by 32 I(mjn2 oo (1 - IV + V5. )dxdy (4.21) n=l L I(mjn,A) is zero when integrated over the domain and only its decomposition in wave number space is of interest. If r, T, H, V'VT, VV:, co, and ~ are written as follows* *All double summations in (4.22), (4.23), and (4.24) have identical limits to j.

91 M M V = 2 M A(m n) eiJ(mx + ny) m=-M n=-M /0 #0 / O ~i(x + ny) T = E Z B(m,n) eiJ(mx+ ny) H = 2 Z C(m,n) eiJ( m + ny) V r = Z D(m,n) eiJ( mx + ny) V V = E Z E(m,n) eJ( + ny) = 2 2 Q(mn i(mxn + ny) (^ = 7, S R(m,n) e v / t = Z Z F(mn) eiJ(m ny) (4.22) then the areal mean of the quantities involved in the energy budget may be written as (* indicates complex conjugate) A = - B(m,n) B*(m,n) e 2 j2 2 2 K = 2 (m + n ) Ai(m,n) A*(m,n) (i = 1,3) f C(A,K ) = _~2 2 B(m,n) Q*(m,n) e e P

92 22 2 2 D(A ) = -2J2 Z Z (m + n2) B(m,n) B*( m,n) e D(K ) = v Z Z (m2 + n ) Ai(m,n) A*(m,n) (i = 1,3) e CK(mjn,l) = Z Z Ai(m,n) E*(m,n) (i = 1,5) K C (mln,2) = -X2 Z Z B(m,n) D*(m,n) G(Ae) = E Z C(m,n) B*(m,n) (4.23) Subscript i indicates summation over pressure surfaces 1 and 3. The enstrophy is C = Z~ F.(m,n) Fp(m,n) (i= 1,5) (4.24) while M Z I(mJn,) = - z Ei(m,n) Fi(m,n) (i = 1,3) m=l 4.4.2 The Diagnostic Calculations The time step used in (4.14) for the time integration was taken equal to 1 hour. The diagnostic calculations of the various terms of ( 4.18)-(4.23), inclusive, were performed every 12 time periods or 12 hours in real time. Integration of the model was attempted for various values of JHI and v. A value of.7 for R equally weights the past heating field and the newly formed component. This value was used in all the experiments discussed in the next section.

93 The value of JIH is particularly influential in the growth rates of the other waves during the early quasi-linear period. For example, a factor of two in JHi results in a factor of about 16 in the growth rate of wave number 32. Phillips (1963) gives the limitation on the heating rate consistent with the quasi-geostrophic approximation to be -1 2 -35 about 10 m sec. The values of the heating rate associated with the values of IHI used in this study are listed in Table 8 and are within the limit given by Phillips. TABLE 8 HEATING AMPLITUDE AND COEFFICIENT OF EDDY VISCOSITY Case A Case B Case C dQh/dt(x 10 cm2 sec ) 70 70 21 v (x 109 cm2 sec-1).85 1.7.85 From our observational calculations, values of the coefficient of 9 2 -1 eddy viscosity near 109 cm sec seem reasonable. The values of v used in the three cases to be discussed are contained in Table 8. Although the continuous equations (4.18) and (4.19) have exact energy integrals the finite difference forms do not because of the spatial and temporal finite difference approximations used in both the prognostic and diagnostic equations. Further, since the diagnostic calculations are performed once for every 12 steps of the prognostic calculations we cannot expect precise correspondence between the energy changes evaluated by the l.h.s. of

94 (4.18) and (4.19) and their r.h.s. A measure of the representativeness of the diagnostic calculations combined with truncation and aliasing errors may be obtained by comparing the energy change as calculated directly from the stream function fields with the budget obtained from the energy generation and dissipation. Table 9 contains these values for the cases under discussion. E denotes total energy. TABLE 9 CHANaES IN TOTAL ENERGY COMPARED WITH RESIDUE OF ENERGY BUDGET CALCULATIONS (x 109 cm2 sec-2) t (hours) 24 500 600 900 1200 1500 Case A E 2.1 20.35 1.1 45.1 56.6 AE 18.1 10.8 12.0 13.5 AE(budget) 14.8 15.9 13.9 14.8 Case B E 2.1 15.3 22.4 30.2 59.2 45.1 AE 13.2 7.2 7.8 8.9 6.0 AE(budget) 12.5 11.8 8.8 11.8 10.0 Case C E.2 1.6 2.6 353 359 4.4 AE 1.4 1.1.7.6.6 AE(budget) 1.3.5.5 1.2 1.1 The comparisons indicate that the diagnostics are representative. Case C has the best agreement with only a 2.5% difference between nE

95 and AE(budget) after 1200 time steps and a 9% difference after 1500 hours. Case B yields the poorest agreement with a 20 difference after 1200 hours and a 27% discrepancy after 1500 hours. Case A shows 9j difference after 1200 hours. There seems to be no particular reason for the variation in the difference between the energy calculated using the r.h.s. of (4.18) and (4.19) and that calculated by the l.h.s. of these two equations. These results may be compared with those of Phillips (1956) (realizing, however, the additional problems which existed at that time) in which the discrepancies reached 100% by about 600 hours and continued to increase thereafter. Smagorinsky (1963) also compares the change in energy with the energy budget for a 22-day period. It seems, however, that he is calculating the diagnostics at every time step and so he avoids the problem of representativeness. He obtains a 3.97% difference per day for the eddy energy and a.18% difference per day for the total energy between the energy change and the budget. Similar quantities for total energy are.355 day,.11% day and.15% day for Cases A, B, and C, respectively. These figures are the averages taken over the entire duration of the runs. If a similar analysis is done on the energy budget as a function of wave number the correspondence is less good. The specific values of w and therefore the conversion of available potential energy to kinetic energy are relatively less reliable than any other process. This is

96 aside from the question of representativeness. The time derivative of the temperature field is evaluated over a 24-hour period while all other factors involved in the calculation of cX are evaluated from their instantaneous values. Although the individual wave number values of C(A,K ) are, therefore, less reliable than the other processes, nevertheless, the shape of the curve is considered to be representative and meaningful. It can also be shown that the changes in kinetic and available potential energy due to the nonlinear transfer processes do not exactly balance the dissipation and energy changes with respect to wave number. Here again, there exists a factor aside from the question of representativeness and that is the response of the finite difference representation of the Jacobian. As mentioned previously and observed in the calculation, the Arakawa second-order Jacobian is conservative if the boundary conditions allow zero net flux into the domain. However, the finite difference Jacobian undoubtably distorts the transfer in wave number space. Evidently this distortion does not obscure the main features of the nonlinear transfers. It would be of interest to compare the nonlinear energy gains as determined by finite differences using the Arakawa Jacobian to that evaluated by the summation of triple products, as done in a previous chapter on the observational calculations. The general distribution of nonlinear energy gain, both for kinetic and available potential energy, resembles observational calcula

97 tions (Yang, 1967) and leads to a reasonable qualitative picture of the energy cycle if not an exact quantitative distribution. 4..4.3 The Energy Cycle and Nonlinear Transfers We shall use the results of Case A to illustrate the following discussion. Figures depicting the results of Cases B and C are contained in Appendix E. Figure 22 contains the generalized energy flow diagram for the model. There is an initial start up period which lasts for 200 to 300 time steps, depending upon the particular case, during which both forms of energy increase rapidly (see Figure 23). Thereafter the available potential energy decreases somewhat and then is quasi-stationary or very slowly decreasing with time. The eddy kinetic energy, on the other hand, continues to increase with time. The relationship between K and e time may be expressed as xt K = Ce 600 < t < 1200 e - - where C and a are constants and are listed in Table 10. The continual increase in K is due to growth for m < 8 due to the nonlinear trapping effect discussed below. effect discussed below.

98 10 - 8 6 K e 4 - 2 CJ o 1010_ Q) 8 D 6 0) 10 2 109 8 240 480 720 960 1200 1400 t (hours or time steps) Figure 23. Ke and Ae as a function of time for Case A.

99 TABLE 10 EVALUATION OF C AND a IN K = Ce,600 < t < 1200, USING K AT t = 600 AND 1200 e Case A B C C(107 erg cm ) 1458 1061 130 (10-4) 4.7 4.5 3.4 The generation of available potential energy is forced, in our model, to be positive at m = 8 and negligible at all other wave numbers. In the atmosphere G(A ) is predominantly negative due to radiational effects almost everywhere and sensible heat transfer near the lower boundary (Lawniczak, 1969). Precipitation effects are usually positive. The source for A in the atmosphere is generally the conversion from e A. Since we exclude zonal forms of energy in our model we allow G(A ) z e to be positive. The energy input to A is through G(A ) which is exclusively of e e wave number 8. The nonlinear transfer of this quantity if predominantly toward large m. Therefore, we do not observe the loss, due to nonlinear effects, of A for m < 8 as observed in the atmosphere (Yang, 1967). The conversion from A to K, C(A,K ), in this spectral e e e e region should not be expected to and does not approximate atmospheric behavior. The nonlinear gain of A for m > 8 is predominantly removed by e C(A,Ke) and to a lesser degree by D(A ) the dissipation effects. In e e <

100 genaeral, for m > 8, C(A,K) K C(mJn/,) + D(A ) with D(Ae).2C(A,K). 9 ~~~~e e A e e e e) The mean values of these processes for 300 < t < 1200 as a function of wave number are exhibited in Figure 24. The eddy kinetic energy, Ke, in the spectral region 1 < m < 7 receives energy by the nonlinear transfer process but our neglect of C(K,K ) and mispresentation of C(A,K ) prevents this spectral region e z eI~) revnt ths ~ reioe e from resembling the atmosphere. Instead, K for 1 < m < 7 in our model should more closely correspond to the two-dimensional theoretical model than the atmosphere since the nonlinear transfer dominates the energy gains. D(K e) the dissipative effect is small and C(Ae,Ke) is small and has alternating sign. As the spectrum of K for m > 8 attains e quasi-stationarity the dominant kinetic energy transfer process is that of trapping energy at m < 8 and chiefly at m = 1. The mean values of CK(mjn,e), C(A,K ) and D(K ) for 300 < t < 1200 are shown in Figure 25. For m > 8, the K spectrum in our model should resemble the atmose pheric spectrum since we account for C(A,K ), C (mln,2), and D(K e) e e K' in a reasonable fashion. C(Ke,K ), which is not accounted for in our model, is relatively smaller in the atmosphere for m > 8. The effect of the nonlinear exchange of K is a loss for 8 < m < 14 e and small gains for m > 14. For 9 < m < 14 a balance is achieved among C (mjn,2), C(A, Ke) and D(K ) with the three quantities of the same order of magnitude. It is noted that the kinetic energy supply for 8 < m < 17 in our model is not C (mln,a) but it is the balance of all the involved

101 18 - - C-(peKe) 16 -1 ---- (m n,l) 14 - -D(A e e 12 - u \ o 84 _ - 6 H Q x 4 0 0; - 6 - ~IJ^ ~ -". 6- I Fz - H F -12q -14 for m = 8 -16 G(Ae) = 29,580 -C(A,Ke) = 21,020 -18" CA(m n,2) = -8336 2 4 6 8 10 12 14 16 18 20 22 24 26 Wave Number Figure 24. Mean values of CA(mln,A), -C(Ae,Ke), and D(Ae) for 500 < t < 1200, Case A.

102 36 - C(A,K )e 32 - --- CK(m n,)l K 28 -.- D(Ke) ~~28 - ~e I 24 - H 20 \ C) o 16 1S Q 12 - i —1 \s~' = -'2 - 8 CM -i -16 -20 <~~~C,~~~~~ ~C(mn) 1-11 460 2 4 6 8 10 12 14 16 18 20 22 24 26 CH(mln,l) -11,460 Wave Number Figure 25. Mean values for CK(mln,~), C(Ae,Ke), and D(Ke) for 300 _ t 1200, Case A.

103 processes which maintain the quasi-stationary spectral slope. In the atmosphere C(A,K ) (Saltzman and Fleisher, 1960; Wiin-Nielsen, 1959) is larger than CK(mln,l) (Yang, 1967) for this spectral range except K perhaps the short wave end where both are small. The nonlinear enstrophy transfer is predicted by the theory to be constant and toward large wave number in the postulated inertial subrange. We find the direction of transfer to be verified in the model. The losses of enstrophy occur for 8 < m < 14 while the gains take place mainly in the spectral region with m > 14. Figure 26 depicts I(mjn,^) as a function of wave number for Case A. The values are the means for 300 < t < 1200. This figure may be compared with Figure 4 which contains the observational results. There is very close correspondence between the action of the nonlinear transfer of enstrophy in the atmosphere and in the model if the wave numbers at which the sign changes are used as reference points. The exception exists for 1 < m < 7 in the model's results. Since the processes in the spectral region are misrepresented (as discussed in connection with the available potential energy spectrum) this region should not and does not have correspondence with the atmospheric case. The results for Cases B and C are contained in Appendix E. 4.4.4'Power Laws' After the initial start up period and overshoot, the kinetic energy at each wave number for m > 8 may be described as quasi

800X -200 - 0; - 4 ~ ~ - 2 0 0 o -00I H -800- I(81n,~) = -5072. I t I! I I I I I I s I I I I 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 Wave Number Figure 26. Mean value of I(mln,e), 300 < t < 1200, Case A.

105 stationary. The kinetic energy for m = 8, 12, 16, and 24, Case A, is exhibited as a function of time in Figure 27. Similar figures for Cases B and C are contained in Appendix E. A small positive trend with time is observed for the kinetic energy at m = 16 and 24 for Case A beginning at t = 1000. This is not observed in Cases B and C in which the quasi-stationary characteristic is more closely adhered to. The kinetic energy spectra at t = 600 and 1200 (Case A) are depicted in Figure 28. It is noted that changes for m > 8 are small whereas the major changes occur for m < 5 and particularly at m = 1, 2, due to the nonlinear trapping effect mentioned in the previous section. This holds true, if anything, to a greater degree in Cases B and C (see figures in Appendix E). Table 11 lists the spectral slopes calculated by least squares of the kinetic energy for 8 < m < 16 as a function of time for the various cases. TABLE 11 KINETIC ENERGY SPECTRAL SLOPES 8 < m < 16 AS A FUNCTION OF TIME (x-l) Time 300 600 900 1200 1500 A 3.1 2.2 1.8 1.7 B 3.6 3.1 3.2 2.8 2.8 C 6.6 3.6 3.6 3.3 3.4 We have obtained three quasi-stationary kinetic energy distribu

106 110 L _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 109 - 1c 0'? _ m=8 b 10 a O a) C) 8 l 10 — 1,, I P/^v,24 107 l I I I I - 240 480 720 960 1200 1440 t (time steps or hours) Figure. 27. The kinetic energy of wave numbers 8, 12, 16, and 24 as a function of time, Case A.

107 1010 - \ t= 12 \ 7- t = Coo\ o \ / hi) so I6 I II X| XH -.... t = 1200 V 2 4 6 8 10 20 30 Wave Number Figure 28. The kinetic energy spectrum number at t = 600 and 1200, Case A.

108 tions in wave number space (Figure 29) by altering the relative intensity of the generation of energy and its dissipation. The nonlinear effects are generally smaller than either of the dissipation or the conversion from available potential energy. Therefore, we must conclude that the arguments of the two-dimensional turbulence theory are not appropriate to deducing the kinetic energy spectral slope in our model. The relative magnitudes of C(A,K ) and CK(mln,e) in the experie e K ments conducted in this study are similar to those calculated from atmospheric data (Yang, 1967; Wiin-Nielsen, 1959; Saltzman and Fleisher, 1960; Saltzman and Teweles, 1964). It is, therefore, unlikely in the atmosphere as well as in our model that the two-dimensional arguments are appropriate, even though the spectral slope it predicts is closely adhered to in the atmosphere and our model. In particular, one should not take the existence of a'-3 power law' in the atmospheric kinetic energy spectrum for 8 < m < 16 as a verification that the assumptions used in the two-dimensional turbulence theory such as isotropy are applicable in the atmosphere. The available potential energy for Cases A, B, and C fit'a -5 power law' fairly well for 8 < m < 20 (Figure 30). In our diagnostic computations the enstrophy is defined as being proportional to the wave number squared times the kinetic energy. If K is represented by a'-2 power law' as in Case A, then the enstrophy e is approximately constant with respect to wave number for 8 < m K 16.

109 1lo-o\ ~ C o \X' \<^ 108\ 9-CaseA \\ \ 7 Case B 10 Figure 29. The kinetic energy spectra at t = 1200 for Cases A, B, and C.

110 8 10 a_ - C A \\ bj IX I Q) 6 1 10 d \ -- Case A — Case B 10 Case C 2 4 6 8 10 20 50 Figure 50. The spectra of A at t =1200 for Cases A, B, and C.

CHAPTER V CONCLUDING REMARKS 5. 1 SUMMARY AND CONCLUSIONS The nonlinear effects on the enstrophy and potential enstrophy were evaluated for a three-month period, February 1963 - April 1963, inclusive. The nonlinear transfer is toward large wave numbers, affirming that prediction of the theories of Kraichnan and Leith. The inertial subrange through which enstrophy should flow quantitatively unaltered, also postulated by the theories, is not verified by the calculations of this study. The enstrophy cascade rate, r, -18 -3 from m < 7 to m > 7 is 172 x 10 sec while the cascade rate for -18 -3 potential enstrophy is 253 x 10 sec. Other factors entering into the enstrophy and potential enstrophy budgets are evaluated. These are fluxes due to the interaction among waves and due to wave-zonal mode interaction. The former process results in loss rates of -18 -3 -18 -3 25 x 10 sec and 24 x 10 sec while the latter process re-18 -3 -18 -3 sults in gain rates of 45 x 10 sec and 290 x 10 sec for enstrophy and potential enstrophy, respectively. The exchange of enstrophy and potential enstrophy between waves and the zonal flow due to their mutual interaction results in -18 -3 -18 -3 40 x 10 sec and -171 x 10 sec rates of change for the waves, respectively. The rates of change due to the n-effect is 111

132 -17 - x -18 -3 -174 x 108 sec and -68 x 10 1sec for the enstrophy and potential enstrophy, respectively. The above results when decomposed in the time domain shows that wave-wave interactions take place dominantly through the fast transient mode while in wave-zonal mode interaction as well as the 0-term effect the slow transient mode is dominant. These results hold for both the exchanges and fluxes of the enstrophy and potential enstrophy. The vertical variation of the absolute values of all evaluated processes indicates maxima at 20 cb < p < 30 cb. P (m) p also has a relative maximum in the layer 60 cb < p < 77.5 cb. When a power law is fitted to the transient enstrophy on a log-log-scale a slope of -1.067 results in the spectral region 8 < m < 16, after a correction factor (Wiin-Nielsen; 1967) is applied to counter the effect of the smoothing due to the finite difference techniques (Ogura; 1957) used in the evaluation of geostrophic winds from the height field. For the transient corrected potential enstrophy the slope is -.997. These values are not significantly different from -1 and correspond to a kinetic energy spectral region following a -3 slope. Using Leith's (1968b)formulation for the enstrophy dissipation, T, the coefficient of eddy viscosity was evaluated to be about 2. 5 8 2 -1 x 10 cm sec. In the absence of intermediate sources or sinks a cut-off wavelength was calculated using Kraichnan's formulation to be 70 km.

113 A transfer coefficient for potential vorticity was sought following Green's suggestion. Whereas Wiin-Nielsen and Sela (1971) found a positive transfer coefficient for yearly and monthly mean fields of potential vorticity this was found not to be so for the daily cases. A quasi-geostrophic two-layer model was integrated, for 1200 to 1500 hours depending on the case. The lateral boundary conditions restrained the fields to be cyclic in both the x - and y - directions. The vertical boundary conditions were X = 0 at p = 0 cb and p = 100 cb. The energy was introduced by a heating field whose harmonic con-17 stituents were zero (~ 10 x H) except those whose larger wave number equalled 8. The total amplitude of the heating field was constant but the relative amplitudes and phase angles of its components were changed at each time step. Every twelve time steps the energies and the various processes involved in the energy budget were evaluated as well as the enstrophy and its nonlinear transfer in wave number space. The energy cycle is similar to that observed in the atmosphere except in the spectral region with m < 8. A, generated at m = 8, is both converted to K and transferred, for the most part, toward e components with m > 8. This is in agreement with cascade direction observed in the atmosphere. The available potential energy attains a quasi-stationary state by t 3 00 and does not change significantly thereafter. The power-law governing the spectral distribution of

ll4 A would require an exponent of about -5. e Available potential energy is converted into kinetic energy for all components with m > 8. The nonlinear transfer of K is from waves "- e with 8 < m < 14 predominantly toward the large scale portion of the spectrum. There is a relatively small transfer toward m > 14. After an initial period of about 300-400 time steps K for m > 8 becomes e quasi-stationary while for m < 8 and, in particular, m = 1, 2 growth continues by virtue of the nonlinear transfer. This portion of the spectrum will resemble a -5/3 slope if the process is allowed to continue. In the atmosphere this growth is countered by the transfer of kinetic energy from the large scale eddies to the zonal mode. In Case A Ck(mln,e), C(A, K ) and D(K ) are all of the same k e e e order of magnitude for 8 < m < 16. In Cases B and C which more closely resemble the atmosphere with respect to spectral slope of kinetic energy, the nonlinear transfer of K contributes the smallest e effect to the kinetic energy budget. The relative effects of CK(mjn,i) and C(A, K ) correspond qualitatively to the results obtained from ate e mospheric data (Yang, 1967; Wiin-Nielsen, 1959). In the three cases discussed in Chapter IV three somewhat different quasi-stationary spectral slopes are obtained. These vary from about -1.7 to -3.4 depending upon both the intensity of energy input and dissipation. In none of the cases is an "enstrophy inertial subrange" observed. In view of these results it is concluded that the closeness with which atmospheric kinetic energy follow a'-3 power law' should not be explained by arguments postulating inertial

115 subranges. The quasi-stationary kinetic energy spectral slope seems to be affected mainly by the relative intensity of energy sinks and sources at each wave number and, to a lesser degree, by the nonlinear effects. It would seem that it is not as important to ascertain the precise spectral slope as was believed 6-12 months ago. Leith (1969) and Fleming (1970) have shown that the -3 slope is not critical for the predictability question as previously postulated by Lorenz. Here, we show, by the simple quasi-geostrophic model, that there is nothing magical about the -3 slope but that other quasi-stationary slopes may be attained by the kinetic energy spectrum as the sinks and sources of energy are modified. Thus, it would seem to be of more importance to determine the distribution of sources and sinks than to measure the resulting spectral slope. The atmosphere may be quasi two-dimensional with respect to the prognostic equations and the magnitudes of the three components of the velocity, however, it is not so with respect to the energy cycle. It is concluded that the two-dimensicnal turbulence theories (Fjortoft, 19553; Kraichnan, 1967; and Leith, 1968a)pertain to the model and the atmosphere only to the extent that they predict the dominant direction of the nonlinear cascades of kinetic energy and enstrophy.

116 5. 2 SUGGESTION FOR FUTURE INVESTIGATION Observational studies should be extended in time so that more reliable estimates of the nonlinear processes may be obtained. These studies should also be extended in wave number space. Since the -1 enstrophy is proportional to m it cannot be said that little enstrophy exists form > 15. The observational studies might very well be carried out in the framework of spherical harmonics and so avoid the filtering effects accompanying finite-differencing the derivatives with respect to the north-south direction. A study incorporating a greater range of wave numbers will also permit better definition of the coefficient of eddy viscosity, v, and especially r, the enstrophy cascade rate. If a study be carried out for an annual cycle then the interactions can be decomposed in frequency space and a more precise statement on the interactions in the space-time domain would result. Further studies simulating realizations of turbulence using the model described in Chapter IV should be carried out for a greater range of a parallel integration carried out in the Fourier domain, as opposed to the grid-point calculation used herein, would be of interest (Orzag, 1969). It would probably require more computer time and coding but would eliminate the filtering effects of the spatial finite difference schemes used to approximate the Jacobian and frictional effects in the prognositc equations. The case of decaying turbulence may be studied by simply

117 setting the heating to zero at some arbitrarily chosen time. The generation of available potential energy will then cease while the dissipative effects will continue. Will the nonlinear effect on the kinetic energy reverse and drain the energy from the spectral range of m < 8 or will the kinetic energy remain trapped due to the inability of the kinetic energy to be carried nonlinearly toward large m? In order for the model to shed light on the similarities and differences between the theories and atmosphere observation its hybrid nature should be retained, as opposed to tuning the model to the atmosphere. I reiterate the suggestion offered by Smagorinsky (1963) which is that experiments should be conducted varying the parameters such as the heating rate and the viscous coefficients. This has not been done to any great extent either with a simple model (as used in this study) or with models of greater complexity, Perhaps, in this way climatological theories may be tested, specifically those involving changes in the solar constant. If a quasi-steady state is arrived at for a given choice of heating then the case would represent a possible realization of the atmospheric circulation. This would not attest to the veracity of the theory, however, it may present means of limiting the number of possibilities.

BIBLIOGRAPHY Arakawa, A., 1966: Computational Design for Long-Term Numerical Integration of the Equations of Fluid Motion: Two-Dimensional Incompressible Flow, Part I. J. Comp. Physics 1, 119-143. Batchelor, G. K., 1969: Computation of the Energy Spectrum in Homogeneous Two-Dimensional Turbulence. Phys. Fluids Suppl. II, 12, II-233-II-239. Bradley, J.H.S., 1967: The Transient Part of the Atmospheric Planetary Waves. The University of Michigan. Technical Report 06372-5-T. Bradley, J.H.S., and A. Wiin-Nielsen, 1968: On the Transient Past of the Atmospheric Planetary Waves. Tellus 20, No. 3, 533-544. Charney, J. G., 1959: On the General Circulation of the Atmosphere. Rossby Memorial Volume. The Rockefeller Institute Press, New York, 178-1953. Charney, J. G., 1966: Some Remaining Problems in Numerical Weather Prediction. Advances in Numerical Weather Prediction. 1965-66 Seminar Series. The Travellers Research Center. Inc. Derome J. F., and A. Wiin-Nielsen, 1966: On the Baroclinic Stability of Zonal Flow in Simple Model Atmospheres. The University of Michigan. Technical Report 06372-2-T. Fj~rtoft, R., 1953: On the Changes in the Spectral Distribution of Kinetic Energy for the Two-Dimensional, Nondivergent Flow. Tellus 5, 225-230. Fleming, R. J., 1970: Concepts and Implications of Stochastic Dynamic Prediction. N.C.A.R. Cooperative Thesis No. 22. Gates, W. L., 1960: Static Stability Measures in the Atmosphere. Scientific Report No. 3, Dynamical Weather Prediction Project. Department of Meteorology, UCLA (AFCRL-TN-60-817). Green, J.S.A., 1970: Transfer Properties of the Large-Scale Eddies and the General Circulation of the Atmosphere. Quart. J. Roy. Meteorol. Soc. 96, 157-185. 118

119 BIBLIOGRAPHY (Continued) Griffith, H. L., H. A. Panofsky, and I. Van der Hoven, 1956: Power Spectrum Analysis Over Large Ranges of Frequency. J. Meteorol. 13, 279-282. Holopainen, E. 0., 1967: On the Mean Meridional Circulation and the Flux of Angular Momentum Over the Northern Hemisphere. Tellus 19, 1-13. Horn, L. H. and R. A. Bryson, 1963; An Analysis of the Geostrophic Kinetic Energy Spectrum of Large-Scale Atmospheric Turbulence. J. Geophys. Res. 68, 1059-1064. Julian, P., W. Washington, L. Hembree, and C. Ridley, 1970: On the Spectral Distribution of Large-Scale Atmospheric Kinetic Energy. J. Atmos. Sci. 27, 376-387. Kao, S. K., 1968: Governing Equations and Spectra for Atmospheric Motion and Transports in Frequency, Wave-Number Space. J. Atmos. Sci. 25, No. 1, 52-58. Kao, S. K., C. Y. Tsag, and L. L. Wendell, 1970: The Meridional Transport of Angular Momentum in Wave Number-Frequency Space. J. Atmos. Sci. 27, No. 3, 619-626. Kasahara, A., and W. M. Washington, 1967: NCAR Global General Circulation Model of the Atmosphere. Mon. Wea. Rev. 95, 389-402. Kraichnan, R. H., 1967: Inertial Ranges in Two-Dimensional Turbulence. Phys. Fluids 10, 1417-23. Kung, E. C., 1967: Diurnal and Long-Term Variations of the Kinetic Energy Generation and Dissipation for a Five-Year Period. Mon. Wea. Rev. 95, 59 5-606. Lawniczak, G. E., Jr., 1969: On a Multiple-Layer Analysis of Atmospheric Diabatic Processes and the Generation of Available Potential Energy. The University of Michigan. Technical Report 08759-5-T, 111 p. Leith, C. E., 1968a: Diffusion Approximation for Two-Dimensional Turbulence. Phys. Fluids 11, 671-673. Leith, C. E., 1968b: Two-Dimensional Eddy Viscosity Coefficients. Proc. WMO/IUGG Syp. Num. Wea. Pred., Tokyo, 41-44.

120 BIBLIOGRAPHY (Continued) Leith, C. E., 1970: Atmospheric Predictabilityand Two-Dimensional Turbulence. N.C.A.R. MS 70-87. 55pp. (Published in J. Atmos. Sci.) Lilly, D. K., 1965: On the Computational Stability of Numerical Solutions to the Time Dependent Nonlinear Geophysical Fluid Dynamics Problems. Mon. Wea. Rev. 93, 11-26. Lilly, D. K., 1969: Numerical Simulation of Two-Dimensional Turbulence. Phys. Fluids Suppl. II, 12, II-240-II-249. Lorenz, E. N., 1969: The Predictability of a Flow which possesses many Scales of Motion. Tellus 21, 289-307. Manabe, S. J., Smagorinsky, J. L. Holloway, Jr., and H. M. Stone, 1970: Simulated Climatology of a General Circulation Model with a Hydrological Cycle. Mon. Wea. Rev. 98, No. 3, 175-212. Ogura, Y., 1957. Spectrum Modification Due to Use of Finite Differences. J. Meteor. 14, 77-80. Ogura, Y., 1958: On the Isotropy of Large-Scale Disturbances in the Upper Troposphere. J. Meteor. 15, 375-382. Ogura, Y., 1962a: Energy Transfer in a Normally Distributed and Isotropic Turbulent Velocity Field in Two Dimensions. Phys. Fluids, 5, 595-401. Ogura, Y., 1962b: Energy Transfer in an Isotropic Turbulent Flow. J. Geophys. Res. 67, 3143-3149. Orzag, S. A., 1969: Numerical Methods for the Simulation of Turbulence. Phys. Fluids Suppl. II, 12, II-250-II-257. Phillips, N. A,, 1956: The General Circulation of the Atmosphere: A Numerical Experiment. Quart. J. Roy. Meteor. Soc. 82, 123-164. Phillips, N. A., 1963: Geostrophic Motion. Rev. Geophys. 1, 123-176. Pinus, N. Z., E. R. Reiter, G. N. Shaw, and N. K. Vinnichenko, 1967: Power Spectra of Turbulence in the Free Atmosphere. Tellus 19, 206-213. Reiter, E. R., 1969a: Atmospheric Transport Processes. U.S. Atomic Energy Commission, Division of Tech. Information, Oak Ridge, Tenn.

121 BIBLIOGRAPHY (Continued) Reiter, E. R., 1969b: Mean and Eddy Motions in the Atmosphere. Mon. Wea. Rev. 97, 200-204. Saltzman, B., and A. Fleisher, 1960: The Modes of Release of Available Potential Energy in the Atmosphere. J. Geophys. Res. 65, 12151222. Saltzman, B., and A. Fleisher, 1961: Further Statistics of the Modes of Release of Available Potential Energy. J. Geophys. Res. 66, 2271-2273. Saltzman, B., and S. Teweles, 1964: Further Statistics on the Exchange of Kinetic Energy Between Harmonic Components of the Atmospheric Flow. Tellus 16, 432-435. Smagorinsky, J., 1963: General Circulation Experiments with the Primitive Equations: I. The Basic Experiment. Mon. Wea. Rev. 91, 991-1064. Smith, P. J., 1970: A Note on Energy Conversions in the Open Atmospheric Systems. J. Atmos. Sci. 27, 518-521. Wiin-Nielsen, A., 1959: A Study of Energy Conversion and Meridional Circulation for the Large-Scale Motion in the Atmosphere. Mon. Wea. Rev. 87, 319-332. Wiin-Nielsen, A., 1959: On Certain Integral Constraint for the Time Integration of Baroclinic Models. Tellus 11, No. 1, 49-59. Wiin-Nielsen, A., 1961a: On the Distribution of Temperature Relative to Height in Stationary Planetary Waves. Tellus 13, No. 2, 127139. Wiin-Nielsen, A., 1961b: A Preliminary Study of the Dynamics of Transient, Planetary Waves in the Atmosphere. Tellus 13, No. 3, 320333. Wiin-Nielsen, A., J. A. Brown, and M. Drake, 1963: On Atmospheric Energy Conversions Between the Zonal Flow and the Eddies. Tellus 15, 261-279. Wiin-Nielsen, A., J. A. Brown, and M. Drake, 1964: Further Studies of Energy Exchange Between the Zonal Flow and the Eddies. Tellus 16, 168-180.

122 BIBLIOGRAPHY (Concluded) Wiin-Nielsen, A., 1965: Some New Observational Studies of Energy and Energy Transformations in the Atmosphere. Proceedings of the Symposium on Research and Development Aspects of Long-Range Forecasting, Boulder, Colorado. W.M.O. Technical Note 66, 177-202. Wiin-Nielsen, A., 1967: On the Annual Variation and Spectral Distribution of Atmospheric Energy. Tellus 19, 540-559. Wiin-Nielsen, A., and J. Sela, 1971: On the Transport of Potential Vorticity. To be published. Yang, C. H., 1967: Nonlinear Aspects of the Large-Scale Motion in the Atmosphere. The University of Michigan. Technical Report 08759-1-T.

APPENDIX A VERTICAL FINITE DIFFERENCE SCHEME FOR THE POTENTIAL VORTICITY The potential vorticity may be written in the form M _ imk p = U + Z U(m)e (1) m=-M The coefficients U(m) may then be written in terms of the vorticity and the stream function as U(m) = F(m)+ f2 A(m) (2) on p a ap / Let 2 1A(m1 Z - (m>)=f ( - ( (X(m) - iY(m)) then Z(-m) = X(m) + iY(m) and X(m) = f <( ) 0 ap c ap 2 6 lC (m~ Y(m) = f2 a ( ap The finite difference form of X(m) in terms of B(m) is 123

124 B(m)-B(m) 1 -1/2 B(m) - B(m) Ai +1/2 X(m) = 2 2 1 a +1/2 i+1/ i - 11/2 2f () Af i+1/2+ i-1/2)2 where the superscripts refer to the pressure surfaces as indicated in Table Al. TABLE Al QUANTITIES USED IN POTENTIAL VORTICITY CALCULATIONS I P(cb)'(10-6 kg2 m4 sA(cb) Winter Spring 2 10.4 95.7 5 3 20 } 54.5 42.8 5 4 }17.9 13.3 1 0 ) 3.73 3.70 20 5 50 3 2.26 2.23 20 6 0 3 2.24 1.77 15 78 10 } 1.88 1.75 15 8 100 A corresponding equation may be written for Y(m) substituting C(m) for B(m) in (3).

APPENDIX B RELATIONSHIP BETWEEN KINETIC ENERGY AND ENSTROPHY Assume the streamfunction 4 to be composed of only one component of wave numbers m and n in the x- and y-directions, respectively. Then we can write t(x,y) = A(m,n)e x ny) where J = 2-. The kinetic energy per unit mass may be written as K(m,n) = 1 f ds ~S S 2 22 2 = J (m + n ) A(m,n) A*(m,n) The enstrophy is o 2 S(m,n) = - f ds 4( 2 2 A(mn) A(mn) = J (m + n ) A(m,n) A*(m,n) 2 2 2 If we define an effective wave number m such that m = m + n e e then we see that 2 g(m,n) a K(m,n) m e We should not expect this relationship to hold exactly for our computations since the derivative in the y-direction is evaluated by finite differences. 125

APPENDIX C MONTHLY MEANS AND STANDARD DEVIATIONS OF THE PROCESSES AFFECTING THE ENSTROPHY AND POTENTIAL ENSTROPHY The large standard deviations associated with the quantities calculated in this study were, to a degree, expected. They indicate that further calculations are a necessity if conclusions are to be drawn with a greater degree of certainty than is possible in the present study. Calculations covering at least an annual cycle should be strongly considered. 126

TABLE C1 MEANS AND STANDARD DEVIATIONS OF PROCESSES AFFECTING ENSTROPHY (x 10-18 sec-3) Function Month m = 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Mean - 83.4 - 31.4 4.1 38.1 -142.6 - 6.0 12.5 - 23.0 - 18.8 - 9.4 - 29.5 42.0 43.9 3.5 118.8 February S.D. 284 256 295 292 218 200 286 299 276 271 210 234 216 144 168 Mean 7.1 - 66.8 - 48.6 - 16.0 44.2 - 5.2 - 62.8 - 2.0 19.8 57. - 7.2 48.2 29.5 1. - S.D. 256 236 293 239 256 306 259 302 281 258 234 21 210 208 196 Mean - 31.5 - 52.5 - 37.3 - 8.7 16.2 -110 - 17.8 3.0 2.3 7.6 27.0 17.0 41.4 100.1 47.5 Ap S.D. 198 263 224 196 202 287 274 288 269 324 254 209 202 173 19 Mean 7.5 12.4 - 9.4 - 37.0 18.6 - 13.9 - 14.5 4.6 16.5 2.8 - 14.4 1.7 4.1 6.2 56 February S.D. 75 90 75 89 81 78 99 83 101 78 90 95 91 64 84 F(mln,l) March Mean 17.8 2.8 17.9 - 14.7 13.0 - 75.5 56.8 14.2 - 29.7 - 31.2 - 20.5 - 30.9 23.8 20.7 0.7 S.D. 89 95 151 123 118 223 147 175 136 188 113 121 103 102 165 Mean 24.6 7.3 - 5.2 - 23.1 6.2 40.3 1.8 - 30.9 - 15.5 - 13.1 - 5.9 17.1 - 8.5 - 13.0 20.0 S.D. 68 88 76 85 84 88 118 131 163 114 154 127 155 121 89 Mean - 18.6 0.1 78.6 - 11.5 34.6 33.4 18.7 3.0 - 15.5 11.1 4.4 - 3.3 2.5 - 1.5 2.4 Februar S.D. 53 66 125 78 81 65 65 45 40 45 18 15 11 8 12 I(o) March Mean - 7.4 - 48.2 38.0 24.8 10.0 12.6 13.3 8.4 - 57 6.4 2.0 - 5.7.4.4 S.D. 32 52 96 81 62 46 52 46 32 26 17 16 9 8 6 Mean -.9 - 14.2 -.5 - 3 12.9 - 28.4 2.1 3.9.7 -.8 - 6.4 - 1.2 -.1 -.6 1.0 S.D. 31 46 58 47 82 66 48 41 27 31 22 16 12 7 8 Mean 15.8 10.9 29.4 24.7 3.4 - 6.9 - 1.4 2.4 2.1 - 3.5 0.0 -.9 -.3 -.9 1.0 Febru S.D. 26 24 29 29 16 23 17 16 13 14 9 65 6 F(o) March Mean 7.4 - 2.9 25.8 12.1 - 8.9 -.1 - 1.5 4. -.4 - 2.9 1.2.0.8.1.6 S.D. 20 20 27 23 24 22 14 25 10 11 9 8 4 4 April Mean 5.0 7.6 4.1 6.7 - 1.4 - 1.5 -.3 - 2.1 1.4 2.3 1.8 2.5 -.2.1 S.D. 21 31 11 17 17 17 15 15 10 9 7 7 4 5

TABLE C2 MEANS AND STANDARD DEVIATIONS OF PROCESSES AFFECTING POTENTIAL ENSTROPHY (x 10-18 sec'3) Function Month m = 1 2 3 4 5 6 7 8 9 10 11 12 1 14 15 Mean - 32.0 - 63.2 - 18.5 23.1 -165.6 - 2.4 28.2 - 26.5 - 20.5 4.0 - 47.0 71.1 65.3 35.4 148.6 February S.D. 379 325 364 363 254 262 363 366 355 340 264 279 259 177 211 I (nn,) Mrch Mean 1.6 - 89.5 - 62.7 16.7 39.4 - 18.3 - 73.0 - 2.0 27.4 83.1 8.3 72.1 43.5 9.0 -.0 S.D. 327 303 369 298 319 376 320 394 45 15 283 295 251 347 2 Mean - 4.0 - 42.8 - 47.4 - 35.0 14.1 -140.0 - 398 - 4.1 9.7 -.2 44.0 25.5 50.9 132.1 75.9 A l S.D. 246 335 388 250 252 363 343 355 335 410 320 242 245 205 236 Mean 8.1 11..8 - 8 - - 33.2 24.9 - 12.6 - 20.0 8.2 28.1 - 1.7 - 14.5 - 8.1 0.0 6.9 5.1 February S.D. 97 109 85 104 88 92 107 87 122 82 108 116 99 66 99 F ( n,) Marc Mean -233 - 6.5 21.4 - 7.2 13.7 - 87.5 69.1 23.7 - 29.2 34.2 - 23.7 - 40.2 20.8 15.2 9.7 n MS.D. 102 119 175 142 128 257 161 199 133 201 119 141 112 107 176 p Mean 33.6 - 19.3 - 5.5 - 29.6 12.7 47.0 - 3 - 5.1 - 18.4 - 12.9 -.2 23.8 - 13.1 - 17.5 19.1 S.D. 95 114 93 93 99 103 137 137 180 120 171 119 159 126 111 e Mean - 63.3 - 62.8 - 24.4 - 83.5 51.9 27.6 39.2 - 1.2 - 21.5 27.9 6.9 -.2 4.9 - 10 - 6.8 S.D. 128 135 158 133 139 137 149 105 92 93 50 33 29 24 6 I lo) March Mean - 6.5 -147.6 - 71.8 - 19.4 41.6 2.7 6.3 - 6.5 5.4 16.5 - 3.7 -.6 2.1 78 - 7 Ip (m S.D. 64 129 171 172 149 133 133 146 66 78 54 46 31 29 25 Mean - 15.6 - 87.5 - 40.2 - 16.5 - 28.3 - 39.5 3.8 31.8 12.2 - 6.6 - 16.6 - 4.9.7 - 1.0 3.1 S.D. 103 96 102 102 201 129 112 97 68 61 59 34 29 30 26 Mean 40.1 47.6 157.7 126.9 10.8 - 35.2 - 3.5 20.9 14.9 - 16.3 - 1.9 - 1.6 5.4 9.6 6.3 February S.D. 93 179 162 154 107 152 118 107 82 88 55 32 39 29 32 F (ml) Marh Mean 3.6 56.4 148.9 82.3 - 35.7 139 22.4 14.2 - 4.6 -10.7 15.4 - 4.3 4.9 - 4.9 Fpmo a S.D. 65 120 182 141 175 151 100 185 59 73 65 53 33 31 22 Mean 4.4 52.6 49.1 30.7 46.9 13.0 93 - 27.6 - 5.1 6.8 10.6 2.7 -.7 -.2 - p S.D. 87 117 82 80 159 137 117 97 75 58 48 34 24 36 25

APPENDIX D DERIVATIVE OF THE STREAMFUNCTION FROM THE POTENTIAL VORTICITY The potential vorticity tendency equation is integrated to produce a new potential vorticity field. This equation is of the form qi F(ti t1' ~ 3' H3'', v) i = 1,3 (1) For the subsequent integration for t+2 the streamfunctions'1 and Jt must be known at t+l. They are calculated from the following equations. V2 k2 =3 =72'V \V1 Ar ) q = V2* 3 + (4 -' v) (2) by decomposing both the potential vorticities and the streamfunctions. If we write ql F q1 F (m,n) q3 M M F3(mn) - z z y, iJ(mx + ny) 1 m=-M n=-M Amn) /0 #0 A3(m,n) 3 3 2f where J = and L is the length of a side of the square domain, then L we get for given values of m and n 129

130 22 2 2 2 Fl(m,n) = -[J2(m + n ) + ] Al(m,n) + x A(m,n) 2 2 2 2 2 F (mn) = -[J2(m + n ) + ] A (m,n) + A (m,n) 222 2 2 Let D =[J (m + n ) + % ], then 2 F (m,n) = -D A (m,n) + 2 A (m,n) 2 F (m,n) = -D A (m,n) + X A (m,n) which gives PF (m,n) + 2 F (m,n) A, (m,n)' 1 4 2 \ -D and D F (m,n) + X Fl(m,n) A7(m,n) =4 A -D

APPENDIX E ILLUSTRATIONS OF RESULTS FOR CASES B AND C Figures E.1 and E.2 are to be compared with Figure 23. Figures E.3 and E.4 are to be compared with Figure 24. Figures E.5 and E.6 are to be compared with Figure 25. Figures E.7 and E.8 are to be compared with Figure 26. Figures E.9 and E.10 are to be compared with Figure 27. Figures E.11 and E.12 are to be compared with Figure 28. 131

132.1.1 8 - 6 4 - KK e 2 oJ S / bo 10- / 8 0o 6. 4A e 2 - 109 - 240 480 720 960 1200 1440 t (hours or time steps) Figure E.1. Ke, Ae as a function of time for Case B.

133 10 11 lolO 8 6 4 K e 2 8 6 A e 4 2 8 10 8 6 240 480 720 960 1200 1440 t (hours or time steps) Figure E.2. Ke, Ae as a function of time for Case C.

134 1800oo -C(A K) e 1600 --- CA(m n,) 1400 D\ -*- D(Ae) \ 1200 - \ C) 1000 cg 800 \ \ (1, r \ F 400- / H/ I 200 - 200 - 600 - \ /. / -1200 4oo -1600 -1600-G for m = 8 G(Ae) = 50,540 -1800 - -C(Ae,Ke) = -19,830 CA(mln,f) = -8292 I I I I I I I I I I I I i 2 4 6 8 10 12 14 16 18 20 22 24 26 Wave Number Figure E.3. Mean values of CA(mln,f), -C(Ae,Ke), and D(Ae) for 300 _ t < 1500, Case B.

135 18018o - -C(A,K ) 160 - \ C(AlA) 140o D(Ae) \ e 120 \'o 100 8o \ C\ 80- 1 \l o! \ 60 / - 640 I 1 \20\ I I- I I 10 1 0' -60- \ 2 4 6 8 10 12 14 16 18 20 22 24 26 Wave Number Figure E.4. Mean values of CA(mfn,), -C(AeKe), and D(Ae) for 500 < t < 1500, Case C. Hio %% \I~~~~~~ —- -1. ~0' -140~ %.~/ ~ ~ ~ ~~~~~~o. ~ -~~60 ~~~~~~(A)=22

136 3600- C(A,K) 3200 - I -— CK(mln, ) I K 2800 -- D(K ) 2400 - 2000 - i- 1600 - / 1''0 I00 / 1 \ I 8 - 400Ioo -_ tJ i 800 - \ 1200- - 400'' I 400 - o -400for m - -3600 C(A,,K) = 19,830; -10CK(mn,) -11,528 | -2800 2 4 6 8 10 12 14 16 18 20 22 24 26 Wave Number Figure E.5. Mean values of C(Ae,K,), CK(m|n,), and D(K,) for 300 < t - 1500, Case B. CK(mlne) -11,o528

137 180 C (A K I \ 16o --- CK(mln, l) I I 140- D \ (K) I 120 I I 100 I I I 80 I o I I I I J J60 i 40 ) 20 - o\ / -.. <100 C C. 20 I C) 7 4.0' -60 - 7 -140 - -160 - for m = 8 C(Ae,Ke) = 1927 -180 - CK(mln,A) = -922 D(Ke) = -232 2 4 6 8 10 12 14 16 18 20 22 24 26 Wave Number Figure E.6. Mean values of C(Ae,Ke), CK(mln,A), and D(Ke) for 300 < t < 1500, Case C.

800 -' 600o -. -200p A v 400 0 m -00 - 0 c o for m = 8 I(m|n,e) = -5259. I I I I' I I I I I I I I I 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 52 Figure E.7. Mean values of I(mln,i), 300 < t < 1200, Case B.

8o cd 6 (D i) eI I -60 0 I(mln,l) = -4704 40 Figure E.8. Mean values o I(m)n,f ) = 00, C 4 2 4 6 8 10 12 14 16 18 20 22 24 26 28 50 52 Wave Number Figure E.8. Mean values of I)mjn,I), 300 < t < 1200, Case C.

14o 10 1010 - 8 6 4 2 m=8 109 - 8 c~ ~~r1~16 O 5 " 6) 2 10 8 24 6 4 2 10 - 240 480 720 960 1200 1440 t (time steps or hours) Figure E.9. The kinetic energy of wave numbers 8, 12, 16, and 24 as a function of time, Case B.

141 8 6 4 2 8 10 8 - Cj |i;12 16 107 8 6 4 2 24 1o6 240 480 720 960 1200 1440 Figure E.10. The kinetic energy of wave numbers 8, 12, 16, and 24 as a function of time, Case C.

142 10 9 I I O 108- ii) O -3 Wave Number Figure E.. K spectra at t = 900, 1500, Case B. t = 1500 6 ---- t = 900 10 2 4 6 8 10 20 30 Wave Number Figure E.11. Ke spectra at t = 900, 1500, Case B.

I I "" 10-0 v - o 109 4-D F 106..- -- t = 900. ry\ 2 4 6 8 1o 20 30 Wave Number Figure E.12. Ke spectra at t = 900, 1500, Case C.

UNIVERSITY OF MICHIGAN 3 901111111111111 111111111 11111103526 6942 3 9015 03526 6942